Analysis of cell size effects in atomistic crack propagation

TL;DR

This paper investigates how cell size influences crack propagation in crystalline materials, revealing bifurcation behavior and providing mathematical proofs for finite-cell approximation convergence.

Contribution

It introduces a bifurcation analysis framework for crack propagation, demonstrating the stress intensity factor as a key parameter and analyzing finite-cell approximation convergence.

Findings

Stress intensity factor acts as a bifurcation parameter

Bifurcation diagram exhibits a periodic snaking curve

Finite-cell approximations converge to the exact bifurcation diagram

Abstract

We consider crack propagation in a crystalline material in terms of bifurcation analysis. We provide evidence that the stress intensity factor is a natural bifurcation parameter, and that the resulting bifurcation diagram is a periodic "snaking curve". We then prove qualitative properties of the equilibria and convergence rates of finite-cell approximations to the "exact" bifurcation diagram.

Figures

Equations

$$E(u):={\int }_{\mathcal{D}\setminus {\Gamma }_{\mathcal{D}}}W(\nabla u)dx,$$

$$-\mathrm{div}\left(\mathbb{C}:\nabla u\right)$$

$$\left(\mathbb{C}:\nabla u\right)\nu$$

$$-\Delta u$$

$$\nabla u\cdot \nu$$

$${\Gamma }_0:=\{(x_1,0)|x_1\le 0\}.$$

$${\hat{u}}_k(x)=k\sqrt{r}\sin \frac{\theta }{2},$$

$$\Lambda:=\big{\{}l-(\tfrac{1}{2},\tfrac{1}{2})\,\big{|}\,l\in\mathbb{Z}^{2}\big{\}}.$$

$$\Gamma_{\pm}:=\big{\{}m\in\Lambda\,\big{|}\,m_{1}<0\,\text{ and }\,m_{2}=\pm\textstyle{\frac{1}{2}}\big{\}}$$

$$\mathcal{R}=\left\{\pm e_1,\pm e_2\right\}.$$

