# A comprehensive study of nonlinear perturbations in the dynamics of planar
crack fronts

Itamar Kolvin  Mokhtar Adda-Bedia  School of Physics, Georgia Institute of
Technology, 837 State St NW, Atlanta, 30332, GA, USA  Laboratoire de Physique,
CNRS, ENS de Lyon, Université de Lyon, F-69342 Lyon, France

###### Abstract

The interaction of crack fronts with asperities is central to the criteria of
fracture in heterogeneous materials and for predicting fracture surface
formation. It is known how dynamic crack fronts respond to small, 1st-order,
perturbations. However, large and localized disturbances to crack motion
induce dynamic and geometric nonlinear effects beyond the existing linear
theories. Because the determination of the 3D elastic fields surrounding
perturbed crack fronts is a necessary step towards any theoretical study of
crack front dynamics, we develop a 2nd-order perturbation theory for the
asymptotic fields of planar crack fronts. Based on previous work, we consider
two models of fracture. In the so-called scalar elastic model, which is
analogous to anti-plane (Mode III) fracture, the stress and displacement
fields are obtained through matched asymptotic expansions. A self-consistent
expansion is used to resolve the fields near tensile (Mode I) crack fronts.
Both methods can be extended to higher perturbation orders. The main results
of this work are the explicit 2nd-order expressions of the local dynamic
energy-release-rates for perturbations of straight fronts. These general
formulae recover the known energy-release-rates of curved quasi-static fronts
and of straight dynamic fronts.

##  1 Introduction

The prediction of crack propagation in heterogeneous media is a central
problem in fracture mechanics (Ponson, 2009; Démery et al., 2014; Steinhardt
and Rubinstein, 2022; Schmittbuhl and Måløy, 1997; Ponson, 2009; Chopin et
al., 2015; Lebihain et al., 2020; Albertini et al., 2021; Roch et al., 2023;
Stanchits et al., 2015; Lubomirsky and Bouchbinder, 2023), in designing
advanced materials (Xia et al., 2012; Mirkhalaf et al., 2014), and in
frictional and earthquake mechanics (Latour et al., 2013; Lebihain et al.,
2021; Gounon et al., 2022; Svetlizky and Fineberg, 2014; Bayart et al., 2018;
Bedford et al., 2022). Computational frameworks for 3D dynamic fracture, such
as the spectral boundary integral method (Geubelle and Rice, 1995; Roch et
al., 2022, 2023), phase-field simulations (Pons and Karma, 2010; Henry and
Adda-Bedia, 2013; Chen et al., 2015; Bleyer and Molinari, 2017; Henry, 2019;
Lubomirsky and Bouchbinder, 2023), and atomistic models (Heizler and Kessler,
2015; Möller and Bitzek, 2015), perform the intensive task of numerically
determining the elastodynamic fields. A complementary approach, pursued in
this work, reduces the complexity of the 3D problem to the motion of 1D crack
fronts. Cracks in brittle materials are described by the theory of Linear
Elastic Fracture Mechanics (LEFM). The cornerstone of LEFM is that the near-
crack region is dominated by a universal $r^{-1/2}$ stress field characterized
by a single intensity factor. The stress intensity factor (SIF), then,
determines the crack dynamics. This approach was adopted in several works
(Rice et al., 1994; Willis and Movchan, 1995; Willis and Movchan, 1997;
Ramanathan and Fisher, 1997; Morrissey and Rice, 1998, 2000; Norris and
Abrahams, 2007; Willis, 2013; Adda-Bedia et al., 2013) which considered small
perturbations to straight crack fronts. The usefulness of this approach was
demonstrated by the discovery of crack front waves (Ramanathan and Fisher,
1997; Morrissey and Rice, 1998, 2000), which were later observed
experimentally (Sharon et al., 2001, 2002; Fineberg et al., 2003; Sagy et al.,
2004). When applied to mixed mode loading configurations, this approach
yielded possible instability mechanisms responsible for the generation of
corrugation waves (Adda-Bedia et al., 2013) and crack front segmentation (Pons
and Karma, 2010; Leblond et al., 2011; Chen et al., 2015). For heterogeneous
materials, in which the toughness contrasts are of order unity, however, the
linear theory is of little use. Understanding nonlinear perturbations is also
needed to determine how heterogeneity affects energy dissipation in fracture,
since 1st-order perturbations have a net zero contribution to dissipation. The
importance of nonlinear effects has been sketched in a recent work (Kolvin et
al., 2017). There, a first-principles nonlinear equation of motion for crack
fronts was derived in the context of scalar elasticity and their response to
externally induced perturbations was computed. It was suggested that nonlinear
front focusing coupled with the rate-dependence of fracture energy dissipation
may govern the transition to micro-branching.

In the framework of LEFM, the modeling of the crack front dynamics
necessitates the knowledge of the variation of dynamic SIF with the crack
front geometry. The mathematical foundations of nonlinear perturbations of
static crack fronts are well-developed (Leblond et al., 2012; Vasoya et al.,
2016; Adda-Bedia et al., 2006). While the general methods for dynamic cracks
were developed by Movchan et al. (1998), Ramanathan and Fisher (1997);
Ramanathan (1997),Norris and Abrahams (2007) and Willis (2013), these works
did not provide explicit, tractable formulas that can be applied to concrete
crack evolution. Morrissey and Rice (2000), however, obtained explicit 1st-
order formula for the energy-release-rate of Mode I crack fronts, and
numerically studied crack propagation in hetereogeneous solids. Based on these
works, we compute the 2nd-order corrections to the energy-release-rate in the
scalar model of elasticity and conventional “vectorial” elastodynamics. The
energy balance between the elastic energy flux into the crack front and the
dissipation can then be used to obtain equations of motion for dynamic planar
crack fronts.

Figure 1: Geometry of the problem. The crack front (red) subtends two semi-
infinite planar crack faces (blue). The front propagates in the $y=0$ plane of
a linear elastic solid at a velocity $V$ with a superimposed spatiotemporal
perturbation $f(z,t)$.

Consider an infinite body containing a semi-infinite planar crack that is
driven by remote loading (Fig. 1). The leading edge of the propagating crack
front experiences perturbations in space and time around a steady motion at a
constant velocity $V$. The coordinate system is defined as follows: $x$ is the
crack propagation direction, the $y$ axis is perpendicular to the fracture
$xz$ plane, and $t$ signifies time. The instantaneous crack front position is
$h(z,t)=Vt+f(z,t)$, where $f$ is the crack front perturbation. The task at
hand is to obtain an expression for the local energy-release-rate $G(z,t)$
along the front. Two different methods are used: matched asymptotic
expansions, in the case of scalar elasticity, and a self-consistent expansion
in the case of Mode I fracture. Both methods can be generalized to higher-
order perturbations.

In the following, we outline the main results by providing the explicit
expressions of the 2nd-order perturbations of the energy-release-rate.
Afterward, we detail, for the interested reader, the calculations of both
scalar and tensile crack problems.

##  2 Summary of the results

###  2.1 Derivation of G in the scalar elasticity model

In section 3, we solve for the asymptotic displacement and stress fields
surrounding a crack in a model scalar elastic solid, extending the
calculations of Norris and Abrahams (2007). Fracture in scalar elasticity is
identical to anti-plane shear (Mode III) fracture in a solid body that is
constrained to displace along one dimension only, which is described by the
displacement $\mathbf{u}=(0,0,u_{z})$. The theory is characterized by a single
wave speed ($b=1$). We assume that the crack is driven at velocity $V$ by
remote loading conditions that are characterized by a distance $l$ from the
crack front. The elastic fields are expanded in the small parameter
$\epsilon\equiv||f||/l\ll 1$ in two separate perturbation series: the terms of
the first series are centered at the perturbed crack front $x=h(z,t)$; the
terms of the second series are centered at the unperturbed, imaginary,
position of the crack front $x=Vt$. Solutions of the equation of elasticity
contain unknown integration constants that are determined by matching the two
expansions. The main result of this calculation is an expression for $G$ that
is exact up to $\order{f^{2}}$ and $\order{\epsilon^{2}}$. When
$\epsilon\rightarrow 0$, the expansion of the energy-release-rate becomes

$G(z,t)=G_{r}g(V_{\perp})\left(1+H^{(1)}[f]+H^{(2)}[f,f]+\order{f^{3}}\right)\,,$
(1)

where the rest energy-release-rate $G_{r}$ is determined by the loading
conditions, $g(V_{\perp})$ is the dynamical contribution given by (Rice et
al., 1994),

$g(V_{\perp})=\sqrt{\frac{1-V_{\perp}}{1+V_{\perp}}}\,,$ (2)

and $V_{\perp}(z,t)=(V+f_{t})/\sqrt{1+f_{z}^{2}}$ is the local normal front
velocity of the crack front. Subscripted functions $f_{t},\,f_{z}$ denote
partial derivatives. The linear functional $H^{(1)}[f]$ and the 2nd-order
functional $H^{(2)}[f,f]$ depend only on the crack front history
$\\{f(z,t^{\prime}<t)\\}$,

$\displaystyle H^{(1)}[f]$ $\displaystyle=$
$\displaystyle-\frac{1}{\alpha^{2}}\Psi[f]\,,$ (3) $\displaystyle H^{(2)}[f]$
$\displaystyle=$
$\displaystyle\frac{1}{4\alpha^{4}}\Psi[f]^{2}+\frac{1}{2\alpha^{4}}\Psi[f\Psi[f]]-\frac{1-2V}{4\alpha^{4}}\Psi_{2}[f^{2}]-\frac{1+2V}{2\alpha^{4}}f\Psi_{2}[f]\,,$
(4)

where $\alpha=\sqrt{1-V^{2}}$, and, with the definition of the Fourier
transform $\bar{f}(k,t)=\int\\!\mathrm{d}z\,e^{-ikz}f(z,t)$,

$\displaystyle\overline{\Psi[f]}(k,t)$ $\displaystyle=$ $\displaystyle\alpha
k\int_{-\infty}^{t}\\!\mathrm{d}t^{\prime}\,\frac{J_{1}(\alpha
k(t-t^{\prime}))}{t-t^{\prime}}\bar{f}(k,t^{\prime})\;,$ (5)
$\displaystyle\overline{\Psi_{2}[f]}(k,t)$ $\displaystyle=$
$\displaystyle\alpha^{2}k^{2}\int_{-\infty}^{t}\\!\mathrm{d}t^{\prime}\frac{J_{2}(\alpha
k(t-t^{\prime}))}{t-t^{\prime}}\bar{f}(k,t^{\prime})\;,$ (6)

where $J_{1}$ and $J_{2}$ are the 1st\- and 2nd-order Bessel functions.

####  2.1.1 Recovery of known limits

To test Eq. (1) against known results, we consider two limiting cases: that of
a straight front ($f_{z}=0$), and that of a steadily propagating front
($f_{t}=0$).

The condition $f_{z}=0$ is equivalent to taking $k\rightarrow 0$ in Eqs.
(5-6), which approach $\Psi[f],\Psi_{2}[f]\rightarrow 0$. Then the energy-
release-rate, Eq. (1), approaches $G\rightarrow G_{r}g(V+f_{t})$, which is
identical to the known equation of motion of a 2D crack (Rice et al., 1994).

When $f_{t}=0$, Eqs. (5-6) become $\overline{\Psi[f]}=\alpha|k|\bar{f}(k)$ and
$\overline{\Psi_{2}[f]}=\frac{1}{2}\alpha^{2}k^{2}\bar{f}(k)$. Then

$\displaystyle\overline{H^{(1)}[f]}=-\frac{1}{\alpha}|k|\bar{f}(k)\,,$
$\displaystyle\overline{H^{(2)}[f,f]}=$ (7)
$\displaystyle\frac{1}{4\alpha^{2}}\int\\!\mathrm{d}k^{\prime}\,\bar{f}(k-k^{\prime})\bar{f}(k^{\prime})\left(|k^{\prime}||k-k^{\prime}|+|k|(|k^{\prime}|+|k-k^{\prime}|)-k^{2}+(1+2V)k^{\prime}(k-k^{\prime})\right)\,.$

Using these results and expanding $g(V_{\perp})\simeq
g(V)-\frac{1}{2}Vg^{\prime}(V)f_{z}^{2}$, Eq. (1) becomes

$\begin{split}\overline{G}(k)&=G_{r}g(V)\Bigg{\\{}2\pi\delta(k)-\frac{1}{\alpha}|k|\bar{f}(k)\\\
&+\frac{1}{4\alpha^{2}}\int\\!\mathrm{d}k^{\prime}\,\bar{f}(k-k^{\prime})\bar{f}(k^{\prime})\left(|k^{\prime}||k-k^{\prime}|+|k|(|k^{\prime}|+|k-k^{\prime}|)-k^{2}+k^{\prime}(k-k^{\prime})\right)\Bigg{\\}}\,,\end{split}$
(8)

after symmetrization $k^{\prime}\leftrightarrow k-k^{\prime}$. In the limit
$V\rightarrow 0$, Eq. (8) becomes

$\begin{split}\overline{G}(k)&=G_{r}\Bigg{\\{}2\pi\delta(k)-|k|\bar{f}(k)\\\
&+\frac{1}{4}\int\\!\mathrm{d}k^{\prime}\,\bar{f}(k-k^{\prime})\bar{f}(k^{\prime})\left(|k^{\prime}||k-k^{\prime}|+|k|(|k^{\prime}|+|k-k^{\prime}|)-k^{2}+k^{\prime}(k-k^{\prime})\right)\Bigg{\\}}\,.\end{split}$
(9)

The same expression is obtained from the 2nd-order approximation of the SIF of
quasi-static crack fronts (Leblond et al., 2012; Vasoya et al., 2016)).

###  2.2 Derivation of G in Mode I fracture

In section 4, we determine the asymptotic fields in the vicinity of a crack
that propagates under remote tensile loading conditions in an infinite body of
Young modulus $E$ and Poisson ratio $\nu$, extending the calculations of
Ramanathan (1997). In the absence of perturbations, the loading conditions
produce an asymptotic Mode I stress field. Ahead of the crack front and at the
fracture plane, the tensile stress is $\sigma_{yy}\sim
K_{0}^{\sigma}/\sqrt{2\pi(x-Vt)}$ with the SIF,
$K_{0}^{\sigma}=K_{r}^{\sigma}k(V)$. $K_{r}^{\sigma}$ is the rest SIF and
$k(V)$ is a universal function of the crack velocity defined in Eq. (6.4.26)
of Freund (1990). The energy-release-rate is then $G=G_{r}g(V)$ where
$G_{r}=\frac{1-\nu^{2}}{E}(K_{r}^{\sigma})^{2}$ is the rest energy-release-
rate, the dynamical prefactor is $g(V)=A_{I}(V)k(V)^{2}$, and
$A_{I}(V)=(1-\nu)^{-1}(V^{2}\alpha_{a})/(R(V)b^{2})$;
$R(V)=4\alpha_{a}\alpha_{b}-(1+\alpha_{b}^{2})^{2}$;
$\alpha_{s}=\sqrt{1-(V/s)^{2}}$; and $s=a,b$ are the longitudinal and the
shear wave speeds (Freund, 1990).

For space-time-dependent crack fronts, we derive a perturbation expansion of
the energy-release-rate

$G(z,t)=G_{r}g(V_{\perp})\left(1+H^{(1)}[f]+H^{(2)}[f,f]+\order{f^{3}}\right)$
(10)

Explicit expressions for the history functionals $H^{(1)}[f]$ and
$H^{(2)}[f,f]$ are found in the space $(k,t)$, which defined via the Fourier
transform $\bar{f}=\int\\!\mathrm{d}z\,e^{-ikz}f(z,t)$. The linear functional
is $\overline{H^{(1)}[f]}=2I_{1}[\bar{f}]$ where
$I_{1}[\bar{f}]=-\int_{-\infty}^{t}\\!\mathrm{d}t^{\prime}\,A_{1}(k,t-t^{\prime})\bar{f}(k,t^{\prime})$
and (Morrissey and Rice, 2000)

$A_{1}(k,t)=k^{2}\left[-\frac{a}{2}\frac{J_{1}(\beta_{a}t)}{\beta_{a}t}+c\frac{J_{1}(\beta_{c}t)}{\beta_{c}t}-\frac{1}{4}\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\left(\frac{\eta^{2}+V^{2}}{\eta^{2}-V^{2}}J_{2}(\beta_{\eta}t)-J_{0}(\beta_{\eta}t)\right)\right]\,,$
(11)

where $\beta_{s}=\sqrt{s^{2}-V^{2}}|k|$. The 2nd-order functional is

$\begin{split}\overline{H^{(2)}[f,f]}&=2I_{1}[\bar{f}*I_{1}[\bar{f}]]-\frac{1}{2}I_{1}[I_{1}[\bar{f}*\bar{f}]]-\bar{f}*I_{1}[I_{1}[\bar{f}]]+I_{1}[\bar{f}]*I_{1}[\bar{f}]\\\
&+\frac{V}{2}I_{2}[\bar{f}*\bar{f}]-V\bar{f}*I_{2}[\bar{f}]\,.\end{split}$

where the convolution operator is defined as
$f*g=\int\\!\mathrm{d}k^{\prime}\,f(k-k^{\prime})g(k^{\prime})$, and
$I_{2}[\bar{f}]=-\int_{-\infty}^{t}\\!\mathrm{d}t^{\prime}\,A_{2}(k,t-t^{\prime})\bar{f}(k,t^{\prime})$
with

$A_{2}(k,t)=|k|^{3}\Bigg{[}-\gamma_{a}\frac{J_{2}(\beta_{a}t)}{\beta_{a}t}+2\gamma_{c}\frac{J_{2}(\beta_{c}t)}{\beta_{c}t}-\frac{1}{4}\int_{b}^{a}\\!\mathrm{d}\eta\,\frac{\Theta(\eta)}{\sqrt{\eta^{2}-V^{2}}}\left(\frac{3\eta^{2}+V^{2}}{\eta^{2}-V^{2}}J_{3}(\beta_{\eta}t)-J_{1}(\beta_{\eta}t)\right)\Bigg{]}\,.$

