Exactly solvable model for a velocity jump observed in crack propagation in viscoelastic solids

Needs to impart appropriate elasticity and high toughness to viscoelastic polymer materials are ubiquitous in industries such as concerning automobiles and medical devices. One of the major problems to overcome for toughening is catastrophic failure linked to a velocity jump, i.e., a sharp transition in the velocity of crack propagation occurred in a narrow range of the applied load. However, its physical origin has remained an enigma despite previous studies over 60 years. Here, we propose an exactly solvable model that exhibits the velocity jump incorporating linear viscoelasticity with a cutoff length for a continuum description. With the exact solution, we elucidate the physical origin of the velocity jump: it emerges from a dynamic glass transition in the vicinity of the propagating crack tip. We further quantify the velocity jump together with slow- and fast-velocity regimes of crack propagation, which would stimulate the development of tough polymer materials.

• A line crack propagates in the positive x-direction.
We consider a two-dimensional model on an N 1 × N 2 rectangular lattice where N 1 is (countably) infinite and N 2 is even. The lattice points are labeled by the index n = (n 1 , n 2 ) where n 1 ∈ Z ≡ {. . . , −1, 0, 1, . . . } and n 2 ∈ {1, 2, . . . , N 2 }. At ε = 0, the lattice constant in the x-and y-direction are a and b, respectively. In other words, under zero strain, the height and width of the sheet are aN 1 (= ∞) and bN 2 (= L), respectively. The position of the lattice point labeled by the index n is given by r 0 (n) = (an 1 , bn 2  FIG. S1: Ordinary and minimal lattice models for the fixed-grip crack propagation. a, Lattice model with a rectangular unit cell under zero strain. The lattice point is specified by two integers n = (n1, n2). The top and bottom boundaries are located at n2 = 1 and n2 = N2, respectively, and the two surfaces of the line crack are located at n2 = N2/2 and n2 = N2/2 + 1. In the box, magnified views on a lattice point and its nearest neighbor points with and without deformation are shown together with the indices of the points, the lattice constants a and b, and associated forces. b, Coarse-grained semilattice model obtained by removing lattice points except for the ones located at n2 = N2/2 and n2 = N2/2 + 1. c, Elementary deformation modes of a unit cell in the rectangular lattice with lattice spacings a and b.
In the model, every lattice point interacts with the nearest neighbor points (at most four points) by a minimal coupling. The tensile and shear stresses are given as F tensile /(ah) = E∆b/b and F shear /(bh) = G∆b/a, respectively, for the deformation characterized by ∆b (See, Fig. S1 for the case in which the unit cell deforms in the y-direction). Here, we have introduced Young's modulus E and the shear modulus G. The elastic energy of the sheet with the minimal coupling is then given by 2 , (I.1) where C 11 = C 22 = Ea/b and C 12 = C 21 = Gb/a. Here,1 ≡x ≡ (1, 0) and2 ≡ŷ ≡ (0, 1), δ µν is the Kronecker delta, and the summation extends over all the bonds on the sheet. We note that Poisson's ratio of this system is zero, i.e., G = E/2 because Eq. (I.1) does not contain the coupling between the displacements in the x-and y-directions. Moreover, assuming that the displacements in the x-direction u 1 (n) for arbitrary n are zero at the initial time, forces in the x-direction are always zero, and thus, every lattice point does not move in the x-direction. Then, in the following, we consider the motions and forces only in the y-direction.
To construct the equation of motion for the system, we consider tensile and shear forces, together with viscous forces. As illustrated in Fig. S1, the tensile and shear forces acting on the lattice point n in the y-direction are given, respectively, as follows: Here, we have introduced the notations, ∆ 2 ∆x 2 f (n) ≡ [f (n +x) − 2f (n) + f (n −x)] /a 2 and ∆ 2 ∆y 2 f (n) ≡ [f (n +ŷ) − 2f (n) + f (n −ŷ)] /b 2 , which correspond to the second-order partial derivatives in the continuum limits, a → 0 and b → 0, respectively. When a Kelvin-Voigt element is employed for the interaction in the y-direction, the following viscous term should be added: In this way, the equation of motion of the lattice point n in the y-direction is given by ] . (I.5) Assuming that the inertial term m ∂ 2 ∂t 2 u 2 (n) is negligible (the overdamped limit), we obtain the equation of motion per unit volume in the following form:

A. Construction of a minimal model
In this section, we construct a minimal model for crack propagation in viscoelastic sheets as shown in Fig. S1b, by ignoring all the lattice points except for the lattice points at n 2 = N 2 /2 and n 2 = N 2 /2+1. The original set of variables {u 2 (n)} is now represented by a much smaller set, For simplicity, we assume that the sheet is always symmetric about the x axis, and thus the lattice points on the upper side, u i ≡ u 2 (i, N 2 + 1) for i ∈ Z, completely describe the dynamics of the present model.
In the following discussion, it is important to distinguish three types of strain in the y-direction (see Fig. S2a; here and hereafter, we set l ≡ b and use l instead of b):   (ii) E i (t): the strain of the i-th upper (or lower) "long spring," i.e., the i-th spring of natural length (L−l)/2 directly connected to the top (or bottom) boundary (the strains of the upper and lower long springs are the same so that in the following we consider only the upper ones and "long springs" indicate the upper ones); N 2 )]/l = 2u i /l: the strain of the i-th "short spring," i.e., the i-th spring of natural length l located at the center in the y-direction (provided that the i-th short spring exists).
Here, ε is a constant, but E i (t) and E (l) i (t) depend on time t. We mainly use E i (t) to describe the dynamics of crack propagation because E (l) i (t) is expressed by E i (t) by the following relation (see Fig. S2b): i (t) provided that the i-th short spring exists because "the spring constant" of long springs is much smaller than that of the short springs (the corresponding Young's moduli are the same for the short and long springs). To realize crack propagation, we remove the i-th short spring when E (l) i (t) ≥ ε c . At places far from the crack tip, E i (t) approaches a constant value and the value depends on whether the i-th long spring is located on the front (i.e., right) or rear (i.e., left) side of the tip of a crack propagating in the positive x-direction (see Fig. S1b): We construct the equation of motion in the y-direction for the i-th lattice point characterized by u i (t) in the minimal model. First, we consider tensile stresses on the basis of Eq. (I.6); as illustrated in Fig. S2c, considering the elongational deformation in the y-direction of the i-th long and short springs, we obtain the tensile stresses acting on the i-th lattice point in the following form: Here, we explicitly show the subscript 0 for the Young's modulus, E 0 , whereas we have omitted the subscript in Sec. I, i.e., E ≡ E 0 . This is because we should distinguish the two springs E 0 and E 1 in a Zener element, which is a generalization of a Kelvin-Voigt element and will be considered in Secs. III and IV. Second, we consider shear stresses; as illustrated in Fig. S2c, considering the shear deformation in the y-direction, we obtain F right,i − F left,i = alhµ ∆ 2 ∆x 2 u i (t). Note that shear modulus G for the ordinary lattice model (Fig. S1a, discussed in Sec. I) is replaced by the "effective" shear modulus µ for the minimal model (Fig. S1b). This is because the shear modulus for the minimal model is considered to effectively represent all the forces acting on the decimated points from the nearest neighbor points located in the x-directions, with a spirit similar to the one employed in renormalization [S1] in statistical physics. Third, combining the tensile and shear stresses together with the viscous terms, we obtain the equations of motion in the y-direction for the i-th lattice point of mass m, by noting that F down,i (t) and F (η) down,i (t) are missing on the rear side of the crack tip: , we can rewrite the above equations of motion in terms of E i (t) by removing the dynamic variable u i (t). Finally, assuming that the inertial term is negligible (the overdamped limit), we obtain the equations of motion for the field E i (t): We here make three comments on the equations of motion (II.5). (i) When we neglect shear force (µ = 0), the two expressions on the front and rear sides reduce to the same equation of motion (except for the equilibrium position), which describes the Kelvin-Voigt dynamics. (ii) The "second-derivative" term E 0 ∆ 2 ∆y 2 u 2 (n) in Eq. (I.6) has been replaced by the strain-field term, The net strains of long springs on the rear and front sides are different and given by E i (t) and E i (t) − ε, respectively, because lim i→−∞ E i (t) = 0 and lim i→−∞ E i (t) = ε on the rear and front sides, respectively (see Eq. (II.2)).
To derive the equations to be solved analytically, we take the continuum limit of Eq. (II.5) in the x-direction, a → 0. In this limit, the finite difference ∆ 2 ∆x 2 is replaced by the derivative ∂ 2 ∂x 2 and the discrete strain field E i (t) by the continuum strain field E ≡ E(τ, χ) where we have introduced the dimensionless parameters, In this way, the above discrete version of the equations motion (II.5) is replaced by the following continuum equations:

B. Derivation of an exact solution of crack propagation with a constant velocity
In this subsection, we solve Eq. (II.7) for the static case (V = 0) and for the dynamic case in which a crack propagates with a constant velocity (0 < V < ∞), by seeking a solution of the form E(τ, χ) = f (χ − ντ ). Here, we have introduced the dimensionless velocity of crack propagation, Substituting f (χ − ντ ) into Eq. (II.7), we have second-order linear ordinary differential equations: Here, since the width (in the x-direction) of the sheet is large enough, no generality is lost by setting the position of crack tip to χ = 0 with τ = 0: the crack exists in the region χ < 0 and is absent in the region χ ≥ 0. We give the boundary conditions for the differential equations (II.11) as follows: the conditions at remote edges, and the matching conditions at the crack tip for the strain field E, As previously mentioned, we remove the short spring (located in the center in the y-direction with natural length l, see Fig. S2) when E (l) (t) ≥ ε c . We rewrite this inequality as We solve the ordinary differential equation (II.11) under the boundary conditions (II.12) and (II.13). First, we solve the differential equation for χ < 0. Substituting the form f (χ) = Ce −χ/ξ into Eq. (II.11), we obtain the characteristic equation (quadratic equation for ξ), (II.15) and the solutions, From the boundary conditions (II.12) and (II.13), only the solution ξ N,− is relevant: Second, we solve the differential equation for χ ≥ 0. Substituting the form f (χ) − ε = Ce −χ/ξ into Eq. (II.11), we obtain the characteristic equation, g 1 (ξ) = ξ 2 + νξ − 1 = 0, and the solutions, . From the boundary conditions (II.12) and (II.13), only the solution ξ 1,+ is relevant: In the following, we determine f (χ) from Eq. (II.17) in the cases of ν = 0 and 0 < ν < ∞.
In the static case (ν = 0), we can rewrite Eq. (II.17) as f (0) N and ξ 1,+ = 1 from Eq. (II.16). This gives the following solution: Note that Eq. (II.19) implies thatε > 1/ √ N in the case of crack propagation, ν > 0. In the case of crack propagation with a constant velocity ν (0 < ν < ∞), since f (0) = f c , the exact relationship between ν andε is obtained from Eq. (II.17) as The expression (II.20) does not exhibit the velocity jump in the physically relevant range ofε, 1/ √ N <ε < 1. By using Eq. (II.20), we obtain the dynamics of the strain distribution: Here, two healing lengths have been introduced (see Fig. S4c below): Equation (II.20) or (II.21) gives the relation between the velocity V and the initially applied energy density w = 1 2 E 0 ε 2 , by virtue of the relation

A. Generalization of viscoelastic interaction: equations of motion in the y-direction
In this section, we generalize the present model by changing the interaction from the one based on Kelvin-Voigt elements to the one on Zener elements (see Fig. S3). A Zener element is a parallel connection of the spring characterized by E 0 (strain E(t)) and a Maxwell element, which is a serial connection of the spring characterized by E 1 (strain E 1 (t)) and the dashpot characterized by η (strain E 2 (t)). The strains of the parallel components are identical: (III.1) For the serial components in the Maxwell element, the stresses on the spring and dashpot are identical: Note here that a Zener element reduces to a Kelvin-Voigt element in the limit E 1 → ∞, in which E 1 (t) = 0. The generalization of the equations of motion in the y-direction is achieved by the following replacements: With this replacement, Eq. (II.5) in the continuum limit in the x-direction, a → 0, is given as As in the case of the Kelvin-Voigt interaction, we use the dimensionless parameters, τ = t/t 0 , χ = x/x 0 , and N ≡ L/l, where t 0 ≡ η/E 0 , and x 0 ≡ l √ (L − l)µ/(2LE 0 ). Then, we rewrite the above equations of motion as whereĖ 2 ≡ ∂ ∂τ E 2 (τ, χ) and E ′′ ≡ ∂ 2 ∂χ 2 E(τ, χ). We eliminate E in Eq. (III.5) by using the following relation obtained by substituting Eq. (III.2) into Eq. (III.1), to have the equations of motion for a single field E 2 in the following form: Here, we have introduced the parameter instead of s. Although the physical meaning of λ is clearer than that of s, we mainly use s in this section for mathematical simplicity. The parameters s and λ vary in the ranges of 0 < s < ∞ and 1 < λ < ∞, respectively. Note that the Kelvin-Voigt interaction (finite E 0 and E 1 = ∞) corresponds to the limit s → 0 or equivalently λ → ∞.

B. Derivation of the exact solution of crack propagation with a constant velocity
In this subsection, we solve Eq. (III.7) in the case of crack propagation with a constant velocity (0 < V < ∞). Note that, in the case of V = 0, a Zener element reduces to a Kelvin-Voigt element, and we have already solved the model consisting of Kelvin-Voigt elements in Sec. II.
We derive the equations of motion with relevant boundary conditions in the present case, following the manner employed in the Kelvin-Voigt case in Sec. II. We substitute the form E 2 (τ, χ) = f (χ − ντ ) into Eq. (III.7) and derive linear ordinary differential equations for f (χ) for χ < 0 and 0 ≤ χ, where we set the position of the crack tip to χ = 0 with τ = 0. Here, V and ν = V /V 0 are the dimensional and dimensionless velocities, respectively, with V 0 ≡ x 0 /t 0 ≃ lE 0 /η as before (see Eq. (II.10)). The result is given by for 0 ≤ χ (front side). (III.10) To determine the boundary conditions, we substitute E 2 (τ, For χ − ντ = 0 which corresponds to the crack tip, we assume that the strain distributions E and E 2 are continuous and differentiable at the crack tip. The values of the functions and their derivatives match at the crack tip: +0) from the others). In addition, as in the case of the Kelvin-Voigt interaction, we assume that the short spring is absent when E (d) ≥ ε c , which yields E(−0) = E(+0) = N ε−εc N −1 . In summary, we give the appropriate boundary conditions as (III.11) We can then establish the exact analytical relation between the initially applied strain ε and the normalized velocity of crack propagation ν ≡ V /V 0 as in the following theorem:
Lemma 1. Let s, ν, and N are positive real numbers. Then, g N (ξ) = 0 has one positive and two negative real solutions for −∞ < ξ < ∞.
We prove Theorem 1 with the aid of Lemma 1 (we give the proof of Lemma 1 at the end of this subsection).
Proof of Theorem 1. First, we solve the differential equation for χ < 0 (the first equation of Eqs. (III.10)). Substituting the form f (χ) = Ce −χ/ξ into the equation, we obtain the characteristic equation (III.14). According to Lemma 1, we can set the three solutions, ξ N,1 , ξ N,2 , and ξ N , of the cubic equation g N (ξ) = 0 to satisfy the relation ξ N,1 < ξ N,2 < 0 < ξ N . To satisfy the boundary conditions at ξ → −∞, the solution of the differential equation for χ < 0 has the form f (χ) = ∑ 2 i=1 C i e −χ/ξ N,i , where C 1 and C 2 will be determined later by using the boundary conditions at the crack tip, i.e., χ = 0.

