# Three-dimensional phase-field model of dislocations for a heterogeneous face-centered cubic crystal

Antoine Ruffini, Yann Le Bouar, Alphonse Finel

► To cite this version:

Antoine Ruffini, Yann Le Bouar, Alphonse Finel. Three-dimensional phase-field model of dislocations for a heterogeneous face-centered cubic crystal. Journal of the Mechanics and Physics of Solids, 2017, 105, pp.95-115. 10.1016/j.jmps.2017.04.008. hal-01632435

HAL Id: hal-01632435

https://hal.science/hal-01632435

Submitted on 5 Jul 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Three-dimensional phase-field model of dislocations for a heterogeneous face-centered cubic crystal

Antoine Ruffini*, Yann Le Bouar, Alphonse Finel

Laboratoire d'Étude des Microstructures, CNRS-ONERA, B.P. 72, F92322 Châtillon Cedex, France

## Abstract

A central aim of current materials studies is to develop a predictive modeling that incorporates dislocation-based plastic activity and microstructural evolution. Phase-field method has emerged as a powerful tool for addressing this issue, providing us with a versatile variational framework able to describe the movement of dislocations in interaction with underlying microstructures. In this article, a three-dimensional phase-field model of dislocations (PFMD) is developed with a discretization scheme that explicitly captures the face-centered cubic (FCC) geometry. Within this framework, continuous fields are discretized in a way that allows to consider strongly heterogeneous materials and sharp interfaces (free surfaces, stiffer precipitates, pores...) without generating numerical artifacts. The PFMD exposed in this work reproduces dislocation activity in FCC geometry, their reactions, and a particular attention is devoted to the dislocation core behaviors in order to remove effects present in prior generic PFMDs that can appear to be spurious for micron-scale applications. This allows us to rigorously reproduce the dislocation's velocity with respect to experimental friction coefficients. The model is discussed and illustrated by applications standing at different space-scales that show how dislocations operate with microstructural heterogeneities such as free-surfaces (cylindrical nanopillar) and voids (pore under isostatic pressure).

**Keywords:** Dislocations, Phase-Field, Numerical Modeling, Nanopillar, Pore

## 1. Introduction

The macroscopic mechanical behavior of heterogeneous metallic alloys results from the evolution of its microstructure (solid precipitates, pore network...) generally coupled to microscopic phenomena (dislocations, microcracks...). Numerically however, at a space and time-scale of the industrial processes, it is impossible to reproduce this behavior without introducing constitutive laws that partially lack physically based justifications [1, 2]. As a consequence, these models can fail to accurately describe and predict the evolution of a given system under specific mechanical loads [3].

<footnote>*Corresponding author</footnote>

<footnote>Email address: antoine.ruffini@onera.fr (Antoine Ruffini)</footnote>
A better understanding of the mechanical behaviors of heterogeneous materials can be achieved by explicitly taking into account the crystalline defects that are at the origin of these behaviors, in particular dislocations. When the typical distance between dislocations is much smaller than the characteristic length-scale of internal heterogeneities, dislocations can be introduced at a continuum level using dislocation densities [4, 5]. However, the development of a crystal plasticity model based on a rigorous transition between discrete dislocations and continuous densities is a challenging task despite recent advances [6, 7].

When the average distance between dislocations is of the order or larger than other internal length-scales (for example, the average distance between internal interfaces), plasticity has to be accounted for by considering individual dislocations. In this context, a relevant model should include the interactions between dislocations and other defects (precipitates, interfaces, voids...) as well as the coupled dynamics of these defects.

Since more than a decade, Phase-Field Models (PFMs) have emerged as very powerful methods to reproduce, within a well controlled variational formalism, the out-of-equilibrium evolution of microstructures. Therefore, it is of great interest to develop a PFM formulation of dislocation glide (PFMD) that will permit, contrary to actual Discrete Dislocation Dynamics (DDD) techniques [8, 9], a straightforward coupling between dislocations and microstructural evolution [10–14]. However, PFMDs can still face technical difficulties to correctly couple the dislocations with strong elastic heterogeneities [15]. In addition, because of their continuous formulations, PFMDs offer multiscale possibilities for which relevant numerical implementations must be set.

For instance, when the issue is to study the core structure of dislocations, a subatomic refined mesh is required in order to properly reproduce the continuous atomic-scale variations of the core profile [16–19]. In one such situation, the underlying symmetry of the mesh can arbitrarily be chosen since, if the grid spacing is small enough, it does not interfere with the physical symmetry of the PFMD. However, due to the required subatomic grid spacing, the size of the simulated systems that can be afforded is necessarily limited. Alternatively, one could consider a PFMD numerical implementation that can potentially stand at the micron-scale. In this case, a non-refined mesh would have to be considered in order to reach relevant space scales. However, its intrinsic symmetry can then detrimentally interfere with that of the crystalline materials.

Because face-centered cubic (FCC) crystallographic symmetry is very common, it is crucial to develop a numerical scheme that incorporates the specificity of this symmetry (glide planes, dislocation line directions and characters...) whatever the scale used for the numerical implementation. Our solution to this constraint is to adopt a numerical grid which is always homothetic to crystalline symmetry such that FCC materials can be described irrespective of the mesh refinement. In this work, this issue is addressed through the elaboration of a three-dimensional PFMD, considering unambiguous free surfaces (with non-interpolated fields) on an explicit FCC grid. Due to its continuous formulation, the present model is expected to propose a numerical framework in which other kinds of continuous fields can also be introduced (diffusion, phase transformations...).

This paper is organized as follows: We first describe the elastic model and its numerical implementation that explicitly captures the FCC underlying symmetry. Then, we implement a PFMD that correctly reproduces the coplanar
dislocation reactions in this geometry. The dynamics of the model are modified which allows us to obtain a better description of the dislocations' physical properties. Finally, we validate the model and illustrate its potentiality when dislocations operate in the presence of heterogeneities made of voids and free surfaces.

## 2. Elasticity

In all PFMDs, the starting point is a continuum material whose constitutive behavior is based on linear elasticity [10–14]. In this section, we begin by recalling this background and describe how it is numerically implemented to account for the underlying FCC symmetry where microstructural heterogeneities may develop.

### 2.1. Elastic model

Within linear elasticity [20], the total free elastic energy of a material is defined by the following functional:

$$F_{el} = \int_V f_{el} \{\varepsilon_{ij}(\mathbf{r})\} dV, \quad (1)$$

with V its total volume and $f_{el}$ the free elastic energy density given by:

$$f_{el}(\mathbf{r}) = \frac{1}{2} C_{ijkl}(\mathbf{r}) \varepsilon_{ij}(\mathbf{r}) \varepsilon_{kl}(\mathbf{r}) - \sigma_{ij}^{A} \varepsilon_{ij}(\mathbf{r}), \quad (2)$$

where $C_{ijkl}$ are the elastic tensor components and $\varepsilon_{ij}$ are the elastic strains defined at position $\mathbf{r}$. The components $\sigma_{ij}^A$ account for the applied stress. In the present work, we consider small strains:

$$\varepsilon_{ij}(\mathbf{r}) = \frac{\beta_{ij}(\mathbf{r}) + \beta_{ji}(\mathbf{r})}{2}, \quad (3)$$

where $\beta_{ij}$ denotes the displacement-gradients given by:

$$\beta_{ij}(\mathbf{r}) = \frac{\partial u_i(\mathbf{r})}{\partial X_j}, \quad (4)$$

with $u_i$ and $X_j$ the displacements and reference positions respectively.

Then, displacements $u_i$ are decomposed into a homogeneous part $\bar{u}_i$ and a heterogeneous part $\Delta u_i$ such that $u_i = \bar{u}_i + \Delta u_i$. The homogeneous part is defined by $\bar{u}_i(X_j) = \varepsilon_{ij}X_j$ and reflects the loading conditions through the average strains $\varepsilon_{ij} = \frac{1}{V} \int_V \varepsilon_{ij}(\mathbf{r}) dV$ (hard driving).

In classical PFMD, dislocations evolve under mechanical equilibrium conditions while their motion is assumed to result from dissipative frictional phenomena (phonon drag, interactions with obstacles...). To do so, the elastic fields (displacements, strains, stresses) are instantaneously relaxed at each time step of the dislocation motion. For a strongly heterogeneous material, the numerical strategy that must be considered for the mechanical equilibrium to be efficiently reached is a study on its own which is not addressed in the present paper. Here, the purpose lies rather in the way we can discretize the continuous fields on a numerical grid whose intrinsic symmetry does not detrimentally
interfere with that of an FCC crystalline material (provided the grid spacing is equal to or higher than the interatomic distance).

Therefore, we decide to consider a simple dissipative scheme on the elastic fields that makes dislocations evolve within a quasi-static environment. For this, the elastic energy is minimized using the following dissipative equation:

$$ \frac{\partial \Delta u_i(\mathbf{r}, t)}{\partial t} = -M \frac{\partial \mathcal{F}_{el}}{\partial \Delta u_i(\mathbf{r}, t)}, \qquad (5) $$

where M is a dissipative coefficient and $\mathcal{F}_{el}$ the functional given in Eq. (1). When applied stresses $\sigma_{ij}^A$ are considered (soft driving), we make average strains $\epsilon_{ij}$ to converge toward the corresponding equilibrium values through the following dissipative scheme:

$$ \frac{\partial \epsilon_{ij}(t)}{\partial t} = -M_{\epsilon} \frac{\partial \mathcal{F}_{el}}{\partial \epsilon_{ij}(t)}, \qquad (6) $$

where $M_{\epsilon}$ denotes another dissipative coefficient. In order for the quasi-static hypothesis to be granted, M and $M_{\epsilon}$ are chosen high enough so that the characteristic time-scale of the dissipation-driven elastic minimization is much smaller than that of the dislocation motion (see subsection 3.2).

In elastically inhomogeneous systems, the way the mechanical equilibrium is set could be improved by the use of more efficient numerical solvers (e.g. Fourier space formulation [21–23] or non-linear conjugate gradient methods [24]). This will be addressed in a separate work.

## 2.2. Numerical implementation

The numerical implementation of the actual model requires two important features: (i) the mesh must satisfy the face-centered cubic (FCC) symmetry, and (ii) the components of a unique physical quantity must be defined at the same point. Therefore, whatever the space-scale used in the PFMD, the first point guarantees that the physical phenomena will necessarily respect the underlying crystallographic symmetry. The second point guarantees that strong microstructural elastic heterogeneities can be easily introduced without requiring a re-meshing procedure or the extrapolation of some physical quantities. This point will be discussed in more detail in subsection 2.3.

A primary FCC grid is built as displayed in Fig. 1.a, by considering four embedded simple cubic (SC) grids that are subsequently referenced with the colorized superscripts I, II, III and IV respectively (see colors on-line). As the calculation of the energy density requires the displacement-gradients (see Eqs. (1–4)), we localize the displacements $u_i$ on the nodes of the primary FCC grid and the displacement-gradients $\beta_{ij}$ on the octahedral sites constituting a secondary staggered FCC grid whose nodes actually correspond to the centers of the four SC grids, as shown in Fig. 1.b. With this scheme, the displacement-gradients $\beta_{ij}$, the strains $\epsilon_{ij}$ and the energy densities $f_{el}$ are all defined on the octahedral sites constituting the secondary FCC grid.

To highlight in more detail the way we calculate the elastic energy density, we focus on the elementary cell of the FCC lattice displayed in Fig. 2. We only reproduce one octahedral site in the center of a cubic cell of the first SC grid. The first diagonal of the octahedron (in red) represents a length element $dX_1$ of length d aligned with the first direction of space. The second diagonal (in green) stands for $dX_2$ and the third (in blue) for $dX_3$. The nodes are
<figure><img src="image_2.png" /><figcaption>Figure 1. a) Primary FCC grid built with four embedded SC grids referenced with the colorized superscripts I, II, III and IV respectively (see colors on-line). b) Secondary staggered FCC grid formed by the octahedral sites of the primary FCC grid. The nodes of the secondary grid actually correspond to the centers of the four SC grids that constitute the primary grid, with consistent nomenclature.</figcaption></figure>

