{'fraction_non_alphanumeric': 0.04714133522727273, 'fraction_numerical': 0.013405539772727272, 'mean_word_length': 3.9651777843079636, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Steady-State Cracks in Viscoelastic Lattice Models II

We present the analytic solution of the Mode III steady-state crack in a square lattice with piecewise linear springs and Kelvin viscosity. We show how the results simplify in the limit of large width. We relate our results to a model where the continuum limit is taken only along the crack direction. We present results for small velocity, and for large viscosity, and discuss the structure of the critical bifurcation for small velocity. We compute the size of the process zone wherein standard continuum elasticity theory breaks down.


II. DESCRIPTION OF THE MODELS
The lattice model we study is identical to that described in our earlier work, [11]. We have a square lattice of mass points undergoing (scalar) displacement out of the plane. The lattice extends infinitely long in the x-direction, with N + 1 rows in the y-direction. The lattice points are connected by linear "springs", with spring constant 1, to their nearest neighbors. The top row is displaced a fixed amount ∆. The bottom row is connected to a fixed line, with piece-wise linear springs. These springs, with spring constant k, "crack" irreversibly if they are stretched an amount ǫ. When k = 2, this model is equivalent to a system of 2N + 2 rows, loaded by ±∆ from top and bottom, with a symmetric crack running down the middle, with extension at cracking of the springs that bridge the middle being 2ǫ. All the (uncracked) springs have a viscous damping η. The equation of motion for the system is then for j = 1 with u i,N +1 ≡ ∆, and Note that in these units, the elastic wave speed is unity, so all velocities are dimensionless, expressed as fractions of the wave speed.
We are interested in steady-state cracks, described by the Slepyan traveling wave ansatz, which implies that every mass point in a given row undergoes the same time history, translated in time. We choose the origin at time such that u 1 (0) = ǫ so that it represents the moment of cracking of the spring attached to the bottom row mass point. The equation of motion is best expressed in terms of the N × N coupling matrix The steady-state equation then reads We will also consider in this paper an x-continuum version of this model, where we replace the nonlocal in time coupling along the crack with its continuum analog

III. SURVEY OF RESULTS
In this section, we survey the major results derived in the bulk of the paper. As the derivations are exceedingly technical, it is useful to present the results first by themselves so that they may be appreciated without getting lost in a welter of technical complications.
We begin by completing the Wiener-Hopf (WH) solution of the continuous x, discrete y model, as the results are simpler and are a useful basis for assimilating the more complicated results of the full lattice model. The key aspect of the solution is the calculation of ∆ as a function of the crack velocity v (in units where the wave speed is unity). We find ∆ ∆ G = √ kN + 1 m q 1,m (1 + ηvQ 1,m ) Q 1,m (1 + ηvq 1,m ) (7) which expresses ∆ (normalized to the Griffith value ∆ G = ǫ √ 2N + 1 (8) at which the uncracked state becomes metastable) in terms of the wave vectors corresponding to the various normal modes of the problem. If we label the normal mode eigenvalues of the y-coupling matrix on the uncracked side M(k) by Λ m , then Q 1,m is the unique positive root of the dispersion relation Similarly, q 1,m is the unique positive root of the dispersion relation using the normal mode eigenvalue λ m of the cracked side M(0). This formula is fairly complicated, but simplifies tremendously for the case of symmetric cracks (k = 2) in the macroscopic limit N >> 1. Then, the product above can be performed analytically, with the simple result where Q 1 (1) is the mode associated with the highest frequency y-mode, with Λ = −4. For typical η's of order 1, Q 1 (1) does not vary much from its zero velocity value of 2. The resulting curve ∆(v)/∆ G starts linearly at v = 0 from 1 with slope η and diverges at v = 1. Thus, in the infinite N limit, the velocity never exceeds the wave-speed. At any finite N , however, the velocity crosses the wave speed at a ∆ of order N 1/6 ∆ G . Since the divergence with N is so weak, crossing the wave-speed barrier may not be as difficult as one would naively think. This is especially true for small dissipation, where the critical ∆ scales as (ηN ) 1/6 . This appears a more likely mechanism for explaining the experimental observation of supersonic cracks than the time-dependent forcing hypothesis of Slepyan. [16] The basic structure is unchanged when we go over to the full lattice model. The essential difference is that the lattice dispersion relation is nonpolynomial and has an infinite number of positive (real-part) solutions for each eigenmode m. The ∆ − v relationship is where now the product extends over all positive real-part roots Q 1,n,m of the lattice dispersion relation for each Λ m (λ m in the case of q 1,n,m ). For a given m, there is one real positive root, Q 1,0,m (q 1,0,m ), and an infinite series of complex-conjugate pairs of complex roots, ordered by increasing imaginary part. For large n, the imaginary part increases by roughly 2π for each successive root. Again, for symmetric cracks we can evaluate analytically the macroscopic (large N ) limit. We obtain where q ∞,n is the root corresponding to the highest frequency Λ = −4 eigenmode and plays the role of Q 1 (1) of the previous x-continuum result. The q 0,n are the roots corresponding to the Λ = 0 eigenmode. These do not have a counterpart in the x-continuum calculation as the n = 0 real solution vanishes, and only the lattice-induced n = 0 modes enter. As indicated by the way we expressed this result, we can consider it as essentially the x-continuum result, Eq. (11), with the real lattice q ∞,0 replacing Q 1 (1), modified by a multiplicative correction factor involving the complex lattice modes. To understand the usefulness of this way of thinking, as well as its limitations, we present in Fig. 1, for η = .5, the exact numerically computed relationship Eq. (13), along with the x-continuum result Eq. (11). In addition, we plot the lattice result, truncated after its n = 0, |n| = 1 and |n| = 5 terms. We see that at larger velocities all these results are close, indicating that the lattice-induced shift in q as well as the additional lattice modes play little role at these velocities. At smaller velocities, the various approximations differ significantly from each other and from the exact curve. We see, in fact, that as v approaches 0, more and more terms must be included in the product to achieve an accurate result. The calculation of the limiting behavior at small velocities requires summing all the terms. The result of the calculation is that for all η, as v → 0, ∆ approaches ∆| 0 + = 1 + √ 2∆ G , the maximal ∆ for which an arrested crack exists. This generalizes the result of Slepyan for infinitesimal dissipation. As v increases, ∆ decreases linearly with the η-independent slope, −∆| 0 + /2 so that the bifurcation from the arrested is subcritical and universal.
More progress can be made in the large η limit. Here, at fixed ∆, the velocity goes to zero as η increases, so that the ratio φ ≡ ηv is fixed. In this limit, we can calculate the infinite product and find This infinite-η result, together with the exact result for various η's, is presented in Fig. 2. We see that this calculation does not reproduce the subcritical bifurcation from the arrested crack at small velocities, which is a higher order effect. We can evaluate this 1/η correction near the bifurcation at small φ, and find which reproduces the small v behavior described above and shows that the η dependent corrections are in fact exponentially small in v. The resulting ∆ − v curve starts at ∆| 0 + at v = 0, heads back linearly for a short distance of order 1/η(ln η) 2 and then sharply veers forward. A last result worth noting is that whereas the Kelvin viscosity model analyzed herein has a nice macroscopic limit when expressed in terms of ∆ G , the model with Stokes viscosity, where the dissipation in put in the masses and not in the bonds, does not have such a limit. There an O(1) Stokes viscosity at the microscopic level changes the continuum elastic fields and requires an ever-increasing ∆/∆ G as the sample is made wider.

IV. THE x-CONTINUUM MODEL
We begin our analysis with the solution of the x-continuum model, Eq. (6), introduced in Kessler, et. al. [11]. It is important to remember that is this model, the lattice structure in the y-direction is left unchanged. The solution of the lattice model is similar in structure to that of the x-continuum model, but the latter is a simpler context in which to develop the necessary techniques. Furthermore, the x-continuum model is an interesting approximation in its own right, which captures a significant amount of the structure of the full lattice problem. In Kessler, et. al. [11], a Wiener-Hopf analysis of the problem was initiated. In this analysis, the key technique is to decompose all the terms in the steady-state equation of motion into terms analytic in the upper-and lower-half planes respectively. However, the analysis was not carried to completion, due to the presence of one term whose decomposition was not evident. Here we use a trick to accomplish the decomposition of this last remaining term, and thereby complete the solution of the problem. We choose not to reproduce the lengthy preliminary stages of this calculation, for which the interested reader is referred to [11]. We do however reiterate the definition of the relevant notations introduced there, so that the current exposition is minimally self-contained.
The problematic term, from [11] Eq. (42), is where the χ, q, and Q are the roots of a certain family of cubic polynomials. In detail, let ℓ l , (l = 1, . . . , N − 1) be the eigenvalues of the (N − 1) × (N − 1) coupling matrix M N −1 (1) defined in Eq. (4) above. Define the polynomial specifying the dispersion relation, P (λ, Q), by Then, P (ℓ m , Q) has, for each m, three roots, one positive which we denote χ 1,m , and two with negative real parts, which we denote by −χ 2,m , −χ 3,m , so that all the χ's have positive real parts. Similarly, denote the eigenvalues of the It is apparent from these brief remarks that the troublesome term in its current state involves singularities and poles in both the upper-and lower-half planes. To proceed, we rewrite the numerator using the following manipulations: The first line of this chain employed an identity from [11], Eq. (40), relating the numerator to the determinant of a certain matrix formed from M N −1 (1) and the identity matrix I together with the function The second line in this chain, claiming this determinant is equivalent, up to a constant factor, to the difference of two N × N determinants can be proven by expanding each of the matrices about the first row. The last line reexpresses each of these two determinants using more identities from [11], Eqs. (38-39). After these manipulations, our term can be written As the outside factor has a pole at −i/ηv in the lower-half plane, the second term only has singularities and poles in the lower-half plane and so is in the desired form. The first term is still mixed and requires further massaging. The idea is to subtract out the unique lower-half plane pole so that what is left has only upper-half plane poles and zeros. Thus, where now g − has only upper-half plane poles and zeros. We will not need the explicit form of g − in the calculation. What we have is now sufficient to solve forũ + , the Fourier transform of the displacement of the bottom masses in the crack region, u 1 (x)θ(x). Using [11], Eq. (42), we find Solving forũ + , we find Fourier transforming and evaluating at x = 0 + yields so that Using u 1 (0) = ǫ, ∆ G = ǫ √ kN + 1 and the relations (see [11] Eq. (43)) we obtain our desired result The primary benefit of this method of solution over the direct approach employed in [11] is that for the symmetric crack (k = 2) we can take the large-N limit. To do this, we break up the N −fold product into two terms. The first is We transform the product into the exponential of a sum over logarithms, a sum which for N large we can approximate by an integral, via the Euler-MacLauren Summation Formula (EMSF) [14]. Thus Now, for k = 2, Λ m = −4 sin 2 ( πm 2N +1 ) and λ m = −4 sin 2 ( π(m−1/2) 2N +1 ). If we define α = m/N , then we see that The integral is now a total derivative, and so When α = 1, m = N and so Λ m ≈ −4, and so Q 1 (1) satisfies Similarly, when α approaches 0, so does m and so also Λ m . This in turn implies that Q 1 (0) = 0. So, finally, The second factor is slightly more difficult to treat, since the numerator and denominator both vanish as m → 0. To handle this, we regularize the product by multiplying and dividing by m λ m /Λ m , which we can perform analytically. Then, the regularized product, Π R 2 , can be transformed to the exponential of an integral of a total derivative, which can be calculated explicitly. In detail, where we have used the fact that for α small, Also, so putting all the pieces together yields the simple result which expresses ∆ in terms of Q 1 (1), the wave-vector at the end of the Brillioun zone.
The most striking lesson of this formula is that ∆ diverges at v = 1, the wave speed. Thus, while at any finite N , there is no upper limit to the velocity, at infinite N the wave speed is an absolute upper bound to the crack velocity. At large ∆, v approaches unity from below as 1/∆ 4 . A second lesson is that at small velocity, ∆ ∼ ∆ G (1 + ηv), so that ∆ approaches ∆ G linearly, as is generally true for this x-continuum model. A third implication is the behavior at large η. For fixed ∆, v decreases as η gets large, so that Q 1 (1) satisfies In this large η limit, of course, ∆ is a function of the scaling variable ηv, which was first introduced in [15]. We have seen how at infinite N , the crack speed v never crosses unity, the wave speed. However, at any finite N , there is a ∆ for which the crack speed crosses unity, which must diverge with N . We now calculate how this threshold scales with N . The key to the calculation is Π 2 , since it is the vanishing of Π 2 which leads to the divergence of ∆ at v = 1 for N infinite. To compute the value of Π 2 at v = 1 for finite large N , we need to choose a different regularization. We now define so that Now, since Q 1 (α)/(−Λ(α)) 1/3 approaches the finite limit 1/η 1/3 as α goes to 0, Π R 2 has a finite limit at v = 1 as N goes to infinity, namely Π R 2 ≈ Q 1 (1)(η/4) 1/3 . Using the infinite N limit of Π 1 from Eq. (33), we get Thus, the threshold ∆ scales as N 1/6 ∆ G , in accord with the numerical evidence discussed in [11]. The coefficient goes to 2 for η large, and vanishes as η 1/6 for small η.
Another manifestation of this same phenomenon, the disappearance of the v = 1 crossing in the infinite-N limit, is the nonuniformity of the the large-N limit as v approaches 1. Working out the corrections to the EMSF, we find that and Thus, the relative error of the infinite-N approximation is O(1/N ), and diverges as v approaches 1 as (1 − v) −3/2 . It is also interesting to note that the relative error vanishes as v goes to zero, son that the infinite-N approximation becomes better at small velocities.
One last interesting piece of information we can derive from our solution is the size of the "process zone", the region where the solution from continuum elastic theory breaks down. The leading-order macroscopic solution was derived in [11], and exhibited the classic square-root singularity at the crack tip, x = 0. This singularity in really present only at infinite N , and is cut off by the upper limit on the Q's, (relative to the smallest Q ∼ 1/N ) at finite N . We can determine the structure of the process zone which replaces the singularity by studying our exact solution forũ + , Eq. (23), for K's of order 1. Using the EMSF to evalute the infinite product, similar to the derivations above, we find Using Q 1 (0) = Q 2 (0) = 0, Q 3 (0) = (1 − v 2 )/ηv, and the result for ∆, Eq. (38), we get Examining this expression for small K, we find the expected K −3/2 singularity, which gives rise to the square-root singularity of the outer solution. The coefficient of the K −3/2 singularity to leading order in N is ∆i 3/2 N −1/2 (1 − v 2 ) −1/4 , which reproduces the same η independent coefficient of the square-root singularity, or equivalently stressintensity factor, found in [11]. The structure of the process zone is governed by the other singularities inũ + that lie off the origin. In particular, the size of the process zone is determined by the singularity nearest the real line. For small η, this is at K = −iQ 2 (1) ≈ −2i/ √ 1 − v 2 , so the process zone is truly microscopic, unless the velocity is very close to 1. For large η, the dominant singularity is at K = −iQ 3 (1) ≈ −i/ηv, so the process zone grows linearly with η in size.

V. THE LATTICE MODEL
In this section, we generalize our solution of the continuum model to the lattice model. For ease of presentation, we will present the derivation only in the N = 1 case. The case of general N follows in a straightforward manner from this derivation and that of the continuum finite N model presented in the previous section.
Our derivation follows directly along the lines of our WH treatment of the continuum N = 1 problem in [11]. The equation of motion of the steady-state crack is Upon Fourier transforming, we find whereũ is the Fourier transform of u andũ ± are the transforms of θ(±t)u(t). We define the function in terms of which This function, R(λ; Q), which is the lattice equivalent of the polynomial P employed in the previous section, has not 3 roots, but in fact an infinite set of zeros in the complex plane. We shall label these zeros according to their real parts, Q 1,n (q 1,n ) are the zeros of R(−(1 + k); Q) (R(−1; Q)) with positive real parts, and Q 2,n ′ (q 2,n ′ ) are their counterparts with negative real parts. The indices n, n ′ run over the entire infinite set of zeros but are otherwise left unspecified for now. We can decompose R in terms of its zeros Using this, we rewrite the equation of motion: As in the last section, the hard part is to decompose the last term. The trick is the same, rewriting the numerator as the difference of R's.
As before, the second term is now fine, but the first term is still mixed. Again we subtract out the unique pole in the lower-half plane which is what we need to findũ + .
where g − only has poles and zeros in the upper-half plane. Separating out the pieces analytic in the upper-half plane yields Solving forũ + yields Fourier transforming and evaluating at x = 0 + , we find so that ∆ = u 1 (0)(1 + k) n q 1,n (1 + ηvQ 1,n ) Q 1,n (1 + ηvq 1,n ) As ∆ G = u 1 (0) √ 1 + k, we obtain our desired result As there is exactly one real positive root of R(λ; Q), it is convenient to assign this the index 0 and to label the complex roots in order of imaginary part, so that for example Q 1,n and Q 1,−n are complex conjugates. It is clear the basic structure of the lattice result is similar to the continuum result Eq. (27) above, with the continuum Q 1,1 , q 1,1 replaced by their lattice counterparts Q 1,0 , q 1,0 , and multiplied by a correction factor due to the additional infinite hierarchy of complex Q, q's which solve the lattice dispersion relation. The generalization to finite N is straightforward and is left as an exercise to the reader. The result is the direct generalization of the N = 1 result. At finite N , there is a set of zeros with positive real part of R(Λ m ; Q), (R(λ m ; Q)), for each m = 1, . . . , N , now labeled Q 1,n,m (q 1,n,m ). Then is the solution to the lattice problem at finite N . It of course reduces in the limit η → 0 + to the result of Marder and Gross [9]. As in the continuum, this rather unwieldy formula simplifies tremendously in the symmetric crack case k = 2 as N goes to infinity. The procedure for evaluating the limit is similar to the continuum calculation and so we do not present the details. What enters again are Q 1,n (α), at the two extremes of the Brillouin zone α = 0, 1. If we label Q 1,n (1) = q ∞,n , Q 1,n (0) = q 0,n then they satisfy the dispersion relations In terms of these q's, the infinite N limit solution is Again, this is very essentially similar to its continuum counterpart, with the real lattice wave-vector q ∞,0 playing the role of the continuum wave vector Q 1 (1), and with a multiplicative correction due to the presence of complex lattice wave-vectors. It should also be noted that this result reduces to that of Slepyan [6] in the η → 0 + limit.

VI. THE SMALL VELOCITY LIMIT
We begin our explorations of the content of our key result, Eq. (64) by examining the η fixed, v → 0 + limit. It is not sufficient to simply set v = 0, since as v gets smaller, more and more terms contribute significantly to the infinite product, as seen in Fig. 1. The proper treatment is to replace the infinite product by an infinite sum of logarithms and then approximate the infinite sum by an integral via the EMSF. Note that for v = 0, q ∞,n satisfies sinh 2 q∞,n 2 = 1, with the solution q v=0 ∞,n = 2πin + ω, where ω is the unique real root of the equation, namely ω = 2 ln(1 + √ 2). Similarly, q v=0 0,n = 2πin.
As we discussed above, we need to consider v → 0, n → ∞, α ≡ 2πηvn fixed. Then, writing q ∞,n ≡ 2πin + ω ∞ , ω ∞ satisfies Similarly, We can now easily approximate the first infinite product, yielding The second product is somewhat trickier, because a naive expansion diverges at small α. We define a regularized product so that (using the product formula for sinh, and the fact that sinh(ω/2) = 1) Our regularized product is now easily approximated, Using the identity we obtain our desired result This result is, as desired, explicitly independent of v, but would appear to depend on η through the very nontrivial η dependence of ω ∞ and ω 0 under the integral. It is possible to explicitly evaluate the integral for small and large η. For large η, ω ∞ − ω ∼ O(1/η 2 ) and ω 0 ∼ O(1/η), so the integral vanishes and so For small η [6], ω ∞ − ω 0 is concentrated at small α ∼ O(η), so it is appropriate to convert the integral into a principal value integral and do the ω integral immediately. In the remaining integral, we change variables to β = α/2η. Then, the denominator in the integrand reduces to 1/β, so only the odd (i.e. imaginary) part of ω ∞ − ω 0 contributes. For β ≥ 0 we find and The integral thus becomes So, again the integral vanishes, and the result for large and small η is the same. One is lead to guess that in fact the integral vanishes for all η, as indeed a numerical computation confirms. This is physically reasonable, since ∆| 0 + should be nothing other than the maximal ∆ for an arrested crack, which was previously found numerically [11] to be approximately 1.55∆ G . This maximal arrested crack ∆ is the result of a static calculation, and is of course completely independent of η. The vanishing of the integral can be demonstrated analytically and is the result of the fact that the integral has no singularities in the lower-half-plane. One can then close the contour there and the result is identically zero.
To see this, one has to study the analtyic structure of the functions ω ∞ (α), ω 0 (α). Consider first ω 0 (α) = 2 sinh −1 (y(α)), where y 2 (α) ≡ −α 2 /(4η 2 (1 + iα)). Since sinh −1 (y) = 2 ln( (y 2 ) + 1 + y 2 , ω 0 has a branch cut singularity along the line y 2 = −r where 1 > r > 0. Working out the algebra, in the complex α plane this works out to be, for η > 1, two separate curves. The first is a segment along the upper imaginary axis from α = 2iη(η − η 2 − 1) up to α = 2iη(η + η 2 − 1). The second is the circle of radius 1 centered at the point η = i. Similarly, ω ∞ has branch cuts for 2 > r > 1, which are two finite segments along the positive imaginary axis extending above and below the ω 0 branch cuts. For √ 2/2η < 1, the branch cut for ω 0 is a sector of the circle, while the branch cut for ω ∞ is the rest of the circle and a finite piece of the the entire imaginary axis extending centered about 2i. For η < √ 2/2, the branch cuts are confined entirely to a part of the circle. Thus, the singularities for all η lie entirely in the upper-half plane and, as advertised, the integrand is analytic in the lower-half-plane and so the integral vanishes.
The next step is to extend this calculation to next order in v. There are two sources for this first correction in v. One comes from the higher-order velocity dependence of q ∞,n , q 0,n . The other comes from the EMSF correction to the replacement of the infinite sum by an integral. The calculation of the first piece is similar in structure to the leading order calculation, just more involved. We expand q s,n ≈ 2πin + ω s + vσ s , where the subscript s = ∞, 0, and find It is straightforward, though tedious, to substitute this in Π 1 , Π R 2 and expand, giving a multiplicative correction factor of The integrand is again a very nontrivial function of η and α, but again a miracle occurs and the integral vanishes identically (as seen by numerical computation) for all η! The same analytic argument as above can be used to prove this point. This leaves us with only the second source for a O(v) correction, namely the first endpoint EMSF correction to the integral. There are only exponentially small corrections to the integral representation of Π 1 , but since Π R 2 does not include an n = 0 term, we need to subtract the n = 0 limit of the summand, namely from ln Π R 2 . This gives a multiplicative correction factor of (1 + v) to Π R 2 , so we find that for small velocity Thus the leading small v behavior of ∆ is completely independent of η. However, it has ∆ as a strictly decreasing function of v. As we shall see in the next section, the turnaround for larger v is a nonperturbative effect. For now, we will conclude this section by showing in Fig. 3 a plot of the small velocity region of the graph for various η's, together with our analytic approximation. We see that the analytic result is confirmed.

VII. LARGE η LIMIT
We now turn to a study of the large η limit. In this limit, as first pointed out by Pla,et. al. [15], the relevant variable is ηv. Thus, we study the limit η → ∞, v → 0, φ ≡ ηv fixed. As we shall see, this calculation will shed much light on the small v results we obtained in the previous section.
To begin the calculation, we need the q ∞,n 's and q 0,n 's at v = 0 that we obtained in the previous section, Eq. (65). Then Similarly, Thus, We can invert this relation, solving for φ in terms of ∆, which yields For large ∆, this approaches 1 2 (∆/∆ G ) 2 − √ 2 . This asymptotic result, which is also presented in Fig. 2, is to be contrasted with the result of our continuum calculation, where we found φ = 1 2 (∆/∆ G ) 2 − 1 . Thus, the continuum infinite-η calculation for all ∆ essentially reproduces the large-∆ limit of the lattice calculation, with the correct functional dependence, but with the graph just shifted down slightly. It is also worth noting that including just the n = 0 term, instead of the whole infinite product, also gives the same result, with an intercept of 2/ω which is intermediate between the continuum calculation and the exact asymptotic result. As ∆ decreases, the true η = ∞ curve falls below the asymptotic result, so as to intercept the ∆-axis at ∆| 0 + . The approach is singular, as can be seen by looking at Eq. (86) for small α. We find with an essential singularity at small φ.
Examining Fig. 2 more carefully, we see that our infinite η result has failed to capture one of the most salient features of the finite η data, namely the subcritical nature of the bifurcation from the arrested state. Instead, it possesses a (very-) marginally supercritical onset of the moving crack. To reproduce the subcritical bifurcation from our analytics, we need to generate the next order correction in 1/η.
We begin by generating the next order correction to the q's. We find that q ∞,n does not change to this order, but now q 0,n ≈ 2πin + ω 0 where This induces a multiplicative correction to ∆ of We are interested in the effect of this correction at small φ, in which case we are again free to replace the sum by an integral. If we add in the n = 0 term to the sum, the error will be exponentially small in 1/φ. So, up to exponentially small terms, the correction for small φ is (defining α = 2πnφ) The integral vanishes, as can be seen by a substitution of variables x = (1 + α 2 ) −1/4 . In fact, the integral is nothing more than the first-order expansion in 1/η of the integral in Eq. (75) which we found vanishes identically in η. We are thus left with a correction factor of simply (1 − φ/2η) = (1 − v/2) up to exponentially small terms. This is precisely the small v correction we found in the previous section. The full behavior to this order for small φ is thus This has the subcritical bifurcation we are seeking. As φ increases from 0, ∆ decreases from ∆| 0 + due to the influence of the second factor, until the exponential kicks in and causes ∆ to turn around and start increasing. The φ at which the turn-around occurs is, for large η, of order 1/ ln η (translating to a velocity of order 1/η ln η) which goes to 0 as η goes to ∞, but very slowly. Thus at infinite η there is no turnaround and ∆ strictly increases with φ as we found in the zeroth-order calculation at the beginning of this section. The minimum ∆ lies, for large η, an amount of order 1/η(ln η) 2 below ∆| 0 + .
Thus we see that it is the subdominant pieces that are responsible for the increase of ∆ with v, while the perturbative pieces give rise to the subcritical bifurcation. Analyzing the subdominant pieces in a little more depth, it is easy to see that for η > √ 2/2 the leading subdominant piece goes as exp(−1/ηv). For smaller η, the subdominant piece falls less rapidly, and has an oscillating component, due to the off-axis branch cut assuming dominance. This picture is consistent with the numerical evidence.

VIII. STOKES VISCOSITY
It is worthwhile to contrast the behavior we have seen for Kelvin viscosity with that which obtains for Stokes viscosity, where the dissipation is associated with the mass points and not the bonds. The calculation in this case is much simpler, since the troublesome η term is not present. For our purposes, it is sufficient to consider our x-continuum theory, as the conclusions we obtain carry over to the full lattice model. The result forũ + is where now the Q's satisfy the dispersion relation (and the q's the parallel form with λ m ) and b is the Stokes viscosity. This can be seen by a simple limiting procedure applied to Eq. (23), or by replaying the derivation leading up to Eq. (41) of [11] with b instead of η. This form of the solution can be shown to be equivalent to that obtained by Marder and Gross [9]. This result leads to the solution for ∆: We are interested in the large-N limit, which we obtain be defining the renormalized product since the Q 2 's are linear in Λ for small Λ. Applying the EMSF, we find that for large N , so that The key difference between this formula and the parallel one for η is that ∆/ǫ is proportional to N , and not N 1/2 as before. The reason for this is that the Stokes viscosity is most effective at damping small wavelengths, and so affects the macroscopic stress fields. The Kelvin viscosity does not damp out small wavelengths and only acts on short wavelengths. Another way to see this is to compute the stress intensity factor, which in the Stokes case is inversely proportional to √ b. The driving force required to propagate the crack is thus much larger in the Stokes case. In particular, in the Stokes case there is no macroscopic scaling limit, where things just scale with the Griffith driving, ∆ G . For these reasons, we feel that the Stokes viscosity is not a good model of dissipation for studying crack propagation.
The only way to obtain a nice macroscopic limit where ∆ scales like ∆ G is to artificially scale b with N so that b = b 0 /N . However, this procedure has no physically satisifying motivation, especially when the Kelvin viscosity model suffers none of these defects.

IX. CONCLUDING REMARKS
We close by making a few comments about this work and prospects for future extensions. First it is important to note that the present work is limited to a consideration of the steady-state crack. Thus, aside from general issues of the size of the process zone, the major output of this problem is the velocity-driving relation. Here the most striking qualitative effect of Kelvin viscosity is near threshold, reducing the extent of the backward bifurcation. Significantly above threshold, the major role of viscosity is to provide a velocity scale, so that the crack velocity becomes inversely proportional to the viscosity. It is important to understand how viscosity impacts on the stability of the crack. It is clear, as Marder and Gross have pointed out [9], that the steady-state crack is unstable in the regime of the backward bifurcation. The more interesting question is in the higher-velocity regime. Here, no systematic studies have been done to examinne the role of viscosity. It is not clear that the piecewise-linear model considered here is altogether appropriate for studies of stability, as instabilities can be masked by inconsistencies of the steady-state solution. We look forward to reporting on work in this direction soon, along with generalization to the problem of Mode I cracking.