An iterative scheme for the generalized Peierls-Nabarro model based on the inverse Hilbert transform

TL;DR

This paper introduces a semi-analytical iterative scheme leveraging the inverse Hilbert transform to efficiently solve the generalized Peierls-Nabarro model, enabling detailed dislocation core structure analysis and external stress incorporation.

Contribution

A novel semi-analytical iterative method based on the inverse Hilbert transform for solving the generalized Peierls-Nabarro model is proposed, reducing complex equations to a local fixed point iteration.

Findings

Validated with 1D Peierls-Nabarro model examples

Applied to dislocation core structures in Aluminium

Discussed external stress incorporation and method limitations

Abstract

A new semi-analytical iterative scheme is proposed in this work for solving the generalized Peierls-Nabarro model. The numerical method developed here exploits certain basic properties of the Hilbert transform to achieve the desired reduction of the non-local and non-linear equations characterizing the generalized Peierls-Nabarro model to a local fixed point iteration scheme. The method is validated with simple examples involving the 1D Peierls-Nabarro model corresponding to a sinusoidal stacking fault energy, and with calculations of the core structure of both edge and screw dislocations on the close-packed $\{111\}$ planes in Aluminium. An approximate technique to incorporate external stresses within the framework of the proposed iterative scheme is also discussed with applications to the equilibration of a dislocation dipole. Finally, the advantages, limitations and avenues for future extension of the proposed method are discussed.

Equations

$$E(s_1,s_2)=\sum _{\alpha ,\beta }\frac{K_{\alpha \beta }}{4\pi }{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\frac{{\rho }_{\alpha }(x^{\prime })s_{\beta }(x)}{x-x^{\prime }}dxd{x}^{\prime }+{\int }_{-\infty }^{\infty }\gamma (s_1(x),s_2(x))dx.$$

$${\rho }_{\alpha }(x)=\frac{d{s}_{\alpha }(x)}{dx}$$

$$\sum _{\beta }\frac{K_{\alpha \beta }}{2\pi }{\int }_{-\infty }^{\infty }\frac{{\rho }_{\beta }(x^{\prime })}{x-x^{\prime }}d{x}^{\prime }={\tau }_{\alpha }(s_1(x),s_2(x)),$$

$${\tau }_{\alpha }({\delta }_1,{\delta }_2)=-{\partial }_{\alpha }\gamma ({\delta }_1,{\delta }_2)$$

$$\gamma (\delta )=\frac{\mu b^2}{4{\pi }^{2}d}\left(1-\cos \frac{2\pi \delta }{b}\right),$$

$${\int }_{-\infty }^{\infty }\frac{\rho (\xi )}{x-\xi }d\xi =-\frac{1-\nu }{d}\sin \frac{2\pi s(x)}{b},$$

$$s(x)=\frac{b}{2}+\frac{b}{\pi }{\tan }^{-1}\frac{x}{\xi },\xi =\frac{d}{2(1-\nu )}.$$

$$\mathcal{H}\left(\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\rho }_{\beta }(x)\right)={\tau }_{\alpha }(s_1(x),s_2(x)).$$

$${\tau }_{\alpha }(s_1(x),s_2(x))=\mathcal{H}(g_{\alpha }(x)),$$

$$\mathcal{H}\left(\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\rho }_{\beta }(x)\right)=\mathcal{H}(g_{\alpha }(x))\Rightarrow \sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\rho }_{\beta }(x)=g_{\alpha }(x)+l_{\alpha },$$

$$\gamma ({\delta }_1,{\delta }_2)=\sum _{m_1=-\infty }^{\infty }\sum _{m_2=-\infty }^{\infty }{c}_{m_1,m_2}\exp \left(i\sum _{\beta }{m}_{\beta }{k}_{\beta }{\delta }_{\beta }\right),$$

$$k_{\alpha }=\frac{2\pi }{b_{\alpha }},$$

$${\tau }_{\alpha }({\delta }_1,{\delta }_2)=-\sum _{m_1,m_2}i{m}_{\alpha }{k}_{\alpha }{c}_{m_1,m_2}\exp \left(i\sum _{\beta }{m}_{\beta }{k}_{\beta }{\delta }_{\beta }\right).$$

$$\begin{array}{rl}{\tau }_{\alpha }(s_1(x),s_2(x))& =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{\hat{d}}_{\alpha }(k)\exp (ikx)dk,\\ {\hat{d}}_{\alpha }(k)& ={\int }_{-\infty }^{\infty }{\tau }_{\alpha }(s_1(x),s_2(x))\exp (-ikx)dx.\end{array}$$

$${\tau }_{\alpha }(s_1(x),s_2(x))=\mathcal{H}\left(\frac{i}{2\pi }{\int }_{-\infty }^{\infty }\text{sgn}(k){\hat{d}}_{\alpha }(k)\exp (ikx)dk\right).$$

$$\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\rho }_{\beta }(x)=\frac{i}{2\pi }{\int }_{-\infty }^{\infty }\text{sgn}(k){\hat{d}}_{\alpha }(k)\exp (ikx)dk+l_{\alpha },$$

$${\tau }_{\alpha }(s_1(x),s_2(x))\simeq \sum _{n=-\infty }^{\infty }{d}_{\alpha ,n}\exp in{k}_{0}x,$$

$$k_0=\frac{\pi }{L},$$

$$d_{\alpha ,n}=\frac{1}{2L}{\hat{d}}_{\alpha }(n{k}_0).$$

$$s_{\alpha }(x)=s_{\alpha ,j}+\frac{s_{\alpha ,j+1}-s_{\alpha ,j}}{h_j}(x-x_j),x\in [x_j,x_{j+1}].$$

$$\begin{array}{rl}{\hat{d}}_{\alpha }(k)& =-\frac{i}{2L}\sum _{j=0}^{N_s-1}\sum _{m_1,m_2}{m}_{\alpha }{k}_{\alpha }{c}_{m_1,m_2}\exp \left(i\sum _{\beta }{m}_{\beta }{k}_{\beta }{\lambda }_{\beta ,j}\right)I_{j,m_1,m_2}^n,\\ I_{j,m_1,m_2}^n& =\{\begin{array}{ll}{h}_j,& \text{ if }{L}_{j,m_1,m_2}^n=0,\\ \frac{\exp \left(i{L}_{j,m_1,m_2}^n{x}_{j+1}\right)-\exp \left(i{L}_{j,m_1,m_2}^n{x}_j\right)}{L_{j,m_1,m_2}^n},& \text{otherwise,}\end{array}\\ L_{j,m_1,m_2}^n& =\sum _{\beta }{m}_{\beta }{k}_{\beta }{\mu }_{\beta ,j}-n{k}_0\\ {\lambda }_{\alpha ,j}& =\frac{x_{j+1}{s}_{\alpha ,j}-x_j{s}_{\alpha ,j+1}}{h_j},{\mu }_{\alpha ,j}=\frac{s_{\alpha ,j+1}-s_{\alpha ,j}}{h_j}.\end{array}$$

$$\mathcal{H}\left(\frac{1}{2}\sum _{\beta }{K}_{\alpha \beta }{\rho }_{\beta }(x)\right)=\sum _{n=-\infty }^{\infty }{d}_{\alpha ,n}\exp (in{k}_{0}x).$$

$$\mathcal{H}\left(\frac{1}{2}\sum _{\beta }{K}_{\alpha \beta }{\rho }_{\beta }(x)\right)=\mathcal{H}\left(\sum _{n}i\text{sgn}(n)d_{\alpha ,n}\exp (in{k}_{0}x)\right),$$

$$\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }\frac{d{s}_{\beta }(x)}{dx}=\sum _{n}i\text{sgn}(n)d_{\alpha ,n}\exp (in{k}_{0}x)+a_{\alpha ,1},$$