characterized by their displacements $u_i$, indexed in space using the alphanumerical notation $a(+1), b(+1)$ and $c(+1)$, and colorized/superscripted according to their corresponding underlying SC grid. It may be noted that $u_i = x_i - X_i$ with $x_i$ and $X_i$ as the current and reference positions of a node respectively. One can also extract the displacements of the octahedral site $u_i^{ave}$ by taking the average displacements of the six face-centered nodes.

<figure><img src="image_3.png" /><figcaption>Figure 2. FCC unit cell where nodes supporting the displacements $u_i$ are identified by discs colorized/superscripted according to their underlying SC grid (see colors on-line). The octahedral site is located at the intersection of its three diagonals that also constitute three length elements $dX_j$ of length $d$. Squares are the actual sites where the displacement-gradients are calculated: $\beta_{i1}$ in red, $\beta_{i2}$ in green, and $\beta_{i3}$ in blue.</figcaption></figure>

To compute the displacement-gradients $\beta_{ij}$ defined in Eq. (4), a centered finite difference calculation can be implemented using the following discretization scheme:

$$ \beta_{i1} = \frac{\partial u_i}{\partial X_1} \equiv \frac{u_i^{II}(a+1, b, c) - u_i^{II}(a, b, c)}{d}, \qquad (7) $$

$$ \beta_{i2} = \frac{\partial u_i}{\partial X_2} \equiv \frac{u_i^{III}(a, b+1, c) - u_i^{III}(a, b, c)}{d}, \qquad (8) $$
$$ \beta_{i3} = \frac{\partial u_i}{\partial X_3} \equiv \frac{u_i^{IV}(a, b, c + 1) - u_i^{IV}(a, b, c)}{d} \quad (9) $$

The issue with this scheme is that it introduces oscillating deformed modes with zero elastic energy. For instance, considering $u_i \propto \cos(\pi(X_j \pm X_k))$ leads with Eqs. (7–9) to $f_{el} = 0$ unlike what it is theoretically (and physically) expected. Actually, this pitfall is not due to the adoption of a centered difference scheme. It is due to the fact that the previous scheme only involves the FCC second-neighbor nodes.

It is well-known that the first-neighbor connectivity of the FCC lattice, which includes first-neighbor equilateral triangles and regular tetrahedrons, may lead to important frustration effects. For example, in the domain of statistical mechanics, the correct modeling of the Ising model for an ordering alloy that lies on the FCC lattice requires the use of a model that incorporates all the many-body correlations included within the tetrahedrons and octahedrons, even in the simplest case when only first-neighbor interatomic interactions are considered [25, 26].

In the present context, where elastic relaxations are involved, we restore the first-neighbor connectivity by splitting the $\beta_{ij}$ components along the diagonals of the octahedron by using the average displacements of the octahedral site $u_i^{ave}$. In Fig. 2, the $\beta_{i1}$ components are identified by the red squares labeled 1 and 2, $\beta_{i2}$ by the green labeled 3 and 4, and $\beta_{i3}$ by the blue labeled 5 and 6. Thus, octahedrons are divided into eight tri-rectangular tetrahedrons leading to eight ways of computing the elastic energy density. For example, when we consider the tetrahedron delimited by edges 1, 3 and 6, we obtain:

$$ \beta_{i1} \equiv \frac{u_i^{ave} - u_i^{II}(a, b, c)}{d/2}, \quad (10) $$

$$ \beta_{i2} \equiv \frac{u_i^{ave} - u_i^{III}(a, b, c)}{d/2}, \quad (11) $$

$$ \beta_{i3} \equiv \frac{u_i^{IV}(a, b, c + 1) - u_i^{ave}}{d/2}, \quad (12) $$

allowing to compute the displacement-gradient tensor and therefore the elastic energy of the corresponding tetrahedron. This scheme allows the first/second-neighbor connectivity to be granted since first/second-neighbor displacements contribute to $u_i^{ave}$ and therefore to the numerical definition of the displacement gradient $\beta_{ij}$. Repeating this procedure for the seven other tri-rectangular tetrahedrons and summing-up the pondered energy contributions of the eight tetrahedrons finally leads to the elastic energy density that is considered on the octahedral site. The present scheme does not lead to the existence of zero-energy deformed modes.

Note that even though the discretized procedure has been highlighted considering the octahedral site constituting the center of the first SC grid, we could also have used octahedral sites supported by SC grids II, III and IV for which, in Eq. (10–12), other nodes supporting the displacement field would have to be involved.

In addition to the initial requirements, this scheme also allows us to use the same number of sites for all the physical quantities considered thus far. Computationally, this enables us to define these quantities with tables of identical size without having to consider empty cells. This provides an efficient way for the computer to access the values, which decreases the numerical cost as well as minimizing the required memory. Having to consider FCC grids is also an
advantage when the time comes to define the discretized forms of the spatial differential operators (see Appendix A).

Finally, the discretized formulation of the total elastic energy $\mathcal{F}_{el}$ given Eq. (1) writes:

$$ \mathcal{F}_{el} = \sum_{octa. sites} f_{el} \times V_{octa}^{eff}, \qquad (13) $$

where $V_{octa}^{eff} = d^3/4$ stands for the effective volume of an octahedral site and $f_{el}$ is the elastic energy density expressed with the discretized form of the $\beta_{ij}$ components that explicitly involve $u_i$. Thus, by taking the functional derivative of $\tilde{\mathcal{F}}_{el} = \mathcal{F}_{el}/(Kd^3)$, with K the bulk modulus, the discretized formulation of Eq. (5) is:

$$ \frac{\partial \Delta \tilde{u}_i(\tilde{X}_j, \tilde{t})}{\partial \tilde{t}} = -\frac{\tilde{M}}{4} \sum_{octa. sites} \frac{\partial \tilde{f}_{el}}{\partial \Delta \tilde{u}_i(\tilde{X}_j, \tilde{t})}, \qquad (14) $$

where tildes denote the reduced quantities $\tilde{u}_i = u_i/a_0$, $\tilde{X}_j = X_j/d$, $\tilde{f}_{el} = s^2 f_{el}/K$, $\tilde{M} = M/(La_0^2)$ and $\tilde{t} = t/t_0$ with $t_0 = s^2/(KL)$. To define these reduced quantities, the scale factor $s = d/a_0$ is introduced and accounts for the ratio between the numerical grid spacing $d$ and the physical FCC lattice parameter $a_0$. It may also be noted that the reduced mobility $\tilde{M}$ and time unit $t_0$ involve a mobility coefficient $L$ which will be related below to the dislocation motion (see subsection 3.1). In order for the quasi-static hypothesis to be granted (the characteristic time-scale of the elastic minimization is much smaller than that of the dislocation motion), we will take $\tilde{M} = M/(La_0^2) = 10 \gg 1$. In this equation, it may be noticed that the sum only refers to the six closest octahedral sites surrounding a node since its displacements only contribute to the energy density calculated on these sites. Using these notations, the discretized formulation of Eq. (6) is:

$$ \frac{\partial \tilde{\epsilon}_{ij}(\tilde{t})}{\partial \tilde{t}} = -\frac{\tilde{M}_{\epsilon}}{4} \sum_{octa. sites} \frac{\partial \tilde{f}_{el}}{\partial \tilde{\epsilon}_{ij}(\tilde{t})}, \qquad (15) $$

where $\tilde{M}_{\epsilon} = M_{\epsilon}s^2/L$ is conveniently taken to be equal to $\tilde{M}$. From the previous equations, it can be noticed that reduced strains $\tilde{\epsilon}_{ij} = s\epsilon_{ij}$ (leading to average strains $\tilde{\epsilon}_{ij} = s\epsilon_{ij}$) have been introduced and consistently enable us to define the reduced applied stresses $\tilde{\sigma}_{ij}^A = s\sigma_{ij}^A/K$.

## 2.3. Pores, cracks and other heterogeneities

Within the present approach, a continuum material is characterized by displacements $u_i$ discretized on nodes that constitute the primary FCC lattice, and by displacement-gradients $\beta_{ij}$ (or strains $\epsilon_{ij}$, or elastic energy densities $f_{el}$) discretized on octahedral sites that constitute the secondary staggered FCC lattice.

In a material whose microstructure is coherent with the FCC lattice, what distinguishes a region from another is the way its elastic energy $f_{el}$ varies with strains $\epsilon_{ij}$. According to Eq. (2), this is controlled by the elastic coefficients $C_{ijkl}$. With our numerical scheme, these values can be chosen to be different from one position $\mathbf{r}$ to another (i.e.: from one octahedral site to another) since it does not affect the way the continuum fields are calculated and does not introduce any singular point in which a physical value would have to be interpolated because of its ambiguous localization.

A solid precipitate can therefore be defined in a specific region of the whole system by changing the elastic coefficients of the corresponding octahedral sites. This results in a domain which is delimited by sharp interfaces that
necessarily follow the FCC symmetry. With this scheme, we can also introduce a non-material region (a void) by canceling its corresponding elastic coefficients. This constitutes the strongest heterogeneity that may be considered. Physically, it allows us to define random shape cavities and even cracks when the non-material regions take the form of thin platelets. Hence, free surfaces of complex shape can naturally be involved without having the risk of generating singular fields.

As an example, we have simulated a $100d^3$ isotropic FCC material. In this case, we only have two independent elastic coefficients which are the bulk modulus $K$ and the shear modulus chosen to be $\mu = 0.40K$ (i.e.: with Poisson ratio $\nu \approx 0.33$). For instance, these coefficients could be aluminium ones by taking $K \approx 80$ GPa. Then we have introduced a spherical pore by canceling these coefficients within a central spherical region of radius $R_0 = 15d$, as it is shown in Fig. 3a. Periodic boundary conditions are considered. The system has been compressed isostatically with pressure $p = 0.00125K$ ($\approx 100$ MPa) and relaxed by dissipating its elastic energy. Some components of its stress-field $\sigma_{ij} = C_{ijkl}\epsilon_{kl}$ have been extracted and displayed in Figs. 3b–3d. This allows us to confirm that the interface of the pore follows the FCC symmetry and that the elastic field (stresses) is not perturbed by the presence of non-physical oscillating modes.

Concerning the characteristic length-scale, one can note that it is given by the scale factor $s = d/a_0$ which accounts for the ratio between the grid spacing $d$ and the FCC lattice parameter $a_0$. In a physical system, when the atomic aspects are deemed as essential (nanometer discontinuities in interfaces, crystalline defects in the bulk...), it is natural to select $s = 1$. In this case, we can adopt a naive representation where the nodes of the FCC grid are atoms interacting through the harmonic elastic energy potentials of the octahedral sites, and where the cancellation of elastic coefficients in one of these sites reproduces the breaking of bonds for the six surrounding atoms. Of course, since the model is a continuum model, higher length-scales can also be considered, where $s > 1$, especially when atomic short-range effects are negligible and all that matters is capturing the FCC symmetry (see subsection 6.2). On the contrary, smaller length-scales are conceivable when the subatomic energy variations must be described [16–19].

### 3. Phase-field model of FCC dislocations

The PFMD exposed in this work is based on the first approaches developed in Refs. [10, 11]. To our knowledge, these models have not significantly been improved since the work of C. Shen and Y. Wang [27] in which the coplanar FCC dislocation reactions are described more rigorously. In this section, this former model is reconsidered with some improvements that allow us to correctly reproduce the dislocation mobility.

#### 3.1. Set-up of the model

