# Prediction of Yield Surface of Single Crystal Copper from Discrete
Dislocation Dynamics and Geometric Learning

Wu-Rong Jian  Mian Xiao  WaiChing Sun<EMAIL_ADDRESS>Wei Cai
<EMAIL_ADDRESS>Department of Mechanical Engineering, Stanford University,
Stanford CA, 94305, USA  Department of Civil Engineering and Engineering
Mechanics, Columbia University, 614 SW Mudd, 4709, New York, NY 10027, USA

###### Abstract

A yield surface of a material is a set of critical stress conditions beyond
which macroscopic plastic deformation begins. For crystalline solids, plastic
deformation occurs by the motion of dislocations, which can be captured by
discrete dislocation dynamics (DDD) simulations. In this paper, we predict the
yield surfaces and strain-hardening behaviors using DDD simulations and a
geometric manifold learning approach. The yield surfaces in the three-
dimensional space of plane stress are constructed for single-crystal copper
subjected to uniaxial loading along the $[100]$ and $[110]$ directions,
respectively. With increasing plastic deformation under $[100]$ loading, the
yield surface expands nearly uniformly in all directions, corresponding to
isotropic hardening. In contrast, under $[110]$ loading, latent hardening is
observed, where the yield surface remains nearly unchanged in the orientations
in the vicinity of the loading direction itself, but expands in other
directions, resulting in an asymmetric shape. This difference in hardening
behaviors is attributed to the different dislocation multiplication behaviors
on various slip systems under the two loading conditions.

††journal: arXiv

##  1 Introduction

Understanding crystal plasticity in terms of fundamental physics has been a
long-standing goal in computational materials science. Since the proposal of
crystal dislocations (Taylor, 1934; Polanyi, 1934; Orowan, 1934) and their
observation by transmission electron microscopy (Hirsch et al., 1956), it has
been well-established that the plastic deformation behaviors of crystalline
materials are controlled by the motion of dislocations (Hirth and Lothe,
1982). The discrete dislocation dynamics (DDD) simulation method (Amodeo and
Ghoniem, 1990; Devincre and Kubin, 1997; Arsenlis et al., 2007) has been
developed to establish the connection between the microscopic motion of
individual dislocations and the macroscopic stress-strain behavior of the
single-crystalline material. A fundamental concept in describing the plastic
deformation behavior of materials is the yield surface (Meyers and Chawla,
2008). When the local stress at a material point is within the yield surface,
the deformation is purely elastic. By contrast, when the local stress reaches
the yield surface, plastic deformation begins. Furthermore, the strain-
hardening behavior corresponds to the change of the yield surface with
increasing plastic strain. For example, isotropic hardening corresponds to the
uniform expansion of the yield surface in all directions, while kinematic
hardening corresponds to the translation of the yield surface in the stress
space without changing its size and shape (Hill, 1998). The Bauschinger effect
(Bauschinger, 1886), in which the plastic deformation causes yield stress in
the reverse loading direction to become lower, is a manifestation of non-
isotropic hardening behavior. Given the central role of yield surface in
plasticity and solid mechanics, it would be a natural goal to predict the
yield surface of single crystals by DDD simulations. However, to our best
knowledge, there have not been predictions of yield surface from DDD
simulations to date.

There have been several challenges that have prevented DDD predictions of
yield surfaces for single crystals. First, DDD simulations have been
computationally very expensive. On the one hand, in order to predict the
plastic deformation behavior of a crystal, the motion of dislocations in a
large enough simulation cell needs to be followed for a long enough time to
accumulate a significant level of plastic strain. On the other hand,
interactions between nearby dislocations require a very small time step to
ensure numerical stability. Second, to construct the yield surface for a given
computational sample (i.e. a dislocation configuration) by brute-force DDD
simulations would require a very large number of DDD simulations, each with a
different loading (stress) orientation. Given that the stress conditions span
a six-dimensional (6D) parameter space, the total number of DDD simulations is
overwhelming if a uniform sampling is used. Dealing with crystal plasticity,
we cannot reduce the stress space to the three-dimensional (3D) space of the
principal stresses, because the orientation of the stress with respect to the
slip systems of the crystal matters for the forces on dislocations and the
yielding behavior. Third, we need a framework to express the yield surface
that is flexible enough to fit to the data from DDD predictions, as well as
providing mathematically well-behaved interpolation in-between. Ideally, the
framework should also inform us where more sampling is needed in the stress
space where additional DDD simulations can be performed to provide the most
valuable information about the yield surface.

In this work, we construct the yield surfaces of single-crystal Cu from DDD
simulations with the following approaches to address the above challenges.
First, recent progress on subcycling integrator (Sills et al., 2016) and GPU
implementation (Bertin et al., 2019) has made DDD simulations much more
efficient than before. For single crystal Cu, plastic strain on the order of
$1\%$ strain can be reached using a single GPU over the time period of a week,
although the strain rate still needs to be kept high (e.g. $10^{3}$ s-1) in
these simulations. In this work, we will keep our effective strain rate to be
about $10^{3}$ s-1. Hence, our predicted yield surface is not the same as the
one typically might have in mind, i.e. for quasi-static loading conditions.
Our work represents a first attempt at this problem, and the strain rate
effect on the yield surface will be examined in the future. In addition, we
show that if we just need to determine the yield point for a given loading
direction, the DDD simulation can be much shorter than what is needed for
determining the strain hardening rate. Second, in order to reduce the total
number of DDD simulations, here we will limit our scope to the plane-stress
condition, where only three stress components are non-zero. This means that we
will construct a 3D “cross-section” of the full yield surface which lives in
the 6D space. The yield surface in the 3D plane-stress conditions is also much
easier to visualize, and we show that quite interesting insights can be gained
by examining how it evolves while the crystal is plastically deformed along
different loading directions. Third, we use the geometric prior method (Xiao
and Sun, 2022) to construct the yield surface from the DDD simulation data,
where patches each having a local coordinate system provide smooth local
descriptions of the yield surface, respectively. Then, the patches overlap
consistently and offer a general description of any surface. The geometric
prior framework can provide ways to evaluate local features on the yield
surface, e.g. local curvature, which can be used to decide where more data
should be collected from additional DDD simulations. Prior to the geometric
prior method, the non-uniform rational B-spline (NURBS) method (Coombs et al.,
2016; Coombs and Motlagh, 2017, 2018) and the level-set function (Vlassis and
Sun, 2021, 2022) had also been developed to construct yield surfaces. Compared
to the geometric prior method that employs various trainable neural networks
in different patches, both NURBS and level-set approaches use a single learned
function to represent the yield surface, and the resultant models are
difficult to update when new data are supplemented from DDD simulations.
Hence, they are not suitable for an active learning environment where the
learned yield surface must be easily updated (i.e., without re-training the
entire model from scratch) whenever new data is presented. While we
demonstrate yield surface in the 3D plane-stress space, our geometric prior
framework to construct yield surface is equally applicable to the full 6D
stress space.

The rest of this paper is organized as follows. First, the DDD simulation
setup and the procedure of yield stress data extraction from DDD simulation
are described in Section 2.1, followed by the details of our yield surface
construction framework using geometric prior in Section 2.2. We present the 3D
yield surfaces of various dislocation configurations of single-crystal Cu upon
plane stress loading in Section 3.1 and the 2D yield loci for uniaxial loading
in Section 3.2. The results show that the yield surface evolves in different
modes during the strain hardening of single-crystal Cu along [100] (isotropic
hardening) and [110] (latent hardening), respectively. Section 4.1 gives an
explanation of this observation based on the different dislocation
multiplication behaviors on various slip systems for these two loading
orientations. The outlooks for DDD simulations and the geometric prior method
for the yield surface construction of single crystals are provided in Section
4.2 and Section 4.3, respectively. The conclusions are given in Section 5.

##  2 Methodology

###  2.1 Data collection from DDD

