% ****** Start of file template.aps ****** %
%
%   This file is part of the APS files in the REVTeX 3.0 distribution.
%   Version 3.0 of REVTeX, November 10, 1992.
%
%   Copyright (c) 1992 The American Physical Society.
%
%   See the REVTeX 3.0 README file for restrictions and more information.
%
%
% This is a template for producing files for use with REVTEX 3.0.
% Copy this file to another name and then work on that file.
% That way, you always have this original template file to use.
%
\documentstyle[preprint,aps]{revtex}
\begin{document}
%\draft command makes pacs numbers print
\draft
\title{Crack Stability and Branching at Interfaces.}
% repeat the \author\address pair as needed
\author{Robb Thomson, NIST Fellow, Emeritus}
\address{Materials Science and Engineering Laboratory\
National Institute of Standards and Technology\
Gaithersburg, MD 20899}
\date{\today}
\maketitle
\begin{abstract}
% insert abstract here
The various events which occur at a crack on an interface are
explored, and described in terms of a simple graphical construction
called the crack stability diagram.  For simple Griffith cleavage in a
homogeneous material, the stability diagram is a sector of a circle in
the space of stress intensity factors, $K_I/K_{II}$. The Griffith 
circle is  limited in both
positive and negative $K_{II}$ directions by nonblunting dislocation 
emission
on the cleavage plane.  For a branching plane inclined at an angle to the
original cleavage plane, both cleavage and emission (which blunts the
crack) can be described as a balance between an elastic driving force
and a lattice resistance for the event.  We use an analytic expression
obtained by Cotterell and Rice for cleavage, and show that it is an
excellent approximation, but show that the lattice resistance includes
a cornering resistance, in addition to the standard surface energy in
the final cleavage criterion.  Our discussion of the lattice
resistance is derived from simulations in 2D hexagonal lattices with
UBER force laws with a variety of shapes. Both branching cleavage and blunting
emission can be described in terms of a stability diagram in the space
of the remote stress intensity factors, and the competition between
events on the initial cleavage plane and those on the branching plane
can be described by overlays of the two appropriate stability
diagrams.   The popular criterion that $k_{II}=0$ on the branching
plane is explored for lattices and found to fail significantly,
because the lattice stabilizes cleavage by the anisotropy of the surface
energy. Also, in the lattice, dislocation emission must always be
considered as an alternative competing event to branching.
\end{abstract}
% insert suggested PACS numbers in braces on next line
\pacs{PACS numbers: 07.05.T, 62.20.M, 61.72.B, 61.72}

% body of paper here

\section {Introduction}
A crack on an interface in a material can take any of four 
different event paths:
It can cleave straight ahead on the interface plane, emit a
dislocation on this plane, or it can branch in cleavage on an inclined
cleavage plane, or it can emit a dislocation on an inclined slip
plane.  Depending on the symmetry of the lattice and the geometry of
the crack line and cleavage plane, a number of slip or branching
planes may be available.  

The subject of this paper is to give a
unified presentation of the criteria for these various events, and
introduce a graphical representation from which it is possible to pick
off what path a crack under a particular (mixed) load will take, given
that the force laws and crystal type are known.  Since it is
impossible to give criteria for the general crystal including all
possible bonding types, we will illustrate the process with a generic
2D hexagonal lattice with a coherent interface, using UBER and
associated pair forces.  The criterion for branching is presumably
associated with a Griffith condition with an appropriate
crack driving force (energy release rate) for the particular plane.
The emission condition is more complicated, and will be analyzed in
terms of Rice's unstable stacking fault energy\cite{R}, $\gamma_{us}$
and its extensions\cite{ZCT,TC,T}.

Both emission and cleavage events must be a balance between a driving
force for the event, and a lattice resistance to the event.  Since
branching is a cleavage process, the driving force must be the
standard energy release rate on the branching plane, written in terms
of the local stress intensity factors, $k_i$, for the branching crack
in the limit of zero length for the kink.
(We will alternatively refer to the branching event as a kinking
of the original crack, and the emerging crack on the branching plane
as a kink.)  The limiting local stress intensity factors for the kink
are not obtainable from the stress intensity factors of the main crack
in an analytic form, but must be obtained
numerically\cite{BC,BCH,DR,L}.  The more complicated interfacial case
has been analyzed by He and Hutchinson\cite{HH}.  Even though analytic
results are not possible, Cotterell and Rice\cite{CR} have shown that
the kink stress intensity factors are closely related to the stress
fields surrounding the main crack, and have written an analytic
approximation for the local $k$'s which is remarkably accurate.  Since
the continuum approximation introduces inaccuacies of its own, the
Cotterell/Rice expressions may be as close as one can come to these
quantities in the continuum approximation, and will be used as
the analytic basis for our lattice comparisons.  