It is possible to use these expressions to evolve crack fronts in time
(Morrissey and Rice, 2000).

It is useful to write $G$ in Fourier space $(k,\omega)$, where
$\hat{f}(k,\omega)=\int\\!\mathrm{d}k\mathrm{d}\omega\,e^{-i\omega
t-ikz}f(z,t)$. We define the convolution
$\hat{f}\otimes\hat{g}=\int\\!\mathrm{d}k^{\prime}\mathrm{d}\omega^{\prime}\,\hat{f}(k-k^{\prime},\omega-\omega^{\prime})\hat{f}(k^{\prime},\omega^{\prime})$.
The perturbation expansion of the energy-release-rate is then

$\widehat{G}(k,\omega)=G_{r}g(V)\left((2\pi)^{2}\delta(k)\delta(\omega)+\widehat{\delta
G^{(1)}}(k,\omega)+\widehat{\delta G^{(2)}}(k,\omega)+\order{f^{3}}\right)\,.$
(12)

The 1st-order correction is (Ramanathan and Fisher, 1997)

$\widehat{\delta
G^{(1)}}(k,\omega)=-2|k|P_{1}(\omega/|k|;V,\nu)\hat{f}(k,\omega)\,,$ (13)

and the 2nd-order correction is

$\begin{split}&\widehat{\delta
G^{(2)}}(k,\omega)=2|k|P_{1}\left\\{\hat{f}\otimes\left(|k|P_{1}\hat{f}\right)\right\\}-\left\\{\frac{1}{2}k^{2}P_{1}^{2}+\frac{i}{2}V\omega|k|P_{2}\right\\}(\hat{f}\otimes\hat{f})\\\
&-\hat{f}\otimes\left\\{\left(k^{2}P_{1}^{2}-iV\omega|k|P_{2}\right)\hat{f}\right\\}+(|k|P_{1}\hat{f})\otimes(|k|P_{1}\hat{f})\,.\end{split}$
(14)

The kernel functions $P_{1}$ and $P_{2}$ are

$\displaystyle P_{1}(s;V,\nu)=$
$\displaystyle-\frac{1}{2}\gamma_{a}\sqrt{1-\frac{\gamma_{a}^{2}s^{2}}{a^{2}}}+\gamma_{c}\sqrt{1-\frac{\gamma_{c}^{2}s^{2}}{c^{2}}}$
$\displaystyle+\frac{1}{2}\int_{b}^{a}\frac{(s^{2}+V^{2})(\eta^{2}+V^{2})-2\eta^{2}V^{2}}{\sqrt{\eta^{2}-(s^{2}+V^{2})}\left(\eta^{2}-V^{2}\right)^{2}}\Theta(\eta)\,\mathrm{d}\eta\,,$
(15)

and

$\displaystyle P_{2}(s;V,\nu)=$ $\displaystyle
2\frac{\gamma_{c}^{3}}{c^{2}}\sqrt{1-\frac{\gamma_{c}^{2}s^{2}}{c^{2}}}-\frac{\gamma_{a}^{3}}{a^{2}}\sqrt{1-\frac{\gamma_{a}^{2}s^{2}}{a^{2}}}+$
$\displaystyle\int_{b}^{a}\frac{(s^{2}+V^{2})(3\eta^{2}+V^{2})-2\eta^{2}(\eta^{2}+V^{2})}{\sqrt{\eta^{2}-(s^{2}+V^{2})}\left(\eta^{2}-V^{2}\right)^{3}}\Theta(\eta)\,\mathrm{d}\eta\,,$
(16)

where $\gamma_{s}=1/\sqrt{1-V^{2}/s^{2}}$, $c$ is the Rayleigh wave speed, and

$\Theta(\eta)=\frac{2}{\pi}\arctan\left[4\sqrt{1-\frac{\eta^{2}}{a^{2}}}\sqrt{\frac{\eta^{2}}{b^{2}}-1}\bigg{/}\left(2-\frac{\eta^{2}}{b^{2}}\right)^{2}\right]\,.$
(17)

The branch cuts of the square roots in Eqs. (2.2,2.2) are defined such that
$\sqrt{1-s^{2}}=i\,\mathrm{sign}(s)\sqrt{s^{2}-1}$ for $|s|>1$. Example curves
of the functions $P_{1}$ and $P_{2}$ are depicted in Fig. 2. Particularly, the
zeros $s_{1}$ and $s_{2}$ of $P_{1}$ and $P_{2}$ do not coincide (Fig. 3). The
roots $s=\pm s_{1}$ of $P_{1}$ signify the existence crack front waves, since
$\delta G^{(1)}=0$ at $\omega/|k|=\pm s_{1}$. Because $s_{1}\neq s_{2}$,
$\delta G^{(2)}\neq 0$ at $\omega/|k|=\pm s_{1}$. Crack front waves are
expected, therefore, to have a non-zero contribution to the energy-release-
rate at the 2nd-order.

Figure 2: Real and imaginary parts of $P_{1}(s;V,\nu)$ (left) and
$P_{2}(s;V,\nu)$ (right) for $V=0.5b$ and $\nu=0.3$.

Figure 3: The real roots $s_{1}$ of $P_{1}$ and $s_{2}$ of $P_{2}$ normalized
by $\sqrt{c^{2}-V^{2}}$. $\nu=0.3$

####  2.2.1 Recovery of known limits

Eq. (10) can be reduced to known expressions at the 2D limit $f=f(t)$ and at
the quasi-static limit $f=f(z)$. In the 2D limit, $H^{(1)}[f]=0$ and
$H^{(2)}[f,f]=0$ since $\bar{f}\propto\delta(k)$. The 2D expression for the
energy-release-rate $G=G_{r}g(V+f_{t})$ is recovered.

In the quasi-static limit, $f=f(z)$. Then, $I_{1}[\bar{f}]=-P_{1}(0;V)|k|$ and
$I_{2}[\bar{f}]=\pi_{1}k^{2}$ where

$\pi_{1}=\frac{1}{2}\frac{a}{a^{2}-V^{2}}-\frac{c}{c^{2}-V^{2}}+\frac{1}{2}\int_{b}^{a}\frac{\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{2}}\Theta(\eta)\mathrm{d}\eta\,.$
(18)

The history functionals become

$\overline{H^{(1)}[f]}=-2P_{1}(0;V)|k|\bar{f}(k).$

and

$\begin{split}\overline{H^{(2)}[f]}=&P_{1}(0;V)^{2}\int\\!\mathrm{d}k^{\prime}\,\bar{f}(k-k^{\prime})\bar{f}(k^{\prime})\left(|k|(|k^{\prime}|+|k-k^{\prime}|)-k^{2}+k^{\prime}(k-k^{\prime})+|k^{\prime}||k-k^{\prime}|\right)\\\
&+V\pi_{1}(k\bar{f})*(k\bar{f})\,.\end{split}$

Since $2\pi_{1}=g^{\prime}(V)/g(V)$, the energy-release-rate expansion becomes

$\begin{split}\bar{G}(k)&=G_{r}g(V)\Bigg{\\{}2\pi\delta(k)-2P_{1}(0;V)|k|\bar{f}(k)\\\
&+P_{1}(0;V)^{2}\int\\!\mathrm{d}k^{\prime}\,\bar{f}(k-k^{\prime})\bar{f}(k^{\prime})\left(|k|(|k^{\prime}|+|k-k^{\prime}|)-k^{2}+k^{\prime}(k-k^{\prime})+|k^{\prime}||k-k^{\prime}|\right)+\order{f^{3}}\Bigg{\\}}\,.\end{split}$
(19)

At the limit $V\rightarrow 0$, Eq. (9) is reproduced.

##  3 Nonlinear perturbation of scalar elastic crack fronts

This section develops a perturbation theory for planar crack fronts that
propagate in a solid whose are described by a scalar displacement field $u$.
In this model, the stress, a vector field, is defined by
$\mathbf{\sigma}=\mu\nabla u$ where $\mu$ is the elastic modulus. This
formulation is identical to the elastic problem of anti-plane shear
deformation when the orthogonal displacement components are artificially set
to zero (Rice et al., 1994). In general, freely propagating anti-plane shear
(Mode III) cracks experience tensile (Mode I) and in-plane shear (Mode II)
stresses when perturbed (Geubelle and Rice, 1995). The scalar model, however,
is useful in obtaining a qualitative understanding of dynamic fracture, as it
is much simpler than “full” elastodynamics (Perrin and Rice, 1994).

Displacements are governed by Newton’s 2nd law

$\rho\partial_{t}^{2}u=\mu\nabla\cdot\mathbf{\sigma}$

where $\rho$ is the density. With the definition of the elastic wave velocity
$b=\sqrt{\mu/\rho}$, the displacement satisfies the wave equation

$\nabla^{2}u-\frac{1}{b^{2}}\partial_{t}^{2}u=0\,.$ (20)

In the following we set $b=1$ and we assume that the crack is driven at
velocity $V$ by stresses that are applied at a large distance $l$ from the
crack front. It is helpful to transform the wave equation to the co-moving
frame, $x\rightarrow(x-Vt)/l,\,y\rightarrow y/l,\,z\rightarrow
z/l,\,t/l\rightarrow t,\,\sigma\rightarrow\sigma/\mu$, such that

$\alpha^{2}u_{xx}+u_{yy}+u_{zz}-u_{tt}+2Vu_{xt}=0\,,$ (21)

where subscripts denote partial derivatives. The equation of motion is
supplemented by the boundary conditions at the fracture plane $y=0$,

$u(x>f(z,t),y=0,z,t)=0;\;\;\sigma_{y}(x<f(z,t),y=0,z,t)=0.$ (22)

For the simple crack, $f(z,t)=0$, the asymptotic elastic fields that solve
this problem are well-known (Eshelby, 1969; Rice et al., 1994; Norris and
Abrahams, 2007). The displacement can be written as a Williams expansion

$u(x,y,z,t)=\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\,\imaginary\left\\{\sqrt{x-Vt+i\alpha
y}\left[1+\frac{m}{3}(x-Vt+i\alpha y)+\frac{n}{5}(x-Vt+i\alpha
y)^{2}+...\right]\right\\}$ (23)

where $i=\sqrt{-1}$, $\alpha=\sqrt{1-V^{2}}$ and
$K^{u}_{0}=K^{\sigma}_{0}/\alpha$. The SIF $K^{\sigma}_{0}$ and the
coefficients $m\propto l^{-1}$ and $n\propto l^{-2}$ are determined by the
loading conditions. The SIF can be written as a product
$K^{\sigma}_{0}=K_{r}k(V)$ where $K_{r}$ is the rest SIF and $k(V)=\sqrt{1-V}$
(Rice et al., 1994). The energy-release-rate for the straight crack is then
$G_{0}=G_{r}g(V)$ where $G_{r}\propto K_{r}^{2}$ is the rest energy release
and $g(V)=\sqrt{(1-V)/(1+V)}$.

The simple crack solution provides a basis for exploring the effect of
arbitrary perturbations. Below, we follow the method of Norris and Abrahams
(2007) where the elastic fields are determined by matching two asymptotic
expansions. In the “inner” expansion the fields are centered at the crack
front position $x=f(z,t)$. In the “outer” expansion, the fields are centered
at $x=0$. The coefficients of the inner and outer solutions are then matched
to yield the displacement and the stress intensity factors, defined as

$\displaystyle K^{u}$ $\displaystyle=$ $\displaystyle\lim_{x-f(z,t)\rightarrow
0^{-}}(x-f(z,t))^{-1/2}\sqrt{\frac{\pi}{2}}\mu u(x,y=0,z,t)$ (24)
$\displaystyle K^{\sigma}$ $\displaystyle=$
$\displaystyle\lim_{x-f(z,t)\rightarrow
0^{+}}\sqrt{2\pi(x-f(z,t))}\sigma_{y}(x,y=0,z,t)\,.$ (25)

The energy release is then $G\propto K^{u}K^{\sigma}$.

###  3.1 The inner solution

We define $\epsilon=\mathrm{max}_{z,t}|f(z,t)|/l$, and the “inner” variables
$X=\epsilon x,Y=\epsilon y$. Writing the displacement field as a function of
the inner variables $U(X,Y,z;t)=u(x,y,z;t)$, Eq. (21) becomes

$\alpha^{2}U_{XX}+U_{YY}+2\epsilon VU_{Xt}+\epsilon^{2}(U_{zz}-U_{tt})=0\,.$
(26)

Expanding
$U=\epsilon^{1/2}U^{(1/2)}+\epsilon^{3/2}U^{(3/2)}+\epsilon^{5/2}U^{(5/2)}+...$
and substituting the expansion in Eq. (26), $U$ is found order-by-order in
$\epsilon$. Since we are seeking 2nd-order corrections to $G$, we will need to
compute the expansion up to the $\epsilon^{5/2}$ term. The zeroth order
solution corresponds to the dominant term in the simple crack field

$U^{(1/2)}=\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\Im{S^{1/2}}\,,$ (27)

where $S=X-f+i\alpha Y$.

To find the next order terms, we will use the identity

$(\alpha^{2}\partial_{X}^{2}+\partial_{Y}^{2})w_{1}(\overline{S})w_{2}(S)=4\alpha^{2}w_{1}^{\prime}(\overline{S})w_{2}^{\prime}(S)\,,$
(28)

where $w_{1}(S),\,w_{2}(S)$ are arbitrary complex functions and $\overline{S}$
is the complex conjugate of $S$. The 1st-order term in the expansion is then

$U^{(3/2)}=\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\Im{A^{(3/2)}S^{3/2}+B^{(3/2)}S^{1/2}+\frac{V}{4\alpha^{2}}f_{t}\overline{S}S^{-1/2}}\;.$
(29)

$A^{(3/2)}$ and $B^{(3/2)}$ are coefficients to be determined by matching to
the outer solution. The 2nd-order term in the expansion is

$\begin{split}&U^{(5/2)}=\\\
&\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\imaginary\\{\left(A^{(5/2)}S^{5/2}+B^{(5/2)}S^{3/2}+C^{(5/2)}S^{1/2}\right)+\left(\frac{mv}{2}f_{t}-vB^{(3/2)}_{t}+\frac{1}{2}f_{zz}-\frac{1}{2\alpha^{2}}f_{tt}\right)\frac{1}{2\alpha^{2}}\overline{S}S^{1/2}\\\
&+\left(\frac{V}{2}B^{(3/2)}f_{t}+\frac{1+V^{2}}{4\alpha^{2}}(f_{t})^{2}-\frac{1}{4}(f_{z})^{2}\right)\frac{1}{2\alpha^{2}}\overline{S}S^{-1/2}-\frac{V^{2}}{16\alpha^{4}}f_{tt}\overline{S}^{2}S^{-1/2}-\frac{V^{2}}{32\alpha^{4}}(f_{t})^{2}\overline{S}^{2}S^{-3/2}\Bigg{\\}}\;.\end{split}$
(30)

$A^{(5/2)}$, $B^{(5/2)}$ and $C^{(5/2)}$ are additional coefficients to be
determined in the matching.

To make the matching to the outer solution straightforward, we expand Eqs.
(27,29,30) around the unperturbed crack front position. Replacing
$S=\epsilon^{-1}(s-\epsilon f)$ where $s=x+i\alpha y$, and gathering terms of
the same order in $\epsilon$ we find the following expressions.

Terms of the order $\epsilon^{0}$:

$\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\Im{s^{1/2}+\frac{1}{3}ms^{3/2}+A^{(5/2)}s^{5/2}}\;.$
(31)

Terms of the order $\epsilon^{1}$:

$\begin{split}&\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\imaginary\\{-\frac{1}{2}fs^{-1/2}+(B^{(3/2)}-\frac{1}{2}mf)s^{1/2}+\frac{Vf_{t}}{4\alpha^{2}}\overline{s}s^{-1/2}+(B^{(5/2)}-\frac{5}{2}A^{(5/2)}f)s^{3/2}\\\
&-\frac{V^{2}f_{tt}}{16\alpha^{4}}\overline{s}^{2}s^{-1/2}+\frac{1}{4}\left(\frac{mV}{\alpha^{2}}f_{t}-\frac{2V}{\alpha^{2}}B^{(3/2)}_{t}+\frac{1}{\alpha^{2}}f_{zz}-\frac{1}{\alpha^{4}}f_{tt}\right)\overline{s}s^{1/2}\Bigg{\\}}\;.\end{split}$
(32)

Terms of the order $\epsilon^{2}$:

$\begin{split}&\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\imaginary\\{-\frac{f^{2}}{8}s^{-3/2}+\left(\frac{m}{8}f^{2}-\frac{fB^{(3/2)}}{2}-\frac{V}{4\alpha^{2}}ff_{t}\right)s^{-1/2}+\frac{Vff_{t}}{8\alpha^{2}}\overline{s}s^{-3/2}\\\
&\left(C^{(5/2)}-\frac{3}{2}B^{(5/2)}f+\frac{15}{8}A^{(5/2)}f^{2}+\frac{V}{2\alpha^{2}}B_{t}^{(3/2)}f-\frac{mV}{4\alpha^{2}}ff_{t}+\frac{ff_{tt}}{4\alpha^{4}}-\frac{1}{4\alpha^{2}}ff_{zz}\right)s^{1/2}\\\
&+\left(\frac{d}{dt}\left(-\frac{1}{16\alpha^{2}}mVf^{2}+\frac{1}{4\alpha^{2}}VB^{(3/2)}f+\frac{\left(1+V^{2}\right)}{16\alpha^{4}}\frac{d}{dt}f^{2}\right)-\frac{1}{16\alpha^{2}}\frac{d^{2}}{dz^{2}}f^{2}\right)\overline{s}s^{-1/2}\\\
&-\frac{V^{2}}{64\alpha^{4}}\frac{d^{2}}{dt^{2}}f^{2}\overline{s}^{2}s^{-3/2}\Bigg{\\}}\;.\end{split}$
(33)

The matching procedure is further simplified by introducing the polar
coordinates $s=re^{i\theta}$ and taking the imaginary part in the above
expressions:
Terms of the order $\epsilon^{0}$:

$\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\left(\sqrt{r}\sin{\frac{\theta}{2}}+A^{(3/2)}r^{3/2}\sin\frac{3\theta}{2}+A^{(5/2)}r^{5/2}\sin\frac{5\theta}{2}\right)\;.$
(34)

Terms of the order $\epsilon^{1}$:

$\begin{split}&\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\Bigg{\\{}\frac{f}{2}r^{-1/2}\sin{\frac{\theta}{2}}+\left(B^{(3/2)}-\frac{3}{2}A^{(3/2)}f\right)r^{1/2}\sin{\frac{\theta}{2}}-\frac{Vf_{t}}{4\alpha^{2}}r^{1/2}\sin\frac{3\theta}{2}\\\
&+\left(\frac{VB^{(3/2)}_{t}}{2\alpha^{2}}-\frac{3A^{(3/2)}Vf_{t}}{4\alpha^{2}}+\frac{f_{tt}}{4\alpha^{4}}-\frac{f_{zz}}{4\alpha^{2}}\right)r^{3/2}\sin{\frac{\theta}{2}}\\\
&+\left(B^{(5/2)}-\frac{5}{2}A^{(5/2)}f\right)r^{3/2}\sin\frac{3\theta}{2}+\frac{V^{2}f_{tt}}{16\alpha^{4}}r^{3/2}\sin\frac{5\theta}{2}\Bigg{\\}}\;.\end{split}$
(35)

Terms of the order $\epsilon^{2}$:

$\begin{split}&\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\Bigg{\\{}\frac{f^{2}}{8}r^{-3/2}\sin\frac{3\theta}{2}-\frac{Vff_{t}}{8\alpha^{2}}r^{-1/2}\sin\frac{5\theta}{2}+\left(\frac{1}{2}B^{(3/2)}f-\frac{3A^{(3/2)}f^{2}}{8}+\frac{Vff_{t}}{4\alpha^{2}}\right)r^{-1/2}\sin{\frac{\theta}{2}}\\\
&+\left(\frac{15}{8}f^{2}A^{(5/2)}-\frac{3}{2}fB^{(5/2)}+C^{(5/2)}+\frac{VfB_{t}^{(3/2)}}{2\alpha^{2}}-\frac{3A^{(3/2)}Vff_{t}}{4\alpha^{2}}+\frac{ff_{tt}}{4\alpha^{4}}-\frac{ff_{zz}}{4\alpha^{2}}\right)r^{1/2}\sin{\frac{\theta}{2}}\\\
&+\left(-\frac{VfB_{t}^{(3/2)}}{4\alpha^{2}}-\frac{Vf_{t}B^{(3/2)}}{4\alpha^{2}}+\frac{3A^{(3/2)}Vff_{t}}{8\alpha^{2}}-\frac{\left(1+V^{2}\right)\left(f_{t}^{2}+ff_{tt}\right)}{8\alpha^{4}}+\frac{ff_{zz}}{8\alpha^{2}}+\frac{f_{z}^{2}}{8\alpha^{2}}\right)r^{1/2}\sin\frac{3\theta}{2}\\\
&+\frac{V^{2}\left(f_{t}^{2}+ff_{tt}\right)}{32\alpha^{4}}r^{1/2}\sin\frac{7\theta}{2}\Bigg{\\}}\;.\end{split}$
(36)

###  3.2 The outer solution

We derive an asymptotic solution for the outer problem, where the fields are
expanded around the position of the unperturbed crack front. Let us consider
an expansion involving a single Fourier component $u=u^{(0)}+\epsilon
u^{(1)}+\epsilon^{2}u^{(2)}+...=u^{(0)}+\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\imaginary\\{(\epsilon
q^{(1)}(x,y)+\epsilon^{2}q^{(2)}(x,y))e^{i(kz+\omega t)}\\}$, where
superscripts mark the order of each term in $\epsilon$. Here

$u^{(0)}=\sqrt{\frac{2}{\pi}}\frac{K^{u}_{0}}{\mu}\imaginary\left\\{s^{1/2}+\frac{m}{3}s^{3/2}+\frac{n}{5}s^{5/2}\right\\}\;,$
(37)

where $s,m$ and $n$ are the same as in the previous sections. The expansion is
substituted in Eq. (21), which is then solved term-by-term. Once a solution is
found, the fields due to arbitrary perturbations can be computed by
superposing Fourier components. To affect the matching of the inner and outer
solutions, explicit expressions of $q^{(i)}$ will be equated with expressions
(34-36). Each $q^{(i)}(x,y)$ satisfies the equation

$\alpha^{2}q_{xx}+q_{yy}+(\omega^{2}-k^{2})q+2iV\omega q_{x}=0\;.$ (38)

A general solution to this equation is

$q=\frac{1}{2\pi}\int_{-\infty}^{\infty}\\!\mathrm{d}\xi\,\hat{q}(\xi)e^{i\omega\xi
x+\omega\gamma(\xi)y}\;,$ (39)

where $\gamma^{2}=\alpha^{2}\xi^{2}+2V\xi+(k/\omega)^{2}-1$. Without loss of
generality, we assume that $\omega>0$. To ensure the convergence of the
integral (39), we must take the branch of $\gamma$ which has a positive real
part. $\gamma$ has two branch points, $\xi=-\lambda_{+}$ and $\xi=\lambda_{-}$
where

$\lambda_{\pm}=\frac{1}{\alpha^{2}}\left(1-\alpha^{2}k^{2}/\omega^{2}\right)^{1/2}\pm\frac{V}{\alpha^{2}}\,.$
(40)

It is convenient, then, to take one branch cut extending from
$\xi=-\lambda_{+}$ to $-i\infty$ and the other from $\xi=\lambda_{-}$ to
$+i\infty$. As we shall see, the exact shape of the branch cuts does not
affect our calculations.

To determine $\hat{q}(\xi)$, we employ the boundary conditions Eq. (22). In
terms of $q$ these boundary conditions are

$\begin{split}q(x,0)&=0;\;x>0\,,\\\ q_{y}(x,0)&=0;\;x<0\,.\end{split}$ (41)

Fourier transforming the boundary conditions results in the relations

$\begin{split}\hat{q}(\xi)&=\int_{-\infty}^{0}\\!\mathrm{d}x\,q(x,0)e^{-i\omega\xi
x}\,,\\\
-\omega\gamma(\xi)\hat{q}(\xi)&=\int_{0}^{\infty}\\!\mathrm{d}x\,q_{y}(x,0)e^{-i\omega\xi
x}\,.\end{split}$ (42)

These relations provide us with information about the analytical domain of
$\hat{q}(\xi)$ and $\gamma(\xi)\hat{q}(\xi)$. The function $\hat{q}(\xi)$ is
analytic for $\Im{\xi}<0$ and the function $\gamma(\xi)\hat{q}(\xi)$ is
analytic for $\Im{\xi}>0$. We decompose $\gamma=\gamma^{+}\gamma^{-}$, where
$\gamma^{\pm}=\alpha^{1/2}(\xi\pm\lambda_{\pm})^{1/2}$ and $\gamma^{+}$
($\gamma^{-}$) is analytic for $\Im{\xi}>0$ ($\Im{\xi}<0$), according to our
choice of the branch of $\gamma$. Then, we readily see that
$W\equiv\gamma^{-}\hat{q}=\gamma\hat{q}/\gamma^{+}$ is an entire function
which can generally be represented as a power series in $\xi$. In fact,
$\gamma^{-}q^{(n)}(x,y)$ should be a polynomial of order $n-1$. This can be
understood from the previous section, where at each order the solution in the
inner variables becomes progressively more singular. As $\xi$ is the Fourier
conjugate of $x$, stronger singularities of $x$ translate into higher powers
of $\xi$.

Let us consider the 1st-order term $q^{(1)}(x,y)$, and set
$\hat{q}^{(1)}(\xi)=q_{0}/\gamma^{-}$. For convenience, we write the $x,y$
coordinates as $\Re{s},\Im{s}$ respectively, and expand the integral in Eq.
(39) in powers of $1/\xi$ (note that $\order{s\xi}\sim\order{1}$)

$\begin{split}&q^{(1)}(x,y)=\frac{1}{2\pi}\int\\!\mathrm{d}\xi\,\frac{q_{0}}{\alpha^{1/2}(\xi-\lambda_{-})^{1/2}}e^{i\omega\xi\Re{s}+\omega\gamma\Im{s}}=\\\
&\frac{q_{0}}{2\pi\alpha^{1/2}}\int\\!\mathrm{d}\xi\,e^{\omega|\xi|\Im{s}+i\xi\omega\Re{s}}\left(\frac{1}{\sqrt{\xi}}+\frac{\lambda_{-}-\omega|\xi|\Im{s}\left(\lambda_{-}-\lambda_{+}\right)}{2\xi^{3/2}}\right.\\\
&\left.+\frac{3\lambda_{-}^{2}+\xi^{2}\omega^{2}\Im{s}^{2}\left(\lambda_{-}-\lambda_{+}\right){}^{2}-\omega|\xi|\Im{s}\left(3\lambda_{-}^{2}+\lambda_{+}^{2}\right)}{8\xi^{5/2}}+\order{\xi^{-7/2}}\right)\;.\end{split}$
(43)

We integrate this expression using the identity

$\begin{split}\frac{1}{2\pi}\int\\!\mathrm{d}\xi\,&\frac{e^{\omega|\xi|\Im{s}+i\xi\omega\Re{s}}}{\xi^{n+1/2}}[1,\mathrm{sgn}\xi]=\\\
&-\frac{2^{2n}n!}{(2n)!}\frac{e^{i\pi/4}}{\sqrt{\pi}}(-i)^{n}|\omega
s|^{n-1/2}[\sin(n-\frac{1}{2})\theta,i\cos(n-\frac{1}{2})\theta]\;,\end{split}$
(44)

where $s=re^{i\theta}$ and $n\geq 0$. A similar identity holds for the
transformation $n\rightarrow-n$,

$\begin{split}\frac{1}{2\pi}\int\\!\mathrm{d}\xi\,&\frac{e^{\omega|\xi|\Im{s}+i\xi\omega\Re{s}}}{\xi^{-n+1/2}}[1,\mathrm{sgn}\xi]=\\\
&\frac{(2n)!}{2^{2n}n!}\frac{e^{i\pi/4}}{\sqrt{\pi}}i^{-n}|\omega
s|^{-n-1/2}[\sin(n-\frac{1}{2})\theta,-i\cos(n-\frac{1}{2})\theta]\,.\end{split}$
(45)

Integrating and substituting $s=re^{i\theta}$, the 1st-order solution is

$\begin{split}&q^{(1)}(x,y)=\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left[\frac{\sin\frac{\theta}{2}}{\sqrt{r}}+\sqrt{r}\left(-\frac{1}{4}i\omega\sin\frac{\theta}{2}\left(3\lambda_{-}+\lambda_{+}\right)+\frac{1}{4}i\omega\sin\frac{3\theta}{2}\left(\lambda_{-}-\lambda_{+}\right)\right)\right.\\\
&+r^{3/2}\left(-\frac{1}{16}\omega^{2}\sin\frac{\theta}{2}\left(5\lambda_{-}^{2}+2\lambda_{-}\lambda_{+}+\lambda_{+}^{2}\right)+\frac{1}{32}\omega^{2}\sin\frac{3\theta}{2}\left(5\lambda_{-}^{2}-2\lambda_{-}\lambda_{+}-3\lambda_{+}^{2}\right)\right.\\\
&\left.\left.-\frac{1}{32}\omega^{2}\sin\frac{5\theta}{2}\left(\lambda_{-}-\lambda_{+}\right){}^{2}\right)\right]\;.\end{split}$
(46)

Substituting the explicit expressions (40) for $\lambda_{\pm}$, the 1st-order
solution becomes

