Field-theoretical description of the formation of a crack tip
process zone

Abstract
The crack tip process zone is regarded as a region where the solid physical properties
are altered due to high stress. They are controlled by the solid degrees of freedom
existing within the zone and vanishing outside, and can be divided into two classes:
(1) zones always existing at the tip and (2) those emerging as soon as certain conditions
are met. We focus on the zones of the second kind and argue that they can be described
analogously to phase transitions taking place locally. We report both a numerical and an
analytical solution for the process zone. We find that the zone can only exist within a
limited domain of the dynamic phase diagram, at one side of the phase transition line. We
describe this domain and establish its dependence on the crack velocity. We show the
existence of a critical crack velocity above which the zone cannot exist.



Mechanics of a propagating crack
To a great extent, the acute interest presently focused on the problem of the fracture of solids is related to high precision experiments that have recently revealed a number of reproducible instabilities in the high-speed dynamics and morphology of an isolated crack. It is especially intriguing that the instabilities often look similar in materials as different as glasses, PMMA and brittle polyacrylamide gel [1]. This has given rise to the hope that the crack dynamics may be considered within the context of the rapidly developing science of pattern formation, and has engendered a number of theoretical attempts, both analytical and numerical, to predict and describe such instabilities. These approaches analyzed the dynamics of a crack accounting for the elastic degrees of freedom of the solid, either linear or nonlinear.
Following the great achievements of the linear elastic fracture mechanics, it came to its limitations, being unable to determine the direction of the crack propagation in isotropic materials, correctly predict the terminal crack velocity, explain the mechanisms of the observed instabilities in the dynamics and morphology of a fast crack, in solids of various origins, both inorganic and organic, thus calling for its modification.
One such modification has been proposed by Buehler and Gao in order to explain the sub-and super-Rayleigh a e-mail: alexei.boulbitch@t-online.de as well as intersonic terminal crack velocities. They put forward the idea of a local material hyperelasticity at the crack tip [2,3]. The same idea explains some dynamic instabilities of a rectilinear crack [4]. Hyperelastic models are often applied to describe the deformation of polymeric networks [5]. It can, however, hardly be a universal model for inorganic solids.
The importance of the elastic nonlinearity near the crack tip has been revealed in brittle polyacrylamide gels, and established by the high precision measurements of the deviations of the actual crack tip profile [6] from the parabolic one predicted by the linear theory. These deviations have been described within the weakly nonlinear elastic approach [7].
Analytical studies of lattice models of crack propagation, such as those pioneered by Slepyan [8] and further developed by other researchers [9,10], have demonstrated the great complexity of the processes following a crack propagation. However, the lattice models admitting analytical solutions suffer from making assumptions that are too specific [9]. In particular, the bond breaking criterion generally employed in these models is usually postulated as a nonphysical deterministic condition of breaking upon exceeding a critical bond length [9,11]. This contradicts the probabilistic nature of the bonds' breaking, and ignores the role of the local atomic structure at the tip. In contrast to such an assumption, computer simulations as well as numerous experiments (see below) show that in general this is not the case. Further, analytical calculations [12] demonstrated a counterintuitively high sensitivity of the results to the choice of model lattice potentials. The latter has already been established within a relatively narrow set. In addition, so far no approach to correctly introduce dissipative processes in a deterministic form at a truly atomistic scale has been proposed [11]. This suggests that the universal description of the experimental facts of the crack dynamics fails within the tractable lattice models.
In a series of papers, the explanation of an apparent universality of dynamic crack instabilities has been suggested using various 2D versions of the phase field approach [13][14][15], where the crack is regarded as a slightly diffused cut in the elastic medium. Indeed, simulations have demonstrated the existence of oscillatory and branching instabilities taking place at the crack velocities exceeding some threshold [13,14,16,17]. It was later revealed, however, that these instabilities differ from those actually observed in experiments [10]. In particular, the phase-field models predicted the wave length of the oscillatory instability to be proportional to the dimension of the system [16,17], while experiments have demonstrated its intrinsic character, independent of the system geometry [18]. An essential difference also exists between the experimentally observed microbranching instability [19,20] exhibiting a 3D geometry [18,21,22] on the one hand, and its predicted 2D phase-field analogue [13,17] on the other hand.
The results mentioned above suggest that elasticity alone, linear or nonlinear, cannot explain all the observed phenomena of crack dynamics. This suggests that some additional degrees of freedom (different from the elastic ones) may be responsible for the anomalies. This implies that such nonelastic degrees of freedom must be explicitly included in the theoretical models. It is reasonable to expect these degrees of freedom to manifest themselves within the process zone (PZ), a nano-to meso-sized domain at the crack tip. The same conclusion follows from a great body of experiments directly demonstrating the activation of nonelastic degrees of freedom within the PZ.

Experimental observations of crack tip process zones
For a long time, the PZ was beyond the reach of experiments. During the last decades, new experimental techniques have emerged, that have significantly changed this situation. It appears to be possible to obtain the PZ structure with a high spatial resolution, combining several techniques in one study, in some cases also appropriate for propagating cracks. One should mention the method of high-angle annular dark-field scanning TEM [23] allowing the direct imaging of atomic locations. To exclude ambiguities, the results can be further combined with electron nanodiffraction patterns [24]. A combination of micromechanical loading with in situ high-resolution X-ray microdiffraction [25] enables the imaging of the PZ in polycrystalline Nitinol. A combination of an in situ SEM with electron backscatter diffraction [26] has made possible the study of the evolution of the emerging phase in the tip vicinity during fatigue experiments. An in situ optical digital image correlation technique allows monitoring the strain field and phase boundary configuration at the tip of a propagating crack [27]. Raman mapping has revealed the local distribution of phases in the tip vicinity [28]. Atomic force microscopy enables one to map the lateral distribution of the surface height directly related to the spontaneous strain, thus making it possible to detect the phase boundary [29,30]. Nanoindentation allows one to study the temperature dependence of the zone size [31]. The above methods together with the more traditional ones [32] have revealed that the solid structure in the close vicinity of the tip undergoes chemical or structural reconstructions. The latter often take the form of conventional phase transformations localized in the vicinity of the tip, referred to below as the "local phase transformation" (LPT).
The great variety of materials in which transformational PZs have been observed as well as the broad range of transformations observed within the PZs suggests that it is a typical phenomenon, occurring for many solids. It should be expected that this phenomenon would have an impact on the static and dynamic properties of a fracture.