The DDD simulations are conducted using the open-source ParaDiS program
(Arsenlis et al., 2007; ParaDiS, 2023), which utilizes the recently developed
sub-cycling time integration scheme (Sills et al., 2016) and its Graphics
Processing Units (GPU) implementation (Bertin et al., 2019). The material
parameters in our DDD simulations correspond to single crystal copper, e.g.,
Burgers vector magnitude $b=0.255$ nm, Poisson’s ratio $\nu=0.324$, and shear
modulus $\mu=54.6$ GPa. During all DDD simulations, a linear mobility law with
drag coefficient $B=1.56\times 10^{-5}\mathrm{~{}Pa}\cdot\mathrm{s}$ is
applied to glissile dislocations. Cross-slip is not allowed, and dislocation
junctions are only allowed to move along their line directions using zipping
or unzipping mode. All the DDD simulation parameters are also summarized in
Table 1. Due to the large amount of simulations necessary to obtain a yield
surface in the general 6D stress space, here we limited our attention to plane
stress conditions, where only $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{xy}$ may
be non-zero.

The dimensions for the simulation model are $\sim$ 15 $\mu$m $\times$ 15
$\mu$m $\times$ 15 $\mu$m. To build the initial dislocation configuration for
the subsequent loading simulations, straight dislocation lines on the
$\frac{1}{2}\langle 110\rangle\\{111\\}$ slip systems are introduced randomly
to the simulation box, where periodic boundary conditions (PBCs) are applied
along all three directions. The dislocation configuration is then relaxed,
resulting in a dislocation density of $\rho_{0}\approx 1.2\times
10^{12}\mathrm{~{}m}^{-2}$. The relaxed configuration, shown in Figure 1 and
also used in our previous works (Akhondzadeh et al., 2020, 2021), is labelled
config-0, and considered as the reference configuration in this work.

Figure 1: The initial dislocation configuration (config-0) with a dislocation
density of $\rho_{0}\approx 1.2\times 10^{12}\mathrm{~{}m}^{-2}$ before
loading in DDD simulation.

To construct the yield surface of a given (relaxed) dislocation configuration,
we perform a set of DDD simulations, each at a different stress orientation,
under a constant stress rate. We choose a coordinate system whose axes are
aligned with the cubic axes of the single crystal, i.e. the $x$, $y$, $z$-axes
are along the $[100]$, $[010]$, $[001]$ directions, respectively. For
simplicity, we focus on the plane stress condition (in the $x$-$y$ plane, i.e.
the $(001)$ crystallographic plane) where the non-zero stress components are
limited to $\sigma_{xx}$, $\sigma_{yy}$ and $\sigma_{xy}$. In this case, the
magnitude and orientation of the stress tensor $\bm{\sigma}$ can be visualized
by considering a 3D vector $\vec{\sigma}$ = ($\sigma_{xx}$, $\sigma_{yy}$,
$\sqrt{2}\sigma_{xy}$). This is because the magnitude of stress as given by
the Frobenius norm

$\sigma=\|{\bm{\sigma}}\|=\sqrt{\sigma_{xx}^{2}+\sigma_{yy}^{2}+2\sigma_{xy}^{2}}\,,$
(1)

is the same as the length of the vector $\vec{\sigma}$. Therefore, for every
plane stress state with non-zero magnitude, we can define a unit vector,
$\hat{\sigma}=\vec{\sigma}/\sigma$, which specifies the stress orientation.
Hence all stress orientations in plane stress correspond to a unit sphere in
the 3D space of ($\sigma_{xx}$, $\sigma_{yy}$, $\sqrt{2}\sigma_{xy}$). To
construct the yield surface under plane stress, we thus need to sample the
stress orientation $\hat{\sigma}$ from this unit sphere, and for each chosen
$\hat{\sigma}$ increase the stress magnitude $\|{\bm{\sigma}}\|$ linearly with
time until yielding is detected.

Figure 2: A total of 92 points distributed nearly uniformly on the unit sphere
are shown as vertices of a polyhedron generated by the icosphere module in
Python (Dahl, 2023). The red circle corresponds to uniaxial tensile loading
along different directions in the $x$-$y$ plane.

Our initial sampling of the plane stress space corresponds to 92 points nearly
uniformly distributed on the unit sphere, which are the vertices of the
polyhedron shown in Figure 2. The coordinates of the sampling points are
obtained using the icosphere module in Python (Dahl, 2023). The uniaxial
loading conditions along different directions in the $x$-$y$ plane are shown
as the red circle on the unit sphere in the 3D space of plane stress.

For each stress orientation shown in Figure 2, we performed a DDD simulation
in which the stress orientation remains fixed and the stress magnitude
increases at a constant rate $\dot{\sigma}=10^{13}\,{\rm Pa}\cdot{\rm
s}^{-1}$. The magnitude of the stress rate is chosen here so that the stress-
strain response of the single crystal upon uniaxial loading condition matches
those at a strain rate of $10^{3}\,{\rm s}^{-1}$, which is a common strain
rate used in DDD simulations of single crystal Cu (Sills et al., 2018; Bertin
et al., 2019; Akhondzadeh et al., 2020, 2021). For example, Figure 1 shows
that for uniaxial tensile loading along $[100]$ crystal orientation, the
stress-strain curves for constant stress rate of $10^{13}\,{\rm Pa}\cdot{\rm
s}^{-1}$ and constant strain rate of $10^{3}\,{\rm s}^{-1}$ are consistent
with each other up to the yielding point.

Figure 3: (a-b) The stress-strain curves for the uniaxial tensile loading at
strain rate $10^{3}\,{\rm s}^{-1}$ along [100] and [110] directions,
respectively. The blue and red dots on the curves represent the unloading
points at the strains of 0.07% and 0.22%, where dislocation configurations are
extracted and relaxed, respectively. The yield point is determined as the
intersection of the curve with the elastic curve offset of 0.022%, as denoted
by the orange arrows. (c-d) show the corresponding relaxed dislocation
configurations from the [100] loading, respectively. (e-f) show the
corresponding relaxed dislocation configurations corresponding from the [110]
loading curve, respectively.

To plot the stress-strain curves for a general plane stress loading, we define
the strain (in the horizontal axis) as the dot product between the stress and
strain tensors,

$\varepsilon^{\rm pl}=\hat{\sigma}_{xx}\varepsilon_{xx}^{\rm
pl}+\hat{\sigma}_{yy}\varepsilon_{yy}^{\rm
pl}+2\hat{\sigma}_{xy}\varepsilon_{xy}^{\rm pl}$ (2)

where $\hat{\sigma}_{xx}=\sigma_{xx}/\|{\bm{\sigma}}\|$,
$\hat{\sigma}_{yy}=\sigma_{yy}/\|{\bm{\sigma}}\|$,
$\hat{\sigma}_{xy}=\sigma_{xy}/\|{\bm{\sigma}}\|$. $\varepsilon_{xx}^{\rm
pl}$, $\varepsilon_{yy}^{\rm pl}$ and $\varepsilon_{xy}^{\rm pl}$ are the
plastic strain components. This definition reduces to the normal plastic
strain along the loading direction for the case of uniaxial loading. In this
work, the yield point of plane stress loading is determined by
$\varepsilon^{\rm pl}=0.022\%$. This yield criterion is the same as using the
strain offset in the case of uniaxial loading, when applied to the stress-
strain curve in which the strain is the total strain, as shown in Figure 3(a)
and (b). The choice of $0.022\%$ offset strain is similar to those used in
previous experimental work (Franciosi et al., 1980).

