## Abstract

We present the full analytical solution for steady-state in-plane crack motion in a brittle triangular lattice. This allows quick numerical evaluation of solutions for very large systems, facilitating comparisons with continuum fracture theory. Cracks that propagate faster than the Rayleigh wave speed have been thought to be forbidden in the continuum theory, but clearly exist in lattice systems. Using our analytical methods, we examine in detail the motion of atoms around a crack tip as crack speed changes from subsonic to supersonic. Subsonic cracks feature displacement fields consistent with a stress intensity factor. For supersonic cracks, the stress intensity factor disappears. Subsonic cracks are characterized by small-amplitude, high-frequency oscillations in the vertical displacement of an atom along the crack line, while supersonic cracks have large-amplitude, low-frequency oscillations. Thus, while supersonic cracks are no less physical than subsonic cracks, the connection between microscopic and macroscopic behaviour must be made in a different way. This is one reason supersonic cracks in tension had been thought not to exist.

## 1. Introduction

One of the conclusions frequently drawn from the continuum theory of fracture mechanics is that cracks in tension cannot travel faster than the speed at which sound travels over a flat surface, the Rayleigh wave speed [1–3]. The reason for this assertion is that the motion of a crack requires energy to break bonds, but energy flux seems to become nonsensical for supersonic cracks. In particular, the rate at which energy flows into a crack tip per time is
1.1
where *v* is the crack speed, *c*_{l} and *c*_{t} are longitudinal and transverse wave speeds, respectively, , , *μ* is a Lamé constant, and *K*_{I}, the mode I dynamic stress intensity factor, is the coefficient of a universal singularity that develops outside of cracks that run as they are pulled symmetrically in tension from above and below. The denominator of this expression vanishes when the crack speed *v* reaches the Rayleigh wave speed, and for slightly higher velocities it becomes negative. Once the crack speed exceeds the transverse wave speed, the expression becomes imaginary. An expression saying that cracks moving above the Rayleigh wave speed need negative energy seems physically impossible, and an expression requiring imaginary energy seems even worse. One resolution of these problems is to conclude that cracks travelling faster than the Rayleigh wave speed in tension are not physically allowed.

However, there is evidence that such cracks exist after all. The first indications came from measurements of earthquakes [4]. These suggested computer simulations [5,6], and laboratory experiments [7] for cracks that move by sliding faces past each other, showing that cracks in shear (mode II) can move faster than the transverse wave speed *c*_{t} (and hence the Rayleigh wave speed as well) and reach speeds close to the longitudinal speed *c*_{l}. Dynamic fracture theory was extended to include these ‘intersonic’ cracks [2,8,9].

Cracks in rubber under tension were found to have a wedge-like tip suggestive of supersonic motion [10]. Additional experimental work confirmed that the cracks do travel faster than the transverse wave speed [11,12]. The experiments led to theoretical descriptions for supersonic cracks in tension [13,14]. Both explicit numerical solutions for atomic equations of motion and the corresponding analytical solutions show the supersonic cracks do in fact exist [15]. This last reference contains an extended discussion of how crack speed depends upon loading for the models studied in this manuscript.

Supersonic cracks were also observed for cracks with hyperelastic constitutive laws that cause stiffening near the tip [16]. As shown, for example, in the current article supersonic cracks exist equally well in materials without hyperelastic constitutive laws. Thus, we are not sure that hyperelasticity is actually needed to promote supersonic cracks.

Furthermore, the original objection from continuum theory about the impossibility of supersonic cracks has not directly been addressed. Doing so is a goal in this article. What we will explain in detail is the specific way that the continuum argument breaks down. In particular, we will show that, as a crack approaches the Rayleigh wave speed, the stress intensity factor *K*_{I}, upon which the continuum theory depends, stops being defined. The deformations around the crack tip take a completely different form from the one they have in the subsonic regime.

Our main technical tool to provide this discussion is the atomic theory of dynamic cracks. Analytical solutions at the atomic level for moving cracks were first found by Slepyan. His original calculations applied to out-of-plane (mode III) cracks in an infinite square lattice [17]. Slepyan and co-workers then generalized the solutions to in-plane (mode I) cracks in an infinite square lattice [18]. Many additional solutions due to Slepyan are found in [19].