$\begin{split}&q^{(1)}(x,y)=\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left[r^{-1/2}\sin\frac{\theta}{2}+\left(\frac{iV\omega}{2\alpha^{2}}+\widehat{\Pi}\right)r^{1/2}\sin\frac{\theta}{2}-\frac{iV\omega}{2\alpha^{2}}r^{1/2}\sin\frac{3\theta}{2}\right.\\\
&\left.+\left(\frac{k^{2}}{2\alpha^{2}}-\frac{\omega^{2}}{2\alpha^{4}}-\frac{V^{2}\omega^{2}}{4\alpha^{4}}+\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{3/2}\sin\frac{\theta}{2}\right.\\\
&\left.+\left(\frac{V^{2}\omega^{2}}{8\alpha^{4}}-\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{3/2}\sin\frac{3\theta}{2}-\frac{V^{2}\omega^{2}}{8\alpha^{4}}r^{3/2}\sin\frac{5\theta}{2}\right]\,.\end{split}$
(47)

Here, we have introduced a new symbol $\widehat{\Pi}$ with the definition

$\begin{split}\widehat{\Pi}(k,\omega)&=-i\alpha^{-2}\mbox{sign}(\omega)\sqrt{\omega^{2}-\alpha^{2}k^{2}};\;\omega^{2}>\alpha^{2}k^{2}\,,\\\
\widehat{\Pi}(k,\omega)&=-\alpha^{-2}\sqrt{\alpha^{2}k^{2}-\omega^{2}};\;\omega^{2}<\alpha^{2}k^{2}\,.\end{split}$
(48)

The 2nd-order field $q^{(2)}(x,y)$ is similarly found by writing
$\hat{q}(\xi)=(q_{1}+q_{2}\xi)/\gamma^{-}$, so that

$q^{(2)}(x,y)=\frac{1}{2\pi}\int\\!\mathrm{d}\xi\,\frac{q_{1}+q_{2}\xi}{\alpha^{1/2}(\xi-\lambda_{-})^{1/2}}e^{i\omega\xi\Re{s}+\omega\gamma\Im{s}}\,.$
(49)

This integral can be separated into two parts, proportional to $q_{1}$ and
$q_{2}$ respectively. We have already calculated the first part in Eq. (47).
The order of the $r^{3/2}$ terms is higher than 2nd-order since the 2nd-order
inner solution in Eq. (36) does not contain them.

It remains therefore to calculate

$\begin{split}&\frac{1}{2\pi}\int\\!\mathrm{d}\xi\,\frac{q_{2}\xi}{\alpha^{1/2}(\xi-\lambda_{-})^{1/2}}e^{i\omega\xi\Re{s}+\omega\gamma\Im{s}}=\\\
&\frac{q_{2}}{2\pi\alpha^{1/2}}\int\\!\mathrm{d}\xi\,e^{\omega|\xi|\Im{s}+i\xi\omega\Re{s}}\left(\sqrt{\xi}+\frac{\lambda_{-}-\omega|\xi|\Im{s}\left(\lambda_{-}-\lambda_{+}\right)}{2\sqrt{\xi}}\right.\\\
&\left.+\frac{3\lambda_{-}^{2}+\omega^{2}|\xi|^{2}\Im{s}^{2}\left(\lambda_{-}-\lambda_{+}\right){}^{2}-\omega|\xi|\Im{s}\left(3\lambda_{-}^{2}+\lambda_{+}^{2}\right)}{8\xi^{3/2}}+\order{\xi^{-5/2}}\right)\,.\end{split}$
(50)

This is done in the same way as the 1st-order and the result in polar
coordinates is

$\begin{split}&-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\left[r^{-3/2}\sin\frac{3\theta}{2}+\left(\frac{3\,i\,V\;\omega}{2\alpha^{2}}+\widehat{\Pi}\right)r^{-1/2}\sin{\frac{\theta}{2}}-\frac{iV\omega}{2\alpha^{2}}r^{-1/2}\sin\frac{5\theta}{2}\right.\\\
&\left.+\left(\frac{k^{2}}{\alpha^{2}}-\frac{\omega^{2}}{\alpha^{4}}-\frac{7V^{2}\omega^{2}}{8\alpha^{4}}+\frac{5iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{1/2}\sin{\frac{\theta}{2}}+\left(-\frac{k^{2}}{2\alpha^{2}}+\frac{\omega^{2}}{2\alpha^{4}}+\frac{3V^{2}\omega^{2}}{4\alpha^{4}}-\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{1/2}\sin\frac{3\theta}{2}\right.\\\
&\left.-\frac{V^{2}\omega^{2}}{8\alpha^{4}}r^{1/2}\sin\frac{7\theta}{2}\right]\,.\end{split}$
(51)

The expression for the outer 2nd-order field becomes

$\begin{split}&q^{(2)}(x,y)=-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\left[r^{-3/2}\sin\frac{3\theta}{2}+\left(\frac{3\,i\,V\;\omega}{2\alpha^{2}}+\widehat{\Pi}\right)r^{-1/2}\sin{\frac{\theta}{2}}-\frac{iV\omega}{2\alpha^{2}}r^{-1/2}\sin\frac{5\theta}{2}\right.\\\
&\left.+\left(\frac{k^{2}}{\alpha^{2}}-\frac{\omega^{2}}{\alpha^{4}}-\frac{7V^{2}\omega^{2}}{8\alpha^{4}}+\frac{5iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{1/2}\sin{\frac{\theta}{2}}+\left(-\frac{k^{2}}{2\alpha^{2}}+\frac{\omega^{2}}{2\alpha^{4}}+\frac{3V^{2}\omega^{2}}{4\alpha^{4}}-\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)r^{1/2}\sin\frac{3\theta}{2}\right.\\\
&\left.-\frac{V^{2}\omega^{2}}{8\alpha^{4}}r^{1/2}\sin\frac{7\theta}{2}\right]\\\
&+\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left[r^{-1/2}\sin\frac{\theta}{2}+\left(\frac{iV\omega}{2\alpha^{2}}+\widehat{\Pi}\right)r^{1/2}\sin\frac{\theta}{2}-\frac{iV\omega}{2\alpha^{2}}r^{1/2}\sin\frac{3\theta}{2}\right]\,.\end{split}$
(52)

###  3.3 Matching the solutions

In the two previous sections, we have derived the asymptotic solutions up till
the 2nd-order in $\epsilon$ for the inner fields Eqs. (34-36) and for the
outer fields Eqs. (37,47,52). The expressions we have derived contain five
inner coefficients $A^{(3/2)},A^{(5/2)},B^{(3/2)},B^{(5/2)},C^{(5/2)}$ and
three outer coefficients $q_{0},q_{1},q_{2}$, which will be found by matching
the inner and outer solutions. We note that, in general, the coefficients may
be functions of $(z,t)$.

To match the terms in the two solutions we equate, order-by-order the
prefactors of the polar functions $r^{-j/2}\sin(j\theta/2)$, where $j$ is an
integer. The matching is made less straightforward by the fact that the inner
solution is expressed in the real space coordinates $(z,t)$ and the outer
solution is expressed in the Fourier coordinates $(k,\omega)$. However, we
will stick to the present notation, with the implicit understanding that the
outer solution terms must be first Fourier transformed before matching to the
corresponding inner solution terms. We also adopt the notation
$\Pi[f](z,t)=\frac{1}{(2\pi)^{2}}\int\\!\mathrm{d}\omega\mathrm{d}k\,e^{i(kz+\omega
t)}\widehat{\Pi}(\omega,k)f(\omega,k)$. The final expressions will always be
given in the $(z,t)$ space.

Term | Inner solution | Outer solution
---|---|---
$\sqrt{r}\sin{\frac{\theta}{2}}$ | $K^{u}_{0}$ | $K^{u}_{0}$
$r^{3/2}\sin\frac{3\theta}{2}$ | $A^{(3/2)}$ | $\frac{m}{3}$
$r^{5/2}\sin\frac{5\theta}{2}$ | $A^{(5/2)}$ | $\frac{n}{5}$
Table 1: Terms of order $\epsilon^{0}$ Term | Inner solution | Outer solution
---|---|---
$r^{-1/2}\sin{\frac{\theta}{2}}$ | $\frac{f}{2}$ | $\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}$
$r^{1/2}\sin{\frac{\theta}{2}}$ | $B^{(3/2)}-\frac{3}{2}A^{(3/2)}f$ | $\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left(\frac{iV\omega}{2\alpha^{2}}+\widehat{\Pi}\right)$
$r^{1/2}\sin\frac{3\theta}{2}$ | $-\frac{Vf_{t}}{4\alpha^{2}}$ | $-\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\frac{iV\omega}{2\alpha^{2}}$
$r^{3/2}\sin{\frac{\theta}{2}}$ | $B^{(3/2)}_{t}\frac{V}{2\alpha^{2}}-A^{(3/2)}\frac{3Vf_{t}}{4\alpha^{2}}+\frac{f_{tt}}{4\alpha^{4}}-\frac{f_{zz}}{4\alpha^{2}}$ | $\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left(\frac{k^{2}}{2\alpha^{2}}-\frac{\omega^{2}}{2\alpha^{4}}-\frac{V^{2}\omega^{2}}{4\alpha^{4}}+\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)$
$r^{3/2}\sin{\frac{3\theta}{2}}$ | $B^{(5/2)}-\frac{5}{2}A^{(5/2)}f$ | $\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left(\frac{V^{2}\omega^{2}}{8\alpha^{4}}-\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\right)$
$r^{3/2}\sin\frac{5\theta}{2}$ | $\frac{V^{2}f_{tt}}{16\alpha^{4}}$ | $-\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\frac{V^{2}\omega^{2}}{8\alpha^{4}}$
Table 2: Terms of order $\epsilon^{1}$ Term | Inner solution | Outer solution
---|---|---
$r^{-3/2}\sin\frac{3\theta}{2}$ | $\frac{f^{2}}{8}$ | $-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}$
$r^{-1/2}\sin{\frac{\theta}{2}}$ | $\frac{1}{2}B^{(3/2)}f-A^{(3/2)}\frac{3f^{2}}{8}+\frac{Vff_{t}}{4\alpha^{2}}$ | $-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\left(\frac{3\,i\,V\;\omega}{2\alpha^{2}}+\widehat{\Pi}\right)+\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}$
$r^{-1/2}\sin\frac{5\theta}{2}$ | $-\frac{Vff_{t}}{8\alpha^{2}}$ | $\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\frac{iV\omega}{2\alpha^{2}}$
$r^{1/2}\sin{\frac{\theta}{2}}$ |  $\frac{15}{8}f^{2}A^{(5/2)}-\frac{3}{2}fB^{(5/2)}+C^{(5/2)}\\\ +\frac{VfB_{t}^{(3/2)}}{2\alpha^{2}}-\frac{3A^{(3/2)}Vff_{t}}{4\alpha^{2}}+\frac{ff_{tt}}{4\alpha^{4}}-\frac{ff_{zz}}{4\alpha^{2}}$ |  $-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\Big{(}\frac{k^{2}}{\alpha^{2}}-\frac{\omega^{2}}{\alpha^{4}}-\frac{7V^{2}\omega^{2}}{8\alpha^{4}}\\\ +\frac{5iV\omega\widehat{\Pi}}{2\alpha^{2}}\Big{)}+\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left(\frac{iV\omega}{2\alpha^{2}}+\widehat{\Pi}\right)$
$r^{1/2}\sin\frac{3\theta}{2}$ |  $-B_{t}^{(3/2)}\frac{Vf}{4\alpha^{2}}-B^{(3/2)}\frac{Vf_{t}}{4\alpha^{2}}\\\ +A^{(3/2)}\frac{3Vff_{t}}{8\alpha^{2}}-\frac{\left(1+V^{2}\right)\left(f_{t}^{2}+ff_{tt}\right)}{8\alpha^{4}}\\\ +\frac{ff_{zz}}{8\alpha^{2}}+\frac{f_{z}^{2}}{8\alpha^{2}}$ |  $-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\Big{(}-\frac{k^{2}}{2\alpha^{2}}+\frac{\omega^{2}}{2\alpha^{4}}+\frac{3V^{2}\omega^{2}}{4\alpha^{4}}\\\ -\frac{iV\omega\widehat{\Pi}}{2\alpha^{2}}\Big{)}-\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}\left[\frac{iV\omega}{2\alpha^{2}}\right]$
$r^{1/2}\sin\frac{7\theta}{2}$ | $\frac{V^{2}\left(f_{t}^{2}+ff_{tt}\right)}{32\alpha^{4}}$ | $-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}\left[-\frac{V^{2}\omega^{2}}{8\alpha^{4}}\right]$
Table 3: Terms of order $\epsilon^{2}$

The zeroth order matching in Table 1 yields $A^{(3/2)}=\frac{m}{3}$ and
$A^{(5/2)}=\frac{n}{5}$. The 1st-order matching is given in Table 2. The first
row yields an expression for $q_{0}$;
$\frac{q_{0}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}=\frac{f}{2}$. The
second row translates into

$B^{(3/2)}=\frac{m}{2}f+\frac{V}{4\alpha^{2}}f_{t}+\frac{1}{2}\Pi[f]\,.$ (53)

Given the expressions for $A^{(3/2)},A^{(5/2)},q_{0}$ and $B^{(3/2)}$, the
third, fourth and sixth rows of Table 2 become identities. This is a
consistency check that our calculations so far have not contained a mistake.
The fifth row yields new information:

$B^{(5/2)}=\frac{n}{2}f-\frac{V^{2}}{16\alpha^{4}}f_{tt}-\frac{V}{4\alpha^{2}}\Pi[f_{t}]\,.$
(54)

The 2nd-order matching in Table 3 provides the three remaining coefficients
$q_{1},q_{2}$ and $C^{5/2}$. The first row of the table translates into
$-\frac{q_{2}\;ie^{i\pi/4}}{2\sqrt{\pi}\alpha^{1/2}\;(-\omega)^{3/2}}=\frac{f^{2}}{8}$.
The determination of $q_{2}$ makes the third and fifth rows of Table 3 an
identity. The second row gives $q_{1}$ through the equation

$\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}=\frac{1}{2}(\frac{m}{2}f+\frac{V}{4\alpha^{2}}f_{t}+\frac{1}{2}\Pi[f])f-\frac{mf^{2}}{8}+\frac{Vff_{t}}{4\alpha^{2}}-\frac{1}{8}\left(\frac{3Vff_{t}}{\alpha^{2}}+\Pi[f^{2}]\right)\,,$
(55)

or

$\frac{q_{1}e^{i\pi/4}}{\sqrt{-\pi\omega}\alpha^{1/2}}=\frac{mf^{2}}{8}+\frac{1}{4}\Pi[f]f-\frac{1}{8}\Pi[f^{2}]\,.$
(56)

The fourth row of Table 3 is the center of this calculation since it would
give the 2nd-order correction for the energy-release-rate $G$. Substituting
the expressions for $A^{(5/2)}$ and $B^{(5/2)}$ the inner solution cell of the
fourth row is

$-\frac{3}{8}nf^{2}+\frac{8+7V^{2}}{32\alpha^{4}}ff_{tt}+\frac{5V}{8\alpha^{2}}f\Pi[f_{t}]+C^{(5/2)}-\frac{ff_{zz}}{4\alpha^{2}}\,,$
(57)

while the outer solution cell is

$\begin{split}&-\frac{1}{4\alpha^{2}}(f_{z}^{2}+ff_{zz})+\frac{(8+7V^{2})}{32\alpha^{4}}(f_{t}^{2}+ff_{tt})+\frac{V\Pi[ff_{t}]}{2\alpha^{2}}+\frac{mVff_{t}}{8\alpha^{2}}\\\
&+\frac{V}{8\alpha^{2}}(f\Pi[f_{t}]+f_{t}\Pi[f])+\frac{m\Pi[f^{2}]}{8}+\frac{1}{4}\Pi[\Pi[f]f]-\frac{1}{8}\Pi[\Pi[f^{2}]]\,.\end{split}$
(58)

Since $\widehat{\Pi}^{2}=-\frac{1}{\alpha^{4}}(\omega^{2}-\alpha^{2}k^{2})$,
equating these two expressions yields a formula for the last remaining
coefficient $C^{(5/2)}$:

$\begin{split}&C^{(5/2)}=\frac{3}{8}nf^{2}+\frac{7V^{2}}{32\alpha^{4}}f_{t}^{2}-\frac{1}{4\alpha^{4}}ff_{tt}-\frac{V}{2\alpha^{2}}f\Pi[f_{t}]+\frac{1}{4\alpha^{2}}ff_{zz}+\frac{V\Pi[ff_{t}]}{2\alpha^{2}}\\\
&+\frac{mVff_{t}}{8\alpha^{2}}+\frac{V}{8\alpha^{2}}f_{t}\Pi[f]+\frac{m\Pi[f^{2}]}{8}+\frac{1}{4}\Pi[\Pi[f]f]\,.\end{split}$
(59)

###  3.4 Calculation of the energy-release-rate

The previous sections dealt with the derivation of the asymptotic field $u$ up
to the 2nd-order in the perturbation to the crack front. The expressions
derived in the last section for the undetermined coefficients of the inner
solution can now be inserted into Eqs. (34-36) to find the explicit dependence
of $u$ on $f$. To calculate the energy-release-rate $G$, it is simpler to use
Eqs. (27-30). Then, the displacement and stress intensity factors are given
respectively by

$\displaystyle K^{u}$ $\displaystyle=\lim\limits_{X\uparrow
f(z,t)}\epsilon^{-1/2}\left(X-f(z,t)\right)^{-1/2}\sqrt{\frac{\pi}{2}}\mu
U(X,0,z,t)\,,$ (60) $\displaystyle K^{\sigma}$
$\displaystyle=\sqrt{2\pi}\lim\limits_{X\downarrow
f(z,t)}\epsilon^{-1/2}\sqrt{2\pi\left(X-f(z,t)\right)}\mu U_{Y}(X,0,z,t)\,.$
(61)

Some algebra leads to the two 2nd-order expressions

$\displaystyle K^{u}=$ $\displaystyle
K^{u}_{0}\left(1+\epsilon\left(B^{(3/2)}+\frac{Vf_{t}}{4\alpha^{2}}\right)\right.$
$\displaystyle\left.+\epsilon^{2}\left(C^{(5/2)}+\frac{VB^{(3/2)}f_{t}}{4\alpha^{2}}+\frac{f_{t}^{2}}{8\alpha^{4}}+\frac{3V^{2}f_{t}^{2}}{32\alpha^{4}}-\frac{f_{z}^{2}}{8\alpha^{2}}\right)\right)\,,$
(62) $\displaystyle K^{\sigma}=$ $\displaystyle\sqrt{\frac{\pi}{2}}\alpha
K^{u}_{0}\left(1+\epsilon\left(B^{(3/2)}-\frac{3Vf_{t}}{4\alpha^{2}}\right)\right.$
$\displaystyle\left.+\epsilon^{2}\left(C^{(5/2)}-\frac{3VB^{(3/2)}f_{t}}{4\alpha^{2}}-\frac{3f_{t}^{2}}{8\alpha^{4}}-\frac{5V^{2}f_{t}^{2}}{32\alpha^{4}}+\frac{3f_{z}^{2}}{8\alpha^{2}}\right)\right)\,.$
(63)

The energy-release-rate is given by their product

$\begin{split}G\propto
K^{u}K^{\sigma}=\alpha(K^{u}_{0})^{2}\left(1+\epsilon\left(2B^{(3/2)}-\frac{Vf_{t}}{2\alpha^{2}}\right)\right.\\\
\left.+\epsilon^{2}\left((B^{(3/2)})^{2}+2C^{(5/2)}-\frac{VB^{(3/2)}f_{t}}{\alpha^{2}}-\frac{(1+V^{2})f_{t}^{2}}{4\alpha^{4}}+\frac{f_{z}^{2}}{4\alpha^{2}}\right)\right)\,.\end{split}$
(64)

Using Eqs. (53,59) derived in the last section for $B^{(3/2)}$ and
$C^{(5/2)}$, we write

$G=G_{r}g(V)(1+\epsilon\delta G^{(1)}+\epsilon^{2}\delta
G^{(2)}+\order{f^{3}})\,,$ (65)

where the 1st-order correction is

$\delta G^{(1)}=mf+\Pi[f]$ (66)

and the 2nd-order correction is

$\begin{split}\delta
G^{(2)}=&\frac{1}{4}m^{2}f^{2}+\frac{3}{4}nf^{2}+\frac{1}{2}mf\Pi[f]+\frac{1}{4}m\Pi\left[f^{2}\right]\\\
&+\frac{1}{4}\Pi[f]^{2}+\frac{1}{2}\Pi[f\Pi[f]]-\frac{Vf\Pi\left[f_{t}\right]}{\alpha^{2}}+\frac{V\Pi\left[ff_{t}\right]}{\alpha^{2}}-\frac{f_{t}^{2}}{4\alpha^{4}}-\frac{ff_{tt}}{2\alpha^{4}}+\frac{f_{z}^{2}}{4\alpha^{2}}+\frac{ff_{zz}}{2\alpha^{2}}\,\end{split}$
(67)

The 1st-order correction (66) was derived by Norris and Abrahams (2007), where
it was shown that the term $mf$ results in wave dispersion. This result is
possibly related to the observations of small amplitude wave dispersion
reported in Fineberg et al. (2003). Re-introducing dimensions, the coefficient
$m$ scales as $1/l$, and therefore the term $mf\sim\order{f/l}$ becomes
negligible when $f\ll l$. The term $\Pi[f]$, however, remains, as it does not
depend on sample geometry or loading conditions. In the following, we will
assume that $f/l\rightarrow 0$ and put $m=n=0$. The expansion of the energy-
release-rate (65) can be rewritten as a product of dynamical and “historical”
contributions, as in Eq. (1). To show that we follow Morrissey and Rice (2000)
and convert $\widehat{\Pi}(k,\omega)$ to the time domain using the identity

$\int\\!\mathrm{d}t\,e^{i\omega
t}\left(\partial_{t}+\frac{p\Theta(t)J_{1}(pt)}{t}\right)=-i\mbox{sign}(\omega)\sqrt{\omega^{2}-p^{2}}\,,$
(68)

where $\Theta$ is the Heaviside function and $J_{1}$ is the 1st-order Bessel
function. Then, $\Pi(k,t)=-\alpha^{-2}(\partial_{t}+\Psi(k,t))$, where
$\Psi(k,t)=\alpha k(\Theta(t)J_{1}(\alpha kt)/t)*$, with $*$ representing a
convolution in time. The appearance of the Heaviside function in $\Psi[f]$
signifies the dependence of real space functional $\Pi[f](z,t)$ on the past
history of the front. With this transformation Eq. (65) becomes (setting
$m=n=0$ and $\epsilon=1$)

$\displaystyle G=$ $\displaystyle
G_{r}g(V)\left\\{1-\frac{f_{t}}{\alpha^{2}}-\frac{1}{\alpha^{2}}\Psi[f]+\frac{1}{4\alpha^{4}}\Psi[f]^{2}+\frac{1}{2\alpha^{4}}\Psi[f\Psi[f]]+\right.$
$\displaystyle\left.\frac{1-2V}{2\alpha^{4}}\Psi[ff_{t}]+\frac{1+2V}{2\alpha^{4}}f\Psi[f_{t}]+\frac{1}{\alpha^{4}}f_{t}\Psi[f]+\frac{1-2V}{2\alpha^{4}}f_{t}^{2}+\frac{1}{4\alpha^{2}}f_{z}^{2}+\frac{1}{2\alpha^{2}}ff_{zz}\right\\}\,.$
(69)

To write this expression under the form given by Eq. (1), we develop
$V_{\perp}=V+f_{t}-\frac{V}{2}f_{z}^{2}+\order{f^{3}}$ so that

$g(V_{\perp})=g(V)\left(1-\frac{f_{t}}{\alpha^{2}}+\frac{V}{2\alpha^{2}}f_{z}^{2}+\frac{1-2V}{2\alpha^{4}}f_{t}^{2}+\order{f^{3}}\right)\,.$
(70)

Inserting this expression into Eq. (1) we have

$G=G_{r}g(V)\left(1-\frac{f_{t}}{\alpha^{2}}+H^{(1)}[f]+\frac{V}{2\alpha^{2}}f_{z}^{2}+\frac{1-2V}{2\alpha^{4}}f_{t}^{2}-\frac{f_{t}}{\alpha^{2}}H^{(1)}[f]+H^{(2)}[f]\right)\,.$
(71)

A direct comparison between Eqs. (3.4) and (71) then yields Eq. (3) for
$H^{(1)}[f]$ and

$H^{(2)}[f]=\frac{1}{4\alpha^{4}}\Psi[f]^{2}+\frac{1}{2\alpha^{4}}\Psi[f\Psi[f]]+\frac{1-2V}{2\alpha^{4}}\Psi[ff_{t}]+\frac{1+2V}{2\alpha^{4}}f\Psi[f_{t}]+\frac{1-2V}{4\alpha^{2}}f_{z}^{2}+\frac{1}{2\alpha^{2}}ff_{zz}\,.$
(72)

The expression for $H^{(2)}[f]$ can be further simplified through integration
by parts in the time domain. Using

$\begin{split}\Psi[f_{t}]&=\alpha k\int_{-\infty}^{t}\frac{J_{1}(\alpha
k(t-t^{\prime}))}{t-t^{\prime}}f_{t}(k,t^{\prime})dt^{\prime}\\\
&=\frac{\alpha^{2}k^{2}}{2}f(k,t)-\alpha^{2}k^{2}\int_{-\infty}^{t}\frac{J_{2}(\alpha
k(t-t^{\prime}))}{t-t^{\prime}}f(k,t^{\prime})dt^{\prime}\,,\end{split}$ (73)

the 2nd-order correction $H^{(2)}[f]$ is then given by Eq. (4).

##  4 Nonlinear perturbation of Mode I elastic crack fronts

In this section, we consider dynamic crack fronts in materials described by 3D
elastodynamics which involves three displacement components $u_{i}(x,y,z,t)$
and a system of three scalar wave equations coupled by boundary conditions.
The wave equations are a consequence of Newton’s 2nd law

$\frac{\partial\sigma_{ix}}{\partial x}+\frac{\partial\sigma_{iy}}{\partial
y}+\frac{\partial\sigma_{iz}}{\partial
z}=\rho\frac{\partial^{2}u_{i}}{\partial t^{2}}\,,$ (74)

where the 3D stress $\sigma_{ij}$ is linearly related to the strain tensor
through the Young modulus $E$ and the Poisson ratio $\nu$. Under pure Mode I
loading, the boundary conditions comprise tensile loads applied at the remote
boundaries and traction-free conditions on the crack faces, i.e.

$\sigma_{ij}n_{j}=0\,.$ (75)

The perturbation scheme of planar crack fronts follows the PhD thesis of
Ramanathan (1997). Parts of the calculation, that appear in the thesis, are
reproduced below for completeness. The explicit 2nd-order expressions and the
time-dependent formulation are the main contributions of this section. To
compute the expansion of $G$ in $f$, Ramanathan’s calculation aims at
obtaining an asymptotic solution for the elastic fields of a running in-plane
crack with a crack front defined by $h(z,t)=Vt+f(z,t)$ (Fig. 1). The solution
is to be obtained as a perturbation series around a straight front $f(z,t)=0$.
More specifically, the $x$-origin of the system of coordinates is locally
translated to the instantaneous position of the crack front $(X\equiv
x-h(z,t),y,z)$. Then the contribution of $f$ in the elastic field components
is “eliminated". This latter condition is fulfilled by introducing $f(z,t)$ at
the functional level through the transformation

$\begin{split}&u_{i}(X,y,z,t)=e^{f(z,t)\partial_{X}}U_{i}(X,y,z,t)\,,\\\
&\sigma_{ij}(X,y,z,t)=e^{f(z,t)\partial_{X}}\Sigma_{ij}(X,y,z,t)\,.\end{split}$
(76)

The advantage of this transformation is that the new fields $U_{i}$ and
$\Sigma_{ij}$ still follow the same constitutive linear elastic relations as
$u_{i}$ and $\sigma_{ij}$. Moreover, Eq. (74) becomes

$\frac{\partial\Sigma_{ix}}{\partial X}+\frac{\partial\Sigma_{iy}}{\partial
y}+\frac{\partial\Sigma_{iz}}{\partial
z}=\rho\left(V^{2}\frac{\partial^{2}U_{i}}{\partial
X^{2}}-2V\frac{\partial^{2}U_{i}}{\partial X\partial
t}+\frac{\partial^{2}U_{i}}{\partial t^{2}}\right)\,.$ (77)

Using the symmetry of Mode I loading, the boundary conditions for the fields
$U_{i}$ and $\Sigma_{ij}$ on the plane $y=0$ are given by

$\displaystyle\Sigma_{xy}(X,0,z,t)$ $\displaystyle=$
$\displaystyle\Sigma_{yz}(X,0,z,t)=0\,,$ (78) $\displaystyle U_{y}(X>0,0,z,t)$
$\displaystyle=$ $\displaystyle 0\,,\;\Sigma_{yy}(X<0,0,z,t)=0\,.$ (79)

The fields $\Sigma_{ij}$ and $U_{i}$ should be computed from the solution of
the elastodynamic equations in the reference frame $(X,y,z)$ with the
prescribed boundary conditions. This is done by using the general relations
obtained by Geubelle and Rice (1995) between the Fourier components of stress
and displacement fields in the plane of fracture. Specifically, the stress
component $\Sigma_{n}(X,z,t)\equiv\Sigma_{yy}(X,y=0,z,t)$ and the displacement
component $U_{n}(X,z,t)=U_{y}(X,y=0^{+},z,t)$ are related through

$\widehat{\Sigma}_{n}(q,k,\omega)=\widehat{\Upsilon}_{yy}(q,k,\Omega)\widehat{U}_{n}(q,k,\omega)\,,$
(80)

where the Fourier transforms are defined such that
$\widehat{F}(q,k,\omega)=\int\,dX\,dz\,dt\,e^{-iqX-ikz-i\omega t}F(X,z,t)$ and

$\hat{\Upsilon}_{yy}(q,k,\Omega)=-\mu\frac{p\,R\left(\frac{\Omega}{p}\right)}{\frac{\Omega^{2}}{b^{2}p^{2}}\sqrt{1-\frac{\Omega^{2}}{a^{2}p^{2}}}}\,,$
(81)

where $\mu$ is the shear modulus, $\Omega=\omega-qV$, $p=\sqrt{q^{2}+k^{2}}$
and $R(\zeta)$ is the Rayleigh function given by

$R(\zeta)=4\sqrt{1-\frac{\zeta^{2}}{a^{2}}}\sqrt{1-\frac{\zeta^{2}}{b^{2}}}-\left(2-\frac{\zeta^{2}}{b^{2}}\right)^{2}\,.$
(82)

Here $a$ and $b$ are the longitudinal and shear wave speeds.

It is left to explicitly compute $\Sigma_{yy}$ and $U_{y}$ by satisfying the
boundary condition (79). Notice that $\Sigma_{n}(X<0,z,t)=0$ (resp.
$U_{n}(X<0,z,t)=0$) is equivalent to $\widehat{\Sigma}_{n}(q,k,\omega)$ (resp;
$\widehat{U}_{n}(q,k,\omega)$) having its analytical domain encompassing the
upper (resp. lower) $q$-half-plane. Then, given a decomposition
$\widehat{\Upsilon}_{yy}=\widehat{\Upsilon}^{+}\widehat{\Upsilon}^{-}$ where
$\widehat{\Upsilon}^{+}$ is analytical for $\mbox{Im}(q)>0$ and
$\widehat{\Upsilon}^{-}$ is analytical for $\mbox{Im}(q)<0$, there exists an
analytical function $\widehat{W}(q,k,\omega)$ for all $q$’s such that
$\widehat{\Sigma}_{n}=\widehat{W}\widehat{\Upsilon}^{+}$ and
$\widehat{U}_{n}=\widehat{W}/\widehat{\Upsilon}^{-}$. Finally, the crux of the
method lies in demanding that the fields have the correct physical singularity
locally at the front such that

$\begin{split}&\sigma_{yy}(X,y=0,z,t)\rightarrow K^{\sigma}(z,t)(2\pi
X)^{-1/2}\;\;\text{ as }\;\;X\rightarrow 0^{+}\,,\\\
&u_{y}(X,y=0^{+},z,t)\rightarrow\frac{1-\nu^{2}}{E}K^{u}(z,t)(-2X/\pi)^{1/2}\;\;\text{
as }\;\;X\rightarrow 0^{-}\,,\end{split}$ (83)

where $K^{\sigma}$ and $K^{u}$ are the stress and displacement intensity
factors respectively. For a straight, unperturbed, crack front,
$K^{\sigma}=K^{\sigma}_{0}=K^{\sigma}_{r}k(V)$ and
$K^{u}=K^{u}_{0}=2A_{I}(V)K^{\sigma}_{0}$ where the rest SIF $K^{\sigma}_{r}$
is determined by the loading conditions (Freund, 1990) (see also section 2.2
for definitions of $k(V)$ and $A_{I}(V)$). For a curved crack front, the
demand, Eq. (83) will be met by eliminating all the higher order singularities
that appear on the right-hand side of Eq. (76). In $q$-space, this demand
imposes that the field $\widehat{u}_{y}$ shall not contain any power of $q$
higher than $q^{-3/2}$.

The energy-release-rate can be directly expressed using $K^{\sigma,u}$.
Following Freund (1990), we construct a small parallelopiped $P=(\delta
x,\delta y,\delta z)$ around the point on the crack front $x=h(z,t)$ that has
two sides of parallel to the $x-z$ plane, two sides parallel to the local
front tangent and $\hat{y}$, and two sides parallel to the $x-y$ plane. The
parallelopiped also travels at velocity $\dot{h}=V+f_{t}$ along the $x$-axis.
With these definitions, the local energy-release-rate is given by

$G(z,t)=\lim_{P\rightarrow 0}\left\\{\frac{1}{\delta
z}\frac{1}{\dot{h}}\int_{P}\\!\mathrm{d}S\,\sigma_{ij}\hat{n}_{j}\dot{u}_{i}\right\\}\,$
(84)

where $\hat{n}_{i}$ is the local normal to $P$ and $\mathrm{d}S$ is a surface
element. Since only $\sigma_{yy}$ and $\dot{u}_{y}$ have singular $X^{-1/2}$
behavior close to the crack front, the above integral simplifies into

$G(z,t)=\lim_{\delta x\rightarrow 0}\lim_{\delta y\rightarrow
0}\left\\{\frac{2}{\dot{h}}\int_{h(z,t)-\delta x}^{h(z,t)+\delta
x}\\!\mathrm{d}x\,\sigma_{yy}\dot{u}_{y}\right\\}=\lim_{\delta x\rightarrow
0}\lim_{\delta y\rightarrow 0}\left\\{2\int_{-\delta x}^{\delta
x}\\!\mathrm{d}X\,\Sigma_{yy}\frac{\partial U_{y}}{\partial X}\right\\}\,$
(85)

which leads to

$G(z,t)=\frac{1-\nu^{2}}{2E}K^{u}K^{\sigma}\,.$ (86)

Alternative formulations of the 3D path-independent energy flux integral give
the same result (Amestoy et al., 1981; Dodds et al., 1988; Eriksson, 2002;
Leguillon, 2014; El Kabir et al., 2018).

###  4.1 Implementation

To compute the energy-release-rate, we begin by applying the Fourier transform
$F(q)=\int\\!\mathrm{d}X\,e^{iqX}F(X)$ to Eqs. (76) ,

$\begin{split}&\widetilde{u}_{y}(q,y=0^{+},z,t)=e^{iqf(z,t)}\widetilde{U}_{n}(q,z,t)\,,\\\
&\widetilde{\sigma}_{yy}(q,y=0,z,t)=e^{iqf(z,t)}\widetilde{\Sigma}_{n}(q,z,t)\,.\end{split}$
(87)

Consider the displacement first. The zeroth order solution in the vicinity of
the straight crack front is given by

${u}_{0}(X,y=0^{+},z,t)=\frac{1-\nu^{2}}{E}K_{0}^{u}\sqrt{\frac{-2X}{\pi}}\Theta(-X)\,,$
(88)

which transforms into

$\widetilde{u}_{0}(q,y=0^{+},z,t)=\kappa^{u}q^{-3/2}\,,$ (89)

in the sense that $q$ has a vanishing positive imaginary part, and
$\kappa^{u}=i^{3/2}(1-\nu^{2})(\sqrt{2}E)^{-1}K_{0}^{u}$. The next orders will
be given by $\widehat{W}/\widehat{\Upsilon}^{-}$, such that

$\widetilde{u}_{y}(q,y=0^{+},z,t)=e^{iqf(z,t)}\left[\widetilde{u}_{0}(q,y=0^{+},z,t)+\int\frac{dk}{2\pi}\frac{d\omega}{2\pi}e^{ikz+i\omega
t}\frac{\widehat{W}(q,k,\omega)}{\widehat{\Upsilon}^{-}(q,k,\Omega)}\right]\,,$
(90)

where $\widehat{W}(q,k,\omega)$ is an analytical function of $q$. Now we
expand all functions in $q$ :

$\begin{split}&e^{iqf(z,t)}=\sum_{j=0}^{\infty}\frac{(iq)^{j}}{j!}f^{j}(z,t)\,,\\\
&\widetilde{W}(q,z,t)=\sum_{r=1}^{\infty}(iq)^{r-1}W_{r}(z,t)\,,\\\
&\frac{1}{\widetilde{\Upsilon}^{-}(q,z,t)}=\kappa^{u}q^{-1/2}\left(1+\sum_{n=1}^{\infty}\frac{\Lambda_{n}^{-}(z,t)}{(iq)^{n}}\right)\,.\end{split}$
(91)

Multiplying (90) by $q^{3/2}/\kappa^{u}$ we obtain

$\frac{q^{3/2}}{\kappa^{u}}\widetilde{u}_{y}(q,y=0^{+},z,t)=\sum_{j=0}^{\infty}(iq)^{j}\frac{f^{j}}{j!}\left[1-i\sum_{r=1}^{\infty}(iq)^{r}W_{r}-i\sum_{r=1}^{\infty}\sum_{n=1}^{\infty}(iq)^{r-n}\Lambda_{n}^{-}\star
W_{r}\ \right]\,,$ (92)

where $\star$ denotes convolution in $z$ and $t$. Eq. (92) can be simplified
into

$\begin{split}&\frac{q^{3/2}}{\kappa^{u}}\widetilde{u}_{y}(q,y=0^{+},z,t)=1-i\sum_{n=1}^{\infty}\sum_{r=1}^{n}\frac{f^{(n-r)}}{(n-r)!}(\Lambda_{n}^{-}\star
W_{r})\\\
&+\sum_{m=1}^{\infty}(iq)^{m}\left[\frac{f^{m}}{m!}-i\sum_{r=1}^{m}\frac{f^{(m-r)}}{(m-r)!}W_{r}-i\sum_{n=1}^{\infty}\sum_{r=1}^{m+n}\frac{f^{(m+n-r)}}{(m+n-r)!}(\Lambda_{n}^{-}\star
W_{r})\right]\,.\end{split}$ (93)

To meet the physical demand given by Eq. (83), we need to equate in Eq. (93)
the coefficients of $(iq)^{m}$ with $m\geq 1$ to zero. Therefore we get

$i\sum_{n=1}^{m}\frac{f^{(m-n)}}{(m-n)!}W_{n}+i\sum_{n=1}^{\infty}\sum_{r=1}^{m+n}\frac{f^{(m+n-r)}}{(m+n-r)!}(\Lambda_{n}^{-}\star
W_{r})=\frac{f^{m}}{m!}\,.$ (94)

Eq. (94) should be solved for $W_{m}$ order by order in $f$. Notice that Eq.
(94) shows that each $W_{r}$ has at least one term of the order $f^{r}$.
Therefore, $W_{r}$ can be written as

$W_{r}(z,t)=-i\sum_{j=r}^{\infty}w_{r,j}(z,t)\,,$ (95)

where $w_{r,j}$ is of the order $O(f^{j})$. Plugging this expansion in Eq.
(94) we find

$\begin{split}&\sum_{n=1}^{m}\frac{f^{(m-n)}}{(m-n)!}w_{n,n}=\frac{f^{m}}{m!}\,\quad
m\geq 1\,,\\\
&\sum_{n=1}^{m}\frac{f^{(m-n)}}{(m-n)!}w_{n,l+n-m}+\sum_{n=1}^{l-m}\sum_{r=1}^{m+n}\frac{f^{(m+n-r)}}{(m+n-r)!}(\Lambda_{n}^{-}\star
w_{r,l+r-m-n})=0\,\quad l>m\geq 1\,.\end{split}$ (96)

It is easy to show that this linear system of equations can be solved
recursively order by order. For the orders of interest, we find

$\begin{split}&w_{1,1}=f\,,\\\
&w_{2,2}=\frac{f^{2}}{2}-fw_{1,1}=-\frac{f^{2}}{2}\,,\\\
&w_{1,2}=-f(\Lambda_{1}^{-}\star w_{1,1})-\Lambda_{1}^{-}\star
w_{2,2}=-f(\Lambda_{1}^{-}\star
f)+\Lambda_{1}^{-}\star\frac{f^{2}}{2}\,.\end{split}$ (97)

It is left to express the SIF as a series in $f$. The expansion of the SIF
comprises the remaining non-zero terms in Eq. (93). Expressing these terms
with $w_{r,j}$,

$\frac{K^{u}}{K^{u}_{0}}=1-\sum_{p=1}^{\infty}{\cal
C}_{p}=1-\sum_{p=1}^{\infty}\sum_{n=1}^{p}\sum_{r=1}^{n}\frac{f^{(n-r)}}{(n-r)!}(\Lambda_{n}^{-}\star
w_{r,p+r-n})\,.$ (98)

For the orders of interest, we find

$\begin{split}&{\cal C}_{1}=\Lambda_{1}^{-}\star f\,,\\\ &{\cal
C}_{2}=-\Lambda_{1}^{-}\star\\{f(\Lambda_{1}^{-}\star
f)\\}+\Lambda_{1}^{-}\star\Lambda_{1}^{-}\star\frac{f^{2}}{2}+f\\{\Lambda_{2}^{-}\star
f\\}-\Lambda_{2}^{-}\star\frac{f^{2}}{2}\,.\end{split}$ (99)

Finally, the expression for the displacement intensity factor is

$\frac{K^{u}}{K^{u}_{0}}=1-\Lambda_{1}^{-}\star
f+\Lambda^{-}_{1}\star\\{f(\Lambda_{1}^{-}\star
f)\\}-\Lambda^{-}_{1}\star\Lambda_{1}^{-}\star\frac{f^{2}}{2}-f(\Lambda_{2}^{-}\star
f)+\Lambda_{2}^{-}\star\frac{f^{2}}{2}\,.$ (100)

A similar calculation shows that the SIF is

$\frac{K^{\sigma}}{K^{\sigma}_{0}}=1-\Lambda_{1}^{+}\star
f+\Lambda^{+}_{1}\star\\{f(\Lambda_{1}^{+}\star
f)\\}-\Lambda^{+}_{1}\star\Lambda_{1}^{+}\star\frac{f^{2}}{2}-f(\Lambda_{2}^{+}\star
f)+\Lambda_{2}^{+}\star\frac{f^{2}}{2}\,.$ (101)

where the functions $\Lambda_{n}^{+}$ are given by the expansion

$\widetilde{\Upsilon}^{+}(q,z,t)=(2i)^{-1/2}K_{0}^{\sigma}q^{1/2}\left(1+\sum_{n=1}^{\infty}\frac{\Lambda_{n}^{+}(z,t)}{(iq)^{n}}\right)\,.$
(102)

To obtain explicit expressions for the SIF and the energy-release-rate, one
should perform the decomposition of $\widehat{\Upsilon}_{yy}(q,k,\omega)$ and
its expansion in powers of $q$. This is done in the next section.

###  4.2 Decomposition of $\widehat{\Upsilon}_{yy}(q,k,\omega)$

According to our plan, we seek a Wiener-Hopf decomposition of
$\widehat{\Upsilon}_{yy}(q,k,\omega)$ (Eq. (81)) into a product of two
functions $\widehat{\Upsilon}^{+}(q,k,\omega)$ and
$\widehat{\Upsilon}^{-}(q,k,\omega)$ that are analytical for $\mbox{Im}\,q<0$
and $\mbox{Im}\,q>0$ respectively. Similarly to section 2.5 in Freund (1990),
a function $T(\Omega/p)\propto\widehat{\Upsilon}_{yy}$ is sought that
satisfies $T\rightarrow 1$ when $|q|\rightarrow\infty$ and $\log T=\log
T^{+}+\log T^{-}$, where $T^{+}$ ($T^{-}$) is analytical and non-zero in the
lower (upper) half of the $q$-plane. To find $T$, let us analyze the Rayleigh
function given by Eq. (82). It is known that $R(\zeta)$ has a double root at
$\zeta=0$ and two simple roots at $\zeta=\pm c$, where $c<b<a$ is the Rayleigh
velocity. To extract these roots from $\widehat{\Upsilon}_{yy}$, we define

$T(\zeta)=\frac{V^{2}}{R(V)\gamma_{c}^{2}}\frac{R(\zeta)}{\zeta^{2}\alpha_{c}(\zeta)^{2}}\,,$
(103)

with $\alpha_{c}(\zeta)=\sqrt{1-\frac{\zeta^{2}}{c^{2}}}$ and
$\gamma_{c}=1/\sqrt{1-V^{2}/c^{2}}$. Since $T(V)=1$, we have
$T(\Omega/p)\rightarrow 1$ as $\absolutevalue{q}\rightarrow\infty$. Then, Eq.
(81) is rewritten as

$\widehat{\Upsilon}_{yy}(q,k,\omega)=-\mu\frac{b^{2}R(V)\gamma_{c}^{2}}{V^{2}}p\frac{\alpha_{c}(\Omega/p)^{2}}{\alpha_{a}(\Omega/p)}T(\Omega/p)\,,$
(104)

with $\alpha_{a}(\zeta)=\sqrt{1-\frac{\zeta^{2}}{a^{2}}}$. Moreover, we notice
that

$p\frac{\alpha_{c}(\Omega/p)^{2}}{\alpha_{a}(\Omega/p)}=\frac{q^{2}+k^{2}-(\omega-
qV)^{2}/c^{2}}{\sqrt{q^{2}+k^{2}-(\omega-qV)^{2}/a^{2}}}\,.$ (105)

Both the denominator and the numerator in Eq. (105) contain similar quadratic
expressions in $q$ and they may be decomposed in the same manner. For example,
one has

$q^{2}+k^{2}-\frac{(\omega-
qV)^{2}}{c^{2}}=\frac{1}{\gamma_{c}^{2}}\left(q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}-iQ_{c}\right)\left(q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}+iQ_{c}\right)\,,$
(106)