To investigate how the yield surface evolves by strain hardening, we also
computed the yield points along the stress orientations defined in Figure 2
for four more dislocation configurations. These configurations are obtained
from DDD simulations starting from config-0 as initial condition, and
subjected to uniaxial tensile loading conditions with strain rate
$10^{3}\,{\rm s}^{-1}$, along the $[100]$ and $[110]$ directions,
respectively. For each loading orientation, two configurations are extracted
at the strains of 0.07% and 0.22%, and unloaded to zero stress and relaxed,
leading to config-1 and config-2, respectively. These four dislocation
configurations are plotted in Figure 3(c), (d), (e) and (f), respectively.
From the normal stress-strain curves in Figure 3(a) and (b), it can be seen
that the $[100]$ loading orientation has a higher strain hardening rate than
the $[110]$ loading orientation. This is because for $[100]$ loading, all 4
slip planes and 8 out of 12 slip systems are active, and their intersections
lead to a high strain hardening rate. In contrast, for $[110]$ loading, only 2
slip planes and 4 slip systems are active. However, the stress-strain curves
in Figure 3(a) and (b) only reveal how the yield point changes along the same
stress orientation as the direction of applied loading. In this work, our goal
is to determine how the yield surface changes with plastic deformation, i.e.
how the yield point changes in all stress orientations most of which do not
coincide with the direction of applied loading.

###  2.2 Construction of yield surface with geometric prior method

We use the method of geometric prior (Xiao and Sun, 2022) to construct the
yield surface ($\mathcal{S}$) under plane stress as a manifold atlas in
$\mathbb{R}^{3}$, based on the yield point data predicted by DDD simulations
along a discrete set of stress orientations ($\mathcal{X}$). The method of
geometric prior also allows us to iteratively improve the yield surface
representation, by providing suggestions on stress orientations where new DDD
simulations can be performed to improve the data density in regions most
needed.

We characterize the mathematical representation of the yield surface as a
manifold atlas rather than the traditional implicit yield function
representation (Hill, 1998). One can draw an analogy between a manifold atlas
and a piece-wise function: (1) the domain of the piece-wise function is
divided into multiple pieces; similarly, the manifold atlas also partitions
the surface into multiple patches. (2) On each piece of the domain, the piece-
wise function is defined continuously, either explicitly or implicitly;
correspondingly, each patch of the manifold atlas is associated with a
coordinate chart, a neural network-trained nonlinear mapping that enables us
to provide local coordinates of the patch within the yielding manifold. An
example of this representation is shown in Figure 4.

Figure 4: Sketch of the geometric prior manifold representation. The Riemann
surface is divided into multiple pieces (patches), and each patch is
associated with a coordinate chart as a mapping from a regular parametric
domain.

After constructing the initial manifold atlas based on DDD simulation data of
92 stress orientations, we can compute local geometric information, e.g.
Gaussian curvature, everywhere on the yield surface. We then identify regions
where local refinement is desirable. These are regions where the yield surface
intersects a bounding sphere of a chosen radius and centered at points of
local maxima of Gaussian curvature. A point cloud is generated in each of
these refinement regions using Poisson disk sampling. Additional DDD
simulations are then performed to predict yield points at stress orientations
corresponding to these new sampling points. The data is combined with the
previously obtained DDD data to obtain a refined manifold atlas representing
the updated yield surface. The process can be iterated to sequentially refine
the yield surface representation, as illustrated in Figure 5.

Figure 5: Major steps in the reconstruction of a yield surface.

####  2.2.1 Geometric prior method for yield surface construction

Here we provide some specifics of the geometric prior method that is applied
to constructing the yield surface $\mathcal{S}$ under plane stress as a
manifold atlas in $\mathbb{R}^{3}$; for more details please refer to (Xiao and
Sun, 2022). To generate the patch discretization of $\mathcal{S}$, we
prescribe a set of anchor points $\vec{x}_{p}$ on $\mathcal{S}$. Then, a
surface patch $\mathcal{S}_{p}$ is defined by the intersection of the
neighborhood ball at $\vec{x}_{p}$ and $\mathcal{S}$:
$\mathcal{S}_{p}=\mathcal{S}\cap\mathcal{B}_{p},\mathcal{B}_{p}=\\{\vec{x}\in\mathbb{R}^{3}|\>\>||\vec{x}-\vec{x}_{p}\|<\epsilon\\}$,
where $\epsilon$ is the radius of the neighborhood ball and $\|\cdot\|$ is the
L2-norm. We then define the parametrization of the local coordinate charts as
$\phi_{p}(\vec{v}):V\rightarrow\mathbb{R}^{3}$, where $V=(0,1)\times(0,1)$.

Subsequently, there are two major learning tasks: (1) fit the local coordinate
charts $\phi_{p}(\vec{v})$; (2) ensure the consistency between patches. The
first task is achieved by approximating each coordinate chart with a multi-
layer-perceptron (MLP) function (Murtagh, 1991) and optimizing the learnable
parameters with the Sinkhorn regularized distance (Cuturi, 2013):

$\mathcal{L}_{p}=\min_{P_{ij}}\>\sum_{i,j\leq
N_{p}}P_{ij}\|\phi(\vec{v}_{j};\mathbf{W}_{p})-\vec{x}_{j}\|^{2}+\chi^{-1}\sum_{i,j\leq
N_{p}}P_{ij}\log P_{ij}$ (3)

where $\mathcal{L}_{p}$ is the training loss function and $\mathbf{W}_{p}$ is
the set of neural network parameters for the $p$th coordinate chart, such that
$\phi_{p}(\vec{v})\approx\phi(\vec{v};\mathbf{W}_{p})$. $N_{p}$ is the number
of sampling points used to calculate the Sinkhorn distance, which is less than
or equal to the total available points within the same coordinate chart.
$\vec{v}_{j}\in V$ indicates the input samples following the Poisson disk
distribution for the reconstructed set. $P_{ij}$ is an $n\times n$ bi-
stochastic matrix. $\chi$ is a regularization parameter such that
$\mathcal{L}_{p}$ approximates the optimal transport distance (Rubner et al.,
1997) as $\chi\rightarrow\infty$.

The second task is performed in two steps: we first find the optimal indices
permutation policy of each patch $p$ via the following objective derived from
(3):

$\min_{\mathbf{W}_{p}}\>\inf_{\pi_{p}}\sum_{i\leq
N_{p}}\|\phi(\vec{v}_{i};\mathbf{W}_{p})-x_{\pi_{p}(i)}\|^{2}$ (4)

where $\pi_{p}$ is a permutation policy, assigning indices of points in
$\mathcal{X}_{p}$ to indices of parametric positions in $V_{p}$. We then
minimize the divergence between embedding functions $\phi_{p},\phi_{q}$ for
all pairs of overlapping patches:

$\min_{\mathbf{W}_{p},\mathbf{W}_{q}}\>\inf_{\pi_{p},\pi_{q}}\>\>\sum_{i\in
T_{pq}}\|\phi(\vec{v}_{i};\mathbf{W}_{p})-\phi(\vec{v}_{\pi^{-1}_{q}(\pi_{p}(i))};\mathbf{W}_{q})\|^{2}$
(5)

where
$T_{pq}=\\{i\,|\,\vec{x}_{\pi_{p}(i)}\in\mathcal{X}_{p}\cap\mathcal{X}_{q}\\}$
indicates the set of indices of parametric points in chart $p$ included within
the intersection of chart $p$ and $q$.

####  2.2.2 Local refinement of yield surface representation

The initial (nearly uniform) sampling of the plane stress conditions provides
an overall representation of the yield surface that may be too coarse. Because
of the substantial computational cost of DDD simulations, it is preferable to
locally refine the yield surface representation in regions with sharp edges or
corners than a uniform refinement. Similar to mesh refinement (Berger and
Oliger, 1984) in numerical simulations, the point cloud (for which DDD
predictions are collected) should be densified at the locations where it is
likely to provide the most information gain. In geometry reconstruction, we
believe curvature is a reasonable choice of quantitative criterion for such
purpose (Tang and Feng, 2018).

There are multiple curvature measures in the application of geometric
learning. In order to efficiently establish a criterion with inequality, we
adopted one of the scalar curvature measures: Gaussian curvature (Kreyszig,
2013), denoted as $\mathcal{K}$ (for mathematical formulation of Gaussian
curvature, please refer to (Xiao and Sun, 2022)). Gaussian curvature can serve
as a smoothness measure for digital images and geometric objects (Tang et al.,
2023), where a large Gaussian curvature usually indicates a relatively abrupt
change in the geometric shape (e.g. non-smoothness). Here we selected regions
for enrichment as locations where $\mathcal{K}>0.07\mathrm{~{}MPa}^{-2}$. As
Figure 5 shows, such locations group together around the “corners” of the
yield surface. For each of the corners, we then defined a bounding sphere
$\mathbf{\mathcal{B}}^{(e)}$ with minimum radius $R^{(e)}$ that includes all
the selected data points (whose $\mathcal{K}>0.07\mathrm{~{}MPa}^{-2}$), where
the superscript $e$ labels the enrichment regions (corners).