$$\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }\frac{d{\tilde{s}}_{\beta }^{(k+1)}(x)}{dx}=\sum _{n}i\text{sgn}(n)d_{\alpha ,n}(\{s_{\beta ,j}^{(k)}\})\exp (in{k}_{0}x)+a_{\alpha ,1}.$$

$$\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\tilde{s}}_{\beta }^{(k+1)}(x)=\sum _n\frac{\text{sgn}(n)}{n{k}_0}{d}_{\alpha ,n}(\{s_{\beta ,j}^{(k)}\})\exp (in{k}_{0}x)+a_{\alpha ,1}x+a_{\alpha ,2},$$

$${\tilde{s}}_{\alpha }^{(k+1)}(-L)=0,{\tilde{s}}_{\alpha }^{(k+1)}(L)=S_{\alpha }.$$

$$\begin{array}{rl}{a}_{\alpha ,1}& =\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }\frac{S_{\beta }}{2L},\\ b_{\alpha ,2}& =\sum _{\beta }\frac{1}{4}{K}_{\alpha \beta }{S}_{\beta }-\sum _n\frac{(-1{)}^n\text{sgn}(n)}{n{k}_0}{d}_{\alpha ,n}(\{s_{\beta ,j}^{(k)}\}).\end{array}$$

$$\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }{\tilde{s}}_{\beta }^{(k+1)}(x)=\sum _n\frac{\text{sgn}(n)}{n{k}_0}\left(d_{\alpha ,n}(\{s_{\beta ,j}^{(k)}\})\exp (in{k}_{0}x)-(-1{)}^n\right)+\sum _{\beta }\frac{1}{2}{K}_{\alpha \beta }\frac{S_{\beta }}{2L}(x+L).$$

$$s_{\alpha ,j}^{(k+1)}=(1-\lambda )s_{\alpha ,j}^{(k)}+\lambda {\tilde{s}}_{\alpha ,j}^{(k+1)},$$

Full text

††articletype: PREPRINT

1 Introduction

It is well-known that the primary mediators for plasticity in crystalline materials are topological defects called dislocations. A proper understanding of the physics of dislocations, which is of crucial importance to many applications in materials science and engineering, entails considerations of processes from sub-atomic to continuum length and time scales. Notwithstanding the interplay of various factors across such an enormous range of scales, numerous simplified models that seek to explain specific aspects of dislocations have been proposed. One of the earliest, yet highly influential, models that throws light on the structure of dislocation cores is the Peierls-Nabarro (PN) model, originally proposed by Peierls [1], and later elaborated upon by Nabarro [2]. The current paper contributes a new iterative scheme to solve a particularly useful generalization of the original Peierls-Nabarro model of dislocations developed by Schoeck [3]. Specifically, the proposed method provides an efficient framework to model the slip distribution that characterizes the core structure of dislocations.

To set the context for the ensuing discussion, it is recalled that the energy of a crystalline solid containing a dislocation can be decomposed into two terms, one involving the elastic energy of the bulk, and another equal to the energy associated with the dislocation itself. Two common strategies are employed in the literature to effect this decomposition. The simpler of the two assigns an ad-hoc core radius to a dislocation and assumes linear elastic behavior outside the core region. The bulk energy associated with the region outside the defect core is computed using well established techniques in linear elasticity theory [4]. The energy of the core region, however, is not amenable to such a computation since the assumptions underlying linear elastic behavior are not valid in the dislocation core. A more serious limitation of the core radius approach is that the core structure of real dislocations is often complex enough (see, for instance, [5], [6]) that assigning a core radius fails to capture the essential features of the core structure. The PN model provides an alternative and more physically justified framework to study dislocation cores by taking advantage of the fact that dislocations are the boundaries between slipped and unslipped regions on well defined slip planes. The core structure is then determined by a competition between the elastic energy associated with the slip distribution and the energetic cost of the misfit between the slipped and unslipped regions. An attractive feature of the PN model is its variational structure, which permits a variety of systematic and controlled numerical approximations of the core structure [7], [8], [9], [10]. In spite of the various assumptions inherent in the PN framework, it has been remarkably successful in explaining the core structure of dislocations in a variety of materials [10], [11], [12], [13], [14], [15].

There are many known limitations of the Peierls-Nabarro model of dislocations, two of which are briefly discussed here. First, the use of the generalized stacking fault energy as the misfit energy, as originally proposed by Vitek [16], [17], assumes a priori that the slip distribution varies slowly enough for the misfit energy to accurately represent the misfit energy. A second limitation of the PN model that is important for this work is that all information about the discreteness of the lattice plane on which the slip plane resides is lost, despite the fact that the generalized stacking fault energy is periodic with the periodicity of the lattice on the slip plane. Thus, there is no energetic cost to moving a dislocation on the slip plane according to the PN model. This is in sharp contrast to real dislocations which experience a lattice resistance to their motion due to the periodic energy landscape on the slip plane, called the Peierls potential. The stress required to overcome the Peierls potential, called the Peierls stress, is thus not captured by the PN model. An approximate means to estimate the Peierls stress was introduced by Nabarro [2] by summing the total energy at discrete locations corresponding to the lattice points on the slip plane. This introduces a dependence of the total energy on the position of the dislocation line with respect to the atomic lattice, thereby recovering the lattice scale periodicity in the energy landscape that is absent in the original Peierls model. A succinct summary of other shortcomings of the PN model can be found in [18], [19].