where $Q_{c}=\gamma_{c}\sqrt{k^{2}-\frac{\omega^{2}\gamma_{c}^{2}}{c^{2}}}$.
Finally, one has

$\widehat{\Upsilon}_{yy}(q,k,\omega)=-\mu\frac{b^{2}R(V)}{V^{2}}\frac{\left(q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}-iQ_{c}\right)\left(q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}+iQ_{c}\right)}{\frac{1}{\gamma_{a}}\sqrt{q+\frac{V\omega\gamma_{a}^{2}}{a^{2}}-iQ_{a}}\sqrt{q+\frac{V\omega\gamma_{a}^{2}}{a^{2}}+iQ_{a}}}T(\Omega/p)\,.$
(107)

and

$\displaystyle\widehat{\Upsilon}^{+}$ $\displaystyle=$
$\displaystyle\frac{1}{(2i)^{1/2}}K_{0}^{\sigma}\frac{q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}-iQ_{c}}{\sqrt{q+\frac{V\omega\gamma_{a}^{2}}{a^{2}}-iQ_{a}}}T^{+}$
(108) $\displaystyle\widehat{\Upsilon}^{-}$ $\displaystyle=$
$\displaystyle\frac{2^{1/2}E}{i^{3/2}(1-\nu^{2})K_{0}^{u}}\frac{q+\frac{V\omega\gamma_{c}^{2}}{c^{2}}+iQ_{c}}{\sqrt{q+\frac{V\omega\gamma_{a}^{2}}{a^{2}}+iQ_{a}}}T^{-}\,.$
(109)