The theoretical analysis of local phase transitions
In the theoretical studies of transformational PZs, three distinct approaches can be pointed out. On the one hand, the atomistic mechanisms of LPTs have been studied by computer simulations and density functional theorylike calculations for several solids, such as iron [81][82][83][84][85][86], silicon [87][88][89][90], tantalum [91], zirconium [92], UO 2 [93], molybdenum [24] and Nitinol [94]. This became possible as the result of the development of computational power on the one hand and the advancement of the molecular dynamics approach on the other hand, as well as by the implementation of hybrid approaches combining molecular dynamics with quantum mechanics [90,[95][96][97][98]. The simulations revealed a strong dependence of the LPT formation upon (i) the loading mode; (ii) the crack plane direction and (iii) the sample geometry [81,82]. They further elucidated the atomistic mechanisms leading to the development of a LPT [82,83,90].
The analytical continuum approach has been developed based on the assumption that the LPT zone only differs from the rest of the solid by its (i) elastic moduli and (ii) the existence of a spontaneous strain [32,[99][100][101].
The above approach exhibits, however, two main deficiencies. The first is related to the proper choice of the criterion of the phase boundary. By analogy with the macroplasticity criterion, most papers assume that the phase boundary is situated in the place in which a stress component (see, e.g., Ref. [102]), or a combination of several stress components (e.g., Refs. [103][104][105][106]), exceeds a critical value σ c , the latter being regarded as a material constant. It has been, however, experimentally shown that accepting this assumption, this parameter is highly anisotropic [27]. In addition, it should exhibit a dependence upon temperature, or more generally, on the position in the solid phase diagram. This shows that the plasticity analogy [102][103][104][105] is an oversimplification of the LPT boundary problem.
It should be noted that there is an alternative to this simplistic criterion. A boundary condition corresponding to a narrow phase boundary has been rigorously derived by Roitburd [107,108] and Grinfel'd [109,110], leading to a nonlinear integral equation on the interface configuration [107,108].
In addition, the phase transformation criterion has been derived thermodynamically at the macroscopic scale for a representative volume rather than for a single interface [105,106], enabling one to efficiently apply a numerical approach. Macroscopic equations for the representative volume based on the local conditions at the interfaces have been derived in references [111,112].
The second deficiency is related to the absence of an explicit treatment of internal degrees of freedom condensing at the tip and becoming excited during the crack motion. The intrinsic peculiarities of the static and dynamic behavior of these degrees of freedom influencing the crack are ignored in such an approach.
The third approach is based on the Landau theory of phase transitions [113] applied to solids with an inhomogeneity introduced by a spatially-dependent stress field. Nabutovskii and Shapiro were the first to describe the equilibrium LPT using the Ginzburg-Landau equation [114]. The possibility of applying this approach to analytically describe the crack tip PZ was briefly indicated in reference [115] in the cases of a motionless and of a steadily propagating crack. The paper [116] described the crack's self-oscillation generated by the PZ of the second kind, [117] described the crack velocity jumps engendered by the PZ that emerged by a second order phase transition, while [118] did so by studying a manifested first order LPT. Numerically, the PZ of the second kind has been addressed in references [119][120][121][122].
In the present paper, we argue that the PZ can be generally regarded and described as an LPT. The description of the PZ requires the explicit introduction of an order parameter, which is an internal degree of freedom distinguishing the zone from the rest of the solid. We introduce such an order parameter subjected to the time-dependent Ginzburg-Landau equation. To be precise, we worked within the example of the local condensation of optical phonons at the crack tip. Such a condensation corresponds to a structural phase transformation localized at the tip. We regard the emergence of a dynamic PZ as the bifurcations of the coupled equations of motion for both the mechanical subsystem and the order parameter. This approach enables one to characterize both the spatial structure and shape of the PZ as well as its dynamic properties.
We give a numerical finite element solution of the equations of motion exhibiting the bifurcation. We analytically derive the exact condition of bifurcation, depending on the crack velocity in a wide range, including high crack speeds. We demonstrate that the PZ does not always exist, but emerges at the crack tip as soon as a certain thermodynamic condition is met. We find an exact expression for this static condition in terms of the temperature and stress intensity factor. We describe the emergence of the zone also in the case of a propagating crack and show how it depends on the crack velocity. We find a critical velocity beyond which no zone can exist at the tip. We further obtain an analytical, asymptotically exact solution for the order parameter describing the PZ in the vicinity of the bifurcation. This is done for the equations of the crack/PZ motion corresponding to a rather wide class of PZs.
We further give numerical estimates for the temperatures corresponding to the emergence of the PZ. These estimates show that the predicted phenomena occur in a wide class of materials, under easily realizable fracture conditions. This paper is organized as follows. In Section 2 we formulate a dynamic equation for the order parameter. In Section 3 we present our numerical results. In Section 4 we find an analytical solution of the PZ equations. In Section 5 we make estimates and comments generalizing our findings. Section 6 summarizes the results.

The order parameter
The very notion "process zone" implies that a small domain in the immediate vicinity of the tip of a brittle crack has specific properties, different from those of the bulk of the solid. The PZ can only be distinguished from the rest of the solid if at least one of its physical properties varies perceptibly across its boundary. This may be either an abrupt quantitative change, such as a steep growth of the elastic nonlinearity [123] or hyperelasticity [2,3] at the tip, or a qualitative change. In the latter case, the PZ differs from the bulk by, e.g., its chemical composition, electronic properties (such as, e.g., metal/isolator, exciton condensate/exciton gas, etc.) or crystal structure.
We parametrize the differences of the solid properties inside the zone from those outside by a field, η = η(r, t), referred to as the "order parameter". Without loss of generality, one may assume η = 0 inside the zone, while vanishing outside. In this respect, our approach is akin to the theory of phase transitions [113] as well as to the popular phase-field approach [124][125][126]. In contrast to the latter, our order parameter explicitly describes the degrees of freedom within the PZ responsible for the variation of its symmetry, crystal structure, or chemical composition.
At present, the terms "phase transition" and "order parameter" unify a crystal structure variation of solids with bifurcations in nonlinear systems, both equilibrium and non-equilibrium, such as, e.g., bifurcations taking place during chemical reactions [127]. This implies that each of these transitions can be described by its inherent degree(s) of freedom responsible for the bifurcation of the nonlinear system or for the structural variation, namely, the order parameter. Let us look at the PZ from this point of view.
With respect to the formation of the PZ, all order parameters split into two classes. Under symmetry transformations, an order parameter belonging to the first class transforms according to the same representation as the strain tensor, ε = ε(r), or some combinations of its component(s). In this case, the equation for the order parameter (see below) admits the absolute term ∼ε, and the PZ exists at all values of temperature and stress intensity. Far from the tip, such an order parameter has the asymptotics η ∼ ε(r). Formally, the crystal structure of such a zone can be obtained from that of the bulk by superimposing a certain strain corresponding to the order parameter.
In the second class, the order parameter and the strain tensor components transform according to different representations of the symmetry group. Therefore, the strain has an indirect effect on the emergence of the order parameter. For this reason, the PZ only emerges as soon as a certain condition is achieved, such as a specific temperature or value of the stress intensity, or both. In other words, the behavior of the PZ exhibits the features of a phase transition or of a bifurcation.
Such a classification can be applied to the cases involving the group-subgroup relation between the PZ and solid bulk symmetries, or if the PZ and the matrix symmetries are subgroups of a common group. The transitions referred to as "reconstructive" exhibit no such relations, and require a special modification of the Landau theory [128] . In particular, the family of martensitic transformations occupies a position intermediate between these two classes. On the one hand, the symmetry changes during such transformations are described by the so-called "transcendental" order parameters, which are different from the strain [129,130]. On the other hand, the transitions are accompanied by spontaneous strains, which are typically large, admitting an approximate treatment of the problem similar to that in proper ferroelastics [131].

The free energy and dissipation function
We study the second class of PZs, in which the order parameter is generally a multicomponent object η = (η 1 , η 2 , . . . , η n ) (where n = 1, 2, . . . is a natural number) representing a set of normal coordinates of an optical phonon describing displacements of the atoms within the crystal's unit cell.
The order parameter η(r,t) = 0 contributes to the solid free energy: (1) while its time variation gives rise to the energy dissipation described by the dissipation function D: is the free energy density and ε ≡ ε ik is the strain tensor defined as follows: and Ω is the spatial domain. Since we consider the case of a thin plate here, Ω represents the infinite plane: dΩ ≡ dxdy. We, thus, assign both F and D to the unit solid thickness in the z direction. Further, κ is the kinetic constant and t is the time.
Here we address the simplest case, a phase transition described by a one-component order parameter (n = 1). The latter catches the basic features of the PZ. In this case, the only possible transformations under the action of the crystal symmetry group of the solid are either η → η or η → −η. Being invariant with respect to the crystal symmetry, the free energy and the dissipation function should only contain even functions of η [113]. The free energy density can be written in the following form: Here the function Φ pt (η) denotes the part responsible for the phase transition: where ∇η is the order parameter gradient, g > 0, and β 0 > 0 are the constant parameters of the Landau potential (1). The case of the first order transition β 0 < 0 will be analyzed elsewhere. Expression (4) coincides with the classical Landau theory describing bulk phase transitions. The latter theory assumes that only the factor α depends on the temperature: α = a(T − T c ), where a > 0 is a constant, T is the temperature, and T c is the Curie temperature [113]. Though in general all these parameters may depend upon the temperature, we stick here to this assumption; its generalization is straightforward. Φ el (ε ik ) is the elastic part of the free energy density, and we use the elastically-isotropic approximation: where E is the Young's modulus and σ is Poisson's ratio. The latter should not be confused with the stress tensor σ ≡ σ ik .
Finally, the contribution Aη 2 ε ii to (3) takes into account the interaction between the strain and order parameter fields, the striction constant A being either positive or negative. It should be noted that the form of the interaction term Aη 2 ε ii in (3) implies that the phase transition only gives rise to a spontaneous dilatation. This is always the case with the one-component order parameter considered in the present paper. It should be noted that in the case of a multicomponent order parameter, this term may have a form admitting deviatoric spontaneous strain. The striction term in this case admits no universal form and should be separately derived for each specific case.
The striction constant is directly related to the slope of the phase transition line k = dT c /dp in the (p, T ) phase diagram. Indeed, one can represent the term ∼η 2 in equations (3) and (4), by It should be mentioned that the constants g, a, β 0 , T c and A (or k) related to the phase transformation together with the elastic constants E and σ are material constants of the solid controlling its fracture behavior. The equations of motion of the order parameter can be built on the basis of the free energy (1) and dissipation function (2) in the standard way [113]: where δ is the sign of variation. One obtains the following system of equations: For simplicity, we have omitted the inertial term in the second equation (8), which is valid for any motion much slower than the speed of sound. Here Δ is the Laplace operator, and σ ik = ∂Φ/∂ε ik is the stress tensor: Here δ ik is the Kronecker symbol. The last term in (9), Aη 2 δ ik , describes the spontaneous stress generated by the phase transition. Equations (8) and (9) constitute the complete system describing the dynamics of the PZ.

Elimination of the elastic variables
Making use of (8) and (9), one can express the displacement vector u i as where u (0) i (r) is the displacement field created by the "undressed" crack (that is, the crack without the transformational PZ) and G ij (r) is the elastic Green's function of the solid with a cut. Using the Fourier transform one finds an integral representation of the dilatation ε ii (r): where the strain ε ii (r) is generated by the "undressed" crack, i.e., the one without the PZ (η ≡ 0). In the elastically-isotropic case, it is given by the the well-known expression where r and θ are the polar coordinates with the origin at the crack tip (see, e.g., [132]) and K I is the stress intensity factor. The second term (12) represents the contribution of the PZ to the strain. Substitution of the strain expression (12) into (8) yields the nonlinear, integro-differential equation of motion for the order parameter: where the nonlinear equation part,N (η), is expressed in terms of an integral: Equations (14) and (15) involve no elastic degrees of freedom and contain only the field η(r). It should be stressed that the derivation of equation (14) involves no approximations.
To the best of our knowledge, no explicit form for a Green's function of a body with a cut is known. We approximate it by the Green's function of the infinite elastically-isotropic body [133]. In more details this approximation and the approach to improve it is analyzed in Discussion (Sect. 5). In this case one finds and the identity holds. Making use of this relation in (12) one finds the strain dilatation in the form Substitution of the strain trace (17) into the free energy (1), (3), (4), and (5) yields the effective free energy with respect to the initial parameter β 0 . Substituting (13) into the first equation of (8) gives a nonlinear differential equation of motion: In this caseN (η) = β 1 η 3 (r) is polynomial. It should be mentioned that equation (20) can also be obtained by the variation procedure (7) using the dissipation function (2) and the effective free energy (18). Equation (14) (or (20) in the elastically-isotropic approximation) describes the dynamics of the PZ of the second kind. Its solution will be addressed in the next sections.

The automodel regime
Assuming a crack tip propagating rectilinearly with velocity V along the Ox axis, one finds η = η(x − V t, y), and passing to the comoving frame, x = x − V t, y = y, the equation of motion for η = η(x , y ) (20) takes the form where One chooses the sign "+" if k > 0 and "−" in the opposite case. Further, r = x 2 + y 2 1/2 , and the Laplace operator is defined by Δ = ∂ 2 /x 2 + ∂ 2 /y 2 . Since in the following we only use the comoving frame, from here on the primes on x, y and r will be omitted. If one describes a PZ η(r) = 0 embedded in the matrix of the bulk phase (α > 0, η = 0), localized at the crack tip, (x, y) = 0, the boundary condition takes the form η(∞) = 0. Because the variation of atomic coordinates must be limited, the order parameter is everywhere finite: To describe a PZ embedded in the matrix of the daughter phase (α < 0; η = (−α/β) 1/2 = 0), one needs to use the boundary condition η(∞) = (−α/β) 1/2 .

Simulations
Below, we study the problem at hand by simulation, using the finite element method.

Rescaling
We report below the simulation of the daughter phase PZ embedded in the matrix of the mother phase adopting the elastically-isotropic approximation (20).
The parameter fixes the characteristic lateral size of the distribution, as will be explicitly shown in Section 6. Equation (21) witĥ N (η) = β 1 η 3 has been rescaled to make the variables dimensionless and minimize the number of control parameters: The rescaled equation (21) takes the form where Δ 1 = ∂ 2 /∂x 2 1 + ∂ 2 /∂y 2 1 , r 1 = (x 2 1 + y 2 1 ) 1/2 , while the dimensionless control parameter q and dimensionless velocity ν are expressed as follows:

Results
We used a pseudo-time stepping approach, representing a version of the iteration method. Its description, along with the technical details of the software and settings, are given in Appendix A. At large positive values of q, equation (25) exhibits a trivial solution u = 0. Below a certain value q c = q c (v) found by trial and error, a non-trivial solution emerges. Figure 1 shows the distributions of the rescaled order parameter u at the tip (x 1 = y 1 = 0) of a motionless crack (ν = 0) for successively increasing values of the difference q c − q > 0 below the bifurcation point q = q c .
The distribution of the rescaled order parameter u(x 1 , y 1 ) at the tip of a propagating crack (ν ≥ 0) at the same value of q is shown in Figure 2. One can see that the motion of the crack transforms the distribution in several ways. First, with an increase in the dimensionless velocity ν, the distribution becomes lower and vanishes as soon as the velocity reaches a certain value. Besides, the order parameter distribution is compressed in front of and stretched out behind the tip in comparison to that at the tip of a motionless crack.
To find the points of bifurcation q c = q c (ν), we made simulations by fixing the value of ν and evaluating  the maximal order parameter value u max at different values of q. One can see that at large values of q, u max vanishes, while emerging beyond a certain threshold (Fig. 3).
The values of these thresholds have been extracted from the above data by fitting to the function where u 0 and q c are the fitting parameters, and q is the variable. Figure 4 shows the obtained dependence q c = q c (ν).
To summarize, we studied numerically the PZ containing the daughter phase, η(r) = 0 embedded into the matrix of the mother phase, η = 0. We have shown that this zone emerges as soon as q ≤ q c (ν), but vanishes at large q and ν. We numerically obtained the q c (ν) dependence.

Bifurcation theory
The analytical results below depend heavily on using the approach of the bifurcation theory of nonlinear equations. For the convenience of the reader, let us first briefly recall the key results of this theory [134], which we use in the following argumentation.
In a close vicinity of the bifurcation point, one can obtain the asymptotically exact, overcritical solution of equation (28) in the form of a series in terms of two small parameters: the amplitude ξ, and the "distance" from the bifurcation point, α − α * . The bifurcation theory ensures, however, that its main term always takes the form where the amplitude ξ will be determined from the nonlinear equation (28) [134].
To obtain ξ one can make use of (29) and the main term of the solution (30) to represent (28) in the form ξ ×[L(α)−L(α * )]Ψ * =N (ξΨ * ). Multiplying both its parts by Ψ * and integrating one finds where we write f, g = f (r)g(r)dΩ for the Hilbert scalar product of two functions, f (r) and g(r). Equation (31) is a nonlinear equation with respect to ξ, only valid if ξ is small, and is referred to as the "ramification equation" [134]. Its solution yields the amplitude ξ, thus finalizing the solution (30) of the bifurcation problem valid at (α − α * )/α * 1. One observes that although the above recipe heavily involves the solution of the linear equation (29), it represents the solution of the nonlinear equation (28).
Full details, theorems, and their proofs can be found in the book of Vainberg and Trenogin [134].

The bifurcation condition
Let us turn to the analytic analysis of the automodel equation (21), of the PZ motion. We start with the case of the crack in the matrix of the mother phase η = 0, and will look for the emergence of the daughter phase η = 0 at the crack tip. Accordingly, let us assume α > 0 and look for the point of instability of the trivial solution η = 0 describing the homogeneous mother phase.
Comparing (28) and (21) one establishes the correspondencê whileN (η) is either given by (15) in the elasticallyanisotropic case, orN (η) = β 1 η 3 in the elasticallyisotropic case. Let us apply the above recipe of bifurcation theory [134] to the case at hand. First of all, one finds that at k > 0 (corresponding to the sign "+"), equation (29) has no discrete spectrum, meaning that the solution η = 0 of equation (21) is stable at α > 0.
At large values of α, equation (21) has the trivial solution η = 0. Below the bifurcation point, in its close vicinity one finds a non-trivial solution η = η(r) = 0 and σ ik (r) = σ (0) ik (r)+O(η 2 ), where η is small. According to its definition (4), α is expressed in terms of the temperature: α = a(T − T c ). The value α * will, thus, yield the temperature T * at which the PZ emerges at the crack tip. The exact solution of equation (29) withL(α) (32) is obtained in the next section.

Exact results for the eigenvalues and eigenfunctions of equation (32)
The solution of (29) plays an important role in the theory, defining both the bifurcation point α * = α 1 and the spatial distribution of the overcritical solution (30) of the nonlinear equation (28). Let us solve it.
The linear part of equation (21) takes the form gΔΨ n + κV × ∂Ψ n /∂x − α n − B cos(θ/2) one proceeds to the equation in terms of ψ n (r): Let us agree that in the case of the daughter PZ embedded in the bulk mother phase, all parameters characterizing the problem, such as the characteristic size R, the bifurcation point, α * , etc., will be written with the subscript "1": R 1 , α * 1 , etc., while the subscript "2" is reserved for the case of the mother PZ embedded into the daughter phase matrix. Passing to dimensionless cylindrical coordinates θ and ρ = r/R 1 , one transforms (35) into the following equation: where λ n (n = 1, 2 . . .) represent the eigenvalues of equation (36). They are related to α n as follows: Equation (36) represents a 2D Schrödinger equation with the anisotropic potential U (ρ, θ) = − cos(θ/2)/ √ ρ shown in Figure 5.

Bifurcation point
The above solution is valid for k < 0. It gives the ground state eigenvalue (41) and the eigenfunction Figure 6A where ν = V/V c is the dimensionless velocity, and the critical velocity V c is defined as Its physical meaning will be discussed below. Note that the expression for ν coincides with the one given in equation (26). Setting the exponent to −1, one finds the size L f of the order parameter distribution in front of the tip: At the tip of a motionless crack (ν = 0), one finds the order parameter distribution size L f 1 ≈ 9.25R 1 decreasing down to L f 1 ≈ 2.96R 1 at ν = 1. The dependence L f 1 = L f 1 (ν) is shown in Figure 7.

The overcritical solution
Substituting the eigenfunction (42) as well as (32) into the ramification equation (31), one obtains where the factors I n (n = 2, 4) are the integrals over the whole plane depending on ν.
It should be noted that substituting (30) and (42) into the effective free energy, (18) yields the PZ free energy: As expected, its minimization with respect to ξ brings one back to the ramification equation (45), thus giving an alternative way to build the overcritical solution.
The solution of the ramification equation has the form The integrals I n (ν) cannot be obtained analytically at ν = 0. We calculated the ratio I 2 /I 4 numerically by using the standard NIntegrate routine of Mathematica 10.
This completes the building of the asymptotically exact overcritical solution (30) of equation (21).

The critical velocity
The result (48) implies that the PZ only takes place at α < α * 1 , while at α > α * 1 , the crack tip is undressed. Since at α < 0 the whole bulk of the solid transforms into the daughter phase η = (−α/β) 1/2 = 0, the domain in which PZ takes place is restricted to the interval 0 < α ≤ α * 1 .
Since α * 1 = α * 1 (V ), the latter result defines the velocity V at which the PZ disappears. Making use of equation (40), one finds that this happens at V = V c . In other words, the PZ can only exist at the crack tip for V < V c , while for V ≥ V c it vanishes. This property is fundamental for any PZ, both of the second and of the first order. It is a direct consequence of the fact that the order parameter has its own dynamics, exhibiting an intrinsic characteristic time, and that as soon as V ≥ V c , the order parameter in front of the crack tip has no time to evolve from η = 0 to η ≈ ξ.
The spatial distribution of the order parameter (30) is shown in Figure 6. Here (A) shows the order parameter in the vicinity of a motionless crack tip, while (B) shows it in the case of a propagating crack. The latter image is obtained with the velocity value V = 0.5V c .
One finds that the order parameter is highly localized in the vicinity r 10R 1 of the tip of a motionless crack. In the case of a moving crack tip (Fig. 6B), the order parameter distribution is compressed in front of the tip and stretched at its back compared to that of a motionless tip. The length L b of the order parameter distribution behind the tip of a propagating crack takes the form The length L b diverges if V → V c . Simultaneously, as V → V c , the amplitude ξ vanishes. Summarizing the results of the present section, we have shown that in a solid with a crack, the trivial solution η ≡ 0 is stable at high temperatures (T > T c ), but loses its stability at the point α = α * 1 > 0 corresponding to a temperature T * 1 somewhat higher than that of the bulk phase transition, T * 1 > T c . This solution describes a region of the phase η = 0 embedded in the matrix η ≡ 0 representing the transformational PZ at the crack tip. In terms of temperature, its existence is limited to the domain T c < T < T * 1 , while in terms of velocity, it is limited to 0 ≤ V ≤ V c . The transformational PZ at the tip of a propagating crack is deformed in comparison to that of a motionless crack: the order parameter distribution is compressed in its front, while stretched out at its back.

The mother phase process zone embedded in the bulk daughter phase
Let us assume that the whole solid is in the daughter phase η = const. = 0. This takes place at α < 0. In this section we study the emergence at the crack tip of the mother phase PZ embedded in the matrix of the daughter phase.
Let us first find a solution for the order parameter η = η(r) in the daughter phase. We will assume here that α is negative and large in its absolute value, so that this solution η(r) is stable. One observes that at r (g/α) 1/2 the terms ∼Δη and ∼∂η/∂x in equation (21) are much smaller than the others and, therefore, one can neglect them. This yields the approximate dependence η = η 0 (r) in the daughter phase: To be specific, from the two solutions of the state equation, we have chosen the positive one. Let us now look for a perturbation of (51) in the form η(r) = η 0 (r) + δη(r).
Substituting it into (21) one finds that the term δη(r) obeys equation (28) witĥ (52) and In this case, at k < 0 -the sign "+" in (52) -equation (52) appears to have no discrete spectrum, implying that the solution (51) is stable. The discrete spectrum indicating the instability, however, exists for k > 0, corresponding to the sign "−" in equation (52). Below we consider this latter case. Let us agree that in the case of the mother PZ embedded in the bulk daughter phase all parameters will be written with the subscript "2": R 2 , α * 2 , etc.
Applying the analysis already described above to the present case, one finds the eigenfunction (42), in which instead of R 1 one should take a characteristic size R 2 expressed as: which is smaller than (23) by the factor 2 −2/3 . The expression for the bifurcation point α * 2 has the form The latter is smaller than α * 1 by a factor of 2. Setting α * 2 = 0 one can make sure that the critical velocity V c has in this case the same value (43) as in the case of the daughter PZ.
Since the overcritical solution (30) for δη is determined by the same eignfunction (42) as that for the hightemperature case, the distribution δη(r) = ξ 2 Ψ * (r) has the same form as that of the order parameter in the hightemperature PZ. The distribution is shown in Figure 6. Analogously to the high-temperature case, the distribution δη(r) at the tip of a propagating crack is compressed in its front, but stretched in its back, and the relation (50) holds. The amplitude ξ 2 should be determined from the ramification equation.
In the low-temperature phase, the ramification equation is most conveniently obtained starting from the effective free energy, which has the form Details on the calculation of the factors I i (i = 2, 3, 4) can be found in Appendix B.
The effective free energy (56) has a cubic term. Since it is positive, one finds that the left minimum of the free energy (56) is more pronounced (Fig. 9).
Analogously to the static state, one concludes that one must only choose the solution of the ramification equation which corresponds to this deeper minimum. In other words, it is the negative solution that takes place. This solution of the ramification equation takes the form In general, the sign of δη is opposite to that of η 0 . Summarizing the findings of this section, the inhomogeneous solution η 0 (r) describing the order parameter distribution in the low-temperature phase is stable at α < α * 2 < 0. At α = α * 2 < 0 this solution becomes unstable and the solution η 0 (r)+δη(r) branches off, the signs of η 0 and δη being different. This means that the zone with the high-temperature phase, η ≡ 0, emerges at the crack tip at α = α * 2 and exists within the domain 0 ≤ α ≤ α * 2 provided 0 ≤ V ≤ V c .

Symmetry aspects
The present theory does not apply to materials belonging to the first class PZ. The latter includes the so-called proper ferroelastics, for which the strain plays the role of the primary order parameter. This requirement means that (i) the symmetry variation affects only the point, rather than space symmetry group and (ii) that the crystal structure of the daughter phase can be obtained as that of the mother phase subjected to some homogeneous strain. These requirements are only fulfilled for a few dozen materials [136,137]. Though rather limited, this class includes several materials very important in applications, such as pseudo-elastic nitinol [138,139], as well as other materials exhibiting proper martensite-austenite transitions.
Instead, the theory developed in this paper applies to the so-called improper ferroelastics, where the primary order parameter differs by its symmetry from strain. The latter class lacks the symmetry limitations of the ferroelastics and, therefore, is much more numerous.

Dilatant versus nondilatant case during elimination
In this paper we only considered the case of (a) a onecomponent order parameter and (b) dilatant spontaneous strain. Our analytical results can be straightforwardly generalized to the case of a multicomponent order parameter and non-dilatant strain.
Indeed, the key point in the calculation is the elimination of the elastic degrees of freedom from the equations of motion, leading to a single nonlinear equation formulated only in terms of the order parameter. This latter equation can be analyzed by means of bifurcation theory. One can see that if non-dilatant terms ∼εη 2 are present in the free energy, one can still eliminate the elastic degrees of freedom by an approach analogous to the one described in Section 3.2 with a different kernel of the integral (15). This does not allow, however, the calculation to be carried out in an explicit form like (17). As a result,N (η) in this case will be a nonlinear integral operator rather than a polynomial. It still will be of the third order in terms of η. This makes the problem technically more intricate, but still tractable.
From the qualitative point of view, taking into account deviatoric terms will change the eigenfunction form yielding the PZ shape, and the numerical factor λ 1 in front of the eigenvalue (39) defining the temperature shift ΔT * . The functional forms of the dependencies of ΔT * and V c on the parameters K I , k, κ, etc., however, are only determined by the dependence of the "potential" of the Schrödinger equation, equation (36), on the distance U ∼ ρ −1/2 , and are, therefore, universal.

Improvement of the Green function approximation
To eliminate the elastic degrees of freedom we used the Green function (16) of a continuous elastically-isotropic medium [133]. This is, however, an approximation, since we apply it to the solid with a crack. The correct Green function is somewhat "softer" then (16). It should account for an additional opening of the crack in response to the stress field generated by the PZ.
In the paper [102] the authors have shown that the stress at the banks of the cut generated by the dilatant phase transformation at the crack tip is strictly zero. This is, however, not the case in our approach. The difference arises due to the approximation of an infinitely thin phase boundary accepted in reference [102]. Our approach is free of this approximation. Therefore, in the present theory the order parameter is distributed backwards along the banks of the cut (x < 0, y = 0) as follows: and according to (17) gives rise to the strain and the corresponding stress there. This stress contributes to the perturbation of the stress intensity factor of the crack. Let us note that the amplitude ξ = ξ(K I ) depends on the stress intensity factor. This gives rise to a self-consistent nonlinear equation on K I yielding the effect of the additional crack opening. The approach formulated above, thus, accounts for the cut and improves the above approximation for the Green function. Its effect is, however, small controlled by the value of the dimensionless parameter R 1 /L c (1 − ν) 1, where L c is the crack length. This calculation, though technically simple, is lengthily, and to simplify reading we will publish it elsewhere.

Small versus high order parameter
In this paper we only addressed incipient PZs, where the order parameter is small. This enabled us to apply bifurcation theory. The theory can be straightforwardly generalized to the case of the so-called first order transitions close to the second order when the order parameter undergoes a small discontinuity in the transition point. This generalization is, however, limited to a close vicinity of the tricritical point. Though such points occur in nature, they are relatively rare. For this reason we did not include these cases in the present paper. On the other hand, phase transformations of the pronounced first order cannot be treated by the approach of bifurcation theory and our approach does not apply to these cases.
Second order phase transitions are generally considered to be soft. The softness of the second order phase transitions only means that the solid continuously passes from the state η = 0 to the state η = 0 at the transition point. This differs from first order transitions, where at the transition point the order parameter exhibits a jump. The impact of the PZ on the crack dynamics may, however, be observed not only in the close vicinity of the bifurcation point, but also far from this point. The effect of the PZ depends on the absolute value that the order parameter can achieve within the PZ zone. In the discussion below, we will argue that the zone is typically very wide. This implies that the PZ order parameter is most often in the saturated state. In this state the effect of the second order PZ exhibits no qualitative difference from that of the first order PZ.

Temperature shift as a material constant
The PZ can be characterized by the temperature shift ΔT c = ΔT * (K IC ) taking place at K I = K IC , at which V = 0: depending only on material constants of the solid. ΔT c can, thus, also be regarded as a solid material constant, representing its important characteristics. In the present paper this expression has been obtained by means of bifurcation theory in the case of a second order phase transition.
Let us now show that it also holds qualitatively for any transitions of any order. The vast majority of solids exhibit first order PTs characterized by phase diagrams, the (p − T )-planes divided into regions of the existence of phases, typically about ∼100 K wide. These regions of existence are separated by the lines of the first order PTs, T 0 = T 0 (p), flanked by hysteresis domains where the merging phases coexist. The latter are usually ∼1 to 10 K wide, only achieving ∼100 K for few a materials, typically belonging to the first class, such as nitinol or ZrO 2 [140].
The PT lines are approximately straight, characterized by their slopes k = dT 0 (p)/dp [140]. The application of a mean stress of Δp = −σ ii shifts the PT temperature by ΔT = kΔp, where ΔT = T − T 0 . Here T is the temperature of the phase transition under the pressure Δp.
The stress σ ii ∼ K IC r −1/2 in the crack tip vicinity yields, thus, a PT temperature shift ΔT ∼ −kK IC r −1/2 at the point specified by the distance r from the tip. Let us introduce the order parameter correlation length, r c = [g/a (T − T c )] 1/2 , the distance over which the variation of the order parameter is small [113]. One observes that the PZ size cannot become smaller than the correlation length r c . Therefore, one finds the temperature "distance" ΔT c = T c − T 0 from the bulk transformation temperature T 0 at which the PZ emerges by setting r ∼ r c in the above expression for ΔT : where r c0 = (g/a) 1/2 is the order parameter correlation radius at 1 K from the transition line. This brings us again to the expression (59).
One concludes that the interval (59) between the bulk transition temperature and that of the emergence of the PZ is universal for the zones of the second kind, independent of the order of the phase transformation.

Typical values of the temperature shift
The universal nature of the temperature shift enables finding the order of magnitude typical for ΔT c . The typical value of the fracture toughness of the inorganic solids is K IC ∼ 1 MPa m 3/2 ∼ 10 8 erg cm −5/21 . The phase diagram slopes can have any values, however, typically most of them are found in the range of k ∼ 1 to 10 K/kbar ∼ (0.1 ÷ 1) × 10 −8 K cm 3 erg −1 [140]. The correlation radius represents one of the most accessible parameters, since it can be extracted from the width of the X-ray spectrum peaks [141], and typically exhibits a value of r c0 ∼ 1 to 10 nm K 1/2 . Making use of (60), one finds the typical temperature shift ΔT c ∼ 10 2 to 10 4 K .
It should be mentioned that the width ΔT c ≈ 300 K of the region of existence of the transformational PZ has been, indeed, observed during glass transition at the tip of a propagating crack in NiTi [142]. Other papers have not revealed the whole region of existence.
Since the typical values of the melting point of inorganic solids is T m ∼ 10 3 K, one concludes that ΔT c typically covers a considerable part, if not the whole, of the phase diagram above or below the line of a bulk phase transition.

On the difficulties and possibilities of detecting a transformation process zone outside of the hysteresis region
In the Introduction, we listed a number of materials for which PZ observation has been reported. A relatively small number of such materials other than those of the martensite type seem to contradict our main findings. This is, however, only an apparent contradiction. Detailed inspection of the papers cited above shows that in most of them, local phase transitions have been detected by the analysis of the fracture surface available after the sample has been broken, the so-called, "post mortem" examination. In the case of the zirconia, for example, the fracture surface exhibited a layer of the monoclinic daughter phase about 1 μm thick on top of the tetragonal mother phase surviving a considerable time after fracture. This requires the zirconia to be deep within the hysteresis region of its phase diagram. Indeed, martensitic transformations in both martensite-austenite metals and zirconia exhibit wide hysteresis regions. This is, however, rare. For most solids the hysteresis does not exceed ∼10 K, while the second order transitions have no hysteresis at all.
In contrast, the transformation zone outside of the hysteresis region is only present under stress. As soon as the stress is removed, it immediately disappears. It cannot, therefore, be detected by the "post mortem" inspection of the fracture surface.
Further, at a high temperature, one may observe no zone at K I = K IC , but since ΔT * ∼ K 4/3 I the zone may show up for a higher stress intensity factor, that is, at the tip of a propagating crack, rather than of a motionless one. The detection of a local phase transition outside the hysteresis region is, therefore, a challenging experimental task.

Summary
We formulated an approach to describe the emergence of a process zone at the tip of a crack. This was achieved by describing the process zone in terms of an order parameter, an internal degree of freedom responsible for the zone formation. It takes place as the response of the solid to a local stress inhomogeneity. The conclusions of the present paper are based on our following fundamental findings.
First, we have found the exact solution for the point of bifurcation both in the case of a motionless and of a propagating crack.
Second, we established the existence of a critical crack velocity controlling the disappearance of the local phase transition zone at high speeds.
Third, we found an analytical, asymptotically exact solution of the nonlinear equation describing the distribution of the order parameter close to the bifurcation point.
Fourth, by numerical simulations, we obtained the solutions away from the bifurcation point, in line with our analytic solution.

Appendix A: Simulation: Technical details
To perform the simulations, we used the software COM-SOL 4.3b. The equations have been simulated in the halfplane y ≥ 0. A semi-circular domain has been defined with the diameter D = 50. By trial and error we find that this is large enough to let the solution vanish well away from the domain boundary. The initial mesh size of 5 was chosen, but the adaptive mesh refinement option has been further used to automatically refine the mesh as appropriate. The no-flux boundary condition has been set at the boundary y = 0 and the condition u = 0 at the rest of its boundary. A straightforward simulation of the static equation, equation (25), with such boundary conditions, however, only returns the trivial solution u = 0 at any value of the control parameter q. To avoid this, instead of (25), we introduced a pseudo-dynamic equation:  (25), represent fixed points of the dynamic system (A.1). As initial condition we used a smoothed step function, only unequal to zero in a vicinity of the point (0, 0). The pseudo-time stepping approach converges away from the points q = 0 and from the bifurcation point, q = q c . In practice, we obtained a good convergence at q > 0.1, above the points of the global bifurcation. Closer to the point q = 0, we were unable to get an equilibrium solution. In the vicinity of the point q = q c , the method produced a small regular error, discussed in detail below. Apart from that, the method exhibited a good convergence, enabling us to study the distribution of the rescaled order parameter u(x 1 , y 1 ) in the vicinity of the crack tip.
The dynamic system has been solved using the direct MUMPS solver with the BDF time stepping. The convergence of the solution to its fixed point has been controlled by the behavior of u max , the maximum value of the function u(x 1 , y 1 , t ps ). By trials we found that 700 pseudo-time steps ensure a good convergence, though sometimes it was necessary to maintain the process for as long as 3000 steps. Figure A.1 shows an example of such a convergence study for a number of simulations in which all parameters except q were fixed, while q varied.
One can see that far from the bifurcation point, the convergence takes place well before 700 pseudo-time steps are completed. As can be expected, the situation is different in the close vicinity of the bifurcation (q = 0.38 and 0.39 corresponding to the vertex-down triangles and open circles in Fig. A.1). Even here, 700 pseudo-time steps guarantee a rather reliable convergence. where η 0 (r) (51) turns into zero. At smaller values of r one finds η = 0. For this reason, in I 2,3,4 one should only integrate over ρ from ρ 0 (θ) to infinity, while the integration over θ runs from −π to π.
The integration has been done numerically using the standard NIntegrate routine of Mathematica 10.1 [135] employing an even-odd subdivision method with a local adaptive strategy. Below, only the ratios I 2 I 4 /I 2 3 and I 3 /I 4 are used. Figure