A lot of effort has gone into remedying the various shortcomings of the Peierls-Nabarro model: see, for instance, [3], [11], [20], for extensions of various aspects of the PN model, [21] for a modern reformulation of the Peierls-Nabarro model as a phase-field model, and [22] for a discrete generalization of the PN model. Among these, the generalized Peierls-Nabarro model developed by Schoeck [3] is of note. The generalized PN model for an infinite dislocation assumes that the slip distribution has components both parallel and perpendicular to the dislocation line, and further uses anisotropic linear elasticity to model the elastic interactions. The generalized PN model has proven to be effective in providing insights into the core structure of dislocations in many materials; see, for instance, [23]. The iterative solution scheme developed in this work is aimed at solving Schoeck’s generalized PN model.

Among the many existing approaches to solve the generalized PN model, two classes of methods are widely used. The first class employs an ansatz for the slip distribution and/or the misfit stress distribution with the desired asymptotic properties and computes the slip distribution in a least squares sense. The ansatz for the slip distribution is chosen typically as a linear combination of the solution of the original PN equation. Instances of this approach can be found in [7], [8], [9], [24]. An attractive feature of this approach is that the various parameters that enter the numerical model have straightforward physical interpretation. The second class of solution techniques is the so-called semi-discrete variational PN model (SVPN) that is based on a piecewise linear finite element discretization of the total energy of a crystalline solid with a dislocation [10]. This is results in a nonlinear and nonlocal equation for the slip distribution which is solved using standard algorithms. It is to be noted that the process of solving these equations involves the use of dense matrices. Unlike the first class of methods mentioned earlier that assume a specific form of slip distribution, the approximation adopted in the SVPN method is ansatz-free. A variety of other methods have also been developed in the literature to solve the PN model: see [25] for a discrete Fourier transform based method to solve the generalized PN model for curved dislocations, and [26] for a solution scheme based on the fast multipole method. The present contribution presents a novel computational scheme to solve the generalized PN model, and may be viewed as a competitive alternative to the aforementioned schemes.

The outline of the paper is as follows. Details regarding the mathematical formulation of the new iterative scheme for solving the generalized PN model are presented first. This is followed by a presentation of various examples that illustrate the applicability of the proposed model. Finally, the advantages, limitations and future extensions of the iterative method are discussed before concluding the paper. Relevant background information is summarized in two appendices.

2 Mathematical formulation

The mathematical formulation of a new iterative scheme to solve the generalized Peierls-Nabarro model is presented in this section. Towards this end, the generalized PN model is reviewed first. This is followed by a reformulation of the same to obtain an inverse Peierls-Nabarro model by exploiting certain properties of the Hilbert transform; relevant properties of the Hilbert transform are summarize in Appendix A. The inverse PN model is subsequently approximated via a simple fixed point iteration scheme to obtain the desired iterative scheme that is the key contribution of this work. Unlike the generalized PN model, there is no straightforward means to incorporate the effect of an external stress within the framework of the inverse PN model; a possible strategy to include external stresses is finally discussed to complete the mathematical formulation of the newly proposed technique.

2.1 Recap of the generalized Peierls-Nabarro model

A summary of the generalized PN model introduced by Schoeck [3] is provided here, primarily for fixing the notation and terminology for the analysis to follow. The generalized PN model aims at describing the core structure of an infinite straight dislocation in an otherwise perfect crystal; a schematic illustration of an infinite (edge) dislocation is shown in Figure 1.

The line direction of the dislocation is assumed to be along the $x_{2}$ direction111The notations $x,y,z$ are used interchangeably with $x_{1},x_{2},x_{3}$., along which it is infinite in extent. The glide plane, also referred to as the slip plane, of the dislocation is the $x_{1}x_{2}$ plane $x_{3}=0$. The Burgers vector of the dislocation is assumed to lie in the slip plane. The crystal is assumed to be infinite along all three directions. The total energy of the crystalline solid with a single dislocation is modeled as the sum of two contributions: the elastic energy of the two halves of the material separated by the slip plane, and the energy associated with the specific slip distribution on the slip plane. Owing to the geometry chosen for this problem, the slip distribution is assumed to be uniform along the $x_{2}$ direction. Furthermore, the slip is assumed to be entirely contained in the $x_{1}x_{2}$ plane. With these assumptions in place, the slip distribution on the slip plane is of the form $(s_{1}(x),s_{2}(x))$. The total energy of the dislocated crystal with the slip distribution $(s_{1}(x),s_{2}(x))$ can thus be written as222The following notation is used: greek indices take the values $1$ and $2$. The Einstein convention regarding repeated indices is not followed; all summations are explicitly indicated, though the range of the variables is often to be inferred from the context.

[TABLE]

Here,

[TABLE]

denotes the component of the dislocation density along the $x_{\alpha}$ direction. The second term in (1) represents the misfit energy associated with a slip distribution $(s_{1}(x),s_{2}(x))$ over the slip plane. A convenient choice for the misfit energy is obtained by choosing $\gamma(\delta_{1},\delta_{2})$ to be the the generalized stacking fault energy obtained as the elastic energy due to a uniform relative slip $(\delta_{1},\delta_{2})$ of the upper ($z>0$) half of the crystal relative to the lower half [17]. $\{K_{\alpha\beta}/4\pi\}$ in (1) are the components of the symmetric anisotropic Stroh tensor [27] which depends on the crystalline lattice and the geometry of the slip plane. The equilibrium slip distribution is then obtained as the minimizer of the energy functional (1). This yields the generalized PN equation:

[TABLE]

where

[TABLE]

is the stress associated with a uniform slip $(\delta_{1},\delta_{2})$, called henceforth the misfit stress. The slip distribution corresponding to the presence of a single dislocation is obtained by solving the generalized PN equation, (3), subject to the boundary conditions $s_{1}(-\infty)=0,s_{1}(\infty)=b,s_{2}(\pm\infty)=0$.

The original PN model ([1], [2]) is obtained by considering an isotropic elastic medium with slip distribution constrained by the condition $s_{2}(x)\equiv 0$. Assuming a simple sinusoidal form of the stacking fault energy,

[TABLE]

the PN equation reduces to the following equation:

[TABLE]

where $s(x)=s_{1}(x)$, $\mu$ is the shear modulus and $\nu$ is the Poisson’s ratio of the isotropic elastic medium. This equation, with the boundary conditions stated earlier admits an exact solution:

[TABLE]

This solution will be used as a basic benchmark for the numerical methods presented later. The quantity $\xi$ serves as a useful characterization of the width of the dislocation.

2.2 Inverse Peierls-Nabarro model

The generalized PN equation (3) is now reformulated in a form that will prove to be useful in the sequel. Using the definition of the Hilbert transform (see Appendix A), the generalized PN equation (3) can be written as

[TABLE]

The key idea behind the numerical strategy proposed in this work is the observation that if the misfit stress distribution $\tau_{\alpha}(s_{1}(x),s_{2}(x))$ can be expressed as the Hilbert transform of some function $g_{\alpha}(x)$, as in

[TABLE]