The decomposition of $\hat{\Upsilon}_{yy}$ reduces to resolving
$T=T^{+}T^{-}$. Let us identify the singularities of $T$ in the $q$ plane.
$R(\zeta)$ has two branch cuts: namely $b<\zeta<a$ and $-a<\zeta<-b$. Since
$\zeta=\frac{\omega-qV}{\sqrt{q^{2}+k^{2}}}$, in the $q$ plane these branch
cuts are transformed into into linear segments from $q_{a}$ to $q_{b}$ in the
upper $q$ plane and from $q_{a}^{*}$ to $q_{b}^{*}$ in the lower $q$ plane,
where $q_{a}=-\frac{\omega V\gamma_{a}^{2}}{a^{2}}+iQ_{a}$,
$q_{b}=-\frac{\omega V\gamma_{b}^{2}}{b^{2}}+iQ_{b}$ and
$Q_{a}=\gamma_{a}\sqrt{k^{2}-\frac{\omega^{2}\gamma_{a}^{2}}{a^{2}}}$,
$Q_{b}=\gamma_{b}\sqrt{k^{2}-\frac{\omega^{2}\gamma_{b}^{2}}{b^{2}}}$.

After Ramanathan (1997), we name the contour circling $[q_{a},q_{b}]$ in the
clockwise direction the branch cut $C_{+}$ and the contour circling
$[q_{a}^{*},q_{b}^{*}]$ in the clockwise direction the branch cut $C_{-}$.
Then, $\log T$ satisfies the conditions for the Wiener-Hopf decomposition: it
decays to $0$ as $\absolutevalue{q}\rightarrow\infty$ everywhere and it has
two finite branch cuts. Therefore, one has

$\displaystyle T^{+}$ $\displaystyle=$ $\displaystyle\exp\left(\frac{1}{2\pi
i}\oint_{C_{+}}\\!\frac{\log T}{\xi-q}\,\mathrm{d}\xi\right)$ (110)
$\displaystyle T^{-}$ $\displaystyle=$ $\displaystyle\exp\left(\frac{1}{2\pi
i}\oint_{C_{-}}\\!\frac{\log T}{\xi-q}\,\mathrm{d}\xi\right)\,.$ (111)

Let’s tackle the first integral. Writing it more explicitly we have

$\frac{1}{2\pi i}\oint_{C_{+}}\\!\frac{\log
T}{\xi-q}\,\mathrm{d}\xi=\frac{1}{2\pi
i}\oint_{C_{+}}\\!\frac{\mathrm{d}\xi}{\xi-q}\,\log\left(R(\zeta)/B(\zeta)\right)\,,$
(112)

where $\zeta=\frac{\omega-\xi V}{\sqrt{\xi^{2}+k^{2}}}$ and
$B(\zeta)=V^{-2}R(V)\gamma_{c}^{2}\zeta^{2}(1-\zeta^{2}/c^{2})$. Then, since
$B(\zeta)$ has no poles or zeros inside the contours $C_{\pm}$, the integral
over its logarithm is trivially zero and we have

$\log T^{\pm}=\frac{1}{2\pi i}\oint_{C_{\pm}}\\!\frac{\log
T}{\xi-q}\,\mathrm{d}\xi=\frac{1}{2\pi
i}\oint_{C_{\pm}}\\!\frac{\mathrm{d}\xi}{\xi-q}\,\log\left(R(\zeta)\right)\,.$
(113)

Changing the integration variable to $\zeta=\zeta(\xi)$ such that

$q_{\pm}(\zeta)=-\frac{\omega V}{\zeta^{2}-V^{2}}\pm
i\absolutevalue{k}\frac{\sqrt{\zeta^{2}}\sqrt{\zeta^{2}-H}}{\zeta^{2}-V^{2}}\,,$
(114)

where $H=V^{2}+\omega^{2}/k^{2}$. Here, the $+$ sign is taken for the $C_{+}$
integral and the $-$ sign is taken for the $C_{-}$ integral. Then

$\frac{\mathrm{d}\xi}{\xi-q}=\frac{q_{\pm}^{\prime}(\zeta)\,\mathrm{d}\zeta}{q_{\pm}(\zeta)-q}\,.$
(115)

This change of variables copies the contour $C_{+}$ in the $q$-plane to a
clockwise contour around the real interval $b<\zeta<a$ and the contour $C_{-}$
to a clockwise contour around $-a<\zeta<-b$. To evaluate the contour
integrals, a choice for the branch cuts of $R(\zeta)$ has to be made. To
ensure that the elastic fields maintain physical behavior far from the crack,
$\alpha_{b}$ must always have a non-negative real part (Geubelle and Rice,
1995; Ramanathan, 1997). For $\zeta=x+i\epsilon$ where $x>b$ and $|\epsilon|$
is arbitrarily small, we have $\alpha_{b}\simeq\pm
i\sqrt{x^{2}/b^{2}-1}(1+i\epsilon x/(x^{2}-b^{2}))$. Hence, when $\epsilon>0$,
we must choose the minus sign, and vice versa when $\epsilon<0$. Accordingly,
the integral which follows the clockwise contour $C_{+}$ is

$\displaystyle\log T^{+}=\frac{1}{2\pi
i}\oint_{C_{+}}\\!\frac{q_{+}^{\prime}(\zeta)\,\mathrm{d}\zeta}{q_{+}(\zeta)-q}\,\log\left(R(\zeta)\right)=$
$\displaystyle\frac{1}{2\pi
i}\int_{b}^{a}\\!\frac{q_{+}^{\prime}(\zeta)\,\mathrm{d}\zeta}{q_{+}(\zeta)-q}\log\left[-4i\sqrt{1-\frac{\zeta^{2}}{a^{2}}}\sqrt{\frac{\zeta^{2}}{b^{2}}-1}-(2-\frac{\zeta^{2}}{b^{2}})^{2}\right]+$
$\displaystyle\frac{1}{2\pi
i}\int_{a}^{b}\\!\frac{q_{+}^{\prime}(\zeta)\,\mathrm{d}\zeta}{q_{+}(\zeta)-q}\log\left[4i\sqrt{1-\frac{\zeta^{2}}{a^{2}}}\sqrt{\frac{\zeta^{2}}{b^{2}}-1}-(2-\frac{\zeta^{2}}{b^{2}})^{2}\right]\,.$
(116)

After a change of variable to $J=\zeta^{2}$ one has

$\log T^{+}=-\frac{1}{2\pi
i}\int_{b^{2}}^{a^{2}}\\!\frac{q_{+}^{\prime}(J)\,\mathrm{d}J}{q_{+}(J)-q}\\{\log(1-iz)-\log(1+iz)\\}\,,$
(117)

where
$z(J)=4\sqrt{1-\frac{J}{a^{2}}}\sqrt{\frac{J}{b^{2}}-1}/\left(2-\frac{J}{b^{2}}\right)^{2}$.
Using trigonometric identities, one can simplify $\log T^{+}$ further to
obtain

$\log
T^{+}=\frac{1}{\pi}\int_{b^{2}}^{a^{2}}\\!\frac{q_{+}^{\prime}(J)\,\mathrm{d}J}{q_{+}(J)-q}\arctan\left[4\frac{\sqrt{1-\frac{J}{a^{2}}}\sqrt{\frac{J}{b^{2}}-1}}{(2-\frac{J}{b^{2}})^{2}}\right]\,.$
(118)

A similar calculation shows that

$\log
T^{-}=-\frac{1}{\pi}\int_{b^{2}}^{a^{2}}\\!\frac{q_{-}^{\prime}(J)\,\mathrm{d}J}{q_{-}(J)-q}\arctan\left[4\frac{\sqrt{1-\frac{J}{a^{2}}}\sqrt{\frac{J}{b^{2}}-1}}{(2-\frac{J}{b^{2}})^{2}}\right]\
$ (119)

Now we are in a position to develop $T^{\pm}$ in powers of $1/(iq)$. Denoting

$\Pi^{\pm}_{n}=\mp\frac{1}{\pi}\int_{b^{2}}^{a^{2}}i^{n}q_{\pm}^{\prime}(J)q_{\pm}(J)^{n-1}\arctan[z(J)]\,\mathrm{d}J\,$
(120)

we find that,

$\log T^{\pm}=\sum_{n=1}^{\infty}\Pi^{\pm}_{n}\frac{1}{(iq)^{n}}\,.$ (121)

Therefore up to the 2nd-order,

$T^{\pm}=\exp\left(\sum_{n=1}^{\infty}\Pi^{\pm}_{n}\frac{1}{(iq)^{n}}\right)\simeq
1+\frac{\Pi_{1}^{\pm}}{iq}+\left(\Pi^{\pm}_{2}+\frac{1}{2}(\Pi^{\pm}_{1})^{2}\right)\frac{1}{(iq)^{2}}\,.$
(122)

Also, note that

$\left(\Pi^{+}_{n}\right)^{*}=(-1)^{n+1}\Pi^{-}_{n}\,,$ (123)

and that

$\Pi^{\pm}_{n}=-\frac{i^{n-1}}{2\pi}\oint_{C_{\pm}}\xi^{n-1}\log
T\\!\mathrm{d}\xi\,\,.$ (124)

We can now substitute the expansions of $T^{\pm}$ back in Eqs. (108). Equating

$\widehat{\Upsilon}^{+}=\frac{1}{(2i)^{1/2}}K_{0}^{\sigma}\frac{q-q_{c}}{\sqrt{q-q_{a}}}T^{+}=\frac{1}{(2i)^{1/2}}K_{0}^{\sigma}q^{1/2}\left(1+\sum_{n}\frac{\hat{\Lambda}_{n}^{+}}{(iq)^{n}}\right)\,$
(125)

where $q_{c}=-\frac{V\omega\gamma_{c}^{2}}{c^{2}}+iQ_{c}$. Since

$\frac{q-q_{c}}{\sqrt{q-q_{a}}}=q^{1/2}\left(1+\frac{iq_{a}/2-iq_{c}}{iq}+\frac{(1/2)q_{a}q_{c}-(3/8)q_{a}^{2}}{(iq)^{2}}+O(iq)^{3}\right)$
(126)

we obtain

$\displaystyle\hat{\Lambda}_{1}^{+}$ $\displaystyle=$ $\displaystyle
i\frac{q_{a}}{2}-iq_{c}+\Pi_{1}^{+}$ $\displaystyle\hat{\Lambda}_{2}^{+}$
$\displaystyle=$
$\displaystyle\frac{1}{2}q_{a}q_{c}-\frac{3}{8}q_{a}^{2}+\Pi_{2}^{+}+\frac{1}{2}(\Pi_{1}^{+})^{2}+\Pi_{1}^{+}\left(i\frac{q_{a}}{2}-iq_{c}\right)\,.$
(127)

The first of these equations correspond to Eq. (A.25) in Ramanathan (1997) and
the second corresponds to Eq. (A.35) in the limit $V=0$.

Now considering

$\frac{1}{\widehat{\Upsilon}^{-}}=\frac{i^{3/2}(1-\nu^{2})K_{0}^{u}}{2^{1/2}E}\frac{\sqrt{q-q_{a}^{*}}}{q-q_{c}^{*}}\frac{1}{T^{-}}=\frac{i^{3/2}(1-\nu^{2})K_{0}^{u}}{2^{1/2}E}q^{-1/2}\left(1+\sum_{n}\frac{\hat{\Lambda}_{n}^{-}}{(iq)^{n}}\right)\,$
(128)

we develop in the same way

$\frac{\sqrt{q-q_{a}^{*}}}{q-q_{c}^{*}}=q^{-1/2}\left(1+\frac{iq_{c}^{*}-iq_{a}^{*}/2}{iq}+\frac{q_{a}^{*}q_{c}^{*}/2+(q_{a}^{*})^{2}/8-(q_{c}^{*})^{2}}{(iq)^{2}}\right)$
(129)

and therefore

