%~Mouliné par MaN_auto v.0.34.0 (b0f8651d) 2025-06-23 17:29:05
\documentclass[CRMECA,Unicode,biblatex,published]{cedram}

\addbibresource{CRMECA_Appolaire_20250129.bib}

\usepackage{siunitx}

\DeclareSIUnit\angstrom{\text{Å}}
\newcommand{\tensor}[1]{\boldsymbol{#1}}
\newcommand{\im}{\mathrm{i}}
\newcommand{\myvec}[1]{\vec{#1}}
\newcommand{\micromegas}{\emph{microMEGAS}}
\newcommand{\strain}{\boldsymbol{\varepsilon}}
\newcommand{\stress}{\boldsymbol{\sigma}}
\newcommand{\polarization}{\boldsymbol{\tau}}
\newcommand{\Green}{\boldsymbol{G}}

\newcommand*{\longrange}{\mathrm{long}}
\newcommand*{\shortrange}{\mathrm{short}}

%%%---------------------------------------------------------------------------

%%% BAR & OVERLINE

%%% Better widebar macro
\makeatletter
\let\save@mathaccent\mathaccent
\newcommand*\if@single[3]{%
  \setbox0\hbox{${\mathaccent"0362{#1}}^H$}%
  \setbox2\hbox{${\mathaccent"0362{\kern0pt#1}}^H$}%
  \ifdim\ht0=\ht2 #3\else #2\fi
}
%The bar will be moved to the right by a half of \macc@kerna, which is computed by amsmath:
\newcommand*\rel@kern[1]{\kern#1\dimexpr\macc@kerna}
%If there's a superscript following the bar, then no negative kern may follow the bar;
%an additional {} makes sure that the superscript is high enough in this case:
\newcommand*\widebar[1]{\@ifnextchar^{{\wide@bar{#1}{0}}}{\wide@bar{#1}{1}}}
%Use a separate algorithm for single symbols:
\newcommand*\wide@bar[2]{\if@single{#1}{\wide@bar@{#1}{#2}{1}}{\wide@bar@{#1}{#2}{2}}}
\newcommand*\wide@bar@[3]{%
  \begingroup
  \def\mathaccent##1##2{%
%Enable nesting of accents:
    \let\mathaccent\save@mathaccent
%If there's more than a single symbol, use the first character instead (see below):
    \if#32 \let\macc@nucleus\first@char \fi
%Determine the italic correction:
    \setbox\z@\hbox{$\macc@style{\macc@nucleus}_{}$}%
    \setbox\tw@\hbox{$\macc@style{\macc@nucleus}{}_{}$}%
    \dimen@\wd\tw@
    \advance\dimen@-\wd\z@
%Now \dimen@ is the italic correction of the symbol.
    \divide\dimen@ 3
    \@tempdima\wd\tw@
    \advance\@tempdima-\scriptspace
%Now \@tempdima is the width of the symbol.
    \divide\@tempdima 10
    \advance\dimen@-\@tempdima
%Now \dimen@ = (italic correction / 3) - (Breite / 10)
    \ifdim\dimen@>\z@ \dimen@0pt\fi
%The bar will be shortened in the case \dimen@<0 !
    \rel@kern{0.6}\kern-\dimen@
    \if#31
      \overline{\rel@kern{-0.6}\kern\dimen@\macc@nucleus\rel@kern{0.4}\kern\dimen@}%
      \advance\dimen@0.4\dimexpr\macc@kerna
%Place the combined final kern (-\dimen@) if it is >0 or if a superscript follows:
      \let\final@kern#2%
      \ifdim\dimen@<\z@ \let\final@kern1\fi
      \if\final@kern1 \kern-\dimen@\fi
    \else
      \overline{\rel@kern{-0.6}\kern\dimen@#1}%
    \fi
}%
  \macc@depth\@ne
  \let\math@bgroup\@empty \let\math@egroup\macc@set@skewchar
  \mathsurround\z@ \frozen@everymath{\mathgroup\macc@group\relax}%
  \macc@set@skewchar\relax
  \let\mathaccentV\macc@nested@a
%The following initialises \macc@kerna and calls \mathaccent:
  \if#31
    \macc@nested@a\relax111{#1}%
  \else
%If the argument consists of more than one symbol, and if the first token is
%a letter, use that letter for the computations:
    \def\gobble@till@marker##1\endmarker{}%
    \futurelet\first@char\gobble@till@marker#1\endmarker
    \ifcat\noexpand\first@char A\else
      \def\first@char{}%
    \fi
    \macc@nested@a\relax111{\first@char}%
  \fi
  \endgroup
}
\makeatother

\let\oldbar\bar
\renewcommand*{\bar}[1]{{\mathchoice{\widebar{#1}}{\widebar{#1}}{\widebar{#1}}{\oldbar{#1}}}}

%%%---------------------------------------------------------------------------

%%% DELIMITERS

\usepackage{mathtools}
\DeclarePairedDelimiter{\parens}{\lparen}{\rparen}
\DeclarePairedDelimiter{\braces}{\{}{\}}
\DeclarePairedDelimiter{\abs}{\lvert}{\rvert}
\DeclarePairedDelimiter{\norm}{\lVert}{\rVert}
\DeclarePairedDelimiter{\bracks}{[}{]}
%\DeclarePairedDelimiter{\bbracks}{\llbracket}{\rrbracket}
%\DeclarePairedDelimiter{\angles}{\langle}{\rangle}
%\DeclarePairedDelimiter{\aangles}{\langle\!\langle}{\rangle\!\rangle}
%\DeclarePairedDelimiter{\pparens}{(\!(}{)\!)}

%%%---------------------------------------------------------------------------

\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}}
\newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel}
\newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel}
\newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel}
\newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel}
\let\oldhat\hat
\renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}}
\let\oldforall\forall
\renewcommand*{\forall}{\mathrel{\oldforall}}

%%%---------------------------------------------------------------------------

%%% Changer la largeur des adresses

% Jouer avec \hsize :
\makeatletter
\def\@setafterauthor{%
  \vglue3mm%
%  \hspace*{0pt}%
\begingroup\hsize=11cm\advance\hsize\abstractmarginL\raggedright
\noindent
%\hspace*{\abstractmarginL}\begin{minipage}[t]{10cm}
   \leftskip\abstractmarginL
  \normalfont\Small
  \@afterauthor\par
\endgroup
\vskip2pt plus 3pt minus 1pt
}
\makeatother

%%%---------------------------------------------------------------------------

\title{Efficient FFT solver for discrete dislocation dynamics with sharp field description of plastic strain in FCC metals}
\alttitle{Un solveur FFT efficace pour la dynamique des dislocations discrètes dans les métaux CFC}

\author{\firstname{Luis} \lastname{Eon}}
\address{Universit{\'e} Paris-Saclay, ONERA, CNRS, Laboratoire d'{\'E}tude des Microstructures (LEM), 92322 Ch{\^a}tillon, France}
\email{luis.eon-ext@onera.fr}

\author{\firstname{Benoît} \lastname{Appolaire}}
\address{Université de Lorraine, CNRS, IJL, France}
\email{benoit.appolaire@univ-lorraine.fr}

\author{\firstname{Riccardo} \lastname{Gatti}}
\address[1]{Universit{\'e} Paris-Saclay, ONERA, CNRS, Laboratoire d'{\'E}tude des Microstructures (LEM), 92322 Ch{\^a}tillon, France}
\email{riccardo.gatti@cnrs.fr}

\keywords{\kwd{Discrete dislocation dynamics} \kwd{fast Fourier transform} \kwd{eigenstrain} \kwd{sharp}}
\altkeywords{\kwd{Dynamique des dislocations discrètes} \kwd{transformée de Fourier rapide (FFT)} \kwd{déformation propre (eigenstrain)} \kwd{champ abrupt}}

\begin{abstract}
Discrete Dislocation Dynamics (DDD) is an established mesoscopic numerical method for simulating dislocation motion to determine the plastic behaviour of metals. However, DDD's reliance on analytical expressions for internal stress fields limits its application to infinite isotropic media. The Discrete-Continuous Model~(DCM), which couples DDD with a finite element elastic solver, has been proposed to handle more complex boundary conditions by regularising plastic strain using the eigenstrain formalism. However, this model still relies on analytical solutions for short-range dislocation interactions, making it unsuitable for anisotropic media. In this work, we present an improved DCM framework that replaces the FE solver with a Fast Fourier transform~(FFT) solver for improved computational efficiency and full numerical calculation of stress fields, eliminating the need for analytical corrections. The proposed FFT solver employs a discrete theory of Green's operators and uses a sharp eigenstrain field to describe dislocations. The solver mesh is aligned with the face-centred cubic (fcc) lattice of the DDD, forming an octahedral cell to address symmetry artefacts around $\{111\}$ slip planes. Our FFT-based approach successfully maintains numerical stability by representing fcc dislocations as sharp fields without generating oscillations. This coupling allows the study of plasticity in anisotropic materials and interactions between dislocations and diffuse inclusions, such as precipitates, without short-range stress corrections.
\end{abstract}

\begin{altabstract}
La dynamique des dislocations discrètes (DDD) est une méthode numérique mésoscopique bien établie pour simuler le mouvement des dislocations et déterminer le comportement plastique des métaux. Toutefois, elle repose généralement sur les expressions analytiques des champs de contraintes obtenues dans des milieux isotropes infinis, ce qui en limite l'application. Le modèle discret-continu~(DCM), qui couple la DDD avec un solveur élastique par éléments finis (EF) grâce au formalisme des déformations libres, avait été proposé pour traiter des conditions aux limites plus complexes. Cependant, ce modèle utilise également des solutions analytiques pour les interactions à courte portée, ce qui le rend inadapté aux milieux anisotropes. Dans ce travail, nous présentons une version améliorée du DCM en remplaçant le solveur EF par un solveur utilisant les transformées de Fourier rapides~(FFT), plus efficace et supprimant la nécessité des corrections analytiques à courte portée. Pour cela, ce solveur s'appuie sur trois ingrédients~: (i) un champ de déformation propre abrupt (discontinu à l'échelle du maillage) pour décrire les dislocations~; (ii) un maillage CFC conforme à celui du code DDD sur lequel un ``stencil'' octaédrique est utilisé pour préserver les symétries des plans de glissement $\{111\}$~; (iii) enfin un opérateur de Green discret. Ce solveur parvient ainsi à gérer des dislocations décrites par des champs abrupts sans générer d'oscillations. Couplé à la DDD, il permet d’étudier la plasticité dans les matériaux anisotropes et les interactions entre dislocations et des inclusions diffuses, sans corriger les contraintes aux cœurs des dislocations comme ce qui est fait habituellement.
\end{altabstract}

\dateposted{2025-09-11}
\begin{document}

\makeatletter
\def\cdr@preabstracthook{\par\smallskip}
\makeatother

\maketitle

\section{Introduction}

Discrete Dislocation Dynamics (DDD) is an efficient numerical method, at the mesoscopic scale, in which dislocation motion is determined via the Peach--Koehler equation, which relates the applied and internal stress fields in the material to the force acting on each dislocation. The applied stress field is generated by an external load, while the internal stress field, which depends on the dislocation-dislocation interaction, is usually based on analytical expressions assuming an infinite isotropic medium~\cite{1967_Hirth, 1995_Devincre, 2006_Cai}. Therefore, although DDD is a reliable tool for predicting plastic strain, anisotropic media, free surfaces or heterogeneous media, including cracks, cannot be easily handled. To overcome these limitations and to consider more complex boundary conditions, the Discrete-Continuous Model~(DCM)~\cite{2014_Vattre, 2016_Jamond} has been proposed by coupling DDD with an elastic finite element solver~(FE).

The DCM is based on the regularization of the plastic strain: dislocation loops are described by the eigenstrain formalism~\cite{1987_Mura}. In the initial formulation of the model, the stress field is divided into two parts: (i) short-range dislocation-dislocation interactions, which are treated analytically with the same formulas as in classical DDD at a distance smaller than the regularization length; (ii) the long-range stress field, including all other interactions due to external loading and free boundaries, which is calculated numerically by FE. To ensure continuity between the short-range and long-range contributions, the plastic slip regularization is computed using the Burgers vector spreading function proposed by Cai et al.~\cite{2006_Cai} to derive non-singular continuum formulas for the dislocation stress field. This decomposition is straightforward for dealing with free surfaces and heterogeneous media, but still depends on analytical expressions to calculate the short-range interaction and is therefore not viable for modeling plasticity in anisotropic media. Consequently, in this work we propose an approach based on the DCM framework to compute the full stress field numerically without the need for analytical corrections.

In addition to physical fidelity, our approach aims to improve numerical performance. For large cumulative plastic strains, finite element~(FE) solvers can be computationally expensive and memory-intensive due to matrix assembly and preconditioning. Fast Fourier transform~(FFT)-based mechanical solvers have gained popularity in recent years due to their computational efficiency and scalability in periodic media~\cite{1998_Moulinec, 2014_Willot, 2015_Willot, 2016_Schneider, 2016_Valdenaire,2021_Schneider}. However, FFT solvers that use continuous differential operators have been shown to suffer from Gibbs oscillations when elastic fields exhibit sharp variations, as is the case near dislocation cores.

This issue is not specific to material heterogeneities; it also appears in early coupling attempts between discrete dislocation dynamics~(DDD) and FFT solvers~\cite{2015_Bertin, 2018_Bertin}, where dislocation cores had to be spread to suppress the spurious oscillations. Some strategies have been proposed to address this issue by modifying the discretization of the equilibrium equation via carefully chosen stencils (see, for example, \cite{2014_Willot,2015_Willot,2016_Valdenaire,2016_Schneider,2017_Ruffini,2017_Schneider,2018_Eloh,2019_Eloh}, or~\cite{2021_Schneider} for a review). Similar ideas have been used in the context of Field Dislocation Mechanics~\cite{BERBENNI20144157,DJAKA2017921}, where centered finite differences schemes expressed on Cartesian grids have been introduced to reduce the oscillations. While these methods can improve numerical robustness, they do not address the fundamental geometric incompatibility between cubic meshes and the $\{111\}$ slip planes that govern dislocation motion in FCC crystals.

In the present study, we develop and apply a DDD-FFT coupling strategy designed specifically to model plasticity in FCC metals. Our approach centers on an FFT solver based on an ad hoc stencil that respects the symmetries of $\langle 110 \rangle \{111\}$ dislocations. Unlike previous methods, our solver can handle sharp eigenstrain fields without regularization, enabling accurate computations of the stress field. This capability relies on several key innovations. First, the FFT grid aligns with the FCC lattice used in the DDD model. Integration cells conform to the geometry of the $\{111\}$ slip planes. This alignment improves symmetry compatibility and avoids artifacts typical of cubic discretizations. Second, we designed a stencil tailored to this grid, enabling the solver to accurately process sharp fields. Third, our method does not require any analytical short-range correction, enabling a fully numerical resolution of dislocation-induced internal stresses. To our knowledge, no existing DDD–FFT coupling combines these three features. Hence, this new framework overcomes the limitations of previous hybrid approaches, such as the discrete-continuous model~(DCM), and opens new possibilities for simulating anisotropic materials and their interactions with microstructural features, such as inclusions or precipitates.

This paper is organized as follows. First, we explain how dislocations are described using sharp eigenstrain fields on an fcc lattice. Second, we introduce the FFT elastic solver based on a suitable template that is able to handle sharp eigenstrain fields on the fcc grid. Next, the coupling between the DDD and the FFT solvers is presented. Finally, the last section is devoted to several examples of increasing complexity that illustrate the advantages of our new coupling, from static loops to complex dislocation microstructures around an inclusion.

\section{Description of dislocations with eigenstrain}\label{section_eigenstrain}

\subsection{Mechanical equilibrium in the small strain framework}

We have considered the small-strain framework of linear elasticity, where stress $\stress$ and strain $\strain$ are symmetric tensors. The strain field is derived from the displacement $\myvec{u}$:
\begin{equation}\label{compatibility}
\varepsilon_{ij}(\myvec{x}) = \frac{1}{2} \parens[\Bigg]{\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}}
\end{equation}
in a standard Cartesian frame. The total strain can be decomposed as the sum of an elastic strain~$\tensor{\varepsilon}^e$ and an inelastic strain $\tensor{\varepsilon}^\star$ (eigenstrain) :
\begin{equation}\label{decomp_strain}
\varepsilon_{ij}(\myvec{x}) = \varepsilon^e_{ij}(\myvec{x}) + \varepsilon^\star_{ij}(\myvec{x})
\end{equation}
Hooke's law relates stress and strain:
\begin{equation}\label{hooke}
\sigma_{ij}(\myvec{x}) = C_{\mathit{ijkl}}(\myvec{x}) \, \varepsilon^e_{kl}(\myvec{x}) = C_{\mathit{ijkl}}(\myvec{x}) \parens[\big]{\varepsilon_{kl}(\myvec{x}) - \varepsilon_{kl}^{\star}(\myvec{x})}
\end{equation}
where $C_{\mathit{ijkl}}$ are the components of the stiffness tensor, and where we use the Einstein notation for the sums over repeated indices.

The mechanical equilibrium writes:
\begin{equation}\label{equilibrium}
\frac{\partial\sigma_{ij}}{\partial x_j} = 0 \quad \forall i=1,2,3.
\end{equation}

\subsection{Discrete-Continuous Model and dislocation core regularization} \label{DCM}

Discrete Dislocation Dynamics is a method for simulating the motion and interactions of dislocations in a continuous elastic medium to study the plasticity of crystals at the mesoscale. In most DDD codes, the dislocation lines are divided into a discrete set of interconnected dislocation segments. As mentioned in the introduction, the stress field generated by each segment is based on the well-known analytical expressions assuming an isotropic, infinite, and homogeneous medium~\cite{1967_Hirth}. These assumptions are cumbersome and are not adapted when we consider an anisotropic material, for example, a finite medium with free surfaces or the presence of any inclusions (interfaces) that modify the boundary conditions. In fact, the interactions between dislocations and free surfaces or interfaces generate a stress field that affects the movement of the dislocations, which is difficult, and in some cases impossible, to calculate analytically. The Discrete-Continuous Model~(DCM) was introduced to address this problem. It was first introduced by Lemarchand and was used by Vattré~\cite{2014_Vattre} to observe the interaction between dislocations and precipitates in $\gamma/\gamma'$ superalloys. Then it was improved by Jamond et al.~\cite{2016_Jamond} and used to study the interaction between dislocations and static cracks by Korzeczek~\cite{2017_Korzeczek}. This model consists of combining the use of analytical formulas and a numerical stress field calculated by the FE method to determine the internal stress field of a dislocation. In continuous theory, dislocations can be considered as homogeneous inclusions that induce a shear (plastic) eigenstrain. The value of this shear eigenstrain can be determined at any point in space by calculating the area swept by the dislocation segments in each of the slip planes.

For a slip plane $g$, the area swept in a time interval $\delta t$ is denoted by $\delta S^g$. The slip increment is:
\begin{equation}
\delta \gamma^p = \frac{b \, \delta S^g}{V}
\end{equation}
where $b$ is the norm of the Burgers vector and $V$ is a reference volume that contains the swept area $\delta S^g$. Usually, $V$ is chosen as the volume of the simulation domain. The eigenstrain field $\delta \varepsilon_{ij}^{\star}$ generated by the dislocations during their movement from $t$ to $t+\delta t$ is then the sum over all slip systems:
\begin{equation}
\delta \varepsilon^{\star}_{ij} = \sum^{12}_{g=1} \frac{1}{2} \parens{m_i^g n_j^g + m_j^g n_i^g} \delta \gamma^p
\end{equation}
with $\myvec{b}^g = b\,\myvec{m}^g$, $\myvec{m}^g$ the unit Burgers vector and $\myvec{n}^g$ the unit normal to the $g$ slip plane.

The code \micromegas{}~\cite{2011_Devincre} used for this work, is a lattice-based DDD code. Dislocations are described by straight connected segments that lie on a homothety of the atomic lattice. In the initial version of DCM, the FE mesh and the DDD lattice do not coincide, so the plastic strain field is regularized at the different nodes of the FE mesh by distributing the Burgers vector, as presented by Cai et al.~\cite{2006_Cai} in their nonsingular continuum theory of dislocations. This is possible by convolution of the plastic eigenstrain with a regular function to obtain a nonsingular field.

Thus, the calculated stress field of a dislocation is affected by regularization. In particular, it underestimates the short-range interaction between dislocations. To solve this problem, the stress field generated by a dislocation is split into two contributions: the long-range stress field $\stress^{\longrange}$, which takes into account the different obstacles and boundary conditions, and the short-range stress field $\stress^{\shortrange}$, which aims to approximate as closely as possible the stress divergence at the core of a dislocation. The $\stress^{\longrange}$ field is the result of the FE calculation~\cite{2016_Jamond}, and the $\stress^{\shortrange}$ field comes from analytical formulas that assume an infinite isotropic homogeneous medium.

\subsection{Sharp field description of dislocations}\label{subsection_sharp_description}

The DCM limitation presented in the previous section can be overcome by considering a sharp description of the eigenstrain field, i.e.\ avoiding regularization. In \micromegas{}, the extremities and centers of the dislocation segments lie in lattice positions similar to fcc for fcc materials. Each node of the DDD simulation lattice can be associated with an elementary volume that forms a Voronoi cell. In an fcc lattice, this volume is a rhombic dodecahedron of volume~$V_d = \frac{d^3}{4}$, where $d$ is the distance from one side of the cubic cell to the other (along $\langle 100 \rangle$ directions, see Figure~\ref{fig1}a). The intersection of this polyhedron with a plane of the family of glide planes $\{111\}$ passing through a simulation node (that is, the center of the dodecahedral cell) results in a regular hexagon of surface $S_h= \frac{\sqrt{3}}{4} d^2$ (Figure~\ref{fig1}b).

Each of the nodes $\myvec{x}$ of the grid has an assigned value of the eigenstrain, depending on the swept area within the voxel, which is written in general terms:
\begin{equation}\label{discrete_eigenstrain}
\varepsilon^{\star}_{ij}(\myvec{x}) = \frac{1}{2} \frac{b \psi S_h}{V_d} \parens{m_i n_j + m_j n_i} = \frac{1}{2} \frac{\sqrt{3}\, \psi b}{d}\parens{m_i n_j + m_j n_i}
\end{equation}
where we see the expression of the elementary shear associated with a dislocation of the Burgers vector $\myvec{b}=b\,\myvec{m}$ in the plane of the normal $\myvec{n}$. Note that this representation implies that the dislocation loop is not a pure 2D object, but has a thickness equal to the height of a dodecahedral voxel in the $\myvec{n}$ direction~($\langle 111 \rangle$ directions).

\begin{figure}[h!]
\includegraphics[width=0.7\textwidth]{Figure_1.pdf} \caption{(a) Rhombic dodecahedral elementary volume of an fcc lattice. (b) Hexagonal section through the $(\overline{1}11)$ plane.} \label{fig1}
\end{figure}

In~\eqref{discrete_eigenstrain}, $\psi$ is the fraction of the hexagonal area $S_h$ that is swept into a voxel of volume~$V_d$ to account for all possible situations:
\begin{itemize}
\item no dislocation swept the voxel, $\psi=0$, the eigenstrain value is zero (Figure~\ref{fig2}a);
\item the node is inside a dislocation loop, i.e.\ a segment has completely swept the hexagonal surface $S_h$ and $\psi=1$ (Figure~\ref{fig2}b);
\item finally, there are many intermediate situations in which the hexagon is partially swept so that $\psi \in\mathopen{]}0,1\mathclose{[}$ (Figure~\ref{fig2}c). Since \micromegas{} has only eight possible segment directions, listing the limited number of possible $\psi$ fractions can be done.
\end{itemize}

\begin{figure}[h!]
\includegraphics[width=0.9\textwidth]{Figure_2.pdf} \caption{\small{Three possible situations of swept areas in a voxel of the fcc mesh.}} \label{fig2}
\end{figure}

\subsection{Eigenstrain field dynamics} \label{subsection_eigenstrain_dynamics}

The eigenstrain field is updated at each time step $t$ of \micromegas{}. The new eigenstrain field at $t+\delta t$ is built by addition to the previous one at time $t$ (Figure~\ref{fig3}). The sign of the increment $\delta\varepsilon^{\star}_{ij}$ in a given voxel depends on the direction of the segment slip, which is determined by the sign of its displacement. In tension, the increment is positive if the loop increases its area, while it is negative if the loop decreases its area.

\begin{figure}[h!]
\includegraphics[width=0.5\textwidth]{Figure_3.pdf} \caption{\small{Eigenstrain increment as a function of swept surfaces in fcc mesh voxels.}}\label{fig3}
\end{figure}

\section{Octahedral solver}\label{section_octahedral_solver}

\subsection{General framework}

To solve the mechanical equilibrium, solvers based on Fast Fourier Transforms~(FFT) are becoming more and more popular thanks to their computation time performance. A coupling between DDD and FFT solver, already realized in the literature by Bertin et al.~\cite{2015_Bertin, 2018_Bertin}, seems indeed to improve the performance of the DCM. In fact, the kinematics and mechanical equilibrium equations~\eqref{compatibility}--\eqref{equilibrium} become simple algebraic expressions in Fourier space:
\begin{gather}\label{discrete_kinematics}
\widehat{\varepsilon}_{ij}(\myvec{q}) = \frac{1}{2}\bracks[\big]{D_j^u(\myvec{q})\,\hat{u}_i (\myvec{q}) + D_i^u(\myvec{q})\hat{u}_j (\myvec{q})}
\\\label{discrete_equilibrium}
D_j^\sigma(\myvec{q})\,\widehat{\sigma}_{ij}(\myvec{q}) = 0
\end{gather}
where $\strain$ is the strain fluctuations with respect to some average $\langle\strain\rangle$ and $u$ is the associated displacement fluctuations. $\widehat{}\,$ denotes the Fourier transform, $\myvec{q}$ is the wave vector and $\myvec{D}(\myvec{q})$ are the differential operators in Fourier space. Since these operators can be different depending on the equation under consideration, they are denoted with a superscript $u$ or $\sigma$.

Considering the case of homogeneous stiffness $C_{\mathit{ijkl}}(\myvec{x}) = C^0_{\mathit{ijkl}}$, and using Hooke's law~\eqref{hooke} to combine~\eqref{discrete_kinematics} and~\eqref{discrete_equilibrium}, one obtains:
\begin{equation}
\widehat{u}_i(\myvec{q}) = -\widehat{G}_{ij}(\myvec{q})\,D_k^\sigma(\myvec{q})\widehat{\tau}_{jk}(\myvec{q}) \label{calcu}
\end{equation}
where $\widehat{\Green}(\myvec{q})$ is the Fourier transform of the Green function for displacements and $\widehat{\polarization}(\myvec{q})$ the Fourier transform of the polarization tensor expressed as:
\begin{equation}\label{polarization}
\widehat{\tau}_{ij} = -C_{\mathit{ijkl}}^0 \widehat{\varepsilon}^{\,\star}_{kl}(\myvec{q})
\end{equation}
The Green function $\widehat{\Green}(\myvec{q})$ is defined as :
\begin{equation}\label{Green_u}
\widehat{G}_{ik}^{-1}(\myvec{q}) = C^0_{\mathit{ijkl}} D_j^{\sigma}(\myvec{q}) D_l^{u}(\myvec{q}) \quad\forall \myvec{q}\ne 0.
\end{equation}
In the case of inhomogeneous elastic constants, the polarization tensor depends on the contrast of elastic constants with respect to some reference stiffness $\Delta C_{\mathit{ijkl}}(\myvec{x}) = C_{\mathit{ijkl}}(\myvec{x}) - C^0_{\mathit{ijkl}}$:
\begin{equation}\label{polarization_tensor_inhomo}
\tau_{ij}(\myvec{x}) = \sigma_{ij}(\myvec{x}) - C_{\mathit{ijkl}}^0 \, \varepsilon_{kl}(\myvec{x}) = \Delta C_{\mathit{ijkl}}(\myvec{x}) \, \varepsilon_{kl}(\myvec{x}) - C_{\mathit{ijkl}}(\myvec{x}) \, \varepsilon^{\star}_{kl}(\myvec{x})
\end{equation}
The problem can be solved using a fixed point algorithm, as proposed in~\cite{1983_Khachaturyan,1998_Moulinec}: $\widehat{\myvec{u}}$ is computed iteratively with~\eqref{calcu} where $\polarization$ uses $\strain$ at a previous iteration.

In the case of linear elasticity with isotropic phases, the reference stiffness tensor ${C}^0$  is chosen as an isotropic tensor with a bulk modulus $K^0$ and a shear modulus $\mu^0$ defined as the average of the extreme values of the heterogeneous medium. This choice has been shown to improve convergence by reducing the spectral radius of the iterative operator, as demonstrated in~\cite{MichelMoulinec01,MoulinecSilva14}. This choice is particularly effective in media with high contrast or strongly differing phase properties. The convergence criterion employed follows the one proposed in the basic scheme of Moulinec and Suquet~\cite{1998_Moulinec}, and is based on a norm of the mechanical equilibrium:
\[
\epsilon = \norm{\operatorname{Div}\sigma/\sigma_{\mathrm{ref}}},
\]
where $\sigma_{\mathrm{ref}}$ is some reference stress. While Moulinec and Suquet have used the applied stress as the reference stress and the $L^2$ norm, we used $\sigma_{\mathrm{ref}}=\mu_0$ and the $L^\infty$  norm. With these choices, we have found that convergence is achieved when $\epsilon < 10^{-4}$. Tests with more drastic values such as $10^{-6}$ did not show significant differences at a much higher computational cost.

In their seminal work, Moulinec and Suquet~\cite{1998_Moulinec} considered the exact partial differential equations such that $\myvec{D}^u(\myvec{q}) = \myvec{D}^\sigma(\myvec{q}) = \im \myvec{q}$. However, oscillations due to the Gibbs phenomenon are observed when considering sharp variations in the eigenstrain field. These oscillations are detrimental to the convergence of the algorithms used to find the equilibrium. To get rid of this problem, different solutions have been proposed so far, based on discretized versions of~\eqref{compatibility} and~\eqref{equilibrium}, for example using finite difference (FD) approximations of the differential operators (e.g.~\cite{1998_Muller,2008_Willot,2016_Schneider,2016_Valdenaire}). However, as discussed in~\cite{2021_Schneider}, the choice of the FD stencil is critical. Indeed, depending on this stencil, oscillations can be generated whose nature is different from the Gibbs phenomenon. They can appear when the differential operators cancel at some $\myvec{q}_{c}\ne 0$ in the first Brillouin zone of the reciprocal grid ($\myvec{D}(\myvec{q_c})=0$). Although some stencils based on staggered grids have been proposed that are free of such problems~\cite{2016_Schneider, 2016_Valdenaire, 2021_Schneider}, they are based on cubic grids and may have other artefacts due to geometry when considering dislocation loops in fcc, as will be discussed later.

\subsection{Grids}

Therefore, we have adopted an ad hoc FD stencil based on two staggered fcc grids, corresponding to the grid used in DDD (Section~\ref{DCM}) and having the symmetries necessary for a proper description of $\{111\}$ dislocation loops by sharp eigenstrain platelets. The grids are defined by their rhombohedral primitive basis:
\begin{equation}\label{fcc_basis}
\begin{cases}
\myvec{a}_1 = \tfrac{1}{2} (\myvec{e}_1 + \myvec{e}_3) = \parens[\big]{\tfrac{1}{2},0,\tfrac{1}{2}}_{\mathcal{E}} 	\\
\myvec{a}_2 = \tfrac{1}{2} (\myvec{e}_1 + \myvec{e}_2) = \parens[\big]{\tfrac{1}{2},\tfrac{1}{2},0}_{\mathcal{E}} 	\\
\myvec{a}_3 = \tfrac{1}{2} (\myvec{e}_2 + \myvec{e}_3) = \parens[\big]{0,\tfrac{1}{2},\tfrac{1}{2}}_{\mathcal{E}}
\end{cases}
\end{equation}
expressed in the canonical Cartesian basis $\mathcal{E} = \{\myvec{e}_1, \myvec{e}_2, \myvec{e}_3 \}$ where $\abs{\myvec{e}_i}=d$ (Figure~\ref{fig4}a), with respect to which all tensors are defined (i.e., the components of the tensors are expressed in $\mathcal{E}$). The first grid (noted P) composed of nodes at $\myvec{r}^{\,p} {= m\,\myvec{a}_1 + n\,\myvec{a}_2 + p\,\myvec{a}_3}$ with $(m,n,p)\in\mathbb{Z}^3$ carries stress, strain, eigenstrain and stiffness fields (red in Figure~\ref{fig4}a) and coincides with the DDD grid. The second grid (noted S), shifted by $\myvec{R} = \tfrac{1}{2}(1,1,1)_\mathcal{E} = \tfrac{1}{2}(1,1,1)_\mathcal{B}$ with respect to the first one, with nodes located at $\myvec{r}^{\,s} = \myvec{r}^{\,p} + \myvec{R}$, carries the displacement field (blue in Figure~\ref{fig4}b).

\begin{figure}[ht]
\includegraphics[width=\textwidth]{Figure_4.pdf} \caption{(a) Left: Cartesian basis $\mathcal{E}$ (black vectors) and rhombohedral primitive basis of fcc $\mathcal{B}$ (red vectors). Right: rhombohedral system built with dodecahedral voxels. (b) The two staggered fcc grids: grid P in red carries strain, eigenstrain, stress and stiffness fields, and grid S in blue carries the displacement field.} \label{fig4}
\end{figure}

Since the fields are defined on the rhombohedral grids, their Fourier transforms are defined on the reciprocal grid built using the dual basis $\mathcal{B}^{\star}=(\myvec{a}^{\star}_1, \myvec{a}^{\star}_2,\myvec{a}^{\star}_3)$, where $\myvec{a}_i\cdot\myvec{a}^{\star}_j=\delta_{ij}$. This gives the discrete Fourier transforms of any components $t$ of strain, eigenstrain, stress and stiffness:
\begin{equation}
\widehat{t}(\myvec{q}) = \frac{1}{N} \sum_{\myvec{r}^{\,p}} t(\myvec{r}^{\,p}) \; \exp\parens{-\im\,\myvec{q}\cdot\myvec{r}^{\,p}}, \label{fourier_eps}
\end{equation}
when the discrete Fourier transforms of the displacement components reads:
\begin{equation}
\widehat{u}_i(\myvec{q}) = \frac{1}{N} \sum_{\myvec{r}^{\,p}} u_i(\myvec{r}^{\,p}+\myvec{R}) \; \exp\parens{-\im\,\myvec{q}\cdot\myvec{r}^{\,p}}, \label{fourier_disp}
\end{equation}
$N$ is the total number of nodes. The wave vectors write $\myvec{q}=q_w\,\myvec{a}_w^{\star}$ ($w\in\{1,2,3\}$) where $q_w=2\pi\alpha_w/N_w$ with $\alpha_w\in\bracks[\big]{-\frac{N_w}{2}+1,\frac{N_w}{2}}$ if $N_w$ is even or $\alpha_w\in\bracks[\big]{-\frac{N_w-1}{2},\frac{N_w-1}{2}}$ if $N_w$ is odd.

\subsection{Discretized kinematics}

Once the grids are defined, the relationship between displacement and strain is obtained by discretizing~\eqref{compatibility}. To ensure that the fields at the dislocation core have the correct symmetries, displacements must be at nodes that form a polyhedron (i) with the same number of vertices on each side of $\{111\}$ planes, (ii) that are equidistant from these planes. According to~\cite{2017_Ruffini}, the octahedron is the simplest choice that satisfies these conditions. The octahedron vertices surrounding a given node at $\myvec{r}^{\,p}$ belongs to the set $O(\myvec{r}^{\,p})$ and are indexed by $o\in[1,6]$ (Figure~\ref{fig5}).

The partial derivatives in~\eqref{compatibility} are approximated by finite differences between the displacement components at neighboring nodes of grid S:
\begin{equation}\label{eij_mnp_oct}
\varepsilon_{ij}(\myvec{r}^{\,p}) = \frac{1}{2d}\sum_{o\,\in\,O(\myvec{r}^{\,p})} \parens{\alpha_{j}^{o}\,u_i^{o} +\alpha_{i}^{o}\,u_j^{o}},
\end{equation}
where the $u_i^{o}$ are the Cartesian components, along the $i$ direction, of the neighbouring nodes $o\in O(\myvec{r^{\,p}})$ (Figure~\ref{fig5}) and where $\alpha_{i}^{o} = -1$, $0$ or $+1$ depending on the component $i$ and the position~$o$ in the octahedron. The choice of the neighboring displacement nodes must comply with the direction of the derivative. For example:
\begin{equation}\label{e12_oct}
\varepsilon_{12}(\myvec{r}^{\,p}) = \frac{1}{2d}\bracks[\Big]{u_1\parens[\big]{\myvec{r}^{\,p}+\tfrac{1}{2}\myvec{e}_2} - u_1\parens[\big]{\myvec{r}^{\,p}-\tfrac{1}{2}\myvec{e}_2} + u_2\parens[\big]{\myvec{r}^{\,p}+\tfrac{1}{2}\myvec{e}_1} - u_2\parens[\big]{\myvec{r}^{\,p}-\tfrac{1}{2}\myvec{e}_1}}
\end{equation}

\begin{figure}[!h]
\centering
\includegraphics[width=0.6\textwidth]{Figure_5.pdf}
\caption{Nodes $o$ of grid S (blue) where displacements are defined forming the octahedron $O(\myvec{r}^{\,p})$ surrounding the node at $\myvec{r}^p$ of grid P (red) where mechanical fields are defined.}\label{fig5}
\end{figure}

\subsection{Discretized mechanical equilibrium}

We consider that the elastic energy is piecewise constant in each dodecahedral voxel (i.e.\ the Wigner--Seitz cell of the fcc) of volume $d^3/4$ (Section~\ref{subsection_sharp_description}) surrounding the nodes of grid P. Hence, assuming that the average strain $\langle\strain\rangle=0$ for simplicity:
\begin{equation}\label{discrete_elastic_energy}
E_{el} = \frac{d^3}{4} \sum_{\myvec{r}^{\,p}} f_{el}(\myvec{r}^{\,p}) = \frac{d^3}{4} \sum_{\myvec{r}^{\,p}} \frac{1}{2} C_{\mathit{ijkl}}(\myvec{r}^{\,p}) \bracks[\big]{\varepsilon_{ij}(\myvec{r}^{\,p}) - \varepsilon^{\star}_{ij}(\myvec{r}^{\,p})} \bracks[\big]{\varepsilon_{kl}(\myvec{r}^{\,p}) - \varepsilon^{\star}_{kl}(\myvec{r}^{\,p})}
\end{equation}
The discretized mechanical equilibrium is obtained in a manner consistent with the discretisazion of the fields by minimizing~\eqref{discrete_elastic_energy} with respect to the discrete displacement of each node of grid S:
\begin{equation}\label{minim_elas}
\frac{\delta E_{el}}{\delta u_j(\myvec{r}^{\,s})} = \sum_{o\,\in\,O(\myvec{r}^{\,s})} \frac{\delta E_{el}}{\delta\varepsilon_{ij}^o} \frac{\delta\varepsilon_{ij}^o}{\delta u_k(\myvec{r}^{\,s})} = 0
\end{equation}
where $o$ denotes the index of a node of grid P belonging to the octahedron $O(\myvec{r}^{\,s})$ surrounding the node at $\myvec{r}^{\,s}$ in grid S. We obtain the following general relation:
\begin{equation}\label{eq_meca_sys_lin_oct}
\frac{1}{d} \sum_{o\,\in\,O(\myvec{r}^{\,s})} \beta_{j}^{o}\,\sigma_{ij}^{o} = 0
\end{equation}
where $\beta_j^o = -1, 0$ or $+1$ depending on the position of the node.

\subsection{Operators in Fourier space} \label{fourier_operators}

When expressed in Fourier space, \eqref{eij_mnp_oct} and~\eqref{eq_meca_sys_lin_oct} write:
\begin{gather}
\widehat{\varepsilon}_{ij}(\myvec{q}) = - \frac{1}{2} \bracks[\big]{D^{\star}_j(\myvec{q}) \; \widehat{u}_i(\myvec{q}) + D^{\star}_i(\myvec{q}) \; \widehat{u}_j(\myvec{q})} \label{discrete_grad_oct}
\\
D_j(\myvec{q}) \,\widehat{\sigma}_{ij}(\myvec{q}) = 0 \label{discrete_equilibrium_oct}
\end{gather}
where $D_i(\myvec{q})$ and its complex conjugate $D_i^{\star}(\myvec{q})$ are the components in canonical Cartesian direction~$i$ of the discrete differential operators. We recover the generic expressions~\eqref{discrete_kinematics} and~\eqref{discrete_equilibrium} with $\myvec{D}^\sigma = \myvec{D}$ and $\myvec{D}^u = \myvec{D}^\star$. Written explicitly in terms of the components $q_w$ of the wave vectors in basis $\mathcal{B}$, the components of the operators read:
\begin{equation}\label{operators_oct}
\begin{aligned}
D_1(\myvec{q}) &= \frac{1}{d} \bracks[\big]{e^{\im (q_1+q_2)} - e^{\im q_3}}, \\
D_2(\myvec{q}) &= \frac{1}{d} \bracks[\big]{e^{\im (q_2+q_3)} - e^{\im q_1}}, \\
D_3(\myvec{q}) &= \frac{1}{d} \bracks[\big]{e^{\im (q_3+q_1)} - e^{\im q_2}}.
\end{aligned}
\end{equation}

In order to analyze the solver stability, one must examine if there exists any wave vector $\myvec{q}_c\ne 0$ in the first Brillouin zone of fcc that cancels $\myvec{D}$. In that case, \eqref{discrete_equilibrium_oct} becomes undetermined and any fluctuations with wavelength $2\pi/\abs{\myvec{q}_c}$ can appear. In general, such $\myvec{q}_c$ are located at the boundary of the Brillouin zone so that checkerboard patterns plague the calculations near sharp variations of the fields. This kind of analysis is simplified when the wave vectors are expressed in the Cartesian dual basis $\mathcal{E}^{\star} = \{\myvec{e}^{\,\star}_1,\myvec{e}^{\,\star}_2, \myvec{e}^{\,\star}_3\}$ with $\myvec{e}_i^{\,\star}\cdot\myvec{e}_j = \delta_{ij}$: $\myvec{q} = q_1'\myvec{e}_1^{\,\star} + q_2'\myvec{e}_2^{\,\star} + q_3'\myvec{e}_3^{\,\star}$. In this basis, the operators read:
\begin{equation}
\begin{aligned}
D_1(q_1', q_2', q_3') = \frac{1}{d}\,e^{\im\frac{1}{2}(q_2' + q_3')} \bracks[\big]{e^{\im q_1'} - 1}, \\
D_2(q_1', q_2', q_3') = \frac{1}{d}\,e^{\im\frac{1}{2}(q_1' + q_3')} \bracks[\big]{e^{\im q_2'} - 1}, \\
D_3(q_1', q_2', q_3') = \frac{1}{d}\,e^{\im\frac{1}{2}(q_1' + q_2')} \bracks[\big]{e^{\im q_3'} - 1}.
\end{aligned}
\end{equation}
The wave vectors in the family $2\pi\langle 1,0,0\rangle_\mathcal{E}$ cancel the three components of the differential operator. Consequently, the octahedral solver is not unconditionally stable. However, in the case of dislocations in $\{111\}$ planes, the solver does not exhibit any oscillations, as will be shown in Section~\ref{results}.


\section{Coupling DDD with the octahedral solver}\label{section_coupling}

The coupling between the FFT solver and \micromegas{} follows the same principle as the DCM~\cite{2016_Jamond}, except that transferring the stress field from the solver to DDD is simpler. Indeed, due to the performance of FFT solvers with respect to FE solvers, the FFT grid nodes and the \micromegas{} simulation lattice, as described above, can coincide. In this way, it is possible to return the stress at any point in the \micromegas{} grid, without needs of any interpolation. Moreover, the performance of the octahedral FFT solver has been increased thanks to a parallelization based on P3DFFT libraries~\cite{2012_Pekurovsky}.

At each time step, the coupling is carried out in four consecutive steps:
\begin{enumerate}\alphenumi
\item Between $t$ and $t+\delta t$, \micromegas{} applies the displacement, computed via the Peach--Koehler force based on the stress field at the time $t$, to dislocations. In particular, it manages the discretization of the loops, connects each discrete segment, and handles junction formations, annihilations and cross-slip events. The areas swept during $\delta t$ and the new positions of the segment midpoints are sent to the FFT solver by MPI communication. The mechanical loading is incremented by \micromegas{} with adaptive algorithm taking into account the applied stress or strain according to the cumulative plastic strain in the simulation volume. This loading is also sent to the FFT solver.
\item The swept surfaces are used by the solver to update the eigenstrain field: $\boldsymbol{\varepsilon}^{\star}(t) + \delta\strain^{\star}$ (Figure~\ref{fig7}). The updated eigenstrain $\strain^{\star}(t+\delta t)$ becomes the initial condition for the FFT solver, as well as the received mechanical loading (stress or strain), assuming periodic boundary conditions.
\item The mechanical equilibrium is computed either in a single step if elastic constants are homogeneous, or by a fixed point algorithm if inhomogeneous. The stress field $\stress(t+\delta t)$ is computed at the nodes of the fcc grid, compliant with \micromegas{} grid. Therefore, the stress at the segment midpoints are directly accessible and are stored in an array sent to \micromegas{} by MPI communication.
\item Once the array in which the stress field received, \micromegas{} computes the Peach--Koehler force on each segment and determines their displacement.
\end{enumerate}

\section{Results} \label{results}

\subsection{Static dislocation loop}

We have first considered the case of a single static $(1\overline{1}1)[\bar{1}01]$ dislocation loop in pure Cu, assuming homogeneous and isotropic elastic constants $\mu = \SI{46.61}{\mega\pascal}$ and $\nu = 0.33$. The amplitude of the Burgers vector is $b = \SI{2.5525}{\angstrom}$, and the applied stress $\sigma^{a}_{33} = \SI{300}{\mega\pascal}$. The grid contains $256^3$ nodes, with grid spacing $d = \SI{2.4366}{\nano\meter}$. The loop, described by a sharp eigenstrain field (cf.\ Section~\ref{subsection_sharp_description}), is rectangular, $\SI{206}{\nano\meter}$ long in the edge direction, $\SI{180}{\nano\meter}$ large in the screw direction, and 1 voxel in thickness (Figure~\ref{fig6}).

\begin{figure*}[h]
\includegraphics[width=0.75\textwidth]{Figure_6.pdf} \caption{Representation of a $(1\overline{1}1)[\overline{1}01]$ dislocation loop by a sharp eigenstrain field with rhombohedral voxels (easier to compute).} \label{fig6}
\end{figure*}

\begin{figure}[ht]
\includegraphics[width=0.85\textwidth]{Figure_7.pdf} \caption{Stress components in different planes associated with a $(1\overline{1}1)[\overline{1}01]$ dislocation loop.} \label{fig7}
\end{figure}

The stress components $\sigma_{11}$ and $\sigma_{13}$ (indices $1,2,3$ corresponding respectively to $[100]$, $[010]$ and $[001]$ cristallographic directions) in the $(1\bar{1}1)$ slip plane are plotted in Figure~\ref{fig7}a-b. The stress component $\sigma_{11}$ around the screw segments in the $(\bar{1}01)$ plane and around the edge segments in the $(121)$ plane is plotted in Figure~\ref{fig7}c-d. It is worth stressing that no oscillation, that plague other FFT solvers~\cite{1998_Moulinec,2015_Willot}, are observed close to the sharp variations of the eigenstrain. We have also plotted $\sigma_{11}$ along $[121]$ (red arrow) across the screw segments in Figure~\ref{fig8}a, and $\sigma_{33}$ along the $[\bar{1}01]$ (green arrow) across the edge segments in Figure~\ref{fig8}b. These profiles are compared to the analytical solution of the non-singular theory of~\cite{2006_Cai} with a regularization width $r_c\approx d$ (Section~\ref{DCM}). The agreement is excellent, whatever the stress component and whatever the discretization. In particular, the characteristic antisymmetric profiles, close to the dislocation line, at the core region, edge or screw, indicate that the octahedral stencil fulfills the fcc symmetries. Although a small difference in the stress peaks can be observed, due to the difference in the distribution of the eigenstrain, the stress value on the dislocation segment, used to compute the Peach--Koehler forces, is the same between the present work and the analytical solution. Finally, it can be noted that the same excellent agreement has been found for different loops for all fcc slip systems.

\begin{figure}[h]
\includegraphics[width=\textwidth]{Figure_8.pdf} \caption{Stress components in (a) edge direction and (b) screw direction compared to the analytical solution of the non-singular theory of~\cite{2006_Cai}.} \label{fig8}
\end{figure}

%\FloatBarrier

\subsubsection{Cross-slip} \label{subsubsec:cross_slip}

In a second step, we have examined the case of two loops with the same Burgers vector joined by their opposite screw segments (red line, Figure~\ref{fig9}a). The elastic constants, grid and applied stress are the same as in the previous calculation. The two loops are rectangular with \SI{138}{\nano\meter} long screw segments and \SI{122}{\nano\meter} long edge segments. This static configuration mimics a cross-slip event, with one loop in the primary glide plane and the second one in the cross-slip plane.

\begin{figure}[ht]
\includegraphics[width=0.87\textwidth]{Figure_9.pdf} \caption{(a) Cross-slip configuration described by two loops joined by their opposite screw segments. (b) $\sigma_{22}$ and $\sigma_{33}$ in plane $(011)$ normal to the screw segments and to the fictitious junction. (c) $\sigma_{22}$ along $[21\overline{1}]$ passing through the fictitious junction.} \label{fig9}
\end{figure}

Along the red segment of Figure~\ref{fig9}a, the total eigenstrain is the sum of the eigenstrains of the two abutting loops weighted by their respective fractions in the voxels containing the red segment (Section~\ref{subsection_sharp_description}):
\begin{equation}
\arraycolsep=4pt
A \strain^\star = \frac{1}{2}
\begin{pmatrix}
\;	0 & 1 & 1 	\;\\
\;	1 & -2 & 0 	\;\\
\;	1 & 0 & 2	\;
\end{pmatrix} + \frac{1}{2}
\arraycolsep=2pt
\begin{pmatrix}
0 & -1 & -1 \\
-1 & -2 & 0 \\
-1 & 0 & 2
\end{pmatrix} =
\arraycolsep=4pt
\begin{pmatrix}
\;	0 & 0 & 0	\;\\
\;	0 & -2 & 0	\;\\
\;	0 & 0 & 2	\;
\end{pmatrix}
\end{equation}
where $A = \frac{2V_d \abs{\myvec{m}}\times \abs{\myvec{n}}}{b S_h}$ (cf.~\eqref{discrete_eigenstrain}).

It can be noted that the diagonal components are continuous in the two loops, whereas the non-diagonal ones vanish in the voxels containing the red segment, corresponding to the annihilation of the two opposite dislocation segments. As a consequence and as expected, the fictitious junction does not generate any stress, as clearly shown in Figure~\ref{fig9}b-c. It is worth stressing that this important feature has been obtained because the octahedral stencil fulfills the necessary symmetries around the fcc slip planes. Indeed, other finite-difference stencils have been tested without the proper symmetries: they have generated spurious local stresses along the red segment that could affect the dynamics of dislocations.

%\FloatBarrier

\subsection{Junction formation}

It has been shown that the formation of junctions between dislocations is a fundamental elementary mechanism to understand and predict plastic deformation. In DDD codes, the junction formation is driven by the elastic interactions. Consequently, we have analyzed the consequence of the coupling on this mechanism.

For that purpose, we have used the same parameters as in the previous calculations, except for a larger system with $512^3$ nodes and a lower constant stress $\sigma^{a}_{33} = \SI{152}{\mega\pascal}$ applied during $t_f = \SI{20}{\nano\second}$. The configuration is described in Figure~\ref{fig10}: it is composed of two rectangular loops with \SI{344.6}{\nano\meter} long screw segments and \SI{298.4}{\nano\meter} edge segments. Their Burgers vectors are, respectively, $\myvec{b}_1=[10\bar{1}]$ and $\myvec{b}_2=[011]$. Only the red segments can move, and the green segments are kept immobile to prevent the planar loops from closing before the formation of the junction between the red segments. A Lomer--Cottrell junction is expected to form with a Burgers vector $\myvec{b}_1 + \myvec{b}_2 = [110]$.

\begin{figure}[hb]
\includegraphics[width=0.4\textwidth]{Figure_10.pdf} \caption{Configuration of two planar loops forming a junction.} \label{fig10}
\end{figure}

As shown in Figure~\ref{fig11}a, the coupled DDD-FFT model (i) properly builds the eigenstrain field with the swept area algorithm discussed in Sections~\ref{subsection_sharp_description}--\ref{subsection_eigenstrain_dynamics}, and (ii) a junction formation is observed, with the line direction $\myvec{\xi}= [\bar{1}10]$, consistent with a Lomer--Cottrell lock formation.

Moreover, Figure~\ref{fig11}b shows that the stress field generated by the junction is perfectly smooth and free of any oscillation (for brevity, only components $\sigma_{11}$ and $\sigma_{12}$ are shown but the other components exhibit similar features): this demonstrates that the octahedral solver is relevant to handle dislocations described by sharp eigenstrain.

\begin{figure}[ht]
\includegraphics[width=0.8\textwidth]{Figure_11.pdf} \caption{(a) Component $\varepsilon^{\star}_{33}$ of the 3D eigenstrain. (b) Stress field in a $(1\overline{1}0)$ cross section perpendicular to the junction line, at $t_f$.} \label{fig11}
\end{figure}

In Figure~\ref{fig12}, we have compared the final shapes of the dislocations obtained with our coupled DDD-FFT model (in blue) and with the classical DDD (in red). The agreement is quite good, although the junction is shorter with the coupled DDD-FFT. Different reasons can be proposed to explain this difference. First, the periodic boundary conditions are applied to different box geometries: cubic for \micromegas{} and rhombohedral for the coupled DDD-FFT. Therefore, in the DDD-FFT coupling replicas of the dislocation loop are closer with respect then in the classical DDD calculation and not uniformly distributed: in some direction, due to the rombohedral cell, the periodicity is quite reduced compared to a cubic box, as we can clearly notice in Figure~\ref{fig12}. This generates an effective Peach--Koehler force acting on the red segment, causing a discrepancy between the two calculations. Second, we have shown that the stress peaks at the dislocation core are slightly different (Figure~\ref{fig8}), affecting, again, the force acting on dislocations.

\begin{figure}[hb]
\includegraphics[width=\textwidth]{Figure_12.pdf} \caption{Junction formation (snapshots at $t=0$ and $t=t_f$): comparison between \micromegas{} (red) and the coupled DDD-FFT calculations (blue).} \label{fig12}
\end{figure}

\subsection{Anisotropic medium}

One major advantage of our coupled DDD-FFT model is that it is straightforward to consider anisotropic elasticity. Indeed, the description of the dislocations as sharp eigenstrain fields discards the need of the short-range correction used in the classical DCM coupled with a Finite Element elastic solver. Moreover, it is worth noting that the collocation of the grids of the DDD code and the FFT solver avoids any data transfer as would necessitate other FFT solvers based on cubic grids, e.g.~\cite{2016_Schneider,2016_Valdenaire}.

To illustrate the capacity of our coupling with respect to anisotropic elasticity, we have computed first the stress field of a static dislocation loop and then the equilibrium shape of such a loop, at constant plastic strain ($\delta\strain^p=0$). We have used the same grid as in the previous section~($512^3$ nodes and grid spacing $d=\SI{2.4366}{\nano\meter}$). Initially, the $(11\bar{1})[101]$ dislocation loop with $b = \SI{2.5525}{\angstrom}$ is rectangular with $\SI{344.6}{\nano\meter}$ long screw segments and $\SI{298.4}{\nano\meter}$ long edge segments. The initial stress is $\sigma^0_{33} = \SI{170}{\mega\pascal}$, and the time step is $\delta t = \SI{1.1e-11}{\second}$. To enhance the effect of elastic anisotropy, the local line tension contribution relying on isotropic elasticity has been discarded.

\begin{table}[h!]
\caption{Cubic elastic constants in \si{\giga\pascal} for three Zener anisotropy ratios $A$.} \label{tab:anisotropic_elastic_constant}
\begin{tabular}{cccc}
\hline
$\boldsymbol{A}$ & $\boldsymbol{1}$ & $\boldsymbol{3}$ & $\boldsymbol{10}$ \\
\hline
$C_{11}$ & 200 & 164 & 144 \\
$C_{12}$ & 100 & 118 & 128 \\
$C_{44}$ & 50 & 68 & 78 \\
\hline
\end{tabular}
\end{table}

Three cubic rigidity tensors have been considered, that all comply with the same equivalent isotropic constants $\lambda$ and $\mu$:
\begin{equation}\label{voigt_average}
\begin{aligned}
C_{11} &= \lambda + 2\mu - \tfrac{2}{5} H, \\
C_{12} &= \lambda + \tfrac{1}{5} H, \\
C_{44} &= \mu + \tfrac{1}{5} H,
\end{aligned}
\end{equation}
where the $C_{ij}$ are the three independent components of the rigidity tensor in Voigt notation, and where the anisotropy factor $H$ and the Zener anisotropy ratio $A$ are respectively~\cite{1967_Hirth} $H = 2C_{44} - C_{11} + C_{12}$ and $A = 2C_{44}/(C_{11}-C_{12})$. The sets of values are reported in Table~\ref{tab:anisotropic_elastic_constant} and correspond to $\lambda=\SI{100}{\giga\pascal}$ and $\mu=\SI{50}{\giga\pascal}$ ($\nu=0.33$).

\begin{figure}[h]
\includegraphics[width=0.5\textwidth]{Figure_13.pdf} \caption{Equilibrium shapes of planar loops for three anistropy ratios $A=1$ (isotropic), $A=3$ and $A=10$.} \label{fig13}
\end{figure}

The equilibrium shapes are shown in Figure~\ref{fig13}. As expected, the larger $A$, the larger the deviation from the reference elliptical shape ($A=1$), with longer screw segments along $[101]$. It seems identical as what was observed in BCC metals in~\cite{2011_Aubry}.

In Figure~\ref{fig14}, we have plotted the non-zero stress components in a $(101)$ section normal to the screw segments. The indices ``$x,y,z$'' refer respectively to crystallographic directions $[101]$, $[1\bar{2}\bar{1}]$ and $[11\bar{1}]$. This choice of basis has been made to compare the results to the literature. The stress loading is retrieved to keep only the contribution of each dislocation segment. The usual two lobes pattern is obtained for the isotropic case ($A=1$). When $A$ increases, the stress lobes can be seen to deform and favour certain directions.

In Figure~\ref{fig15} different stress components are plotted in a $(\bar{1}21)$ section normal to the edge segments. The indices ``$x,y,z$'' refer this time respectively to cristallographic directions $[\bar{1}21]$, $[101]$ and $[11\bar{1}]$. Components $\sigma_{xz}$ and $\sigma_{yz}$ are not shown here, we observe zero values in the isotropic case, in accordance with the literature, and non-zero values in both anisotropic cases. For both screw and edge segments, the shape of stress field is consistent with analytical solutions expressed in~\cite{2024_Lazar} for Al (isotropic) and Cu ($A \approx 3$).

\begin{figure*}[!t]
\includegraphics[width=\textwidth]{Figure_14.pdf} \caption{Non-zero components of stress across two opposite screw segments expressed in glide system basis for three anisotropy ratios $A=1$ (isotropic), $A=3$ and $A=10$.} \label{fig14}
\end{figure*}

\begin{figure*}[h!]
\includegraphics[width=\textwidth]{Figure_15.pdf} \caption{Stress components $\sigma_{xy}$, $\sigma_{yy}$ and $\sigma_{zz}$ of stress across two opposite edge segments in glide system basis for three anisotropy ratios $A=1$ (isotropic), $A=3$ and $A=10$.} \label{fig15}
\end{figure*}

%\FloatBarrier

\subsection{Interactions with an inclusion}

Finally, we illustrate the ability of our DDD-FFT model to handle the interaction of dislocations with microstructures. We have considered the case of a soft inclusion $\mathcal{I}$ with isotropic elastic constants, two orders of magnitude lower than in the matrix, i.e.\ $\lambda = \SI{1}{\giga\pascal}$ and $\mu = \SI{0.5}{\giga\pascal}$. As discussed in Section~\ref{fourier_operators}, the chosen stencil is not unconditionally stable, which limits the convergence of the iterative solver. For this reason, we do not consider the infinite contrast case, such as voids or rigid inclusions. The inclusion is a cube with $43.66~\si{\cubic\nano\meter}$ edges aligned with the $\langle 100\rangle$ fcc directions. The inclusion is described by a field $\phi(\myvec{x})$ such that $\phi(\myvec{x})=1, \ \forall\myvec{x}\in\mathcal{I}$ and $\phi(\myvec{x})=0, \ \forall\myvec{x}\notin\mathcal{I}$. To avoid unstable modes $\langle 100\rangle_\mathcal{E}$, the field has been smoothed by convolution with an appropriate Gaussian function so that the interface is diffused and spans 6 node spacings. This strategy lies at the heart of phase-field modeling of microstructure evolution with internal stresses~\cite{AmmarAppolaire09b,AmmarAppolaire21,CotturaAppolaire21}. It has also been proposed for numerical homogenization of composite material, e.g.~\cite{MorinBrenner21}.

The elastic properties are simply weighted by $\phi$. It is an arbitrary choice, and other possibilities based on homogenization schemes (Voigt, Reuss, rank-1 or laminates) have been proposed in different contexts, such as the phase field approach~\cite{AmmarAppolaire09b,SchneiderTschukin15} and numerical homogenization studies~\cite{MareauRobert17}. Depending on the scheme, stresses can differ significantly at the interfaces. The rhombohedral system is discretized using $512^3$ nodes with a grid spacing $d = \SI{2.4366}{\nano\meter}$. We have considered that the precipitate is contained in a single rhombohedral grain whose boundaries are impenetrable barriers to dislocations. The system is submitted to tension along $[001]$-axis, with an adaptative increment of stress $\langle\sigma_{33}\rangle$ depending on dislocation velocity for $t_f = 1400\,\delta t$ where $\delta t = \SI{1e-11}{\second}$.
       	 
The pattern of dislocations around the inclusion achieved at $t_f$ is shown in Figure~\ref{fig16}a. Some segments are stuck at the boundaries of the system due to the impermeable grain boundaries. Moreover, the inclusion is free of dislocation because the stress almost vanishes in the inclusion due to the very low elastic constants. Therefore, in this case, choosing to interpolate the elastic constants should not affect the results.

Finally, in Figure~\ref{fig16}b we have plotted the component $\sigma_{33}$ of the stress field in a $(100)$ cross section containing the center of the cubic inclusion. The accumulation of dislocations close to the impermeable boundaries results locally in larger stresses. The benefit of using sharp eigenstrain fields to describe the dislocations is clearly visible with the local stress patterns around each dislocation. The number of dislocations that can be handled with the present coupling is significantly larger than that with the previous versions of the DCM. In addition, the fields are free of the oscillations generated by the simpler~\cite{1998_Moulinec} or more popular stencils~\cite{2008_Willot}.

\begin{figure}[!h]
\includegraphics[width=\textwidth]{Figure_16.pdf} \caption{(a) Dislocations pattern around a soft cubic inclusion (blue). The different colors indicate the different slip systems. (b) Stress component $\sigma_{33}$ in a $(100)$ section crossing the inclusion at its center (isosurface $\phi=0.5$).} \label{fig16}
\end{figure}

\section{Conclusion}

In this work we have detailed an FFT solver working with a discrete theory of Green's operators that fulfills this requirement. Firstly, the mesh of the FFT solver is chosen to coincide exactly with the fcc network of the DDD, in the context of the study of fcc materials. The integration points presented here form an octahedral cell, so, the solver is called octahedral solver. It is used to solve the various artefacts due to symmetry around $\{111\}$ slip planes. The discrete differential operators of this solver do not ensure unconditional numerical stability. Nevertheless, it is possible to solve the mechanical equilibrium of fcc dislocations in $\{111\}$ planes represented as a sharp field without generating numerical oscillations. The coupling between the DDD and the latter solver has been developed. Because of the lack of mathematical stability, other types of inclusions than fcc dislocations may produce oscillations. However, this coupling is interesting in many ways since it allows the study of plasticity mechanisms in anisotropic materials, without short-range correction of the internal stress field. The results for anisotropic elasticity are very promising, allowing fine analysis of the stress field near the dislocation core and coinciding very well with the literature and known analytical formulas. The interaction between dislocations and diffuse inclusions such as precipitates can also be handled. This work has also established good prospects for the study of dislocation-crack interaction~\cite{2022_Eon}.

\section*{Declaration of interests}

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography

\end{document}