To refine the yield surface inside the bounding sphere
$\mathbf{\mathcal{B}}^{(e)}$, we first projected the bounding sphere onto the
unit sphere $\tilde{S}^{2}$ centered at the origin. The projection is defined
as $\mathcal{P}_{\tilde{S}^{2}}(\mathbf{\mathcal{B}}^{(e)})$. We used Poisson
disk sampling (Bowers et al., 2010; Williams, 2022) to generate a random but
relatively uniform distribution of points over the entire unit sphere
$\tilde{S}^{2}$ and extracted sampling points that lie within
$\mathcal{P}_{\tilde{S}^{2}}(\mathbf{\mathcal{B}}^{(e)})$. Each one of these
sampling points on the unit sphere represents a stress orientation
$\hat{\sigma}$ for which a DDD simulation is performed under constant stress
rate loading to determine the yield point.

To determine how many sampling points are needed in the refined regions, we
monitor the convergence of the local Gaussian curvature predicted by the
geometric prior at a given stress orientation,
$\hat{\sigma}=\left(-\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},0\right)$. Figure
6 plots the predicted Gaussian curvature at this point when 33, 67 and 97
additional data points are combined with the initial 92 data points to
construct the yield surface manifold. We observe a gradual increase of the
predicted local Gaussian curvature and an indication of convergence. Due to
the limitations of computational resources, we limit the total number of
sampling points on the unit sphere to 189, with each sampling point
corresponding to a DDD simulation at a different stress orientation.

Figure 6: Gaussian curvature of the learned yield surface manifold for
config-0 at the stress orientation
$\hat{\sigma}=\left(-\frac{\sqrt{2}}{2},-\frac{\sqrt{2}}{2},0\right)$ as a
function of the total number of stress orientations used to construct the
yield surface manifold.

The geometric learning process is then applied to the enriched dataset to
produce a relatively more accurate representation of the DDD yield surface.
This iterative procedure is illustrated in Figure 7. The original and refined
constructions of the yield surface for config-0 are shown in Figure 7(b) and
(d), respectively. Limited by the computing power and the high computational
cost of DDD simulations, the data enrichment step cannot be continued beyond
the first iteration, and we cannot obtain thousands of data points, a typical
number training data points for surface reconstruction (Williams et al.,
2019). In this work, we perform one step of data enrichment and our final
yield surface is trained on 189 data points from DDD simulations.

Figure 7: The geometric refinement of yield surface for config-0. (a) Yield
surface is generated by interpolated triangulation from the initial 92 data
points of yielding stress. (b) Based on (a), the yield surface is constructed
by the fitting of geometric prior. (c) Based on (b), more sampling points are
selected where DDD simulations are performed for data enrichment, resulting in
189 data points on the yield surface. (d) Based on (c) the yield surface is
constructed by the fitting of geometric prior.

##  3 Results

###  3.1 Yield surfaces from plane stress loading

Figure 8: Along three different viewing angles, yield surfaces of (a-c)
config-0, [100] config-1 and [100] config-2, and (d-f) config-0, [110]
config-1 and [110] config-2 under plane stress loading. The green straight
line in (d-f) denotes the stress orientation of uniaxial loading along [110]
crystal orientation, i.e., $\hat{\sigma}=\left(1/2,1/2,\sqrt{2}/2\right)$.

Figure 8 shows the yield surface constructed from the DDD simulations under
plane stress loading from the initial dislocation configuration (config-0),
together with how the yield surface evolves by plastic deformation along [100]
and [110] directions, respectively. Figure 8 (a-c) show the yield surfaces
from different viewing angles for config-0, [100] config-1, [100] config-2. It
can be seen that with increasing plastic strain along [100] crystal
orientation, the yield surface appears to expand isotropically in all stress
orientations. This indicates that plastic deformation along [100] direction
results in isotropic hardening. Figure 8 (d-f) show the yield surfaces from
different viewing angles for config-0, [110] config-1, [110] config-2. It can
be seen that with increasing plastic strain along [110] crystal orientation,
the expansion of the yield surface is not isotropic in all stress directions.
Specifically, the yield surface appears to be not expanding appreciably along
the stress orientations neighboring that of the [110] uniaxial loading itself.
(The stress orientation corresponding to [110] uniaxial loading is indicated
by the green straight lines.) This is consistent with the low strain-hardening
rate in the stress-strain curve of [110] loading shown in Figure 3(b).
However, in stress orientations “perpendicular” to that of the [110] uniaxial
loading in the space of
$\left(\sigma_{xx},\,\sigma_{yy},\sqrt{2}\sigma_{xy}\right)$, the expansion of
yield surface (i.e. strain hardening) is significant. This indicates that
plastic deformation along [110] direction results in latent hardening,
representing the phenomenon that yield stress and hardening rate are higher in
previously unactivated slip systems than in the previously activated primary
slip systems (Edwards et al., 1953; Nakada and Keh, 1969; Kocks and Mecking,
2003). Latent hardening has been reported in the experiments of FCC single
crystals where the initial loading orientations are single-slip orientations,
which lie in the central region of the stereographic triangle (Kocks and
Brown, 1966; Jackson and Basinski, 1967; Wessels and Jackson, 1969). In
comparison, the $[110]$ loading direction is at a corner of the stereographic
triangle and activates multiple slip systems.

Figure 9: Hardening ratios of (a)-(c) [100] config-1 and [100] config-2, and
(d)-(f) [110] config-1 and [110] config-2 in different perspectives. The
innermost unit sphere is plotted as a reference. The stress orientation
$\hat{\sigma}=\left({1}/{2},{1}/{2},{\sqrt{2}}/{2}\right)$ of the [110]
uniaxial loading is shown as a green line.

To better illustrate the strain hardening behaviors, in Figure 9 we plot the
hardening ratio, which is defined as the yield surfaces at different plastic
strains over the initial yield surface (i.e. config-0). (The innermost shape
is a unit sphere plotted as a reference.) For each configuration, the
orientation-dependent hardening ratio data is fitted to an ellipsoid for
smoothing. Figure 9 (a-c) shows that for plastic deformation along [100]
crystal orientation, the hardening ratio is nearly isotropic in all stress
orientations. In contrast, Figure 9 (d-f) shows that for plastic deformation
along [110] crystal orientation, the hardening ratio is close to 1 (i.e. no
hardening) along the stress orientation corresponding to the [110] loading
orientation itself, and is greater than 1 and nearly isotropic in the stress
orientations perpendicular (in the $\hat{\sigma}$ space) to that of [110]
uniaxial loading.

Figure 10: Cut-plane views of the yield surfaces of (a-c) config-0, [100]
config-1 and [100] config-2, and (d-f) config-0, [110] config-1 and [110]
config-2 under plane stress loading. The green and black arrows denote the
stress orientations of [100] and [110] uniaxial loadings, respectively.

For a more detailed examination of how the yield surface evolves with plastic
strain, Figure 10 plots the cross-section views of the yield surfaces for
different cut-planes in the space of stress orientations. When viewed on the
cut-plane $\left(-\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2},0\right)$, the yield
surface for the original configuration (config-0) has approximately a square
shape. Figure 10(a) shows that with [100] uniaxial loading, the cross-section
view of the yield surface expands while keeping the square shape. In contrast,
Figure 10(d) shows that with [110] uniaxial loading, the cross-section view of
the yield surface expands asymmetrically into a rectangular shape. Figure
10(b-c) shows the cross-section views of yield surfaces under [100] uniaxial
loading on cut-planes perpendicular to that in (a); the yield surface expands
nearly isotropically in all directions. In comparison, Figure 10(e-f) shows
corresponding views of yield surfaces under [110] uniaxial loading, where the
expansion appears less in magnitude and is asymmetric.

