Multiscale modeling of dislocations: Combining peridynamics with gradient elasticity

Modeling dislocations is an inherently multiscale problem as one needs to simultaneously describe the high stress fields near the dislocation cores, which depend on atomistic length scales, and a surface boundary value problem which depends on boundary conditions on the sample scale. We present a novel approach which is based on a peridynamic dislocation model to deal with the surface boundary value problem. In this model, the singularity of the stress field at the dislocation core is regularized owing to the non-local nature of peridynamics. The effective core radius is defined by the peridynamic horizon which, for reasons of computational cost, must be chosen much larger than the lattice constant. This implies that dislocation stresses in the near-core region are seriously underestimated. By exploiting relationships between peridynamics and Mindlin-type gradient elasticity, we then show that gradient elasticity can be used to construct short-range corrections to the peridynamic stress field that yield a correct description of dislocation stresses from the atomic to the sample scale.


Introduction
Peridynamics as introduced by Silling (2000) is a nonlocal continuum theory based on the formulation of integrodifferential equations for the displacement field, without directly involving classical concepts such as stress and strain.Like other generalized continuum approaches such as nonlocal elasticity (Eringen 1983) or gradient elasticity (Mindlin 1965;Mindlin and Eshel 1968;Kröner 1967), it can serve to regularize problems that, in classical elasticity theory, are associated with singular or discontinuous solutions.Peridynamics shares this regularization property with higher-order continuum theories such as gradient elasticity (see e.g.Silling and Lehoucq (2008)), and several authors have proposed matching schemes between peridynamics and higher-order continua (Dayal 2017;Chen and Chan 2020).
Here we use bond-based peridynamics to explore how peridynamics can be used to model dislocations -a paradigmatic multiscale problem as, in classical elasticity, the elastic fields of dislocations exhibit singularities both at infinity (where they are regularized by surface boundary conditions) and at the dislocation line (where they are regularized by atomic-scale non-local interactions).One thus needs to establish a consistent description of dislocation stress fields from the atomic to the sample scale.In the section "Structure of bond-based peridynamics and relations with gradient elasticity", we give a brief introduction into bond-based peridynamics and derive its relations with bulk gradient elasticity.In particular, we show how to relate the micromodulus and length scale parameter ('horizon') of peridynamics to the elastic constants and higher-order elastic constants of Mindlin's first strain gradient elasticity.In section "Constructing dislocations in peridynamics", the introduction of dislocations into a peridynamic framework is discussed, using concepts of continuum dislocation theory which relate dislocations to spatial derivatives of the plastic distortion.We first consider singular dislocations which, in generalization of Volterra's construction, can be considered in terms of displacement discontinuities across slip planes; in this formulation, the plastic distortion localized on a slip plane leads to force-free changes in bond vectors that cross this plane.Similar ideas have been formulated previously by other authors, see in particular Zhao et al. (2021); the formulation here differs from theirs in detail as the horizon of peridynamic interactions is defined in an intermediate (plastically deformed) configuration rather than the material reference configuration.This allows to correctly account for the physical implications of changes in atomic neighborhood relations.The formulation is then generalized to dislocations of distributed Burgers vector, where accordingly the plastic distortion is spatially distributed around the dislocation slip plane and the force-free change of the bond vector must be evaluated as a line integral over the plastic distortion along the bond.For both singular and distributed dislocations, it is shown that the dislocation related fields exhibit a peculiar scaling form which is common to the classical, singular solutions obtained from linear elasticity, to solutions computed in the framework of gradient elasticity, and to the results obtained numerically from peridynamics.Section "Simulations and comparison with gradient elasticity" provides a brief discussion of numerical methodology and then presents numerical results obtained for singular as well as distributed edge dislocations in the bond-based peridynamic framework.For singular dislocations where the plastic distortion is localized on the slip plane, it is found that the peridynamic solution removes the stress and strain singularities at the dislocation line, but still exhibits non analytic behavior as the asymptotic solution, that is approached when the number of collocation points in the horizon is increased, exhibits a discontinuous jump across the dislocation line.It is then shown that this jump is removed when dislocations with distributed Burgers vector are considered, and numerical results obtained with appropriate Burgers vector distribution are shown to be in excellent agreement with theoretical results for dislocation stress fields in gradient elasticity.Section "Construction of a hybrid multiscale dislocation model" shows how the results can be used to construct a multiscale description of dislocations.Several such models have been proposed in the literature (Huang and Li 2015;Jamond et al. 2016;Lu et al. 2019Lu et al. , 2022;;Bertin and Capolungo 2018;Bertin 2019).These models use numerical continuum approaches, in the form of FEM or spectral methods, to evaluate stresses associated with the dislocation network (or the corresponding eigenstrain) on the sample scale.This is combined with analytical results that describe the stress fields close to the individual dislocation lines on the level of dislocation segments.In a similar spirit, we propose in the present study to use peridynamics for describing the elastic field far from the dislocation and the associated boundary value problems, while gradient elasticity results are used to evaluate corrections that provide an accurate description of the stresses close to the dislocation line.Finally, in "Discussion and conclusions" we compare our approach with other multiscale dislocation models and outline potential applications.

