text:
# Adsorption of benzene on Si(100) from first principles

\[

## Abstract

Adsorption of benzene on the Si(100) surface is studied from first principles. We find that the most stable configuration is a tetra-$`\sigma `$-bonded structure characterized by one C-C double bond and four C-Si bonds. A similar structure, obtained by rotating the benzene molecule by 90, lies slightly higher in energy. However, rather narrow wells on the potential energy surface characterize these adsorption configurations. A benzene molecule impinging on the Si surface is most likely to be adsorbed in one of three different di-$`\sigma `$-bonded, metastable structures, characterized by two C-Si bonds, and eventually converts into the lowest-energy configurations. These results are consistent with recent experiments.

\]

Adsorption of benzene on the Si(100) surface is a topic of great current interest  both because it represents a prototype system for the study of molecular adsorption (and desorption) of hydrocarbons on semiconductor surfaces, and because it is considered a promising precursor for technologically relevant processes, such as the growth of Si-C and CVD diamond thin films on Si surfaces. However, despite many experimental and theoretical investigations, the adsorption mechanism is not yet well understood. In particular, at present there is no consensus about the lowest-energy structure of benzene on Si(100): results obtained from surface science experimental techniques, semiempirical methods, and first-principles approaches provide a number of different predictions.

Benzene is known from experiments to adsorb exclusively on top of the Si(100) surface dimer rows, thus avoiding energetically disfavored structures with unsaturated, isolated Si dangling bonds. Even so, since the size of the benzene molecule is comparable to the spacing between two adjacent dimers on the same row, many different bonding configurations are possible. Among the structures proposed in the literature as the lowest-energy configurations, the 1,4-cyclohexadiene-like (“butterfly”) configuration, in which the benzene molecule is di-$`\sigma `$ bonded to the two dangling bonds of the same Si surface dimer, is supported by thermal desorption and angle-resolved photoelectron spectroscopy, STM, vibrational IR spectroscopy and near-edge X-ray absorption fine structure techniques, and first-principles cluster calculations. Instead, other STM experiments suggest the 1,3-cyclohexadiene-like (“tilted”) structure. Finally, semiempirical calculations, STM and IR spectroscopy experiments  favor a tetra-$`\sigma `$-bonded configuration where benzene is bonded to two adjacent surface dimers.

Another open issue concerns the occurrence and nature of metastable adsorption states. In fact, the results of STM and IR spectroscopy support the hypothesis that benzene is initially chemisorbed in a metastable, “butterfly”-like state, and then slowly converts (within minutes) to a lower-energy final state, which is a “tilted” structure according to Ref. , or a tetra-$`\sigma `$-bonded one according to Ref. . Moreover, recent IR experiments  suggest that, at room temperature, benzene is predominantly adsorbed in the “butterfly” configuration, while the existence of a less stable structure, consistent with a tetra-$`\sigma `$-bonded configuration, is proposed.

Previous theoretical calculations on benzene on Si(100) have been restricted to semiempirical or ab initio cluster-model methods. In the latter approach the Si surface is modeled with a cluster of Si atoms, thus considerably reducing the cost of a first-principles calculation. However, the effects of such an approximation can be relevant. It is well known, for instance, that the characteristic buckling of the Si dimers on the clean Si(100) surface can only be obtained by using models with a slab geometry and periodic boundary conditions. As shown in the following, the details of the surface reconstruction (i.e. buckling and periodicity of the surface dimers) are crucial ingredients in determining the adsorption structure of benzene. Moreover, the convergence of different properties, such as the binding energies of adsorbed molecules, is rather slow as a function of the cluster size.

In order to overcome these limitations and to clarify the open issues discussed above we have performed a full ab initio study of benzene adsorption on Si(100). Total-energy calculations and molecular dynamics (MD) simulations have been carried out within the Car-Parrinello approach in the framework of the density functional theory, in the local spin density approximation. Tests have been also performed using gradient corrections in the BLYP implementation. The calculations have been carried out considering the $`\mathrm{\Gamma }`$-point only of the Brillouin zone, and using norm-conserving pseudopotentials, with $`s`$ and $`p`$ nonlocality for C and Si. Wavefunctions were expanded in plane waves with an energy cutoff of 35 Ry.

The Si(100) surface is modeled with a periodically repeated slab of 5 Si layers and a vacuum region of 7 Å (tests have been also carried out with a vacuum region of 10 Å , without any significant change in the results). A monolayer of hydrogen atoms is used to saturate the dangling bonds on the lower surface of the slab. We have used a supercell with $`p(\sqrt{8}\times \sqrt{8})R45^{}`$ surface periodicity, corresponding to 8 Si atoms/layer; however, in order to check finite-size effects, the geometry optimizations have been repeated using a larger $`p(4\times 4)`$ supercell with 16 atoms/layer.

Structural relaxations of the ionic coordinates are performed using the method of direct inversion in the iterative subspace . During ionic relaxations and MD simulations the lowest Si layer and the saturation hydrogens are kept fixed. We verified that, by starting with the unreconstructed, clean Si(100) surface, the structural optimization procedure correctly produces asymmetric surface dimers, with a dimer bond length and buckling angle in good agreement with previous, highly converged ab initio calculations . We have considered different surface periodicities for the dimer reconstruction which may occur on the Si(100)surface, i.e. $`(2\times 1)`$, $`p(2\times 2)`$ and $`c(4\times 2)`$. A single benzene molecule is added on top of the slab and the system is then fully relaxed towards the minimum energy configuration. To better explore the complex potential energy surface of this system, in most of the cases the optimization procedure was repeated using a simulated-annealing strategy and also starting from different initial configurations.

We find that the lowest-energy configurations are given by two tetra-$`\sigma `$-bonded structures, characterized by the presence of one C-C double bond, which we refer to as “tight bridge” (TiB) and “twisted bridge” (TwB) (see Fig. 1). TwB is similar to TiB but the benzene molecule is rotated by 90 with respect to the Si surface and is slightly higher in energy (see Table I). This result is in agreement with the findings of Ref.  and turns out to be independent on the size of the supercell used in the simulation and on the different reconstructions of the Si(100) surface. It remains true also using BLYP gradient corrections, as can be seen in Table I.

We also find, at somewhat higher energies, three different, metastable “butterfly” structures, characterized by two C-Si bonds, which are shown in Fig. 1. One of them (“standard butterfly”, SB) is the well-known configuration with the benzene molecule adsorbed on top of a single Si dimer. The others (“tilted-bridge butterfly”, TB, and “diagonal-bridge butterfly”, DB), which bridge two adjacent surface dimers, have not been reported in any previous study.

The Si(100) reconstruction crucially affects the occurrence and energetic ordering of the three “butterfly” structures. In fact, in the $`(2\times 1)`$ reconstruction (with parallel buckled dimers), SB and TB are the most stable (almost isoenergetic) “butterfly” configurations, while DB is considerably less favored; in contrast, with reconstructions involving alternating buckled Si dimers, such as the $`p(2\times 2)`$ and the $`c(4\times 2)`$, SB and DB are the lowest-energy configurations, while the binding energy of TB is significantly smaller. This clearly happens because the two C-Si bonds of the TB structure are more easily created when the benzene molecule is adsorbed onto Si(100) $`(2\times 1)`$, while the formation of the DB structure is favored by the presence of alternating buckled Si dimers.

The other configurations proposed in the literature, that is the “tilted” (T) and the “pedestal” (P) ones, lie higher in energy for all the Si(100) reconstructions considered (see Table I). In particular, the P structure is only found to be stable in the $`(2\times 1)`$ reconstruction; however, even in this case, a MD simulation performed at 300 K shows that the structure converts very rapidly (in less than 1 ps) into a DB structure. Although the P structure has four C-Si bonds, it is energetically disfavored because it involves the presence of two radical centers.