then the generalized PN equation (3) can be replaced by an equation of the form

[TABLE]

where $l_{\alpha}$ is a constant. Such an equation is much simpler to handle numerically compared to the generalized PN model (3), as will become evident shortly. This is accomplished through a sequence of steps as follows. To begin with, the periodicity of the lattice implies that the generalized stacking fault energy $\gamma(\delta_{1},\delta_{2})$ is a periodic function of both its arguments, with periodicity $b_{1}$ and $b_{2}$ along the $x_{1}$ and $x_{2}$ directions, respectively. It is thus possible to express the generalized stacking fault energy in terms of a Fourier series of the form

[TABLE]

where

[TABLE]

and $\{c_{m_{1},m_{2}}\}$ are complex constants. It is noted that in practice only a finite number of such coefficients $c_{m_{1},m_{2}}$ are non-zero. A further reduction in the number of coefficients $c_{m_{1},m_{2}}$ that need to be computed is achieved by noting that $\gamma(\delta_{1},\delta_{2})$ is real. Using (4), the misfit stresses $\tau_{\alpha}(\delta_{1},\delta_{2})$ are obtained from the Fourier expansion of the generalized stacking fault energy (9) as

[TABLE]

Note that $\tau_{\alpha}(\delta_{1},\delta_{2})$ is also a periodic function with periodicities $b_{1}$ and $b_{2}$ along the $x_{1}$ and $x_{2}$ directions. However, given an arbitrary slip distribution $(s_{1}(x),s_{2}(x))$, the misfit stress distribution $\tau_{\alpha}(s_{1}(x),s_{2}(x))$ is in general not a periodic function of $x$. It is helpful to express the misfit stress distribution using Fourier transform as

[TABLE]

Using (40) and (38), it is easily seen that

[TABLE]

Since (13) expresses the misfit stress distribution as the Hilbert transform of some function, the arguments given in the beginning of this section can be used to invert the generalized Peierls-Nabarro model. Specifically, using (13) in (8) and using the inverse Hilbert transform yields the following equation:

[TABLE]

where $l_{\alpha}$ is a constant determined by the boundary conditions. Since the process of deriving (14) from (8) is effected with the aid of the inverse Hilbert transform, the foregoing equation (14) is henceforth referred to as the inverse Peierls-Nabarro model. It is emphasized that the inverse PN model (14) is essentially exact, and it is easy to check that solutions of the inverse PN model are identical to those of the generalized PN model (3).

The inverse PN model (14) does not admit an analytical solution in general. Nevertheless, an analytically tractable approximation of the inverse PN model (14) can be obtained as follows. Since the lattice structure is almost perfectly restored far away from the dislocation, it is clear that $s_{\alpha}(L)=b_{\alpha}$ and $s_{\alpha}(-L)=0$ asymptotically for sufficiently large $L$. This in turn implies that $\tau_{\alpha}(s_{1}(\pm L),s_{2}(\pm L))\simeq 0$. Thus, if the domain of slip is modeled as $[-L,L]$ along the $x$ axis, for some sufficiently large $L$, then the Fourier transform (12) of $\tau_{\alpha}(s_{1}(x),s_{2}(x))$ can be approximated using a Fourier series of the form

[TABLE]

where

[TABLE]

and the Fourier coefficients $d_{\alpha,n}$ in (15) are related to the coefficients $d_{\alpha}(k)$ in (12) as

[TABLE]

In practice, the summation in (15) is carried over a finite range $n\in[-N_{\tau},N_{\tau}]$. To obtain analytical expressions for $\hat{d}_{\alpha}(k)$ (and thereby the coefficients $\{d_{\alpha,n}\}$), the slip distribution is further approximated as a piecewise linear function over a suitable discretization of the domain $[-L,L]$: for $j=0,1,\ldots,(N_{s}-1)$,

[TABLE]

Here $(x_{0},\ldots,x_{N_{s}})$ is a partition of $[-L,L]$, with $x_{0}=-L$ and $x_{N_{s}}=L$, $h_{j}=x_{j+1}-x_{j}$, and $s_{\alpha,i}=s_{\alpha}(x_{i})$. Using the piecewise linear approximation (18) for the slip distribution, the Fourier expansion (9) of the stacking fault energy, and the definition of the misfit stress (4), the coefficients $\hat{d}_{\alpha}(k)$ of the Fourier transform of $\tau_{\alpha}(s_{1}(x),s_{2}(x))$ in (12) can be computed analytically as follows:

[TABLE]

With these approximations in place, the generalized Peierls-Nabarro equation (8) can be written as

[TABLE]

It is emphasized that the foregoing equation is nonlinear on account of the dependence of the coefficients $d_{\alpha,n}$ on the piecewise linear slip distribution $\{s_{\alpha,i}\}$. Using (41), the right hand side of (20) can be expressed using Hilbert transforms, thereby resulting in the following equation:

[TABLE]

where sgn is the sign function. Using (39) and taking the inverse Hilbert transform of (21) yields the following equation:

[TABLE]

where $a_{\alpha,1}$ is a constant that is fixed by the boundary conditions. Since (22) is obtained by discretizing the slip distribution using the piecewise linear approximation (18), it is henceforth referred to as the discrete inverse Peierls-Nabarro model. A simple algorithm to solve the discrete inverse PN model is presented next.

2.3 Iterative scheme for the inverse Peierls-Nabarro equation

The discrete inverse Peierls-Nabarro equation (22) is a nonlinear and non-local equation. One of the advantages of this equation over the original PN equation (3) is that integrals with the singular ‘$1/r$’ kernel have been eliminated. A further simplification can be achieved by solving the discrete inverse Peierls-Nabarro equation (22) in an iterative fashion using a relaxed Picard iteration scheme: if $s_{\alpha}^{(k)}$ denotes the piecewise linear slip distribution at the $k^{\text{th}}$ iteration, then a predictor $\tilde{s}_{\alpha}^{(k+1)}$ of the solution of (22) at iteration $k+1$ can be computed by solving

[TABLE]

Equation (23) can be integrated exactly to yield

[TABLE]

where $a_{\alpha,2}$ is an arbitrary constant. To fix the constants $a_{\alpha,1}$ and $a_{\alpha,2}$ in the predictor (24), the following boundary conditions are used

[TABLE]

Here $S_{\alpha}$ is the total slip across the slip plane along the $x_{\alpha}$ direction. Using these boundary conditions, the constants $a_{\alpha,1}$ and $a_{\alpha,2}$ in the predictor (24) are computed as

[TABLE]

The use of these expressions for the constants in (24) yields the following analytical expression for the the predictor $\tilde{s}^{(k+1)}_{\alpha}$:

[TABLE]