The correct driving force for cleavage has not been clearly
identified until lattice studies\cite{AR} showed it to be
given by the standard mechanics expression
\begin{equation}
{\cal G}_c={1\over 2\mu'}(k_I^2+k_{II}^2),
\label{eqn1}
\end{equation}
where $\mu'$ is an appropriate elastic shear modulus. It has
sometimes been claimed\cite{T1} that the cleavage condition should be
predominantly  a
matter of the Mode I loading, because the opening mode is the only one
to actually force the cleavage planes apart; otherwise the crack should
simply reclose itself.
But the lower limit for $k_I$ is actually determined
by the fact that the crack will break down in shear and emit a
dislocation when the shear load is sufficiently great on the Griffith
circle,  and for shear loads below this
limit, the criterion is given by the $\cal G$ formula above.  Thus,
contrary to He and Hutchinson\cite{HH}, we will always take the local $\cal
G$ on the kinking crack as the appropriate driving force.  When this
quadratic expression is transformed back to the remote loading stress
intensity factors, of course, a more complex quadratic form in $k_I$
and $k_{II}$ results.


The appropriate driving force for emission has been proposed by Rice\cite{R}
to be the Mode II component of $\cal G$, or 
\begin{equation}
{\cal G}_e={1\over 2\mu'}k_{II}^2.
\label{eqn2}
\end{equation}
Again, when this equation is transformed back to the remote or lab system of
stress intensity factors, a more complex form results, to be discussed
in \S III.  

The lattice resistance for an event will be computed in the
simulations to be reported.  One would expect that the lattice
resistance for cleavage should simply be the intrinsic surface energy
of the lattice, but we will show that a special resistance due to
turning the corner of the kink apparently comes into play, so that the
standard Griffith relation is not satisfied exactly in this case. 

For emission, the lattice resistance has been the subject of intense
study, and it has been shown\cite{ZCT2} 
that when the crack is not blunted by the
emission, the lattice resistance is simply given by the unstable
stacking fault, $\gamma_{us}$, as proposed by Rice\cite{R}.  When the
crack is blunted by the emission, however, the lattice resistance is
more complex, and involves a product of $\gamma_s\gamma_{us}$, as shown
by Zhou, Carlsson and Thomson\cite{ZCT} and by Thomson and
Carlsson\cite{TC}. Blunting emission for the interface was studied by
Thomson\cite{T}. 

The lattice modeling will be done with the methodology described
elsewhere\cite{T,TZCT} for the 2D hexagonal lattice with nearest
neighbor central force laws, and the reader is referred to the
references for the details of the simulations.  The reason for
choosing the hexagonal lattice is that it is elastically isotropic in
the continuum limit, so anisotropic corrections need not be considered
for the analysis to be carried out.  Also, our purpose is to
illucidate the generic physics of cracks in lattices, and a method
which gives quick results for a variety of different force laws in a
simulation where the load mix and lattice mismatch can be varied
easily is desired.  

\section{Analysis.}

The Cotterell/Rice\cite{CR} prescription for the local stress
intensity factor of a kinking crack when applied to an interface is
given by\cite{RSW}
\begin{eqnarray}
k_I&&=\sigma_{\theta\theta}\sqrt{2\pi r}={\cal K}_I\Sigma^I_{\theta\theta}
+{\cal K}_{II}\Sigma^{II}_{\theta\theta}\nonumber\\
k_{II}&&=\sigma_{r\theta}\sqrt{2\pi r}={\cal K}_I\Sigma^I_{r\theta}
+{\cal K}_{II}\Sigma^{II}_{r\theta}\nonumber\\
{\cal K}&&={\cal K}_I+i{\cal K}_{II}=|K|e^{i(\psi-\eta)}\nonumber\\
e^{i\eta}&&=\left({2a\over r}\right)^{i\epsilon}\nonumber\\
\epsilon&&={1\over
2\pi}\ln{\kappa_1\mu_2+\mu_1\over\kappa_2\mu_1+\mu_2}
={1\over 2\pi}\ln{11(c_1/c_2)+5\over 11+5(c_1/c_2)}.
\label{eqn3}
\end{eqnarray}
$K$ is the (complex) remote load stress intensity factor of the bulk 
material (i.e. without the interface) written as
$K=K_I+iK_{II}=|K|\exp{i\psi}$.  $\psi$ is the phase angle of the 
remote load. The stress intensity factor of the unkinked interfacial crack is
written as $\cal K$ with the connection given above to the local
stress intensity, $k_{II}$, and to the shear stress,
$\sigma_{r\theta}$. This definition for the stress intensity factor
differs slightly from that in common use\cite{RSW}, but is appropriate for the
crack and load geometry in use in this work. 
In these equations, the additional phase angle at
the crack tip generated by the elastic mismatch at the interface is
given by $\eta$.  $\eta$ is a singular logarithmic function
of the distance, $r$, from the crack tip---the mode mixing anomaly
characteristic of interfacial cracks in the continuum limit.
$\epsilon$ is a constant which depends on the
elastic mismatch, where $\kappa$ and $\mu$ are the standard isotropic
elastic parameters for the two materials.  The second form for
$\epsilon$ in the last equation in (\ref{eqn3}) refers to the 2D
hexagonal lattice, with the two spring constants, $c_1$ and $c_2$.

The approximation by Cotterell and Rice\cite{CR} is to 
recognize that the shear
stress when normalized by the square root of the radius has the
dimensions of the stress intensity factor, and when written in terms
of the angular variables, $r$ and $\theta$, has the actual character
of the appropriate tensile or shear stress intensity factor for the
branching crack.  This approximation is rigorous, of course, in the
limit of small $\theta$.

Explicit expressions for the $\Sigma$'s are given by\cite{RSW}
\begin{eqnarray}
\Sigma_{\theta\theta}^I&&={1\over\cosh\pi\epsilon}
\left[\sinh((\pi-\theta)\epsilon)\cos{3\theta\over 2}+
e^{-(\pi-\theta)\epsilon}\cos{\theta\over2}(\cos^2{\theta\over 2}
-\epsilon\sin\theta)\right]\nonumber\\
\Sigma_{\theta\theta}^{II}&&=-{1\over\cosh\pi\epsilon}
\left[\cosh((\pi-\theta)\epsilon)\sin{3\theta\over 2}+
e^{-(\pi-\theta)\epsilon}\sin{\theta\over 2}
(\sin^2{\theta\over 2}+\epsilon\sin\theta)\right]\nonumber\\
\Sigma^I_{r\theta}&&={1\over\cosh\pi\epsilon}\left[\sinh((\pi-\theta)\epsilon)
\sin{3\theta\over 2}+
e^{-(\pi-\theta)\epsilon}\sin{\theta\over 2}(\cos^2{\theta\over 2}
-\epsilon\sin\theta)\right]\nonumber\\
\Sigma^{II}_{r\theta}&&={1\over\cosh\pi\epsilon}
\left[\cosh((\pi-\theta)\epsilon)
\cos{3\theta\over 2}+
e^{-(\pi-\theta)\epsilon}\cos{\theta\over 2}(\sin^2{\theta\over 2}
+\epsilon\sin\theta)\right].
\label{eqn4}
\end{eqnarray}

It is more useful to write these complicated expressions in terms
of a set of angular phase shifts in the form
\begin{eqnarray}
k_I'&&=
\Sigma_{\theta\theta}K\,\cos(\psi-\eta_e+\lambda_{\theta\theta})\nonumber\\
(\Sigma_{\theta\theta})^2&&=(\Sigma_{\theta\theta}^I)^2
+(\Sigma_{\theta\theta}^{II})^2\nonumber\\
\lambda_{\theta\theta}&&=-\tan^{-1}(\Sigma_{\theta\theta}^{II}/
\Sigma_{\theta\theta}^I)\nonumber\\
k_{II}'&&=
\Sigma_{r\theta}K\sin(\psi-\eta_e+\lambda_{r\theta})\nonumber\\
(\Sigma_{r\theta})^2&&=(\Sigma^{II}_{r\theta})^2+(\Sigma^I_{r\theta})^2
\nonumber\\
\lambda_{r\theta}&&=\tan^{-1}(\Sigma^I_{r\theta}/\Sigma^{II}_{r\theta})
\nonumber\\
e^{i\eta_e}&&=\left({2a\over r_0}\right)^{i\epsilon}.
\label{eqn5}
\end{eqnarray}
These last equations are written for a critical value of the local stress
intensity factor (with a prime) corresponding to emission, cleavage,
etc. The total local phase shift for the crack is composed of the
separate contributions from the phase shift for the remote load, $\psi$; the
geometrical contribution due to the kink angle, $\lambda$; and the
interfacial phase shift, $\eta$, written in terms of the core size,
$r_0$.   $\lambda$ is a function of the interface mismatch, as well 
as the kink angle, $\theta$. 
Earlier\cite{T}, we have found that the core size, $r_0$,
can simply be set equal to the range parameter in the force
law.  The reader should note that the stress field of the
kinked crack is not phase shifted by the length of the kink, because the
crack is now off the interface in homogeneous material.  However, the
stress intensity factor in (\ref{eqn3}) does contain the phase 
shift\cite{HH}. The physical reason why the
local stress intensity of the kinked crack is phase shifted by the
interfacial mismatch is that the kinked crack is ``loaded'' by
the phase shifted stress field of the parent interface crack.  

The criteria for critical events in terms of the local stress
intensity factors are written with the subscripts corresponding to 
cleavage or emission.  Thus
\begin{equation}
{\cal G}_c={(k'_{Ic})^2+(k'_{IIc})^2\over 2\mu'_{eff}}
\label{eqn6}
\end{equation}
for cleavage, and 
\begin{equation}
{\cal G}_{IIe}={(k'_{IIe})^2\over 2\mu'_{eff}},
\label{eqn7}
\end{equation}
for emission.
In these equations, an appropriate shear modulus must be taken.  If
the crack is kinked off the interface into material 2, then
$\mu'_{eff}=\mu'_2$, while if the crack is on the interface, then
$\mu'_{eff}=2\mu'_1\mu'_2/(\mu'_1+\mu'_2)$, and the prime on the
$\mu$ represents taking the appropriate plane strain or plane stress
elastic modulus.  In plane stress, which must be used in the 2D
simulations, $\mu'=\mu(1+\nu)$, where $\mu$ and $\nu$ have their
standard elasticity definitions.  The appropriate ``right hand sides''
for the lattice resistance in these equations are to be determined 
by  computer simulations.

The accuracy of these Cotterell/Rice expressions are easily checked in
the zero mismatch case from their paper, where comparisons are graphed
as a function of the branching angle\cite{CR}.  The main deviation
appears in the value for $\Sigma_{r\theta}^{II}$, where the analytic
expression is about 10\% higher than the correct numerical value
at $\theta=60^\circ$.  All other $\Sigma$'s are accurate to
within a percent or so. 


\section{Crack Stability Diagram.}

The crack stability diagram is a generalization of a graphical
description of the Griffith criterion for simple homogeneous cracks.
That is, for a straight crack in homogeneous
material, the Griffith condition is 
\begin{equation}
{\cal G}_c={k_{Ic}^2+k_{IIc}^2\over 2\mu'}=2\gamma_s,
\label{eqn8}
\end{equation}
and is graphed in $k_I/k_{II}$ space as a simple circle.  In continuum
elasticity, the crack is stable so long as the loads correspond to
stress intensity factors lying on the Griffith circle.  In the
lattice, a crack will be stable over a region in $k_I/k_{II}$--space
containing the Griffith circle, because of lattice trapping.  More
important, the stability of the crack will be limited by the ability
of the core of the crack to sustain shear stresses, and beyond some
critical shear, the lattice will break down with emission of a
dislocation. A diagram showing such a stable
region of loading is shown in Fig. \ref{f1} (Use axes $k_I/k_{II}$).  
As explained in the
introduction, when the crack is on an interface, the critical
parameter for cleavage or emission is not the bare remote or 
lab stress intensity
factor, but the local quantity for the core, Eqn. (\ref{eqn5}).  

A generalization of the stability diagram valid for the interface, 
for simple straight cleavage on the initial cleavage plane, or for
nonblunting emission on that plane, will now look like the full
Fig. \ref{f1}. Here, the
local stress intensity factors are rotated relative to the remote axes by
the core phase angle, $\eta_e$, where the core phase angles for
emission and cleavage are assumed to be the same.  It is important to note 
that, except for a $\cosh\pi\epsilon$ term, 
the $\cal G$--force for a straight interface crack is
an invariant in the crack length, and does not contain the phase 
shift\cite{MS}.

When the crack branches or kinks off the interface ($\theta=60^\circ$
is assumed), the more general
expressions for the local stress intensity factors must be used, and
now a contour plot for the values of constant $\cal G$ plotted in the
remote or lab system of stress intensity factors, is shown in Fig. \ref{f2}.
In the figure, the simple case of zero lattice mismatch (no interface)
is shown contrasted with that where the spring constants of the two
lattices differ by a factor of 10 ($c_1/c_2=10$).  Note that
even when the interface mismatch disappears, the maximum gradient of
$\cal G$ is rotated off the axes, and the presence of a remote Mode II
enhances the tendancy to branch.  

The crack stability diagram is completed when the limiting values of
cleavage corresponding to dislocation emission and shear breakdown are
provided by means of emission criteria for the particular force law.
These criteria and actual points for specific cases  are plotted on 
these diagrams in the next section.  

Finally, events on the inclined (kinking) plane will compete 
with events on the
initial interface, and by plotting event loci on the stability
diagram, it is possible to determine when, for example, emission of a
dislocation on the inclined plane will occur before cleavage on the
interface, etc.

\section{The Lattice Resistance.}

The lattice modeling uses the same lattice Green's function techniques
used in previous work\cite{T}, and the general methodology is given in
Thomson, etal.\cite{TZCT}  Also we use the same 2D hexagonal lattice
with the same set of nearest neighbor forces.
Figure  \ref{f3} shows the lattice with a crack on an interface
between atoms of one kind below and a second kind above.  In the
lattice case, we are at liberty to make the bonding between the layers
different from that of either bulk, as in the real physical situation.
The elastic constant of material 2 is arbitrarily normalized to unity,
since all the physical results scale with the elastic constant.  
Also, we allow the crack to emit a
dislocation into material 2, only, and assume that material 1 is 
strong, incapable of deforming or cleavage.  We also make the 
assumption that dislocation emission takes place on the ``forward'' 
slip plane, $\theta=60^\circ$,
because the shear stress is largest on this plane.  Subsequent
emission could conceivably take place on the second slip plane,
$\theta=120^\circ$,  but multiple
emission of dislocations is not explored in this work.  The crack is
loaded at the center of the crack with a concentrated load, see
Fig. \ref{f3}.   This
method of loading means that the crack system is stable, because the
$K$ at the crack tip decreases as the crack length increases by the
elastic equation
\begin{equation}  
K={F_0\over\sqrt{\pi a}},
\label{eqn9}
\end{equation}
where $F_0$ is the point load on upper and lower crack planes, and $a$
is the half length of the parent crack.

As in previous work, our simulations are for a bimaterial slab.  
The slab is $2\times 10^3$ atoms thick,
with the interface running down the center.  The slab has periodic
boundary conditions in the lateral direction, again with repeat
distance of $2\times 10^3$ atom spacings.  The crack itself is 201
atom spacings in total length, the cohesive zone is 12 atom
spacings long on the cleavage plane to the right, and the inclined slip plane
is 16 atom spacings long.  Thus we have no worries about short crack
effects, or interactions with neighboring cracks in the repeating
cells, or with the free surfaces.

In contrast with our previous work, we have restricted the force law
in this work to be the UBER of Rose,
etal.\cite{RSF} derived from the energy expression,
\begin{eqnarray}
E(u)&&=-{c\over 1-e^{-u_0/\alpha}}\left[\alpha^2(1+{u\over 
\alpha})e^{-u/\alpha}
+ {u^2-e_0\over2}e^{-u_0/\alpha}\right]\nonumber\\
e_0&&=u_0^2+2u_0\alpha+2\alpha^2.
\label{eqn10}
\end{eqnarray}
$c$ is the lattice spring constant, $u$ is the displacement from the
equilibrium distance between two atoms, and $\alpha$ is the range
parameter. The lattice parameter will be normalized to unity. 
This equation has been modified from the standard UBER
expression so it has zero energy and force at the second neighbor
position, $r=u_0+1=\sqrt{3}$. It will be convenient to normalize the
lattice spring constant in sublattice 2 to unity, $c_2=1$.


We have performed a set of simulations for a variety of elastic
mismatches, Figs. \ref{f4},\ref{f5}.  In all these simulations,  
the range parmeters are 
held fixed ($\alpha_{12}=0.2,\ \alpha_2=0.13$),  
and $c_{12}=c_2=1.0$.  Thus, for all 
cases, the interface energy and unstable stacking fault in material 2
are held constant, $\gamma_{us}=0.0116$ and $\gamma_s=0.0284$, in 
the natural units of the simulation (lattice parameter and spring
constant, $c_2$, normalized to unity). Note that with these choices,
the bonds in the interface are
stronger than in the matrix.  This
choice amounts to an assumption about the kind of chemical interaction
taking place between the two interfaces to the right of the crack tip, and
the dislocation emission criterion depends on this chemistry.
Nevertheless, the choice we make allows us to explore the form of the
stability diagram, which is our purpose.  We will return at a later
point to discuss the implications of the interfacial chemistry
problem. 


In these simulations, we have
arbitrarily cut the bonds to the left of the base of the spur, so that
the crack does not have to grow through ``good material'' to reach the
branching plane.  This corresponds to a chemical knife which has acted
over the cleavage plane up to the crack tip, in the time honored
manner of the continuum mechanics community.  That is, a chemical
adsorbate is assumed to have formed over the cleavage plane which has
zero interaction with the atoms on the opposite cleavage plane--highly
unlikely physically, but perfectly possible, mathematically.  It
has the advantage of allowing the branching crack to grow out of the
tip, without regard for the problem of how it got there.
Specifically, the first active bond at the tip is the one at
120$^\circ$ from the base of the right hand spur line where it meets
the lower cleavage plane, and is depicted as the first connected bond 
shown in Fig. \ref{f3}.

In Figs. \ref{f4},\ref{f5}, as the solid line, we plot the 
$\cal G$ appropriate to a
simple Griffith law, where $\gamma=\gamma_s$ for material 2.  The
individual points are the simulation results, and they span the 
complete range of possible mixed loads, $K$, for which
branching cleavage is possible.  The upper limits correspond to
lattice shear breakdown and dislocation emission on the inclined
plane, while the lower values correspond to the lowest loads where the
crack is stable on the inclined plane.  
 

The results show that the expected Griffith condition is not satisfied
for the branching crack, but that an extra ``cornering'' energy is
required to turn the crack out of its initial cleavage plane.  This
cornering energy is increased by increasing the lattice mismatch.  
The line of points, however, is parallel to the expected $\cal
G$--curve, showing that the driving force for cleavage is correctly
given by the sum of the squares of the local $k$'s in (\ref{eqn6}).
The emission point obtained is consistent with the analytic emission
criterion given in the previous paper\cite{T}.
We note again, that the lower limit in the
line of simulation points simply corresponds to the load where the
first bond of the branching crack is not broken.  For lower loads, the
crack opens up to the artificially created crack tip, but cannot go
further. 

In the case of the zero misfit, the numerical results for ${\cal
G}_c$ are known, and can be compared with the actual stability diagram
in Fig. \ref{f4}.  As noted in \S II, the analytic expression is
larger than the numerical value, so the cornering energy is slightly
larger than the 30\% shown in Fig. \ref{f4}.  Unfortunately, the
numerical results of He and Hutchinson\cite{HH} are given for 
misfit parameters which differ completely
from those of the present work, so we cannot comment on what
portion of the cornering energy in Fig. \ref{f5} might be due to an
error in the analytic Cotterell/Rice approximation for the interface.
However, we see no
physical reason why the error should be significantly larger than in
the zero misfit case, and believe that  increasing the misfit
increases the cornering energy, as shown in Fig. \ref{f5}.  

There is another possibly important effect for the interface, however, 
associated with
the dependence of ${\cal G}_e$ on the length of the kink
when the elastic mismatch is nonzero. (See the discussion following
Eqn. (\ref{eqn5}).) Because of the length
dependence of the local $k$ for the kinking crack, the crack may
require a larger ${\cal G}_c$ to nucleate than to propagate.
This effect could masquerade for the increased cornering energy we
observe for the high elastic mismatch case.  We would not be able to tell
the difference between the oscillation in ${\cal G}_c$ at the growing
kink and an increase in cornering resisitance.  
In He and Hutchinson\cite{HH}, the oscillations in ${\cal G}_c$ are 
noted, but not
quantified, because they have no way to infer the actual phase at the
nucleating kink.  (They simply suggest that the phase be ignored,
which is not a good approximation in our results.  In our work,
we calculate it in terms of the force law range parameter.)
In our simulations, the kinking crack
always jumps forward, once it has been nucleated, which would be
consistent with an oscillatory ${\cal G}_c$ similar to that shown in Fig. 2
of their paper, but it is also consistent with the existence of a 
local cornering energy
for nucleating the kink. So, it is possible that the additional cornering
energy we observe for large elastic mismatch case could be an effect
due to oscillations in ${\cal G}_c$.  That there is a physical cornering
resistance, however,
is clear, because it occurs for zero elastic mismatch, where the phase
oscillation in ${\cal G}_c$ would disappear.  


Figs. \ref{f4},\ref{f5} represent the competition between cleavage on
the spur plane and emission on that plane.  But we noted in the
Introduction that there is also competition with events on the
cleavage plane to the right of the spur.  Depending on circumstances,
it is possible for the crack to either cleave or emit on that plane,
as well, as shown in Fig. \ref{f1}.  The total picture is obtained by
overlaying events on the interface plane with events on the spur, as
shown scematically in Fig. \ref{f6}. The
radius of the Griffith circle for cleavage on the interface is proportional to
$\sqrt{\gamma_s}$, and by changing the chemical bonding between the two
interfaces, this circle can be expanded or contracted, relative to
branching, as shown.  If the
events on the spur lie below the Griffith circle for cleavage
(or for nonblunting emission) on the interface plane, then
events on the spur plane will dominate over those on the interface
plane. Alternatively, if the interface bonding is weak, then events
on the spur plane can dominate.  
(The reader will now appreciate why it was
necessary to choose a strong force law on the interface to explore 
events on the spur.)  

In the Fig. \ref{f6}, the branching plane stability
line is shown crossing the smaller Griffith circle.  In this case,
events on the upper portion of the Griffith circle for the interface 
 will be stable
relative to branching, while for loads on the lower portion, branching
will be stable relative to events on the interface plane.
Such a degenerate situation can easily be set up in the
computer with precisely the results predicted by the diagram, but it
requires a careful choice of relative bonding between matrix and
interface. 

Note that ``crossover'' from the interface plane to the
spur plane takes place over a range of loads, but in the
particular situation depicted in 
Fig. \ref{f6}, there is no combination of loads where the crack
will emit a dislocation on the spur plane for the small Griffith
circle case.  For this particular choice of
force laws, then, the crack may cleave or emit on the interface plane, or
it may branch onto the spur plane as a cleavage crack.  But it
will not ever be possible to emit a blunting dislocation, and in this
sense, the material is intrinsically brittle.  However, there exists a
combination of force laws such that the situation depicted in
Fig. \ref{f7} is realized, and for a particular choice of loads, the
crack stability is degenerate---it can either emit or cleave,
corresponding to a crossover from brittle to ductile behavior.  In
previous papers\cite{RT,ZCT}, we have worked with pure Mode I loads,
and this is seen as a special case of the more general criterion for
crossover, where the loading mode is not specified. In Fig. \ref{f7},
zero elastic mismatch is assumed for clarity, and
the upper limit of the branching curve where emission
occurs lies precisely on the $K_I$ axis.  In general, of course, the
emission point does not lie on the $K_I$ axis.  But it is true that 
negative Mode II enhances cleavage on the branching plane,
so the upper limit of the branching stability line will tend to be 
in the vicinity of the $K_I$ axis, if not precisely on it. (This
statement, of course, does not apply to the interface case.)

Crossover from cleavage to emission can also be described in terms
of the ductility parameter, $\cal D$, defined as 
\begin{equation}
{\cal D}={{\cal G}_c\over{\cal G}_e}
\label{eqn11}
\end{equation}
When ${\cal D}>1$ the material is ductile, brittle otherwise.  Because
both ${\cal G}_c$ and ${\cal G}_e$ are quantities depending on the 
mixing of the modes, the
crossover condition defined in this way is not unique, 
but depends somewhat on the
mode mixity, $\psi$, as described graphically, above. 


\section{Crack Path in a Homogeneous Crystalline Material.}

In a separate paper\cite{T2}, we have shown that in an
amorphous material, where the surface energy is completely isotropic
in direction, a crack under mixed modes will deviate out of the
initial crack path into the path where $k_{II}=0$ on the branching
crack. In such cases, crack deviation is, of course, easy, and widely
observed.   But in a lattice,
where $\gamma_s$ has cusps on the high density planes, only certain
special planes are allowed for cleavage. In this paper, we have shown
that significant Mode II can be sustained by a crack, up to the limits
governed by the stability diagram.  Indeed, for any given lattice,
branching may be quite difficult, because not only must the Griffith
relation on the branching plane be satisfied (with any additional
cornering resistance which may be necessary),
but dislocation emission on one of the allowed planes might intervene 
before the relevant Griffith relation is satisfied.  

We performed a simulation to explore the case of branching in a
completely homogeneous lattice where all interfacial and chemical
effects were taken out.  For 
the UBER force law, we found that branching was possible only for a
range parameter corresponding to the largest possible unstable
stacking fault,  for the reason quoted above---sufficiently large
(negative) $K_{II}$ values were
necessary, such that dislocation emission tended to occur on the
cleavage plane. When branching was finally
induced, however, it became possible only for a very narrow,
essentially unique, load range, for which it was found that $k_{II}$
was not zero, but rather $k_{II}/k_I=.25$ when branching occurred.
We consider this finding to be a significant violation of the
expectation of He and Hutchinson\cite{HH} that $k_{II}=0$ at branching.
Presumably the reason for the violation is that the stability diagram
on the branching plane covers a range of $k_{II}$ values, and in the
UBER law, only the extreme lower limit of the branching stability
curve is reachable, and then only for a specific value of the range
parameter. 

Thus, in the lattice, in spite of the continuum
prediction that ${\cal G}_c$ is maximized when $k_{II}=0$, the
stability line on the branching plane encompasses a range of $k_{II}$
values. Therefore, a complete  analysis of the stability diagram 
for the given
force law, including competing dislocation emission on all possible
planes is necessary, in order to predict when
branching will occur.

Hutchinson Mear and Rice\cite{HMR} also apply the $k_{II}=0$ criterion
for predicting the plane on which cracking will occur when the crack
runs parallel to an interface.  Presumably similar difficulties will
arise in that case when the material is a crystal and not amorphous.
Any plane parallel to the interface where the crack satisfies the
stability criterion will be a possible plane of cracking.  Since
$k_{II}/k_I$ can reach values in the approximate range of 0.25 to 0.5
for reasonable force laws for stable cracks, and a crack, once
established on a stable plane will not deviate from it, there should
be a range of possible planes where a crack will be stable on a plane
parallel to an interface.



\section{Discussion.}

In this paper, we have presented a simple graphical construction in
the remote K--space of the crack, with which one can
describe the competition between emission and cleavage on two 
nonequivalent slip/cleavage planes.  In general, there
may be more planes where the lattice can either cleave or emit
dislocations.  These
more complex situations can be described by simply overlaying the
additional diagrams on the ones described here.  

For events  restricted to the initial cleavage plane, 
the stability diagram is a sector of a
Griffith circle, where the circle limits represent lattice breakdown
in shear and dislocation emission.  For those force laws where lattice
trapping is significant, then the Griffith circle fuzzes and develops a
finite width, within which the lattice is stable under the relevant loads.  

When events are possible on an inclined slip plane, then a similar
(noncircular) contour line for the critical ${\cal G}_c$ can be
constructed in the same space. Dislocation emission corresponds to
one limit on this curve, and the other limit corresponds to a critical
local Mode I value below which the crack is no longer able to
break the relevant bonds in the branching crack core.  
By overlaying these two
plots, one can determine which events are favored, and which are not
in a given load range.  Moreover, the crossover between simple
cleavage on the initial plane and blunting emission on an inclined
plane can be determined by varying the interfacial bond strength.
This crossover corresponds to the crossover between materials which
are intrinsically brittle, and those which are intrinsically ductile.
In this paper, we have used a simple lattice and simple force laws to
illustrate the physical principles involved.  In actual cases, one
must carry out calculations for realistic force laws and lattices. 

For the case of branching, the criterion is quite simple.  Branching
is expected to occur whenever the extended Griffith relation
(including the cornering resistance) with the bulk $\gamma_s$ in the
stability diagram drops below the Griffith relation
for simple straight crack extension with the interfacial $\gamma_s$.
Since negative Mode II translates into a strong ${\cal G}_c$ in the branching
plane, mixed modes on the interface strongly influence branching. The
magnitude of the cornering resistance can range from about 30\% of the
normal Griffith $\gamma_s$ to a factor of about 2 for the highest
elastic mismatch interface (where branching would be quite difficult).
However, in a crystal lattice, where large Mode II may be necessary in
order to induce branching, the competing dislocation emission may
intervene.  Thus, depending on the details of the force law,
increasing Mode II in a lattice may branch the crack, or it may break
it down in shear with dislocation emission.

It also follows that a crack which emits a
dislocation on an inclined glide plane will develop a large mode
mixity from the shielding by the emitted dislocation, 
and this could cause the crack to branch, as is sometimes seen
in the computer simulations of Holian, etal.\cite{HZ}.

We explored the prospects that branching will take place on the plane
where the local $k_{II}=0$.  While this criterion does hold for an
amorphous material\cite{T2}, we found that significant deviations from
this popular criterion take place for crack branching in a crystal
lattice. 





% now the references. delete or change fake bibitem. delete next three
%   lines and directly read in your .bbl file if you use bibtex.
\begin{references}
\bibitem{R} J. R. Rice, J. Mech. Phys. Solids, {\bf 40}, 239 (1992).
\bibitem{ZCT} S. J. Zhou, A. E. Carlsson, and R. Thomson,
Phys. Rev. Lett., {\bf 72}, 852 (1994).
\bibitem{TC} R. Thomson and A. E. Carlsson, Phil. Mag. {\bf 70}, 893
(1994).
\bibitem{T} R. Thomson, Submitted Phys. Rev.
\bibitem{BC} B. A. Bilby and G. E. Cardew, Int. J. Fract., {\bf 11},
708 (1975).
\bibitem{BCH} B. A. Bilby, G. E. Cardew and I. C. Howard, Fracture
1977, vol. 3, Univ. Waterloo Press, (1977) p 197.
\bibitem{DR} V. V. Dudukelanko and N. B. Romalis, Izv. An. SSR, {\bf
8}, 129 (1973).
\bibitem{L} K. K. Lo, J. App. Mech., {\bf 45}, 797 (1978).
\bibitem{HH} M.-Y. He and J. W. Hutchinson, J. Appl. Mech., {\bf 56}, 
270 (1989).
\bibitem{CR} B. Cotterell and J. R. Rice, Int. J. Fract., {\bf 16},
155 (1980).
\bibitem{AR} P. Anderson and R. Thomson, J. Appl. Phys., {\bf 76}, 843 (1994).
\bibitem{T1} R. Thomson, Solid State Physics, {\bf 39}, 1 (1986),
ed. D. Turnbull and H. Ehrenreich, Academic Press, New York.
\bibitem{ZCT2} S. J. Zhou, A. E. Carlsson and R. Thomson, Phys. Rev. B,
{\bf 47B}, 7710 (1993).
\bibitem{TZCT} R. Thomson, S. J. Zhou, A. E. Carlsson, and
V. K. Tewary, Phys. Rev. B, {\bf 46B}, 10613 (1992). 
\bibitem{RSW} J. R. Rice, Z. Suo and J.-S. Wang, Metal--ceramics 
Interfaces, ed. M. Rulhe, A. G. Evans, M. F. Ashby and J. Hirth, 
Acta Scripta Met. Proc. Series {\bf 4}, 269 (1990)
\bibitem{MS} B. M. Malyshev and R. L. Salganik, Int. J. Fract. Mech.,
{\bf 5}, 261 (1965).
\bibitem{RSF} J. Rose, J. Smith and J Ferrante, Phys. Rev., {\bf 28},
1835 (1983).
\bibitem{ZT} S. J. Zhou and R. Thomson, Phys. Rev. B, {\bf 49B}, 44
(1994).
\bibitem{RT} J. R. Rice and R. Thomson, Philos. Mag., {\bf 29}, 73
(1974).
\bibitem{T2} R. Thomson, to be published.
\bibitem{HMR} J. W. Hutchinson, M. E. Mear, and J. R. Rice,
J. Appl. Mech., {\bf 54}, 828, (1987).
\bibitem{HZ} B. Holian and S. J. Zhou, Unpublished report.
\end{references}


% figures follow here
%
% Here is an example of the general form of a figure:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
%

\begin{figure}
\caption{Stability diagram.  (For interpretation of the simple
homogeneous lattice, consider only the axes $k_I/k_{II}$.) Any load
satisfying the Griffith condition will lie on the circle.  Dislocation
emsission corresponds to the limits of the curve in the $k_{II}$
direction, where the lattice breaks down in shear.  The locus of
points on the Griffith line between these two limiting points
corresponds to possible stable loadings of the lattice. When an interface is
present, the criteria for cleavage or emission are given in the local
core stress intensity factors ($k_I/k_{II}$ axes), which are rotated 
from the lab axes ($K_I/K_{II}$ axes) by the core phase angle.}
\label{f1}
\end{figure}

\begin{figure}
\caption{Contour plot of $\cal G$ for a branching crack.  When the
crack branches, then the simple circular Griffith line becomes a more
general form, as shown for two cases.  In (a), the
interface mismatch is zero, while in (b), $c_1/c_2=10$.  A stability
diagram would correspond to a segment of a critical $\cal G$ contour
between allowed values of $k_{II}$.}
\label{f2}
\end{figure}

\begin{figure}
\caption{The interfacial crack in a hexagonal lattice.  The atoms
below the centerline lie in sublattice 1 and those above the line in
sublattice 2.  The crack line lies on the interface between the two.
Sublattice 1
has spring constant $c_1$, sublattice 2 has spring constant $c_2=1$,
and the interfacial bonds have a third spring constant, $c_{12}$.
Nonlinear bonds are formed at all atoms on the slip plane
or the cleavage plane in the cohesive zone of the crack.  Different
nonlinear bonds can form on lower atoms, upper atoms (including the
slip plane), and between atoms facing one another across the
interface. The nonlinear bonds are depicted by the zig-zag lines.}
\label{f3}
\end{figure}

\begin{figure}
\caption{Simulations plotted as points on the stability diagram for
zero mismatch between sublattices ($c_1=c_2$).  The
solid curve corresponds to the ${\cal G}_c=\gamma_s$ in material 2.
If the Griffith relation were satisfied, then the points would lie on
the curve.  That they lie parallel to it corresponds to the fact that
branching the crack off its initial cleavage plane requires an
additional energy from that for the straight crack, which we term a
cornering energy.  The upper limit for the simulation points
corresponds to dislocation emission, and the lower limit corresponds
to the critical driving force required to break the first bond on the
spur plane.}

\label{f4}
\end{figure}

\begin{figure}
\caption{Same as previous figure, except $c_1=10c_2$.}
\label{f5}
\end{figure}

\begin{figure}
\caption{Superposition of stability diagrams for initial plane and
inclined plane for remote stress intensity factors, $K_I$ and
$K_{II}$.   Two Griffith circles are drawn with different radii,
corresponding to different $\gamma_s$ values on the interface plane.  
The third curve corresponds to the stability line on the branch plane.
In this case, only the upper limit of the curve corresponds to
dislocation emission. In the case shown, the branching events are all
more stable than any events on the large Griffith circle.  However, the
branch line crosses the small Griffith circle, and the upper part of
the small Griffith circle is more stable than the upper part of the
branch line, and vice versa for the lower parts of those curves.  Note
that for the case shown, there is no dislocation emission on the
branch plane.}
\label{f6}
\end{figure}

\begin{figure}
\caption{Same as above, for ductility crossover.  Here, blunting
emission is degenerate with the cleavage Griffith circle.  This choice
of force laws thus corresponds to the crossover between ductile and
brittle behavior for the material.  This figure would correspond to a
homogeneous lattice without interface.  It is not expected that
blunting emission should always lie exactly on the $K_I$ axis, 
but it probably does lie close to that axis in most cases.}
\label{f7}
\end{figure}



% tables follow here
%
% Here is an example of the general form of a table:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Insert the column specifiers (l, r, c, d, etc.) in the empty braces of the
% \begin{tabular}{} command.
%
% \begin{table}
% \caption{}
% \label{}
% \begin{tabular}{}
% \end{tabular}
% \end{table}

\end{document}
%
% ****** End of file template.aps ******