Inspection of the C-C distances for the various stable structures reveals the existence of two kind of bonds: a long one (”single”) and a short one (”double”), of length $`1.491.59`$ and $`1.341.36`$ Å, respectively. These values should be compared with the C-C bond length in the isolated benzene molecule, 1.39 Å. One double bond characterizes the TiB and TwB structures, while two double bonds are found in the “butterfly” structures. In contrast, in the P configuration all the C-C bonds are single ones. These conclusions are confirmed by a more quantitative analysis of the electronic orbitals, which we performed by using both the notion of Mayer bond order and the method of the Localized Wannier functions . In the three “butterfly” configurations (SB, TB, DB), the bond angles ($`119^{}`$-$`122^{}`$) at the C atoms not involved in the Si-C bonds are close to that ($`120^{}`$) of the isolated benzene molecule, while those ($`103^{}`$-$`113^{}`$) at the 4-fold coordinated C atom are closer to the ideal tetrahedron ($`109.5^{}`$) angle. This clearly indicates $`sp^2`$ and $`sp^3`$ hybridization, respectively. After benzene chemisorption, although the Si-Si dimers are preserved, the Si dimer buckling angle is almost reduced to zero, with the exception of the TB and DB structures. In the lowest-energy TiB structure the angle between the double bond and the Si(100) surface is $`45^{}`$ in good agreement with the experimental estimate , $`43^{}`$.

The structural parameters do not change appreciably when a larger $`p(4\times 4)`$ surface supercell is used. Use of BLYP gradient corrections makes bond lengths about 1-2 % longer, while binding energies are significantly reduced (see Table I). Moreover, in the $`(2\times 1)`$ reconstruction, the P configuration is no longer stable and, among the three “butterfly” structures, BLYP favors SB, while the binding energy of DB is even smaller than that of the T structure. Note, however, that TiB and TwB remain the lowest-energy configurations.

According to the results of some experiments and theoretical calculations, adsorbed benzene predominantly forms a “butterfly” (SB) configuration, while the TiB one (and perhaps TwB) appears in detectable amounts on relatively long timescales only, thus indicating the existence of an energy barrier between the two structures.

In order to identify possible metastable states, occurring in the early stages of adsorption, we have tried to find, in the simplest way, the most probable structure of a benzene molecule impinging on the Si(100) surface. If we place the molecule at some distance from the surface we observe that, regardless of the initial position and orientation of the molecule, after full relaxation the final structure is almost invariably one of the three “butterfly” configurations. This happens because the dimers are tilted, favoring the formation of the di-bonded “butterfly” structures rather than the tetra-bonded ones. The specific “butterfly” configuration which is actually formed depends critically on the type of reconstruction of the Si surface that is considered, as already discussed above. On the contrary, there are only very few initial positions which lead to the low-energy TiB and TwB configurations.

We have tried to characterize the energy barrier which must be overcome to relax from the “butterfly” configurations to the lower-energy TiB and TwB structures. To this aim we started with the benzene molecule in the SB configuration. Let Cd be one of the C atoms involved in the Si-C bonds. Many calculations have been performed in which the ionic coordinates of both the molecule and the substrate were optimized under the constraint that the $`x,y`$ coordinates of the two Cd atoms are held fixed. A particular pathway, connecting the SB to the DB structure, is shown in Fig. 2, where the reaction coordinate is defined as the distance between the Cd-Cd axis of the initial configuration and that of the displaced structure. The pronounced energy minimum corresponds to the occurrence, during the transformation, of the lowest-energy Tib structure. Note however that this is characterized by a very narrow well. From Fig. 2 a lower bound of $`0.5`$ eV can be inferred for the energy barrier, to be compared with the experimental estimates, $``$ 0.9-1.0 eV. A similar calculation for the TB$``$TwB transition gives a smaller value of $`0.4`$ eV. As a consequence the conversion from TB to TwB is expected to be somewhat faster than that from SB to TiB.

A large fraction of experiments on benzene on Si(100) is based on STM techniques. However, different interpretations of similar STM images led to contradictory conclusions about the adsorption sites and geometry of the adsorbed molecules. For each of the structures reported in Table I we have produced “theoretical” STM images to be compared with the experimental ones, following the recipe of Ref. . Charge density iso-surfaces have been obtained by including electron states in an energy range down to $`2`$ eV below the highest occupied state, which corresponds to typical STM bias voltages. The simulated images are obtained by viewing these iso-surfaces at typical tip-surface distances (a few Å above the benzene molecule).

Our computed STM image for the TiB structure exhibits a density maximum above one of the two Si dimers involved in bonding with benzene, while the TwB configuration produces a similar image but rotated by 90. These images resemble those obtained by Lopinski et al.. The theoretical STM image for the SB structure is characterized by a bright two-lobe protrusion centered symmetrically above a single Si dimer unit and oriented orthogonal to the dimer axis, in qualitative agreement with the experimental findings. Instead, the STM images of the TB and DB structures are quite different from that of SB. In fact the TB image is qualitatively similar to that of TwB (and the experimental STM resolution could be insufficient to distinguish between the two configurations), while DB gives rise to a much fainter feature, bridging in diagonal two Si dimers, which is probably hardly visible in experiments. These observations could explain why the DB and TB structures have not been detected in STM experiments. The T configuration produces an asymmetric (with respect to Si dimers) image, appearing as a bright region (placed between two Si dimers) adjacent to a dark region. Finally the P structure is characterized by two spots corresponding to the dangling bonds of benzene; this result supports the conjecture  which rules out the presence of a significant fraction of benzene molecules adsorbed in the P structure because of the absence of such spots in the STM images.

We have also computed the vibrational spectra for a representative “butterfly” structure, SB, and for the lowest-energy TiB configuration, by performing Car-Parrinello MD simulations at room temperature. Our results for TiB show a slightly more quantitative agreement with the experimental results  than those for the SB structure, although the main features of the spectra are similar in the two structures. Let C$`^{\prime \prime }`$ (C$`^{}`$) denote a C atom which shares a double (single) bond with another C atom. The C$`^{}`$-H and C$`^{\prime \prime }`$-H frequencies (2880 and 3010 cm-1) are in agreement with the $`sp^3`$ and $`sp^2`$ stretching modes observed in recent IR spectroscopy experiments  (2945 and 3044 cm-1), and semiempirical cluster calculations . Note that the C-H vibrations for the isolated benzene molecule are characterized by a single frequency of 3040 cm-1. For the C$`^{}`$-C$`^{\prime \prime }`$ and C$`^{\prime \prime }`$-C$`^{\prime \prime }`$ frequencies we find 1230 and 1520 cm-1, respectively, to be compared with the EELS experimental values , 1170 and 1625 cm-1. The C-H bending modes are found at 900 and 1100 cm-1, whereas experimentally  they are at 910 and 1075 cm-1.

In conclusion, using state-of-the-art ab initio simulations, we have shown that a tetra-$`\sigma `$ bonded structure is the most stable configuration for benzene adsorbed on Si(100). However, this structure and a very similar one, lying only slightly higher in energy, correspond to very narrow wells in the potential energy surface for a benzene molecule impinging on the surface. Therefore it is more likely for the molecule to be adsorbed into one of three different, metastable “butterfly” configurations, and eventually convert into the lowest-energy structures. Our study provides detailed information about structural, electronic, and vibrational properties of the system, and allows a critical comparison with results obtained from different experimental techniques and previous theoretical calculations.

We thank M. Boero and A. Vittadini for useful discussions. This work is partially supported by INFM through the Parallel Computing Initiative.


format: text

text:
# Steady-state mode I cracks in a viscoelastic triangular lattice

## I Introduction

In the previous paper (Pechenik, et al., 2000), we showed how one could extend the methodology originally devised by Slepyan (1981) and obtain a closed form solution for mode III cracks propagating steadily in a square lattice in the presence of dissipation in the form of a Kelvin viscosity. This approach allowed us to determine the effect of dissipation on the crack velocity as a function of driving. Similarly, we found the critical velocity (and the associated critical time) at which other bonds in the lattice have displacements larger than the assumed breaking criterion, signaling an instability in the steadily-propagating crack and the onset of more complex spatio-temporal behavior.

In this second paper, we extend this analysis to the more interesting (and more experimentally relevant) case of mode I cracks. We now work on a triangular lattice and assume that there are central force ideally brittle viscoelastic springs connecting particles at the lattice sites. Again, we use ideas borrowed from Kulamekhtova (1984) to solve this problem exactly on an infinite lattice. We also solve the same system numerically on a lattice of finite transverse extent so as to provide an independent check on some of the results.

Perhaps the most interesting finding reported here concerns the aforementioned onset of additional bond-breaking at large enough velocity. At small dissipation, diagonal bonds that are offset from the assumed crack line in the vertical direction are the most “dangerous”. The critical velocity at which one of these bonds exceeds critical displacement is more or less independent of the dissipation at low values thereof and eventually rises sharply as the system becomes increasingly damped. This is similar to the behavior obtained for mode III. At moderate to large damping, however, horizontal bonds become more relevant and eventually give rise to a critical velocity which decreases with increasing damping. At the end, we will comment on the possible relevance of our results for trying to make sense of recent experiments (Fineberg, et al., 1991, 1992; Sharon, et al., 1995, 1996) and simulations on instabilities during dynamic fracture.