For well understood reasons related to numerical stability, the actual update that yields the slip distribution $s_{\alpha}^{(k+1)}$ at iteration $k+1$ from the slip distribution $s_{\alpha}^{(k)}$ at iteration $k$ is formulated as

[TABLE]

where $\lambda\in(0,1]$ is a (small) parameter (see Appendix B). The iterations are carried out until $\sum_{\alpha}\left(\sum_{j}(s^{k+1}_{\alpha,j}-s^{k}_{\alpha,j})^{2}\right)^{\frac{1}{2}}<\text{TOL}$, where TOL is a specified tolerance. It is emphasized that the foregoing iterative strategy is local at every iteration, and is thus very efficient. Notably, this iterative scheme is matrix-free and does not require the formulation of a linear algebraic system of equations, unlike existing techniques to solve the generalized PN model.

In practice, the fixed point iteration scheme developed above converges slowly. A variety of techniques have been developed to accelerate fixed point iterations; these are briefly reviewed in Appendix B. The Anderson acceleration scheme with $n_{AA}$ stages is used in this work to accelerate the convergence of the proposed iterative scheme to solve the generalized PN model (3).

2.4 Effect of applied stresses

The effect of an external stress on the slip distribution, while straightforward in the generalized PN model, is not trivial within the framework of the inverse PN model proposed here. This is due to the fact that the Hilbert transform of a constant is zero (see (39)). A numerical strategy that circumvents this problem is proposed herein.

Taking advantage of the fact that the iterative scheme proposed in this work is defined over the domain $[-L,L]$, the effect of an externally applied stress $F_{\alpha}$ is modeled as

[TABLE]

where $\chi_{[a,b]}(x)$ is the characteristic function of the interval $[a,b]$ that is equal to $1$ when $x\in[a,b]$ and [math] otherwise. Using (43), the external stress can be modeled as follows:

[TABLE]

The generalized PN model (3) in the presence of an applied stress $F_{\alpha}$ takes the form

[TABLE]

Following the same procedure outlined earlier in obtaining (21) from (3), (31) can be written as

[TABLE]

Taking the inverse transform, the resulting inverse PN model can be solved analytically as before at each iteration of the iterative scheme.

While this strategy provides a way to handle external stresses within the current framework, it has its limitations owing to the fact that the domain size $L$ explicitly figures in (32). The specific choice (29) used to model external stresses requires careful handling, especially close to the boundaries. This and related details are discussed in a later section.

3 Numerical examples

A few numerical examples are now presented to illustrate the utility of the proposed iterative scheme for the inverse Peierls-Nabarro model. This includes simple test cases using the original PN model (6), dissociated dislocations in one dimension, the equilibration of a dipole to illustrate the handling of external stresses, and a study of the core structure of edge and screw dislocations on the close-packed $\{111\}$ planes in Aluminium.

In all the simulations presented below, $h$ refers to the mesh discretization size, and is defined as $h=\text{max}(h_{j})$. When reporting the results, two choices of discretization, namely a fine discretization $h\to 0$ and a coarse discretization with the atomic scale periodicity $h_{j}=b$ are chosen. The dislocation density in each case is computed numerically from the slip distribution data by employing a central difference approximation.

3.1 Original PN model

The availability of an exact analytical solution for the simple case of a sinusoidal misfit stress, as in (7), provides a simple benchmark to test the proposed numerical strategy to solve the PN equation. A comparison of the slip distribution and the dislocation density obtained by the iterative solution of the inverse PN model and the exact analytical solution is shown in Figure 2.

The numerical solution was obtained with the following choice of parameters: $b=1.0$, $\xi=1.0$, $L=50.5$, $N_{\tau}=202$ and $\lambda=0.01$. Two approximate solutions, one corresponding to a fine discretization with $N_{s}=801$, and a coarse discretization with $N_{s}=101$, were considered. $N_{s}=801$ corresponds to a mesh size of approximately $b/8$, while $N_{s}$ corresponds to the physically relevant mesh size of $b$. The initial guess for the slip distribution was chosen as a step function with slip equal to [math] in $[-L,0)$, and equal to $b$ in $[0,L]$. The iterations were carried out until the error in the slip update is less than a specified tolerance level of $\text{TOL}=10^{-6}$. It is seen that both the fine and coarse meshes produce results in good agreement with the exact solution, with the finer mesh being more accurate than the coarser mesh - this will be discussed in more detail shortly.

The naive fixed point iteration used in obtaining Figure 2 requied $1002$ iterations when $h=b$, $1154$ iterations when $h\simeq b/8$. To accelerate the convergence, the Anderson acceleration technique for fixed point iteration (see Appendix B) with $n_{AA}$ stages was used. The number of iterations required to reach a tolerance level $\text{TOL}<10^{-6}$ is shown in Table 3.1 and Figure 3. The other parameters were kept fixed at the same values as before: $b=1.0$, $\xi=1.0$, $L=50.5$, $N_{\tau}=202$, and $N_{s}=101,801$. The starting guess for the slip distribution was set as the same step distribution as before. It is seen that there is an optimal value of $n_{AA}$, in this case close to $10$, for which the number of iterations for convergence is minimum. This value, however, is dependent on the problem at hand. The key observation is that the number of iterations with Anderson acceleration is $O(10)$ - a significant improvement compared to the $O(10^{3})$ iterations required for convergence in the naive iteration scheme. For this reason, all numerical results presented henceforth are those obtained using the Anderson accelerated fixed point iterations.

To monitor the effect of the domain size $L$ on the accuracy of the iterative solution, the error $\epsilon$ between the iterative and exact solutions, defined as

[TABLE]

is plotted as a function of $L$ in Figure 4, and listed in Table 3.1, for both the fine ($h\simeq b/8$) and coarse ($h=b$) mesh discretizations. Here $I_{\tilde{L}}$ is the set of indices with the property that $x_{i}\in[-\tilde{L},\tilde{L}]$ whenever $i\in I$, and $N_{I_{\tilde{L}}}$ is the number of such indices. For concreteness, $\tilde{L}=10b$ was chosen in computing the errors shown in Figure 4. The various constants are chosen as follows: $b=1.0$, $\xi=1.0$, $N_{\tau}=2(2L/b)$, and $n_{AA}=10$. $N_{s}=2L/b$ for the coarse mesh, and $N_{s}=(16L/b)-7$ for the fine mesh. The number of iterations required for convergence is also listed in Table 3.1. As can be seen from Figure 4 and Table 3.1, the error associated with the fine discretization decreases rapidly with the domain size, while the error in the coarse mesh discretization initially decreases rapidly, but eventually decreases slowly with increase in domain size.