$\displaystyle\hat{\Lambda}_{1}^{-}$ $\displaystyle=$ $\displaystyle
iq_{c}^{*}-i\frac{q_{a}^{*}}{2}+\Pi_{1}^{-}$
$\displaystyle\hat{\Lambda}_{2}^{-}$ $\displaystyle=$
$\displaystyle\frac{1}{2}q_{a}^{*}q_{c}^{*}+\frac{1}{8}(q_{a}^{*})^{2}-(q_{c}^{*})^{2}+\Pi_{2}^{-}+\frac{1}{2}(\Pi_{1}^{-})^{2}+\frac{1}{2}\Pi_{1}^{-}\left(iq_{c}^{*}-i\frac{q_{a}^{*}}{2}\right)\,.$
(130)

###  4.3 Evaluation of the energy-release-rate

To compute $G$ we can now use Eq. (86) and substitute Eqs. (100,101). Since
explicit expressions for $\Lambda_{n}^{\pm}$ were obtained in $(k,\omega)$
space, we will compute the Fourier transform $\widehat{G}(k,\omega)$. First,
note that from Eqs. (100,101)

$\begin{split}\frac{\widehat{K}^{\sigma,u}(k,\omega)}{K^{\sigma,u}_{0}}=(2\pi)^{2}\delta(k)\delta(\omega)-\hat{\Lambda_{1}}^{\pm}\hat{f}+\hat{\Lambda}_{1}^{\pm}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{\pm}\hat{f}\right)\\}\\\
-\frac{1}{2}(\hat{\Lambda}_{1}^{\pm})^{2}(\hat{f}\otimes\hat{f})-\hat{f}\otimes\left(\hat{\Lambda}_{2}^{\pm}\hat{f}\right)+\frac{1}{2}\hat{\Lambda}_{2}^{\pm}(\hat{f}\otimes\hat{f})\,,\end{split}$
(131)

where

$\widehat{(f\cdot
g)}(k,\omega)=\widehat{f}\otimes\widehat{g}=\int\frac{d\omega^{\prime}}{2\pi}\frac{dk^{\prime}}{2\pi}\hat{f}(k-k^{\prime},\omega-\omega^{\prime})\hat{g}(k^{\prime},\omega^{\prime})\,.$
(132)

Then, the Fourier transformed Eq. (86) becomes

$\begin{split}\frac{\widehat{G}(k,\omega)}{G_{r}g(V)}=(2\pi)^{2}\delta(k)\delta(\omega)-(\hat{\Lambda}_{1}^{+}+\hat{\Lambda}_{1}^{-})\hat{f}+\hat{\Lambda}_{1}^{+}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{+}\hat{f}\right)\\}+\hat{\Lambda}_{1}^{-}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{-}\hat{f}\right)\\}\\\
+\frac{1}{2}\left(\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-}-(\hat{\Lambda}_{1}^{+})^{2}-(\hat{\Lambda}_{1}^{-})^{2}\right)(\hat{f}\otimes\hat{f})-\hat{f}\otimes\\{(\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-})\hat{f}\\}+(\hat{\Lambda}_{1}^{+}\hat{f})\otimes(\hat{\Lambda}_{1}^{-}\hat{f})\,.\end{split}$
(133)

Let us write explicit expressions for each of these terms. First, using the
expressions of $\hat{\Lambda}_{1}^{+}$ and $\hat{\Lambda}_{1}^{-}$, we obtain

$\widehat{\delta
G^{(1)}}(k,\omega)=\hat{\Lambda}_{1}^{+}+\hat{\Lambda}_{1}^{-}=2|k|P_{1}\left(\frac{\omega}{|k|};V,\nu\right)\,,$
(134)

where

$\begin{split}P_{1}(s;V,\nu)&=-\frac{1}{2}\gamma_{a}\sqrt{1-\frac{\gamma_{a}^{2}s^{2}}{a^{2}}}+\gamma_{c}\sqrt{1-\frac{\gamma_{c}^{2}s^{2}}{c^{2}}}\\\
&+\frac{1}{2\pi}\int_{b^{2}}^{a^{2}}\frac{(s^{2}+V^{2})(J+V^{2})-2JV^{2}}{\sqrt{J}\sqrt{J-(s^{2}+V^{2})}\left(J-V^{2}\right)^{2}}\arctan[z(J)]\,\mathrm{d}J\,.\end{split}$
(135)

The branch cuts of the square roots are defined so that
$\sqrt{1-s^{2}}=i\,\mathrm{sign}(s)\sqrt{s^{2}-1}$ for $|s|>1$. This
definition is consistent with the physical demand that elastic waves are
radiated away from the crack front (Geubelle and Rice, 1995; Ramanathan,
1997).

Second, separating real and imaginary parts of the term
$\hat{\Lambda}_{1}^{+}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{+}\hat{f}\right)\\}+\hat{\Lambda}_{1}^{-}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{-}\hat{f}\right)\\}$
yields

$\begin{split}&\hat{\Lambda}_{1}^{+}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{+}\hat{f}\right)\\}+\hat{\Lambda}_{1}^{-}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{-}\hat{f}\right)\\}=\\\
&2\mathrm{Re}(\hat{\Lambda}_{1}^{+})\\{\hat{f}\otimes\left(\mathrm{Re}(\hat{\Lambda}_{1})^{+}\hat{f}\right)\\}-2\mathrm{Im}(\hat{\Lambda}_{1}^{+})\\{\hat{f}\otimes\left(Im(\hat{\Lambda}_{1}^{+})\hat{f}\right)\\}\,.\end{split}$
(136)

Notice that the imaginary part of $\hat{\Lambda}_{1}$ is linear in $\omega$

$\displaystyle\mathrm{Im}(\hat{\Lambda}^{+}_{1})=VQ_{1}(V)\omega\,,$ (137)

with

$Q_{1}(V)=\frac{1}{2(a^{2}-V^{2})}-\frac{1}{c^{2}-V^{2}}+\frac{1}{\pi}\int_{b^{2}}^{a^{2}}\frac{1}{\left(J-V^{2}\right)^{2}}\arctan[z(J)]\,\mathrm{d}J\,,$
(138)

or (Adda-Bedia et al., 2013; Freund, 1990)

$Q_{1}(V)=-\frac{1}{2V}\frac{\mathrm{d}}{\mathrm{d}V}\left(\log\left(\frac{V^{2}}{R(V)\gamma_{a}b^{2}}\right)\right)\,,$
(139)

where
$R(V)=4\sqrt{1-\frac{V^{2}}{a^{2}}}\sqrt{1-\frac{V^{2}}{b^{2}}}-\left(2-\frac{V^{2}}{b^{2}}\right)^{2}$.
Hence,

$\hat{\Lambda}_{1}^{+}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{+}\hat{f}\right)\\}+\hat{\Lambda}_{1}^{-}\\{\hat{f}\otimes\left(\hat{\Lambda}_{1}^{-}\hat{f}\right)\\}=2|k|P_{1}\\{\hat{f}\otimes\left(|k|P_{1}\hat{f}\right)\\}-2V^{2}Q_{1}(V)^{2}\omega\\{\hat{f}\otimes\left(\omega\hat{f}\right)\\}\,.$
(140)

Third, to simplify the terms
$\frac{1}{2}\left(\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-}-(\hat{\Lambda}_{1}^{+})^{2}-(\hat{\Lambda}_{1}^{-})^{2}\right)(\hat{f}\otimes\hat{f})$
and $f\otimes\\{(\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-})\hat{f}\\}$, we
notice that

$\hat{\Lambda}_{2}^{+}=\frac{1}{2}(\hat{\Lambda}^{+}_{1})^{2}+\Pi_{2}^{+}+\frac{1}{2}\left(q_{c}^{2}-\frac{1}{2}q_{a}^{2}\right),\,\hat{\Lambda}_{2}^{-}=\frac{1}{2}(\hat{\Lambda}^{-}_{1})^{2}+\Pi_{2}^{-}-\frac{1}{2}\left(q_{c}^{2}-\frac{1}{2}q_{a}^{2}\right)^{*}$
(141)

Since $\hat{\Lambda}^{-}_{1}=(\hat{\Lambda}^{+}_{1})^{*}$ and
$(\Pi^{+}_{2})^{*}=-\Pi^{-}_{2}$, we find

$\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-}=\real[(\hat{\Lambda}_{1}^{+})^{2}]+2i\imaginary[\Pi_{2}^{+}]+i\imaginary[q_{c}^{2}-\frac{1}{2}q_{a}^{2}]\,.$
(142)

Using the identity
$\real[(\hat{\Lambda}_{1}^{+})^{2}]=\real[\hat{\Lambda}_{1}^{+}]^{2}-\imaginary[\hat{\Lambda}_{1}^{+}]^{2}=k^{2}P_{1}^{2}-V^{2}Q_{1}^{2}\omega^{2}$
and defining
$iV\omega|k|P_{2}(\omega/|k|;V,\nu)=-2i\imaginary[\Pi_{2}^{+}]-i\imaginary[q_{c}^{2}-\frac{1}{2}q_{a}^{2}]$
we write

$\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-}=k^{2}P_{1}^{2}-V^{2}Q_{1}^{2}\omega^{2}-iV\omega|k|P_{2}\,.$
(143)

Thus, we can write

$\frac{1}{2}\left(\hat{\Lambda}_{2}^{+}+\hat{\Lambda}_{2}^{-}-(\hat{\Lambda}_{1}^{+})^{2}-(\hat{\Lambda}_{1}^{-})^{2}\right)=-\frac{1}{2}\left(k^{2}P_{1}^{2}-V^{2}Q_{1}^{2}\omega^{2}\right)-\frac{i}{2}V\omega|k|P_{2}\,,$
(144)

where

$\begin{split}P_{2}(s;V,\nu)=&2\frac{\gamma_{c}^{3}}{c^{2}}\sqrt{1-\frac{\gamma_{c}^{2}s^{2}}{c^{2}}}-\frac{\gamma_{a}^{3}}{a^{2}}\sqrt{1-\frac{\gamma_{a}^{2}s^{2}}{a^{2}}}\\\
&+\frac{2}{\pi}\int_{b^{2}}^{a^{2}}\frac{(s^{2}+V^{2})(3J+V^{2})-2J(J+V^{2})}{\sqrt{J-(s^{2}+V^{2})}\left(J-V^{2}\right)^{3}}\arctan[z(J)]\,\frac{\mathrm{d}J}{2\sqrt{J}}\,.\end{split}$
(145)

Finally, using the distributive law, one obtains

$(\hat{\Lambda}_{1}^{+}\hat{f})\otimes(\hat{\Lambda}_{1}^{-}\hat{f})=(|k|P_{1}\hat{f})\otimes(|k|P_{1}\hat{f})+V^{2}Q_{1}^{2}(\omega\hat{f})\otimes(\omega\hat{f})\,.$
(146)

Collecting all the above expressions, the 2nd-order correction to the energy-
release-rate can be written as

$\begin{split}&\widehat{\delta G^{(2)}}(k,\omega)=\\\
&2|k|P_{1}\\{\hat{f}\otimes\left(|k|P_{1}\hat{f}\right)\\}-2Q_{1}^{2}V^{2}\omega\\{\hat{f}\otimes\left(\omega\hat{f}\right)\\}-\left\\{\frac{1}{2}\left(k^{2}P_{1}^{2}-V^{2}Q_{1}^{2}\omega^{2}\right)+\frac{i}{2}V\omega|k|P_{2}\right\\}(\hat{f}\otimes\hat{f})\\\
&-\hat{f}\otimes\left\\{\left(k^{2}P_{1}^{2}-V^{2}Q_{1}^{2}\omega^{2}-iV\omega|k|P_{2}\right)\hat{f}\right\\}+(|k|P_{1}\hat{f})\otimes(|k|P_{1}\hat{f})+V^{2}Q_{1}^{2}(\omega\hat{f})\otimes(\omega\hat{f})\,.\end{split}$
(147)

This expression can be further simplified using the identity
$\hat{f}\otimes(\omega\hat{f})=\frac{1}{2}\omega(\hat{f}\otimes\hat{f})$.
Then,

$\begin{split}\widehat{\delta
G^{(2)}}(k,\omega)=&2|k|P_{1}\\{\hat{f}\otimes\left(|k|P_{1}\hat{f}\right)\\}-\left\\{\frac{1}{2}k^{2}P_{1}^{2}+\frac{i}{2}V\omega|k|P_{2}\right\\}(\hat{f}\otimes\hat{f})\\\
&-\hat{f}\otimes\left\\{\left(k^{2}P_{1}^{2}-iV\omega|k|P_{2}\right)\hat{f}\right\\}+(|k|P_{1}\hat{f})\otimes(|k|P_{1}\hat{f})\,.\end{split}$
(148)

###  4.4 Transformation to $(k,t)$ space

To obtain a real space expression of the form
$G=g(V_{\perp})(1+H^{(1)}[f]+H^{(2)}[f,f]+\order{f^{3}})$ where the
functionals $H^{(i)}$ only depend on the history of the front
dynamics$\left\\{f(z,t^{\prime});\;t^{\prime}<t\right\\}$, we transform
$\widehat{\delta G^{(1)}}(k,\omega)$ and $\widehat{\delta G^{(2)}}(k,\omega)$
from $\omega$ to $t$. Following Morrissey and Rice (2000), we use the Bessel
identities

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}\sqrt{\beta^{2}-\omega^{2}}=\delta^{\prime}(t)+p^{2}\frac{J_{1}(pt)}{pt}\mathcal{H}(t)$

and

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}\frac{1}{\sqrt{\beta^{2}-\omega^{2}}}=J_{0}(pt)\mathcal{H}(t)$

where $\mathcal{H}(t)$ is the Heaviside function, to write

$\begin{split}&\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}|k|P_{1}\left(\frac{\omega}{|k|}\right)=\left(-\frac{1}{2}\frac{a}{a^{2}-V^{2}}+\frac{c}{c^{2}-V^{2}}\right)\delta^{\prime}(t)\\\
&-\partial_{t}^{2}\frac{1}{2}\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{2}}J_{0}(\beta_{\eta}t)\mathcal{H}(t)\\\
&+k^{2}\left(-\frac{a}{2}\frac{J_{1}(\beta_{a}t)}{\beta_{a}t}+c\frac{J_{1}(\beta_{c}t)}{\beta_{c}t}-\frac{1}{2}\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{V^{2}}{\eta^{2}-V^{2}}J_{0}(\beta_{\eta}t)\right)\mathcal{H}(t)\end{split}$
(149)

where $\beta_{s}=\sqrt{s^{2}-V^{2}}|k|$. From the Bessel identities
$J_{0}^{\prime}(x)=-J_{1}(x)$ and
$J_{0}^{\prime\prime}(x)=(J_{2}(x)-J_{0}(x))/2$

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}|k|P_{1}\left(\frac{\omega}{|k|}\right)=C_{1}\delta^{\prime}(t)+B_{1}(k,t)\delta(t)+A_{1}(k,t)\mathcal{H}(t)$
(150)

$\displaystyle C_{1}(k,t)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\frac{a}{a^{2}-V^{2}}+\frac{c}{c^{2}-V^{2}}-\frac{1}{2}\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{2}}\
J_{0}(\beta_{\eta}t)$ $\displaystyle B_{1}(k,t)$ $\displaystyle=$
$\displaystyle-|k|\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{3/2}}\
J_{1}(\beta_{\eta}t)$ $\displaystyle A_{1}(k,t)$ $\displaystyle=$
$\displaystyle
k^{2}\left[-\frac{a}{2}\frac{J_{1}(\beta_{a}t)}{\beta_{a}t}+c\frac{J_{1}(\beta_{c}t)}{\beta_{c}t}-\frac{1}{4}\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\left(\frac{\eta^{2}+V^{2}}{\eta^{2}-V^{2}}J_{2}(\beta_{\eta}t)-J_{0}(\beta_{\eta}t)\right)\right]\,.$

Then,

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}|k|P_{1}\left(\frac{\omega}{|k|}\right)\hat{f}(k,\omega)=-\pi_{1}\bar{f}_{t}-I_{1}[\bar{f}]$
(151)

where $\pi_{1}$ is given by Eq. (18), and
$I_{1}[\bar{f}]=-\int_{-\infty}^{t}\\!\mathrm{d}t^{\prime}\,A_{1}(k,t-t^{\prime})\bar{f}(k,t^{\prime})$.

Since $2\pi_{1}=g^{\prime}(V)/g(V)$, to the 1st-order in $f$ we have
$G=g(V_{\perp})(1+H^{(1)}[f]+\order{f^{2}})$ where
$\overline{H^{(1)}[f]}=2I_{1}[\bar{f}]$ Morrissey and Rice (2000).

To extend this result to the 2nd-order, we use the Bessel identity
$J_{0}^{\prime\prime\prime}(x)=3J_{1}(x)/4-J_{3}(x)/4$ to find

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}i\omega|k|P_{2}\left(\frac{\omega}{|k|}\right)=D_{2}\delta^{\prime\prime}(t)+C_{2}\delta^{\prime}(t)+B_{2}\delta(t)+A_{2}\mathcal{H}(t)\,,$
(152)

where