## II Mode I on a triangular lattice

Our model consists of central force springs between mass points located on a triangular lattice. We introduce unit vectors in the lattice directions as follows: $`\widehat{d}_1`$ points in the direction of the $`x`$-axis and each of the subsequent unit vectors $`\widehat{d}_i`$ with $`2i6`$ rotated by an additional $`\pi /3`$ in the counter-clockwise direction. This means that $`\widehat{d}_{i+3}=\widehat{d}_i`$. We will always assume that the steady-state crack breaks bonds between the rows with $`y=0`$ and $`y=\sqrt{3}/2`$, as it moves in the direction of the $`x`$-axis.

As already discussed in our previous paper on mode III cracks, each bond is taken to be linear until a critical displacement at which point it completely breaks. Also, dissipation is incorporated as a Kelvin viscosity term. Given these assumptions, the equation describing the motion of the masses in the lattice is

$$(1+\eta \frac{}{t})\underset{i=1}{\overset{6}{}}Q_i(t,\stackrel{}{x})\widehat{d}_i\frac{^2\stackrel{}{u}(t,\stackrel{}{x})}{t^2}=\stackrel{}{\sigma }(t,\stackrel{}{x}).$$
(1)
Here, $`\stackrel{}{x}=(x,y)`$ is a lattice vector in the plane and $`Q_i(t,\stackrel{}{x})`$ is the elongation of bond $`i`$ emanating from the specific lattice site, in the approximation that this elongation is much smaller than the unstretched bond length,

$$Q_i(t,\stackrel{}{x})=\left(\stackrel{}{u}(t,\stackrel{}{x}+\widehat{d}_i)\stackrel{}{u}(t,\stackrel{}{x})\right)\widehat{d}_i.$$
(2)
Finally, the driving term $`\stackrel{}{\sigma }(t,\stackrel{}{x})`$ consists of forces which compensate for the forces from broken bonds (which are included on the left-hand side) as well as any external forces acting on the crack boundaries.

Next we suppose that the crack is moving with a constant speed $`v`$. This allows us to change variables to $`\tau =xvt`$, $`\stackrel{}{z}=(\tau ,y)`$, which gives

$$(1\eta v\frac{}{\tau })\underset{i=1}{\overset{6}{}}Q_i(\stackrel{}{z})\widehat{d}_iv^2\frac{^2\stackrel{}{u}(\stackrel{}{z})}{\tau ^2}=\stackrel{}{\sigma }(\stackrel{}{z}).$$
(3)
Because of the form of $`Q_i(\stackrel{}{z})`$, we have the obvious symmetries

$$Q_2(\stackrel{}{z})=Q_5(\stackrel{}{z}+\widehat{d}_2),Q_3(\stackrel{}{z})=Q_6(\stackrel{}{z}+\widehat{d}_3).$$
(4)
We will impose one additional symmetry condition on the $`Q_i`$. As the crack propagates, it alternately breaks bonds in the $`\widehat{d}_2`$ and $`\widehat{d}_3`$ directions. Let us define our coordinates such that at $`\tau =0`$ and $`t=0`$, (so that $`x=0`$), the crack breaks a bond in the direction of $`\widehat{d}_2`$. This means that for $`\tau >0`$ this bond ($`Q_2(\stackrel{}{z}_0)`$ with $`\stackrel{}{z}_0=(\tau ,0)`$) is unbroken and for $`\tau <0`$ it is always broken. If we shift one lattice spacing to the left, the $`\widehat{d}_2`$ bond at $`x=1`$ will break at $`t=1/v`$ (here also $`\tau =0`$ in $`Q_2(\stackrel{}{z}_0)`$). It is clear that the bond in the $`\widehat{d}_3`$ direction at $`x=0`$ breaks at a time-point between these two events; we will assume that this time is exactly at the midpoint, $`t=1/(2v)`$ or equivalently $`\tau =1/2`$ in $`Q_3(\stackrel{}{z}_0)`$. This means that for $`\tau >1/2`$ this bond is unbroken and for $`\tau <1/2`$ it is always broken. We will assume an even more strict symmetry condition for these bonds, namely

$$Q_2(\stackrel{}{z}_0)=Q_3(xv(t\frac{1}{2v}),y=0)=Q_3(\stackrel{}{z}_0+\frac{1}{2}\widehat{d}_1)$$
(5)
which will need to be checked from our solution later.

Let us denote $`Q_2(\stackrel{}{z}_0)Q(\tau )`$. Then, we can write down an explicit expression for the forces $`\stackrel{}{\sigma }(z)`$ on the right hand side of Eq. 3. For the row with $`y=0`$

$`\stackrel{}{\sigma }(\stackrel{}{z}_0)`$ $`=`$ $`\theta (\tau )\left(\stackrel{}{\sigma }_e(\stackrel{}{z}_0)(1\eta v{\displaystyle \frac{}{\tau }})Q_2(\stackrel{}{z}_0)\widehat{d}_2\right)`$
$`+\theta (\tau +{\displaystyle \frac{1}{2}})\left(\stackrel{}{\sigma }_e^{}(\stackrel{}{z}_0)(1\eta v{\displaystyle \frac{}{\tau }})Q_3(\stackrel{}{z}_0)\widehat{d}_3\right)`$
$`=`$ $`\theta (\tau )\left(P_0(\tau )(1\eta v{\displaystyle \frac{}{\tau }})Q(\tau )\right)\widehat{d}_2`$
$`+\theta (\tau +{\displaystyle \frac{1}{2}})\left(P_0(\tau {\displaystyle \frac{1}{2}})(1\eta v{\displaystyle \frac{}{\tau }})Q(\tau {\displaystyle \frac{1}{2}})\right)\widehat{d}_3.`$ (7)
Here $`\stackrel{}{\sigma }_e`$ and $`\stackrel{}{\sigma }_e^{}`$ are the assumed external forces. In the second form, these have been taken to match the vectorial nature and also the symmetry of the bond-breaking terms. The (scalar) $`P_0`$ will be specified later. Note that the second theta-function in (II) reflects the already mentioned fact that the bond from point $`\stackrel{}{z}_0`$ in the $`\widehat{d}_3`$ direction breaks earlier (by time interval $`1/2v`$) than the bond from the same point in the $`\widehat{d}_2`$ direction.

Let us write down the corresponding force for the row $`y=\sqrt{3}/2`$. For the point $`\stackrel{}{x}=\widehat{d}_2`$, the bond in the $`\widehat{d}_5`$ direction is physically the same as the $`\widehat{d}_2`$ bond at $`\stackrel{}{x}=(0,0)`$ and therefore breaks at $`\tau =0`$ in $`Q_5(\stackrel{}{z}_0+\widehat{d}_2)=Q_2(\stackrel{}{z}_0)`$. Now, the $`\widehat{d}_6`$ bond at $`\stackrel{}{x}=\widehat{d_2}`$ breaks later, at $`\tau =1/2`$ in $`Q_6(\stackrel{}{z}_0+\widehat{d}_2)`$. Thus we have

$`\stackrel{}{\sigma }(\stackrel{}{z}_0+\widehat{d}_2)`$ $`=`$ $`\theta (\tau )\left(\stackrel{}{\sigma }_e(\stackrel{}{z}_0+\widehat{d}_2)(1\eta v{\displaystyle \frac{}{\tau }})Q_5(\stackrel{}{z}_0+\widehat{d}_2)\widehat{d}_5\right)`$
$`+\theta (\tau {\displaystyle \frac{1}{2}})\left(\stackrel{}{\sigma }_e^{}(\stackrel{}{z}_0+\widehat{d}_2)(1\eta v{\displaystyle \frac{}{\tau }})Q_6(\stackrel{}{z}_0+\widehat{d}_2)\widehat{d}_6\right)`$
$`=`$ $`\theta (\tau )\left(P_0(\tau )(1\eta v{\displaystyle \frac{}{\tau }})Q(\tau )\right)\widehat{d}_5+`$
$`+\theta (\tau {\displaystyle \frac{1}{2}})\left(P_0(\tau +{\displaystyle \frac{1}{2}})(1\eta v{\displaystyle \frac{}{\tau }})Q(\tau +{\displaystyle \frac{1}{2}})\right)\widehat{d}_6.`$ (9)
where we used the symmetry in Eq. 5 to write $`Q_6(\stackrel{}{z}_0+\widehat{d}_2)=Q_3(\stackrel{}{z}_0+\widehat{d}_2\widehat{d}_3)=Q_2(\stackrel{}{z}_0+\widehat{d}_2\widehat{d}_3\widehat{d}_1/2)=Q(\tau +1/2)`$. Finally, we can rewrite all these forces in a more compact form,