In the micro-elasticity theory, a dislocation is a linear defect defined by the front delimiting two coplanar crystallographic regions, one sheared and one not (the sheared region reflecting more fundamentally the displacement of one atomic plane to another). The sheared region can be viewed as a platelet-like inclusion labeled $p$ characterized by the
<figure><img src="image_4.png" /><figcaption>Figure 3. Spherical pore in a 100d³ compressed isotropic FCC material. In Fig. (a) the (111) clipping plane reveals the removed non-material part defining the pore. In Figs. (b), (c) and (d), the (001) clipping plane allows us to display the stress field components σ₁₂, σ₁₁ and σ₂₂, respectively (given in the unit of the bulk modulus K, see colors on-line).</figcaption></figure>

stress-free displacement-gradients $\beta_{ij/p}^0$ [28, 29]. For dislocations, these components write:

$$ \beta_{ij/p}^0 = \frac{b_{i/p} n_{j/p}}{h_p}, \qquad (16) $$

where $b_{i/p}$ is a component of the Burgers vector $b_p$, $n_{j/p}$ is a component of the normal unit vector of the dislocation slip-plane $n_p$ and $h_p$ is the interlayer spacing of the slip-plane.

In an FCC crystal, one can consider three 1/2⟨110⟩ slip-directions in four {111} slip-planes of identical spacing $h = d/\sqrt{3}$. This leads to twelve platelet-like types of inclusion that characterize the twelve $p \equiv 1/2\langle 110 \rangle\{111\}$ dislocation slip-systems. In Fig. 4, these systems are displayed using the nomenclature of Tab. 1.

A platelet-like inclusion characterizing one slip-system $p$ is described as a phase of order parameter $\eta_p$. This order parameter is chosen at position $r$ to be an integer that reflects the shear state of the corresponding element of volume. According to this point, and because $\beta_{ij/p}^0$ given in Eq. (16) stands for an elementary shear of the slip-system $p$, the
<figure><img src="image_5.png" /><figcaption>Figure 4. The twelve $p = 1/2\langle 110 \rangle\{111\}$ FCC slip-systems with nomenclatures of Tab. 1.</figcaption></figure>

<table><caption>Table 1. Normal unit vectors **n**<sub>p</sub> and Burgers vectors **b**<sub>p</sub> related to the twelve $p = 1/2\langle 110 \rangle\{111\}$ FCC slip-systems.</caption><thead><tr><th><b>n</b><sub>p</sub></th><th colspan="3"><b>b</b><sub>p</sub></th></tr></thead><tbody><tr><td><b>n</b><sub>1,2,3</sub> = <sup>1</sup>&frasl;<sub>&radic;3</sub> (1<span style="text-decoration: overline">11</span>)</td><td><b>b</b><sub>1</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [110]</td><td><b>b</b><sub>2</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [1<span style="text-decoration: overline">10</span>1]</td><td><b>b</b><sub>3</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [0<span style="text-decoration: overline">1</span><span style="text-decoration: overline">1</span>]</td></tr><tr><td><b>n</b><sub>4,5,6</sub> = <sup>1</sup>&frasl;<sub>&radic;3</sub> (1<span style="text-decoration: overline">11</span>)</td><td><b>b</b><sub>4</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [101]</td><td><b>b</b><sub>5</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [1<span style="text-decoration: overline">1</span>10]</td><td><b>b</b><sub>6</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [01<span style="text-decoration: overline">1</span>]</td></tr><tr><td><b>n</b><sub>7,8,9</sub> = <sup>1</sup>&frasl;<sub>&radic;3</sub> (111)</td><td><b>b</b><sub>7</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [1<span style="text-decoration: overline">1</span>0]</td><td><b>b</b><sub>8</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [1<span style="text-decoration: overline">1</span>01]</td><td><b>b</b><sub>9</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [01<span style="text-decoration: overline">1</span>]</td></tr><tr><td><b>n</b><sub>10,11,12</sub> = <sup>1</sup>&frasl;<sub>&radic;3</sub> (1<span style="text-decoration: overline">1</span>1)</td><td><b>b</b><sub>10</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [101]</td><td><b>b</b><sub>11</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [1<span style="text-decoration: overline">1</span>0]</td><td><b>b</b><sub>12</sub> = <sup>a<sub>0</sub></sup>&frasl;<sub>2</sub> [0<span style="text-decoration: overline">1</span><span style="text-decoration: overline">1</span>]</td></tr></tbody></table>

total stress-free displacement-gradients at position **r** are:

$$ \beta_{ij}^{0}(\mathbf{r}) = \sum_{p=1}^{12} \beta_{ij/p}^{0} \eta_{p}(\mathbf{r}), \qquad (17) $$

where the dislocation field $\eta_p(\mathbf{r})$ represents the amount of cumulative shear experienced by point **r** due to the slip-system $p$. The dislocation loops therefore introduce a sheared region characterized at position **r** by the total eigenstrain components:

$$ \varepsilon_{ij}^{0}(\mathbf{r}) = \frac{\beta_{ij}^{0}(\mathbf{r}) + \beta_{ji}^{0}(\mathbf{r})}{2}. \qquad (18) $$

Inside one such region, the variation of the elastic energy density $f_{el}$ must be computed taking the sheared state as the new reference state. To do so, the expression of $f_{el}$ given in Eq. (2) is modified by a change in the strain origin with eigenstrains given in Eq. (18). The strain components $\varepsilon_{ij}$ are thus replaced by $(\varepsilon_{ij} - \varepsilon_{ij}^0)$ and the next mentions of the elastic energy density will now refer to this new $\eta_p$-dependent definition: $f_{el}\{\varepsilon_{ij}(\mathbf{r}), \eta_p(\mathbf{r})\}$.

Then, we consider kinetic equations that control the temporal evolution of dislocations. For this, we adopt the time-dependent Ginzburg-Landau approach that requires us to define a phenomenological non-elastic energy written as [30]:

$$ \mathcal{F}_{GL} = \int_V f_{\text{cryst}}\{\eta_p(\mathbf{r})\} + f_{\text{grad}}\{\eta_p(\mathbf{r})\} dV. \qquad (19) $$
In Eq. (19), $f_{cryst}$ and $f_{grad}$ are the crystalline and gradient energy densities that stand for the Landau and Ginzburg terms respectively. We will also use the Ginzburg-Landau energy density $f_{GL} = f_{cryst} + f_{grad}$.

In the simplest formulation of the model, we only consider perfect dislocations with a diffuse realistic core structure. In this case, the crystalline energy density $f_{cryst}$ is chosen to stabilize the order parameter $\eta_p$ at values that are integers, since $\eta_p$ reflects the amount of elementary shear as defined in Eq. (17). We also chose to implement a model in which the coplanar dislocation interactions are correctly reproduced. To do so, the crystalline energy must reflect the fact that a dislocation of one slip-system can have several equivalent representations. For example, as it can be noticed in Fig. 4, a dislocation with Burgers vector $\mathbf{b} = -\mathbf{b}_3$ characterized by $(\eta_1 = \eta_2 = 0, \eta_3 = -1)$ is equivalent to $\mathbf{b} = \mathbf{b}_1 + \mathbf{b}_2$ with $(\eta_1 = \eta_2 = 1, \eta_3 = 0)$. This is achieved in the model developed in Ref. [27] where the crystalline energy writes:

$$ 
\begin{aligned}
f_{cryst} = & 2A \{4 - \cos[\pi(\eta_1 - \eta_2)] \cos[\pi(\eta_2 - \eta_3)] \\
& \qquad \cos[\pi(\eta_3 - \eta_1)] \\
& - \cos[\pi(\eta_4 - \eta_5)] \cos[\pi(\eta_5 - \eta_6)] \\
& \qquad \cos[\pi(\eta_6 - \eta_4)] \\
& - \cos[\pi(\eta_7 - \eta_8)] \cos[\pi(\eta_8 - \eta_9)] \\
& \qquad \cos[\pi(\eta_9 - \eta_7)] \\
& - \cos[\pi(\eta_{10} - \eta_{11})] \cos[\pi(\eta_{11} - \eta_{12})] \\
& \qquad \cos[\pi(\eta_{12} - \eta_{10})]\}.
\end{aligned} 
\quad (20) $$

In Eq. (20), A corresponds to the half amplitude of the crystalline energy density. Generally speaking, $f_{cryst}$ is related to the generalized stacking fault energy (or $\gamma$-surface) that can be extracted from first principle calculations [16–18] and A can be related to the unstable stacking fault energy [31].

In this formulation, the gradient term $f_{grad}$ is chosen to penalize the energy when the order parameters $\eta_p$ vary in space. This allows us to obtain a diffuse front of $\eta_p$ at its transition region, where the order parameter passes from one integer value to another (i.e.: where a crystallographic plane passes from one shear state to another). Here, we use the formula of the gradient energy density introduced in Ref. [27]:

$$ f_{grad} = \frac{B}{2} \sum_{p,q=1}^{12} \frac{\mathbf{b}_p \cdot \mathbf{b}_q}{b^2} (\mathbf{n}_p \times \nabla \eta_p) \cdot (\mathbf{n}_q \times \nabla \eta_q). \quad (21) $$

In Eq. (21), B is an amplitude scaling as a stress multiplied by a unit area which is chosen to be the same for each slip-system p. The symbol $\nabla$ stands for the gradient vector, and $\times$ is the cross product.

For a non-conserved order parameter, as it is the case for dislocations, the spatial-temporal evolution of $\eta_p$ is controlled by the Allen-Cahn equation given by [32]:

$$ \frac{\partial \eta_p(\mathbf{r}, t)}{\partial t} = -L \left( \frac{\delta \mathcal{F}_{el}}{\delta \eta_p(\mathbf{r}, t)} + \frac{\delta \mathcal{F}_{GL}}{\delta \eta_p(\mathbf{r}, t)} \right), \quad (22) $$
where L is a mobility coefficient. In subsection 3.3, it is shown how this parameter can be related to the mobility of dislocations.

## 3.2. Numerical details

For the numerical implementation, from Eq. (17), we introduce a reduced total stress-free displacement gradient $\tilde{\beta}_{ij}^0 = s\beta_{ij}^0$ with $s = d/a_0$ such that $\eta_p$ still reflects the amount of cumulative shear experienced by point r whatever the working space-scale. This is the same as defining $\tilde{\beta}_{ij}^0$ with reduced Burgers vectors $\tilde{b}_p = b_p/a_0$ and interlayer spacing $\tilde{h} = h/d$. At small strain, from Eq. (18), this procedure leads to the definition of reduced eigenstrains $\tilde{\epsilon}_{ij}^0 = s\epsilon_{ij}^0$. By using these reduced quantities, the reduced elastic energy density $\tilde{f}_{el}$ already introduced in Eqs. (14) and (15) writes:

$$ \tilde{f}_{el} = \frac{1}{2} \tilde{C}_{ijkl} (\tilde{\epsilon}_{ij} - \tilde{\epsilon}_{ij}^0) (\tilde{\epsilon}_{kl} - \tilde{\epsilon}_{kl}^0) - \tilde{\sigma}_{ij}^A \tilde{\epsilon}_{ij}, \quad (23) $$

with $\tilde{C}_{ijkl} = C_{ijkl}/K$ and, as above, $\tilde{\epsilon}_{ij} = s\epsilon_{ij}$ and $\tilde{\sigma}_{ij}^A = s\sigma_{ij}^A/K$.

Since the order parameter reflects the shear state of an element of volume, we decide to localize $\eta_p$ on the octahedral site of the FCC grid. As a consequence, the Ginzburg-Landau energy density $f_{GL} = f_{cryst} + f_{grad}$ as well as the amplitudes A and B are also defined on the octahedral sites. In this case, the discretized formulation of the non-elastic Ginzburg-Landau energy $\mathcal{F}_{GL}$ given Eq. (19) writes:

$$ \mathcal{F}_{GL} = \sum_{\text{octa. sites}} f_{GL} \times V_{\text{octa}}^{\text{eff}}, \quad (24) $$

where $V_{\text{octa}}^{\text{eff}} = d^3/4$. Thus, by taking the functional derivative of $\tilde{\mathcal{F}}_{GL} = \mathcal{F}_{GL}/(Kd^3)$, Eq. (22) becomes:

$$ \frac{\partial \eta_p(\tilde{X}_i, \tilde{t})}{\partial \tilde{t}} = -\frac{1}{4} \left( \frac{\partial \tilde{f}_{el}}{\partial \eta_p(\tilde{X}_i, \tilde{t})} + \frac{\partial \tilde{f}_{GL}}{\partial \eta_p(\tilde{X}_i, \tilde{t})} - \frac{\partial}{\partial \tilde{X}_j} \frac{\partial \eta_p(\tilde{X}_i, \tilde{t})}{\partial \tilde{X}_j} \right), \quad (25) $$

or by explicitly expanding the gradient energy contributions given Eq. (21):

$$ \frac{\partial \eta_p(\tilde{X}_i, \tilde{t})}{\partial \tilde{t}} = -\frac{1}{4} \left( \frac{\partial \tilde{f}_{el}}{\partial \eta_p(\tilde{X}_i, \tilde{t})} + \frac{\partial \tilde{f}_{cryst}}{\partial \eta_p(\tilde{X}_i, \tilde{t})} \right. \\ \left. - \tilde{B} \sum_{q=1}^{12} \frac{\tilde{b}_p \cdot \tilde{b}_q}{\tilde{b}^2} \left\{ n_{p/j} n_{q/j} \frac{\partial^2 \eta_q(\tilde{X}_i, \tilde{t})}{\partial \tilde{X}_j \partial \tilde{X}_q} \right. \right. \\ \left. \left. - n_{p/j} n_{q/k} \frac{\partial^2 \eta_q(\tilde{X}_i, \tilde{t})}{\partial \tilde{X}_j \partial \tilde{X}_k} \right\} \right\}, \quad (26) $$

where tildes denote the aforementioned reduced quantities with $\tilde{f}_{cryst} = s^2 f_{cryst}/K$, $\tilde{f}_{grad} = s^2 f_{grad}/K$ and considering $\tilde{A} = s^2 A/K$ and $\tilde{B} = B/(K a_0^2)$ with, as defined above, $t_0 = s^2/(KL)$. In this equation, it may be noticed that the sum is removed since the order parameter (and its spatial derivatives) only contribute to the energy density calculated on the same octahedral site. A quantitative estimation of $t_0$ for a typical situation (nickel-based superalloys) is given in the next subsection.
In Eqs. (21) and (26), the first and second spatial derivatives of the order parameters $\eta_p$ (gradient vector and Hessian matrix respectively) must also be introduced. These differential operators discretized on the FCC grid are calculated in Appendix A. As a matter of fact, these operators correspond to the FCC stencil used in the acoustics community to perform the finite difference approximations of the Laplacian (Hessian trace). Interestingly, the FCC stencil is numerically shown to provide good precision and an optimized computational efficiency compared to other cubic stencils [47]. In the present work, this FCC stencil emerges naturally because of the use of an explicit FCC grid.

### 3.3. Connection with dislocation’s physical properties

Let us consider a single FCC straight edge dislocation of Burgers vector **b** located at the origin of the Cartesian frame ($\mathbf{e}_1$, $\mathbf{e}_2$, $\mathbf{e}_3$) such that $\mathbf{b} = b\mathbf{e}_1 = a_0/\sqrt{2}\mathbf{e}_1$ and the normal vector is $\mathbf{n} = \mathbf{e}_3$.

In this situation, the total eigenstrain tensor given in Eq. (18) has two non-zero components $\epsilon_{13}^0(x) = \epsilon_{31}^0(x) = \frac{b}{2h}\eta(x)$, and the order parameter that characterizes the amount of shear in the slip-plane is such that $\eta(x \to -\infty) = 0$ and $\eta(x \to +\infty) = 1$ as it is shown in Fig. 5. Formally, its kinetic evolution is given by Eq. (22) that here becomes:

$$ \frac{\partial \eta(x, t)}{\partial t} = -L \left( -2\sigma_{13}(x) \frac{\partial \epsilon_{13}^0(x, t)}{\partial \eta(x, t)} + 2\pi A \sin[2\pi\eta(x, t)] - B \frac{\partial^2 \eta(x, t)}{\partial x^2} \right), \quad (27) $$

with the stress component $\sigma_{13} = C_{1313}(\epsilon_{13} - \epsilon_{13}^0)$.

<figure><img src="image_6.png" /><figcaption>Figure 5. Profile of the order parameter $\eta_0$ from Eq. (29) in full line and normalized Ginzburg-Landau energy density $f_{GL}/f_{GL}^0$ in dashed line – calculated with the solution of $\eta_0$. In this graph, the values are normalized to unity taking $\xi = d$ and $f_{GL}^0 = 4A$. The total Ginzburg-Landau energy given Eq. (19) is the surface area of the filling zone. This quantity also corresponds to the core energy of the straight edge dislocation.</figcaption></figure>

To determine the analytic function that describes the profile of $\eta$, we analyze the static case by taking $\partial\eta/\partial t = 0$, assuming that the elastic contributions can be neglected ($\sigma_{13} \approx 0$). In this case, Eq. (27) reduces to the sine-Gordon equation [33]:

$$ 2\pi A \sin[2\pi\eta(x)] - B \frac{\partial^2\eta(x)}{\partial x^2} = 0. \quad (28) $$
whose solution is the soliton $\eta_0$ given by:

$$ \eta_0(x) = \frac{2}{\pi} \arctan \left( \exp \left( \frac{x}{\xi} \right) \right), \qquad (29) $$

with:

$$ \xi = \frac{1}{2\pi} \sqrt{\frac{B}{A}}, \qquad (30) $$

the length characterizing the spreading of the order parameter at point $x = 0$, as it is shown in Fig. 5. Since $\eta$ reflects the amount of shear in the slip-plane, $\xi$ also characterizes the size of the dislocation core.

By introducing this solution in $\mathcal{F}_{GL}$ given in Eq. (19), and integrating within the slip-plane of thickness $h$, one can determine the Ginzburg-Landau non-elastic energy associated with the profile of $\eta$ which is nothing but the core energy (per unit length) of the straight edge dislocation:

$$ W_{core} = \frac{4h}{\pi} \sqrt{AB}. \qquad (31) $$

The amplitudes $A$ and $B$ are thus related to $\xi$ and $W_{core}$, the size and the energy of the dislocation core respectively. In a very generic way, the transition region described by a PFMD is diffuse at the space-scale of the numerical grid while the other parameters are chosen so that the core energy and dislocation velocity are physically correctly reproduced. In this work, we decide to consider $W_{core} = 0.1\mu b^2 = 0.05\mu a_0^2$ as a realistic physical value [34, 35] and $\xi = 0.5d$. With this parametrization, the dislocation core is diffuse enough so that the FCC numerical grid does not introduce lattice resistance on dislocations. In other words, Peierls stress is assumed negligible. This leads to the reduced quantities:

$$ \tilde{A} = \frac{s^2 A}{K} = \frac{\sqrt{3}}{80} \tilde{\mu} \quad \text{and} \quad \tilde{B} = \frac{B}{K a_0^2} = \pi^2 \tilde{A}, \qquad (32) $$

with $\tilde{\mu} = \mu/K$. Interestingly, we note that these value are independent of the scale factor $s$ which only modifies the definition of reduced energies, strains, stresses, displacement-gradients and time unit.

Let us now consider the situation where an external stress is applied such that it results in a resolved shear stress $\sigma_r \equiv \sigma_{13}$ in the dislocation slip-plane. In this case, the dislocation glides and reaches a steady-state motion which is characterized, in the PFMD, by the displacement along $x$ of the profile of the order parameter $\eta$. The dislocation velocity $v$ can be calculated by following a standard asymptotic analysis of the phase-field equations (see for example [36, 37]) where the profile is assumed to be unchanged and still well described by the soliton $\eta_0$. In one such situation, $\eta_0$ is shifted over a length $-vt$ such that $\eta_0(x) = \eta_0(x + vt)$. In the Cartesian frame $(e_1, e_2, e_3)$, since $\eta_0$ describes the static equilibrium profile, Eq. (27) reduces to:

$$ v \frac{\partial \eta_0(x)}{\partial x} = 2 L \sigma_r \frac{\partial \varepsilon_{13}^0(x)}{\partial \eta(x)} \Big|_{\eta_0(x)}. \qquad (33) $$

By introducing the eigenstrain $\varepsilon_{13}^0 = \frac{b}{2h}\eta(x)$, multiplying each term with $(\partial\eta_0(x)/\partial x)$ and taking the integral from $-\infty$ to $+\infty$, one obtain the dislocation velocity at the so-called asymptotic limit:

$$ \frac{2h}{L\pi^2\xi} v = b \sigma_r. \qquad (34) $$
This calculation allows us to show that the dissipative formalism of the PFMD naturally leads to a linear relationship between dislocation velocity and the resolved shear stress. Therefore, we can use the model to reproduce any situation where a linear relationship between dislocation velocity and stress is observed. Such a dependence is generally characterized by the relationship:

$$B_{fr} v = b \sigma_r, \quad (35)$$

where $B_{fr}$ is a friction coefficient and b the Burgers vector norm. A direct comparison between Eqs. (34) and (35) yields to:

$$L = \frac{2h}{B_{fr}\pi^2\xi}, \quad (36)$$

which provides us with a way to link the phenomenological mobility coefficient L and the material parameter $B_{fr}$ [38]. The present development also leads to a link between the time unit $t_0$ and material constants. More precisely, using Eq. (36) and the previous parameters, we obtain $t_0 = s^2/(KL) \approx 4.27s^2B_{fr}/K$.

In pure metals, the linear relationship given Eq. (35) is experimentally related to phonon drag mechanisms and $B_{fr}$ is nothing but the phonon drag coefficient [34, 35]. That said, we can also take the example of the γ-phase in nickel-based superalloys where a linear relationship is also experimentally observed. In this case, dislocations glide in a multi-component solid solution and $B_{fr}$ is viewed as a more general friction term that accounts for all the operative dissipative phenomena. In this context, we can obtain an order of magnitude of $t_0 \approx 4.27s^2B_{fr}/K$ by considering that the Burgers vector of the dislocation is about $b \approx 10^{-10}$ m. In Ref. [39], dislocation velocities $v \approx 10^{-7}$ m/s can be extracted with experienced stresses of about $\sigma_r \approx 10^7$ Pa. Using Eq. (35), we can thus estimate $B_{fr} \approx 10^4$ Pa.s. By taking a bulk modulus of about $K \approx 10^{11}$ Pa as in Ref. [3], a value $t_0 \approx s^20.427 \mu s$ is finally obtained. When atomistic space-scales are considered, $s = 1$ and $t_0 \approx 0.427 \mu s$. That said, when atomic-range effects are negligible and larger grid spacing are conceivable ($s > 1$), the characteristic time-scale of the dislocation motion is opportunistically increased with a square dependence.

# 4. Modification of the PFMD

In this section, we point out some limitations in the prior formulations of the model and make some improvements that provide a more physical coupling between the different fields.

## 4.1. On the coupling between $\eta_p$ and $\sigma_{ij}$

According to the way the PFMD is developed, a dislocation loop is characterized by a region in which the order parameter $\eta_p$ is an integer whose value reflects the amount of elementary shear $\beta_{ij/p}$ given Eq. (16). When the dislocation is static, this requirement is fulfilled by considering an eigenstrain $\epsilon_{ij}^0$ which varies linearly with $\eta_p$, as written in Eqs. (17) and (18). However, a problem arises when considering its temporal variations.