Bond-based peridynamics
In bond-based peridynamics model, as originally defined by Silling (2000), the deformation of a D dimensional continuous body B of density ρ(x) is described by the displacement field u(x) where x are material coordinates.The force balance equation for the material point x is given by where f x, x ′ is a pair force density between x ′ and x , b is a body force field, and inter- actions are restricted to an interaction domain H x which we take to be a D dimensional sphere of radius δ p around x , the so-called horizon: |x − x * | ≤ δ p ∀ x * ∈ H x .The pair force density is specified constitutively.To this end, the bond stretch is defined as where ξ = x ′ − x is the bond vector and η = u ′ − u is the relative bond displacement.The pair force density is then taken to be linearly proportional to the bond stretch s, and pointing in the direction of the vector e u (ξ , η) connecting both points in the current configuration: Here c(ξ ) is the so-called bond micro-modulus which for an isotropic bulk material depends on the bond length ξ = |ξ | only.The bond energy density can then be written in terms of the bond length ξ and bond stretch s as The total elastic energy of the body B is obtained by integrating over all pairs of interact- ing material points: Note that the factor 1 2 is needed in order not to count bonds twice in the double integration.In Eq. ( 5), we have for comparison with classical and gradient elasticity theories introduced the strain energy density W p (x) associated with the point x , defined as ( For the purpose of comparing peridynamics with classical or gradient elasticity theories, it is convenient to express the energy density in the limit of small deformations, |η| ≪ |ξ | , to obtain a force balance equation that is linear in the displacement field u , and an energy functional that is a quadratic form of u (Silling and Lehoucq 2010).This gives where and e ξ = ξ /ξ .The strain energy density W p (x) associated with the point x is then written as This expression serves as our starting point for relating peridynamics to gradient elasticity.

Relation with gradient elasticity
To establish the relationship between bond-based peridynamics and higher-order strain gradient elasticity, we expand the displacement u into a Taylor series around the point x .We define the first, second and third-order displacement gradients via where we use the notation u j,i := ∂u j /∂x i and similarly for higher-order derivatives.The expansion of the displacement then reads where we use the standard summation convention.Upon insertion in Eq. ( 9), we obtain an expansion of the energy density in terms of higher-order displacement gradients: We separate radial and angular parts of the integrals by writing ξ i = ξ e i where the e i are components of a unit vector in ξ direction and depend only on angular coordinates: where D is the solid angle in D dimensions and d the corresponding angle element.Now it is easy to see that the angular integrations give zero unless the indices of the e i factors are pair-wise equal.At lowest order we find: where For comparison with strain-based elasticity theory formulations it is convenient to split the deformation gradient into symmetric and antisymmetric parts, The antisymmetric part of the first-order deformation gradient does not contribute to the energy, and we write the lowest-order elastic energy contribution as which is the elastic energy of a classical linear-elastic material with Lamé constants and µ = , hence Poisson number ν = 1/4.
We turn to the first-order contribution, following a similar line of argument.We write and introduce the characteristic length l via The first-order coupling constants can then be written as ( 15) e i e j e k e l d�, (16) 1 16 e i e j e k e l e m e n d�, (17) where we have grouped terms that, after summation over equal indices, yield equal energy contributions.Using the symmetry f ′ ijk = f ′ jik and re-labeling we can write the gradient contribution to the elastic energy as Again we may perform a transition to strain gradients.Defining the strain gradient via η ijk := ǫ jk,i , we obtain Either way, we find that up to first-order strain gradients the energy functional corresponds to a special case of Mindlin's energy functional for isotropic first gradient elasticity (Mindlin and Eshel 1968), where the five strain gradient coefficients a 1 . . .a 5 are given by For standard parameterizations of the micro-modulus functions, the length l is proportional to the horizon radius δ p , e.g. in 2D, for a constant micro-modulus c(ξ ) = c 0 , we obtain l 2 = 3 5 δ 2 p , and for a 'conical' micro-modulus function c(ξ

General formalism
The approach to constructing dislocations in peridynamics builds upon concepts from the classical continuum theory of dislocations.This theory relates dislocations to gradients of the plastic distortion, which in turn is envisaged in terms of slip on crystallographic slip planes.In the following a single dislocation line is considered as schematically illustrated in Fig. 1 for an edge dislocation.We take the slip plane without loss of generality to be the plane y = 0 with normal vector n = e y , and the Burgers vector as b = be x .The dislocation line is parameterized as r d (s) , where r d • n = 0 ∀ s and s is the line length.The unit tangent vector to the dislocation line is denoted by t(s) = ∂ s r d and g = t × n is the dislocation glide vector. (23) In the spirit of a continuum theory, we allow for the dislocation to be a continuously distributed object, by considering the Burgers vector to be distributed around the dislocation line.To describe this distribution, we use a scalar function φ(r(s), r ⊥ ) where r ⊥ • t(s) = 0 , i.e., r ⊥ is the position vector in the plane locally perpendicular to the dislocation line, with the dislocation line located at r ⊥ = 0 .The function φ is normal- ized, i.e. φdr ⊥ = 1 .The corresponding dislocation density tensor field is given by This tensor is linked to the plastic distortion β pl via So far, these are general relations used in the continuum theory of dislocations.The connection with peridynamics is established by noting that the plastic distortion represents a stress-free deformation which, in bond-based peridynamics, can be understood in terms of the force-free deformation of bonds.We parameterize a bond ξ connecting to x as r(ξ , x, s) = x + e ξ s, 0 ≤ s ≤ ξ .The bond vector in the plastically deformed configuration is then evaluated as with length ξ ′ = ξ ′ .The plastic bond deformation does not create a force but simply changes the bond vector and bond vector length entering the force calculation, which still is conducted according to Eqs. (3) or ( 7). ( 27) Once the modified bond vectors are computed, one more step is needed which amounts to a re-definition of the interaction sphere ('horizon').This needs now to be defined in the bond-deformed configuration: In other words, the plastic bond deformation will convect some material points out of the original sphere of influence, while some new ones enter.This is in line with the idea that plastic deformation implies a change of atomic neighborhood relations and atomic interactions, while the physical nature of the environment of an atom remains unchanged (here: the influence region remains spherical, hence material isotropy is preserved if deformation is homogeneous).

Singular edge dislocation
To illustrate the above ideas, we first consider the case of a singular dislocation which induces a discontinuity of slip along its glide plane.We consider a straight edge dislocation running along the z-axis with Burgers vector b = be x , defining a plane strain prob- lem in the perpendicular plane.The corresponding density function is φ(r ⊥ ) = δ(r ⊥ ) , where here δ(•) is the Dirac delta distribution.In the following, the subscript ⊥ is dropped since a 2D, plane strain implementation in peridynamics is envisaged, where only space dependence in the xy-plane needs to be considered.The dislocation tangent vector is t = e z , and the slip plane normal vector is n = e y .The singular dislocation den- sity tensor follows as and the plastic distortion field which fulfills α = ∇ × β pl is, up to a constant offset, given by where H is Heaviside's unit step function and the second, equivalent term on the righthand side has been added for later use.From the non analytic function given by Eq. ( 32), the bond displacement is calculated via Eq.( 29).The result can be stated in simple words: If the bond vector crosses the slipped parts of the slip plane (i.e. in 2D, the negative x-axis) in +y-direction, we add b , if it crosses the slip plane in −y-direction, we sub- tract b.

Edge dislocation with distributed Burgers vector
The concept of a singular point-like dislocation sits uneasily with the general spirit of peridynamics, where interactions are considered of finite range.A straightforward generalization is to consider a Burgers vector distribution φ of finite range.For illustration, we consider the same edge dislocation as in the previous section but now assume that the Burgers vector is evenly distributed over a circle of radius δ p around the origin: (30) where r = |r| .Alternatively, we consider a Burgers vector distribution which is 'designed' to match the field of a dislocation in gradient elasticity.To this end, the Green's function of the Helmholtz equation is used as Burgers vector distribution: with l H = δ p √ 3/5 .The distributed Burgers vector induces a distributed dislocation den- sity tensor and a spatially distributed plastic distortion which is not explicitly given but computed numerically, with results shown in Fig. 2. Again, the bond deformation follows by using Eq. ( 29).

Scaling property of the dislocation related fields
We now assume that the dislocation related fields can be computed using a small strain approximation where b ≪ δ p .If, moreover, the same length scale δ p governs the range of elastic interactions and the non-local distribution of strain, as in our above example, then the ensuing fields exhibit scaling properties which are well-known in the context of classical dislocation theory (Zaiser and Sandfeld 2014;Berdichevsky 2023).In particular, the stress field σ , plastic distortion β pl and elastic distortion β el all have the scaling form and the displacement field has the form Equation ( 35) holds for the singular distortion field given by Eq. ( 32) as demonstrated by the second term on the right hand side of that equation.Accordingly, the plastic bond displacements computed from Eq. ( 29) obey Eq. ( 36).Moreover, the elastic displacements that arise in response to the bond deformation must be proportional to the same owing to the linearity of Eq. ( 7).The scaling property for the stress field then follows by noting that, in the small strain limit, η ≈ ξ ≈ ξ ′ is used in the calculation of the virial stress via Eq.( 45).It is noted that the stress, distortion and displacement fields of a dislocation in classical elasticity fulfill Eqs. ( 35) and (36) for arbitrary δ p .Finally it is noted that, in the small strain limit and in an infinite body, the field of a distributed dislocation follows from that of a singular distribution by convolution with the Burgers vector distribution provided this field again has the above formulated scaling form.

Dislocations in gradient elasticity theory
The elasticity theory of dislocations in Mindlin's first gradient elasticity has been discussed extensively by Lazar (2021), considering both the general theory (Lazar 2021) and simplified versions with a single length scale ('Helmholtz-type gradient elasticity' (Lazar and (33) Fig. 2 Plastic distortion fields of distributed dislocations in peridynamics, Burgers vector points in x-direction and the glide plane trace corresponds to the x-axis, left graph: constant Burgers vector distribution given by Eq. ( 33); right graph: Burgers vector distribution given by Eq. ( 34) Maugin 2005; Lazar 2013)).Here, we focus on the case of edge dislocations which represents a plane strain problem.At first, the main results of the general theory, as applied to edge dislocations by Lazar (2021), are summarized.
In the full theory, the stress field of an edge dislocation is controlled by four length scales which relate to the gradient coefficients a 1 . . .a 5 , Eq. ( 26), via (Lazar 2021) The displacement, elastic and plastic distortion fields as well as the stress field of the dislocation can be expressed in terms of these four characteristic lengths in conjunction with the elastic constants of the material and the dislocation Burgers vector.Explicit expression for the stress field are given by Lazar (2021), considering an edge dislocation with Burgers vector b = be x : (37) (38) Here the K i are modified Bessel functions of the second kind.It is noted that, as in the classical theory, the sum σ xx + σ yy has the same functional form as the yz-stress field component of a screw dislocation.
As discussed by Lazar (2021), a simplified theory emerges if all four characteristic lengths l 1 . . .l 4 are equal, l 1 = l 2 = l 3 = l 4 = l H , leading to so-called gradient elasticity of Helm- holtz-type (Lazar and Maugin 2005;Lazar 2013).This theory has the advantage that, just as in the classical theory of dislocations, the stress fields of general dislocation configurations can be expressed analytically in terms of integrals over the dislocation lines.For Helmholtztype gradient elasticity, the relationships between the classical and gradient elasticity fields are exceedingly simple (Lazar and Maugin 2005): where is the Laplace operator, X denotes the field in gradient elasticity and X 0 its clas- sical counterpart.In other words, the gradient elasticity fields can be obtained by convolution of the classical fields with the Green's function of the Helmholtz operator, where K 0 is again a modified Bessel function of the second kind.
As to the scaling properties discussed in the previous section, we observe that the stresses given by Eq. ( 38) exhibit the same scaling properties as the peridynamic stress field.This can be seen by noting that, when the gradient elasticity is matched to peridynamics, according to Eqs. (26, 37), all four characteristic lengths of the gradient theory are proportional to δ p , and the scaling form of Eq. ( 35) is then, for the stress field, immediately evident from Eq. ( 38).As to the other fields, the scaling behavior of Eqs. ( 35) and (36) immediately follows from Eq. ( 39) by noting that the Helmholtz operator is scaling invariant if l H ∝ δ p and the classical fields exhibit the same scaling behavior as their peridynamic counterparts.

Simulation method
In our simulations, we consider a square system deforming in 2D plane strain.For discretization a regular square grid of 200 × 200 collocation points is used.Space is measured in units of the collocation point spacing.The peridynamic horizon is varied in the range 3.01 ≤ δ ≤ 12.01 .The origin of the coordinate system is set in the centre of the square which also defines the location of the dislocation, c.f. Fig. 1.
In this type of peridynamics implementation, Eq. ( 1) is discretized into a system of coupled equations for the mass-endowed collocation points akin to a molecular dynamics or molecular mechanics simulation.In the present study, we focus on the quasi-static properties of dislocations and accordingly, in the spirit of molecular mechanics, we consider Eq.
(1) in the static limit where the system is in force equilibrium such that the physical value of the mass contained in the volume associated with each collocation point does not matter.To find the force equilibrium configurations, we use energy minimization via an adaptive dynamic relaxation scheme as proposed by Kilic and Madenci (2010).To this end Eq.( 29) is first used to create the force-free bond deformation associated with the plastic distortion (39 due to the dislocation.Equation ( 1) is then solved for this disturbed system using an explicit leap-frog type time integration scheme and adding artificial damping.Mass, time step and damping coefficient are chosen to achieve optimal convergence.Exploiting Gershgorin's theorem provides an estimate for the optimal fictitious mass density matrix: where d ii are the diagonal components of the mass density matrix, V i the volume of col- location point i, t is the chosen relaxation time step, and k ij the ij-component of the stiffness matrix.The stiffness matrix for the system is obtained by numerical differentiation.The damping of the system is set to be proportional to the mass matrix of the system, with the damping factor c n adaptively computed for each iteration n by (Kilic and Madenci 2010) where u n is the displacement at iteration n and 1 k n an iteration dependent diagonal local stiffness matrix that is computed through (Kilic and Madenci 2010) Here un−1/2 i is the velocity at the intermediate step n − 1/2 and d ii the mass density of the collocation point i and r i is the corresponding force density vector including internal and external forces for the iteration n and n − 1 , respectively.The conver- gence criterion used, requires that the norm of an estimation for the displacement field for the next intermediate step in the explicit integration scheme, computed from the velocity field u for the iteration n + 1/2 multiplied with the time step t and scaled by the norm of the current displacement field u n+1 , is below a thresh- old which we set to 10 −8 (Sauvé and Metzger 1995), i.e.
We emphasize that this approach is not geared towards efficiently implementing dislocation dynamics, rather, we use it to reliably establish, in a first step, the properties of dislocations in the quasi-static regime.
For comparison with classical or gradient elasticity theories, it is useful to compute stress fields.Stresses are however not a natural concept in peridynamics which is based upon forces and displacements, rather than stresses and strains.Therefore, we use the analogy with particle systems where stresses are normally computed on a particle level using the virial theorem.Thus we evaluate the virial stress in the deformed configuration as a proxy of the Cauchy stress, (41) We note that omitting η would instead lead to the second Piola-Kirchhoff stress in the intermediate configuration.For this stress measure it was recently shown (Li et al. 2022) that it is equivalent to the peridynamic stress measure proposed by Silling and Lehoucq (2008).

Singular edge dislocation
We first consider the case of an edge dislocation with singular distribution of the plastic distortion, Eq. ( 32).Resulting stress fields are shown in Fig. 3 (top) for different values of the horizon δ p as indicated in the graphs, which show sections of the σ xy shear stress component along the x and of the σ xx component along the y-axis.It is noted that the functional shape of σ yy (x, y) is almost identical with σ xy (y, x) , thus no separate plot is shown for this component.
Different from the classical solution, the peridynamic stress fields do not exhibit a singularity in the origin.The deviation from the classical solution is evident at distances of the order of δ p from the singularity (Fig. 3, top).At the same time, there is a kind of convergence to the classical solution as decreasing the horizon δ p leads to an increase of the maximum stress near the dislocation and to a decrease of the distance over which differences between the peridynamic and classical solutions are evident.This qualitative (45) observation was reported previously by Zhao et al. (2021), who however fail to realize the underlying mathematical scaling structure.This scaling structure is illustrated in Fig. 3 (bottom), where components of the stress tensor, multiplied with δ p , vs. the re-scaled space coordinates x/δ p and y/δ p are plotted.According to Eq. ( 35), the result should be independent of δ p , which is approximately correct for δ p ≥ 9 nm (9 collocation point spacings).This result depends only on the number of collocation points within the horizon, thus, deviations at smaller values result from the fact that the number of collocation points within the horizon is too small.It is noted in this context, that the actual solution is difficult to represent because it is non analytic in the origin where both stress components σ xx and σ xy exhibit a jump, as can be seen from Fig. 3. Accordingly, a finite spacing of the collocation points used for stress evaluation may lead to the erroneous suggestion of smooth behavior in the origin, a problem that may be exacerbated by injudicious use of interpolation functions.The angular dependence of the stress fields matches their classical counterparts as illustrated in Fig. 4. We next analyze to which extent the observed behavior is contingent on the choice of the micro-modulus function.
Figure 5 compares results obtained for constant and 'conical' micro-modulus distribution, all other parameters being identical.It can be seen in Fig. 5, left, that the 'conical' distribution, which gives a higher weight to short bonds and leads to a smaller length parameter l, results in more pronounced positive and negative stress peaks located at the collocation points adjacent to the origin.However, if we exploit the scaling invariance expected for both gradient elasticity and peridynamics simulations, and plot lσ xx vs. r/l where l = δ p √ 3/5 for constant and l = δ p √ 2/5 for 'conical' micro-modulus dis- tribution, one finds a near-perfect collapse of both curves as seen in Fig. 5, right.This indicates, on the one hand, that the non analytic behavior in the origin is not just an artefact of a particular choice of micro-modulus distribution, and demonstrates on the other hand that dislocation fields obtained from different distributions can be analyzed in a common mathematical framework which measures the range of non-local forces in terms of the length scale l as given by Eq. ( 22).

Distributed edge dislocation
Next, the stress fields of dislocations evaluated with distributed Burgers vector are considered.All simulations reported in the following use a micro-modulus that is constant over the influence sphere.First, the case where the Burgers vector is evenly distributed over the influence sphere is considered, Eq. ( 33).
Figure 6 (left) gives profiles of the σ xx component along the y-axis, evaluated for dif- ferent values of δ p ∈ {3, 6, 9, 12} nm.Comparison with Fig. 3 demonstrates that the peak stress values are reduced as compared to the singular dislocation, and that the non analytic behavior in the origin disappears.The scaling collapse in Fig. 6 (right) now indicates excellent superposition of the curves for all δ p values, indicating much improved convergence of the dislocation fields towards the scaling solution, even for small numbers of collocation points within the horizon.This is important from the viewpoint of efficient numerical implementation since, for fixed physical values of the system size and the horizon, the computation time increases in proportion to the square of the number of collocation points within the horizon.

Comparison with dislocation in gradient elasticity
To compare with gradient elasticity, a Burgers vector distribution is used that matches the one underlying the gradient elasticity solution, Eq. ( 34).Results obtained with this Burgers vector distribution are shown for δ p = 3 nm in Fig. 7.
For comparison, two variants of gradient elasticity are shown.First, the solution of the full Mindlin equations as given by Eq. ( 38) is used, where the four length parameters l 1 . . .l 4 were obtained from the value δ p = 3 nm using Eqs.(26, 37) in conjunction with Fig. 5 Effect of micro-modulus distribution on stress of a singular dislocation; all parameters as in Fig. 3; left: σ xx (x = 0, y) for δ p = 12nm, simulated using a constant micro-modulus and a 'conical' micro-modulus distribution; right: scaling collapse obtained by plotting lσ xx vs. y/l where l = δ p √ 3/5 for data evaluated with constant and l = δ p √ 2/5 for data evaluated with 'conical' micro-modulus l = δ p √ 3/5 for the constant micro-modulus distribution that was used in the peridy- namics calculations, to obtain l 1 = 3.35 nm, l 2 = 2.60 nm, l 3 = 3.07 nm, and l 4 = 3.67 nm.This solution is represented by the blue lines in Fig. 7.A simplified calculation is also performed, setting l 1 = l 2 = l 3 = l 4 = l H = 3 nm to obtain the Helmholtz gradi- ent elasticity solution, which is shown in red in Fig. 7.Both gradient elasticity solutions are quite similar, and both are in good agreement with the peridynamic solution in the ; the open circles indicate virial stresses evaluated at the respective collocation points, the red line gives the matching solution evaluated using Helmholtz gradient elasticity, and the blue line the solution evaluated using (full) Mindlin gradient elasticity central part of the simulated sample.Agreement is slightly better for the Helmholtz gradient elasticity solution, which is to be expected since the Burgers vector distribution has been chosen to match the Helmholtz gradient elasticity result.Different from previous figures, where we have focused on the region near the dislocation, Fig. 7 shows the stress tensor components for the entire domain from surface to surface.While the gradient elasticity results represent solutions in an infinite body which go to zero only asymptotically, the peridynamic solution automatically accounts for surface effects.Thus, near the surfaces, both gradient elasticity solutions differ from the peridynamic solution: On the upper and lower surfaces with normal vector e y , the hydrostatic stress component σ yy is approximately zero for the peridynamic dislocation, as is the shear stress σ xy on the left and right surfaces with normal vector e x .On the other hand, the hydrostatic stress component σ xx is non zero on the upper and lower surface.This behavior is in line with the surface boundary conditions imposed in classical elasticity theory on a free surface.The gradient elasticity stress fields, computed for a dislocation in an infinite body, assume finite values on all surfaces.

Construction of a hybrid multiscale dislocation model
It has been shown that the fields associated with a dislocation in peridynamics have certain scaling properties, and that they can be matched with results for dislocations in gradient elasticity if an appropriate Burgers vector distribution is chosen.Moreover, the peridynamic dislocation satisfies surface boundary conditions for the stress fields while gradient elasticity solutions calculated for a bulk dislocation do not.The question arises what these findings are good for.
First it is noted that peridynamics as such is not a tool of choice for dislocation simulations.If the horizon δ p is taken much larger than the atomic spacing, as it should for reasons of computational cost, then the stress fields close to the dislocation core are underestimated by a factor of the order of b/δ p and accordingly the interactions of close dislocations are ill represented.These interactions are however crucial for our understanding of work hardening (Kubin et al. 2008); the problem is well known from attempts to describe dislocation interactions by FEM in so-called discrete-continuum schemes (Lemarchand et al. 2001) where the element size plays a similar role.If, on the other hand, the horizon is taken of the order of the atomic spacing, then the simulations become computationally very expensive and molecular dynamics is a both cheaper and more accurate option.However, the close match between numerical solutions obtained by peridynamics and analytical solutions obtained by gradient elasticity still allows us to devise a useful hybrid scheme.
The method is based on using a combination of peridynamic and gradient elasticity for evaluating the stress field of a dislocation in a multiscale model as the sum of three contributions: • On the peridynamic side, a horizon δ p is chosen which is computationally conveni- ent and may be much larger than the physical extension of the dislocation core.With this horizon the methodology described above is used to construct a peridynamic dislocation with desired line direction and Burgers vector, with the Burgers vector distribution given by Eq. ( 34).The resulting peridynamic stress field is denoted as σ p r, δ p .It is noted that, since this stress field is free of singularities or disconti- nuities near the dislocation core, it can be safely interpolated between the collocation points to evaluate stresses at any desired spatial resolution.• Using the relations of (Helmholtz) gradient elasticity, matching gradient elasticity fields with length scale parameter l H ∼ δ p are evaluated.The corresponding stress field is denoted as σ g,δ p (r) .Near the dislocation core, it has the same behavior as the peridynamic stress field, whereas at distances r ≫ δ p from the dislocation it behaves like the classical bulk solution for the dislocation stress.• Again using Helmholtz gradient elasticity, a gradient elasticity field with a length scale parameter l b ∼ b is evaluated.The parameter l b is chosen such to match the core properties of the dislocation on the atomistic level as determined e.g. by molecular dynamics.The corresponding length scale parameter will in general be proportional to, but not identical with the Burgers vector length b.The corresponding stress field is denoted as σ g,b (r) .As δ p ≫ l b , this field will typically exhibit much larger val- ues near the dislocation core, while its behavior for r ≫ δ p is the same as for the field σ g,δ p (r).
The global solution is now constructed by simple superposition of these three fields in the form Thus, the peridynamic stress field is corrected by a stress field evaluated using gradient elasticity, which is the difference between the fields σ g,b (r) and σ g,δ p (r) .The rationale behind this procedure is as follows: (i) As can be seen from Eq. ( 38), taking the difference of two gradient elasticity stress fields with different length scales eliminates the slowly decaying classical elasticity solution.The difference is composed of modified Bessel functions only which decay exponentially as functions of r/b and r/δ p .As a conse- quence, on scales well above the horizon δ p , the field σ g is exponentially small and can be neglected.Thus, unless the dislocation is within about one horizon from the surface, the near-surface stress field is dominated by the peridynamic solution which correctly accounts for the surface boundary conditions.(ii) Within a distance of the order of δ p from the dislocation core, the fields σ p and σ g,δ p are practically identical (see Fig. 7) and their difference is negligibly small.The overall stress field is thus completely controlled by σ g,b , which has the correct behavior at the dislocation core.The construction is illus- trated in Fig. 8, showing in double-logarithmic scale the behavior of the σ xy shear stress component on the positive x-axis.On the left we see the three component fields with their different asymptotic behavior at large and small x, on the right the peridynamic stress field and the overall gradient correction together with the hybrid stress field that arises from the superposition.The angular dependence of the gradient correction field σ g (r) is illustrated in Fig. 9.It is clearly seen that the correction field becomes negligibly small above the scale of a few horizons of the peridynamic field (here: δ p = 3 nm).While the overall morphology matches that of the classical solution, it is noted that the angular dependency of the correction field exhibits some deviations from the classical expectation, for instance, the location of the zero-field level lines in Fig. 9 indicates that the ( 46

Comparison with approaches of other authors
There are a few recent papers where peridynamics has been used for evaluating dislocation stress fields in finite samples.In generalization of the classical approach of Needleman and Van der Giessen (1995) where analytically obtained bulk stress fields of dislocations are corrected by the FEM solution of a modified elastic boundary value problem, Dong et al. (2022) replace the finite elements by a peridynamic simulation, using a layer of virtual collocation points outside the surface to impose the required surface boundary conditions.This framework differs from ours inasmuch as the peridynamic grid experiences no plastic distortion.On the other hand, Zhao and Shen (2021) implemented a dislocation within a state-based peridynamic framework by imposing a discontinuous, force-free bond deformation along the sheared part of the slip plane.Their dislocation model can be considered an implementation of our Eq.( 29) with the plastic distortion given by the singular expression, Eq. (32).At the same time, their state based peridynamic framework has the advantage that it can be generalized to describe dislocations in anisotropic materials.However, similar to other models that rely exclusively on continuum methodology to numerically calculate the stress fields of discrete dislocations, the model of Zhao and Shen suffers from the generic drawback that the core size is controlled by the characteristic scale inherent in the model (for peridynamics: the horizon).This poses the conundrum that either the stresses near the dislocation core are grossly under-estimated or the horizon must be taken to be on the scale of the Burgers vector length b, which would incur a prohibitive computational cost.
Our generalization of the formulation of Shen and Zhao to dislocations with distributed Burgers vector allows us to combine peridynamics and gradient elasticity in order to construct a two-scale model of a dislocation which uses peridynamics for evaluating dislocation fields that correctly account for surface boundary conditions, and combines this evaluation with a short-range correction to account for the high stress fields close to the dislocation.Our two-scale model is similar in spirit to other multiscale dislocation models proposed in the literature, all of which use a decomposition of the dislocation stress field analogous to our Eq.46: into a system-scale long-range stress field, which is computed numerically using a continuum model (in our case the peridynamic stress field σ p ), and short-range corrections.In the multiscale approaches of Jamond et al. (2016) and of Huang (2015), the system-scale stresses are evaluated from a coarse grained dislocation eigenstrain using a standard elastic-plastic FEM framework suitable for dealing with boundary value problems in finite samples.Bertin amd co-workers (Bertin and Capolungo 2018;Bertin 2019) use a framework based on the coarse grained dislocation density tensor (field dislocation mechanics, FDM) in conjunction with spectral methods (FFT), which makes their approach particularly suited for dealing with bulk samples using periodic boundary conditions.Both Jamond et al. and Bertin et al. construct the short-range corrections on the level of dislocation segments, using a formalism based on a distributed Burgers vector introduced by Cai et al. (2006), as the difference of two fields with different radii of the Burgers vector distribution.This procedure is very similar to our Eq.( 46), and has the same rationale.Our approach differs in the numerical method used (peridynamics vs. FFT/FDM and FEM) and also in the analytical framework used for constructing the short-range correction (gradient elasticity vs. classical elasticity with distributed Burgers vector).When comparing the respective advantages and disadvantages, we note that peridynamics is numerically more expensive than FEM and even more so than FFT/FDM, but offers important advantages when dealing with problems involving non stationary cracks and damage processes.At the same time, we note that the computational cost is still manageable: a quasi-static or dynamic 3D simulation with 200 3 collocation points and a horizon of 3 collocation point spacings (a comparatively small value which is however acceptable for our multiscale scheme), would be comparable in numerical effort to a molecular mechanics or molecular dynamics simulation with about 10 7 atoms using a standard EAM potential, respectively.We further note that the correction fields obtained via Eq.( 46) from gradient elasticity are mathematically expressed as differences of Bessel functions and thus decay exponentially, whereas the formalism of Cai et al. leads to an algebraic decay of the correction function.This is important in practice, since a cut-off to the short-range interactions is required in order to make the problem of their evaluation of order N, N being the number of dislocations/dislocation segments in the system.In short, using gradient elasticity for the short-range correction may reduce potential hiccups associated with cutting it off at a finite, and preferentially short, distance from the dislocations.

Perspectives and outlook
The presented results provide all necessary 'ingredients' for constructing a plasticity model.Calculating the evolution of the plastic distortion is straightforward once the dislocation velocity is known: For the geometry of the edge dislocation, the plastic distortion rate fulfils the relation βpl = γ e x ⊗ e y , where the scalar slip rate is related to the Burgers vector distribution via γ = bvφ(r) where v is the scalar dislocation glide veloc- ity.The temporal change of the plastic distortion gives, in turn, according to Eq. ( 29) the evolution of the bond structure due to the plastic deformation.The driving forces that make the dislocation move and cause multiple dislocations to interact are straightforward to evaluate since, as shown by Lazar et al. (2007;2014), in Helmholtz-type gradient elasticity the Peach-Koehler force retains its classical form in the sense that it involves only the Cauchy stress acting on the dislocation line -which we have shown how to compute.Moreover, if Helmholtz gradient elasticity is used to compute the short-distance correction fields, the approach is straightforward to generalize to threedimensional systems of curved dislocation lines.Such an approach may, in conjunction with the well known advantages of peridynamics in dealing with incipient cracks and crack branching, be of particular interest for studying dislocation-crack interactions in a dynamic setting.
There are limitations.For a dislocation within distance δ p or less from the surface, the computation becomes unreliable because the correction stress σ g (r) does not account for the surface boundary conditions; as a consequence, dislocation interactions with the surface are under-estimated at short distances from the surface.This problem exists, in one form or another, in all approaches which use a continuum simulation method to provide surface corrections to bulk dislocation fields, e.g., in the classical approach of Needleman and Van der Giessen (1995) where bulk stress fields are corrected by the FEM solution of a modified surface boundary value problem, the problem becomes manifest on scales below the resolution of the finite element mesh.Similarly, in the recently proposed method of Dong et al. (2022), where the finite element mesh is replaced by a peridynamic simulation approach, inaccuracies arise once the dislocationsurface spacing gets on the order of the peridynamic horizon.How to correct potential near-surface artefacts remains an important question for further work on developing the present approach.

Fig. 1
Fig.1Schematic of a single edge dislocation with slip plane normal vector n.The origin of the coordinates is centred in the considered domain and the grey area indicates the domain where the plastic distortion is provided by the gradient theory solution; to match peridynamics and Helmholtz type gradient elasticity, the horizon δ p must be taken proportional to the length scale parameter l H of the gradient theory

Fig. 3
Fig. 3 Stress field components of a singular dislocation in peridynamics, the Burgers vector points in x-direction and the glide plane trace corresponds to the x-axis; left: σ xy (x, y = 0) , right: σ xx (x = 0, y) ; the top graphs show the original results for different values of the horizon where δ p = 3nm, 6nm, 9nm and 12nm (marked H3, 6, 9, 12 in the graph labels); the bottom graphs show the same data after rescaling σ → δ p σ , r → r/δ p ; in these graphs the locations of the collocation points at which virial stresses are evaluated are marked by open circles

Fig. 6
Fig. 6 Stress field components of a distributed dislocation in peridynamics, the Burgers vector points in x-direction and the glide plane trace corresponds to the x-axis, Burgers vector distribution given by Eq. (33); left: σ xy (x = 0, y) for different values of the horizon, right: same data after rescaling σ → δ p σ , r → r/δ p ; in these graphs the locations of the collocation points at which virial stresses are calculated are marked by open circles

Fig. 8
Fig. 8 Construction of a hybrid dislocation model, Burgers vector b = 0.3nm, horizon of the peridynamic solution δ p = 3nm, the length scale of the 'inner' gradient solution σ g,b is l b = 0.3nm, all

Fig. 9
Fig. 9 Angular dependence of the short-range correction function σ g as evaluated for the parameters of Fig. 8; left: stress component σ g xx , right: stress component σ g xy