###  3.2 Yield loci for uniaxial tension

Experimentally uniaxial loading tests are much easier to carry out than the
general plane-stress loading (which may be performed on a tube with combined
tension and torsion (Taylor and Quinney, 1931)). Therefore we examine the
predicted yield conditions under uniaxial tension. The uniaxial tensile
loading conditions correspond to a sub-set of stress orientations within the
plane-stress conditions and are shown as the red circle on the unit sphere in
Figure 2. Note that the center of the circle is not at the origin in the
$\hat{\sigma}$ space. The central inversion of this circle with respect to the
origin produces another circle, which corresponds to uniaxial compressive
loading conditions.

Figure 11: The 2D view of the yield loci along the
$\hat{\sigma}=\left(\sqrt{2}/2,\sqrt{2}/2,0\right)$ direction for uniaxial
tensile loading applied to (a) config-0, [100] config-1 and [100] config-2,
and (b) config-0, [110] config-1 and [110] config-2.

Figure 11 plots the yield loci for the uniaxial tensile loading conditions, in
terms of $\sqrt{2}\sigma_{xy}$ v.s. $(\sigma_{xx}-\sigma_{yy})/\sqrt{2}$. This
figure corresponds to viewing the yield loci corresponding to uniaxial tension
on the yield surface along the
$\hat{\sigma}=\left(\sqrt{2}/2,\sqrt{2}/2,0\right)$ direction. Figure 11(a)
shows the evolution of yield loci from config-0 during plastic deformation
along the [100] direction. The yield loci appears to expand isotropically
while maintaining an oval shape. We note that the Von Mises (i.e. $J_{2}$)
yield criterion (widely used for polycrystals) would predict the yield loci
(for uniaxial tension) in the shape of a circle when plotted using these axes.
Hence the deviation of the yield loci from a circle is a manifestation of the
plastic anisotropy of single crystals. Figure 11(b) shows the evolution of
yield loci from config-0 during plastic deformation along the [110] direction.
There is little expansion of the yield loci around the [110] loading direction
itself. But a substantial expansion is observed in the opposite end of the
loop, corresponding to the $[\bar{1}10]$ loading orientation. The predictions
of yield loci evolution from plastic deformation in the $[110]$ loading
direction, especially the pronounced latent hardening in the $[\bar{1}10]$
direction, provide new insight to the strain hardening behavior of FCC
crystals. Existing experimental studies on latent hardening in FCC crystals
were usually conducted where the initial uniaxial loading is along a direction
that favors single slip, i.e. only a single slip system is activated
(Franciosi et al., 1980; Franciosi, 1985). Subsequently, a “daughter” sample
is cut from the “parent” sample and then loaded in a different orientation.
The goal was to let the primary slip system in the “daughter” sample interact
strongly with the primary slip system previously activated in the “parent”
sample, to cause latent hardening. We hope our prediction of the latent
hardening behavior for initial loading along the $[110]$ direction will
motivate new experiments for its verification.

##  4 Discussion

###  4.1 Strain hardening mechanisms

The flow stress $\tau$ of a metal is known to correlate with the dislocation
density $\rho$. In particular, the Taylor relation (Taylor, 1934),
$\tau=\alpha\mu b\sqrt{\rho}$, is well supported by experimental data on
polycrystals (Wang et al., 2021), where $\mu$ is the shear modulus, $b$ is the
Burgers vector magnitude, and $\alpha$ is a dimensionless constant. Here the
total dislocation density for config-0 is $\sim 1.2\times
10^{12}\mathrm{~{}m}^{-2}$. After uniaxial deformation along $[100]$ crystal
orientation, the toal dislocation density becomes $\sim 1.3\times
10^{12}\mathrm{~{}m}^{-2}$ and $\sim 2.0\times 10^{12}\mathrm{~{}m}^{-2}$, for
config-1 and config-2, respectively. After uniaxial deformation along $[110]$
crystal orientation, the total dislocation density becomes $\sim 1.3\times
10^{12}\mathrm{~{}m}^{-2}$ and $\sim 1.6\times 10^{12}\mathrm{~{}m}^{-2}$, for
config-1 and config-2, respectively. Although in both deformation
orientations, the strain hardening is accompanied by an increase of total
dislocation density, this does not explain the difference between the
isotropic hardening and latent hardening behaviors, for which an examination
of dislocation density on individual slip systems is necessary.

Table 1: Schmid factors and dislocation densities (in $10^{11}\mathrm{~{}m}^{-2}$) for the 12 slip systems in various dislocation configurations (Config-0, Config-1 and Config-2) extracted from the $[100]$ uniaxial loading. The active slip systems are highlighted in gray. Slip system index: $i$, Slip plane: n, Burgers vector: b, Schmid-Boas notation: SB, Schmid factors: SF $i$ |  n |  b |  SB |  SF |  Config-0 |  Config-1 |  Config-2
---|---|---|---|---|---|---|---
1 |  $(\overline{1}11)$ |  $\frac{1}{2}[0\overline{1}1]$ |  A2 |  0 |  0.6945 |  0.6354 |  1.0221
2 |  $(\overline{1}11)$ |  $\frac{1}{2}[101]$ |  A3 |  0.4082 |  0.8599 |  0.8833 |  1.4480
3 |  $(\overline{1}11)$ |  $\frac{1}{2}[110]$ |  A6 |  0.4082 |  0.8882 |  1.0106 |  1.7339
4 |  $(111)$ |  $\frac{1}{2}[0\overline{1}1]$ |  B2 |  0 |  0.7711 |  0.5372 |  1.0167
5 |  $(111)$ |  $\frac{1}{2}[\overline{1}01]$ |  B4 |  0.4082 |  0.6238 |  0.8580 |  1.3295
6 |  $(111)$ |  $\frac{1}{2}[\overline{1}10]$ |  B5 |  0.4082 |  0.8691 |  1.2066 |  1.9012
7 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[011]$ |  C1 |  0 |  0.8838 |  0.8195 |  1.0448
8 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[101]$ |  C3 |  0.4082 |  0.7518 |  0.8251 |  1.2895
9 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[\overline{1}10]$ |  C5 |  0.4082 |  0.8541 |  0.7792 |  1.1405
10 |  $(1\overline{1}1)$ |  $\frac{1}{2}[011]$ |  D1 |  0 |  0.7214 |  0.5979 |  0.8278
11 |  $(1\overline{1}1)$ |  $\frac{1}{2}[\overline{1}01]$ |  D4 |  0.4082 |  0.7460 |  0.7744 |  1.4479
12 |  $(1\overline{1}1)$ |  $\frac{1}{2}[110]$ |  D6 |  0.4082 |  0.6779 |  0.6584 |  0.9163

Table 1 shows the dislocation densities on individual slip systems during
uniaxial deformation along $[100]$ direction. In FCC metals, there are 4 slip
planes, each having 3 slip directions, with a total of 12 slip systems. For
uniaxial loading along $[100]$ direction, 8 slip systems (2 on every slip
plane) are activated with the same Schmid factor. (The remaining 4 slip
systems have zero Schmid factor.) Table 1 shows that the dislocation densities
on the active slip systems increase with strain for both config-1 (strained to
0.07%, except the slight decrease for C5 and D6 slip systems) and config-2
(strained to 0.22%). On the inactive slip systems (A2, B2, C1, D1), the
dislocation densities in config-1 is actually lower than those in config-0;
however, the dislocation densities on these slip systems in config-2 exceed
those in config-0. This phenomenon is termed “slip-free multiplication”
(Akhondzadeh et al., 2021), because dislocations multiply on these slip
systems with zero Schmid-factor and zero shear strain rate. Slip-free
multiplication has been linked to coplanar interactions between active slip
systems on the same slip plane (e.g. dislocations on inactive slip system A2
are produced because both A3 and A6 are active). We conclude that the
isotropic hardening produced by plastic deformation along $[100]$ is mainly
because all 4 slip planes have 2 slip systems activated resulting in increased
dislocation density on all slip planes. The slip-free multiplication is deemed
not essential for the isotropy of strain hardening, because even though
dislocation densities on inactive slip systems (A2, B2, C1, D1) are lower in
config-1 than in config-0, the yield surface of config-1 appears to have
expanded isotropically from that of config-0 nonetheless. In other words,
having dislocations multiply on 2 out of 3 of slip systems on all 4 slip
planes appears to be sufficient for isotropic hardening. Existing experiments
also suggest that the difference of hardening behaviors in FCC single crystal
metals depends on the activated slip planes instead of the slip directions or
slip systems (Kocks and Brown, 1966).