To understand the origin of the slower decrease in the error associated with the coarse discretization $h=b$, as seen in Figure 4, a plot of the error $\epsilon_{10b}$ is plotted as a function of the domain size $L$ for various values of the dislocation width $\xi$ in Figure 5. The various parameters used in obtaining Figure 5 are the same as those used for the simulations in Figure 4. It can be seen that the error corresponding to wider dislocation profiles is lower than that of narrower dislocation profiles when the domain size is sufficiently large. Stated differently, the error associated with the coarse mesh discretization is lower when $b/\xi$ is small, or, equivalently, when the slip distribution varies slowly at the scale of the atomic lattice.

3.2 Dissociated dislocations in 1d

The slip distribution corresponding to a 1d PN model with a more complex slip distribution is considered in this section. Specifically, the misfit stress is chosen to be of the form

[TABLE]

In this case, there is no exact analytical solution available. Further, the presence of a local minimum in the misfit stress distribution is expected to cause a dissociation of the dislocation into partial dislocations separated by a finite width stacking fault. To validate the solution obtained by the proposed method, the 1D PN model was also solved using the semi-discrete variational Peierls-Nabarro (SVPN) model [10]. A comparison of the solution obtained by these two methods is shown in Fig. 6.

The following parameters were used for the reference SVPN solution: $b=1.0$, $\xi=2.5$, $L=50.0$, $N_{s}=2000$. For the inverse PN model, the following parameters were used: $L=50.5$, $N_{s}=201,1601$, $N_{\tau}=402$, $\lambda=0.01$ and $n_{AA}=20$. It can be seen that the proposed method agrees well with the dissociated core structure that is computed using the SVPN method. This example is included to highlight the fact that the proposed iterative scheme for the inverse PN model is competitive with existing techniques like the SVPN method for non-trivial misfit stresses.

3.3 Effect of external stress

As mentioned earlier, the fact that the Hilbert transform removes the DC component of the function it acts on requires that external stresses are handled in a non-standard manner. Specifically, the effect of an external stress is represented by a term of the form $F\chi_{[-L,L]}$. However, it is to be noted that since the current formulation is based on the original PN model, it inherits one its fundamental limitations: the energy of an infinite linear elastic medium with a single dislocation is independent of the position of the dislocation. The fact that the present formulation enforces the slip to have definite values at the boundaries of a finite domain $[-L,L]$ introduces, however, an artificial dependence of the energy on the slip distribution over the finite domain $[-L,L]$ via a confining potential, whose origin lies in the choice of the approximation (29) for external stresses within the current framework.

This is illustrated in Figure 7, where the center of the dislocation, defined as

[TABLE]

is plotted as a function of the external force $F$. The graph shown in Figure 7 was obtained by solving (32) with sinusoidal misfit stress, and with the following set of parameters: $b=1$, $\xi=1.0$, $L=60.0$, $N_{s}=1201$, $N_{\tau}=240$, and $n_{AA}=20$. The slope of the initial part of the force-displacement curve in Figure 7 is $1.0$ (within numerical errors). Thus, the dislocation experiences no lattice resistance, as expected in the original PN model. The slope decreases towards zero as the dislocation approaches the boundary, since the boundary conditions require the slip to be zero on the boundaries, and hence prevent further dislocation motion. The non-linearity in the force-displacement curve in Figure 7 is thus a direct consequence of the finiteness of the domain considered. The fact that the dislocation has a definite position for a given valued of the external stress illustrates the artificial confining potential introduced by modeling the external stress as in (29). The foregoing result also shows that the Peierls stress cannot be calculated within the present approach. The Peierls stress may be calculated as in the conventional procedure by using the discrete solution obtained by the present formulation in a discrete approximation of the energy which recovers the effect of lattice periodicity [28]. The computation of the Peierls stress along these lines is not carried out in this work.

To illustrate the foregoing comments on handling external stresses within the current formulation, as in (32), the equilibration of a dislocation dipole is considered. Consider two edge dislocations of equal and opposite Burgers vector $b$ separated by a distance $2w$ on a glide plane. In the absence of an external stress, the linear elastic theory predicts that this dipole will annihilate itself. The stress $\sigma_{e}$ required to stabilize this dipole can be easily computed using the linear elastic theory as

[TABLE]

The equilibration of the dipole is now reconsidered using the inverse PN model (32) using the sinusoidal stacking fault energy (5). Fixing the width of the dipole to $w=40.0$, the dependence of the equilibration stress as a function of the domain size is shown in Figure 8 and Table 3.3. The following parameters were chosen: $N_{s}=10(2L/b)+1$, $N_{\tau}=4L/b$, $n_{AA}=10$ and $\lambda=0.01$. It is seen that the stress required for equilibration converges if the domain is large enough.

A comparison of the (scaled) external stress required for equilibration of the dipole as a function of its width is shown in Figure 9 and Table 3.3. The following parameters were used: $b=1$, $\xi=1$, $L=200$, $N_{s}=4001$, $N_{t}=800$, $\lambda=0.01$ and $n_{AA}=20$. It can be seen that the predictions of the sinusoidal PN model, computed using the iterative scheme proposed in this work, are in qualitative agreement with the predictions of linear elasticity. The expected deviations from the linear elastic theory that is seen in Figure 9 arise due to two distinct sources. First, there is a modeling error due to the fact that the PN solution is computed assuming that the stacking fault energy is sinusoidal in nature, as in (5). Second, referring to Figure 8, the error associated with smaller dipole widths are smaller since $L/w$ is larger than that for dipoles with larger width. As a final comment, it is important to place the dipole at the center of the domain since boundary effects that arise due to the specific manner in which external stresses are modeled here become important as the dipole approaches the boundaries.

3.4 Core structure of edge and screw dislocations in Aluminium

While the foregoing examples involved toy models with one-component slip distributions, a more realistic problem concerning the core structure of both edge and screw dislocations in Aluminium is now analyzed using the inverse Peierls-Nabarro model. The core structure of dislocations in Al is well studied (see, for instance, [24], [29], [23], [6], [11]), and hence serves as a good benchmark to validate the proposed iterative method to solve the generalized PN model. As is well known, perfect dislocations on the $\{111\}$ planes in FCC materials split into Shockley partials for energetic reasons [4]. One of the goals of this study is to reproduce the core structure of these dissociated dislocations and compare them with results obtained earlier.

For concreteness, dislocations on the close-packed $\{111\}$ planes in Aluminium were considered. Specifically, edge dislocations, with line direction $t=[\bar{1}12]$, and screw dislocations, with line direction $t=[110]$, were studied. The symmetric anisotropic Stroh tensor [27] is used to model the elastic interactions. If $n$ denotes the normal to the $\{111\}$ plane on which the dislocation resides, the direction $m$ perpendicular to the dislocation line is given by the vector cross product of $n$ and $t$: $m=n\times t$. The components of the symmetric anisotropic Stroh tensor on the $\{111\}$ plane for Aluminium are tabulated in Table 3.4; these values are based on those reported in [30].

