The Peierls--Nabarro FE model in two-phase microstructures -- a comparison with atomistics

This paper evaluates qualitatively as well as quantitatively the accuracy of a recently proposed Peierls--Nabarro Finite Element (PN-FE) model by a direct comparison with an equivalent molecular statics simulation. To this end, a 2D microstructural specimen subjected to simple shear is considered, consisting of a central soft phase flanked by two hard-phase regions. A hexagonal atomic structure with equal lattice spacing is adopted, the interactions of which are described by the Lennard--Jones potential with phase specific depths of its energy well. During loading, edge dislocation dipoles centred in the soft phase are introduced, which progress towards the phase boundary, where they pile up. Under a sufficiently high external shear load, the leading dislocation is eventually transmitted into the harder phase. The homogenized PN-FE model is calibrated by an atomistic model in terms of effective elasticity constants and glide plane properties as obtained from simple uniform deformations. To study the influence of different formulations of the glide plane potential, multiple approaches are employed, ranging from a simple sinusoidal function of the tangential disregistry to a complex model that couples the influence of the tangential and the normal disregistry. The obtained results show that qualitatively the dislocation structure, displacement and strain fields, and the dislocation evolution are captured adequately. The simplifications of the PN-FE model lead, however, to some discrepancies within the dislocation core structure. This plays a dominant role in the dislocation transmission process, which thus cannot quantitatively be captured properly. Despite its simplicity, the PN-FE model proves to be an elegant tool for a qualitative study of edge dislocation behaviour in two-phase microstructures, including dislocation transmission, although it may not be quantitatively predictive.