For this, we only have to look at Eq. (27) to note that there exists no stationary solution of $\eta$ that can be an integer as long as the stress $\sigma_{13}$ is different from zero. More generally, the order parameters $\eta_p$ are always different from zero
as long as there is a stress field $\sigma_{ij}$. This causes dissipation in the elastic regions which have not yet been sheared by the passage of dislocations, and generates a non-physical dependence between the amount of shear introduced by a dislocation and the stress which is locally experienced.

This spurious dissipation mechanism is automatically removed within the *loopon* formalism proposed in [14], in which the order parameters $\eta_p$ associated to the plastic activity are subjected to automaton-like dynamics only controlled by elastic driving forces with no recourse to a Ginzburg-Landau core energy. In the present model, the dynamics of $\eta_p$ are partially controlled by a Ginzburg-Landau energy which confers to these order parameters a spatially continuous character. This property is at the origin of the artificial dissipation mechanism mentioned above.

To remove this artifact, one strategy is to replace the order parameter $\eta_p$ by a function $f(\eta_p)$ in the definition of the stress-free displacement-gradients:

$$ \beta_{ij}^{0}(\mathbf{r}) = \sum_{p=1}^{12} \beta_{ij/p}^{0} f(\eta_{p}(\mathbf{r})). \quad (37) $$

Within this formulation, the function $f$ must be such that $f(\eta_p) = \eta_p$ when $\eta_p$ is an integer. Then, in Eq. (27), we note that the elastic term becomes proportional to $\partial f(\eta_p)/\partial\eta_p$. The coupling between $\eta_p$ and $\sigma_{ij}$ can therefore be canceled by taking a function so that $\partial f(\eta_p)/\partial\eta_p = 0$ when $\eta_p$ is an integer. In some previous works, we found two functions that possess the aforementioned properties which have been used for this purpose. The first, introduced by Levitas et al. in Refs. [40–42], is given by:

$$ f(\eta_p) = [\eta_p] + \{\eta_p\}^2 (3 \operatorname{sgn}(\eta_p) - 2 \{\eta_p\}), \quad (38) $$

where $[\cdot, \cdot]$ and $\operatorname{sgn}(\cdot)$ denote the integer part, the fractional part and the sign functions, respectively. The second, introduced by Hu et al. in Ref. [43], writes:

$$ f(\eta_p) = \eta_p - \frac{\sin(2\pi\eta_p)}{(2\pi)}. \quad (39) $$

In both cases, it can be checked that the use of $f(\eta_p)$ does not affect the analytical expressions of the phase-field parameters A, B and L determined in subsection 3.3. As discussed in detail in Ref. [42], the first function given Eq. (38) allows us to nucleate dislocations beyond a critical stress, while the function given Eq. (39) definitively cancels the spontaneous generation of any dislocation. Even though the first approach seems to be more realistic, the critical stress of the dislocation nucleation cannot be easily controlled. For this reason, and because nucleation aspects are not considered in the present article, we chose to use the second function $f(\eta_p) = \eta_p - \sin(2\pi\eta_p)/(2\pi)$ which is also numerically more tractable.

## 4.2. On the coupling between several $\eta_p$

Another contradiction arises when we focus on the concomitant temporal evolution of the different order parameters $\eta_p$ that describe three slip-systems of a common slip-plane. To understand this, let us consider the slip-systems $p = 1, 2, 3$ with the slip-plane (111) described in Fig. 4. In this case, only $\eta_1, \eta_2$ and $\eta_3$ can take values different from
zero. As a consequence, the crystalline energy density given Eq. (20) becomes:

$$f_{\text{cryst}} = 2A \{1 - \cos[\pi(\eta_1 - \eta_2)] \cos[\pi(\eta_2 - \eta_3)] \\ \cos[\pi(\eta_3 - \eta_1)]\}, \quad (40)$$

and the temporal evolution of the different $\eta_p$, controlled by the Allen-Cahn equation Eq. (26), involves the following crystalline driving forces:

$$\frac{\partial f_{\text{cryst}}}{\partial \eta_1} = 2\pi A \cos[\pi(\eta_2 - \eta_3)] \sin[\pi(2\eta_1 - \eta_2 - \eta_3)], \quad (41)$$

$$\frac{\partial f_{\text{cryst}}}{\partial \eta_2} = 2\pi A \cos[\pi(\eta_3 - \eta_1)] \sin[\pi(2\eta_2 - \eta_3 - \eta_1)], \quad (42)$$

$$\frac{\partial f_{\text{cryst}}}{\partial \eta_3} = 2\pi A \cos[\pi(\eta_1 - \eta_2)] \sin[\pi(2\eta_3 - \eta_1 - \eta_2)]. \quad (43)$$

Now we imagine a system with only one moving dislocation characterized by $\eta_1$. The other parameters $\eta_2$ and $\eta_3$ are initially equal to zero. Almost everywhere, the order parameter is $\eta_1 = 0$ or 1, except in the diffuse transition regions in which $\eta_1 \in [0, 1]$. According to Eq. (41), almost everywhere, the crystalline driving force of $\eta_1$ is equal to zero, except in the transition region. Herein, since the temporal derivative of $\eta_1$ is proportional to $\partial f_{\text{cryst}}/\partial\eta_1$, the value is changed, which actually contributes to the dislocation motion.

The problem in this region is that, according to Eqs. (42) and (43), values of $\eta_1$ different from zero also contribute to the crystalline driving forces of $\eta_2$ and $\eta_3$. As a consequence, the two other order parameters take a value different from zero in the transition region of the first one. This leads to the spontaneous generation of the coupled slip-systems and for one order parameter it is impossible to proceed without the two others.

Although the formulation of crystalline energy is relevant when it is to consider the interaction of preexistent dislocations, it appears that it can generate nonphysical situations when we consider their dynamics. Notably: (i) kinematics are perturbed by the unexpected changes in the profile of a given $\eta_p$, and (ii) the activation of other slip-systems can spontaneously generate dislocations without justification. Even though this mechanism can be physically relevant, its correct reproduction within a continuous PFMD formulation would require at least a mesh with a sub-atomic grid spacing. Conversely, it may lead to spurious effects when using a coarser grid, which is the aim of the present PFMD formulation.

Therefore, since we expect that the dynamics of a given $\eta_p$, controlled by Eq. (22), be executed only if an integer value preexists in the system, we decide to introduce a test algorithm that activates it only when the order parameter on a given octahedral site is different from the values of its neighbors. Formally, this is the same as performing a test on the gradient of $\eta_p$ and verifying that its absolute value is higher than zero. Because values very close to zero are expected to be computed as long as a the corresponding slip-systems take place somewhere in the system, we chose to compare the gradient to a threshold value $\nabla\eta_p^{\text{thrs}}$ not much higher than zero and set the order parameter to the closest integer when the gradient is below. In our simulations, the absolute value of the gradient is compared to its theoretical
value calculated at a distance $x = 3\xi$ from the center of the profile. Taking the spatial derivative of $\eta$ as it is given in Eq. (29), one obtains $\partial\eta(x)/\partial x = 1/(\pi\xi\cosh(x/\xi))$ and thus $\nabla\eta_p^{\text{thrs}} = 0.0316/\xi = 0.0632d^{-1}$ with $\xi = 0.5d$.

It is worth mentioning that this test only activates the dynamics within the core transition regions, as displayed in Fig. 5, instead of the whole system. As a consequence, in a material that contains only a few dislocations, the computational time required for this test is compensated by the time which is saved by avoiding the calculation of the dynamics in every octahedral sites.

# 5. Validation of the present model

To validate the model and show its consistency with the analytical developments, we have simulated a single dislocation loop collapsing under its own stress field, two interacting dislocations forming junctions and three-dimensional loops emerging from a Frank-Read source. In this section, simulations are performed on the atomic-scale considering $d = a_0$ which corresponds to $s = 1$.

## 5.1. One dislocation loop

We start with the simulation of a single dislocation loop collapsing under its own stress field. To do so, we have considered a $100d^3$ FCC isotropic material of shear modulus $\mu = 0.40K$ with periodic boundary conditions. In the $(1\bar{1}1)$ slip-plane, we have introduced a circular dislocation loop of initial radius $R_0 = 40d$ with the Burgers vector $\mathbf{b}_1 = a_0[110]/2$. The phase-field amplitudes are chosen to be $\tilde{A} = \sqrt{3}\tilde{\mu}/80 = 0.0087$ and $\tilde{B} = \pi^2\tilde{A} = 0.085$ in agreement with the static theoretical expressions that give $W_{\text{core}} = 0.1\mu b^2$ and $\xi = 0.5d$ for the energy and the size of the dislocation core respectively.

In Fig. 6, we have extracted at a given time the simulated values of the order parameter $\eta_1$ and the Ginzburg-Landau energy density $f_{GL}$ as a function of the normalized coordinate $x/d$ along the [110] direction. These values are compared to the static theoretical ones obtained in subsection 3.3 considering $\tilde{A} = 0.0087$ and $\tilde{B} = 0.085$.

First, as expected, $\eta_1$ is found to be strictly equal to an integer value as long as it does not describe the diffuse transition region. This is due to the function $f(\eta)$ given in Eq. (39), introduced in subsection 4.1 to cancel the dependence between the amount of plastic shear and the stress field. Similarly, in consistency with the modified algorithm proposed in subsection 4.2, it has been checked (but not shown here) that the coupled order parameters $\eta_2$ and $\eta_3$ remain equal to zero during the kinetic evolution of $\eta_1$.

Inside the transition region, the order parameter $\eta_1$ is found to be slightly different from $\eta_0$ as given in Eq. (29). This difference is due to the fact that $\eta_0$ is determined in the static case by neglecting the elastic contributions. That said, it is interesting to note that the spreading of the core is close in both cases, confirming that its order of magnitude is about a few $\xi = 0.5d$.

Note that even though the profile of $\eta_1$ in the transition region may reflect a dissociated core, it still describes a perfect core in consistency with the present PFMD. In fact, the core structure would have to be characterized by the
<figure><img src="image_7.png" /><figcaption>Figure 6. Order parameter $\eta_1$ and Ginzburg-Landau energy density $f_{GL}$ as a function of $x/d$. The values extracted from the simulations are compared to $\eta_0$ given in Eq. (29) and $f_{GL}$ that can be calculated in subsection 3.3. The Ginzburg-Landau energy densities are normalized considering the maximal static theoretical value $f_{GL}^0 = 0.0348Kd^2$.</figcaption></figure>

disregistry vector [44] which here, is not proportional to $\eta_1$ but to the function $f(\eta_1)$ given Eq. (39). Its effect is to make the disregistry profile more regular than the $\eta_1$ profile.

If we now focus on the simulation's value of the Ginzburg-Landau energy density $f_{GL}$, we note that the pic is concomitantly flattened and widened compared to the one given by the static approximated analytical solution. However, the area below the curve that gives the core energy of the dislocation is not significantly modified. Therefore, we found a simulated value of about $W_{core} = 0.13\mu b^2$ which is actually close to the static theoretical one ($W_{core} = 0.1\mu b^2$).

To investigate the kinematic aspects of our model, we chose to study the temporal evolution of the surface area of the dislocation loop following [14]. For this, we start to write its elastic energy $W$ as it can be calculated within the framework of linear elasticity [34]:

$$W(R) = \frac{\mu b^2 R}{4} \frac{2 - v}{1 - v} \left( \ln \left( \frac{4R}{r_0} \right) - 2 \right), \quad (44)$$