Table 2: Schmid factors and dislocation densities (in $10^{11}\mathrm{~{}m}^{-2}$) for the 12 slip systems in various dislocation configurations (Config-0, Config-1 and Config-2) extracted from the $[110]$ uniaxial loading. The active slip systems are highlighted in gray. Slip system index: $i$, Slip plane: n, Burgers vector: b, Schmid-Boas notation: SB, Schmid factors: SF $i$ |  n |  b |  SB |  SF |  Config-0 |  Config-1 |  Config-2
---|---|---|---|---|---|---|---
1 |  $(\overline{1}11)$ |  $\frac{1}{2}[0\overline{1}1]$ |  A2 |  0 |  0.6945 |  0.5536 |  0.4540
2 |  $(\overline{1}11)$ |  $\frac{1}{2}[101]$ |  A3 |  0 |  0.8599 |  0.6544 |  0.5682
3 |  $(\overline{1}11)$ |  $\frac{1}{2}[110]$ |  A6 |  0 |  0.8882 |  0.7212 |  0.5563
4 |  $(111)$ |  $\frac{1}{2}[0\overline{1}1]$ |  B2 |  0.4082 |  0.7711 |  1.0985 |  1.4954
5 |  $(111)$ |  $\frac{1}{2}[\overline{1}01]$ |  B4 |  0.4082 |  0.6238 |  1.0217 |  1.3904
6 |  $(111)$ |  $\frac{1}{2}[\overline{1}10]$ |  B5 |  0 |  0.8691 |  0.8092 |  1.1442
7 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[011]$ |  C1 |  0.4082 |  0.8838 |  1.5671 |  2.4345
8 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[101]$ |  C3 |  0.4082 |  0.7518 |  0.9866 |  1.6820
9 |  $(\overline{1}\overline{1}1)$ |  $\frac{1}{2}[\overline{1}10]$ |  C5 |  0 |  0.8541 |  0.5988 |  0.9298
10 |  $(1\overline{1}1)$ |  $\frac{1}{2}[011]$ |  D1 |  0 |  0.7214 |  0.5783 |  0.4366
11 |  $(1\overline{1}1)$ |  $\frac{1}{2}[\overline{1}01]$ |  D4 |  0 |  0.7460 |  0.5217 |  0.3388
12 |  $(1\overline{1}1)$ |  $\frac{1}{2}[110]$ |  D6 |  0 |  0.6779 |  0.5183 |  0.4730

Table 2 shows the dislocation densities on individual slip systems during
uniaxial deformation along $[110]$ direction, in which only 4 slip systems
(B2, B4, C1, C3) are activated with the same Schmid factor. The remaining 8
slip systems have zero Schmid factor. All slip systems on slip planes A and D
are inactive. Table 1 shows that the dislocation densities on the active slip
systems increase with strain for both config-1 (strained to 0.07%) and
config-2 (strained to 0.22%). On the inactive slip systems B5 and C5, the
dislocation densities first decrease (in config-1) and eventually increase (in
config-2), due to slip-free multiplication (because their coplanar slip
systems B2, B4 and C1, C3 are active). On the other hand, all slip systems on
slip planes A and D are inactive; consequently slip-free multiplication does
not occur for these slip systems and their dislocation densities only decrease
with strain. We conclude that the latent hardening produced by plastic
deformation along $[110]$ crystal orientation is due to the highly uneven
dislocation densities on different slip planes. For example, if a sample first
deformed along the $[110]$ direction is unloaded and then reloaded along the
$[\bar{1}10]$ direction, then the slip systems A2, A3, D1, D4 will be
activated; dislocations on these slip systems will interact strongly with
existing ones on B2, B4, C1, C3 through glissile and collinear junctions
(although no Lomer junctions).

###  4.2 Next steps in DDD simulations

In this work, we have demonstrated that DDD simulations can be used to
construct the yield surface of single-crystal Cu under plane-stress conditions
and how the yield surface evolves under plastic deformation along uniaxial
tensile deformation along $[100]$ and $[110]$ directions. Due to computational
limits on the simulation time scale, the yield surface calculations are
calculated under a stress rate of $10^{13}\,{\rm Pa}\cdot{\rm s}^{-1}$, which
is equivalent to the strain rate of $10^{3}\,{\rm s}^{-1}$ for the uniaxial
deformation. However, what we usually have in mind when discussing yield
surface or strain-hardening behaviors as well as most of the existing
experimental data pertain to much lower, quasi-static strain rates, ranging
from $10^{-4}\,{\rm s}^{-1}$ to $10^{-1}\,{\rm s}^{-1}$ (Meyers, 1994; Meyers
and Chawla, 2008). Since strain-rate effect on the plastic deformation exists
in most of the materials (Meyers, 1994), including polymers (Farrokh and Khan,
2010), metals (Butcher and Karnes, 1966; Meyers, 2012) and composites
(Sierakowski, 1997), it is expected that the yield surface of single crystal
copper is strain-rate dependent. At present, lowering the strain-rate of DDD
simulations to the level of, say $10^{-3}\,{\rm s}^{-1}$, by brute force
appears out of reach. Fortunately, existing experimental data (Edington, 1969)
suggests that the flow stress of single-crystal Cu appears to be strain-rate
independent in a wide range of strain-rates from $10^{-4}\,{\rm s}^{-1}$ to
$10^{1}\,{\rm s}^{-1}$. DDD simulations of single-crystal Cu at the strain
rate of $10^{2}\,{\rm s}^{-1}$ have been reported earlier (Sills et al.,
2018). Hence DDD simulations at the strain rate of $10^{1}\,{\rm s}^{-1}$
appear feasible in the near future, with a bigger commitment of computational
resources. A challenge that one may need to face is that the size of the DDD
simulation cell is likely to increase in order to reach convergence when the
strain-rate is lowered, which would further increase the computational cost.

In addition, more DDD simulations would be needed if we wish to go beyond the
confines of plane-stress loading conditions and construct the yield surface in
the general 6-dimensional stress space. It also appears necessary to repeat
DDD simulations under the same stress conditions for different initial
dislocation configurations to improve statistics.

###  4.3 Next steps in geometric prior method

We have demonstrated that the geometric prior method can be used to construct
the yield surface as a manifold in stress space based on data obtained from
DDD simulations. So far, the manifolds corresponding to each dislocation
configuration (e.g. [100] config-1, [100] config-2, etc.) are constructed
independently of each other, even though they are related to each other
through plastic deformation starting from a common configuration (config-0).
As a natural next step, it would be of interest for the geometric prior method
to learn not only the yield surfaces separately, but its evolution with
plastic deformation. Such a model is what is ultimately needed in a continuum
constitutive model of crystal plasticity.

Intuitively, a starting point for modeling the evolution of the yield surface
is to generate enough snapshots of the yield surface at different plastic
deformation levels. However, the number of yield surfaces that can be obtained
from DDD simulations is limited due to the high computational cost. Hence an
interpolation scheme is needed between sparse data points along the strain
axis. Based on differential manifold theory (Kreyszig, 2013), the
interpolation between two yield surfaces should be rigorously defined by a
transformation mapping from one manifold to another one. To derive an
interpolation policy from a manifold transformation, it is necessary that the
two underlined surfaces are isomorphic, i.e., it is required to establish the
one-to-one correspondence between two points on these two surfaces. This
appears to be a valid assumption at the strain level accessible to DDD
simulations in the near future.