$`\stackrel{}{\sigma }(\stackrel{}{z}_0)`$ $`=`$ $`N(\tau )\widehat{d}_5+N(\tau {\displaystyle \frac{1}{2}})\widehat{d}_6,`$ (10)
$`\stackrel{}{\sigma }(\stackrel{}{z}_0+\widehat{d}_2)`$ $`=`$ $`N(\tau )\widehat{d}_2+N(\tau +{\displaystyle \frac{1}{2}})\widehat{d}_3,`$ (11)
where

$$N(\tau )=\theta (\tau )\left(P_0(\tau )+(1\eta v\frac{}{\tau })Q(\tau )\right).$$
(12)
Next, we proceed by doing a Fourier transformation of Eq. (3), over the continuous variable $`\tau `$ and the discrete variable $`y=n\sqrt{3}/2`$, with integer $`n`$. In detail,

$`u^F(\tau ,s)`$ $`=`$ $`{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}u(\tau ,{\displaystyle \frac{\sqrt{3}}{2}}n)e^{i\frac{\sqrt{3}}{2}sn},`$ (13a)
$`u(\tau ,y)`$ $`=`$ $`{\displaystyle \frac{\sqrt{3}}{4\pi }}{\displaystyle _{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}}u^F(\tau ,s)e^{isy}𝑑s;`$ (13b)
$`u^{FF}(q,s)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^+\mathrm{}}u^F(\tau ,s)e^{iq\tau }𝑑\tau ,`$ (14a)
$`u^F(\tau ,s)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}u^{FF}(q,s)e^{iq\tau }𝑑q.`$ (14b)
We thereby obtain

$$v^2q^2\stackrel{}{u}^{FF}(1+i\eta vq)\underset{i}{}Q_i^{FF}\widehat{d}_i=\stackrel{}{\sigma }^{FF},$$
(15)
with

$$Q_i^{FF}=(e^{i\stackrel{}{r}\widehat{d}_i}1)\stackrel{}{u}^{FF}\widehat{d}_i$$
(16)
and $`\stackrel{}{r}=(q,s)`$.

To find $`\stackrel{}{\sigma }^{FF}`$ we first transform Eqs. (10, 11) over $`\tau `$

$`\stackrel{}{\sigma }^F(q,n=0)`$ $`=`$ $`N(q)\widehat{d}_5+e^{i\frac{q}{2}}N(q)\widehat{d}_6,`$ (17a)
$`\stackrel{}{\sigma }^F(q,n=1)e^{i\frac{q}{2}}`$ $`=`$ $`N(q)\widehat{d}_2+e^{i\frac{q}{2}}N(q)\widehat{d}_3,`$ (17b)
with

$$N(q)=P_0^{}+(1+i\eta vq)Q^{}\eta vQ(0).$$
(18)
The superscript “$``$” refers to the part of the Fourier transform which is analytic in the lower half plane of the variable $`q`$; our notation follows that in our previous mode III paper. Finally, using (13a) we get

$$\stackrel{}{\sigma }^{FF}=N(q)\widehat{d}_5+e^{i\frac{q}{2}}N(q)\widehat{d}_6+\left(e^{i\frac{q}{2}}N(q)\widehat{d}_2+N(q)\widehat{d}_3\right)e^{is\frac{\sqrt{3}}{2}}.$$
(19)
This then completes the derivation of the basic equation that needs to be solved to determine the steady-state crack field.

## III Wiener-Hopf Solution

We wish to use the Wiener-Hopf technique to solve the equation derived above for the elongation field $`Q`$. To start, we note that Eq. (15) can be written in matrix form

$$𝐌\stackrel{}{u}^{FF}=\stackrel{}{\sigma }^{FF}$$
(20)
with the 2x2 matrix

$$𝐌=\left(\begin{array}{cccc}𝐀_\mathrm{𝟏}& 𝐀_\mathrm{𝟑}\hfill & & \\ 𝐀_\mathrm{𝟑}& 𝐀_\mathrm{𝟐}\hfill & & \end{array}\right)$$
(21)
where

$`A_1`$ $`=`$ $`v^2q^2+(1+i\eta vq)\left(4\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_1}{2}}+\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_2}{2}}+\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_3}{2}}\right),`$ (22a)
$`A_2`$ $`=`$ $`v^2q^2+3(1+i\eta vq)\left(\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_2}{2}}+\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_3}{2}}\right),`$ (22b)
$`A_3`$ $`=`$ $`\sqrt{3}(1+i\eta vq)\left(\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_2}{2}}\mathrm{sin}^2{\displaystyle \frac{\stackrel{}{r}\widehat{d}_3}{2}}\right),`$ (22c)
and

$$\stackrel{}{\sigma }^{FF}=\frac{N(q)}{2}\left(\begin{array}{cccc}(1e^{i\frac{q}{2}})(1+e^{i\frac{\sqrt{3}}{2}s})& & & \\ \sqrt{3}(1+e^{i\frac{q}{2}})(1e^{i\frac{\sqrt{3}}{2}s})& & & \end{array}\right).$$
(23)
The solution of (20) is therefore

$$\stackrel{}{u}^{FF}=\frac{1}{det𝐌}\left(\begin{array}{cccc}det𝐌_1& & & \\ det𝐌_2& & & \end{array}\right)$$
(24)
where

$$𝐌_1=\left(\begin{array}{cccc}\sigma _x^{FF}& \sigma _y^{FF}\hfill & & \\ A_3& A_2\hfill & & \end{array}\right),𝐌_2=\left(\begin{array}{cccc}A_1& A_3\hfill & & \\ \sigma _x^{FF}& \sigma _y^{FF}\hfill & & \end{array}\right).$$
(25)
Now, we want to extract from the solution (24) an equation for the variable we are interested in, namely $`Q^F(q)`$. From the definition of $`Q`$ we have

$$Q^F(q)=Q_2^F(q,n=0)=\frac{\sqrt{3}}{4\pi }_{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}(e^{i\stackrel{}{r}\widehat{d}_2}1)\stackrel{}{u}^{FF}\widehat{d}_2𝑑s.$$
(26)
From the basic solution above, we can evaluate $`\stackrel{}{u}^{FF}\widehat{d}_2`$.

$$\stackrel{}{u}^{FF}\widehat{d}_2=\frac{\left|\begin{array}{cccc}d_{2x}A_3+d_{2y}A_1& d_{2x}A_2+d_{2y}A_3\hfill & & \\ \sigma _x^{FF}& \sigma _y^{FF}\hfill & & \end{array}\right|}{|𝐌|}.$$
(27)
Before proceeding to evaluate the integral, it is worthwhile to pause and check the basic symmetry conditions. We find that

$$Q_3^F(q,n=0)=\frac{\sqrt{3}}{4\pi }_{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}(e^{i\stackrel{}{r}\widehat{d}_3}1)\stackrel{}{u}^{FF}\widehat{d}_3𝑑s,$$
(28)
where the new dot product is given by

$$\stackrel{}{u}^{FF}\widehat{d}_3=\frac{\left|\begin{array}{cccc}d_{2x}A_3+d_{2y}A_1& d_{2x}A_2+d_{2y}A_3\hfill & & \\ \sigma _x^{FF}& \sigma _y^{FF}\hfill & & \end{array}\right|}{|𝐌|}.$$
(29)
Now, we can change variables of integration in (28) from $`s`$ to $`s`$. Under this transformation $`A_1|_{ss}=A_1`$, $`A_2|_{ss}=A_2`$, $`A_3|_{ss}=A_3`$, $`\sigma _x^{FF}|_{ss}=e^{i\sqrt{3}s/2}\sigma _x^{FF}`$, $`\sigma _y^{FF}|_{ss}=e^{i\sqrt{3}s/2}\sigma _y^{FF}`$, which gives that $`det𝐌|_{ss}=det𝐌`$ and $`\stackrel{}{u}^{FF}\widehat{d}_3|_{ss}=e^{i\sqrt{3}s/2}\stackrel{}{u}^{FF}\widehat{d}_2`$. Therefore (28) changes into