$\displaystyle D_{2}(k,t)$ $\displaystyle=$ $\displaystyle
2\frac{c}{(c^{2}-V^{2})^{2}}-\frac{a}{(a^{2}-V^{2})^{2}}-\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{3\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{3}}J_{0}(\beta_{\eta}t)$
$\displaystyle C_{2}(k,t)$ $\displaystyle=$ $\displaystyle
3|k|\int_{b}^{a}\\!\mathrm{d}\eta\,\Theta(\eta)\frac{3\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{5/2}}J_{1}(\beta_{\eta}t)$
$\displaystyle B_{2}(k,t)$ $\displaystyle=$ $\displaystyle
k^{2}\Bigg{[}2\frac{c}{c^{2}-V^{2}}\frac{J_{1}(\beta_{c}t)}{\beta_{c}t}-\frac{a}{a^{2}-V^{2}}\frac{J_{1}(\beta_{a}t)}{\beta_{a}t}$
$\displaystyle+$
$\displaystyle\frac{1}{2}\int_{b}^{a}\\!\mathrm{d}\eta\,\frac{\Theta(\eta)}{(\eta^{2}-V^{2})^{2}}\left((5\eta^{2}+V^{2})J_{0}(\beta_{\eta}t)-(9\eta^{2}+3V^{2})J_{2}(\beta_{\eta}t)\right)\Bigg{]}$
$\displaystyle A_{2}(k,t)$ $\displaystyle=$
$\displaystyle|k|^{3}\Bigg{[}-\gamma_{a}\frac{J_{2}(\beta_{a}t)}{\beta_{a}t}+2\gamma_{c}\frac{J_{2}(\beta_{c}t)}{\beta_{c}t}$
$\displaystyle-$
$\displaystyle\frac{1}{4}\int_{b}^{a}\\!\mathrm{d}\eta\,\frac{\Theta(\eta)}{\sqrt{\eta^{2}-V^{2}}}\left(\frac{3\eta^{2}+V^{2}}{\eta^{2}-V^{2}}J_{3}(\beta_{\eta}t)-J_{1}(\beta_{\eta}t)\right)\Bigg{]}\,.$

Then,

$\frac{1}{2\pi}\int\\!\mathrm{d}\omega\,e^{i\omega
t}i\omega|k|P_{2}\left(\frac{\omega}{|k|}\right)\hat{f}(k,\omega)=-\pi_{2}\bar{f}_{tt}-\pi_{1}k^{2}\bar{f}-I_{2}[\bar{f}]$
(153)

where

$\pi_{2}=\frac{a}{(a^{2}-V^{2})^{2}}-2\frac{c}{(c^{2}-V^{2})^{2}}+\int_{b}^{a}\frac{3\eta^{2}+V^{2}}{(\eta^{2}-V^{2})^{3}}\Theta(\eta)\mathrm{d}\eta\,.$
(154)

The second term results from the identity
$\partial_{t}^{2}D_{2}(0)-\partial_{t}C_{2}(0)+B_{2}(0)=-\pi_{1}k^{2}$, and
the third term is
$I_{2}[\bar{f}]=-\int_{-\infty}^{t}\\!\mathrm{d}t\,A_{2}(k,t-t^{\prime})\bar{f}(k,t^{\prime})$.

The 2nd-order contribution to $G$ is then,

$\begin{split}\overline{\delta
G^{(2)}}(k,t)&=(2\pi_{1}^{2}+V\pi_{2})\bar{f}_{t}*\bar{f}_{t}+2\pi_{1}\bar{f}_{t}*\overline{H^{(1)}[f]}+V\pi_{1}(k\bar{f})*(k\bar{f})+\overline{H^{(2)}[f,f]}\,,\end{split}$
(155)

where the convolution operator is defined as
$f*g=\int\\!\mathrm{d}k^{\prime}\,f(k-k^{\prime})g(k^{\prime})$ and

$\begin{split}\overline{H^{(2)}[f,f]}&=2I_{1}[\bar{f}*I_{1}[\bar{f}]]-\frac{1}{2}I_{1}[I_{1}[\bar{f}*\bar{f}]]-\bar{f}*I_{1}[I_{1}[\bar{f}]]+I_{1}[\bar{f}]*I_{1}[\bar{f}]\\\
&+\frac{V}{2}I_{2}[\bar{f}*\bar{f}]-V\bar{f}*I_{2}[\bar{f}]\,.\end{split}$

Using the identities $\pi_{1}^{\prime}(V)=V\pi_{2}$ and
$g^{\prime\prime}(V)/2g(V)=2\pi_{1}^{2}+\pi_{1}^{\prime}(V)$ we obtain

$G(z,t)=G_{r}g(V)\left(1+\frac{g^{\prime}(V)}{g(V)}f_{t}+\frac{g^{\prime\prime}(V)}{2g(V)}f_{tt}-\frac{g^{\prime}(V)}{g(V)}\frac{V}{2}f_{z}^{2}+H^{(1)}[f]+H^{(2)}[f,f]+\order{f^{3}}\right)\,,$
(156)

which is equivalent to

$G(z,t)=G_{r}g(V_{\perp})\left(1+H^{(1)}[f]+H^{(2)}[f,f]+\order{f^{3}}\right)\,.$
(157)

## Acknowledgments

This work was supported by the International Research Project “Non-Equilibrium
Physics of Complex Systems” (IRP-PhyComSys, France-Israel).

## References

  * Adda-Bedia et al. (2013) Adda-Bedia, M., Arias, R.E., Bouchbinder, E., Katzav, E., 2013\.  Dynamic stability of crack fronts: Out-of-plane corrugations.  Physical Review Letters 110, 14302\.  doi:10.1103/PhysRevLett.110.014302.
  * Adda-Bedia et al. (2006) Adda-Bedia, M., Katzav, E., Vandembroucq, D., 2006.  Second-order variation in elastic fields of a tensile planar crack with a curved front.  Physical Review E 73, 35106\.
  * Albertini et al. (2021) Albertini, G., Lebihain, M., Hild, F., Ponson, L., Kammer, D.S., 2021.  Effective Toughness of Heterogeneous Materials with Rate-Dependent Fracture Energy.  Physical Review Letters 127, 035501\.  doi:10.1103/PhysRevLett.127.035501.
  * Amestoy et al. (1981) Amestoy, M., Bui, H.D., Labbens, R., 1981.  On the definition of local path independent integrals in three-dimensional crack problems.  Mechanics Research Communications 8, 231–236.  doi:10.1016/0093-6413(81)90058-6.
  * Bayart et al. (2018) Bayart, E., Svetlizky, I., Fineberg, J., 2018.  Rupture Dynamics of Heterogeneous Frictional Interfaces.  Journal of Geophysical Research: Solid Earth 123, 3828–3848.  doi:10.1002/2018JB015509.
  * Bedford et al. (2022) Bedford, J.D., Faulkner, D.R., Lapusta, N., 2022.  Fault rock heterogeneity can produce fault weakness and reduce fault stability.  Nature Communications 13, 326\.  doi:10.1038/s41467-022-27998-2.
  * Bleyer and Molinari (2017) Bleyer, J., Molinari, J.F., 2017\.  Microbranching instability in phase-field modelling of dynamic brittle fracture.  Applied Physics Letters 110, 151903\.
  * Chen et al. (2015) Chen, C.H., Cambonie, T., Lazarus, V., Nicoli, M., Pons, A.J., Karma, A., 2015\.  Crack front segmentation and facet coarsening in mixed-mode fracture.  Physical review letters 115, 265503\.
  * Chopin et al. (2015) Chopin, J., Boudaoud, A., Adda-Bedia, M., 2015.  Morphology and dynamics of a crack front propagating in a model disordered material.  Journal of the Mechanics and Physics of Solids 74, 38–48.  doi:10.1016/j.jmps.2014.10.001.
  * Démery et al. (2014) Démery, V., Rosso, A., Ponson, L., 2014.  From microstructural features to effective toughness in disordered brittle solids.  EPL (Europhysics Letters) 105, 34003\.  doi:10.1209/0295-5075/105/34003.
  * Dodds et al. (1988) Dodds, R.H., Carpenter, W.C., Sorem, W.A., 1988.  Numerical evaluation of a 3-D J-integral$<$math$><$mtext$>$J-$<$/mtext$><$mtext$>$integral$<$/mtext$><$/math$>$ and comparison with experimental results for a 3-Point bend specimen.  Engineering Fracture Mechanics 29, 275–285.  doi:10.1016/0013-7944(88)90017-3.
  * El Kabir et al. (2018) El Kabir, S., Dubois, F., Moutou Pitti, R., Recho, N., Lapusta, Y., 2018.  A new analytical generalization of the J and G-theta integrals for planar cracks in a three-dimensional medium.  Theoretical and Applied Fracture Mechanics 94, 101–109.  doi:10.1016/j.tafmec.2018.01.004.
  * Eriksson (2002) Eriksson, K., 2002.  A domain independent integral expression for the crack extension force of a curved crack in three dimensions.  Journal of the Mechanics and Physics of Solids 50, 381–403.  doi:10.1016/S0022-5096(01)00059-X.
  * Eshelby (1969) Eshelby, J.D., 1969.  The elastic field of a crack extending non-uniformly under general anti-plane loading.  Journal of the Mechanics and Physics of Solids 17, 177–199.
  * Fineberg et al. (2003) Fineberg, J., Sharon, E., Cohen, G., 2003.  Crack front waves in dynamic fracture.  International Journal of Fracture 121, 55–69.
  * Freund (1990) Freund, L.B., 1990.  Dynamic Fracture Mechanics.  Cambridge University Press, Cambridge; New York.
  * Geubelle and Rice (1995) Geubelle, P.H., Rice, J.R., 1995\.  A spectral method for three-dimensional elastodynamic fracture problems.  Journal of the Mechanics and Physics of Solids 43, 1791–1824.  doi:10.1016/0022-5096(95)00043-I.
  * Gounon et al. (2022) Gounon, A., Latour, S., Letort, J., El Arem, S., 2022\.  Rupture Nucleation on a Periodically Heterogeneous Interface.  Geophysical Research Letters 49, e2021GL096816.  doi:10.1029/2021GL096816.
  * Heizler and Kessler (2015) Heizler, S.I., Kessler, D.A., 2015\.  Microbranching in mode-I fracture using large-scale simulations of amorphous and perturbed-lattice models.  Physical Review E 92, 12403\.
  * Henry (2019) Henry, H., 2019.  Limitations of the modelling of crack propagating through heterogeneous material using a phase field approach.  Theoretical and Applied Fracture Mechanics 104, 102384.  doi:10.1016/j.tafmec.2019.102384.
  * Henry and Adda-Bedia (2013) Henry, H., Adda-Bedia, M., 2013\.  Fractographic aspects of crack branching instability using a phase-field model.  Physical Review E 88, 60401\.
  * Kolvin et al. (2017) Kolvin, I., Fineberg, J., Adda-Bedia, M., 2017.  Nonlinear Focusing in Dynamic Crack Fronts and the Microbranching Transition.  Physical review letters 119, 215505\.
  * Latour et al. (2013) Latour, S., Voisin, C., Renard, F., Larose, E., Catheline, S., Campillo, M., 2013\.  Effect of fault heterogeneity on rupture dynamics: An experimental approach using ultrafast ultrasonic imaging.  Journal of Geophysical Research: Solid Earth 118, 5888–5902.  doi:10.1002/2013JB010231.
  * Lebihain et al. (2020) Lebihain, M., Leblond, J.B., Ponson, L., 2020.  Effective toughness of periodic heterogeneous materials: The effect of out-of-plane excursions of cracks.  Journal of the Mechanics and Physics of Solids 137, 103876.  doi:10.1016/j.jmps.2020.103876.
  * Lebihain et al. (2021) Lebihain, M., Roch, T., Violay, M., Molinari, J.F., 2021\.  Earthquake Nucleation Along Faults With Heterogeneous Weakening Rate.  Geophysical Research Letters 48, e2021GL094901.  doi:10.1029/2021GL094901.
  * Leblond et al. (2011) Leblond, J.B., Karma, A., Lazarus, V., 2011.  Theoretical analysis of crack front instability in mode I+{}III{}.  Journal of the Mechanics and Physics of Solids 59, 1872–1887.  doi:10.1016/j.jmps.2011.05.011.
  * Leblond et al. (2012) Leblond, J.B., Patinet, S., Frelat, J., Lazarus, V., 2012\.  Second-order coplanar perturbation of a semi-infinite crack in an infinite body.  Engineering Fracture Mechanics 90, 129–142.
  * Leguillon (2014) Leguillon, D., 2014.  An attempt to extend the 2D coupled criterion for crack nucleation in brittle materials to the 3D case.  Theoretical and Applied Fracture Mechanics 74, 7–17.  doi:10.1016/j.tafmec.2014.05.004.
  * Lubomirsky and Bouchbinder (2023) Lubomirsky, Y., Bouchbinder, E., 2023\.  Quenched disorder and instability control dynamic fracture in three dimensions.  doi:10.48550/arXiv.2311.11692, arXiv:2311.11692.
  * Mirkhalaf et al. (2014) Mirkhalaf, M., Dastjerdi, A.K., Barthelat, F., 2014.  Overcoming the brittleness of glass through bio-inspiration and micro-architecture.  Nature Communications 5, 3166\.  doi:10.1038/ncomms4166.
  * Möller and Bitzek (2015) Möller, J.J., Bitzek, E., 2015\.  On the influence of crack front curvature on the fracture behavior of nanoscale cracks.  Engineering Fracture Mechanics 150, 197–208.  doi:10.1016/j.engfracmech.2015.03.028.
  * Morrissey and Rice (1998) Morrissey, J.W., Rice, J.R., 1998\.  Crack front waves.  Journal of the Mechanics and Physics of Solids 46, 467–487.
  * Morrissey and Rice (2000) Morrissey, J.W., Rice, J.R., 2000\.  Perturbative simulations of crack front waves.  Journal of the Mechanics and Physics of Solids 48, 1229–1251.  doi:10.1016/S0022-5096(99)00069-1.
  * Movchan et al. (1998) Movchan, A.B., Gao, H., Willis, J.R., 1998.  On perturbations of plane cracks.  International Journal of Solids and Structures 35, 3419–3453.
  * Norris and Abrahams (2007) Norris, A.N., Abrahams, I.D., 2007\.  A multiple-scales approach to crack-front waves.  Journal of Engineering Mathematics 59, 399–417.
  * Perrin and Rice (1994) Perrin, G., Rice, J.R., 1994\.  Disordering of a dynamic planar crack front in a model elastic medium of randomly variable toughness.  Journal of the Mechanics and Physics of Solids 42, 1047–1064.
  * Pons and Karma (2010) Pons, A.J., Karma, A., 2010\.  Helical crack-front instability in mixed-mode fracture.  Nature 464, 85–89.
  * Ponson (2009) Ponson, L., 2009.  Depinning Transition in the Failure of Inhomogeneous Brittle Materials.  Physical Review Letters 103, 055501\.  doi:10.1103/PhysRevLett.103.055501.
  * Ramanathan (1997) Ramanathan, S., 1997.  Crack Propagation through Heterogeneous Media.  Ph.D. thesis. Harvard University.
  * Ramanathan and Fisher (1997) Ramanathan, S., Fisher, D.S., 1997\.  Dynamics and instabilities of planar tensile cracks in heterogeneous media.  Physical Review Letters 79, 877–880.  doi:10.1103/PhysRevLett.79.877.
  * Rice et al. (1994) Rice, J.R., Ben-Zion, Y., Kim, K.S., 1994.  Three-dimensional perturbation solution for a dynamic planar crack moving unsteadily in a model elastic solid.  Journal of the Mechanics and Physics of Solids 42, 813–843.
  * Roch et al. (2022) Roch, T., Barras, F., Geubelle, P.H., Molinari, J.F., 2022\.  cRacklet: A spectral boundary integral method library for interfacial rupture simulation.  Journal of Open Source Software 7, 3724.  doi:10.21105/joss.03724.
  * Roch et al. (2023) Roch, T., Lebihain, M., Molinari, J.F., 2023.  Dynamic Crack-Front Deformations in Cohesive Materials.  Physical Review Letters 131, 096101\.  doi:10.1103/PhysRevLett.131.096101.
  * Sagy et al. (2004) Sagy, A., Fineberg, J., Reches, Z., 2004.  Shatter cones: Branched, rapid fractures formed by shock impact.  Journal of Geophysical Research: Solid Earth 109\.
  * Schmittbuhl and Måløy (1997) Schmittbuhl, J., Måløy, K.J., 1997\.  Direct Observation of a Self-Affine Crack Propagation.  Physical Review Letters 78, 3888–3891.  doi:10.1103/PhysRevLett.78.3888.
  * Sharon et al. (2001) Sharon, E., Cohen, G., Fineberg, J., 2001.  Propagating solitary waves along a rapidly moving crack front.  Nature 410, 68–71.
  * Sharon et al. (2002) Sharon, E., Cohen, G., Fineberg, J., 2002.  Crack front waves and the dynamics of a rapidly moving crack.  Physical Review Letters 88, 85503\.
  * Stanchits et al. (2015) Stanchits, S., Burghardt, J., Surdi, A., 2015.  Hydraulic Fracturing of Heterogeneous Rock Monitored by Acoustic Emission.  Rock Mechanics and Rock Engineering 48, 2513–2527.  doi:10.1007/s00603-015-0848-1.
  * Steinhardt and Rubinstein (2022) Steinhardt, W., Rubinstein, S.M., 2022\.  How Material Heterogeneity Creates Rough Fractures.  Physical Review Letters 129, 128001\.  doi:10.1103/PhysRevLett.129.128001.
  * Svetlizky and Fineberg (2014) Svetlizky, I., Fineberg, J., 2014\.  Classical shear cracks drive the onset of dry frictional motion.  Nature 509, 205–208.
  * Vasoya et al. (2016) Vasoya, M., Unni, A.B., Leblond, J.B., Lazarus, V., Ponson, L., 2016.  Finite size and geometrical non-linear effects during crack pinning by heterogeneities: An analytical and experimental study.  Journal of the Mechanics and Physics of Solids 89, 211–230.
  * Willis (2013) Willis, J.R., 2013.  Crack front perturbations revisited.  International Journal of Fracture 184, 17–24.
  * Willis and Movchan (1995) Willis, J.R., Movchan, A.B., 1995\.  Dynamic weight functions for a moving crack. I. Mode I loading.  Journal of the Mechanics and Physics of Solids 43, 319–341.
  * Willis and Movchan (1997) Willis, J.R., Movchan, A.B., 1997\.  Three-dimensional dynamic perturbation of a propagating crack.  Journal of the Mechanics and Physics of Solids 45, 591–610.
  * Xia et al. (2012) Xia, S., Ponson, L., Ravichandran, G., Bhattacharya, K., 2012\.  Toughening and asymmetry in peeling of heterogeneous adhesives.  Physical Review Letters 108, 1–5.  doi:10.1103/PhysRevLett.108.196101, arXiv:1203.3634.