$$\tilde{\mathcal{R}}(m):=\{\begin{array}{ll}\mathcal{R}& \text{for }m\notin ({\Gamma }_{+}\cup {\Gamma }_{-}),\\ \mathcal{R}\setminus \left\{\mp e_2\right\}& \text{for }m\in {\Gamma }_{\pm }.\end{array}$$

$$\big{(}Du(m)\big{)}_{\rho}{\,:=\,}D_{\rho}u(m)\text{ and }\big{(}\tilde{D}u(m)\big{)}_{\rho}:=\begin{cases}D_{\rho}u(m)\quad&\text{if }\rho\in\tilde{\mathcal{R}}(m),\\ 0\quad&\text{if }\rho\not\in\tilde{\mathcal{R}}(m).\end{cases}$$

$$\begin{array}{c}{\dot{\mathcal{H}}}^1:=\left\{u:\Lambda \to \mathbb{R}|\tilde{D}u\in {\ell }^2\text{and}u(\frac{1}{2},\frac{1}{2})=0\right\},\\ \begin{array}{rl}\text{which has associated norm}& \Vert u{\Vert }_{{\dot{\mathcal{H}}}^1}:=\Vert \tilde{D}u{\Vert }_{{\ell }^2}={\left(\sum _{m\in \Lambda }|\tilde{D}u(m){|}^2\right)}^{\text{nicefrac}12}\\ \text{ and inner product}& (u,v{)}_{{\dot{\mathcal{H}}}^1}:=\sum _{m\in \Lambda }\tilde{D}u(m)\cdot \tilde{D}v(m).\end{array}\end{array}$$

$${\mathcal{H}}^{\mathrm{c}}:=\{u:\Lambda \to \mathbb{R}|\mathrm{supp}(Du)\text{ is compact}\}.$$

$$\mathcal{E}(u,k)=\sum_{m\in\Lambda}V\big{(}D\hat{u}_{k}(m)+Du(m)\big{)}-V\big{(}D\hat{u}_{k}(m)\big{)}.$$

$$V(Du(m))=\sum_{\rho\in\mathcal{R}}\phi\big{(}D_{\rho}u(m)\big{)},$$

$$\varphi (0)=0,{\varphi }^{\prime }(0)=0,{\varphi }^{\prime \prime }(0)=1,\text{and}{\varphi }^{\prime \prime \prime }(0)=0.$$

$$\text{there exists }{R}_{\varphi }>0\text{ such that }{\varphi }^{\prime }(r)=0\forall r\text{ with }|r|\ge R_{\varphi }.$$

$$S:=\big{\{}(u,k)\in\dot{\mathcal{H}}^{1}\times\mathbb{R}\,\big{|}\,\delta_{u}\mathcal{E}(u,k)=0\in(\dot{\mathcal{H}}^{1})^{*}\big{\}},$$

$$⟨{\delta }_u\mathcal{E}(u,k),v⟩=\sum _{m\in \Lambda }\delta V(D{\hat{u}}_k(m)+Du(m))\cdot Dv(m).$$

$$H_0:={\delta }_{uu}^2\mathcal{E}({\bar{u}}_0,{\bar{k}}_0):{\dot{\mathcal{H}}}^1\to ({\dot{\mathcal{H}}}^1{)}^{\ast }\text{ is an isomorphism.}$$

$$H_s:={\delta }_{uu}^2\mathcal{E}({\bar{u}}_s,{\bar{k}}_s):{\dot{\mathcal{H}}}^1\to ({\dot{\mathcal{H}}}^1{)}^{\ast }.$$

$${\delta }_{uk}^2\mathcal{E}({\bar{u}}_b,{\bar{k}}_b)[{\gamma }_b,1]$$

$${\delta }_{uuu}^3\mathcal{E}({\bar{u}}_b,{\bar{k}}_b)[{\gamma }_b,{\gamma }_b,{\gamma }_b]$$

$$\mathcal{B}(s):=({\bar{u}}_s,{\bar{k}}_s),$$

$$⟨H_{s}v,v⟩\ge c\Vert v{\Vert }_{{\dot{\mathcal{H}}}^1}^2$$

$${\mathcal{I}}_{\mathrm{pt}}:=\bigcup _{k=1}^{M/2}{I}_k\subset [0,1]\text{and}{\mathcal{I}}_{\mathrm{pos}}:=[0,1]\setminus {\mathcal{I}}_{\mathrm{pt}},$$

$${\mathcal{B}}_{\mathrm{pt}}:=\mathcal{B}({\mathcal{I}}_{\mathrm{pt}})\text{and}{\mathcal{B}}_{\mathrm{pos}}:=\mathcal{B}({\mathcal{I}}_{\mathrm{pos}})$$

$$s\in {\mathcal{I}}_{\mathrm{pos}}⟹U_s={\dot{\mathcal{H}}}^1.$$

$$H_s{\gamma }_s={\mu }_{s}J{\gamma }_s,$$

$$⟨H_s{\gamma }_s,v⟩={\mu }_s({\gamma }_s,v{)}_{{\dot{\mathcal{H}}}^1}\text{for all }v\in {\dot{\mathcal{H}}}^1.$$

$$\big{|}\tilde{D}\bar{u}_{s}(l)\big{|}\leq C|l|^{-3/2+\delta}.$$

Full text

MB is supported by EPSRC as part of the MASDOC DTC, Grant No. EP/HO23364/1.

TH is supported by the Leverhulme Trust through Early Career Fellowship ECF-2016-526

CO is supported by EPSRC Grant EP/R043612/1

1. Introduction

A fundamental task of materials modelling is to understand the process of failure, which is often facilitated by crack propagation. Cracks (and other defects) initiate and propagate via atomistic mechanisms, which renders the task of creating accurate and efficient simulations of this phenomenon on a large scale particularly difficult [BKG15]. In addition, many of the simulation techniques in operation today rely on simplifying assumptions which are often phenomenological; for example, it is generally unclear under which conditions continuum models become invalid [GSHY10].

There is thus a need for a robust mathematical theory of crack propagation at an atomistic scale, providing a rigorous grounding for a subsequent study of bottom-up multiscale and coarse-grained models. In [BHO19] we began to lay the foundation of such a theory by formulating the equilibration problem on a lattice in the presence of a crack as a variational problem on an appropriate discrete Sobolev space, and establishing existence, local uniqueness and stability of equilibrium displacements for small loading parameters. Crucially, we also established decay properties of the lattice Green’s function in crack geometry, which enabled us to prove qualitatively sharp far-field decay estimates of the atomistic core contribution to the equilibrium fields in order to quantify the “range” of atomistic effects. This work relied on and extended the recent rigorously formalised atomistic theory of single localised defects in crystalline structures [EOS16, HO14, BBO19].

The purpose of the present work is to go beyond the small-loading regime and introduce a key component missing in [BHO19], demonstrating that crack propagation, facilitated by bond-breaking events, can be described in the framework of [BHO19]. The mathematical tools we exploit to do so are taken from bifurcation theory in Banach spaces [CST00]. While this idea has already been explored numerically in [Li13, Li14], a key new conceptual insight is that the stress intensity factor (SIF), which acts as a measure of stability in continuum fracture can be interpreted as the “loading parameter” on the atomistic crack through the far-field boundary condition allowing us to obtain rigorous results about cell size effects.

More specifically, we model the equilibration of an atomistic crack embedded in an infinite homogeneous crystal as a variational problem with the continuum linearised elasticity (CLE) solution as the far-field boundary condition. The SIF enters the model as a scaling parameter multiplying the CLE solution, and so varying this naturally leads to a bifurcation diagram. Moreover, the fact that the CLE crack equilibrium displacement does not belong to the energy space suggests that the bifurcation diagram consists solely of regular points and quadratic fold points, at which the equilibria found transition from being linearly stable to linearly unstable (or vice versa).

This observation and the numerical evidence we obtain together motivate structural assumptions on the bifurcation diagram: we assume (and confirm numerically) that it is a ‘snaking curve’ [TD10] with the stability of solutions changing at each bifurcation point. In particular, under our assumptions, a jump from one stable segment to another captures the propagation of the crack through one lattice cell, with the unstable segment that is crossed in that jump being a corresponding saddle point, which represents the energetic barrier which must be overcome for crack propagation to occur at a given value of the SIF. This allows us to capture the phenomenon of ‘lattice trapping’ [THR71, GC00], a term which refers to the idea that in discrete models of fracture there can exist a range of values of SIF for which the crack remains locally stable despite being above or below the critical Griffith stress.

As in [BHO19], we avoid significant technicalities by restricting the analysis to a two–dimensional square lattice with nearest neighbour pair interaction. The notable difference in the models considered is that in [BHO19], in order to prove that the variational problem is well-posed, the bonds crossing the crack were explicitly removed from the interaction; by contrast, in the present paper they are included in the interaction range, and instead, the fact that they are effectively broken is encoded in the interatomic potential. This gives rise to a physically realistic periodic bifurcation diagram, for which we subsequently prove regularity results both in terms of its smoothness as a submanifold of an appropriate space, as well as uniform spatial regularity of the equilibria along the corresponding solution path.

Our results for the infinite lattice model naturally lead to an investigation of the numerical approximation of these solutions on a finite-domain, and we use the technical tools established in [BHO19] to establish sharp convergence rates as the domain radius tends to infinity. A notable novelty is that our results apply uniformly to finite segments of the bifurcation diagram; moreover, we establish a superconvergence result for the critical values of the SIF at which fold points occur. Since the unstable segments of the bifurcation diagram correspond to index–1 saddle points of the energy, our work in this regard also extends the convergence results of [BO18] for saddle point configurations of point defects and suggests possible future extensions to a full transition state analysis [Eyr35, HTB90, Wig97, Ber13, BDO18].

1A Outline:

In Section 2 we provide a detailed motivation for our work, introducing the model for crack propagation in the anti-plane setup, describing the underlying assumptions, and providing a statement of the main results about the model and its numerical approximation. In Section 2.1 we give a brief overview of the continuum mechanics context and describe how it motivates our work, and in Section 2.2 we discuss the discrete kinematics of the atomistic model. Then, in Section 2.3 the key assumptions are presented and discussed, and the novel components of the theory, in particular the role of the stress intensity factor, are highlighted. The main results of the paper are also stated. Section 2.4 is dedicated to the finite-domain approximation of the problem, with sharp convergence results stated, including the superconvergence result for the critical values of the stress intensity factor. In Section 2.5 we present a numerical setup employed to compute bifurcation paths, enabling us to numerically verify the sharpness of our results with respect to regularity and rate of convergence. Section 3 then provides a discussion about the significance of our results, and the proofs of the main results are given in Section 4.

2. Main results

2.1. Motivation

The principal motivation for our work stems from the following limitation of the continuum elasticity approaches to static crack problems. Consider a domain $\mathcal{D}\subset\mathbb{R}^{2}$, representing a cross-section of a three-dimensional elastic body, with a crack set $\Gamma_{\mathcal{D}}\subset\mathcal{D}$. Given a material-specific strain energy density function $W\,:\,\mathbb{R}^{2\times k}\to\mathbb{R}\cup\{+\infty\}$ (where $k\in\{1,2,3\}$ depending on the loading mode), c.f. [LLS89], one can hope to find a non-trivial equilibrium displacement $u\,:\,\mathcal{D}\to\mathbb{R}^{k}$ accommodating the presence of a crack by minimising the continuum energy given by

[TABLE]

over a suitable function space. In line with CLE, one can approximate $W$ by its expansion around zero to second order and obtain the associated equilibrium equation

[TABLE]

supplied with a suitable boundary condition coupling to the bulk [Fre90]. Here $\mathbb{C}$ is the elasticity tensor with entries $\mathbb{C}^{j\beta}_{i\alpha}:=\partial_{i\alpha j\beta}W(0)$.

It is well-known that regardless of the details of the geometry of $\mathcal{D}$ and $\Gamma_{\mathcal{D}}$, near the crack tip, the gradients of solutions to (2.1)-(2.2) exhibit a persistent $1/\sqrt{r}$ behaviour, where $r$ is the distance from the crack tip, c.f. [Ric68]. The singularity at the crack tip implies the failure of CLE to accurately describe a small region around it where atomistic (nonlinear and discrete) effects dominate. This near-tip nonlinear zone is argued to exhibit autonomy [Fre90, Bro99], meaning that the state of the system in the vicinity of the singular field is determined uniquely by the value of the SIF, and therefore systems with the same SIF but different geometries will behave similarly within the near-tip nonlinear zone.

In order to better understand the microscopic features of this zone, we may exploit the spatial invariance of elasticity and zoom in on the region near the crack tip by performing a spatial rescaling $R\,u(x/R)$, which leads to a simplified geometry of an infinite domain with a half-infinite straight crack line, as illustrated in Figure 1.

In what follows we focus on Mode III cracks, restricting to anti-plane displacements $u:\mathbb{R}^{2}\to\mathbb{R}$. Under natural assumptions on the stored energy density $W$ which result from coupling it with a corresponding interatomic potential, as discussed in [BHO19], the set of equations (2.1)-(2.2) reduces in the simplified geometry to

[TABLE]

where

[TABLE]

This PDE has a canonical solution, as discussed in e.g. [SJ12], given by

[TABLE]

with $(r,\theta)$ representing standard cylindrical polar coordinates centred at the crack tip. The scalar parameter $k$ corresponds to the (rescaled) stress intensity factor (SIF) [Law93].

As we increase the spatial rescaling parameter $R$, we eventually approach the atomic lengthscale at which the hypothesis that the material behaves as a continuum breaks down. We must thus speak of atomic displacements and finite differences rather than differential operators, and consider an atomistic model supplied with the function $\hat{u}_{k}$ as a far-field boundary condition. As such, in the model we describe below, the function $\hat{u}_{k}$ will act as a CLE ‘predictor’, representing the behaviour of the material in the far field away from the crack tip.

2.2. Discrete kinematics

The atomistic setup is similar to the one introduced in [BHO19]; here, we recall it in detail and highlight new concepts. Let $\Lambda$ denote the shifted two dimensional square lattice defined as

[TABLE]

We consider a crack opening along $\Gamma_{0}$ defined in (2.5), and distinguish the lines that include lattice points directly above and below $\Gamma_{0}$. These are defined as

[TABLE]

and we refer to Figure 2 for a visualisation of the setup.

For the purposes of our analysis, it is helpful to consider two notions of interaction neighbourhood for lattice points. First, the nearest neighbour (NN) directions of the homogeneous square lattice are given by

[TABLE]

Second, we modify these interaction neighbourhoods by disregarding the directions across the crack, as these bonds are effectively already broken; for any $m\in\Lambda$, we therefore define

[TABLE]

For an anti–plane displacement defined on the lattice $u\,\colon\,\Lambda\to\mathbb{R}$, we define the finite difference operator as $D_{\rho}u(x):=u(x+\rho)-u(x)$ and introduce two notions of the discrete gradient, denoted by $Du(m),\,\tilde{D}u(m)\,\in\mathbb{R}^{\mathcal{R}}$ and defined as

[TABLE]

We note that since $|\mathcal{R}|=4$, $\mathbb{R}^{\mathcal{R}}$ is a $4$-dimensional space indexed by each member of $\mathcal{R}$. It therefore follows from (2.8) that $Du$ corresponds to homogeneous NN interactions, whereas $\tilde{D}u$ reflects a defective lattice, as when $m\in\Gamma_{\pm}$, the components of $\tilde{D}u(m)$ which correspond to erased lattice directions are always zero.

The removal of NN bonds and subsequent introduction of the discrete gradient operator $\tilde{D}$ allows us to define the appropriate discrete energy space (discrete Sobolev space) for handling arbitrarily large differences in the far–field displacements across the crack,

[TABLE]

The choice to restrict $\textstyle{u(\frac{1}{2},\frac{1}{2})}=0$ ensures that only one constant displacement lies in the space, making $\|\cdot\|_{\dot{\mathcal{H}}^{1}}$ a norm.

It is also helpful to introduce the space of compactly supported displacements,

[TABLE]

Remark 2.1*.*

To avoid future confusion, we note that compared to [BHO19] the definitions of $D$ and $\tilde{D}$ have been swapped to accommodate the differing nature of the two papers. In [BHO19] the interactions across the crack are always explicitly excluded, leaving little need for this explicit distinction. In the present work the distinction is crucial and we opted for $D$ to denote the usual intuitive notion of the discrete gradient. Furthermore, a similar change in notation occurs for $\mathcal{R}$ and $\tilde{\mathcal{R}}$. Note, however, that in both papers the definition of $\dot{\mathcal{H}}^{1}$ remains the same.

2.3. Analysis of the model

We modify the theory developed in [BHO19] for a small-load anti-plane crack to frame it in the context of bifurcation theory. We consider the energy difference functional $\mathcal{E}\,:\,\dot{\mathcal{H}}^{1}\times\mathbb{R}\to\mathbb{R}$ supplied with CLE solution as a far–field boundary condition:

[TABLE]

Here $V\,\colon\,\mathbb{R}^{\mathcal{R}}\to\mathbb{R}$ is a suitable interatomic site potential and $\hat{u}_{k}\,\colon\,\Lambda\to\mathbb{R}$ is the CLE predictor, introduced in (2.6). The function $u\in\dot{\mathcal{H}}^{1}$ is a core correction, thus the total displacement is given by $\hat{u}_{k}+u$.

We assume the site potential to be a NN pair-potential of the form

[TABLE]

with $\phi\in C^{\alpha}(\mathbb{R})$ for $\alpha\geq 5$. Assuming that the lattice is reflection symmetric in the anti–plane direction, we assume without loss of generality that

[TABLE]

The first and third assumptions may be made by subtracting an appropriate constant and rescaling the potential, while the second and fourth follow from the assumption of anti-plane symmetry, as discussed in [BHO19].

Note that we employ the homogeneous discrete gradient operator $Du$ in the definition of $\mathcal{E}$, while the ‘crack-aware’ gradient $\tilde{D}u$ is used to define the space $\dot{\mathcal{H}}^{1}$. This is helpful in the context of capturing crack propagation, since it enables us to consider displacements with arbitrarily large strains across the crack, but raises the issue that for any $m\in\Gamma_{\pm}$ and $\rho\not\in\tilde{\mathcal{R}}(m)$ crossing the crack surface, we have $D_{\rho}\hat{u}(m)\sim|m|^{1/2}$. Thus, in order for such $\mathcal{E}$ to be well-defined on $\dot{\mathcal{H}}^{1}\times\mathbb{R}$, we further assume that the pair-potential satisfies

[TABLE]

Such an assumption is sufficient and simplifies the exposition, but is by no means a necessary condition and can be easily replaced by an appropriate decay property (e.g. exponential or sufficiently high algebraic decay). In particular, under this assumption, we firstly prove the following result.

Theorem 2.2**.**

The energy difference functional $\mathcal{E}$ expressed in (2.10) is well-defined on $\dot{\mathcal{H}}^{1}\times\mathbb{R}$ and is $\alpha$-times continuously differentiable.

The proof is given in Section 4.2 and mostly relies on the analogous result in [BHO19], with the extra work needed to handle the now-included bonds across the crack.

The inclusion of the stress intensity factor $k$ as a variable in the definition of $\mathcal{E}$ allows us to employ bifurcation analysis to describe the propagation of the crack as a series of bifurcations, which we view as corresponding to bond-breaking events.

The primary task of our analysis is to characterise the set of critical points of the energy, $S$, defined as

[TABLE]

where $\delta_{u}\mathcal{E}\,:\,\dot{\mathcal{H}}^{1}\times\mathbb{R}\to(\dot{\mathcal{H}}^{1})^{*}$ is the partial Fréchet derivative given by

[TABLE]

For future reference, we summarize our notation for linear and multi-linear forms, in particular defining the meaning of $\langle\delta_{u}\mathcal{E}(u,k),\,v\rangle$. For any $n$-linear form $L$, we write $L[v_{1},\dots,v_{n}]$ to denote its evaluation at $v_{1},\dots,v_{n}$ and if $m<n$, then $L[v_{1},\dots,v_{m}]$ is the $(n-m)$-linear form $(w_{1},\dots,w_{n-m})\mapsto L[v_{1},\dots,v_{m},w_{1},\dots,w_{n-m}]$. For the sake of readability and only when there is no risk of confusion, we often write $\langle L,\,v_{1}\rangle$ for linear forms and $\langle L_{1}v_{1},\,v_{2}\rangle$ as well as $L_{1}v_{1}=L_{1}[v_{1}]$ for bilinear forms.

It is of particular interest to compute continuous paths contained in $S$, as it allows to characterise the response of the model to variations in SIF. This is often possible if we are able to identify one particular pair, say $(\bar{u}_{0},\bar{k}_{0})\in S$ and it can be further shown that it is a regular point, by which we mean

[TABLE]

In this case, a standard application of the Implicit Function Theorem [Lan99] yields existence of a locally unique path of solutions $(\bar{u}_{s},\bar{k}_{s})$ in the vicinity of $(\bar{u}_{0},\bar{k}_{0})$ which we will assume to be parametrised with an index $s\in\mathbb{R}$; exactly this strategy was used in [BHO19] to show existence of solutions in a static crack problem with crack bonds removed from the definition of $\mathcal{E}$, for $k$ small enough. We will set

[TABLE]

As we will see in the numerical examples of Section 2.5, beyond some critical value of $k$, bifurcations of the following type begin to occur.

Definition 2.3**.**

A (simple quadratic) fold point occurs at $(\bar{u}_{b},\bar{k}_{b})\in S$ if there exists $\gamma_{b}\in\dot{\mathcal{H}}^{1}$ such that ${\rm Ker}(H_{b})={\rm span}\{\gamma_{b}\}$,

[TABLE]

with formulae for these variations of energy given in (4.12) and (4.13), respectively.

A schematic representation of the idea behind Definition 2.3 is shown in Figure 3(a). As already discussed in the introduction, the fact that $\hat{u}_{k}\not\in\dot{\mathcal{H}}^{1}$ is key to (2.16) holding true, and suggests that a full bifurcation diagram is an infinite non-self-intersecting snaking curve [TD10], consisting solely of regular and fold points as shown in Figure 3(b). Our functional setup is well-suited to considering an arbitrary finite segment of it, so we begin with the following set of assumptions. We emphasise that all our subsequent results rely on the validity of these assumptions which are natural (see discussion below) but likely difficult to prove rigorously. Moreover, it is not guaranteed that Assumption 1 in particular is generic, but different potentials and loading geometries may indeed give rise to qualitatively different bifurcation diagrams.

Assumption 1**.**

There exists a bifurcation diagram in the form of an injective continuous path $\mathcal{B}\,:\,[0,1]\to\dot{\mathcal{H}}^{1}\times\mathbb{R}$ given by

[TABLE]

where ${\rm Im}(\mathcal{B})\subset S$ (defined in (2.13)) is compact and for each $s\in[0,1]$, $\mathcal{B}(s)$ is either a regular point, as in (2.14), or a fold point, as in Definition 2.3. We further assume that there are finitely many fold points occuring at $s\in\{b_{1},\dots,b_{M}\}\subset(0,1)$. In particular, this implies that ${\rm Im}(\mathcal{B})$ is a non–self–intersecting curve.

For future reference, if $f:\mathcal{B}\to X$, where $X$ is a Banach space, is differentiable, then we write $f^{\prime}_{s}:=\frac{d}{ds}f_{s}$.

Assumption 2**.**

There exists $c>0$ such that for each $s\in[0,1]$ there exists a subspace $U_{s}$ of $\dot{\mathcal{H}}^{1}$ of codimension at most 1 for which it holds that

[TABLE]

for all $v\in U_{s}$.

The fact that a succession of fold points occurs is assumed to be an inherent feature of the lattice and the potential in place, much as the existence of a solution to a static dislocation problem is assumed in [EOS16]. Assumption 2 ensures that each $\mathcal{B}(s)=(\bar{u}_{s},\bar{k}_{s})$ represents either a bifurcation point, a stable solution or an unstable solution which is an index–1 saddle point. This assumption is motivated by the fact that the anti-plane setup and lattice symmetry naturally binds the crack propagation to the $x_{1}$-axis, leaving little room for any more involved bifurcating behaviour. Moreover, this is also supported by numerical evidence presented in Section 2.5, in particular with Figure 5 clearly exhibiting the snaking curve structure of the bifurcation diagram. We also refer to Section 3.2 for a discussion about the periodicity of the bifurcation diagram which further justifies Assumptions 1 and 2. Finally we note that in [Li13] a similar numerical evidence is presented for a vectorial Mode I fracture model posed on a triangular lattice under Lennard-Jones potential.

As will be shown in Proposition 2.5, requiring that (2.17) holds ensures that a change in the stability of the solution occurs at each fold point. This implies that near bifurcation points and on the unstable segments the infimum of the spectrum of $H_{s}$ is an eigenvalue, which motivates the following decomposition of the parametrisation interval $[0,1]$: since we look at a finite segment of the full bifurcation diagram, we will assume for notational convenience that it starts on a stable segment and the number of fold points $M$ lying in ${\rm Im}(\mathcal{B})$ is even and define sets

[TABLE]

where $I_{k}:=(b_{2k-1}-\xi,b_{2k}+\xi)$ with $\xi>0$ small enough. The cases where $M$ is odd or we start on an unstable segment can be handled in an entirely analogous way. We refer to

[TABLE]

as the collection of segments of the bifurcation diagram with $\sigma_{p}(H_{s})\neq\emptyset$ (non–empty point spectrum) and $\sigma(H_{s})\subset[c,\infty)$ (positive spectrum, $c$ from Assumption 2), respectively. We note that both the unstable segments and neighbourhoods of the bifurcation points belong to $\mathcal{B}_{\rm pt}$, thus the constant $c$ in Assumption 2 can be chosen to be small enough so that

[TABLE]

We now establish some initial results about the model. First, a regularity result.

Proposition 2.4** (Regularity of the diagram).**

The set ${\rm Im}(\mathcal{B})\subset\dot{\mathcal{H}}^{1}\times\mathbb{R}$ is a one–dimensional $C^{\alpha-1}$ manifold.

This result will be proven in Section 4.2 and in particular entails that, without loss of generality, we may make the following assumption concerning the parametrisation $\mathcal{B}$.

Assumption 3**.**

The function $\mathcal{B}:[0,1]\to\dot{\mathcal{H}}^{1}\times\mathbb{R}$ is a constant speed, $C^{\alpha-1}$ parametrisation of the manifold ${\rm Im}(\mathcal{B})\subset\dot{\mathcal{H}}^{1}\times\mathbb{R}$.

We next prove a result concerning the existence of linearly unstable directions and corresponding negative eigenvalues for some sections of the bifurcation diagram.

Proposition 2.5** (Existence of an eigen-pair).**

Under Assumptions 1, 2 & 3, there exist $C^{\alpha-2}$ functions $\gamma\,:\,\mathcal{I}_{{\rm pt}}\to\dot{\mathcal{H}}^{1}$ and $\mu\,:\,\mathcal{I}_{{\rm pt}}\to\mathbb{R}$ such that

[TABLE]

where $H_{s}$ was defined in (2.15) and $J$ represents the Riesz mapping [Rud66], i.e. an isometric isomorphism between $\dot{\mathcal{H}}^{1}$ and $(\dot{\mathcal{H}}^{1})^{*}$, thus we can equivalently say that

[TABLE]

Furthermore, for $j=1,\dots,M$, we have $\mu_{b_{j}}=0$ with the corresponding eigenvector $\gamma_{b_{j}}$ introduced in Definition 2.3 and also $\mu_{b_{j}}^{\prime}\neq 0$, implying that a change of stability occurs at $s=b_{j}$.

This will be proven in Section 4.2.

We subsequently establish the following decay and regularity results for the atomistic core corrector, which rely on the precise characterisation of the lattice Green’s function for the anti-plane crack geometry developed in [BHO19].

Theorem 2.6** (Decay properties of solutions and eigenvectors).**

For any $s\in[0,1]$ and $l\in\Lambda$ with $|l|$ large enough it holds that for any $\delta>0$ the atomistic correction $\bar{u}_{s}$ satisfies

[TABLE]

If $s\in\mathcal{I}_{\rm pt}$, then the eigenvector $\gamma_{s}\in\dot{\mathcal{H}}^{1}$ from Proposition 2.5 satisfies

[TABLE]

In both cases $C$ is a generic constant independent of $s$.

As in the case of Theorem 2.2, we note that (2.24) can be proven as in [BHO19], except for an extra dificulty arising from the fact that bonds across the crack are now included. The estimate in (2.25) follows from a two-step argument that is similar in nature. Both these estimates will be proven in Section 4.2.

Remark 2.7*.*

The arbitrarily small $\delta>0$ in Theorem 2.6 appears due to a technical limitation of the method employed in [BHO19] to estimate the mixed second discrete derivative of the lattice Green’s function in the anti-plane crack geometry. We expect the result to hold for $\delta=0$ too, but this cannot be achieved with the current bootstrapping argument, which saturates at the known decay rate of the corresponding continuum Green’s function.

2.4. Approximation

As numerical simulations are naturally restricted to a computational domain of finite size, we now consider and analyse a finite-dimensional scheme that approximates the solution path $\mathcal{B}$ defined in (2.18) and establish rigorous convergence results.

The starting point is a computational domain $\Omega_{R}$ with $B_{R}\cap\Lambda\subset\Omega_{R}\subset\Lambda$ (where $B_{R}$ is a ball of radius $R$ centred at the origin) and the boundary condition prescribed as $\hat{u}$ on $\Lambda\setminus\Omega_{R}$. The approximation to (2.13) can thus be stated as a Galerkin approximation, that is we seek to characterise

[TABLE]

where

[TABLE]

We prove the following.

Theorem 2.8**.**

Under Assumptions 1, 2 & 3, there exists $R_{0}>0$, such that for all $R\geq R_{0}$, there exists a $C^{\alpha-1}$ approximate bifurcation path $\mathcal{B}_{R}\,:\,[0,1]\to\mathcal{H}_{R}^{0}\times\mathbb{R}$ given by

[TABLE]

where ${\rm Im}(\mathcal{B}_{R})\subset S_{R}$, such that for any $\beta>0$

[TABLE]

and

[TABLE]

where $\mathcal{B}(s)=(\bar{u}_{s},\bar{k}_{s})$ as in (2.18).

For the proof of this result, we refer to Section 4.3.

While the estimate in (2.27) appears to be almost sharp (our numerical results in Section 2.5 indicate that this estimate holds with $\beta=0$), more can be said about the approximation of the critical values of the stress intensity factor for which fold points occur.

Theorem 2.9**.**

For $R$ suffiently large, the approximate bifurcation path $\mathcal{B}_{R}$ from Theorem 2.8 contains $M$ fold points in the sense of Definition 2.3 occuring at $s\in\{b_{1}^{R},\dots,b_{M}^{R}\}\subset(0,1)$, and for each of these we have

[TABLE]

For the proof, we again refer to Section 4.3.

This superconvergence result is well-known in bifurcation theory [CST00], and is observed numerically in our tests in Section 2.5.

2.5. Numerical investigation

In this section we present results of numerical tests that confirm the rate of decay of $|\tilde{D}\bar{u}_{s}|$ and $|\tilde{D}\gamma_{s}|$ established in Theorem 2.6, as well as the convergence rates from Theorems 2.8 and 2.9. The computational setup is similar to the one described in [BBO19, Section 3], with $\Lambda$ and $\mathcal{R}$ as specified in Section 2.2 and the pair-potential given by

[TABLE]

We employ a pseudo-arclength numerical continuation scheme to approximate $\mathcal{B}_{R}$ [BCD*+*02]. To compute equilibria we employ a standard Newton scheme, terminating at an $\ell^{\infty}$-residual of $10^{-8}$.

Theorem 2.6 suggests that $\lvert\tilde{D}\bar{u}_{s}(l)\rvert\lesssim\lvert l\rvert^{-3/2}$ and $\lvert\tilde{D}\gamma_{s}(l)\rvert\lesssim\lvert l\rvert^{-3/2}$. This is verified in Figure 4. Theorem 2.8 suggests that in the supercell approximation of $\mathcal{B}$ in (2.18) we expect $\|\bar{u}_{s}^{R}-\bar{u}_{s}\|_{\dot{\mathcal{H}}^{1}}+|\bar{k}_{s}^{R}-\bar{k}_{s}|\sim\mathcal{O}(R^{-1/2})$, where $R$ is the size of the domain. To verify this numerically, we first compute $\mathcal{B}_{R}$ for $R=32,\dots,256$ via a pseudo–arclength continuation scheme. The results are shown in Figure 5, with stable segments plotted as solid lines and unstable segments as dashed lines. To measure the distance between the segments of the bifurcation diagram, we compute the Hausdorff distance [RW98] with respect to $\|\cdot\|_{\dot{\mathcal{H}}^{1}}$-norm between the critical points on $\mathcal{B}_{R}$ (for $R=32,\dots,90.51$) and on $\mathcal{B}_{R_{*}}$, where $R_{*}=256$. The result is shown in Figure 6(a).

Finally, we test the superconvergence result for the bifurcation points from Theorem 2.9, which predicts that $|\bar{k}^{R}_{b_{j}^{R}}-\bar{k}_{b_{j}}|\sim\mathcal{O}(R^{-1})$. To this end we accelerate the convergence of the sequence $\{\bar{k}^{R}_{b_{j}^{R}}\}$ (for $R=32,\dots,256$) by employing Richardson extrapolation [Ric11], thus giving us an approximate limit value for $\bar{k}_{b_{j}}$. The values of the approximate limits, as well as the $R^{-1}$ convergence is exhibited in Figure 6(b).

Remark 2.10*.*

The pair-potential $\phi$ defined in (2.29) does not satisfy the strong assumption of compact support of $\phi^{\prime}$ introduced in (2.12), but has the slightly weaker property of exponential decay in first derivative. It thus emphasises the already discussed point that (2.29) is by no means a necessary condition. It also emphasises the point that our results are of more practical applicability, despite the severe non-convexity of the energy landscape forcing us to introduce structural assumptions on the bifurcation diagram as opposed to proving them.

3. Conclusion and discussion

The results obtained here, in tandem with those of our previous paper [BHO19], introduce a mathematical framework in which a rigorous formulation and study of atomistic models of cracks and their propagation is possible. In particular, we have shown how the theory of atomistic modelling of defects developed in [EOS16, HO14, BBO19] can be combined with classical results from bifurcation theory [Kel77, BRR80] to study this problem, and a key insight is the identification of the stress intensity factor as a suitable bifurcation parameter which allows us to explore the energy landscape. While our results are of conditional nature in that they rely on assumptions that are reasonably justified and numerically verified, our analysis nonetheless sets earlier numerical work of [Li13, Li14] into a rigorous framework and provides a comprehensive explanation as to why the bifurcation diagram is a snaking curve.

While further work is needed to extend our analytical results to more general models (and particularly to the case of other crack modes), from a numerical perspective several aspects of our theory are of universal applicability. We therefore conclude by pointing out a series of interesting conclusions which arise from our analysis.

3.1. Applicability of the results to other crack models

Due to reasons explored in the concluding section of [BHO19], at present we do not see an easy way of rigorously extending our results beyond a 2D model for scalar Mode III crack with nearest neighbours pair interactions on a square lattice. On a practical level, however, it can be numerically verified that our theory is entirely applicable to any 2D model for scalar Mode III crack with arbitrary finite range interactions under an arbitrary feasible interatomic potential and the resulting bifurcation diagrams is indeed a snaking curve.

For vectorial models of other modes of crack the numerical method described is still feasible, but depending on the potential employed, numerical tests indicate that one can expect a more complicated bifurcating behaviour, not least because of surface relaxation effects. The numerical evidence in [Li13] in particular indicates that the structure of the bifurcation diagram described in this paper does apply to some vectorial models. We hope to investigate this in greater detail in the future.

3.2. Periodicity of the bifurcation diagram

In an infinite lattice, shifting the crack tip by one lattice spacing results in a physically identical configuration. Therefore, it is reasonable to conjecture that $\mathcal{B}_{R}$ in the limit as $R\to\infty$ generates a bifurcation diagram in which the critical points exist for values of the SIF $k$ within a fixed finite interval of admissible values. In Section 2.5 we have exploited the superconvergence result in Theorem 2.9 to test this hypothesis numerically and the results summarised in Figure 6 confirm this intuition, as the extrapolated limit values of SIF as $R\to\infty$ for every second bifurcation point are numerically identical, occuring at

[TABLE]

A translation invariance in the critical points further implies that, if we denote the centre of the CLE predictor by

[TABLE]

for some $\lambda\in\mathbb{Z}$, and define $\mathcal{E}_{\lambda}\,:\,\dot{\mathcal{H}}^{1}\times\mathbb{R}\to\mathbb{R}$ by

[TABLE]

then assuming $\mathcal{B}(s)=(\bar{u}_{s},\bar{k}_{s})\in\dot{\mathcal{H}}^{1}\times\mathbb{R}$ is a parametrisation as described in Section 2.3, we naturally have

[TABLE]

We further notice for any $s\in[0,1]$ the total displacement $y_{s}=\hat{u}_{\bar{k}_{s}}+\bar{u}_{s}$ can be rewritten as $y_{s}(m)=\hat{u}_{\bar{k}_{s}}(m-x_{\lambda})+w_{s,\lambda}(m)$, where $\lambda\in\mathbb{Z}$ and

[TABLE]

Crucially,

[TABLE]

and for any choice of $\lambda\in\mathbb{Z}$

[TABLE]

In other words, no matter which $x_{\lambda}$ we choose to centre the crack predictor $\hat{u}_{k}$ at, the same configuration $y_{s}=\hat{u}_{\bar{k}_{s}}+\bar{u}_{s}$ always remains an equilibrium and $w_{s,\lambda}$ exactly captures the resulting changes to the atomistic correction.

To be precise, let us fix some $s\in[0,1]$ thus giving us a pair $\mathcal{B}(s)=(\bar{u}_{s},\bar{k}_{s})$, let $K:=\bar{k}_{s}$ and further consider

[TABLE]

Equations (3.2) and (3.4) indicate that for any $s^{\prime}\in I_{K}$ for which $\bar{u}_{s^{\prime}}$ is a solution of the same type as $\bar{u}_{s}$ (either stable, or unstable or a bifurcation point), we can find a unique $\lambda\in\mathbb{Z}$ such that $\bar{u}_{s^{\prime}}=w_{s,\lambda}$.

In particular we note that the above strongly suggests that in some cases one may be able to prove results about periodicity and boundedness in $k$ of the bifurcation diagram, which would also have the additional benefit of proving that Assumptions 1 and 2 hold true. The particular difficulty, however, lies in the fact that without these structural assumptions we cannot easily conclude that $\bar{u}_{s^{\prime}}=w_{s,\lambda}$ for some unique $\lambda\in\mathbb{Z}$. This can potentially be proven under suitable (prohibitively restrictive) technical asssumptions on the potential, as explored for anti-plane screw dislocations in [HO15]. We hope to investigate this in future work.

3.3. Interplay between the stress intensity factor and the domain size and its effect on lattice trapping.

The tilt of bifurcation diagrams seen in Figure 5 indicates that the size of the domain heavily impacts the shape of the corresponding solution curve. Notably, each successive bond-breaking event has a different interval of admissible values of SIF associated to it and the corresponding unstable segments are much shorter than stable ones for small domain sizes. The fact that the influence of such finite-domain effects can still be observed for a fairly large $R$ can be explained by the very slow rate of convergence in Theorem 2.9.

In practice, one hopes to investigate crack propagation and associated energy barriers for a fixed value of SIF and subsequently compare it against other admissible choices of SIF to measure the strength of lattice trapping [THR71, GC00], measured by the relative height of the energy barrier. Our work indicates that such investigations are particularly challenging due to the extent to which finite–size effects dominate, an effect we observe to be strong even in the simple model considered here. Only a very large choice of truncation radius $R$ ensures that the resulting solution paths are close to the periodic results one expects in the full lattice case. It may be possible to overcome such difficulties by prescribing a more accurate predictor describing the far–field behaviour, in line with the idea of development of solutions introduced in [BBO19]: this is a clear direction for future investigation.

3.4. Identification of the correct bifurcation parameter.

It is interesting to note that varying the intuitively natural bifurcation parameter $\lambda$ introduced in (3.1) to reflect the crack tip at which the continuum prediction is centred in fact fails to capture the bifurcation phenomenon. This can be seen by considering $\tilde{\mathcal{E}}\,:\,\dot{\mathcal{H}}^{1}\times\mathbb{R}\to\mathbb{R}$ given by

[TABLE]

where $\hat{v}_{\lambda}=\hat{u}_{K}(\cdot-x_{\lambda})$ for some fixed SIF $K>0$. This fundamentally differs from the energy defined in (2.10), as in that case we have linear dependence of the displacement on $k$, $D\hat{u}_{k}(m)=kD\hat{u}_{1}(m)$, which in turn leads to a particular form of derivatives with respect to $k$; in particular, this ensures we have quadratic fold points. In $\tilde{\mathcal{E}}$, however, the dependence is inside $\hat{v}_{\lambda}$, thus the crucial derivative with respect $\lambda$ is given by

[TABLE]

where $f_{\lambda}(m)=\nabla\hat{u}_{K}(m_{1}-\lambda,m_{2})\cdot e_{1}$. In this case, as in (3.3), we can conclude $f_{\lambda}\in\dot{\mathcal{H}}^{1}$. This implies that a fold point cannot occur, as that would require that there exists $\gamma\in\dot{\mathcal{H}}^{1}$ such that

[TABLE]

for all $w\in\dot{\mathcal{H}}^{1}$, which further implies that

[TABLE]

since $f_{\lambda}\in\dot{\mathcal{H}}^{1}$. This breaches the defining property of a fold point given in Definition 2.3; in fact, it is not possible to drive a bifurcation in this way precisely because of the observation made in (3.3) and the resulting periodicity.

3.5. Parameter-driven analysis for other models and defects

An overarching idea of this paper is that a careful analysis of a crucial parameter involved in the model can reveal the energy landscape of the problem, and this can be particularly fruitful in the study of defect migration. In the case of a crack, using the SIF as a driving parameter naturally generalises to more complex fracture models, though in general there may be multiple SIFs.

It would be interesting to undertake future study to see whether such an analysis is applicable more widely to other defects. In the particular case of dislocations, nucleation and motion have been studied in the atomistic context in a number of recent studies, including [CFR06, ADLGP14, GAM15, ADLGP17, Hud17]. Since these defects are the carriers of plastic deformation, the study of their mobility is important, and natural candidates for the parameters in this case are the shear modulus and externally–applied stress [LLS89].

4. Proofs

4.1. Preliminaries

Our approach is based on two classic results from bifurcation analysis on Banach spaces, cf. [CST00], which we state in this section for convenience. The first result is known as ’ABCD Lemma’ and is adapted from [Kel77].

Lemma 4.1** (ABCD Lemma).**

Let $H$ be a Hilbert space with the dual $H^{*}$ and consider the linear operator $M\,:\,H\times\mathbb{R}\to H^{*}\times\mathbb{R}$ of the form

[TABLE]

where $A\,:\,H\to H^{*}$ is self-adjoint in the sense that $\langle Av,\,w\rangle=\langle Aw,\,v\rangle$ for all $v,w\in H$, $b\in H^{*}\setminus\{0\}$, $c\in H\setminus\{0\}$ and $d\in\mathbb{R}$. Then

1. (i)

if $A$ is an isomorphism from $H$ to $H^{*}$, then $M$ is an isomorphism between $H\times\mathbb{R}$ and $H^{*}\times\mathbb{R}$ if and only if $d-(c,A^{-1}b)_{H}\neq 0$; and 2. (ii)

if ${\rm dim}\,{\rm Ker}(A)={\rm codim}\,{\rm Range}(A)=1$ with ${\rm Ker}(A)={\rm span}\{\gamma\}$, then $M$ is an isomorphism if and only if $\langle b,\,\gamma\rangle\neq 0$ and $(c,\gamma)_{H}\neq 0$.

To state the second result we introduce the following setup: let $X$, $Y$ and $Z$ be real Banach spaces and $F\in C^{k}(U\times Y;Z)$ for some $k\geq 1$, where $U$ is a bounded open subset of $X$. The total derivative of $F$ at $(x,y)\in X\times Y$ is denoted $DF(x,y)\in\mathcal{L}(X\times Y,Z)$, with partial derivatives denoted $D_{x}F(x,y)\in\mathcal{L}(X,Z)$ and $D_{y}F(x,y)\in\mathcal{L}(Y,Z)$. We now state a version of [BRR80, Theorem 1] tailored to our setting.

Theorem 4.2**.**

Suppose a function $y\,:\,U\to Y$ is Lipschitz continuous with Lipschitz constant $c_{2}$, and there exist constants $c_{0}$ and $c_{1}$ and a monotonically increasing function $L_{1}:\mathbb{R}\to\mathbb{R}$ such that the following hypotheses are satisfied:

1. (i)

for any $x_{0}\in U$, D_{y}F\big{(}x_{0},y(x_{0})\big{)} is an isomorphism of $Y$ onto $Z$ with

[TABLE] 2. (ii)

we have the uniform bound

[TABLE] 3. (iii)

for any $x_{0}\in U$ and all $(x,y)$ satisfying $\|x-x_{0}\|+\|y-y(x_{0})\|\leq\xi$, we have

[TABLE]

It then follows that there exist constants $a,d>0$ depending only on $c_{0},c_{1},c_{2}$ and $L_{1}$ so that whenever it holds that

[TABLE]

then there exists a unique function $g$ defined from $\bigcup_{x_{0}\in U}B(x_{0},a)$ (union of balls centered at $x_{0}$ of radius $a$) to $Y$ such that

[TABLE]

Moreover, $g$ is a $C^{k}$ function on its domain of definition, and for all $x_{0}\in U$ and all $x\in B(x_{0},a)$

[TABLE]

where $K_{0}>0$ depends only on the constants $c_{0}$ and $c_{1}$.

For future reference we note that a suitable choice $a$ is given by

[TABLE]

where $\hat{a}=\frac{b}{4M}$, $b$ is such that $bL_{1}(b)\leq\frac{1}{2M}$, $M=\max\{c_{0},1+c_{0}c_{1}\}$ and $\epsilon>0$ is sufficiently small to ensure that $a$ is positive.

4.2. Proofs about the model

We begin with a technical lemma that is required to prove Theorem 2.2.

Lemma 4.3**.**

If $v\in\dot{\mathcal{H}}^{1}$, then, for any $l\in\Gamma_{\pm}$,

[TABLE]

Proof.

The argument in [OS12a, Proposition 12(ii)] proves the result for the case without a crack present. In that setting, the proof follows directly from [OS12b, Theorem 2.2]. In a crack geometry the constructions has to avoid crossing the crack $\Gamma_{0}$, hence we modify the argument. We distinguish two cases, depending on whether $l\in\Gamma_{+}$ or $l\in\Gamma_{-}$. We recall that $\hat{x}=\left(\tfrac{1}{2},\tfrac{1}{2}\right)$ and that by definition $v\in\dot{\mathcal{H}}^{1}\implies v(\hat{x})=0$.

Case 1: Let $l\in\Gamma_{+}$, which implies that $(l-\hat{x})\cdot e_{1}=0$. We consider a sequence of squares $(Q_{i})_{i=0}^{N}\subset\mathbb{R}^{2}\setminus\left(\mathcal{D}_{\Gamma}\cup\Gamma_{0}\right)$, where $\mathcal{D}_{\Gamma}$ denotes the continuum region enclosed by $\Gamma_{+}$ and $\Gamma_{-}$. The squares are aligned in the direction $e_{2}$, which is possible due to the assumption on $l$, and defined as follows. $Q_{0}$ and $Q_{N}$ are unit squares corresponding to sites $\hat{x}$ and $l$, defined in such a way that $\hat{x}$ (respectively $l$) is the midpoint of the side of $Q_{0}$ (resp. $Q_{N}$) which borders $\Gamma_{+}$. The squares $Q_{1},\dots,Q_{N-1}$ are defined to fill the space between $\hat{x}$ and $l$ in such a way that they have disjoint interiors and are such that their side-lengths differ by at most a factor of 2, with one side of the smaller square contained in one side of the larger square. It is easy to see that there is at most

[TABLE]

squares in the sequence. See Figure 7.

For any two neighbouring squares $Q_{j},\,Q_{j+1}$ it follows from a special case of [JN61, Lemma 2] that

[TABLE]

where

[TABLE]

and $I$ denotes the crack domain P1 interpolation operator employed in [BHO19, Theorem 2.3]. The final inequality in (4.6) follows from the fact that both components of $\nabla Iv$ are piecewise constant and each corresponds to $D_{\rho}v(l)$ for some bond $(m,m+\rho)$.

As a result

[TABLE]

Furthermore, it is naturally true that

[TABLE]

and since on each $Q_{i}$, the piecewise linear interpolant $v$ belongs to a finite-dimensional space, we obtain

[TABLE]

where the first inequality follows from the equivalence of norms for finite-dimensional spaces and the second from extending the domain from $Q_{i}$ to the whole of $\mathbb{R}^{2}\setminus\Gamma_{0}$.

With (4.7) and (4.8)-(4.9) in hand, we obtain

[TABLE]

which concludes the proof for $l\in\Gamma_{+}$.

Case 2: Let $l\in\Gamma_{-}$. The fact that $l$ is on the other side of the crack relative to $\hat{x}$ deems the previous argument invalid, as we can no longer define the sequence of squares aligned with $\hat{x}$ and $l$ which will be a subset of $\mathbb{R}^{2}\setminus\left(\mathcal{D}_{\Gamma}\cup\Gamma_{0}\right)$. Thus we first ’jump’ to the other side, that is we $\hat{a}=\left(\tfrac{1}{2},-\tfrac{1}{2}\right)$ and conclude that

[TABLE]

where the second inequality follows from applying (4.10) to a sequence of squares between $l$ and $\hat{a}$ (in fact there are only two such squares) and the fact that a bound on $|v(\hat{a})-v(\hat{x})|$ can be incorporated into the general form . ∎

We now show that the model is well-defined, with the particular emphasis on two new elements of the analysis that are distinct from previous arguments of this kind, e.g. [EOS16, BHO19].

Proof of Theorem 2.2.

We can decompose the energy into a bulk part and a crack surface part by writing

[TABLE]

where

[TABLE]

[TABLE]

We notice that $\mathcal{E}_{\rm bulk}$ excludes the bonds across the crack and thus is well-defined on $\dot{\mathcal{H}}^{1}\times\mathbb{R}$, as shown in the first part of Theorem 2.3 in [BHO19].

To establish the same for the $\mathcal{E}_{\Gamma}$, we note that we have the symmetry

[TABLE]

using this observation, for $m\in\Gamma_{\pm}$ and $\rho\in\mathcal{R}\setminus\tilde{\mathcal{R}}(m)$ we have

[TABLE]

Furthermore, it follows from Lemma 4.3 that $u\in\dot{\mathcal{H}}^{1}$ implies that $|Du(m)|\lesssim\log|m|$, which in particular implies that

[TABLE]

for suitable constants $C_{0},C_{1}$. Therefore, using assumption (2.12), it follows that for any $m\in\Gamma_{\pm}$ with $|m|$ sufficiently large,

[TABLE]

this entails that for each $u\in\dot{\mathcal{H}}^{1}$ we effectively only sum over a finite domain, and implies that $\mathcal{E}_{\Gamma}$ is indeed well-defined over $\dot{\mathcal{H}}^{1}\times\mathbb{R}$.

The differentiability properties of the functional follow from a standard argument, see [OT13]; here we simply provide formulae for derivatives of relevance to our subsequent arguments. In particular, we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

As stated, the foregoing expressions are valid for $v\in\mathcal{H}^{\rm c}$. To define them for $v\in\dot{\mathcal{H}}^{1}$ one requires an extension argument relying on showing that $\delta_{u}\mathcal{E}(0)\in(\dot{\mathcal{H}}^{1})^{*}$, which is proven in [BHO19, Theorem 2.3], and analogous results for the remaining terms.

We now turn to the analysis of the bifurcation path $\mathcal{B}$.

Proof of Proposition 2.4.

This is a standard result and follows from the fact that $\mathcal{E}$ is a $C^{\alpha}$ functional, so we may apply the local uniqueness of the function $g$ whose existence was asserted in Theorem 4.2. We therefore only outline the proof. Define sets corresponding to neighbourhoods of fold points

[TABLE]

Since ${\rm Im}(\mathcal{B})$ is compact, then so is ${\rm Im}(\mathcal{B})\setminus\mathcal{B}_{\rm f}$, thus the latter can be covered with a finite collection of neighbourhoods of points $\{(\bar{u}_{s_{i}},\bar{k}_{s_{i}})\}_{i=1,\dots,N}$. It will be shown in the proof of Theorem 2.8 that $H_{s_{i}}=\delta^{2}_{uu}\mathcal{E}(\bar{u}_{s_{i}},\bar{k}_{s_{i}})$ at each such point is an isomorphism, thus rendering Theorem 4.2 applicable to $\delta_{u}\mathcal{E}(\bar{u}_{s_{i}},\bar{k}_{s_{i}})$, giving us a locally unique $C^{\alpha-1}$ graph of critical points $k\mapsto u(k)$, which by its uniqueness together with injectivity of $\mathcal{B}$ has to coincide with $(\bar{u}_{s},\bar{k}_{s})$, thus ${\rm Im}(\mathcal{B})\setminus\mathcal{B}_{\rm f}$ is a piecewise $C^{\alpha-1}$ manifold. To establish the same in $\mathcal{B}_{\rm f}$, for each fold point $b_{i}$, one considers an extended system $\tilde{F}\,:\,(\dot{\mathcal{H}}^{1}\times\mathbb{R})\times\mathbb{R}$ given by \tilde{F}(u,k,t)=\big{(}\delta_{u}\mathcal{E}(u,k),(u-\bar{u}_{b_{i}},\gamma_{b_{i}})_{\dot{\mathcal{H}}^{1}}-t\big{)}, where $\gamma_{b_{i}}$ was introduced in Definition 2.3. The ABCD Lemma is applicable to this extended system evaluated at $(\bar{u}_{b_{i}},\bar{k}_{b_{i}},b_{i})$, thus ensuring that Theorem 4.2 is also applicable, giving us a locally unique $C^{\alpha-1}$ graph $t\mapsto(u(t),k(t))$, where in particular $k(0)=\bar{k}_{b_{i}}$. Again, due to uniqueness this must coincide with $(\bar{u}_{s},\bar{k}_{s})$, and hence this finishes the argument. ∎

Likewise, the existence of an eigen-pair can be established.

Proof of Proposition 2.5.

In what follows we always consider a system $G\,:\,B\times Y\to Z$ where $B\subset[0,1]$, $Y=\dot{\mathcal{H}}^{1}\times\mathbb{R}$ and $Z=(\dot{\mathcal{H}}^{1})^{*}\times\mathbb{R}$, given by

[TABLE]

where $c\in\dot{\mathcal{H}}^{1}$ will be chosen appropriately. We consider two subsets of $\mathcal{B}_{pt}$ seperately.

Throughout this proof we endow the product spaces with their canonical norms, for example, $\|(u,k)\|_{\dot{\mathcal{H}}^{1}\times\mathbb{R}}=\|u\|_{\dot{\mathcal{H}}^{1}}+|k|$.

(a) Vicinity of a bifurcation point: We let $c=\gamma_{b_{i}}$ and, in order to simplify the notation for derivatives of $G$, we introduce $y:=(\gamma,\mu)$ and observe that

[TABLE]

Thus, Lemma 4.1 and Theorem 4.2 together imply that, for $s\in(b_{i}-\xi,b_{i}+\xi)$, where $\xi>0$ is small enough (cf. (2.20)), there exists an eigen-pair $(\mu_{s},\gamma_{s})$. To show that $\mu_{s}^{\prime}:=\frac{{\rm d}\mu_{s}}{{\rm d}s}\neq 0$ at $s=b_{i}$, we differentiate both sides of (2.23) with respect to $s$ to obtain

[TABLE]

By definition, at a fold point we have $\bar{k}_{b_{i}}^{\prime}=0$ and since we can differentiate $\delta_{u}\mathcal{E}(\bar{u}_{s},\bar{k}_{s})=0$ with respect to $s$ to get $H_{s}\bar{u}_{s}^{\prime}+\bar{k}_{s}^{\prime}b_{s}=0$, it follows that $\bar{u}_{b_{i}}^{\prime}=\alpha\gamma_{b_{i}}$ for some $\alpha\neq 0$ (constant speed of parametrisation). Testing (4.15) at $s=b_{i}$ with $\gamma_{b_{i}}$ and simplifying, we obtain

[TABLE]

which is nonzero by Assumption (2.17). This completes case (a).

*(b) Unstable segment away from bifurcations: * We assume without loss of generality that at the bifurcation point $b_{i}$ we switch from a stable segment to an unstable segment. The result in (a) establishes existence of an eigenvector for $s\in(b_{i}-\xi,b_{i}+\xi)$, so we let $t_{1}:=b_{i}+\xi-\epsilon$, where $0<\epsilon<\xi$ and thus are able to set $c=\gamma_{t_{1}}$ in (4.14). A subsequent application of Theorem 4.2 to system $G$ with this newly chosen $c$ yields existence of a new interval $(t_{1}-\xi_{1},t_{1}+\xi_{1})\subset\mathcal{I}_{\rm pt}$ for which the premise of the theorem is true. This procedure can be iterated, for example by incrementing $t_{2}:=t_{1}+\xi_{1}/2$ and repeating the argument. To cover the entire unstable segment in this way we need to bound $\xi_{j}$ from below, independently of $t_{j}$.

Due to (4.16) we know that $\mu_{t}<0$, which implies that the subspace $U_{t}$ from Assumption 2 can be characterised as

[TABLE]

With this in hand we consider any $(u,k)\in\dot{\mathcal{H}}^{1}\times\mathbb{R}$, decompose $u$ as $u=\alpha\gamma_{t}+v$, where $\alpha\in\mathbb{R}$ and $v\in U_{t}$, and aim to uniformly bound

[TABLE]

from below.

To do so, we observe that

[TABLE]

Together, (4.18), (4.19) imply that

[TABLE]

Moreover, we trivially have

[TABLE]

Let $\tilde{c}_{0}(s)^{-1}=\min\{\frac{1}{2}(c-\mu_{t}),\frac{1}{2}\}>\min\{\frac{1}{2}(c),\frac{1}{2}\}=:c_{0}^{-1}$, then combining (4.20) and (4.21) yields

[TABLE]

In a similar vein, we observe that

[TABLE]

where $c_{1}:=\max_{t\in\mathcal{I}_{pt}}\tilde{c}_{1}(t)$, which is guaranteed to exist due to $\mathcal{E}$ being a $C^{\alpha}$ functional and that $\mathcal{B}$ is a $C^{\alpha-1}$ function of $s$, where $\alpha\geq 5$ by assumption.

It is also evident that Condition (iii) from Theorem 4.2 is satisfied with

[TABLE]

where $S(\tilde{s},\gamma_{\tilde{s}},\mu_{\tilde{s}},\xi)=\{(s,\gamma,\mu)\in\mathbb{R}\times\dot{\mathcal{H}}^{1}\times\mathbb{R}\,\colon\,|s-\tilde{s}|+\|\gamma-\gamma_{\tilde{s}}\|+|\mu-\mu_{\tilde{s}}|\leq\xi\}$.

Guided by (4.4), we define $M:=\max\{c_{0},1+c_{0}c_{1}\}$, choose $b$ such that $bL_{1}(b)\geq\frac{1}{2M}$, let $\hat{a}=\frac{b}{4M}$ and recall from (4.4) that $\xi=\min\left\{\frac{\hat{a}}{2},\frac{b}{2c_{2}}\right\}-\epsilon$ with sufficiently small $\epsilon>0$ is an admissible choice for $\xi$ which is in particular independent of $t$.

This completes the proof of part (b).

*(c) Regularity: * It remains to establish the $C^{\alpha-2}$-regularity of $s\mapsto(\mu_{s},\gamma_{s}),s\in\mathcal{I}_{\rm pt}$. To that end, we note that $G$ is a smooth function of $\gamma$ and $\mu$ and $C^{\alpha-2}$ function with respect to $s$, thus uniqueness and regularity parts of Theorem 4.2 immediately imply that $\gamma$ and $\mu$ are $C^{\alpha-2}(\mathcal{I}_{\rm pt})$ functions. ∎

We are now in a position to prove the spatial regularity of $\bar{u}_{s}$ and $\gamma_{s}$.

Proof of Theorem 2.6.

We begin by defining

[TABLE]

where $\tau\in\tilde{\mathcal{R}}(l)$ and $\mathcal{G}$ is the lattice Green’s function for the anti-plane crack geometry, as introduced in [BHO19, Theorem 2.6] and proven to satisfy decay property

[TABLE]

where $\omega$ is the complex square root mapping defined in polar coordinates as

[TABLE]

and $\delta>0$ is arbitrarily small. Here, and throughout this proof, $\lesssim$ should be read as $\leq C_{\delta}$ where $C_{\delta}$ is a constant that may depend on $\delta$.

*Proof of Theorem 2.6: estimate (2.24): * We first prove the decay estimate for $\bar{u}_{s}$. We can write

[TABLE]

The term (4.23) can be estimated by $|l|^{-3/2+\delta}$ due to the argument given in [BHO19, Theorem 2.4]; that is,

[TABLE]

The additional term (4.24) appears because we define the energy with the homogeneous discrete gradient operator and concerns the bonds crossing the crack surface. The estimate in (4.22) thus does not apply, but this term can be estimated as follows. Using (4.11) we see that we only sum over at most $2R_{\phi}^{3}$ lattice sites in (4.24) and thus we can decompose $D_{\rho}v(m)$ into a sum of finite differences along bonds that go around the crack. There will be at most $2R_{\phi}^{3}$ many of them and each separately decays like $|D_{1\rho}D_{2\tau}\mathcal{G}(m,l)|\lesssim|l|^{-3/2+\delta}$, thus ensuring that (4.24) can be bounded by

[TABLE]

where $C<4R_{\phi}^{6}$ (though an optimal $C$ is much smaller). This concludes the proof of (2.24).

*Proof of Theorem 2.6: estimate (2.25): * To estimate $\gamma_{s}$ we employ an analogous argument. We begin by writing

[TABLE]

Using precisely the same argument as for (4.24) we can bound (4.26) by $|l|^{-3/2+\delta}$.

To estimate (4.25) we Taylor-expand $\phi^{\prime\prime}$ around 0 and, noting that we assume that $\phi^{\prime\prime}(0)=1$, we observe that

[TABLE]

where $|R(m)|\lesssim|m|^{-1}$ is the remainder of the expansion. This readily implies that $|D_{\tau}\gamma_{s}(l)|\lesssim|l|^{-1}$. Thus looking again at (4.27), instead of applying Cauchy-Schwarz inequality, we directly observe that

[TABLE]

since $|R(m)|{|}\tilde{D}\gamma_{s}(m)|\lesssim|m|^{-2}$. As before, $\delta>0$ is arbitrarily small. This completes the proof of the second bound (2.25). ∎

Remark 4.4*.*

It is interesting to note that while the model includes a full interaction between nearest-neighbour atoms, even across the crack, it is nonetheless the lattice Green’s function for the fractured domain that is employed to estimate the atomistic solutions. The homogeneous lattice Green’s function fails because the finite differences of $\hat{u}_{k}(m)$ across the crack grow like $\sim|m|^{1/2}$.

4.3. Convergence Proofs

In tandem with the results from bifurcation theory stated in Section 4.1, in order to prove the results from Section 2.4 we rely on the following auxiliary result from [EOS16] that was adapted to domain with cracks in [BHO19].

Lemma 4.5**.**

There exists a truncation operator $T_{R}\,:\,\dot{\mathcal{H}}^{1}\to\mathcal{H}_{R}^{0}$ such that $T_{R}v=0$ in $\Lambda\setminus B_{R}$ and which satisfies

[TABLE]

We can now prove the main result of this section.

Proof of Theorem 2.8.

We consider an extended system $F\,:\,B\times Y\to Z$ where $B=[0,1]$, $Y=\mathcal{H}^{0}_{R}\times\mathbb{R}$ and $Z=(\mathcal{H}^{0}_{R})^{*}\times\mathbb{R}$ given by

[TABLE]

where $\bar{u}_{s}^{\prime}=\frac{d\bar{u}(s)}{ds}$ and $y=(u_{y},k_{y})$. We further introduce a mapping $y_{R}\,:\,B\to Y$ given by $y_{R}(s)=(T_{R}\bar{u}_{s},\bar{k}_{s})$. We shall now show that, with the help of ABCD Lemma, $F$ satisfies the conditions of Theorem 4.2.

One can easily obtain that

[TABLE]

We further define $A_{s}:=H_{s}=\delta^{2}_{uu}\mathcal{E}(\bar{u}_{s},\bar{k}_{s})$ (renamed to keep intuitive notation) and $b_{s}:=\delta^{2}_{uk}\mathcal{E}(\bar{u}_{s},\bar{k}_{s})$ and observe that

[TABLE]

Here, we treat $A_{s}$ as a restriction to $\mathcal{H}_{0}^{R}$ and $b_{s}$ as an element of $(\mathcal{H}_{0}^{R})^{*}$. Since $T_{R}\bar{u}_{s}\to\bar{u}_{s}$ as $R\to\infty$ strongly in $\dot{\mathcal{H}}^{1}$ (a consequence of (4.28), the decay estimate from Theorem 2.6 and $\mathcal{E}\in C^{\alpha}$),

[TABLE]

as $R\to\infty$. Our strategy will therefore be to apply the ABCD Lemma to $M_{s}^{1}$, interpreted as an operator from $\dot{\mathcal{H}}^{1}\times\mathbb{R}$ to $(\dot{\mathcal{H}}^{1})^{*}\times\mathbb{R}$, and show that $M_{s}^{2}$ is a small perturbation.

To carry out this strategy we begin by differentiating $\delta_{u}\mathcal{E}(\bar{u}_{s},\bar{k}_{s})=0$ with respect to $s$ to obtain

[TABLE]

along the bifurcation path $\mathcal{B}$. At a fold point, when $s=b_{i}$, due to (2.16), we have $\bar{k}_{s}^{\prime}=0$, thus revealing that $\bar{u}_{b_{i}}^{\prime}=\alpha\gamma_{b_{i}}$ for some non-zero $\alpha\in\mathbb{R}$. For $s\neq b_{i}$, the operator $A_{s}$ is invertible and thus

[TABLE]

We can now show that $M_{s}^{1}$ satisfies the conditions of ABCD Lemma.

*Case 1, $s\in\mathcal{I}_{\rm pos}$: * Suppose that $s\in\mathcal{I}_{\rm pos}$ from (2.20). In this case $A_{s}$ is an isomorphism due to (2.19) and (2.22). Thus, to apply ABCD Lemma to $M_{s}^{1}$, we have to check that $(\bar{u}^{\prime}_{s},(A_{s})^{-1}b_{s})_{\dot{\mathcal{H}}^{1}}\neq 0$, which is true since

[TABLE]

since by definition at a regular point we have $\bar{k}_{s}^{\prime}\neq 0$ and $\bar{u}_{s}^{\prime}\neq 0$.

*Case 2, $s\in\mathcal{I}_{\rm pt}$: * Now suppose $s\in\mathcal{I}_{\rm pt}$ but $s\neq b_{i}\;\forall i\in\{1,\dots,M\}$. It can be shown that $A_{s}$ remains an isomorphism as follows. Proposition 2.5 tells us that we have an eigen-pair $(\mu_{s},\gamma_{s})$ satisfying (2.23). Any $v\in\dot{\mathcal{H}}^{1}$ can be decomposed into $v=\alpha\gamma_{s}+w$, where $w\in U_{s}$ with $U_{s}$ given by (4.17) and $\alpha\in\mathbb{R}$. Thus,

[TABLE]

and we can further estimate

[TABLE]

where $c$ is the stability constant from Assumption 2. This, together with the fact that $\|v\|\leq|\alpha|+\|w\|$ readily implies that we can set $\tilde{c}:=\frac{1}{2}\min\{|\mu_{s}|,c\}$ and conclude that for all $v\in\dot{\mathcal{H}}^{1}$

[TABLE]

Thus, as in the case $s\in\mathcal{I}_{\rm pos}$, (4.32) ensures that we can apply the ABCD lemma and deduce again that $A_{s}$ is an isomorphism.

*Case 3, $s=b_{i}$: * Finally, suppose $s=b_{i}$ for some $i\in\{1,\dots,M\}$. Due to Assumption 2 we know that the kernel of $A_{s}$ is one-dimensional at a fold point and thanks to Proposition 2.5 we know that it is spanned by $\gamma_{s}$, which means that (4.31) implies that $\bar{u}_{s}^{\prime}=\gamma_{s}$. By Definition 2.3 we know that $\langle b_{s},\,\gamma_{s}\rangle\neq 0$, which implies that the ABCD Lemma is again applicable.

*Uniform Stability of $M_{s}^{1}$: * We have shown so far that, for all $s\in[0,1]$, $M_{s}^{1}$ is an isomorphism from $\dot{\mathcal{H}}^{1}\times\mathbb{R}$ to $(\dot{\mathcal{H}}^{1})^{*}\times\mathbb{R}$. In particular, this implies that for any $x=(u_{x},k_{x})\in\dot{\mathcal{H}}^{1}\times\mathbb{R}$ we have

[TABLE]

where $\tilde{c}_{s}>0$.

Since $s\mapsto M_{s}^{1}$ is continuous in operator-norm due to smoothness of $\mathcal{E}$ established in Theorem 2.2 and smoothness of $s\mapsto(\bar{u}_{s},\bar{k}_{s})$ established in Proposition 2.4, it follows that the infimum $\inf\tilde{c}_{s}$ is attained on $[0,1]$ and must therefore be positive. In summary, we have established the existence of $\tilde{c}>0$ such that

[TABLE]

*Uniform Stability: * Next, using the definition of $M_{s}^{2}$ we can bound

[TABLE]

for $R$ large enough, thus ensuring that $D_{y}F(s,y_{R}(s))$ is an isomorphism from $\mathcal{H}_{0}^{R}\times\mathbb{R}$ to $(\mathcal{H}_{0}^{R})^{*}\times\mathbb{R}$, thus satisfying condition (i) from Theorem 4.2, with uniform bound

[TABLE]

that is $c_{0}$ from (4.1) is given by $c_{0}=\frac{2}{\tilde{c}}$.

*Conclusion: * So far we have confirmed Condition (i) of Theorem 4.2. To conclude the proof, we now need to also check conditions (ii, iii).

It can be readily checked that $\|D_{s}F(s,y_{R}(s))\|=|-1+(T_{R}\bar{u}_{s}-\bar{u}_{s},\bar{u}_{s}^{\prime\prime})_{\dot{\mathcal{H}}^{1}}|\leq 2$ for $R$ large enough. Thus the condition (ii) in Theorem 4.2 is satisfied with $c_{1}$ in (4.2) given by $c_{1}=2$.

The condition (iii) from Theorem 4.2 is satisfied with

[TABLE]

where $S(s,y_{R}(s),\xi)=\{(s_{0},y_{0})\in\mathbb{R}\times(\mathcal{H}_{0}^{R}\times\mathbb{R})\,\colon\,|s-s_{0}|+\|T_{R}\bar{u}_{s}-u_{y}\|+|\bar{k}_{s}-k_{y}|\leq\xi\}$.

Finally, we observe that

[TABLE]

as $R\to\infty$, which implies that no matter how large constants $c_{0},c_{1}$ and $c_{2}$ were and how badly behaved $L_{1}$ was, we would still fall within the regime where the result of Theorem 4.2 was applicable for $R$ large enough.

We can thus conclude that there exists $\mathcal{B}_{R}\,:\,[0,1]\to\mathcal{H}_{0}^{R}\times\mathbb{R}$ given by $\mathcal{B}_{R}(s):=(\bar{u}^{R}_{s},\bar{k}_{s}^{R})$, such that $F(s,(\bar{u}^{R}_{s},\bar{k}_{s}^{R}))=0$, which in particular implies $\delta_{u}\mathcal{E}(\bar{u}^{R}_{s},\bar{k}_{s}^{R})=0$. Furthermore, using (4.3) we can conclude that

[TABLE]

for arbitrarily small $\beta>0$. Crucially, $K_{0}$ depends only on $c_{0}$ and $c_{1}$, which are independent of $s$ and the last inequality follows from Lemma 4.28 and the regularity estimate from Theorem 2.6. This concludes the result, since trivially $\|\bar{u}_{s}^{R}-\bar{u}_{s}\|_{\dot{\mathcal{H}}^{1}}\leq\|\bar{u}_{s}^{R}-T_{R}\bar{u}_{s}\|_{\dot{\mathcal{H}}^{1}}+\|T_{R}\bar{u}_{s}^{R}-\bar{u}_{s}\|_{\dot{\mathcal{H}}^{1}}$.

Finally, we note that (2.28) follows as an immediate collorary, arguing exactly as in the proof of [EOS, Theorem 2.4]. ∎

To prove our final result, the superconvergence of critical values of SIF, we first need to quote two intermediate technical steps. The first lemma, which highlights the origin of this superconvergence, is taken from [CST00, Theorem 4.1], restated in our notation for the sake of convenience.

Lemma 4.6**.**

Let $(\bar{u}_{b_{i}},\bar{k}_{b_{i}})\in\mathcal{B}$ be a simple quadratic fold point. Under Assumptions 1,2 & 3, for $R$ large enough, the approximate bifurcation diagram $\mathcal{B}_{R}$ has a quadratic fold point at $s=b_{i}^{R}$, where $|b_{i}^{R}-b_{i}|\to 0$ as $R\to\infty$. Furthermore,

[TABLE]

To exploit the inequality from Lemma 4.6, we adapt [BRR80, Theorem 2], which is a follow-up result to Theorem 4.2 for derivatives.

Lemma 4.7**.**

Assume the hypotheses of Theorem 4.2 and in addition that

[TABLE]

Then there exists a continuous function $K\,:\,\mathbb{R}_{+}\to\mathbb{R}_{+},$ which depends only on $c_{0},c_{1},L_{1}$ such that for all $x_{0}\in U$ and all $x\in B(x_{0},a)$ it holds that

[TABLE]

Proof of Theorem 2.9.

Lemma 4.6 implies that, for domain radius $R$ large enough, we have exactly $M$ approximate fold points $b_{j}^{R}\to b_{j}$ as $R\to\infty$. By assumption, $b_{j}\in(0,1)$ and hence also $b_{j}^{R}\in(0,1)$. Arguing analogously as in the proof of Theorem 2.8, it is not difficult to show that $F$ defined in (4.29) satisfies the conditions of Lemma 4.7, thus

[TABLE]

for arbitrary small $\beta>0$. Furthermore,

[TABLE]

with the first inequality following from the obvious fact that $T_{R}\gamma_{b_{i}}\in{\mathcal{H}^{0}_{R}}$ and the second from the regularity result for $\gamma_{b_{i}}$ in Theorem 2.6.

Applying Lemma 4.6 we therefore obtain the desired result that

[TABLE]