IV. MINIMAL MODEL CONSISTING OF ZENER ELEMENTS 2: LOW-AND HIGH-VELOCITY REGIMES AND VELOCITY JUMP
In this section, we use Eq. (III.12) to investigate the dependences of the initially applied energy density w on the velocity ν ≡ V /V 0 in low-and high-velocity regimes, and derive the existence condition of the velocity jump. Note that Eq. (III.12) gives the relation between the initially applied strain ε and ν and that ε is related to w simply as w =ε 2 N w 0 (see Eq. (II.24)). In the following, the present model will be analyzed for arbitrary positive real number s ≡ E 0 /E 1 > 0, i.e., for λ ≡ E ∞ /E 0 = 1 + 1/s > 1. Note that the relations derived below are further simplified for elastomers, for which the relation 1 ≪ λ ≪ N is valid (typically λ ≃ 10 2 -10 3 and N ≃ 10 6 -10 9 ).

Theorem 2 (Asymptotic forms in low-and high-velocity regimes).
If λ > 1 and N > 1, thenε Comparing the first and second terms on the right-hand side of Eq. (IV.5), we can estimate the range of ν in which is satisfied, then w(ν) ≃ w 0 (Here, we have omitted the factor 2 in the second term on the right-hand side of Eq. (IV.5)). When √ N ≫ 1, the condition (IV.8) is simplified as Similarly, Eq. (IV.6) enables us to evaluate the w-V relation in the high-velocity regime. From the first term on the right-hand side of Eq. (IV.6), we have lim ν→∞ε (ν) = λ/( √ N + λ − 1), or equivalently We can estimate the range of ν in which w(ν) ≃ λ 2 N w 0 /( √ N + λ − 1) 2 holds, following the manner we employed in the low-velocity case. However, we should pay attention to the fact that Eq. (IV.6) is insufficient for the estimation in the case of (unrealistic) large λ. In fact, for example, in the Kelvin-Voigt limit (λ → ∞) the second term on the right-hand side of Eq. (IV.6) goes to zero. Thus, in this discussion, we assume λ ≲ N , which includes the case of typical viscoelastic materials (1 ≪ λ ≪ N ). Comparing the first and second terms on the right-hand side of Eq. (IV.6), we have the range of ν in which w(ν) ≃ λ 2 N w 0 /( √ N + λ − 1) 2 holds: if the crack propagation velocity is sufficiently high, i.e., In the case of √ N ≫ 1 and λ ≫ 1, which includes the case of typical viscoelastic materials, the above statement can be simplified, i.e., if is satisfied, thenε(ν) ≃ λ √ N +λ , or equivalently We summarize characteristic scales derived from Lemma 2 and Theorem 2 for 1 ≪ λ ≪ N in Fig. S4, together with the scales associated with the jump derived in the next subsection.