The generalized stacking fault energy surface is computed using the fitting procedure outlined in [24]. The generalized stacking fault energy on $\{111\}$ planes is taken to be of the form

[TABLE]

where

[TABLE]

$b=a/\sqrt{2}$, and $a=4.05$Å is the lattice constant of Al. The various constants in (35) are listed333Some of the constants are listed incorrectly in [24]. The corrected version of these constants are reported here. in Table 3.4.

These constants are obtained by requiring that the expression (35) for the stacking fault energy agree with known values of the stacking fault energy (computed, for instance, from lower scale models) for specific stacking configurations. These configurations, referred to as $A,T,T_{1},G,G_{1},G_{2},G_{3}$ following [24], are shown in the contour plot of the stacking fault energy on $\{111\}$ planes in Al is shown in Figure 10. The stacking configurations used to fit the constants in (35) are listed in Table 3.4 - these values are taken from [24].

The slip and dislocation density distributions corresponding to an edge dislocation on the $(1\bar{1}1)$ plane oriented along the $[\bar{1}12]$ direction are shown in Figure 11. The corresponding distributions for a screw dislocation on the $(1\bar{1}1)$ plane with line direction $[110]$ are shown in Figure 12. The distributions in Figure 11 and Figure 12 are shown for different domain sizes and different mesh sizes to highlight the convergence of the numerical procedure for various discretization choices. The parameters used to obtain these figures are as follows: $L=50.5b^{\prime},75.5b^{\prime}$, $N_{s}=101,801,1201$, $N_{\tau}=202,302$, $n_{AA}=25$. Here, $b^{\prime}=\sqrt{3}b$ for screw dislocations and $b^{\prime}=b$ for edge dislocations on $\{111\}$ planes. The simulations for all these cases converge within a few hundred iterations, and hence are competitive in terms of computational efficiency with matrix based techniques like the SVPN method.

The dislocation density profile is shown in Figure 11 and Figure 12 for an edge and screw dislocation, respectively. It is seen that the edge dislocations splits while the screw dislocation has a wide core that may be thought of as two very closely spaced partials. These results agree with those present in [24], [29], [23]. The separation of the partials in the case of edge dislocations is approximately 8.0Å, which is in good agreement with experimental results as reported in [29]. It is pointed out that unlike the SVPN method which predicts the splitting of the core for screw dislocations [31], the current method does not; this is in line with the predictions made in [23].

The slip distribution of both the edge and screw dislocations on the $(1\bar{1}1)$ plane is more conveniently visualized by plotting the components of slip along the $[110]$ and $[\bar{1}12]$ directions over the stacking fault energy surface, as shown in Figure 13. This clearly shows that the core of the dislocation is characterized by the splitting of the dislocation into Shockley partials, as expected. The iterative scheme proposed in this work for solving the generalized PN model is thus useful to study the core structure of dislocations in real materials, given appropriate generalized stacking fault energy data.

4 Discussion and Summary

The results presented so far indicate that the iterative scheme proposed in this work is a reliable and computationally efficient alternative to solve the generalized Peierls-Nabarro model. It is emphasized that the matrix-free nature of the present approach allows for efficient parallelization of the iterative solution scheme. A number of issues pertinent to the iterative scheme proposed in this work to solve the generalized PN model and topics for future investigations are now discussed, followed by a summary of the key ideas presented in this work.

4.1 Handling the constant term in the misfit stress distribution

It is to be noted that the constant term $d_{0}$ in the Fourier representation (12) of the misfit stress distribution $\tau_{\alpha}(s_{1}(x),s_{2}(x))$ on the slip plane is handled in the same manner as an external stress. Specifically, the misfit stress distribution is modeled as

[TABLE]

In the 1D case where the stacking fault energy is of the form (5), it is expected that $d_{0}=0$. This has been verified numerically.

4.2 Choice of initial guess for slip distribution

The convergence of the iterations depends on a good initial guess for the slip distribution. In particular, an initial guesses based on the approach developed in [32] can be used, as is now illustrated in a 1D setting. Assuming that the stacking fault energy is piecewise quadratic of the form

[TABLE]

where $G$ is a constant, $m$ is related to the amount of plastic slip, the inverse PN model reduces to the following equation

[TABLE]

Assuming further that $s$ and $m$ admit a Fourier expansion of the form

[TABLE]

where $k_{0}=\pi/L$, it is straightforward to verify that

[TABLE]

This approach has been used successfully to develop a phase field model for dislocation dynamics in [21]. Within the context of the iterative scheme proposed in this work, the foregoing solution can be used as an initial guess with $G$ computed from the curvature of the stacking fault energy at zero slip. While this provides a good initial guess for the iterative scheme proposed here, a simple piecewise constant slip distribution was found to be sufficient in practice.

4.3 Computing the Peierls stress

As remarked earlier, the current framework doesn’t provide a direct means to compute the Peierls stress. A lower bound for the Peierls stress may, however, be obtained by estimated by considering a coarser mesh discretization. In the (unphysical) limit when the mesh size $h=\text{max}(h_{j})\to 0$, the discretization gets closer to the continuum PN model, as illustrated by the numerical examples presented earlier. However, as is the case with the original PN model, information about the discreteness of the atomic lattice is lost in the limit $h\to 0$. A possible solution to introduce lattice discreteness directly into the current model is to make the choice $h_{j}=b$ (for every $j$) and consider the cases when $L/b$ is large. This is a physically meaningful choice since lattice misfits due to the presence of a dislocation are only truly defined at discrete locations corresponding to the perfect lattice points on the slip plane. This is similar to the procedure followed in the conventional PN model to compute the Peierls stress, where the energy functional is evaluated as a discrete sum of its values at lattice points, thereby introducing the necessary periodicity in the energy landscape. It is to be noted that the PN model is no longer valid in such a setting, and alternative approaches like the discrete dislocation equation [22] that generalize the PN model to a corresponding discrete model may need to be used. Approximating the PN model directly with the coarse mesh whose size is equal to the lattice spacing thus introduces errors, as can be seen from Figure 2, Figure 4, Figure 5 and Figure 6. An estimate of the Peierls stress can be obtained by following a procedure similar to that adopted in the SVPN method [10] as the external stress at which the generalized PN model fails to converge, or using the procedure describe in [33]. The solution of the inverse Peierls-Nabarro model in the presence of an external stress requires additional care since the domain size $L$ directly figures in the approximate model of external stress (29), and consequently in (32). Preliminary simulations showed domain size effects in the estimates of the Peierls stress. Alternatively, the approach adopted in [33] to compute the Peierls stress could be employed here. The continuum energy of a crystal with a dislocation is modified in [33] to include terms that reflect the periodicity of the lattice on the slip plane. This is accomplished by taking advantage of the specific ansatz introduced in [8] and using it in conjunction with the approach developed by Nabarro [2] to estimate the Peierls stress. The use of the discrete Hilbert transform right from the start is another modeling choice that could prove to be helpful. These and related issues pertaining to the computation of the Peierls stress within the framework of the inverse Peierls-Nabarro model will be investigated in a future work.