$$Q_3^F(q,n=0)=\frac{\sqrt{3}}{4\pi }_{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}(e^{iqd_{2x}+isd_{2y}}1)e^{id_{2y}s}\stackrel{}{u}^{FF}\widehat{d}_2𝑑s=e^{i\frac{q}{2}}Q_2^F(q,n=0).$$
(30)
This is exactly the symmetry condition (5); hence, our solution is consistent with the assumed symmetry. We also note here that a similar derivation gives

$$Q_2(\tau ,y)=Q_3(\tau +\frac{1}{2},y)=Q_6(\tau +\frac{1}{2},y+\frac{\sqrt{3}}{2}).$$
(31)
More generally, we can see what these symmetries mean for the displacement field $`\stackrel{}{u}(z)`$. Using $`det𝐌_1|_{ss}=e^{i\sqrt{3}s/2}det𝐌_1`$, $`det𝐌_2|_{ss}=e^{i\sqrt{3}s/2}det𝐌_2`$, we obtain $`u_x^{FF}|_{ss}=e^{i\sqrt{3}s/2}u_x^{FF}`$ and $`u_y^{FF}|_{ss}=e^{i\sqrt{3}s/2}u_y^{FF}`$. These give $`u_x(\tau ,y)=u_x(\tau ,y+\sqrt{3}/2)`$ and $`u_y(\tau ,y)=u_y(\tau ,y+\sqrt{3}/2)`$.

We now proceed to calculate the integral in (26). The explicit form of the integrand in (26) is

$$\frac{N(q)}{|𝐌|}\left|\begin{array}{cccc}3\left(\frac{v^2q^2}{2}(1+i\eta vq)(2\mathrm{sin}^2\frac{q}{2}+\mathrm{sin}^2\frac{\stackrel{}{r}\widehat{d}_3}{2})\right)& \frac{v^2q^2}{2}+3(1+i\eta vq)\mathrm{sin}^2\frac{\stackrel{}{r}\widehat{d}_3}{2}\hfill & & \\ (\mathrm{cos}\stackrel{}{r}\widehat{d}_21)(\mathrm{cos}\frac{\sqrt{3}}{2}s\mathrm{cos}\frac{q}{2})& (\mathrm{cos}\stackrel{}{r}\widehat{d}_21)+(\mathrm{cos}\frac{\sqrt{3}}{2}s\mathrm{cos}\frac{q}{2})\hfill & & \end{array}\right|.$$
(32)
As $`|𝐌|`$ is an even function of $`s`$, only the even part of the other determinant contributes to the integral. Then, by the change of variable to $`w=e^{i\sqrt{3}s/2}`$ we can rewrite (26) as an integral over the unit circle in the complex $`w`$ plane. After some tedious algebra, we find that we can write

$$det𝐌=3\xi ^2+\alpha \xi +\beta $$
(33)
where

$`\xi `$ $`=`$ $`(1+i\eta vq){\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{w}}+w\right),`$ (34a)
$`\alpha `$ $`=`$ $`2\mathrm{cos}{\displaystyle \frac{q}{2}}\left(2v^2q^2+3(1+i\eta vq)(1+2\mathrm{sin}^2{\displaystyle \frac{q}{2}})\right),`$ (34b)
$`\beta `$ $`=`$ $`3(1+i\eta vq)^2(1+3\mathrm{sin}^2{\displaystyle \frac{q}{2}})4(1+i\eta vq)v^2q^2(1+\mathrm{sin}^2{\displaystyle \frac{q}{2}})+v^4q^4.`$ (34c)
Similarly, the even part of the determinant in the numerator becomes

$$\frac{F(\xi )}{1+i\eta vq},$$
(35)
with

$`F(\xi )`$ $`=`$ $`3\left[(\xi (1+iq\eta v)\mathrm{cos}{\displaystyle \frac{q}{2}})^2+2(1+iq\eta v)\mathrm{sin}^2{\displaystyle \frac{q}{2}}(1+\mathrm{cos}{\displaystyle \frac{q}{2}})(1+iq\eta v\xi )\right]`$ (36)
$`v^2q^2\left[(1+iq\eta v\xi )(1+\mathrm{cos}{\displaystyle \frac{q}{2}})\xi \mathrm{cos}{\displaystyle \frac{q}{2}}+1+i\eta vq\right].`$
We write down here also the odd part of this determinant for future reference,

$$G(w)=i(\frac{1}{w}w)\mathrm{sin}\frac{q}{2}(v^2q^23(1+iq\eta v)\mathrm{sin}^2\frac{q}{2}).$$
(37)
Inserting (33) and (35) into (32) we can see that (26) becomes

$$Q^F(q)=\frac{N(q)}{2\pi i(1+i\eta vq)}\underset{|w|=1}{}\frac{F(\xi )}{3\xi ^2+\alpha \xi +\beta }\frac{dw}{w}$$
(38)
with integration around the unit circle in the counter-clockwise direction. The integrand in Eq. (38) has three poles inside the unit circle; i.e., $`w=0`$ and one $`w`$ for each of the two roots of the expression (33),

$$\xi _{1,2}=C\pm \sqrt{C^2D},$$
(39)
with

$`C`$ $`=`$ $`\mathrm{cos}{\displaystyle \frac{q}{2}}\left({\displaystyle \frac{2}{3}}v^2q^2+(1+i\eta vq)(1+2\mathrm{sin}^2{\displaystyle \frac{q}{2}})\right),`$ (40a)
$`D`$ $`=`$ $`(1+i\eta vq)^2(1+3\mathrm{sin}^2{\displaystyle \frac{q}{2}}){\displaystyle \frac{4}{3}}(1+i\eta vq)v^2q^2(1+\mathrm{sin}^2{\displaystyle \frac{q}{2}})+{\displaystyle \frac{1}{3}}v^4q^4.`$ (40b)
Then calculating the residues at these poles we get

$$Q^F(q)=\frac{N(q)}{1+i\eta vq}(1S)$$
(41)
where

$$S(q)=\frac{F(\xi _1)\sqrt{\xi _2^2(1+i\eta vq)^2}F(y_2)\sqrt{\xi _1^2(1+i\eta vq)^2}}{3(\xi _1\xi _2)\sqrt{\xi _1^2(1+i\eta vq)^2}\sqrt{\xi _2^2(1+i\eta vq)^2}}.$$
(42)
Note that we are supposed to take those branches of the square roots in (42) which ensure that

$$|w_{1,2}|<1$$
(43)
with $`w_{1,2}`$ defined by

$$w_{1,2}=\frac{\xi _{1,2}\sqrt{\xi _{1,2}^2(1+i\eta vq)^2}}{1+i\eta vq}$$
(44)
If we now substitute $`N(q)`$ in Eq. (41) into Eq. (18), we get exactly the same equation as was obtained for the mode III case in the first paper (Pechenik, et al., 2000)

$$Q^++SQ^{}+\frac{\eta v}{1+i\eta vq}(1S)Q(0)=\frac{(1S)P_0^{}}{1+i\eta vq}$$
(45)
with $`Q`$ having the same meaning as $``$ in (Pechenik, et al., 2000), namely the elongation of bonds between rows where the crack propagates.

Therefore, we do not need to repeat any of the discussion from our previous paper regarding the choice of $`P_0^{}`$ and the factorization into pieces analytic in the upper and lower half planes. Instead, we can just write down the final answer. For the relationship between the velocity and the driving displacement, we have

$$\frac{\mathrm{\Delta }}{\mathrm{\Delta }_G}=\sqrt{\frac{1+\varphi \eta v}{A\varphi }}\mathrm{exp}\frac{1}{4\pi i}_{\mathrm{}}^+\mathrm{}\frac{\mathrm{ln}K(\chi )}{\chi (1+i\eta v\chi )}𝑑\chi $$
(46)
where for this mode I case

$$A=\frac{1}{\sqrt{2}v^2}\left(\frac{\sqrt{\frac{3}{8}v^2}}{2}\frac{(\frac{3}{4}+v^2)^2}{\sqrt{3}\sqrt{\frac{9}{8}v^2}}\right)$$
(47)
and