Introduction
Over the past decades, the Peierls-Nabarro (PN) model [2][3][4] has gained popularity in the dislocation community due to its ability to model dislocations on the atomistic scale while using a continuum framework. In its classical form, an infinite and homogeneous crystal is split into two linear-elastic regions connected by a glide plane. Along this glide plane, a relative tangential displacement (or tangential disregistry) is allowed, which is mapped onto an intrinsic misfit energy. In the classical PN model, the adopted glide plane potential is a periodic function, which is based on the Frenkel sinusoidal [5]. Minimising the total free energy, consisting of the elastic strain energy of the bulk and the misfit energy of the glide plane, leads to an arctan-type dislocation disregistry profile with satisfactory properties in terms of non-singular stress field, total energy, and Peierls energy as well as Peierls stress.
The PN model furthermore excels through its versatility when solved numerically. In this fashion, the limits of the classical approach can be breached. For instance, the influence of heterogeneous crystals on dislocation obstruction at phase boundaries, as well as on dislocation transmission, can readily be modelled, see, e.g., [6,7,1,8]. It has, however, been pointed out by some authors that, with the utilisation of the simple sinusoidal, the PN model is limited in its quantitative description of dislocations in real crystals [9]. A variety of extensions has therefore been suggested in the literature, as follows.
A significant increase in the accuracy of the misfit energy compared to the simple 1D sinusoidal can be achieved through the Generalized Stacking Fault Energy (GSFE) surface, introduced initially by Vítek [10]. It provides a full 2D-energy landscape (in contrast to the 1D sinusoidal), intrinsic to the crystal considered. With this extension, the PN model is capable of describing mixed dislocations, splitting of dislocations into Shockley partials, and the related recombination energy of partials, see, e.g., [9,[11][12][13]. The GSFE surface is commonly obtained from Molecular Statics (MS) calculations or directly from density functional theory (DFT) [14,15]. It has been shown recently that for a bilayer system, the solutions of the PN model with the GSFE extension are asymptotically close to the full atomistic model [16]. Yet, the GSFE is limited to tangential (i.e., in-plane) disregistry only. Lifting this constraint and introducing an additional normal (i.e., out-of-plane) disregistrydependency improves the dislocation description even further [17][18][19]. As the misfit energy resides between the atoms located above and below the glide plane (recall that the glide plane often reduces to a zero-thickness interface), it has been suggested to subtract the linear elastic part of the misfit energy [17,20]. This correction leads to a larger activation energy (as discussed in [21]), and to a decrease of the Peierls stress, cf. [22]. The inclusion of elastic anisotropy is an important feature for dislocations in anisotropic crystals [19,[23][24][25]. Further PN model extensions include non-local formulations reflecting the discreteness of the underlying atomic lattice [18,[26][27][28][29], or an additional gradient energy term [30] motivated by DFT calculations [31].
In previous research, various authors have studied the accuracy and predictability of the PN model in comparison with atomistic simulations. Von Sydow et al. [32] analysed Shockley partial dislocations in Pd and found a relatively good agreement of the PN model. More recently, Dai et al. [33] unveiled a high accuracy of the generalized PN model in capturing the structure and energy of low-angle grain boundaries and near-twin grain boundaries for (111) twist boundaries in Al, Cu and Ni. Mianroodi et al. [34] compared the generalised PN model of Xiang et al. [19], the phase-field dislocation dynamics model of Hunter et al. [35,36], the phase-field based models of Shen and Wang [37], and Mianroodi and Svendsen [38] with MS simulations and found a rather good agreement in terms of the dissociation of a dislocation dipole in fcc single crystals of Al and Au. Considering a non-local formulation of the PN model, Liu et al. [29] presented a good predictive capability for the dislocation core structure and the Peierls stress in Fe and Cu. Up to this point, all of the previous studies have been based on single crystal systems. For a two-phase microstructure, however, no such comparison has been performed yet.
This paper aims to partially bridge this gap, by the qualitative and quantitative comparison of the predictions made by the recently proposed Peierls-Nabarro finite element (PN-FE) model [1] and by MS simulations. The focus of this comparison is on the pile-up evolution of edge dislocations in a two-phase microstructure under an increasing external shear load and the eventual dislocation transmission across the phase boundary. The considered problem is limited to a 2D lattice to reduce the complexity of the 3D reality and to facilitate a sharp comparison. Different complexities of the glide plane potential are considered, ranging from a simple sinusoidal up to a fully coupled potential as a function of the tangential disregistry ∆ t as well as the normal disregistry ∆ n . Their performance, and predictive capabilities, will be assessed critically.
The problem considered employs a 2D microstructure consisting of two phases. A soft central phase (Phase A) is flanked by two hard-phase regions (Phase B). Both phases have a hexagonal atomic lattice structure with homogeneous spacing. The specimen is subjected to a shear deformation applied through the external boundary, which induces in a defectfree configuration a state of uniform shear stress. Edge dislocation dipoles are generated in the soft Phase A on a glide plane perpendicular to the coherent and non-damaging phase boundaries. Due to the phase contrast, the dislocations pile up at the phase boundary where eventually, under a sufficiently high external shear load, the leading dislocation is transmitted into the harder Phase B.
The remainder of this paper is divided into four sections as follows. In Section 2, the considered problem is described in detail and the basics of the MS and PN-FE models are briefly recalled. The different glide plane potentials employed for the purpose of comparison are subsequently introduced in Section 3. Individual solutions of the PN-FE model for a single dislocation dipole and a dipole pile-up are compared with the results of the MS model in Section 4. The paper closes with a summary and discussion in Section 5.

Problem statement
The geometry of the test example employed throughout this manuscript is first specified in Section 2.1, whereas basics of the atomistic MS calculations along with prescribed boundary conditions and used potentials are detailed in Section 2.2. The definition and numerical values of effective quantities required for the homogenized PN-FE model are specified in Section 2.3 in terms of constitutive linear-elastic constants and glide plane potential. The continuum PN-FE framework itself is finally described in Section 2.4.

Atomistic system
A two-phase microstructure, as illustrated in Figure 1, is considered. It consists of a soft Phase A that is flanked by a harder Phase B. Both phases have an identically oriented hexagonal lattice structure of equal spacing. The crystal orientation is chosen such that one set of glide planes Γ gp is oriented perpendicular to the coherent phase boundaries Γ pb . Through the external boundary, a shear deformation is induced, which corresponds in a defect-free and linear-elastic configuration to the state of homogeneous shear stress τ . Stable dislocation dipoles are initialised in the centre of the Phase A, which under increasing shear load move towards the phase boundary. Due to the phase contrast, the dislocations get obstructed and gradually pile up. Eventually, under a sufficiently high external shear load, the leading dislocation is transmitted into the harder Phase B. Figure 1: A sketch of the atomistic representation of a two-phase microstructure used for simulations of edge dislocation dipoles interacting with a coherent phase boundary. The entire domain (Ω, all atoms) consists of two phases, Phase A (Ω A , grey atoms) and Phase B (Ω B , black atoms), and of the boundary domain (Ω BC , white atoms). Positions of all atoms situated inside Ω BC are prescribed to induce shear deformation.

Molecular statics problem
Full discrete MS simulations are performed, which serve two purposes. First, they provide the reference solution of the above-specified problem to which the PN-FE model is compared in Section 4. Second, atomistic simulations of simple, uniform, single-phase problems provide the homogenised material properties of the corresponding continuum PN-FE model such as linear elasticity constants and glide plane properties.
For studying edge dislocation vs. phase boundary interaction, the problem domain Ω = Ω A ∪ Ω B is considered, as sketched in Figure 1, with the regions occupied by the two Phases A and B specified as where x = x e x + y e y is the position vector. The hexagonal lattice with spacing a 0 is initiated such that a stress free state results (i.e., σ = 0, cf. Eq. (11) below). Each atom α located inside the domain Ω is stored in an index set N = N A ∪ N B , where and r α 0 = r α 0x e x + r α 0y e y denotes the spatial position of atom α in the reference configuration. The mechanical behaviour of the system is governed by its total potential energy V, defined as the sum of all pair potentials φ αβ (r αβ ), where r αβ = r αβ 2 denotes the Euclidean distance between a pair of atoms α and β. A shifted Lennard-Jones potential is employed, with the cutoff radius r cut . Standard notation is adopted for the unshifted Lennard-Jones potential, i.e., the well-depth is ε whereas the distance to the potential minimum is r m ; for further details see, e.g., [39]. Introducing a material contrast ratio ρ, the pair potential of Phase B is set to φ BB = ρφ AA . The constitutive parameters associated with the individual phases are summarised in Table 1. The interaction potential across the phase boundary is considered as the average of that of the two individual phases, i.e., φ AB = 1 2 (φ AA + φ BB ).
The prescribed boundary conditions applied on Ω BC are chosen such that in a defect-free and linear-elastic configuration a state of constant shear stress τ would result. Assuming a phase-wise homogeneous and linear-elastic response, the corresponding shear strains read tτ /µ A in Phase A and tτ /ρµ A in Phase B. Here, τ is the target shear load, µ A the homogenised shear modulus of the Phase A (cf. Section 2.3), and t ∈ [0, 1] a pseudo-time parametrising the evolution of the otherwise rate-independent system (i.e., zero temperature is assumed). In order to introduce the shear deformation, a layer of atoms N BC ⊂ N is displaced accordingly. The thickness of the layer N BC is chosen in accordance with the cutoff radius r cut to eliminate any surface effects. Atom positions are prescribed as In order to predict the mechanical behaviour of the specimen, the total potential energy is minimised at each time step t k , The configuration space at a time step t k is denoted as Q k ⊆ R 2nato , and reflects any prescribed atom displacements (recall Eqs. (5) and (6)). The minimisation problem is solved using the Trust-region methodology (cf., e.g., [39]), which has been implemented within an in-house code. Because multiple glide planes may be activated due to random perturbations and roundoff errors, the numerical solver is initialised towards the preferred glide plane Γ gp = { x ∈ R 2 : y = 0}. In particular, the current equilibrium position x = ± e x of the dislocation dipole is estimated by means of the analytical stress-field of a Volterra dislocation in front of a phase boundary, cf. [40]. Considering all possible dislocation dipole positions, an unstable and a stable equilibrium position may be distinguished. The former relates to the position x where an infinitesimal perturbation of the position leads to either the annihilation of the dipole or the propagation towards the phase boundary into the stable equilibrium position. The latter is hence the position x in which the dislocation is introduced. A scaling of the Volterra displacement field (see, e.g., [41]) is used to perturb the current relaxed configuration at a time step t k , i.e., This perturbed configuration is used as an initial guess for the minimization algorithm in a consecutive time step t k+1 . In Eq. (8), C m is the magnitude of the perturbation (chosen as C m = a 0 /2), whereas ν is the Poisson's ratio (obtained as a part of the homogenized lattice properties of Section 2.3).

Homogenization of the atomistic properties
The PN-FE model necessitates effective material properties matching those of the atomistic system. They comprise the linear-elastic behaviour and the glide plane potential. To reconcile the atomistic model and the PN-FE model under plane strain condition, a virtual thickness of r A m is considered for the atomistic model. First, the linear-elastic constants are extracted from the homogenised stiffness tensor 4 C, following the definition (see, e.g., [41]) where V = A · r A m denotes the (virtual) volume of the simulation cell in the deformed configuration with the according in-plane sectional area A, and r αβ i the components of the r αβ vectors. The (isotropic) stiffness tensor is obtained for a doubly periodic simulation cell of The obtained parameters, computed for Phase A (cf. Table 1), are summarised in Table 2. The parameters of Phase B are obtained by simply multiplying those of Phase A with the material contrast ratio ρ, i.e., 4 C B = ρ 4 C A . For later reference, the definition of the Virial stress σ is also recalled: A second set of homogenised quantities specifies the glide plane behaviour. The glide plane is characterised by the potential ψ * gp that represents the misfit energy between two rigidly shifted bulks of atoms, as shown in Figure 2a. ψ * gp is computed by considering a lattice with the stress free spacing a 0 in a simulation box (again of size 20a 0 × 12a 0 √ 3). The lattice is first equilibrated with periodic boundary conditions applied in the x-direction while the top and bottom edge are left free. The glide plane potential is subsequently computed by shifting the upper part of the lattice rigidly according to ∆ = ∆ t e t + ∆ n e n , where ∆ t and ∆ n denote the tangential and normal disregistry (relative displacements) across the considered interface. These disregistry components are varied in the ranges 0 ≤ ∆ t ≤ a 0 and −0.2 a 0 ≤ ∆ n ≤ 2 a 0 . For each of these states the energy of the system is evaluated and the energy of the initial stress free system is subtracted. The obtained result for Phase A is shown in Figure 2b. Reference values of ψ * gp are the work of separation G c = lim ∆n→∞ ψ * gp , the unrelaxed unstable stacking fault energy γ is the normal disregistry of zero normal traction T n (∆ t ) = ∂ψ * gp (∆ t )/∂∆ n = 0, and the characteristic length l c , i.e., the normal disregistry ∆ n at ∆ t = 0 where ∂ 2 ψ * gp (∆ t )/∂∆ 2 n = 0. The values obtained for the landscape of Figure 2b are listed in Table 3. In analogy to the elasticity parameters, the potential of Phase B is scaled to the potential of the Phase A by the material contrast ratio ρ, i.e., ψ * B gp = ρψ * A gp .   Figure 2: (a) Sketch of the applied displacements for the computation of the glide plane potential ψ * gp (∆ n , ∆ t ); (b) shape of the glide plane potential corresponding to Phase A, i.e., ψ * A gp (∆ n , ∆ t ).

Peierls-Nabarro finite element (PN-FE) model
The atomistic problem as specified in Figure 1 is translated to the PN-FE framework, introduced in [1], in a straightforward manner. Let Ω be the two-phase microstructure under consideration, cf. Figure 3. Similar to the atomistic model, it consists of two regions Ω A and Ω B , that are separated by perfectly bonded and coherent phase boundaries Γ pb normal to e x , i.e., In Ω lies a single glide plane Γ i gp , normal to e y , that splits each phase Ω i , i ∈ {A, B}, into two subdomains Ω i ± as follows: Assuming that the (non-linear) dislocation-related deformation is considered only within the glide plane, the total free energy (per unit thickness) consists of two contributions of the elastic strain energy density ψ e (considered inside Ω i ± ) and of the glide plane potential ψ gp (non-zero along Γ i gp ): Here ψ e is calculated based on linear elasticity under a plane strain condition and with the phase specific fourth-order elasticity tensor 4 C i , which is employed according to Section 2.3. The elastic strain tensor is The linear elastic stress tensor follows accordingly with The glide plane potential ψ gp is a function of the displacement jump (or disregistry) between Ω i + and Ω i − , expressed as ψ gp is a periodic function in ∆ t to capture the effect of lattice periodicity. The different glide plane potentials, which are employed in this paper, are introduced and calibrated in Section 3. Note that the elasticity and glide plane properties are homogeneous in each phase. This results in a jump in material properties across Γ pb , in contrast to the atomistic model in which the interaction between both phases results in a smooth transition. The perfectly bonded and fully coherent phase boundary Γ pb is modelled by enforcing displacement and traction continuity, i.e., Similarly to the atomistic model, the shear deformation is applied on the external boundary ∂Ω in the form of Dirichlet boundary conditions following Eqs. (5) and (6). Under such constraints, the mechanical equilibrium inside Ω is established by minimising the total potential energy Ψ of Eq. (18) subject to the phase boundary constraints of Eq. (24). Due to the periodicity of ψ gp this formulation results in a non-convex minimisation problem, which is solved using the finite element method by discretising Ω i ± and Γ i gp in space. The resulting total potential energy is minimised with the help of the truncated Newton optimization algorithm, see, e.g., [42]. Individual dislocations are initialized using an approach similar to the one described at the end of Section 2.2, where in addition a condition on the resolved shear stress around x = 0, which needs to exceed the critical nucleation stress τ nuc before a dislocation is introduced, is added. A more detailed description on the dislocation introduction can be found in [1].

Glide plane models and their calibration
To study the influence of the glide plane potential ψ gp on the accuracy of the PN-FE framework outlined above in Section 2.4, different complexities of ψ gp are considered in this section. The calibration of each potential is discussed for Phase A only, the properties of the Phase B potential follow accordingly using the contrast ratio ρ. First, the classical sinusoidal PN potential, which depends only on the tangential opening ∆ t , is described in Section 3.1.
The accuracy of this model can be improved by considering a Fourier series, as described in Section 3.2. A coupling of the tangential ∆ t as well as the normal ∆ n disregistry is introduced through an analytical model outlined in Section 3.3, whereas the direct use of the numerical data obtained from the MS calculations (i.e., ψ * gp (∆ n , ∆ t ), recall Figure 2b) is described in Section 3.4.

Classical Peierls-Nabarro potential
The standard choice for the glide plane potential is based on the Frenkel sinusoidal as a function of the tangential disregistry ∆ t , and is expressed as where γ us is the unstable stacking fault energy and b the Burgers vector magnitude. The normal disregistry is constrained to zero, i.e., ∆ n = 0. To calibrate this potential to the MS simulations, b is chosen to be equal to the atomic spacing a 0 to match the lattice periodicity. The unstable stacking fault energy, γ us , is determined from the atomically calculated energy landscape ψ * gp (∆ n , ∆ t ), recall Section 2.3. Keeping in mind that the system evolves towards the state of minimum energy, two potentials are considered to investigate the influence of γ us . The first potential, ψ PN,1 , is calibrated by setting γ us equal to the relaxed unstable stacking fault energy of ψ * gp , i.e., γ us = γ (r) us . Note that this choice introduces with the constraint ∆ n = 0 a slight inconsistency with the atomistic response. The second potential, ψ PN,2 , employs γ us = γ (u) us instead. Both potentials ψ PN,i along with their tractions T t,i = ∂ψ PN,i /∂∆ t , i ∈ {1, 2}, are shown in Figure 4 as a function of the tangential disregistry ∆ t .

Fourier potential
Assuming that for any ∆ t the glide plane potential always occupies the state of minimum energy in ψ * gp for that ∆ t , i.e., where T * n (∆ t ) = ∂ψ * gp /∂∆ n = 0, the corresponding potential ψ gp (∆ t ) cannot be represented with sufficient accuracy by the simple sinusoidal function of Eq. (25). Although the amplitudes of the energy and traction are captured accurately, their profiles are inaccurate (cf. Figure 4, in which the dashed curves represent the atomistic data). A higher complexity is therefore required, which can be achieved through a Fourier series expansion for the glide plane tractions. The corresponding potential is obtained through integration and reads The Fourier coefficients γ us,k in Eq. (26) are determined by a least-squares fit of the shear traction T t,F = ∂ψ F /∂∆ t to T * t (∆ † n , ∆ t ). In analogy to the classical PN potential, the normal disregistry is neglected. The first four obtained Fourier coefficients are listed in Table 4, whereas higher order terms (i.e., k > 4) are dropped because of their negligible influence. The fitted potential ψ F (∆ t ) and its traction T t,F (∆ t ) are shown in Figure 4. As can be observed, it fits the atomistic data nearly perfectly along the relaxed path (∆ † n , ∆ t ) in the diagrams.

Coupled analytical potential
Instead of restricting the potential ψ gp to represent only one section ψ * gp (∆ † n , ∆ t ), the full potential ψ * gp (∆ n , ∆ t ) can be approximated by introducing the additional dependency of ψ gp on the normal disregistry ∆ n . As a result, also the out-of-plane reaction force of the glide plane is included, as opposed to the previous models which neglected this contribution. Following Sun et al. [17], the coupled potential can be expressed as with the work of separation G c , the characteristic length l c (i.e., such a normal disregistry ∆ n at ∆ t = 0 for which ∂ 2 ψ * gp /∂∆ 2 n = 0), and p = ∆ † n (∆ t = b/2)/l c . Sun et al. [17] further suggested to set q equal to the ratio γ us to be a key quantity of the glide plane potential, q is to be calibrated such that γ (r) us is represented correctly. This relationship yields, cf. [17], The correspondingly calculated parameters for ψ C are listed in Table 5. The energy ψ C and shear traction T t,C are plotted in Figure 4 along ∆ † n , ∆ t , to allow for an adequate comparison with the other glide plane potentials.

Coupled numerical potential
The highest possible accuracy of the glide plane potential for the PN-FE model is obtained if the data points of ψ * gp are directly used as input for the glide plane potential, which is denoted as ψ N . Intermediate values are calculated numerically via a cubic interpolation scheme to ensure sufficient smoothness of the glide plane traction ∂ψ N /∂ ∆ and stiffness ∂ 2 ψ N /∂ ∆ 2 . Figure 4 shows the energy ψ N and corresponding traction T t,N along ∆ † n , ∆ t .

Comparative analysis
In this section a two-fold objective is pursued. i) The predictive capability of the PN-FE methodology is assessed by a comparison against the MS simulation. As a benchmark of the optimum accuracy the coupled numerical glide plane potential ψ N (∆ n , ∆ t ) (cf. Section 3) is used. ii) The influence of the various glide plane potentials (cf. Section 3) on the obtained results is studied.
First, the employed geometric properties of the two-phase microstructure are specified in Section 4.1. The results for a single dislocation are discussed in Section 4.2 in terms of the disregistry and traction profiles, as well as the displacement and strain fields. Similarly, the results (without traction profiles) for a 3-dislocation pile-up configuration follow in Section 4.3. Lastly, the full evolution process, from zero dislocations to dislocation transmission, is discussed in Section 4.4 by means of the dislocation positions and the total free energy.

Parameter set
The two-phase microstructure, as introduced in Section 2, is considered, with 512 × 516 atoms in Phase A and two times 258

Single dislocation 4.2.1. Disregistry profile
Consider first an externally applied shear load τ = 0.00473 µ A , at which in both models (MS and PN-FE) a single dislocation exists, obstructed by the phase boundary. The local dislocation configuration of the MS model and the PN-FE model with its different glide plane potentials is illustrated in Figure 5 in terms of the disregistry profiles ∆ t (Figure 5a) and ∆ n (Figure 5b). While for the PN-FE model, the disregistries are straightforwardly calculated from the nodal displacements through Eq. (22), the atomistic model requires first the interpolation of the atomic displacements above and below the glide plane before Eq. (22) can be applied.
As a characteristic of the PN model, the presence of the dislocation is indicated by the drop in disregistry from ∆ t = b (fully slipped) to ∆ t = 0 (non-slipped), which in the PN-FE model is established through the minimisation of the total free energy of Eq. (18) -without the requirement of additional criteria. Naturally, the dislocation is taken to be situated at ∆ t = b/2.
The comparison of the PN-FE model (with the coupled numerical glide plane potential ψ N ) with the MS model shows that, despite a deviation within the dislocation core, the rather simple PN-FE model is able to approximate the disregistry profiles (∆ t and ∆ n ) of the atomistic model (and dislocation position) reasonably well. The deviations within the dislocation core are assumed to origin from the local and linear-elastic continuum formulation used in the PN-FE model. Note that the relatively large deviations in ∆ n in the tailing part of the dislocation is an artefact. It origins from the MS model where ∆ is calculated in the reference configuration (at t = 0). In this context, atoms which are increasingly separated horizontally by ∆ t experience in relation with the external shear load an increasing difference in vertical displacement, that artificially leads to an increase in ∆ n .
The study of the different analytical glide plane potentials ψ gp unveils a significant influence on the local tangential dislocation disregistry profile ∆ t . Consider the PN-FE model with the coupled numerical glide plane potential ψ N as the reference solution. The classical PN potentials (ψ PN,1 and ψ PN,2 ) exhibit certain deviations in ∆ t , which can be related to the difference in ψ gp as follows (cf. also Figure 4a). With ψ PN,1 small deviations of ∆ t from the energy-free state (∆ t = i b with i = 0, 1, 2, . . . ) are energetically less penalised, as compared to ψ N . Consequently, the surrounding bulk outside of the dislocation core relaxes more to reduce elastic strain energy, leading to the slow decay/increase towards ∆ t = i b. On the contrary, ψ PN,2 penalises any tangential disregistry ∆ t to a significantly larger extent. The dislocation core is therefore more compressed to reduce the contribution of the glide plane potential to the total energy, resulting not only in a fast decay/increase towards ∆ t = i b, but also in a higher gradient within the dislocation core. Comparing the glide plane potentials ψ F and ψ C with ψ N shows very closely matching tangential disregistry profiles. This supports the in section 3 stated assumptions that: (i) in the PN-FE model, γ (r) us is the key quantity for the proper description of the dislocation behaviour; (ii) the disregistry profile ∆ of a dislocation follows the path of minimum energy in ψ * gp , i.e., along ∆ † n . Only the coupled numerical and coupled analytical glide plane potentials ψ N and ψ C include the dependency on ∆ n , and are thus capable of approximating the normal disregistry ∆ n of the atomistic model. The Fourier potential, however, while neglecting the normal disregistry, nevertheless does a good job at capturing the tangential disregistry by virtue of its calibration on the relaxed energy Ψ * gp .

Glide plane tractions
The predictions made with the PN-FE model in terms of the dislocation induced glide plane tractions T t , and the implication on the dislocation behaviour are discussed here. An illustrative comparison of the tractions is presented in Figure 6, corresponding to the disregistry profiles in Figure 5. The atomistic tractions are taken as the Virial shear stress σ xy (recall Eq. (11)) at atom positions below the glide plane. Note that the Virial stress, since it is an averaged quantity, becomes questionable within the dislocation core. Hence, no comparison between the atomistic model and the PN-FE model in the core region (±5b) follows. Outside of the dislocation core, nevertheless, it can be seen that the traction distribution in the PN-FE model (with ψ N ) agrees to the MS model. A minor irregularity of the traction with a discontinuity at the phase boundary is present in the PN-FE model, which relates to the discontinuity in material properties across the phase boundary. This similarity outside of the dislocation core exemplifies the good representative capability of the PN-FE model for dislocation induced long-range stresses in two-phase microstructures. For sufficiently separated dislocations in a pile-up (> 5b), similar dislocation induced repulsive shear stresses can thus be expected.
A comparison of the glide plane potentials on the shear tractions T t , shows no influence on the tractions outside of the dislocation core (±5b). Hence, no difference in the pile-up configuration is expected if dislocations are no closer than ≈ 5b. Within the dislocation core, however, and similar to the tangential disregistry profile, a strong influence on the tractions is unveiled. Thus, differently strong repulsive shear stresses are assumed for sufficiently close dislocations (i.e., ≤ 5b). In such a circumstance the classical PN potentials would exhibit a substantially larger (PN2 with γ (u) us ) or slightly lower (PN1 with γ (r) us ) repulsion as compared to the coupled numerical potential. The Fourier potential and the coupled analytical potential, on the contrary, would invoke a similar repulsion as the coupled numerical potential -and thus exhibit similar pile-up configurations.

Displacement and strain fields
The local displacement and strain fields are next evaluated for the MS model and the PN-FE model with selected glide plane models. Considered are the coupled numerical potential (ψ N ) and the Fourier potential (ψ F ), which shows relative good agreement in ∆ t and T t but is constrained in ∆ n . For the atomistic model, the displacements are plotted in terms of the change of atom positions. The strains correspond to the Green-Lagrange strain tensor E G obtained from the deformation gradient, which is calculated by a local least squares fit of the displacements of neighbouring atoms relative to the central one with respect to a homogeneous deformation gradient, as presented in [43,44]. The post-processing is performed in Ovito [45]. Because the strain fields are averaged quantities in MS, unreliable values are obtained within the dislocation core region. The strains are, therefore, to be evaluated with caution within the core regions. To this end, the used colour scales are cropped. For the PN-FE model, either directly the nodal quantities (for displacements) or the element averaged nodal quantities (for infinitesimal strains E) are shown. Due to (i) the questionable accuracy of strains within the core region in the atomistic model, and (ii) the difference between the geometrical non-linearity of the atomistic model vs. linearity of the PN-FE model, the obtained results can only be compared qualitatively. All figures span the window [100 a 0 ≤ x ≤ 320 a 0 , −60 a 0 ≤ y ≤ 60 a 0 ]. The resulting quantities are shown in Figure 7.
The comparison shows a good agreement between the PN-FE model (ψ N and ψ F ) and the atomistic model. Notwithstanding its simplicity, the PN-FE model is able to capture the quite complex displacement and strain fields rather well. In both models, MS and PN-FE, similar discontinuities at the phase boundary are noticeable, e.g., in E xx , E G,xx , E xy , and E G,xy . These discontinuities are due to the change in material properties across the phase boundary, and hence are more pronounced in the PN-FE model, in which this boundary is sharper. Note that the large strains along the glide plane in the MS model (E G,yy and E G,xy ) relate to the local least squares fit and occur in this context due to the large relative displacement of atoms across the glide plane.
Despite the relatively good agreement, some differences between the atomistic and the PN-FE model with ψ N are present, especially for E yy and E G,yy below the glide plane. The Fourier glide plane potential ψ F , on the contrary, achieves better qualitative agreement with the atomistic model. This, however, occurs in relation with the glide plane constraint ∆ n = 0, which limits the dislocation relaxation out of plane, and hence may not be misinterpreted as a higher accuracy of ψ F . Both glide plane potentials show minor deviations for E xx and E G,xx .

Dislocation pile-up 4.3.1. Disregistry profile
To study the model responses under a pile-up configuration, an externally applied shear load of τ = 0.0106 µ A is considered, at which in the MS model and the PN-FE model a 3-dislocation pile-up has formed. The disregistry profiles for the pile-up configuration is plotted in Figure 8. A mismatch of the dislocation positions between the MS model and the PN-FE model is apparent. In consideration of the equal long-range stresses (recall Figure 6), this difference can be related to the Peierls barrier. In the atomistic system, the lattice discreteness gives rise to the Peierls barrier which poses a resistance against dislocation motion and hence to the larger pile-up length than in the continuum PN-FE model, where no Peierls barrier is present. Again, an artificial contribution is included in the calculation of the normal disregistry ∆ n which naturally increases with a larger ∆ t (and a larger number of dislocations).
The comparison of the different potentials unveils an indifference of the adopted glide plane potential on the dislocation positions, which is in alignment with the equality of the long-range repulsive shear stresses for sufficiently separated dislocations (> 5b). An evaluation of the glide plane shear tractions, as done for the single dislocation case (recall Figure 6), is here omitted as it does not provide any additional insight.

Displacement and strain field
The corresponding displacement and strain fields in [100 a 0 ≤ x ≤ 320 a 0 , −60 a 0 ≤ y ≤ 60 a 0 ] are plotted in Figure 9 for the MS model and the PN-FE model with the coupled numerical potential (ψ N ) and the Fourier potential (ψ F ). In addition to the insight, obtained for the single dislocation case (recall Section 4.2.3), the additional discrepancy in dislocation position (between the atomistic and the PN-FE models) is here apparent through the slight difference in the displacement and strain fields.

Pile-up evolution
In this section, the pile-up formation is studied for the MS and the PN-FE model in detail, i.e., from initially dislocation-free up to the transmission of the leading dislocation. The PN-FE model with ψ N is first compared with the atomistic model; a discussion on the influence of the adopted glide plane potential follows thereafter.
By subjecting the considered microstructure to a constantly increasing external shear deformation (corresponding to the shear load τ ), stable dislocations are successively nucle- ated in both models. As a consequence of that applied shear load, the nucleated dislocations tend to move towards the phase boundary where, as a result of the present phase contrast, they pile up. To this end, the parametrisation time t is updated by increments of 1/450 to reach the target strain τ = 0.03 µ A , until dislocation transmission is recorded. The detailed evolution of the dislocation positions as a function of the externally applied shear load τ is shown for the MS and the PN-FE model (ψ N ) in Figure 10a, where x j represents the horizontal coordinate of dislocation j. The corresponding evolution of the total free energy Ψ is illustrated in Figure 10b, with Ψ = V(r ) − V(r 0 ) for the MS model. The specific model evolutions, in terms of dislocation positions and total free energy can be understood as follows.
Initially, both models are defect free and behave (nearly) linear elastic. Only after the external shear deformation has increased sufficiently, dislocation nucleation is triggered, at τ = 0.0012 µ A for the atomistic system; in the PN-FE model, the criterion for dislocation nucleation is not fulfilled yet. By nucleating the dislocation, a large amount of energy is introduced into the system, which results in the energy increase (cf. Figure 10b). At the instance of nucleation (i.e., Point 1 in Figure 10), the stable and unstable equilibrium positions of the dislocation lie close to each other, meaning also that multiple local energy minima exist. The individual local minima are separated by energy barriers that need to be overcome to propagate the dislocation further -the so-called Peierls barriers -originating from the lattice discreteness. This phenomenon of multiple local minima, and thus Peierls barrier, does not exist in the continuum formulation of the PN-FE model. Hence, no stable dislocations are nucleated in the PN-FE model yet. It is also to be noted, that in an MD simulation the dislocation may overcome the barrier due to thermal activation and annihi-late (recall the dipole configuration). Nevertheless, due to the presence of that barrier in the MS model the externally applied shear load τ needs to be increased sufficiently to propagate the dislocation further towards the phase boundary into its new equilibrium position (τ = 0.0033 µ A -Point 2). Subsequently, a sudden motion occurs, which is reflected in the energy curve by a large drop in the energy to a slightly higher level compared to the PN-FE model. Shortly thereafter (τ = 0.0035 µ A -Point 3), the shear deformation suffices to trigger dislocation nucleation in the PN-FE model close to the phase boundary. Here, the energy introduced to nucleate the dislocation equals the amount of energy by which the system relaxes due to the presence of the dipole. Not a jump, but rather a kink is therefore observed in Ψ(τ ). At this point, a slight difference in Ψ between the MS and the PN-FE model reveals a minor underestimation of the dislocation induced energy by the PN-FE model. With further increasing externally applied shear load, both models exhibit similar behaviour in terms of the energy and position of the first dislocation. The minor differences in dislocation position is explained by the presence of the Peierls barrier, which leads to a step-wise motion of the dislocations in the atomistic model with jumps of approximately one Burgers vector b -a phenomenon that is not present in the PN-FE model due to its continuous nature.
A second dislocation is nucleated in the atomistic model at τ = 0.0051 µ A (Point 4). An energy increase is again required to nucleate it, and the dislocation remains temporarily in its nucleation position until the applied shear load reaches the value of τ = 0.0065 µ A . Again, in MD simulations this nucleated dislocation may not be stable. By further increasing the applied shear load, the dislocation moves towards the phase boundary where it piles up, leading to a relaxation of the system and a small drop in energy (Point 5). At approximately τ = 0.0065 µ A (Point 6), the second dislocation is nucleated in the PN-FE model, which now requires some energy, resulting in a small jump. With the second dislocation the difference in energy between both models increases, as opposed to where only one dislocation was present, due to the marginal underestimation of the dislocation induced energy in the PN-FE model. The comparison of the dislocation positions (Figure 10a) furthermore reveals a small deviation between both models which, similar to the step-wise dislocation motion, can be attributed to the Peierls barrier.
With the nucleation of the third dislocation the mechanism changes. Not only does the dislocation nucleate in the atomistic model at a higher shear load (τ = 0.0103 µ A -Point 8) compared to the PN-FE model (τ = 0.0098 µ A -Point 7), it is also immediately mobile. Both the later nucleation and mobility result from the already existing dislocation structure.
In the atomistic model, the dislocation pile-up is at all stages less compressed and it hence induces a higher long-range backstress compared to the PN-FE model, triggering a later nucleation of subsequent dislocations. With the nucleation of the third dislocation a jump in the position of the second dislocation and the absence of the motion barrier for the third dislocation (in the MS model) is apparent, which can be as follows. At the beginning of the nucleation step the third dislocation is introduced at position x 3 where it initially exerts a repulsive shear stress on the 2-dislocation pile-up. Hence, in the course of minimising the total free energy (at this step), the pile-up is compressed which in turn lowers the repulsive shear stress of the pile-up on the nucleated dislocation. As a result, the third dislocation Eventually, the leading dislocation is transmitted in the PN-FE model at an externally applied shear load of τ = 0.011 µ A . In the atomistic model, a significantly higher stress level of τ = 0.0167 µ A is required to reach transmission, at which the size of the pile-up has already increased by one dislocation. For the sake of conciseness and to facilitate a sharp comparison between both models this is not included in Figure 10. The cause of the different obstruction against dislocation transmission lies within the complex interaction of the dislocation core with the phase boundary during its transmission across the phase boundary. Due to the strong approximations made on the core behaviour in the PN-FE model (i.e., small deformations, linear elasticity, local glide plane potential), it is not able to capture the process of dislocation transmission accurately.
A comparison of the energy evolution obtained with the different glide plane potentials, plotted in Figure 11, shows a small difference up to a shear load of τ ≈ 0.011 µ A . Similar to the dislocation glide plane shear tractions (recall Figure 6), PN1 (with γ (r) us ) relates to the lowest total free energy and PN2 (with γ (u) us ) to the highest, whereas the energy for the Fourier, the coupled analytical and coupled numerical potentials exhibit a (nearly) matching total free energy. Most glide plane potentials invoke dislocation transmission at τ = (0.0109 ± 0.00017) µ A . Only with ψ PN,2 (fitted to γ (u) us ) dislocation transmission in triggered at a strongly increased shear load with τ = 0.0214 µ A , where already six dislocations are nucleated. This illustrates the significant influence of γ us on dislocation obstruction at the phase boundary.

Summary and discussion
This paper presented a qualitative and quantitative comparison of the PN-FE model against equivalent atomistic simulations in terms of edge dislocation behaviour using a 2D two-phase microstructure. Different glide plane potentials were considered to study their influence on the dislocation structure (i.e., disregistry profile) and on dislocation transmission across the phase boundary. Although the continuous PN-FE methodology is based on numerous simplifications, it is able to capture the atomistic response to a good extent in terms of dislocation positions, disregistry profile and strain field. This exemplifies that the PN-FE model is able to capture the long-range influence of a second phase on the dislocation behaviour fairly well. Some quantitative differences between the atomistic and PN-FE model are present originating from a distinct dislocation core behaviour in terms of the disregistry profile, the pile-up length (≈ 15%) and the shear load which is required for dislocation transmission. These differences can be related to the simplicity of the PN-FE model as follows: (i) linear elasticity describes the behaviour of the bulk above and below the glide plane; (ii) the problem is solved in an infinitesimal strain framework; (iii) the material properties are homogeneous within the individual phases, with a jump across the interface; (iv) the behaviour of a discrete (atomic) system is adopted in a continuum formulation, which thus cannot capture the discrete nature of the Peierls barrier; and (v) the effect of non-locality is neglected. Although it is possible to advance the PN-FE model to higher complexity, its simplicity and computational efficiency would be compromised.
The comparison of the different glide plane potentials of the PN-FE model showed that the amplitude of the glide plane potential has a significant influence on dislocation transmission. Glide plane potentials fitted to the relaxed unstable stacking fault energy γ (r) us invoked dislocation transmission at a similar externally applied shear deformation compared to atomistic simulation. The glide plane model fitted to the unrelaxed unstable stacking fault energy γ (u) us , on the contrary, required a much higher shear deformation to trigger dislocation transmission. In terms of tangential disregistry profile ∆ t , the highest accuracy (compared to the reference numerical potential ψ N (∆ n , ∆ t )) was obtained by the coupled potential ψ C (∆ n , ∆ t ) and the Fourier series based potential ψ F (∆ t ), although the latter disregards the normal disregistry ∆ n . The normal disregistry ∆ n is only incorporated in the coupled model ψ C and matches adequately with the reference model ψ N .
Although numerous simplifications are adopted in the PN-FE model, it has proven to be an elegant and computationally efficient approach to study edge dislocation interactions with phase boundaries. It enables to study isolated mechanisms of the underlying physics, since only the key mechanisms employed in the model are activated.