where $R$ is the radius of the loop, $b$ its Burgers vector norm, $\mu$ the shear modulus and $v$ is the Poisson ratio. In this equation, $r_0 \ll R$ denotes the cut-off distance enclosing the core region below which the continuous description of elasticity does not remain valid. In the present model, $r_0$ is expected to be equivalent to the core size of the dislocation, which is characterized in Fig. 6 by the length of the diffuse transition region.

During the collapse, the stress experienced by a line element of the loop is given by $\sigma = -W'(R)/(2\pi b)$ [34]. Since the velocity of the line element is given by Eq. (35), we can show, using Eq. (44), that the temporal evolution of $S(t) = \pi R^2(t)$ must finally satisfy:

$$\dot{S}(t) = - \frac{\mu b^2}{4B_{fr}} \frac{2 - v}{1 - v} \left( \ln \left( \frac{4\sqrt{S(t)/\pi}}{r_0} \right) - 1 \right), \quad (45)$$

where $B_{fr} = 2h/(L\pi^2\xi)$ is the viscous coefficient given in Eq. (36). In this equation, the only free parameter is the cut-off $r_0$ since $b = a_0/\sqrt{2}$, $v = 0.33$ and $B_{fr} = 2.34Kt_0$. Therefore, we can estimate $S$ by the numerical integration
of Eq. (45) considering a reasonable value of $r_0 \ll R$.

In Fig. 7, the normalized surface area $S/S_0$ of the loop has been extracted from the simulation and plotted as a function of the reduce time $t/t_0$. The value is compared to the theoretical one resulting from Eq. (45) considering $r_0 = 2.05d$, which has been found to be the best fitting parameter.

<figure><img src="image_8.png" /><figcaption>Figure 7. Theoretical and simulated normalized surface area of the loop $S/S_0$ versus time $t/t_0$. The insert shows different snapshots of the dislocation loop of Burgers vector $\mathbf{b}_1 = a_0[0\bar{1}\bar{1}]/2$ during its collapse.</figcaption></figure>

First, we note that $r_0$, which characterizes the size of the core, is consistent with the one that would be found in Fig. 6. Then, we confirm that the temporal evolution of $S/S_0$ extracted from the simulation reproduces the theoretical normalized value found with linear elasticity. The remaining discrepancies can be explained by the fact that, during its collapse, the dislocation loop does not remain circular, as it is shown in the insert of Fig. 7. This is because the line tension of the loop is not really isotropic, as supposed in the analytical integration, but higher along its screw component.

## 5.2. Self and coplanar interactions

In the simulation box considered previously, we introduce two dislocation loops of Burgers vectors $\mathbf{b}_3 = a_0[0\bar{1}\bar{1}]/2$ in the same slip-plane $(1\bar{1}1)$. They are initially located and oriented as exposed in Fig. 8a. By applying the stress:

$$ \sigma = \frac{\tau}{3\sqrt{2}} \begin{pmatrix} 2 & \bar{2} & \bar{1} \\ \bar{2} & 2 & 1 \\ \bar{1} & 1 & \bar{4} \end{pmatrix}, \qquad (46) $$

we obtain in the plane $(1\bar{1}1)$ a shear stress $\tau > 0$ (aligned with $[1\bar{1}2]$) that provokes the expansion of the loops. As shown on Fig. 9b, the loops react and annihilate. This reaction is called the self annihilation. It results in a larger loop of Burgers vectors $\mathbf{b}_3 = a_0[0\bar{1}\bar{1}]/2$.

Then, reviving the work performed in Ref. [27], we consider two dislocation loops with dissimilar Burgers vectors $-\mathbf{b}_3 = a_0[011]/2$ and $-\mathbf{b}_2 = a_0[10\bar{1}]/2$, introduced and orientated as displayed in Fig. 9a. They grow under the
<figure><img src="image_9.png" /><figcaption>Figure 8. Self-reaction annihilation of dislocation loops with Burgers vector **b**₃ = a₀[0**1**0]/2 in the slip-plane (1**1**1). (a) Initial configuration. (b) Reaction under the planar stress **τ**.</figcaption></figure>

applied stress σ given in Eq. (46). As shown in Fig. 9b, the dislocations meet themselves and react to form a coplanar junction of Burgers vector −**b**₂ − **b**₃ = **b**₁ = a₀[110]/2. Because τ is perpendicular to this new line segment, the latter remains static. In consistency with the formulation of the crystalline energy density given Eq. (20), the core energy of the junction is the same as the original dislocation segments (not the double as in the first historical models [10, 11]). We therefore confirm that the model correctly reproduces the coplanar dislocation reactions in an FCC material.

<figure><img src="image_10.png" /><figcaption>Figure 9. Coplanar reaction of dislocation loops with Burgers vector −**b**₃ = a₀[011]/2 and −**b**₂ = a₀[10**1**]/2 in the slip-plane (1**1**1). (a) Initial configuration. (b) Reaction under the planar stress **τ** resulting in the formation of a junction of Burgers vector **b**₁ with identical core energy.</figcaption></figure>

5.3. *Interaction in two intersecting slip-planes*

We now consider the interaction of dislocations gliding in two intersecting slip-planes. First, a dislocation loop of Burgers vector −**b**₃ = a₀[011]/2 is introduced in its (1**1**1) slip-plane while another of Burgers vector −**b**₄ = a₀[**1**0**1**]/2 is first localized in its own (1**1**1) slip-plane. They are both oriented as shown in Fig. 10a. Under an applied stress
$\sigma_{33} < 0$, the loops undergo the same positive resolved shear stress $\sigma_r = \sigma_{33} / \sqrt{6}$ that induces their expansion. In Fig. 10b, the dislocations react to form a junction of Burgers vector $-\mathbf{b}_3 - \mathbf{b}_4 = -\mathbf{b}_7 = \mathbf{b}_{11} = a_0[\overline{1}10]/2$. Since this Burgers vector does not correspond to any of the original slip-systems, the junction remains stationary, constituting a sessile Lomer lock [35].

<figure><img src="image_11.png" /><figcaption>Figure 10. Sessile Lomer lock reaction of dislocation loops with Burgers vector $-\mathbf{b}_3 = a_0[011]/2$ and $-\mathbf{b}_4 = a_0[\overline{1}0\overline{1}]/2$ in the slip-plane $(1\overline{1}1)$ and $(\overline{1}11)$ respectively. (a) Initial configuration. (b) Reaction under a stress $\sigma_{33} < 0$ resulting in the formation of a junction of Burgers vector $-\mathbf{b}_7 = \mathbf{b}_{11} = a_0[\overline{1}10]/2$.</figcaption></figure>

Then, two dislocation loops of Burgers vectors $-\mathbf{b}_1 = a_0[\overline{1}\overline{1}0]/2$ and $-\mathbf{b}_5 = a_0[110]/2$ lying in their own slip-planes $(1\overline{1}1)$ and $(\overline{1}11)$ respectively, have been introduced and orientated as shown in Fig. 11a. By applying a stress $\sigma_{11} > 0$, the loops expand under the positive resolved shear stress $\sigma_r = \sigma_{11} / \sqrt{6}$. At their intersection, the loops annihilate ($-\mathbf{b}_1 - \mathbf{b}_5 = 0$). This reaction is called the colinear annihilation.

<figure><img src="image_12.png" /><figcaption>Figure 11. Colinear annihilation reaction of the dislocation loops with Burgers vector $-\mathbf{b}_1 = a_0[\overline{1}\overline{1}0]/2$ and $-\mathbf{b}_5 = a_0[110]/2$ in the slip-plane $(1\overline{1}1)$ and $(\overline{1}11)$ respectively. (a) Initial configuration. (b) Annihilation under a stress $\sigma_{11} > 0$.</figcaption></figure>
Physically, the colinear annihilation leads to a perfect crystal that does not generate any kind of energy. In the simulations however, as in every PFMD considered thus far [10–14], the intersecting line resulting from the annihilation reaction still exhibits a dipole of dislocations that generates a small non-zero stress field, as well as a core energy due to the formulation of the Ginzburg-Landau free energy. Because the range of the dipole stress field is much lower than the one generated by an isolated dislocation, its effects are assumed to be negligible. As for the residual core energy, one should propose a non trivial development of the PFMD which is beyond the scope of the present article.

5.4. Frank-Read source

To simulate a Frank-Read source, we decide to introduce a dislocation segment of type $\mathbf{b} = a_0[\bar{1}\bar{1}0]/2$ with pinned extremities in the simulation box. For this, a static rectangular loop of Burgers vector $-\mathbf{b}_5 = a_0[110]/2$ has been considered in the (\bar{1}11) slip-plane, supplemented by a small embryo of loop of Burgers vector $-\mathbf{b}_1 = a_0[\bar{1}\bar{1}0]/2$ in the (1\bar{1}1) plane. Due to the annihilation of the two loops along the MN segment, the dynamics on the order parameter $\eta_1$ allow us to reproduce the bowing of the screw dislocation segment MN in the (1\bar{1}1) slip-plane. By applying a stress $\sigma_{11} > 0$, the MN dislocation line of screw character bows and forms an extended loop of Burgers vector $-\mathbf{b}_1 = a_0[\bar{1}\bar{1}0]/2$ as is shown in Fig. 12. After one cycle, another loop is spontaneously nucleated, still activating the Frank-Read source. One such simulation allows us to confirm that the model can reflect realistic dislocation multiplication mechanisms.

<figure><img src="image_13.png" /><figcaption>Figure 12. Superimposed snapshots showing the Frank-Read source during subsequent equidistant time moments under the applied stress $\sigma_{11} = 0.08\text{kN/m}$.</figcaption></figure>
# 6. Applications

In this section, we present some applications that demonstrate the ability of the model to describe, at different space-scales, dislocations operating in the presence of heterogeneities made of voids and free surfaces.

## 6.1. Traction of a cylindrical nanopillar

First, we simulate the traction of a cylindrical nanopillar [45] that contains a source of spiral dislocations provoking the shear of the system. To be more concrete, we assume that the material is an aluminium single crystal with bulk modulus $K \approx 80$ GPa, shear modulus $\mu = 0.40K$ and lattice parameter $a_0 = d = 0.405$ nm (corresponding to $s = 1$).

Initially, a $200 \times 100 \times 100d^3$ system is generated. Then, by taking the pillar aligned with the direction of traction [100], we define the non-material parts (of zero elastic and phase-field coefficients) outside a cylindrical region generated by a disc of radius $R_0 = 49d$ and normal vector [100]. By doing so, we obtain a cylindrical pillar with lateral free surfaces oriented perpendicularly to the direction [100]. Along this direction, periodic boundary conditions are still considered.

A spiral source is subsequently introduced by reproducing the procedure described in subsection 5.4. Because of the presence of the lateral free surfaces, the rectangular loop of Burgers vector $-b_5$ stands partially outside the nanopillar, forming a step made of the intersecting section between its slip-plane and the free surface. Inside the material, the upper screw segment of the loop bows under a uniaxial stress $\sigma_{11} > 0$, forming a spiral dislocation of Burgers vector $-b_1 = a_0[\bar{1}\bar{1}0]/2$ in the $(1\bar{1}1)$ slip-plane. One such loop appears as a spiral because its second extremity ends on the free surface which does not constitute a pinned point. Through its $N$ rotations, the spiral dislocation induces the shear of the nanopillar, forming a step of amplitude $Nb$ along its entire surface section.

In Fig. 13, different snapshots of the systems are exposed during subsequent equidistant time moments. It is shown how the spiral dislocation moves and generates the step on the free surface (visible in Fig. 13.d). During its motion, the dislocation line intersects its initial position (in red, see color on-line) that causes the formation of a small discontinuity. This is actually due to the diffuse formulation of the PFMD which does not allow to entirely remove the Ginzburg-Landau energy along this line. That said, as shown in Fig. 13, this does not modify the global shape of the spiral and does not change the consequence of its movement. Above all, we confirm that the coupling between dislocations and free surfaces is naturally considered within the present model.

