{'fraction_non_alphanumeric': 0.044717678486478576, 'fraction_numerical': 0.015954916383064368, 'mean_word_length': 4.186562915164168, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Steady-state mode III cracks in a viscoelastic lattice model

We extend the Slepyan solution of the problem of a steady-state crack in an infinite ideally brittle lattice model to include dissipation in the form of Kelvin viscosity. As a demonstration of this technique, based on the Wiener-Hopf method, we apply the method to mode III cracks in a square lattice. We use this solution to find the critical velocity at which the steady-state solution becomes inconsistent due to additional bond-breaking; this point signaling the onset of complex dynamical behavior.


I. INTRODUCTION
Recent experiments on dynamical fracture (Fineberg, et al., 1991(Fineberg, et al., , 1992 have shown that cracks in brittle materials become unstable above a critical velocity. This instability is associated with the formation of a roughened fracture surface, with an increased dissipation of elastic energy, and with complicated velocity dynamics -for a review, see (Fineberg and Marder, 1999). Although there are hints of such an instability in the traditional continuum formulations of ideally brittle cracks (Yoffe, 1951), a systematic treatment does not appear possible. Indeed, recent attempts (Langer and Lobkovsky, 1998) to utilize a cohesive stress modification (Barenblatt, 1959) of the continuum elastic equations in order to address this problem have been rather unsuccessful.
In the absence of a compelling continuum approach, one must deal in some manner with discrete dynamics on the "atomic scale". In this regard, Slepyan (Kulamekhtova, 1984;Slepyan, 1981Slepyan, , 1982 pioneered the idea of studying lattice models in which atoms interact (via piece-wise linear springs) with their neighbors in a predetermined lattice geometry. This conceptual framework was further developed by Marder and collaborators (Marder and Gross, 1995;Marder and Liu, 1993). Here, one can directly obtain the relationship between the crack tip velocity and the imposed driving and furthermore one can search for conditions which do not allow stable steady crack motion. Clearly, lattice models are not fully realistic, especially at displacements that are not small compared to the underlying lattice spacing; to be more realistic, one must resort to molecular dynamics simulations including all possible atomic interactions (Abraham, et al., 1994;Gumbsch, et al., 1997;Marder, 1997, 1998;Omeltchenko, et al., 1997;Zhou, et al., 1997). However, the analytic tractability of these models as well as indications (Marder and Gross, 1995;Pla, et al., 1998;Sander and Ghasias, 1999) that they do contain the essential mechanism responsible for the observed dynamical instability make them well worthy of serious attention.
Slepyan recognized that in an ideally brittle material (where each spring is linear until a displacement at which it completely breaks) one can use Fourier methods to find the lattice analog of a traveling wave solution. In this solution, each point on the lattice (at the same transverse coordinate) undergoes the same time history of motion as any other, albeit with some time delay. Once one obtains the solution, one must check that all springs assumed to be linear have displacements that are in fact below the breaking threshold. The violation of this assumption by the steadily propagating solution signals the onset of more complex dynamical behavior, at least in qualitative accord to what is seen experimentally (Fineberg and Marder, 1999;Marder and Gross, 1995). Using the Wiener-Hopf technique, Slepyan solved for the velocity-driving curve in the limit where the width of the lattice transverse to the crack direction goes to infinity.
One important consideration in all models of brittle fracture concerns dissipation mechanisms. From our perspective, it makes sense to put in dissipation at the lattice scale in such a way so as not dominate the large-scale continuum elastic field. If in addition one demands an equation that is local in time, one is led (Langer, 1992;Pla, et al., 1998) to the introduction of a Kelvin viscosity which dissipates energy proportional to the rate of change of spring lengths. In the naive continuum limit, this gives rise to a third derivative (two space, one time) term. That means that if the viscosity is chosen to be O(1) on the lattice scale, it will scale to zero as far as the macroscopic dynamics is concerned. We saw this explicitly in a previous set of papers (Kessler, 2000;Kessler and Levine, 1998) which solved this problem for finite width lattices and considered the nature of the solution as the width became large. Nevertheless, the viscosity can have a considerable effect on aspects of the solution which depend on the microscopic details, namely the crack speed (for a fixed stress intensity factor) and the self-consistency of the traveling wave ansatz.
The rest of this paper is organized as follows. In the next section, we formally define the model, employing Slepyan's idea (Kulamekhtova, 1984;Slepyan, 1981) of replacing the driving at some external boundary with some local forcing on the crack surface. This allows for the solution in Section III of the strictly infinite lattice model; the connection to a finite problem with some displacement driving is made by appealing to the universality of the microscopic crack solution given a fixed stress intensity factor. In the following section, we discuss the numerical evaluation of the velocity curve from its formal expression. Following that, we turn to the aforementioned self-consistency condition and study the effect of dissipation on the critical velocity. We conclude with some observations in the final section.

II. SQUARE LATTICE MODEL
We wish to study mode III cracks and the effects thereupon of Kelvin viscosity. Our model thus consists of a square lattice of mass points undergoing (scalar) displacements out of the plane. The lattice extends infinitely long in both the x (along-crack) and y (transverse) directions. The lattice points are connected to their nearest neighbors by ideally brittle springs which behave linearly with spring constant 1 until some threshold elongation 2ǫ at which point they irreversibly crack. All the (uncracked) springs have a viscous damping η. Let us assume that the crack corresponds to a sequential breaking of the vertical bonds between the mass points at rows y = 1 and y = 0. The equation of motion for the masses in rows y > 1 reads We allow the vertical bonds along the crack surface to have a different spring constant k and damping parameterη, to model the effect of having weak links along the crack surface. The equation for y = 1 then reads, using the assumed symmetry u(x, y) = −u(x, 1 − y), It is easily verified that the equations for y ≤ 0 are consistent with the above symmetry. We can, without loss of generality choose ǫ = 1/2 (so that the threshold elongation is unity). 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 of time such that u(0, 1) = ǫ, so that it represents the moment of cracking of the spring in column 0. If we define τ = x − vt, the equations of motion for y > 1 become 1 − ηv ∂ ∂τ u(τ + 1, y) + u(τ − 1, y) + u(τ, y − 1) For y = 1, we separate out the terms proportional to θ(−τ ), giving 1 − ηv ∂ ∂τ u(τ + 1, 1) + u(τ − 1, 1) + u(τ, 2) − 3u(τ, 1) The driving term for this equation has the form In this expression, we have inserted an external force F 0 which acts on the crack boundaries. This force is an artifice which serves as an inhomogeneous source term, in place of introducing the driving through the boundary condition on the top and bottom surfaces. We will see later how to choose this function. We Fourier transform with respect to τ (see the appendix for our conventions regarding Fourier transforms and Fourier integrals). If we assume that u(q, y) = χ q [ξ(q)] y−1 for y ≥ 1, we obtain from Eqn. (4) the dispersion relationship This can be rewritten as which therefore implies that ξ 1/2 = ±h/2 ± (h/2) 2 + 1. If we define Note that we have must choose the sign that corresponds to a decaying solution in the y direction, i.e. |ξ| < 1. If we substitute this ansatz into Eqn. (5), we get , we obtain with This equation is implicit, since the driving term σ F on the right hand side still depends on the displacement field u(τ, 1). In the next section, we solve this equation via the Wiener-Hopf technique.

III. WIENER-HOPF SOLUTION
Let us define ∨(τ ) = 2u(τ, 1), the bond elongation between rows y = 0 and y = 1. In the appendix, we define ∨ + (q) and ∨ − (q) as the decomposition of ∨(τ ) into terms which are analytic in the upper and lower half planes respectively. This allows us to write Eq. (6) as where we have employed a similar breakup for the driving term F 0 , and where we have utilized Eq. (62). Substituting the formula for σ F into Eqn. (10) for χ = (∨ + + ∨ − )/2 leads to a closed form equation for ∨; namely Defining S by we can rewrite Eq. (14) as We will solve this equation using the Wiener-Hopf technique. The key is to factor S into a product of two pieces, S(q) = S + (q)S − (q), each of which is regular in the upper (S + ) or lower (S − ) half plane. Doing so, we get We then have to deal with the extra singularity at q = ik ηv . We do this by a simple subtraction We have now to choose F 0 , our external forcing. In taking the width to infinity, in essence we are solving for the "inner" solution in the sense of boundary-layer theory; i.e., only on the scale of the lattice spacing. The external forcing in the finite width problem only acts on the large scale, and does not vary on the lattice scale. Thus, F 0 must have support only at q = 0. Following Slepyan, then, we choose F 0 such that where B is a constant. This then gives us the the two separate equations As shown in Eqs. (59) and (61) we can solve for ∨(0) from Eqs. (20a) or (20b) by multiplying the first by −iq → +∞ or the second by iq → +∞. This yields Substituting this back into Eqns. (20), and defining η k ≡η/k, we find These equations directly determine the displacements given the strength of the external driving, B. To get a physical handle on the meaning of this constant, we note that Eqns. (20) can be directly solved for the leading behavior of ∨ ± as q goes to 0, using the fact that in this limit, , and κ is a constant, reflecting the arbitrariness in how we perform the decomposition of S. We find This is just the expected (Fourier transform of) the singular solution of continuum elasticity as one approaches the crack tip. Thus, as in all boundary-layer problems, the large-distance limit of the "inner" solution matches on to the short-distance limit of the "outer" solution, confirming the justification of the choice of F 0 above. We need to relate the strength of the short-distance singularity of the continuum solution to the driving displacement ∆. This can be done directly by solving the continuum elastic equations. It is easier, however, to follow Slepyan and derive this relation by energy considerations. One can calculate (Slepyan, 1981) the flux of energy into the boundary-layer, or "process" zone, obtaining Now consider the flux of energy T 0 being used to break bonds, which is simply related to ∨(0) Then, using the idea that the Griffith's displacement is determined by exactly the condition that all the energy flux is needed just to break the bonds, we have Notice that the factor of κ guarantees that the result is invariant with respect to how the decomposition of S is performed.
In order to make contact with other treatments, we consider the standard case k = 1 and η k =η = η, in which all bonds are equivalent. In this case, the function S has a particularly simple form, with H 2 = 4(1 + iqηv) sin 2 q 2 − q 2 v 2 and R 2 = 4(1 + iqηv)(sin 2 q 2 + 1) − q 2 v 2 . We can represent these functions in term of the products over all their roots where Q + n (q + n ) are the (nonzero) roots of H 2 (R 2 ) in the lower half plane and Q − n (q − n ) are corresponding roots of H 2 (R 2 ) in the upper half plane. This decomposition allows us to write down explicitly the factors S + , S − : Examining the small q behavior, we find the κ = 1, so that substitution of S + into Eqn. (27) gives If we explicitly take out of the product in Eq. (32) the one imaginary q + , which we will denote as q + 0 , we get This formula was obtained by Kessler (2000), who started from a finite lattice of transverse size N y and considered the N y → ∞ limit. Note that this can be done for k = 1 on a square lattice, (the case being considered here), but not for general k or for modes I or III on a triangular lattice; the method here works for these cases as well (Pechenik, et al., 2000). Next, it is easy to see that in the limit η → 0 + , we have S| q= i ηv = 1 and that therefore Furthermore, in this limit q + 0 = −q − 0 . All other complex roots cancel out, leaving only real roots. Hence, a result originally obtained by Slepyan. Note that in our treatment, the presence of one of these real roots in either the numerator or the denominator of the product expression depends entirely on the half-plane from where it came as we took the limit η → 0 + . To see what this condition implies, let us denote by H 2 (or R 2 ) as f ; then the small imaginary part of a specific real root, iδ, satisfies the equation which gives as the determining factor This condition is exactly what appears in Slepyan's work.

IV. NUMERICAL CALCULATIONS
In this section we derive an integral form for the basic result Eq. (32). We then present a numerical procedure for finding the displacement-velocity curve. The basic notion is to utilize Eq. (63) to find S + . It is more convenient for the numerical work to isolate explicitly the singularity at q = 0, defining K(q) by where φ is an arbitrary positive constant. Because of the factor A 2 q 2 in the denominator, K(0) = 1; the factor (q 2 + φ 2 )/φ 2 does not alter the desired asymptotic behavior at infinity: K → const as q → ±∞. Applying Eq. (63), we get for Im q > 0 and for Im q < 0 where the upper sign is for Im q > 0 and the lower sign for Im q < 0. There is one point regarding these expressions that must be clarified. We must make the correct choice of the branch cuts for the square roots in h and r appearing in S. It is convenient to express things in terms of S * (q) ≡ h/r = H/R, the simple k = 1 and η k = η, form of S, so that from Eqn. (15) and from Eqn. (9) We must choose S * such that |ξ| < 1 as noted below Eqn. (9). It is easy to see from the above expression for ξ that it is sufficient to choose S * to lie in the first and fourth quadrants to guarantee this. For S * , we can then generate S directly via Eqn. (42), and from this K.
As K(0) = 1, there is no pole in the integral (41) at γ = 0 when q = 0. This means that K + (0) = (K − (0)) −1 and the asymptotic expressions for S ± are Comparing this to Eqn. (23), we have so that Explicitly writing out the integrals, we get The extra degree of convergence of the integrand at infinity is very desirable, as it helps control the fact that ln K(ξ) is a rapidly oscillating function. Also, the fact that the answer must be independent of φ provides a check on the numerical routines. A change of variables γ → −γ allows the integral to be rewritten as To actually compute this integral, we divided the region of integration into two parts from 0 to 1 and from 1 to +∞. The second integral was further transformed so as to have in the numerator the factor ln A 2 φ 2 K(ξ) which goes to zero at infinity; the subtracted integral with the numerator ln 1/(Aφ) 2 was performed analytically. After all these manipulations, the integrals were successfully computed by standard mathematical library subroutines. Results from our numerical calculations for k = 1 and η k =η = η, are presented in Figs. 1, 2, which reproduce, of course, those of (Kessler, 2000). At small η, there are large oscillations in the v vs ∆ curve, as found initially by Marder and Gross for lattices with small transverse sizes. The curves all hit the v = 0 line at the point which marks the end of the band of lattice-trapped static cracks (Kessler and Levine, 1998). For larger damping, the oscillations disappear and the moving crack branch bifurcates smoothly (in the backward direction) from this point, eventually turning around and becoming stable. In Fig. 3, we contrast the standard result to that for smaller k, keeping η k = η. We see that the v = 0 limit has moved closer to unity, in accord with the finding in (Kessler and Levine, 1998) that the window of arrested cracks narrows as k decreases (note for purposes of comparison that our k is a factor of 2 smaller than that defined in (Kessler and Levine, 1998). The velocity at the minimum ∆/∆ G however is hardly affected, so the minimum stable velocity is essentially unchanged. Past this point, the velocity is higher at smaller k, even when the change in ∆ G is factored in. This is further evidence of how the microscopic details influence the velocity versus driving displacement curve.

V. CONSISTENCY OF THE SOLUTION
As mentioned in the introduction, cracks become unstable above a critical velocity. In the framework of the model being considered here, one might assume that this velocity is associated with the solution as found by the Wiener-Hopf method becoming inconsistent. As we shall see, there are two kinds of inconsistencies. The first sets in for small v, for small enough η. Here the vertical bond between y = 0 and y = 1 first achieves critical extension prematurely, at some positive τ , contradicting the assumption that the bond breaks at τ = 0. In the second kind of inconsistency, a horizontal bond achieves critical extension. If the model is defined so that only the vertical bonds between y = 0 and y = 1 are breakable, then this does not present a problem. If one the other hand, all the springs are chosen identical, so that the horizontal springs also break at the same critical extension as the vertical ones, then the solution is indeed inconsistent at this point. This is what occurs at large enough velocity, for all η.
To proceed, we must calculate the bond lengths, back transforming our Fourier space solutions so as to obtain the physical displacements. We find the time dependence of the function ∨(τ ). From Eqs. (20a), (20b) and (21), we have Performing the inverse Fourier transform, we get for τ > 0 and for τ < 0. Here we have used the expressions for S ± derived in the last section.
Clearly, these integrals must be done numerically. To proceed, we transform the pieces containing the factors K ± by dividing both the numerator and denominator by (K + (0)) 1/2 . In the denominator, we get ∆ G /∆; we have already discussed how this can be computed. We now have integrands containing factors of the form The ± refers on the right hand side to the sign of the imaginary part of q. As we will see later on, we need to calculate this function only for positive q. Then, the integral in the exponent can be written as We break up the region of integrations into three parts,(0, 1/2), (1/2, 3/2) and (3/2, +∞). In the first interval, we numerically calculate the integral directly as written. For the third interval, we proceed as discussed in the last section for a similar semi-infinite integral; this leaves us with having to integrate numerically a function which behaves as 1/ξ 4 near infinity. Finally, for the integral near 1 we add an integral with K(qξ) replaced by K(q); this added integral is done analytically and the resultant subtracted integral with the integrand now containing ln K(qξ)/K(q) is done numerically.
Applying the described procedure, we can evaluate the functions (K ± (q)) ± 1 2 / (K + (0)) 1 2 for any real positive q. We must now perform the final integration over the variable q. First, we note that the requirement of having a real displacement necessitates ∨ ± (−q) = ∨ ± (q) * . This allows us to change the region of integration to (0, +∞) in (51); we will return shortly to 52. The trickiest part of this calculation concerns the behavior near q = 0. For ∨ + , there is a 1/q 1/2 behavior. Here the leading order term can be integrated analytically over the interval (0, 1) and then the overall function with this term subtracted can be used for a numerical integration. This works in a straightforward manner.
We need to be more careful with the integral (52), because here the integrand behaves near zero as 1/q 3/2 . To proceed, we first divide the interval of integration into three parts: (−∞, −1), (−1, 1), and (1, +∞). The integrals over (−∞, −1) and (1, +∞) can again be combined to yield an integral twice the real part of the integrand over the latter range. For the range spanning zero, we first subtract from the integrand the leading term After this subtraction, the integrand becomes of order 1/q 1/2 . Now, the one-sided integral converges, and thus we can again transform the range to (0, 1). This subtracted integral is evaluated numerically from a very small value of q (typically q = 10 −3 ) to 1. In the remaining part of the range, i.e. from 0 to 10 −3 , we evaluate the integral analytically after approximating it by its leading small q behavior, α/q 1/2 ; the value of the constant α is determined by fitting the the behavior of the function near the point 10 −3 . Finally we need to perform analytically the integral of the subtracted piece 1/( To accomplish this, we deform the contour of integration from (−1, 1) on the real axis to the curve from −1 to 1 going over the lower half of the unit circle. Note that the branch-cut for this function must be taken as before over negative real axis; this choice guarantees the fact that f (−q) = f * (q). Note that once all our transformations are complete, we only need values for the integrand for positive (real) values of q; this was used earlier in the technique for calculating K ± . A good check of this complex numerical technique of computing the Fourier transformation is applying it to the Fourier transform ∨ + (q) for τ < 0 and ∨ − (q) for τ > 0, as in this case result must be equal to zero.
We first examine what happens in the standard case, k = 1 and η k =η = η. Figure 4 shows the behavior of ∨(τ ) for η = .2 and various values of the velocity. As discussed above, we see the two different kinds of inconsistencies exhibited by these solutions. For example, the solution for v = 0.2 is inconsistent as ∨(τ )/ ∨ (0) > 1 for τ > 0. This type of inconsistency, seen at small velocities, corresponds to the region of backward behavior on v − ∆ graph. As already claimed by Marder and Gross, this region is therefore unphysical.
More interesting is the inconsistency which sets in at large velocity, and has been argued to be related to the experimentally observed microbranching. What happens here is that the elongation of the horizontal bond, given by ∨ hor (τ ) = (∨(τ − 1) − ∨(τ ))/2 increases as a function of increasing driving. At some critical speed it reaches the value ∨(0) and the bond breaks. Figure 5 shows the dependence of the critical speed on η. We note in passing that these graphs agree with what can be gotten by extrapolating a similar graph for the finite N y lattice, a graph that can be obtained using the methods of (Kessler and Levine, 1998). The actual location in τ at which the horizontal break occurs at the critical velocity is shown in Fig. 6.
Several details of our findings are worthy of note. First, there is always a critical speed for any dissipation, but this speed gets very close to the maximal crack speed for large η. The value of the critical speed for the dissipation-less limit is v cr ≃ .725; this is above what has been seen in some experiments but recall that we are doing mode III on a square lattice, a far cry from mode I in an amorphous system. The is a curious break in the curve at η ≃ .665, which corresponds in Fig. 6 to the point where the break location varies most strongly with η.
It should be noted that this inconsistency does not appear related in any obvious way to the Yoffe (1951) criterion. First, the Yoffe criterion, concerning the direction of maximal stress, is calculated for the continuum elastic field, and is thus completely independent of the viscosity η (Kessler and Levine, 1998). The onset of the inconsistency is strongly η dependent however. Also, the Yoffe criterion is stated for Mode I cracks, and its applicability to Mode III is subject to question.
We turn now to the case of k < 1, with weakened bonds along the crack path, again keeping η k = η. Figures  5, 6 show our results. As is reasonable, the inconsistencies persist, but are pushed to higher speeds as k decreases. It is interesting to speculate that since for finite width systems, velocities greater than the wave speed are possible, sufficiently small k might enable stable steady-state fracture at supersonic velocities. Also, intersonic wave speeds are possible in the case of mode II fracture, and therefore weakened bonds on the path of the crack would in this case as well allow steady-state propagation at higher speeds (Rosakis, et al., 1999).

VI. DISCUSSION
We have shown how to extend the Slepyan approach to cracks in ideally brittle infinite lattice systems to the case of having a Kelvin viscosity for each lattice spring. In addition, we have examined the onset of the additional, off-axis cracking at high velocity. We showed how dissipation delays but does not eliminate the microbranching instability, an effect impossible to see within the context of any continuum elastic treatment.
In the following work, we will present the results of similar calculations for mode I cracks on a triangular lattice. That system is much closer to those that have been studied experimentally and one can therefore hope for more direct lessons to emerge. Also, we should mention complementary studies on nonlinear lattice models which relax the assumption of ideally brittle springs and thereby allow one to view the onset of inconsistent behavior as the limiting case of more standard bifurcations (Kessler andLevine, 1999, 2000).

Acknowledgments
DAK acknowledges the support of the Israel Science Foundation. The work of HL and LP is supported in part by the NSF, grant no. DMR94-15460. DAK and LP thank Prof. A. Chorin and the Lawrence Berkeley National Laboratory for their hospitality during the inital phase of this work.

VII. APPENDIX
We discuss here some of our conventions and notations regarding Fourier transformations. We use Fourier transformation in the form We define also u + and u − as This gives which means that u + has poles only in the lower half plane and u − only in the upper half plane. Now let us define p = −iq. Then Similarly we can get for p = iq The Fourier transform of θ(−τ ) d dτ u(τ ) is given by Finally, we consider separating an arbitrary function S(q) into the product of two pieces S + and S − , with poles and zeroes in lower and upper half-planes respectively. These can be done using the identity, f or Im q > 0, ln S + (q) f or Im q < 0, − ln S − (q) .