4.4 Nonlocal misfit energy

To overcome the limitations of the use of the generalized stacking fault energy to represent the misfit energy due to a specified slip distribution on the slip plane, nonlocal terms in the misfit energy have been introduced in [20]. The present formulation motivates in a natural manner how such non-local terms may be introduced. To see this, the property of the inverse of the Hilbert transform (38) and the property (39) are used to transform (8) to get

[TABLE]

upto an additive constant that is determined by boundary conditions as before. The form of this equation suggests a natural non-local extension of the inverse PN model as follows:

[TABLE]

where $A_{\alpha}(s_{1}(x),s_{2}(x);s_{1}(x^{\prime}),s_{2}(x^{\prime}))$ is a non-local term that can be modeled using techniques like those developed in [20]. An iterative strategy for such non-local models can be developed in a manner analogous to the one presented in this work. These non-local extensions that extend the generalized PN model using data from lower scale models will be investigated in a future work.

4.5 Extension to two dimensions

The special properties of the Hilbert transform for functions of one variable have been exploited in the current work to develop the iterative scheme. To extend the current approach to two dimensions, as will be required for studying slip distributions on the slip plane, an extension of the Hilbert transform to handle the ‘$1/r$’ kernel in (3) in dimensions larger than one that simultaneously preserve all the nice properties of the 1D Hilbert transform is desired. A variety of extensions of the Hilbert transform to higher dimensional equivalents have been proposed in the past: [34], [35]. The extension of the iterative scheme developed in this work to the two-dimensional setting is also planned for a future work. It is worth reiterating, however, that the generalized PN model (3) has been widely used to study the dislocation core structure in realistic materials despite the fact that it requires the slip distribution to be a function of one variable only. The iterative scheme developed in this work can thus be useful to study dislocation core structures in materials of practical interest.

4.6 Conclusion

To summarize, a novel numerical method that exploits certain properties of the Hilbert transform to solve the generalized Peierls-Nabarro model by reducing it to a fixed point iteration scheme has been developed in this work. The key advantage of the new method in comparison with standard approaches to solve the PN model is that it is matrix-free since it is of the form of a local fixed point iteration which is computationally efficient and is therefore amenable to parallelization schemes. The form of the inverse PN model further suggests in a natural manner extensions to non-local models that is compatible with the iterative structure developed here. The utility of the proposed method has been illustrated with a variety of examples. In particular, it is demonstrated that the proposed numerical scheme can handle realistic dislocation cores by illustrating its application to study the core structure of edge and screw dislocations in Aluminium on the close-packed $\{111\}$ planes. The dissociation of edge dislocations, the width of the stacking fault, and the diffuse core structure of screw dislocations that have been reported in previous works are captured satisfactorily by the present model. A possible strategy to include external stresses within the current framework and its limitations are discussed. The incorporation of external stresses in the inverse PN model so as to predict the Peierls stress and the change in the core structure of dislocations as it moves across the lattice will be carried out in a future work.

Appendix A Hilbert transforms

The Hilbert transform [36] of a function $f:\mathbb{R}\to\mathbb{R}$ is defined as the function $\mathcal{H}(f):\mathbb{R}\to\mathbb{R}$ given by

[TABLE]

where

[TABLE]

is the Cauchy principal value of the improper integral in the Hilbert transform (37). The Hilbert transform is widely used for applications in signal processing. Some useful properties of the Hilbert transform that are relevant to this work are summarized here.

The inverse Hilbert transform of a given function is equal to the negative of the Hilbert transform of the function: if $\mathcal{H}(f)(x)=g(x)$, then $\mathcal{H}(g)(x)=-f(x)$. Thus

[TABLE] 2. 2.

The Hilbert transform of a constant $k$ is zero:

[TABLE] 3. 3.

The Hilbert transform of complex exponentials have an especially simple form:

[TABLE]

where $\text{sgn}(\omega)$ is the sign of $\omega$. Thus, the Hilbert transform rotates a sinusoidal signal by $\pi/2$. As a particular case of the foregoing result, note that for any $\omega>0$,

[TABLE] 4. 4.

If $\chi_{[a,b]}$ denotes the characteristic function for the set $[a,b]$, then

[TABLE]

Using (42) and (38), it is easy to see that

[TABLE]

In the main text, $\mathcal{H}(f)(x)$ is often written as $\mathcal{H}(f(x))$.

Appendix B Acceleration schemes for fixed point iteration

The convergence of the iterative scheme (27) presented earlier to solve the Peierls-Nabarro model could be slow, as is typically the case with fixed point iteration strategies. A variety of acceleration schemes have been developed in the past to speed-up fixed point iterations. These are briefly outlined here in the context of solving the nonlinear fixed point equation

[TABLE]

where $x\in\mathbb{R}^{n}$ and $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$. The naive fixed point iteration, also known as Picard iteration, to find solutions of (44) takes the form

[TABLE]

A relaxed form the iterative scheme (45) is obtained, for a given choice of $\lambda\in(0,1]$, as

[TABLE]

A useful means to rewrite (46) is as follows:

[TABLE]

which can be thought of as the forward Euler discretization of the differential equation

[TABLE]

Restrictions on the size of $\lambda$ thus correspond to well known restrictions on the choice of time step for the numerical solution of ordinary differential equations. Indeed, methods that accelerate the solution of the nonlinear equation (48) using specially designed Runge-Kutta methods have been proposed (see [37]). The iterative scheme (28) for the inverse Peierls-Nabarro equation (27) proposed in this work is accelerated using an alternative strategy known as Anderson acceleration [38]. The $n_{AA}$-stage Anderson acceleration for the nonlinear fixed point equation (44) is computed using $n_{AA}$ successive iterates of a fixed point iterative update, like (45). The basic idea behind the $n_{AA}$-stage Anderson acceleration scheme is to choose constants $\{\alpha_{1},\ldots,\alpha_{n_{AA}}\}$ such that they minimize the residual

[TABLE]

subject to the constraint

[TABLE]

Thus, the current iterate is chosen as the particular linear combination of the $n_{AA}$ preceding iterates that minimizes the joint residual. Specific details regarding the numerical implementation of the foregoing minimization problem can be found in [38]. For the present work, the first $n_{AA}$ iterations are carried out using the relaxed Picard iteration scheme (46) with a small valued of $\lambda$ to bootstrap the process. Subsequent iterations are computed using the $n_{AA}$-stage acceleration scheme as outlined above.