$$K(q)=\frac{S^2(q)(q^2+\varphi ^2)}{A^2q^2\varphi ^2}.$$
(48)
This expression for $`A`$ contains the transverse speed $`c_t=\sqrt{3/8}`$ and longitudinal speed $`c_l=\sqrt{9/8}`$. The zero of $`A`$ gives the limiting speed of crack propagation $`c_R=\frac{\sqrt{3\sqrt{3}}}{2}=0.563`$, which is as expected equal to the Raleigh wave speed.

In our previous paper, we discussed the details of how to evaluate this equation so as to obtain numerical results for the crack velocity response curve. These methods can be applied to this case as well and hence we do not repeat any of the details here. The results of these calculations are presented in Figs. 1 and 2. In fact, these findings are quite similar to the analogous mode III results. Later, we will show that the upper limiting value of $`\mathrm{\Delta }/\mathrm{\Delta }_G`$ for an arrested (i.e. non-moving) crack is around $`1.94`$. As was the case for mode III, the $`\mathrm{\Delta }v`$ curves in Figures 1, 2 approach this value as $`v`$ goes to 0 at any value of the damping $`\eta `$. Slepyan’s original calculations, done for the dissipationless limit $`\eta =0`$, show the same asymptotic value.

## IV consistency of the steady-state solution

In this section we calculate the bond displacements in the vicinity of the crack. This is crucial, as we need to check whether all bonds assumed to be unbroken in the derivation in fact have elongations less than unity, the value at which bonds break.

Let us start with bonds along the crack path, where the elongation is just determined by transforming $`Q^F(q)`$ back to physical space. This part follows exactly the analogous mode III calculation, and we easily obtain for $`Q^\pm `$

$`Q^+`$ $`=`$ $`{\displaystyle \frac{Q(0)}{(iq+0)(1+iq\eta v)}}{\displaystyle \frac{S^+}{S^+|_{q=\frac{i}{\eta v}}}}{\displaystyle \frac{\eta vQ(0)}{1+iq\eta v}},`$ (49a)
$`Q^{}`$ $`=`$ $`{\displaystyle \frac{Q(0)}{S^+|_{q=\frac{i}{\eta v}}S^{}}}{\displaystyle \frac{1}{(iq+0)(1+iq\eta v)}}+{\displaystyle \frac{\eta vQ(0)}{1+iq\eta v}},`$ (49b)
and for $`Q(\tau )`$

$`Q(\tau )`$ $`=`$ $`Q(0){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{\sqrt{A\varphi }e^{iq\tau }}{\sqrt{(iq0)(iq\varphi )}(1+iq\eta v)}}{\displaystyle \frac{(K^+)^{\frac{1}{2}}}{S^+|_{q=\frac{i}{\eta v}}}}{\displaystyle \frac{dq}{2\pi }},\text{ for }\tau >0`$ (50a)
$`Q(\tau )`$ $`=`$ $`Q(0){\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{e^{iq\tau }}{\sqrt{A\varphi }(iq+0)(1+iq\eta v)}}\left({\displaystyle \frac{iq+\varphi }{iq+0}}\right)^{\frac{1}{2}}{\displaystyle \frac{(K^{})^{\frac{1}{2}}}{S^+|_{q=\frac{i}{\eta v}}}}{\displaystyle \frac{dq}{2\pi }}`$ (50b)
$`+Q(0)e^{\frac{\tau }{\eta v}},\text{ for }\tau <0.`$
The numerical evaluation of these expressions follows the same methodology as described for mode III. Typical results are shown in Fig. 3. We note that the elongation along the crack line is rather similar to the same object in the mode III case. In that case, however, this function also determined the horizontal bond elongations by simple subtraction. This is no longer true for mode I, since the vectorial nature of the problem requires that we take a different component of the displacement (different than the ones which goes into $`Q`$) to evaluate this elongation. In detail, we now have for

$`Q_h(q)Q_1^F(q,n=0)={\displaystyle \frac{\sqrt{3}}{4\pi }}{\displaystyle _{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}}(e^{ir\widehat{d}_1}1)\stackrel{}{u}^{FF}\widehat{d}_1𝑑s=`$ (51)
$`(e^{iq}1){\displaystyle \frac{\sqrt{3}}{4\pi }}{\displaystyle _{\frac{2\pi }{\sqrt{3}}}^{\frac{2\pi }{\sqrt{3}}}}{\displaystyle \frac{det𝐌_1}{det𝐌}}𝑑s`$ (52)
where these matrices were defined in the last section.

We now proceed as before to change variables to $`w`$ and re-write this expression in terms of the auxiliary variable $`\xi `$ defined in Eq. (34a).

$$det𝐌_1=\frac{N(q)H(\xi )}{1+i\eta vq}+\text{term odd in }s$$
(53)
where

$`H(\xi )`$ $`=`$ $`{\displaystyle \frac{3}{2}}(1e^{i\frac{q}{2}})(1+i\eta vq+\xi )(1+i\eta vq\xi \mathrm{cos}{\displaystyle \frac{q}{2}}{\displaystyle \frac{v^2q^2}{3}})`$ (54)
$`{\displaystyle \frac{3i}{2}}\mathrm{sin}{\displaystyle \frac{q}{2}}(1+e^{i\frac{q}{2}})(\xi ^2(1+i\eta vq)^2)`$
As before, the integrand in (52) has three poles. Dropping the irrelevant odd term and performing the integral leads to the expression

$$Q_h(q)=\frac{N(q)}{1+i\eta vq}(1e^{iq})\psi _h(q)$$
(55)
with

$$\psi _h(q)=\frac{1e^{i\frac{q}{2}}}{2}\frac{H(\xi _1)\sqrt{\xi _2^2(1+i\eta vq)^2}H(\xi _2)\sqrt{\xi _1^2(1+i\eta vq)^2}}{3(\xi _1\xi _2)\sqrt{\xi _1^2(1+i\eta vq)^2}\sqrt{\xi _2^2(1+i\eta vq)^2}}.$$
(56)
Branches of the square root satisfy the same conditions (43), (44), as earlier. The function $`N(q)`$ can be found from Eq. (41)

$$N(q)=\frac{(1+i\eta vq)Q^F(q)}{1S}.$$
(57)
Our final answer is

$$Q_h(q)=\frac{Q(0)}{S^+|_{q=\frac{i}{\eta v}}S^{}}\frac{(1e^{iq})\psi _h(q)}{(iq+0)(1+iq\eta v)}.$$
(58)
For $`q`$ close to $`0`$, the function $`(1e^{iq})\psi _h(q)`$ behaves as $`q`$, giving a divergence $`1/q^{1/2}`$ for $`Q_h(q)`$. This is similar to the divergence in $`Q^+(q)`$. Thus, the numerical calculation of $`Q_h(\tau )`$ is similar to the calculation of $`Q(\tau )`$ for $`\tau >0`$. Figure 4 displays $`Q_h(\tau )`$ for several values of $`\eta `$ and $`v`$. Generally speaking, the function $`Q_h(\tau )`$ has maxima in two different places, one somewhere in the vicinity of $`1`$ and a second for $`\tau >0`$.

We also need to find the bond elongation between the layers with $`n=1`$ and $`n=2`$ and between the layers with $`n=1`$ and $`n=0`$. Due to the symmetry (31) it is sufficient to consider only $`Q_2^F(q,n=\pm 1)`$, which we will denote as $`Q_\pm (q)`$. We can derive for it an expression similar to that in Eq. (38) with the major difference being that in this case the odd parts of the explicit determinant in Eq. (32) do not cancel out because of the additional factor $`w^{\pm 1}`$.

$`Q_\pm (q)`$ $`=`$ $`{\displaystyle \frac{N(q)}{2\pi i}}{\displaystyle \underset{|w|=1}{}}w^{\pm 1}{\displaystyle \frac{\frac{F(\xi )}{1+i\eta vq}+G(w)}{3\xi ^2+\alpha \xi +\beta }}{\displaystyle \frac{dw}{w}}`$ (59)
$`=`$ $`{\displaystyle \frac{N(q)}{2\pi i}}{\displaystyle \underset{|w|=1}{}}{\displaystyle \frac{\frac{F(\xi )}{1+i\eta vq}\pm G(w)}{3\xi ^2+\alpha \xi +\beta }}𝑑w.`$
To derive the second equality, we used a change of variables $`w1/w`$, along with the symmetries of $`F`$ and $`G`$. Again the integration in (59) is in the counter-clockwise direction. In this case, the integrand in the second form of the integral has just two poles inside the unit sphere, at $`w_{1,2}`$. After the by now familiar tedious calculations, we find