For the geometric prior method used in this study, we use a few patches to
describe the local features of a yield surface. The patches can overlap and
then provide a complete description of the entire yield surface. To achieve an
accurate description of yield surface in this work, different numbers of
patches are used for various yield surfaces and the patch number varies from
30 to 37. If two yield surfaces are represented by different number of
patches, it is difficult to build point-to-point mapping between them for
interpolation. Therefore, in the next step, the yield surfaces from two
configurations along the same deformation trajectory should be described by
the same number of patches. One promising approach is to put effort to the
parameterization of the motions of individual patches (translation, scaling,
rotation, …). This would allow us to construct a plasticity evolution law by
establishing functional relationships between the patch motion parameters and
the cumulative plastic strain (Xiao and Sun, 2022). Efforts are needed to
maintain consistency between patches in the overlapping regions in this new
scheme (Xiao and Sun, 2022).

To extend the yield surface to the full stress space (six-dimensional), one
promising approach would be turning back to the implicit yield function
representation but fitting it with neural kernel method (Xiong et al., 2023),
so that the yield surface is constituted by a smooth basis extracted from
kernel regression. This generally produce a smoothed yield surface with good
generalizability, i.e., the yield surface will not deform to complex shapes
locally due to insufficient training data, but some accuracy will be lost at
locations with real geometric complexities. Essentially, in more-than-three-
dimensional spaces, the hyper-surface reconstruction is still a challenging
task, and the trade-off between local accuracy and global generalizability is
more difficult to make than the cases in three-dimensional spaces.

##  5 Conclusions

In this work, we demonstrate a framework to construct the yield surface of
single crystals using DDD simulations and the geometric prior method. DDD
simulations under constant stress-rate conditions have been performed to
identify the yielding conditions and the data are used by the geometric prior
method to construct a cross-section of the yield surface as a manifold in the
3-dimensional sub-space of plane-stress. An iterative workflow is adopted,
where the geometric prior method identifies regions of interest (based on
local curvature) where further sampling of the stress conditions by DDD
simulations is performed.

We found that the yield surface evolves differently by plastic deformation
along the $[100]$ and $[110]$ loading directions, respectively. Isotropic
hardening is observed for $[100]$ deformation, in which the yield surface
expands with strain by nearly the same ratio in all stress orientations. This
is traced to dislocation multiplication on all four slip planes during $[100]$
deformation. In contrast, latent hardening is observed for $[110]$
deformation, in which the yield surface does not expand much at all in the
vicinity of the stress orientations corresponding to $[110]$ tension, but
expands significantly in other orientations. This is traced to dislocation
multiplication on only two out of the four slip planes during $[110]$
deformation.

##  6 Acknowledgements

This work was supported by the National Science Foundation, under Award Number
DMREF 2118522 (W.J. and W.C.). W.C.S. and X.M. is supported by the UPS
Foundation Visiting Professorship from Stanford University, with additional
support from the National Science Foundation under grant contracts
CMMI-1846875 and the Dynamic Materials and Interactions Program from the Air
Force Office of Scientific Research under grant contracts FA9550-21-1-0391 and
FA9550-21-1-0027.

##  Appendix A Constant stress rate vs constant strain rate loading

The comparison between the stress-strain curves for uniaxial tensile loading
with constant stress rate of $10^{13}\,{\rm Pa}\cdot{\rm s}^{-1}$ and constant
strain rate of $10^{3}\,{\rm s}^{-1}$ along the $[100]$ orientation is shown
in Figure 1.

Figure 1: The stress-strain curves for uniaxial tensile loadings with constant
stress rate of $10^{13}\,{\rm Pa}\cdot{\rm s}^{-1}$ and constant strain rate
of $10^{3}\,{\rm s}^{-1}$ along the $[100]$ orientation, respectively.

##  Appendix B DDD simulation parameters

The parameters of our DDD simulations are summarized in Table 1.

Table 1: Summary of DDD simulation parameters for single crystal copper. Property |  Parameter |  Value
---|---|---
Shear modulus |  $\mu$ |  54.6 GPa
Poission’s ratio |  $\nu$ |  0.324
Burgers vector magnitude |  $b$ |  0.255 nm
Drag coefficient |  $B$ |  15.6 $\mathrm{\mu~{}Pa}\cdot\mathrm{s}$
Core radius |  a |  6 $b$
Relative tolerance |  $f_{\rm tol}$ |  0.1 $b$
Absolute tolerance |  $r_{\rm tol}$ |  10 $b$
Threshold tolerance |  $r_{\rm th}$ |  1 $b$
Subcycling group 1 radius |  $r_{g1}$ |  0
Subcycling group 2 radius |  $r_{g2}$ |  100 $b$
Subcycling group 3 radius |  $r_{g3}$ |  600 $b$
Subcycling group 4 radius |  $r_{g4}$ |  1600 $b$
Max segment length |  $l_{\rm max}$ |  2000 $b$
Min segment length |  $l_{\rm min}$ |  $\sqrt{\frac{4}{\sqrt{3}}A_{\min}}$
Max area between segments |  $A_{\rm max}$ |  $\frac{1}{2}\left(4A_{\min}+\frac{\sqrt{3}}{4}l_{\max}^{2}\right)$
Min area between segments |  $A_{\rm min}$ |  $\min\left(2r_{\mathrm{tol}}l_{\max},\frac{\sqrt{3}}{4}l_{\min}^{2}\right)$
Collision radius |  $r_{\rm col}$ |  10 $b$
Simulation cell size |  $L$ |  15 $\mu$m
Strain rate for uniaxial loading |  $\dot{\varepsilon}$ |  $10^{3}\,{\rm s}^{-1}$
Stress rate for plane stress loading |  $\dot{\left\|\sigma\right\|}$ |  $10^{13}\,{\rm Pa}\cdot{\rm s}^{-1}$