Instead of directly considering a non-material part outside the cylinder, one could also imagine an intermediate material with different elastic coefficients, playing the role of a shell with a cylindrical ring shape. In this case, the present PFMD would be able to reproduce a nanopillar whose surface properties have been changed (oxidation, elaboration processes...) and its effects could be systematically investigated. A fracturing criterion could also be introduced, by following the strategy exposed in Ref. [38] for instance, which would provide a way of studying the brittle to ductile transition. Of course, these examples are beyond the scope of the present article.
<figure><img src="image_14.png" /><figcaption>Figure 13. Traction of cylindrical nanopillar under uniaxial stress $\sigma_{11} > 0$. (a) Initial stages where the upper segment of type b₁ bows and generates an external step highlighted by a dashed line. The presence of another static step formed by the partial rectangular loop of Burgers vector $-b_5$ may also be noted. (b,c,d) The bowed segment propagates, extending the step and taking the form of a spiral due to the presence of the free-surface. (e) External appearance of the pillar after the spiral has rotated three times.</figcaption></figure>

## 6.2. Dislocations interacting with a spherical pore

For the last application, we drew on the work of Refs. [3, 46] which consists in studying the effect of plasticity on the pore closure in nickel-based superalloys during hot isostatic pressing (at $T = 1288^{\circ}$C). In this case, the material is elastically strongly anisotropic. Of course, in the present article, we restrict our analysis to a simple example whose purpose is to illustrate the potentiality of the model instead of achieving a complete investigation on the pore-dislocations interactions.

Therefore, we have reconsidered the system exposed in subsection 2.3 that consists of a $100d^3$ FCC material under periodic boundary conditions, containing a central spherical pore of radius $R_0 = 15d$. Due to the cubic elastic behavior, the elastic energy density now writes:

$$ f_{el} = \frac{K}{2}(\varepsilon_{11} + \varepsilon_{22} + \varepsilon_{33})^2 + \frac{C'}{3}[(\varepsilon_{11} - \varepsilon_{22})^2 + (\varepsilon_{11} - \varepsilon_{33})^2 + (\varepsilon_{22} - \varepsilon_{33})^2] + \frac{\mu}{2}[(2\varepsilon_{12})^2 + (2\varepsilon_{13})^2 + (2\varepsilon_{23})^2], \quad (47) $$

where $K = 164.9$ GPa is the bulk modulus and $C' = 0.09K$, $\mu = 0.42K$ are the two shear modulus found in [3]. The strong anisotropy can be highlighted by the Zener anisotropy ratio $\mu/C' = 4.6$. For this material, we consider the lattice parameter of nickel $a_0 = 0.352$ nm in a context where short-range dislocation interactions are not involved.
This enables us to work at the micron-scale by considering $s = d/a_0 = 690$ which corresponds to $d = 234$ nm.

Then, as it is shown in Fig. 14.a, nine dislocation loops of radius $R = 45d$ and Burgers vector $b_1 = a_0[\bar{1}\bar{1}0]/2$ are introduced into nine different (1\bar{1}1) slip-planes. They are stacked and centered along the diagonal of direction [1\bar{1}1].

<figure><img src="image_15.png" /><figcaption>Figure 14. [0\bar{1}0] and [1\bar{1}1] views of a spherical pore of radius $R_0 = 15d \approx 3.6 \ \mu\text{m}$ in a $100d^3 \approx 24.3 \ \mu\text{m}^3$ FCC material surrounded by nine dislocations of Burgers vector $b_1$ lying in nine different (1\bar{1}1) slip-planes. They are stacked and centered along the green diagonal of direction [1\bar{1}1] (see colors on-line). The dislocation loop whose slip-plane intersects the equatorial plane of the pore is colorized in red. (a) Initial configuration. (b,c) Evolution during subsequent moments in time. (d) Final compressed equilibrium configuration.</figcaption></figure>

The effect of the pore on the dislocations' behavior is observed under isostatic pressure $p = 6.10^{-4} K \approx 100$ MPa which corresponds to the loading involved in Refs. [3, 46]. The snapshots exposed in Fig. 14 show the evolution of dislocations as the system reaches its compressed equilibrium configuration.

It is found that the dislocation loops adopt a non-trivial equilibrium shape and become stacked in a way that is influenced by the anisotropic stress-field of the pore (and the other dislocations). Notably, for each of the loops, we observe that one of its two edge parts moves from its initial position to a position consisting in the intersection between
its slip-plane and the diagonal of direction [11̄1]. At the equilibrium, this results in a pattern of dislocations described in Fig. 14.d for which the center of the pore is a point of symmetry and where edge parts are aligned with the diagonal.

For our purpose, it is interesting to focus on the equatorial dislocation, colorized in red in Fig. 14, which is observed to shrink more rapidly than the other loops while a part of its line is captured by the pore through the spherical free surface. At the end of the simulation, this loop results in two smaller loops partially attached to the pore, leading to the formation of a surface step (not shown in the figure). To our knowledge, this typical behavior can neither be expected with simple analytical calculations nor easily modeled with already existing numerical approaches.

Within the present description, all mechanical effects are automatically considered. They range from the stress-field generated by the pore inside the material domain to the image force experienced by the dislocations due to the presence of the surface. They also include the effects of the other dislocations and those of the steps generated on the pore surface. As a consequence, one such model would be of great interest in the study of the microscale mechanisms involved during hot isostatic pressing of FCC materials. This is an ongoing study.

7. Conclusion

In this article, we first expose a new finite difference scheme that allows us to rigorously discretize any continuous
fields on an explicit face-centered cubic (FCC) grid. In this scheme, multi-component physical quantities (vectors,
tensors...) are all localized on the same computational node which enables us to consider any kind of elastic hetero-
geneity (precipitates, voids, cracks...) and any complex specimen shape.

Then, we introduce a classical phase-field model of dislocations (PFMD) considering the slip-systems of the FCC
geometry. A particular attention is paid to the dislocation core behaviors in order to remove potentially spurious effects
that are present in previous PFMDs when the chosen numerical resolution is not subatomic. We also establish the link
between the mobility parameter *L* associated to the phase-field dynamics of the dislocation fields η<sub>p</sub> and experimental
friction coefficients *B*<sub>fr</sub> that control the viscous glide of dislocations.

The model is tested in a few simple cases (dislocation collapse, junctions, Franck-Read source) and illustrated
in more concrete applications, standing at different space-scales, where dislocations operate in the presence of mi-
crostructural heterogeneities, such as free-surfaces (cylindrical nanopillar) and voids (pore under isostatic pressure).
As a perspective of the present work, we plan to use this model to investigate pore closure by plastic flow under
hot isostatic pressing [3]. Concerning the numerics, we also plan to improve the mechanical solver in the context
of heterogenous materials by considering a Fourier space formulation and/or minimizing the elastic fields by using
non-linear conjugate gradient methods.

Finally, thanks to the variational formulation of the PFMD framework, its extension to the coupling with an
evolving microstructure is straightforward (provided all the phases display an FCC crystallography), notably because
the present discretization scheme does not depend on whether the phase boundaries are mobile are not. Of course, the
present approach could also be the base of a lot of numerical models where considering the intrinsic FCC symmetry
is of paramount importance.

Acknowledgments and founding

We gratefully acknowledge the German Scientific Research Foundation DFG and the financial support from the French National Research Agency ANR (grants ANR15-MERA-0003-03). We would also like to thank A. Epishin, B. Fedelich, T. Feldmann, A. Charmi, B. Viguier and D. Poquillon for fruitful discussions concerning the pore– dislocations coupling in subsection 6.2 of this paper.

Appendix A. Discretization of the spatial differential operators

In Fig. A.15, we have reproduced the FCC unit cell displayed in Fig. 2 by conserving the octahedral sites only. At the center, the order parameter η ( a , b , c ) is the one defined on the first SC grid, indexed in space with the integers a , b and c (removed in what follows, for clarity). Its neighboring values η ( a$_{n}$ , b$_{n}$ , c$_{n}$ ) are defined on the twelve closest octahedral sites such that a$_{n}$ = a ( ± 1), b$_{n}$ = b ( ± 1) and c$_{n}$ = c ( ± 1). They are colorized according to their underlying SC grid but indexed independently from it (see colors on-line). In the reference cubic frame, the coordinates of the nearest neighbor vectors are Δ X$_{1}$ = ( a$_{n}$ − a ) / 2, Δ X$_{2}$ = ( b$_{n}$ − b ) / 2 and Δ X$_{3}$ = ( c$_{n}$ − c ) / 2, respectively, where d still designates the FCC lattice parameter.

<figure><img src="image_16.png" /><figcaption>Figure A.15. FCC unit cell where the octahedral sites supporting the order parameters η are identified by octahedrons colorized according to their underlying SC grid (see colors on-line).</figcaption></figure>

To calculate the spatial differential operators of η on its own octahedral site, the order parameters η ( a$_{n}$ , b$_{n}$ , c$_{n}$ ) are expressed as a function of η using the Taylor series:

$$
\begin{align}
\eta(a_n, b_n, c_n) &= \eta + \sum_i \Delta X_i \frac{\partial \eta}{\partial X_i} + \sum_i \sum_j \frac{\Delta X_i \Delta X_j}{2} \frac{\partial^2 \eta}{\partial X_i \partial X_j} \nonumber \\
&\quad + O(d^3). \tag{A.1}
\end{align}
$$
The next step consists in performing a linear combination of the Taylor series that conserves the components of interest and removes the others. This is how we obtain the discretized form of the gradient components $\partial\eta/\partial X_i$:

$$ \frac{\partial \eta}{\partial X_1} = \frac{1}{4d} [\eta(a+1, b+1, c) + \eta(a+1, b-1, c) \\ - \eta(a-1, b+1, c) - \eta(a-1, b-1, c) \\ + \eta(a+1, b, c+1) + \eta(a+1, b, c-1) \\ - \eta(a-1, b, c+1) - \eta(a-1, b, c-1)] + O(d^2), \quad (A.2) $$

$$ \frac{\partial \eta}{\partial X_2} = \frac{1}{4d} [\eta(a+1, b+1, c) + \eta(a-1, b+1, c) \\ - \eta(a+1, b-1, c) - \eta(a-1, b-1, c) \\ + \eta(a, b+1, c+1) + \eta(a, b+1, c-1) \\ - \eta(a, b-1, c+1) - \eta(a, b-1, c-1)] + O(d^2), \quad (A.3) $$

$$ \frac{\partial \eta}{\partial X_3} = \frac{1}{4d} [\eta(a+1, b, c+1) + \eta(a-1, b, c+1) \\ - \eta(a+1, b, c-1) - \eta(a-1, b, c-1) \\ + \eta(a, b+1, c+1) + \eta(a, b-1, c+1) \\ - \eta(a, b+1, c-1) - \eta(a, b-1, c-1)] + O(d^2), \quad (A.4) $$

and the discretized form of the Hessian components $\partial^2\eta/(\partial X_i \partial X_j)$:

$$ \frac{\partial^2 \eta}{\partial X_1^2} = \frac{1}{d^2} [\eta(a+1, b+1, c) + \eta(a+1, b-1, c) \\ + \eta(a-1, b+1, c) + \eta(a-1, b-1, c) \\ + \eta(a+1, b, c+1) + \eta(a+1, b, c-1) \\ + \eta(a-1, b, c+1) + \eta(a-1, b, c-1) \\ - \eta(a, b+1, c+1) - \eta(a, b+1, c-1) \\ - \eta(a, b-1, c+1) - \eta(a, b-1, c-1) \\ - 4\eta] + O(d^2), \quad (A.5) $$

$$ \frac{\partial^2 \eta}{\partial X_2^2} = \frac{1}{d^2} [\eta(a+1, b+1, c) + \eta(a+1, b-1, c) \\ + \eta(a-1, b+1, c) + \eta(a-1, b-1, c) \\ + \eta(a, b+1, c+1) + \eta(a, b+1, c-1) \\ + \eta(a, b-1, c+1) + \eta(a, b-1, c-1) \\ - \eta(a+1, b, c+1) - \eta(a+1, b, c-1) \\ - \eta(a-1, b, c+1) - \eta(a-1, b, c-1) \\ - 4\eta] + O(d^2), \quad (A.6) $$
$$ \frac{\partial^2 \eta}{\partial X_3^2} = \frac{1}{d^2} [\eta(a, b + 1, c + 1) + \eta(a, b + 1, c - 1) \\ + \eta(a, b - 1, c + 1) + \eta(a, b - 1, c - 1) \\ + \eta(a + 1, b, c + 1) + \eta(a + 1, b, c - 1) \\ + \eta(a - 1, b, c + 1) + \eta(a - 1, b, c - 1) \\ - \eta(a + 1, b + 1, c) - \eta(a + 1, b - 1, c) \\ - \eta(a - 1, b + 1, c) - \eta(a - 1, b - 1, c) \\ - 4\eta] + O(d^2), \qquad (A.7) $$

$$ \frac{\partial^2 \eta}{\partial X_1 \partial X_2} = \frac{1}{d^2} [\eta(a + 1, b + 1, c) + \eta(a - 1, b - 1, c) \\ - \eta(a + 1, b - 1, c) - \eta(a - 1, b + 1, c)] \\ + O(d^2), \qquad (A.8) $$

$$ \frac{\partial^2 \eta}{\partial X_1 \partial X_3} = \frac{1}{d^2} [\eta(a + 1, b, c + 1) + \eta(a - 1, b, c - 1) \\ - \eta(a + 1, b, c - 1) - \eta(a - 1, b, c + 1)] \\ + O(d^2), \qquad (A.9) $$

$$ \frac{\partial^2 \eta}{\partial X_2 \partial X_3} = \frac{1}{d^2} [\eta(a, b + 1, c + 1) + \eta(a, b - 1, c - 1) \\ - \eta(a, b + 1, c - 1) - \eta(a, b - 1, c + 1)] \\ + O(d^2). \qquad (A.10) $$

As a matter of fact, these formulas correspond to the FCC stencil used in the acoustics community to perform the finite difference approximations of the Laplacian (Hessian trace) [47].

## References

[1] M. Yaguchi, M. Yamamoto, T. Ogata, A viscoplastic constitutive model for nickel-base superalloy, part 1: kinematic hardening rule of anisotropic dynamic recovery, Int. J. Plasticity 18 (8) (2002) 1083–1109.

[2] F. Roters, P. Eisenlohr, L. Hantcherli, D. Tjahjanto, T. Bieler, D. Raabe, Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications, Acta Mater. 58 (4) (2010) 1152–1211.

[3] A. Epishin, B. Fedelich, T. Link, T. Feldmann, I. L. Svetlov, Pore annihilation in a single-crystal nickel-base superalloy during hot isostatic pressing: Experiment and modelling, Mat. Sci. Eng. A 586 (2013) 342–349.

[4] J. Kundin, D. Raabe, H. Emmerich, A phase-field model for incoherent martensitic transformations including plastic accommodation processes in the austenite, J. Mech. Phys. Sol. 59 (10) (2011) 2082–2102.

[5] M. Cottura, B. Appolaire, A. Finel, Y. Le Bouar, Coupling the phase field method for diffusive transformations with dislocation density-based crystal plasticity: Application to Ni-based superalloys, J. Mech. Phys. Sol. 94 (2016) 473–489.

[6] P.-L. Valdenaire, Y. Le Bouar, B. Appolaire, A. Finel, Density-based crystal plasticity: From the discrete to the continuum, Phys. Rev. B 93 (2016) 214111.

[7] I. Groma, M. Zaiser, P. D. Ispánovity, Dislocation patterning in a two-dimensional continuum theory of dislocations, Phys. Rev. B 93 (2016) 214110.
[8] A. Vattré, B. Devincre, A. Roos, Orientation dependence of plastic deformation in nickel-based single crystal superalloys: Discrete-continuous model simulations, Acta Mater. 58 (6) (2010) 1938–1951.

[9] A. M. Hussein, S. I. Rao, M. D. Uchic, T. A. Parthasarathay, J. A. El-Awady, The strength and dislocation microstructure evolution in superalloy microcrystals, J. Mech. Phys. Sol. 99 (2017) 146–162.

[10] A. Finel, D. Rodney, Phase field methods and dislocations, in: MRS Fall Meeting Y4.9, Boston, 2000, pp. 652–.

[11] Y. U. Wang, Y. M. Jin, A. M. Cuitiño, A. G. Khachaturyan, Nanoscale phase field microelasticity theory of dislocations: model and 3d simulations, Acta Mater. 49 (10) (2001) 1847–1857.

[12] M. Koslowski, A. M. Cuitiño, M. Ortiz, A phase-field theory of dislocation dynamics, strain hardening and hysteresis in ductile single crystals, J. Mech. Phys. Sol. 50 (12) (2002) 2597–2635.

[13] L.-Q. Chen, Phase-field models for microstructure evolution, Annu. Rev. Mater. Res. 32 (1) (2002) 113–140.

[14] D. Rodney, Y. Le Bouar, A. Finel, Phase field methods and dislocations, Acta Mater. 51 (1) (2003) 17–30.

[15] P.-L. Valdenaire, Y. Le Bouar, B. Appolaire, A. Finel, Variational scheme for full field micromechanics of inhomogeneous anisotropic media, J. Mech. Phys. Sol. (To be submitted).

[16] C. Denoual, Dynamic dislocation modeling by combining Peierls Nabarro and Galerkin methods, Phys. Rev. B 70 (2) (2004) 024106.

[17] C. Shen, Y. Wang, Incorporation of $\gamma$-surface to phase field model of dislocations: simulating dislocation dissociation in fcc crystals, Acta Mater. 52 (3) (2004) 683–691.

[18] N. Zhou, C. Shen, M. J. Mills, J. Li, Y. Wang, Modeling displacive – diffusional coupled dislocation shearing of precipitates in Ni-base superalloys, Acta Mater. 59 (9) (2011) 3484–3497.

[19] J. Mianroodi, A. Hunter, I. Beyerlein, B. Svendsen, Theoretical and computational comparison of models for dislocation dissociation and stacking fault/core formation in fcc crystals, J. Mech. Phys. Sol. 95 (2016) 719–741.

[20] L. D. Landau, E. M. Lifshitz, Theory of Elasticity, Pergamon Edition, Vol. 7 of Course of theoretical physics, 1959.

[21] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and nonlinear properties of composites, C. R. Acad. Sci. Paris, Ser II 318 (1994) 1417–1423.

[22] J. Zhu, L. Q. Chen, J. Shen, Morphological evolution during phase separation and coarsening with strong inhomogeneous elasticity, Modelling Simul. Mater. Sci. Eng. 9 (2001) 499–511.

[23] G. Boussinot, Y. Le Bouar, A. Finel, Phase-field simulations with inhomogeneous elasticity: Comparison with an atomic-scale method and application to superalloys, Acta Mater. 58 (2010) 4170–4181.

[24] J. Zeman, J. Vondřejc, J. Novák, I. Marek, Accelerating a fft-based solver for numerical homogenization of periodic media by conjugate gradients, J. Comp. Phys. 229 (21) (2010) 8065–8071.

[25] A. Finel, F. Ducastelle, On the phase diagram of the fcc Ising model with antiferromagnetic first-neighbour interactions, Europhys. Lett. 1 (3) (1986) 135 ; erratum 1 (10) (1986) 543.

[26] R. Tétot, A. Finel, F. Ducastelle, Superdegenerate point in fcc phase diagram: CVM and Monte Carlo investigations, J. Stat. Phys. 61 (1) (1990) 121–141.

[27] C. Shen, Y. Wang, Phase field model of dislocation networks, Acta Mater. 51 (9) (2003) 2595–2610.

[28] F. R. N. Nabarro, The synthesis of elastic dislocation fields, Phil. Mag. 42 (334) (1951) 1224–1231.

[29] T. Mura, Micromechanics of defects in solids, Vol. 3, Springer Netherlands, 1987.

[30] V. L. Ginzburg, L. D. Landau, On the theory of superconductivity, Sov. Phys. JETP. 20.

[31] J. R. Rice, Dislocation nucleation from a crack tip: An analysis based on the Peierls concept, J. Mech. Phys. Sol. 40 (2) (1992) 239–271.

[32] W. J. Cahn, M. S. Allen, A microscopic theory for domain wall motion and its experimental verification in fe-al alloy domain growth kinetics, J. Phys. Colloq. 38 (C7) (1977) C7–51.

[33] J. Frenkel, T. Kontorova, On the theory of plastic deformation and twinning, Izv. Akad. Nauk, Ser. Fiz. 1 (1939) 137–149.

[34] J. P. Hirth, J. Lothe, Theory of dislocations, John Wiley and Sons Edition, 1982.

[35] L. Kubin, Dislocations, Mesoscale Simulations and Plastic Flow, Oxford Series on Materials Modelling, OUP Oxford, 2013.
[36] P.-A. Geslin, B. Appolaire, A. Finel, A phase field model for dislocation climb, Appl. Phys. Lett. 104 (1) (2014) 011903–011903–4.

[37] P.-A. Geslin, B. Appolaire, A. Finel, Multiscale Theory of Dislocation Climb, Phys. Rev. Lett. 115 (26) (2015) 265501.

[38] A. Ruffini, A. Finel, Phase-field model coupling cracks and dislocations at finite strain, Acta Mater. 92 (0) (2015) 197–208.

[39] F. Pettinari-Sturmel, A. Coujou, N. Clément, The fluctuation of short-range order evidenced by mobile dislocations in the γ-phase of a nickel-based superalloy, Mat. Sci. Eng. A 400 (2005) 114–117.

[40] V. I. Levitas, D. L. Preston, D.-W. Lee, Three-dimensional landau theory for multivariant stress-induced martensitic phase transformations. III. alternative potentials, critical nuclei, kink solutions, and dislocation theory, Phys. Rev. B 68 (13) (2003) 134201–.

[41] V. I. Levitas, M. Javanbakht, Advanced phase-field approach to dislocation evolution, Phys. Rev. B 86 (14) (2012) 140101.

[42] V. I. Levitas, M. Javanbakht, Thermodynamically consistent phase field approach to dislocation evolution at small and large strains, J. Mech. Phys. Sol. 82 (2015) 345–366.

[43] S. Y. Hu, Y. L. Li, Y. X. Zheng, L. Q. Chen, Effect of solutes on dislocation motion – a phase-field simulation, Int. J. Plast. 20 (3) (2004) 403–425.

[44] M. J. Demkowicz, J. Wang, R. G. Hoagland, Chapter 83 Interfaces Between Dissimilar Crystalline Solids, in: A Tribute to F.R.N. Nabarro, Vol. Volume 14, Elsevier, 2008, pp. 141–205.

[45] S. H. Oh, M. Legros, D. Kiener, G. Dehm, In situ observation of dislocation nucleation and escape in a submicrometre aluminium single crystal, Nat. Mater. 8 (2) (2009) 95–100.

[46] A. Epishin, B. Bokstein, I. Svetlov, B. Fedelich, T. Feldmann, A. Ruffini, A. Finel, B. Viguier, D. Poquillon, Vacancy model for pore annihilation during hip of single-crystal nickel-base superalloys, Materialovedeniye (To be published).

[47] B. Hamilton, S. Bilbao, On finite difference schemes for the 3-d wave equation using non-cartesian grids, in: Proc. Sound and Music Computing (SMC) Conf., Stockholm, Sweden, 2013, pp. 592–599.