$$Q_\pm (q)=\frac{N(q)}{1+i\eta vq}\psi _\pm (q)$$
(60)
where

$`\psi _\pm (q)={\displaystyle \frac{F(\xi _1)w_1}{3(\xi _1\xi _2)\sqrt{\xi _1^2(1+i\eta vq)^2}}}+{\displaystyle \frac{F(\xi _2)w_2}{3(\xi _1\xi _2)\sqrt{\xi _2^2(1+i\eta vq)^2}}}`$
$`\pm {\displaystyle \frac{g(q)(w_1w_2)}{3(\xi _1\xi _2)}}`$ (61)
with

$$g(q)=2i\mathrm{sin}\frac{q}{2}(v^2q^23(1+iq\eta v)\mathrm{sin}^2\frac{q}{2})$$
(62)
and $`N(q)`$ is again given by (57). Finally we obtain

$$Q_\pm (q)=\frac{Q(0)}{S^+|_{q=\frac{i}{\eta v}}S^{}}\frac{\psi _\pm (q)}{(iq+0)(1+iq\eta v)}$$
(63)
Just as was the case for $`Q_h(q)`$ and $`Q^+(q)`$, $`Q_\pm (q)`$ behaves as $`1/q^{1/2}`$ near $`0`$. Thus numerical calculations can be performed just as in the previous cases. Figure 5 shows several curves $`Q_\pm (\tau )`$ for differing parameters.

Given the elongations of these “vulnerable” bonds, we can investigate the critical speed at which one of these bonds should be broken. Figure 6 shows the results of our calculations of this critical speed. In Figure 7 we plotted $`\tau _{cr}`$ at which these bonds break; this will allow us to identify spatial and time coordinates of braking and make contact with our finite lattice calculations later. We found that the maximum value of $`Q_{}(\tau )`$ always reaches 1 before $`Q_+(\tau )`$ and that this is the most dangerous bond for small dissipation. This curve turns around at an $`\eta `$ of around $`1.2`$, so that $`Q_{}(\tau )`$ is always subcritical for larger $`\eta `$. This is a result of the fact that the maximum value of $`Q_{}(\tau )`$ is, surprisingly enough, not monotonic with velocity, but reaches a maximum and then decreases with increasing velocity. The maximum extension of this type of bond occurs, for large velocity, at some distance from the crack surface. For small $`\eta `$, where the maximum $`Q_{}(\tau )`$ does exceed critical extension, the decrease of the maximum $`Q_{}`$ with $`v`$ restabilizes the bond beyond some velocity. The dominant threshold for $`\eta `$ above about $`1.1`$, then, comes from the horizontal bond breaking. Note that as $`\eta `$ increases, the critical velocity decreases. Also, for small $`\eta `$, a crossover occurs as the relative importance of the two different maxima in the horizontal bond elongation reverse; this region is anyway irrelevant as the next-layer vertical bond breaks at a lower velocity. Whenever the horizontal bond dominates, the maximum at $`\tau 1`$ is the governing one. As we see from Figure 7 even for large $`\eta `$ this critical $`\tau `$ stays near $`1`$, what means that horizontal bond always breaks near the tip of the crack, contrary to the Mode III case, where the point of breaking drifts backward with increasing $`\eta `$.

Many of these features are special to the Mode I problem, and have no analog in the Mode III calculation. For Mode III, the horizontal bond breaking always dominates. Furthermore, the maximum off-crack surface bond extension is a strictly increasing function of the velocity. This points to the possibility that the dynamics of Mode I may be much richer than the Mode III dynamics.

## V Finite lattice model.

For our mode III calculations, it proved interesting to compare the exact results derived for an infinite lattice with the numerical determination of crack propagation properties in lattices of small transverse size (Kessler and Levine, 1998; Kessler, 2000). We now discuss similar calculations for our mode I model. In addition to providing details regarding the size needed to attain answers relevant to the macroscopic limit (a question of direct relevance for direct molecular dynamics simulations (Abraham, et al., 1994; Zhou, et al, 1997; Gumbsch, et al., 1997; Holland and Marder, 1997, 1998; Omeltchenko, et al., 1997), for example), finite lattice results can be used as a rough check on some of the predictions obtained above. Given the complexity of the analysis, having such a check is quite useful.

### V.1 Arrested crack