## References

  * Akhondzadeh et al. (2021) Akhondzadeh, S., Bertin, N., Sills, R.B., Cai, W., 2021\.  Slip-free multiplication and complexity of dislocation networks in fcc metals.  Mater. Theory 5, 1–24.
  * Akhondzadeh et al. (2020) Akhondzadeh, S., Sills, R.B., Bertin, N., Cai, W., 2020\.  Dislocation density-based plasticity model from massive discrete dislocation dynamics database.  J. Mech. Phys. Solids 145, 104152\.
  * Amodeo and Ghoniem (1990) Amodeo, R.J., Ghoniem, N.M., 1990\.  Dislocation dynamics. i. a proposed methodology for deformation micromechanics.  Phys. Rev. B 41, 6958\.
  * Arsenlis et al. (2007) Arsenlis, A., Cai, W., Tang, M., Rhee, M., Oppelstrup, T., Hommes, G., Pierce, T.G., Bulatov, V.V., 2007\.  Enabling strain hardening simulations with dislocation dynamics.  Modelling Simul. Mater. Sci. Eng. 15, 553.
  * Bauschinger (1886) Bauschinger, J., 1886.  On the change of the position of the elastic limit of iron and steel under cyclic variations of stress.  Mitt. Mech.-Tech. Lab., Munich 13\.
  * Berger and Oliger (1984) Berger, M.J., Oliger, J., 1984\.  Adaptive mesh refinement for hyperbolic partial differential equations.  J. Comput. Phys. 53, 484–512.
  * Bertin et al. (2019) Bertin, N., Aubry, S., Arsenlis, A., Cai, W., 2019\.  Gpu-accelerated dislocation dynamics using subcycling time-integration.  Modelling Simul. Mater. Sci. Eng. 27, 075014.
  * Bowers et al. (2010) Bowers, J., Wang, R., Wei, L.Y., Maletz, D., 2010.  Parallel poisson disk sampling with spectrum analysis on surfaces.  ACM Trans. Graph. (TOG) 29, 1–10.
  * Butcher and Karnes (1966) Butcher, B., Karnes, C., 1966\.  Strain-rate effects in metals.  J. Appl. Phys. 37, 402–411.
  * Coombs and Motlagh (2017) Coombs, W.M., Motlagh, Y.G., 2017\.  Nurbs plasticity: yield surface evolution and implicit stress integration for isotropic hardening.  Comput. Methods Appl. Mech. Eng. 324, 204–220.
  * Coombs and Motlagh (2018) Coombs, W.M., Motlagh, Y.G., 2018\.  Nurbs plasticity: non-associated plastic flow.  Comput. Methods Appl. Mech. Eng. 336, 419–443.
  * Coombs et al. (2016) Coombs, W.M., Petit, O.A., Motlagh, Y.G., 2016.  Nurbs plasticity: Yield surface representation and implicit stress integration for isotropic inelasticity.  Comput. Methods Appl. Mech. Eng. 304, 342–358.
  * Cuturi (2013) Cuturi, M., 2013.  Sinkhorn distances: Lightspeed computation of optimal transport.  Adv. Neural Inf. Process. 26.
  * Dahl (2023) Dahl, V.A., 2023.  Icosphere 0.1.3.  https://pypi.org/project/icosphere/.
  * Devincre and Kubin (1997) Devincre, B., Kubin, L.P., 1997\.  Mesoscopic simulations of dislocations and plasticity.  Mater. Sci. Eng. A 234, 8–14.
  * Edington (1969) Edington, J., 1969.  The influence of strain rate on the mechanical properties and dislocation substructure in deformed copper single crystals.  Philos. Mag. 19, 1189–1206.
  * Edwards et al. (1953) Edwards, E., PARKER, E., WASHBURN, J., 1953.  Some observations on the work hardening of metals.  Trans. Metall. Soc. AIME 197, 1525–1529.
  * Farrokh and Khan (2010) Farrokh, B., Khan, A.S., 2010\.  A strain rate dependent yield criterion for isotropic polymers: low to high rates of loading.  Eur. J. Mech. A/Solids 29, 274–282.
  * Franciosi (1985) Franciosi, P., 1985.  The concepts of latent hardening and strain hardening in metallic single crystals.  Acta Metall. 33, 1601–1612.
  * Franciosi et al. (1980) Franciosi, P., Berveiller, M., Zaoui, A., 1980.  Latent hardening in copper and aluminium single crystals.  Acta Metall. 28, 273–283.
  * Hill (1998) Hill, R., 1998.  The mathematical theory of plasticity. volume 11.  Oxford university press.
  * Hirsch et al. (1956) Hirsch, P., Horne, R., Whelan, M., 1956.  LXVIII. Direct observations of the arrangement and motion of dislocations in aluminium.  Philos. Mag. 1, 677–684.
  * Hirth and Lothe (1982) Hirth, J., Lothe, J., 1982.  Theory of Dislocations.  Wiley.
  * Jackson and Basinski (1967) Jackson, P., Basinski, Z., 1967\.  Latent hardening and the flow stress in copper single crystals.  Can. J. Phys. 45, 707–735.
  * Kocks and Brown (1966) Kocks, U., Brown, T., 1966.  Latent hardening in aluminum.  Acta Metall. 14, 87–98.
  * Kocks and Mecking (2003) Kocks, U., Mecking, H., 2003\.  Physics and phenomenology of strain hardening: the fcc case.  Prog. Mater. Sci. 48, 171–273.
  * Kreyszig (2013) Kreyszig, E., 2013.  Differential geometry.  Courier Corporation.
  * Meyers (2012) Meyers, M., 2012.  Shock waves and high-strain-rate phenomena in metals: concepts and applications.  Springer Science & Business Media.
  * Meyers (1994) Meyers, M.A., 1994.  Dynamic behavior of materials.  John wiley & sons.
  * Meyers and Chawla (2008) Meyers, M.A., Chawla, K.K., 2008\.  Mechanical behavior of materials.  Cambridge university press.
  * Murtagh (1991) Murtagh, F., 1991.  Multilayer perceptrons for classification and regression.  Neurocomputing 2, 183–197.
  * Nakada and Keh (1969) Nakada, Y., Keh, A., 1969.  Latent hardening in rock-salt type crystals.  Phys. Status Solidi B 32, 715–730.
  * Orowan (1934) Orowan, E., 1934.  Plasticity of crystals.  Zeit. Fur Phys. 89, 605–659.
  * ParaDiS (2023) ParaDiS, 2023.  Parallel dislocation simulator.  https://gitlab.com/opendis/ParaDiS-2.7.
  * Polanyi (1934) Polanyi, M., 1934.  Lattice distortion which originates plastic flow.  Zeit. Fur Phys. 89, 660–662.
  * Rubner et al. (1997) Rubner, Y., Guibas, L.J., Tomasi, C., 1997.  The earth mover’s distance, multi-dimensional scaling, and color-based image retrieval, in: Proceedings of the ARPA image understanding workshop, p. 668\.
  * Sierakowski (1997) Sierakowski, R.L., 1997.  Strain rate effects in composites.  Appl. Mech. Rev. 50, 741–761.
  * Sills et al. (2016) Sills, R.B., Aghaei, A., Cai, W., 2016.  Advanced time integration algorithms for dislocation dynamics simulations of work hardening.  Modelling Simul. Mater. Sci. Eng. 24, 045019.
  * Sills et al. (2018) Sills, R.B., Bertin, N., Aghaei, A., Cai, W., 2018\.  Dislocation networks and the microstructural origin of strain hardening.  Phys. Rev. Lett. 121, 085501\.
  * Tang et al. (2023) Tang, W., Lin, Z., Gong, Y., 2023.  Gc-net: An unsupervised network for gaussian curvature optimization on images.  J. Signal Process. Syst 95, 77–88.
  * Tang and Feng (2018) Tang, Y., Feng, J., 2018.  Multi-scale surface reconstruction based on a curvature-adaptive signed distance field.  Comput. Graph. 70, 28–38.
  * Taylor (1934) Taylor, G.I., 1934.  The mechanism of plastic deformation of crystals. part I.—theoretical.  Proc. R. Soc. Lond. A 145, 362–387.
  * Taylor and Quinney (1931) Taylor, G.I., Quinney, H., 1931\.  The plastic distortion of metals.  Phil. Trans. Roy. Soc. A 230, 323–362.
  * Vlassis and Sun (2021) Vlassis, N.N., Sun, W., 2021\.  Sobolev training of thermodynamic-informed neural networks for interpretable elasto-plasticity models with level set hardening.  Comput. Methods Appl. Mech. Eng. 377, 113695.
  * Vlassis and Sun (2022) Vlassis, N.N., Sun, W., 2022\.  Component-based machine learning paradigm for discovering rate-dependent and pressure-sensitive level-set plasticity models.  J. Appl. Mech. 89.
  * Wang et al. (2021) Wang, L., Lu, Z., Li, H., Zheng, Z., Zhu, G., Park, J.S., Zeng, X., Bieler, T.R., 2021.  Evaluating the taylor hardening model in polycrystalline ti using high energy x-ray diffraction microscopy.  Scr. Mater. 195, 113743\.
  * Wessels and Jackson (1969) Wessels, E., Jackson, P., 1969\.  Latent hardening in copper-aluminium alloys.  Acta Metall. 17, 241–248.
  * Williams (2022) Williams, F., 2022.  Point cloud utils.  https://www.github.com/fwilliams/point-cloud-utils.
  * Williams et al. (2019) Williams, F., Schneider, T., Silva, C., Zorin, D., Bruna, J., Panozzo, D., 2019\.  Deep geometric prior for surface reconstruction, in: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 10130–10139.
  * Xiao and Sun (2022) Xiao, M., Sun, W., 2022.  Geometric prior of multi-resolution yielding manifolds and the local closest point projection for nearly non-smooth plasticity.  Comput. Methods Appl. Mech. Eng. 400, 115469.
  * Xiong et al. (2023) Xiong, Z., Xiao, M., Vlassis, N., Sun, W., 2023\.  A neural kernel method for capturing multiscale high-dimensional micromorphic plasticity of materials with internal structures.  Comput. Methods Appl. Mech. Eng. 416, 116317.