Additional solutions for cracks in a finite strip formed from a triangular lattice were found by Marder & Gross [20]. This paper also considers linear and nonlinear instabilities of the solutions. Many additional observations on the experimental implications of the solutions in discrete lattices are contained in [21].

In all these analytical solutions, the central quantity obtained from the analytical methods is a relationship between the external load applied to make a crack move and the speed *v* at which the crack travels. In fact, the speed *v* of the crack enters the theory as an input parameter, and the calculations give the load as an output. For subsonic cracks, one can equivalently say that one obtains *K*_{I}(*v*), the stress intensity factor as a function of crack speed. For supersonic cracks, the end result is the external system strain *e*_{yy} as a function of crack speed [14].

The theory provides more than just this relation. In principle, the theory makes it possible to find the behaviour of every atom as a function of time for a crack moving in steady state. This information is relatively easy to obtain for cracks moving in anti-plane shear [20]. However, for two-dimensional cracks where atoms move freely in the plane, the expressions are so lengthy that to our knowledge no one has completed the process of extracting analytically the motion of every atom.

Thus, in this paper, we have two aims. The first is to present for the first time analytical solutions for every atom in a two-dimensional lattice where a crack moves in plane. The algebra is extremely lengthy, so we just present an outline of the work. Then, as an application of these exact solutions, we turn to the transition between subsonic and supersonic cracks. We examine the precise way that the stress intensity factor vanishes across the transition point and show what the cracks look like once they pass the Rayleigh wave speed. Because the methods produce solutions very rapidly, we are also able to discuss how these phenomena depend upon system size.

The structure of the article is as follows: in §2, we describe briefly the means we have used to solve for atomic motions around steady-state in-plane lattices. In §3, we examine the structure of the solutions in the vicinity of the Rayleigh wave speed, and explain how the stress intensity factor disappears. In §4, we make final remarks and conclude.

## 2. Explicit in-plane solution

### (a) Lattice system

We carry out our computations in a triangular lattice of atoms with 2(*N*+1) rows. The motion of each atom is described by the displacement **u**_{m,n} from its equilibrium position, with *m* the column and *n* the row. The index *m* takes integer values while *n* takes values of the form (*k*+1/2), with *k* an integer ranging from −(*N*+1) to *N*.

The atoms in the lattice interact linearly with their nearest neighbours. The interaction is a function of both parallel and perpendicular displacements, with two spring constants *k*_{∥} and *k*_{⊥}. We must make the restriction *k*_{⊥}=0 along the crack line, otherwise the mathematical formalism fails. We can make up for this by setting for the two rows along the crack. The bonds between atoms snap when they are displaced a distance 2*u*_{f}, and the interaction becomes zero.

We note that, in the long-wavelength limit, the motions of atoms in our lattice are described by isotropic continuum elasticity. In particular, longitudinal and transverse wave speeds are given, for unit mass and lattice spacing *a*, by
2.1a
and
2.1b
In turn, this correspondence implies that crack motion in this lattice contains all results of dynamic linear elastic fracture mechanics as a special case.