Let us start with trying to find arrested cracks, expected to exist in some range of the driving $`\mathrm{\Delta }`$ around the Griffith’s displacement. This involves looking for a solution of Eq. (1) with $`\stackrel{}{u}(x,y)`$ independent of time and with the forces on the right hand side arising solely from the broken bonds. If we have a system with a finite number of rows, the natural boundary condition requires that the top and bottom rows have fixed vertical displacements of $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }`$ respectively. Since the entire system is linear, we can choose to use $`\mathrm{\Delta }=1`$ and rescale the breaking criterion accordingly.

Since we are doing a numerical calculation, we must introduce artificial boundaries in the direction along the crack, $`\widehat{x}`$. What we do is cut off the range of points whose coordinates are variables and outside of this range impose fixed asymptotic displacements. On the cracked side, these are just $`\stackrel{}{u}=\pm \mathrm{\Delta }\widehat{y}`$ for positive and negative $`y`$. For the uncracked side, the asymptotic displacement involves constant strain. At each of the sites that has a variable displacement, we impose the two components of the (vectorial) equation of motion. We also impose the equations of motion at the boundaries of the system. This gives us more equations than we have variables. The coupling of sites which contain variable displacements with the ones with fixed displacements gives rise to inhomogeneous terms. Combining these observations, we can write our system in the schematic form $`Mub=0`$, with a non-square matrix $`M`$. The field $`u`$ is then determined by the requirement that the error be minimal. In this manner, the small errors introduced at the boundary by having to have a fixed box size are prevented from causing any large errors (by modes which grow exponentially away from the edges) in the bulk of the lattice.

The solution of this linear system determines the elongation of all bonds. In Figure 8 we show a typical lattice given by this solution. In fact we need to know elongation of just two bonds right at the tip of the crack; the one which in the case of the moving crack would break next ($`Q_n`$) and the other which would have been the last one broken ($`Q_l`$). Now, recall that we have scaled our displacement to equal unity and the only remaining displacement scale is the elongation at which the springs break, which we can call $`ϵ`$. The solution we have found is consistent as long as $`ϵ`$ is between the lower limit set by the $`Q_n`$ and the upper limit set by $`Q_l`$. Since $`\mathrm{\Delta }_G=ϵ\sqrt{(2N+1)/3}`$, where $`2N`$ is the number of rows in the $`y`$-direction (excluding the boundary rows whose displacement is constrained), we directly obtain the upper and lower limits of the arrested crack band

$$\frac{1}{\sqrt{\frac{1}{3}(2N+1)}Q_l}\mathrm{\Delta }/\mathrm{\Delta }_G\frac{1}{\sqrt{\frac{1}{3}(2N+1)}Q_n}$$
(64)
The results for these thresholds as a function of lattice size are shown in Figs. 9, 10. In the limit of infinite lattice, the lower limiting value $`\mathrm{\Delta }_{}/\mathrm{\Delta }_G`$ approaches $`0.515`$, while the upper limiting value $`\mathrm{\Delta }_+/\mathrm{\Delta }_G`$ approaches $`1.94`$.

### V.2 Stable moving cracks

We now turn to the moving crack problem. Again, we need to solve Eq. (1) with the forces in the right part due to the broken bonds. Now, the displacement field $`\stackrel{}{u}(x,y,t)`$ is of course time dependent. Therefore, we need to introduce a time-step $`\mathrm{\Delta }t`$ so as to define the time-points at which we will obtain a numerical solution. Given some specific speed $`v`$ of the crack, we choose to divide half of the time interval that it takes for the crack to propagate one lattice spacing, $`1/(2v)`$, into $`n`$ equal time intervals; thus $`\mathrm{\Delta }t=1/(2vn)`$. Now we can discretize our system, obtaining equations at times $`t=0,\mathrm{\Delta }t,\mathrm{},(n1)\mathrm{\Delta }t`$. We use a symmetric discretization of first and second derivatives $`(\stackrel{}{u}(t+\mathrm{\Delta }t)\stackrel{}{u}(t\mathrm{\Delta }t))/(2\mathrm{\Delta }t)`$ and $`(\stackrel{}{u}(t+\mathrm{\Delta }t)+\stackrel{}{u}(t\mathrm{\Delta }t)2\stackrel{}{u}(t))/\mathrm{\Delta }t^2`$ correspondingly. These equations depend on displacements $`\stackrel{}{u}`$ outside of this time interval, because of the temporal derivatives. These can be found via use of the assumed symmetries of the moving crack; in fact, it is easy to see that we thereby trade in displacements outside of the modeled interval with displacements inside this interval albeit at a different spatial location. The boundaries are treated the same way as described for the arrested crack; the displacements outside some range are replaced by asymptotic values for all of the time-points in our interval. Including the equations at these boundary points again gives us a system in which the number of equations exceeds the number of variables and again a least squared error algorithm is used to find the required solution. Fig. 11 shows a snapshot of the lattice near the tip of the crack at one particular time for one set of parameters.

We used these finite lattice calculations to provide an independent check on our analytic infinite lattice results. We were specifically interested in checking the critical speed estimates. The difficulty is that convergence in the parameter $`1/N`$ is fairly slow; and, using large lattices is rather time-consuming and memory-demanding. In practice, we limited our calculations to $`N`$’s ranging from $`6`$ to $`24`$. It turns out that for small dissipation where the next-row vertical bond vulnerability is the thing which determines the critical velocity, we can do a credible job of verifying that the infinite lattice results are consistent with the finite lattice ones. For example, in Fig. 12, we show the critical speed as a function of lattice size for several different values of $`\eta `$. These numbers, if extrapolated to infinite $`N`$ are clearly consistent with the results given earlier in Fig. 6, especially considering that having a finite time-step does introduce a small numerical error on its own. We can also check the critical time $`\tau `$ at which the bond breaking occurs. From Fig. (7), we see that in the infinite lattice the break occurs for $`\tau 0.5`$. For the finite lattice with $`n=5`$, this mean that the bond $`a`$ highlighted in Fig. 11 should go above breaking threshold at the first time-point of the modeled interval. This is exactly what occurs.

A different strategy is necessary for studying the onset of horizontal bond breaking. This is due to the fact that quite large $`N`$’s are required in order to for the results to quantitatively approach the infinite $`N`$ limit. This is another striking difference between the Mode I results and those of Mode III. For Mode III, the qualitative nature of the extraneous bond breaking was always similar to that of the infinite width limit, and the results were quantitatively accurate even for fairly small $`N`$’s. Here, however, the picture of which bonds are critical changes dramatically between finite $`N`$ and the infinite $`N`$ limit. In Fig. (13), we present the “phase diagram” for critical bonds for $`N=38`$, plotted together with the infinite $`N`$ result, and in Fig. (14), we compare the $`N=38`$ results to those for $`N=26`$. Since the steady-state code is difficult to run at these large values of $`N`$, the data in these two plots were obtained from time-dependent simulations, with only the central vertical bonds allowed to break. After running to times of $`t=100`$ to eliminate transients, the extensions of the dangerous bonds were checked to see if they exceeded criticality. Note that, in contradistinction to the “phase diagram” plotted in Fig. (6), we now plot $`\mathrm{\Delta }/\mathrm{\Delta }_G`$ on the vertical axis, since this is the input control parameter for the simulations. We see that qualitative behavior of the “$`Q_{}`$” bonds is similar to the infinite-$`N`$ result, with these bonds being below critical extension for both large $`\eta `$ and large $`\mathrm{\Delta }/\mathrm{\Delta }_G`$. The horizontal bonds also behave qualitatively like their infinite-$`N`$ counterparts, but the quantitative agreement, is as we noted above, significantly worse. In fact, for large $`\eta `$, the most dangerous bond is in fact of “$`Q_+`$” type. However, the threshold at which the “$`Q_+`$” type bond becomes dominant is pushed to larger $`\eta `$ as the system size is increased, presumably going to infinity with $`N`$. Also, there is a small region of $`\eta `$ around 1.1 where the inconsistency is re-entrant, so that for intermediate values of $`\mathrm{\Delta }`$, no bonds of the crack surface are critical. Thus, to see dynamics qualitatively representative of the macroscopic system requires very large system sizes at moderate to large $`\eta `$. Presumably, this is connected to the process zone increasing in size with $`\eta `$. Again, it is worth emphasizing that this is a feature not present in Mode III, and leads to the conclusion that in molecular dynamics simulations, the width of the material should be taken very large to accurately study the micro-branching instability.

## VI Discussion

In this work, we have discussed steady-state mode I crack propagation in a viscoelastic lattice model. Our primary method of analysis utilizes the Wiener-Hopf technique to write down closed form expressions for both the $`v\mathrm{\Delta }`$ curves (for various values of the dissipation $`\eta `$) and for the bond elongation field. The latter enables us to define a limit of consistency for the solution past which some bonds not along the crack path have elongations greater than the assumed breaking criterion. These results are new, and correspond to non-trivial extensions of the work of Kulamekhtova (1984) on the dissipationless limit and the work of Marder and Gross (1995) on small lattices.

The most interesting results, in our opinion, concern the dependence of the critical velocity (where the aforementioned inconsistency sets in) on the amount of dissipation. For small dissipation, this threshold is relatively insensitive to $`\eta `$, as already suggested in direct numerical simulations. This threshold, occurring at roughly .73 of the Rayleigh speed is in some ways reminiscent of the Yoffe (1951) branching criterion which suggests that straight crack propagation will become unstable once the largest stress direction shifts away from being straight ahead. It is unfortunately hard to be more precise regarding this correspondence since the crack dynamics in our model is fundamentally tied to the lattice scale, not the macroscopic scale - in fact, the latter is completely invisible in the Slepyan approach aside from providing a driving force in the form of a stress intensity factor.

At larger $`\eta `$, the instability picture changes. Now, it is a horizontal bond breaking which signals the onset of more complex crack dynamics. Also, the threshold goes down with increasing dissipation. This was not the case for our mode III calculations. This instability has nothing to do with Yoffe, as it is strongly dissipation dependent and in any case is not associated with crack branching.

Much of the recent theoretical work on mode I cracks has been motivated by experiments which show clearly that instabilities limit the range of steady-state crack propagation. These instabilities introduce more complex spatio-temporal dynamics to the fracture process, causing additional dissipation and leaving behind a roughened crack surface. It has been tempting to associate these results with the onset of inconsistencies in lattice models, although the details of this correspondence remains uncertain. First, most of the experimental work has been carried out in amorphous materials, making the idea of an ordered lattice model somewhat suspect. Next, the instabilities seen experimentally typically occur at smaller speeds than the ones seen in lattice systems with small dissipation. Finally, the experiments seem to show a typical frequency for micro-branching which is not connected in any obvious way to dynamics at a lattice scale.

Notwithstanding all these issues, we remain optimistic that the study of this class of models will lead to insight into dynamic fracture. There are several intriguing possibilities that need to be investigated in the future. First, we have shown that for mode I cracks, increasing the dissipation (in the form of a Kelvin viscosity) eventually results in a decrease of the instability threshold. If our model is really applied at the atomic scale, it is hard to see why there should be a large linear dissipation; on the other hand, it is well-known that lattice models miss an essential non-linear dissipation mechanism, namely the creation of dislocations which remain pinned to the crack. Would inclusion of these effects also lower the threshold? On the other hand, applying the model on a large scale (possibly for a disordered system) would naturally require a large dissipation and recent numerical simulations indicate that proper inclusion of the thermal fluctuations might also push the model into better agreement with experiment (Pla, et al., 1998; Sander and Ghasias, 1999). Finally, the exact nature of the state which occurs past the instability onset has not been addressed in our work to date, and in fact cannot be addressed by the elegant but ultimately limiting analytic methods utilized for the steady-state problem. Instead, we plan to study a generalized force law in which the sharp breaking criterion is replaced by an analytic nonlinear force versus displacement (Kessler and Levine, 1999). In this formulation, the inconsistency found here becomes a linear instability of the steady-state crack (Kessler and Levine, 2000) and one can use some of the methods developed in the field of nonequilibrium pattern formation to unravel the dynamics past onset. These studies, together with additional experimental data regarding the differences between brittle fracture in crystalline versus non-crystalline materials, will hopefully lead to a better understanding of dynamic fracture.

###### Acknowledgements.

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 initial phase of this work.