B. Existence condition of the velocity jump
As illustrated in Fig. S4a, Eq. (III.12) guarantees the existence of the velocity jump under an appropriate condition. In this subsection, we rigorously derive the existence condition of the velocity jump by analyzing Eq. (III.12). Before we give the rigorous derivation, we roughly derive the existence condition of the velocity jump from Theorem 1 with the aid of Lemma 2.
We evaluate the positive solution of the characteristic equation (III.14) under the condition (λ−1)/ √ λ ≪ ν ≪ √ N . According to Lemma 2, the positive solution of the characteristic equation g N (ξ) = 0 is expressed as ξ ≃ ξ Similarly, the positive solution of g 1 (ξ) = 0 is expressed as ξ ≃ ξ (IV.14) By noting the relation 1 for the numerator in Eq. (IV.14) and the relation for the denominator in Eq. (IV.14), we obtaiñ (IV.16) Therefore, we obtain the following rough estimation for the velocity jump.
which is the position of the velocity jump shown in Fig. S4a.
We have two remarks for Theorem 3: (i) the present model consisting of Kelvin-Voigt elements (which corresponds to λ → ∞) never satisfies the existence condition of velocity jump in Eq. (IV.33) because N is finite; (ii) if λ ≫ 1, then the existence condition (IV.33) reduces to √ λ ≪ ν ≪ √ N .
where we have used We evaluate the lower boundε lower as follows.
where we have introduced a positive number α: .