To study crack motion, one must impose loading on the system. We do this by displacing the top row of atoms in our strip by a fixed vertical amount *U*_{N}. The crack consists in a separation between the middle rows of the lattice, where *n*=±1/2. We study the system when the crack moves at constant speed, in a steady state. Steady crack motion at velocity *v* in the continuum means that elastic fields are functions of *x*−*vt*. For cracks in a lattice, one cannot employ this definition. However, there is a related symmetry, which is that in a lattice of spacing *a*
2.2
This means that what an atom at some height does now, its neighbour to the right will repeat exactly at a time *a*/*v* later. Employing this relation lets us eliminate the *m* index and relate the components of different columns in the lattice to each other by a simple time symmetry. In particular, since
2.3
we will be able to drop the index *m* and set about finding **u**_{n}(*t*).

### (b) Analytical steps

The lattice in which we solve for crack motion is shown in figure 1. The displacement, for example, in the first direction is **Δ**_{1}=**u**_{m−1,n+1}−**u**_{m,n}. To find the force on an atom, define **e**_{∥j} and **e**_{⊥j} to be unit vectors in the unstretched lattice. The unit vectors **e**_{∥j} point from atom 0 to atoms *j*=1…6 in figure 1. Each of the unit vectors **e**_{⊥j} is perpendicular to the corresponding **e**_{∥j}. That is,
2.4
The total force on any given atom, except those along the crack line where some bonds may be broken, is
2.5
and the motion of atoms, of unit mass, is given by
2.6
The addition of the Stokes dissipation through the second term on the right-hand side is necessary to break time reversal invariance and tell the crack whether it should advance or retreat as time runs forward.

Define *x* to be the horizontal direction (the direction along which the crack moves) and *y* to be the vertical direction (along the width of the strip). The forces on atoms above the crack line, *n*>1/2, take the explicit component form
2.7
and
2.8

The number *g*_{n} helps keep track of the location on the triangular lattice. It is defined by
2.9

Much of the algebra needed to solve this system appears in detail in [20]. There is a missing factor of 12 in the equations of motion (VI.21, p. 48) that appear in [20], although the results reported in the body of the paper are correct. The corrected forces along the crack line, *n*=1/2, are
2.10
and
2.11
where *U*(*t*) is a linear combination of *u*^{x}(*t*) and *u*^{y}(*t*) multiplying the Heaviside *θ*(*t*) functions
2.12

Fourier transforming the equations of motion gives 2.13 and 2.14 where 2.15

These equations can be combined and rearranged to give finally
2.16
Here the function *Q*(*ω*) can be calculated explicitly from the Fourier transformed equations of motion for the lattice, but its explicit expression is so lengthy that it would serve no purpose to type it out here; we compute it from MAXIMA or Mathematica scripts, and for the purposes of rapid numerical evaluation find it in Fortran. In addition, *Q*_{0}=*Q*(0), and *α* is a small positive constant introduced when rewriting the boundary condition as *U*_{N} e^{−α|t|}. The Fourier transform of this function is easier to deal with than *δ*(*ω*). We let *α* tend to zero at the end of the calculation.

The Wiener–Hopf technique [22] lets us write
2.17
where *Q*^{−} is analytic in the lower half-plane, and *Q*^{+} is analytic in the upper half-plane. This decomposition is given by the explicit formula:
2.18

We can split court (2.16) into two pieces, one which is analytic only in the upper half-plane, and one which is analytic only in the lower half-plane: 2.19

The left- and right-hand sides of the above equation are analytic in opposite sections of the complex plane; therefore, they must both equal a polynomial *P*(*ω*). We must have *P*(*ω*)=0, otherwise *U*^{−}(*t*) and *U*^{+}(*t*) will include unphysical *δ*(*t*) terms. This lets us write
2.20
From the two functions *U*^{+}(*ω*) and *U*^{−}(*ω*), one can obtain the linear combination of horizontal and vertical motions defined in equation (2.12) from
2.21
Finally, one can obtain the horizontal and vertical motions of any atom in the system from
2.22
and
2.23
The functions and , like *Q*(*ω*), are very lengthy to write down here, but can be obtained from essentially straightforward algebra resulting from the equations of motion. Again, we obtain explicit expressions using symbolic algebra programs, and implement the results for rapid evaluation in Fortran.

### (c) Addition of the Kelvin dissipation

The solution we have described so far included the Stokes dissipation , which was present in equation (2.6). However, the Stokes dissipation is not physically realistic. It acts as though atoms are embedded in some sort of ether that slows them down according to the magnitude of their velocity in some arbitrarily chosen reference frame. A better form of dissipation from the physical point of view is the Kelvin dissipation, which produces a force opposing atomic motion according to the relative motion of adjacent atoms. To preserve generality, we keep both the Stokes and Kelvin dissipation in the equation of motion, which becomes almost, but not quite, 2.24 The reason this expression is not quite right is that when bonds break the time derivative produces unphysical discontinuities on the right-hand side. A solution of this problem was laid out in [23]. The solution is to define 2.25

Then the equations of motion on the crack line are modified to read 2.26 and 2.27

The functions *W*^{±} are defined just like *U*^{±} in equation (2.15). From this point, the analysis proceeds as before, finding a new function *Q*(*ω*) that is not appreciably more complicated because of the presence of *β*. From this function, one determines *W* using the Wiener–Hopf technique as
2.28
and
2.29
Once *U*(*ω*) is in hand, the motion of every atom can again be determined. There is also a compact result relating the velocity *v* of a crack (which goes into the computation of *Q*) and the boundary extension *U*_{N}. The boundary extension can be expressed in terms of a dimensionless load Δ,
2.30
Then, defining *ω*_{0}=1/*β*, Δ can be computed from *Q*(*ω*) directly as
2.31

### (d) Numerical calculations

All atomic positions in time can be evaluated from rapid numerical operations involving nothing but algebra and fast Fourier transforms. First we sample the function *Q*(*ω*) over an interval of width 2*ω*_{max} and use fast Fourier transforms to compute *U*^{−}(*ω*) and **u**_{n}(*ω*) from equation (2.18). Before taking the inverse Fourier transform, any singular behaviour of the form 1/(*iω*) must be subtracted from . The coefficient is determined from the asymptotic behaviour of and is given by
2.32

The subtracted function is then added back analytically after the inverse Fourier transform. This completes the numerical calculation of **u**_{n}(*t*). We can check the accuracy of the components knowing the exact asymptotic behaviour along the crack line,
2.33

Also, in the limit, the load is given by . This can be compared to direct numerical integration of equation (2.31). We can choose an *ω*_{max} that minimizes the difference between these two Δ values. Figure 2 shows components of the atomic position calculated in this manner. We choose the boundary condition *U*_{N} in such a manner that at the moment bonds break they have increased in length by *u*_{f}=1. Thus, we say that all displacements in the plots are measured in units of *u*_{f}. Keeping the breaking length fixed means in turn that all the computations correspond to physical systems with the same fracture energy. As the number of rows 2(*N*+1) increases, the extension *U*_{N} needed to bring the system to the point of fracture increases as . For supersonic cracks, the extensions become large multiples of *u*_{f}.

It is worth pausing to ask how figure 2 might be produced were one not using the Wiener–Hopf technique. It would be possible. It would require integrating the equations of motion for a crack in a system 1002 rows high. In order for atomic motions to reach steady state, the crack would need to run for a distance around 10 times the height of the system [24]. Although this could be sped up with a cutting and pasting procedure, accurate results would require a system 3000 columns long. Thus, one would have to run to steady state a system with around 3 million atoms; this would require a supercomputer. By contrast, producing figure 2 required just a bit over a second using a single processor on a laptop.

## 3. Behaviour of solutions

### (a) Supersonic transition as a function of velocity

Having developed a tool that allows us very quickly to find the time history of atoms in the vicinity of a running crack, we now turn to the question with which we began. We ask what happens near a crack tip as the motion of the crack moves from subsonic to supersonic motion. The important wave speed is the Rayleigh wave speed *c*. For central forces (*k*_{⊥}=0), it is given explicitly in our model by
3.1
which is the root of the denominator in equation (1.1), and the limiting speed for cracks in continuum fraction mechanics.

Supersonic solutions in the lattice model look very different from subsonic solutions. Figure 3 shows how the vertical displacement of an atom on the crack line varies as the crack speed increases through the Rayleigh wave speed from *v*=0.9*c* to *v*=1.05*c*. We can see two distinct behaviours for subsonic and supersonic cracks. In both cases, the atom is nearly motionless until the crack approaches. After the crack passes, in subsonic solutions, the vertical displacement approaches the boundary condition *U*_{N} with small-amplitude, high-frequency oscillations that continue for a time of the order of 1/*β* as shown in figure 3*a*. These are phonons carrying energy left over after the bonds along the crack line have snapped [25].

In supersonic solutions, the vertical displacement also approaches the boundary displacement, but with large-amplitude, low-frequency oscillations as seen in figure 3*h*. The phonons of the subsonic solutions are no longer present in supersonic solutions. Interestingly, the subsonic phonons and supersonic large oscillations appear simultaneously in the subsonic solutions just below the Rayleigh wave speed, as we see in figure 3*c*–*e*.

The continuum solutions for in-plane fracture are graphed in [2]. Using the same notation as equation (1.1), the vertical displacement is
3.2
where *R*=4*αβ*−(1+*β*^{2})^{2}, and *x*=*vt* is the horizontal position of the crack tip as it moves in steady state. This same behaviour of *u*^{y} is present in the subsonic lattice solutions shown in figures 2*b* and 3*a*. We could use this square root displacement to extract a stress intensity factor *K*_{I} from the discrete system, as a complement to methods involving energy balance used in the past. However, as the crack speed increases towards the Rayleigh wave speed, the square root profile becomes increasingly indistinct, and it is completely lost at the Rayleigh speed and above.

As an illustration of the difference between subsonic and supersonic solutions, in figure 4 we plot the asymptotic behaviour of subsonic and supersonic solutions for atomic motion near the crack tip. For the subsonic solutions, the displacement rises as as expected. For supersonic solutions, it rises instead as .

### (b) Dependence upon system size

The subsonic lattice solutions have a vertical displacement that matches the continuum result equation (3.2), superposed with small-amplitude, high-frequency oscillations that result from periodic bond breaking, and carry off all the energy flux to the crack tip not absorbed by the bond-breaking process.

The supersonic solutions feature large-amplitude, low-frequency oscillations. There is a simple explanation for these oscillations, which take place at what we call the supersonic block frequency *ω*_{b}. This frequency is what we find if the lattice oscillates vertically as a block, with the top row held fixed. We can solve for this frequency using the equations of motion along the crack line equations (2.10) and (2.11). These oscillations occur after the crack has passed, so we can eliminate the *θ*(*t*) functions multiplying *U*^{−}(*t*). Fourier transforming these gives equations (2.13) and (2.14) with *U*^{−}(*ω*) replaced by zero. The components **u**_{3/2}(*ω*) can be written as linear combinations of the **u**_{1/2}(*ω*) components, which gives us a system to solve. The determinant of this system for **u**_{1/2}(*ω*) must vanish, which lets us solve for *ω*_{b}, the lowest normal mode frequency. The coefficients in equations (2.22) and (2.23) are so lengthy that this procedure is extremely cumbersome and not very informative.

It is more useful to consider the lattice as a continuous block and use the equations of linear elasticity [1]. The top of the block is fixed while the bottom is free to oscillate. If the relaxed length of the block is *L*, the lowest normal mode frequency can be shown to be
3.3
In our lattice model, the result becomes . This expression is approximate to the extent that our lattice system is not actually a continuum, and because the waves are not infinite in horizontal wavelength. Figure 5 shows how the vertical displacement along the crack line for a fixed supersonic speed varies with the system size. As the system grows from 402 rows high to 3202 rows high, the oscillation frequency decreases in proportion. A fit to the calculated results gives *ω*_{b}∼3.1742/*N*. This result agrees with the result in equation (3.3) within 5%.

We have seen how the vertical displacement along the crack line changes as the crack speed becomes supersonic for a fixed system size (figure 3). The subsonic solutions just below the Rayleigh wave speed have both small-amplitude, high-frequency oscillations as well as large-amplitude, low-frequency oscillations. The low-frequency oscillations as seen in figure 3*h* begin to appear in figure 3*c*,*d* alongside the high-frequency ones, which disappear completely at the Rayleigh speed (figure 3*f*) and above.

In figure 6, we examine the vertical motion of atoms along the crack line for a fixed subsonic speed as the system size varies. For the small system with 202 rows (*N*=100), the solution shares the character of both subsonic and supersonic solutions. Phonon oscillations are visible, but at the same time there is a long-wavelength oscillation with the frequency of the supersonic block frequency *ω*_{b}. Now we increase the system size and monitor the location of the first long-wavelength peak. As one can see in figure 6 it slides to the right and diminishes in amplitude. The location of the peak is given approximately by *t*_{b}∼0.7150*N*, just as the period of the large-amplitude oscillations goes as *T*∝*N* in figure 5. The height of the peak decreases roughly as and approaches the boundary condition *U*_{N} asymptotically from above as seen in figure 6*d*. Thus for any given subsonic crack speed below the Rayleigh wave speed, there exists a system size sufficiently large to produce behaviour like the continuum solution equation (3.2), and this makes it possible to find a stress intensity factor.

## 4. Conclusion

We have presented the full analytical solution for atomic motions accompanying steady-state in-plane crack motion in a brittle triangular lattice. This solution allows rapid numerical evaluation of atomic motion for very large systems, which is necessary to make comparisons with results from continuum fracture mechanics.

Equation (1.1) has been interpreted to mean that cracks cannot propagate faster than the Rayleigh wave speed. However, supersonic solutions do exist in lattice systems. Subsonic solutions are characterized by small-amplitude, high-frequency oscillations in the time dependence of the vertical displacement of an atom along the crack line. The overall envelope of the displacement agrees with the continuum result in equation (3.2). Supersonic solutions are characterized by large-amplitude, low-frequency oscillations, which can be attributed to the lattice oscillating as a block after the crack has passed. We call this frequency the supersonic block frequency *ω*_{b} and show that its value can be approximated using linear elasticity. We have also shown that subsonic solutions close to the Rayleigh wave speed exhibit small-amplitude, high-frequency and large-amplitude, low-frequency oscillations at the same time. This mixing of the character of subsonic and supersonic solutions can be eliminated by increasing the system size *N* sufficiently.

There are three natural objections to what we have presented. The first is that systems 2000 atoms high are computationally challenging, but they are a far cry from macroscopic systems 10^{12} atoms high. Thus, we have not really reached the macroscopic limit. The second is that we are drawing general conclusions about crack behaviour, but we have solved only a particular atomic model with very special and particularly simple interatomic forces. The third is that since the experiments of Schardin [26] it has been known that cracks do not reach the Rayleigh wave speed, let alone exceed it, so discussions of supersonic cracks are not relevant.

In response to the first objection, we have focused on scaling our solutions in such a manner that the large *N* behaviour becomes apparent. For example, in figure 6, we assert that in a system with *N*=10^{12} the plot of vertical displacement versus time for an atom just above the crack tip would look indistinguishable from figure 6*d*. In figure 5, we assert that, in a system with *N*=10^{12}, the period of the large oscillations would continue to increase in proportion to *N*, and in order for the oscillations to remain visible the dissipation would need to decrease accordingly.

In response to the second objection, we invoke the universality of fracture mechanics. Although our atomic force laws are much simpler than any forces that genuinely would arise from the quantum mechanics of interactions, they capture the essence of brittle solids, and include all of dynamic linear elastic fracture mechanics as a special case in the long-wavelength limit. Thus, we believe that conclusions we draw about the continuum (large *N*) limit of our theory should have general validity. The details of what happens on small scales should be indicative of what would happen in real brittle systems, but would ultimately depend upon realistic microscopic details. For examples of what more realistic systems might look like, see [24].

In response to the third objection, we note that cracks in brittle isotropic materials are usually limited by dynamic instabilities to terminal velocities below the Rayleigh wave speed [25]. However, in some cases, the instabilities are suppressed and supersonic cracks become possible. Rubber provides one instance [11], and materials with highly anisotropic fracture energy might provide others. Supersonic cracks are not inevitable just because instabilities have been suppressed [27]: sufficient energy must be supplied to drive cracks above the Rayleigh wave speed. But they do become possible.

The calculations presented here are just the first application of this method. There are some additional questions still pending. Because the stress intensity factor stops being defined as the crack passes the Rayleigh wave speed, the conventional continuum view of energy flux stops applying. This observation provides an unsatisfying answer to the question of why the crack speed is not limited by energy flux. Once the crack travels faster than the Rayleigh wave speed, it is clear that it must extract the energy needed for bond breaking from a local environment of finite size.

Thus, we intend to investigate the energy flux described by (1.1) more closely from an atomic point of view, and ask whether the seemingly pathological behaviour it predicts above the Rayleigh wave speed can be given any physical interpretation, or whether there exists some other quantity that describes energy flux and can be employed in both the subsonic and supersonic regimes.

## Funding statement

Support from the National Science Foundation Condensed Matter and Materials Theory Program (DMR 1002428) is gratefully acknowledged.

## Footnotes

One contribution of 19 to a theme issue ‘Fracturing across the multi-scales of diverse materials’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.