# Interfaces in crystalline materials

## Aurélien Vattré

► To cite this version:

Aurélien Vattré. Interfaces in crystalline materials. Engineering Sciences [physics]. Sorbonne Paris Nord, 2023. tel-04149115

HAL Id: tel-04149115

https://hal.science/tel-04149115

Submitted on 3 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
HABILITATION À DIRIGER DES RECHERCHES

AURÉLIEN VATTRÉ

Maître de Recherche à l'ONERA

**Interfaces in crystalline materials**

Université Sorbonne Paris Nord
Habilitation soutenue le 10 mai 2023 devant le jury composé de :

Ioan Ionescu

Brigitte Bacroix

Stéphane Berbenni

Marc Fivel

Renald Brenner

Président

Rapporteure

Rapporteur

Rapporteur

Examinateur
Interfaces in crystalline materials

Interfaces such as grain boundaries in polycrystalline as well as and heterointerfaces in mul-
tiphase are ubiquitous in materials science and engineering with wide-ranging properties and
applications. Therefore, understanding the basics of interfaces is key in optimization of ceramics
for a wide range of applications including electrochemical energy conversion and storage, optical,
magnetic, and mechanical applications, thermal applications including thermal and environmen-
tal barrier coatings in automobile and aeronautical industries.

Far from being featureless dividing surfaces between neighboring crystals, elucidating fea-
tures of solid-solid interfaces is challenging and requires theoretical and numerical strategies to
describe the physical and mechanical characteristics of these internal interfaces. The first part
of this manuscript is concerned with interface-dominated microstructures emerging from poly-
morphic structural (diffusionless) phase transformations. Under high hydrostatic compression
and shock-wave conditions, the pressure-driven phase transitions and the formation of internal
diffuse interfaces in iron are captured by a thermodynamically consistent framework for combin-
ing nonlinear elastoplasticity and multivariant phase-field approach at large strains. The calcu-
lations investigate the crucial role played by the plastic deformation in the morphological and
microstructure evolution processes under high hydrostatic compression and shock-wave condi-
tions. The second section is intended to describe such imperfect interfaces at a finer scale, for
which the semicoherent interfaces are described by misfit dislocation networks that produce a
lattice-invariant deformation which disrupts the uniformity of the lattice correspondence across
the interfaces and thereby reduces coherency. For the past ten years, the constant effort has been
devoted to combining the closely related Frank-Bilby and O-lattice techniques with the Stroh sex-
tic formalism for the anisotropic elasticity theory of interfacial dislocation patterns. The structures
and energetics are quantified and used for rapid computational design of interfaces with tailored
misfit dislocation patterns, including the interface sink strength for radiation-induced point de-
fects and semicoherent interfaces.
# Contents

ii

1 Introduction 2

2 Crystalline interfaces during solid-solid phase transitions in iron 4

2.1 Motivation 4

2.2 A phase-field model coupled with finite elastoplasticity 6

2.2.1 Kinematics 6

2.2.2 Balance laws 7

2.2.3 The Clausius-Duhem inequality 7

2.2.4 Constitutive equations 9

2.2.5 Multiple reaction pathways and energy landscape 12

2.2.6 Computational framework 19

2.3 Pure hydrostatic compression 20

2.3.1 Material and model inputs 20

2.3.2 Analysis of the pressure-volume responses 21

2.3.3 Microstructure and variant selection 23

2.4 Shock wave propagation 25

2.4.1 The internal structure of shock waves 25

2.4.2 Effect of plasticity in shock-loaded iron 28

2.4.3 Residual stresses in the plastically-deformed microstructure 31

2.4.4 Dynamical instability in structural phase transitions 34

2.5 Limitations 34

3 Dislocation structures and energetics at heterophase interfaces 36

3.1 Motivation 37

3.2 Determining the Burgers vectors of interface dislocation arrays 38

3.2.1 Planar interfaces in linear elastic bicrystals 38

3.2.2 Volterra dislocations in the reference state 39

3.2.3 Crystallographic constraints on interface dislocations 40

3.2.4 Solution strategy 41

3.2.5 Elastic fields of interface dislocation arrays 42

3.2.6 Interface elastic strain energy 47

3.3 Symmetric example applications 47

3.3.1 Pure tilt grain boundary 48

3.3.2 Twist grain boundary 49

3.3.3 Pure misfit interface 50

3.4 Partitioning of elastic distortions at fcc/bcc interfaces 54

3.4.1 Mapping between states in the Nishiyama-Wassermann orientations 55

3.4.2 Far-field strains and rotations 57

3.4.3 Spurious fields from incorrect reference states 59

3.4.4 Orientations differing from the Nishiyama-Wassermann relations 60

3.4.5 Short-range elastic fields 61
<table>
  <tbody>
    <tr>
      <td>3.4.6</td>
      <td>Comparison with atomistic simulations</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Application to the sink strength of semicoherent interfaces</td>
      <td>66</td>
    </tr>
    <tr>
      <td>3.5.1</td>
      <td>Computational multi-model strategy</td>
      <td>66</td>
    </tr>
    <tr>
      <td>3.5.2</td>
      <td>Kinetic Monte Carlo simulations with elastic interactions</td>
      <td>68</td>
    </tr>
    <tr>
      <td>3.5.3</td>
      <td>Effect of elastic interactions on interface sink strength</td>
      <td>71</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Elastic strain relaxation in interfacial dislocation patterns</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.6.1</td>
      <td>General considerations on hexagonal-shaped dislocation patterns</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.6.2</td>
      <td>Solution methodology for strain-relaxed rearrangements</td>
      <td>77</td>
    </tr>
    <tr>
      <td>3.6.3</td>
      <td>Parametric energy-based framework</td>
      <td>78</td>
    </tr>
    <tr>
      <td>3.6.4</td>
      <td>Boundary conditions with surface/interface constitutive relations</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.6.5</td>
      <td>Application to Au/Cu heterosystems</td>
      <td>87</td>
    </tr>
    <tr>
      <td>3.6.6</td>
      <td>Comparison with atomistic simulations</td>
      <td>95</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>Interaction with extrinsic dislocations in bimaterials</td>
      <td>100</td>
    </tr>
    <tr>
      <td>3.7.1</td>
      <td>Extrinsic dislocation arrays and loops</td>
      <td>100</td>
    </tr>
    <tr>
      <td>3.7.2</td>
      <td>Internal forces on intrinsic and extrinsic dislocations</td>
      <td>105</td>
    </tr>
    <tr>
      <td>3.7.3</td>
      <td>On the piled-up dislocations in the (111)Cu/(011)Nb bimaterial</td>
      <td>107</td>
    </tr>
    <tr>
      <td>3.7.4</td>
      <td>Limitations</td>
      <td>114</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Extension to non-singular fields in multilayered magneto-electro-elastic plates</td>
      <td>115</td>
    </tr>
    <tr>
      <td>3.8.1</td>
      <td>Boundary-value problem and singularity-free field solutions</td>
      <td>116</td>
    </tr>
    <tr>
      <td>3.8.2</td>
      <td>A primary case: 2D bilayered composites</td>
      <td>127</td>
    </tr>
    <tr>
      <td>3.8.3</td>
      <td>Energy-based criterion for interlayers in A/B/A trilayers</td>
      <td>129</td>
    </tr>
    <tr>
      <td>3.8.4</td>
      <td>Dislocation-induced response under applied external loading</td>
      <td>138</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Conclusion and future works</td>
      <td>141</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Concluding remarks</td>
      <td>141</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Perspectives</td>
      <td>143</td>
    </tr>
    <tr>
      <td>4.2.1</td>
      <td>Thermoelasticity of semicoherent interfaces</td>
      <td>143</td>
    </tr>
    <tr>
      <td>4.2.2</td>
      <td>Distributed dislocations for periodic networks of cracks</td>
      <td>144</td>
    </tr>
    <tr>
      <td>4.2.3</td>
      <td>Towards a general treatment for {interfaces, dislocations, cracks}</td>
      <td>145</td>
    </tr>
  </tbody>
</table>
Remerciements

Je tiens à remercier les membres du jury, qui ont accepté d'évaluer ce mémoire : Brigitte Bacroix, Stéphane Berbenni, Renald Brenner, Marc Fivel et Ioan Ionescu. Je remercie plus particulièrement les rapporteurs d'avoir pris le temps si précieux de rapporter ce travail dans les moindres détails.
Merci pour vos retours positifs si encourageants !

Bien entendu, le contenu de ce travail aurait été réduit à une peau de chagrin sans les échanges constants avec mes anciens collègues de la Direction des Applications Militaires du Commissariat à l'Énergie Atomique, Christophe Denoual, Jean-Lin Dequiedt, Yves-Patrick Pellegrini et Ronan Madec. Outre-Atlantique, je mesure la chance d'avoir fait des rencontres inspirantes, en particulier avec Robert Balluffi, David Barnett, Michael Demkowicz, John Hirth, Ernian Pan, et j'en oublie, Niaz Abdorrahim, Tom Arsenlis, Sylvie Aubry, Nicolas Bertin, Wei Cai, Christian Brandl, Kedar Kolluri, Enrique Martinez, Ryan Sills, et j'en oublierai encore ! Je tiens aussi à remercier mes collègues les plus proches de l'Office, Christophe Bovet, Jean-Didier Garaud, Serge Kruch, Johann Rannou, sans désir d'exhaustivité. Je remercie Anne Tanguy d'avoir régulièrement soutenu cette habilitation, et ce depuis le début de l'aventure. Une attention particulière et amicale se tourne vers Vincent Chiaruttini, la seule personne disponible à 3h au mat' pour discuter, en partie, des correspondances théoriques et numériques entre une fissure et une dislocation... bienvenue dans le monde de cette dernière, mais arrêtons d'échanger si tard (quoique, continuons, mais n'en parlons ni à Aurélie, ni à Aurélie...). Je souhaite chaleureusement remercier toute l'équipe du secrétariat du Département Matériaux et Structures : votre aide à résoudre quotidiennement des problèmes administratifs est précieuse.

Merci enfin à un groupe spécial $A^3 = \{ \text{Achille} (4 \text{ mois}), \text{Anton} (2 \text{ ans}), \text{Aurélie} \}^*$ pour le bonheur non borné qu'il m'apporte au quotidien. Cette habilitation, qui contient les « mille-feuilles » et autres « Rubik's cubes » déjà contemplés, est aussi la votre !

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

<footnote>* Solutions régularisées et anisotropes de la contrainte normale, de cisaillement et de la densité d'énergie d'une boucle prismatique de dislocation simplement connexe, plongée dans un « mille-feuille »</footnote>
# Chapter 1

# Introduction

Interfaces in polycrystalline as well as multiphase solids of natural and synthetic origin have found their places in various applications, ranging from semiconductor devices to advanced multifunctional coatings in automobile and aeronautical industries. Remarkably, the behavior of polycrystalline materials is often reduced to the analysis of their inherent grain boundaries, while the most recent roadmaps on photonics and phononics propose to design on-demand bandgaps by tailoring the topological interface states in metamaterials. As claimed by Wolfgang Pauli, however, because "God made the bulk; the surface was invented by the devil!", the interface engineering of solid-state materials inevitably requires specific experimental and numerical contributions to describe the physical and mechanical characteristics of these internal interfaces. Far from being featureless dividing surfaces between neighboring crystals, the study of the structure and properties of homo- and hetero-phase interfaces has thus become as a central area in a broader field of the materials science and engineering.

The manuscript is divided into two chapters, considering first the thermodynamics of diffuse interfaces in chapter 2, which was developed more than a hundred years ago by Gibbs. The description of the structures and energetics of imperfect interfaces, namely semicoherent interfaces, is then treated in chapter 3. These semicoherent interfaces are also described by misfit dislocation networks that produce a lattice-invariant deformation which disrupts the uniformity of the lattice correspondence across the interfaces and thereby reduces coherency. This topic has more recently received considerable attention due to the development of high-resolution techniques and increased computational resources in recent decades.

The first introductory chapter 2 is thus concerned with the internal interfaces emerging from polymorphic structural (diffusionless) phase transformations. The formation of these solid-solid interfaces during the pressure-driven phase transitions in iron is captured by a thermodynamically consistent framework for combining nonlinear elastoplasticity and multivariant phase-field approach at large strains. Treatments of thermodynamics and kinetic relations of the phase transitions are formulated by the free energy landscape that involves the concept of reaction pathways with respect to the point group symmetry properties of both low- (cubic) and high- (hexagonal) pressure crystal lattices of iron. The phase-field formalism coupled with finite elastoplastic deformations is implemented into a three-dimensional finite element scheme and is applied to the body-centered cubic into hexagonal close-packed phase transitions under high hydrostatic compression and shock-wave conditions. The calculations exhibit the crucial role played by the plastic deformation in the morphological and microstructure evolution processes. However, the coexistence over a wide range of pressure of both cubic and hexagonal lattice structures in the interface-dominated microstructure leads, in general, to the loss of lattice coherence at the interfaces, for which the lattice correspondence across the grain boundaries and heterophase interfaces require a fine dislocation-based description of internal interfaces. It is the last objective that is covered by the main chapter 3.

Chapter 3 is therefore dedicated to the structures and energetics of heterophase interfaces. Although the simplest interface is a single isolated planar interface separating two adjacent crystals, also viewed as a planar interface in bimaterials, such an idealized interface between two dissimilar
crystals provides the essential basis for understanding the properties of interface-dominated ma-
terials. For the past ten years, the constant effort has been devoted to combining the closely related
Frank-Bilby and O-lattice techniques with the Stroh sextic formalism for the anisotropic elasticity
theory of interfacial dislocation patterns. The key formalism is used by means of a Fourier series-
based analysis to determine the reference states of semicoherent interfaces that gives rise to dis-
location arrays whose far-field elastic fields meet the condition of vanishing far-field strains and
prescribed misorientations. In accordance with the quantized Frank-Bilby equation, these inter-
face dislocation structures, which are also viewed as Volterra dislocations that have been inserted
into the reference state, generate persistent short-range elastic stresses near the interfaces. The cor-
responding energetics have been quantified and used for rapid computational design of interfaces
with tailored misfit dislocation patterns. In particular, a coupled approach with an object kinetic
Monte Carlo code has revealed that elastic interactions between radiation-induced point defects
and semicoherent interfaces lead to significant increases in interface sink strength, compared to
the case with no defect-interface interactions. The original work has also been extended to bi-
layers of finite thickness terminated with free surfaces, layered superlattices with differing layer
thicknesses as well as multilayered magneto-electro-elastic plates for semicoherent interfaces with
relaxed dislocation patterns at semicoherent interfaces including core-spreading effects. Overall,
the elastic full-field solutions have been compared with atomistic calculations for many specific
lattice structures, which provide an opportunity for rigorous validation of the anisotropic elas-
ticity theory of interfacial dislocations as well as for collaborations with individuals outside the
home laboratory.

Although the reader may be disappointed (I understand it...) not to find the content of the two
chapters combined together in a unified formalism, chapter 4 provides concluding remarks and
further directions for near future developments.
## Chapter 2

# Crystalline interfaces during solid-solid phase transitions in iron

### Selected peer-reviewed articles

[P26] N. Bruzy, C. Denoual, **A. Vattré.** *Polyphase crystal plasticity for high strain rate: application to twinning and retwinning in tantalum*. Journal of the Mechanics and Physics of Solids, 104921, 2022.

[P22] **A. Vattré**, C. Denoual. *Continuum nonlinear dynamics of unstable shock waves induced by structural phase transformations in iron*. Journal of the Mechanics and Physics of Solids, 131, 387-403, 2019.

[P12] C. Denoual, **A. Vattré.** *A phase field approach with a reaction pathways-based potential to model reconstructive martensitic transformations with a large number of variants*. Journal of the Mechanics and Physics of Solids, 90, 91-107, 2016.

[P11] **A. Vattré**, C. Denoual. *Polymorphism of iron at high pressure: a 3D phase-field model for displacive transitions with finite elastoplastic deformations*. Journal of the Mechanics and Physics of Solids, 92, 1-27, 2016.

### 2.1 Motivation

The high-pressure and high-deformation states of iron (Fe) are of vital importance in many technological and sociological applications [33] as well as in geophysics due to the role of Fe properties in the Earth and telluric exoplanet internal structure [233]. Fundamental understanding of the physical and mechanical properties of Fe under extreme conditions, where the deformation state is caused by various irreversible processes (e.g. plasticity and polymorphic structural (difusionless) solid-solid phase transformations), is therefore crucial in both materials science and condensed matter physics.

The first indirect evidence of polymorphic phase transitions in iron has been discovered by [17] under shock compression. The authors reported a series of three discontinuous jumps in the velocity of the free surface and postulated that the three shock-wave structure is produced by a compressive elastic precursor (Ep wave) followed by a plastic wave (P wave), and, a third wave attributed to a phase transformation (PT wave). Wave profile measurements indicate that the onset of the phase transition occurred at a pressure of ~ 13 GPa and room temperature on the Hugoniot. Since the pioneering experiments, efforts succeeded in acquiring static high pressure X-ray diffraction analysis, where the stable ferromagnetic body-centered cubic ground state (bcc α-Fe) has shown a magnetic and structural transition to the nonmagnetic hexagonal close-packed phase (hcp ε-Fe) at about 13 GPa, revealing the same transition as in shock experiments. Therefore, both bcc and hcp phases have been observed to coexist over a wide range of pressure, which captures
the signature of a diffusionless solid-to-solid martensitic transition in iron. While the phase diagram of iron under hydrostatic pressure is well established [217], detailed in situ observations via dynamic X-ray diffraction techniques during shock-loading have supported unambiguously that the high pressure phase has hcp crystal structure [139, 287]. However, due to the considerable experimental difficulties of quantifying plasticity with respect to the polymorphic phase transformations during shock wave propagation in solids, the complete irreversible deformation mechanism still remains poorly investigated.

The high pressure-induced transition in iron has been intensively described using ab-initio electronic structure calculations, where some simulation results remain debated. Although the broad outline of the transition has been settled by crystallographic considerations [49, 179, 23], a major problem deals with the accuracy in determining the energy landscape for the bcc-to-hcp transition [84, 170]. Furthermore, ab-initio computational resources are limited to small system sizes, for which plasticity-induced effects in iron cannot be captured by first-principles calculations. Alternative approaches are based on large-scale molecular dynamics simulations that give insight into the motion of multi-million-atoms. Shock waves have also been simulated by employing embedded atom method potentials and varying initial shock strength [137, 138, 136]. For low particle velocities, an elastic shock wave of uniaxially compressed bcc was observed. With increasing shock strength, a two-wave shock structure was identified with an elastic precursor followed by a slower phase-transition wave. No direct evidence of plastic wave profile was observed, certainly due to the small time scale compared to experiments that exhibit a three-wave structure at the nanosecond scale [17, 19]. While further work is needed to understand the detailed mechanisms of plasticity under shock conditions, phase-field models provide a companion approach to shock response of crystalline materials at higher time and length scales.

Various continuum mechanics approaches to simulate martensitic phase transitions in the context of plasticity theory have been developed and can be categorized by the nature of the scale description of the constitutive relations. A first micromechanical class of models aims to deliver predictions of macroscopic observables, e.g. stress-strain curves, by including microstructural aspects via homogenization and averaging techniques. In a multiscale strategy, relevant approaches track the volume fraction of martensite phase in the small [129, 215] and large [148, 177] strain formulations. However, these models are generally unable to predict detailed microstructural changes and spatial arrangements of parent-product interfaces during phase transformations at the nanometer scale. A second class of models for displacive transformations has pushed toward smaller scales in an effort to capture transformational processes by tracking the kinetics of interface orientations and variants with respect to the associated configurational forces. Thus, structural phase-field approaches have been successfully applied to model microstructure evolution by formulating thermodynamic driving forces for martensitic transitions between stable states [157, 10, 142, 77, 291]. Treatments of thermodynamics and kinetic relations in phase-field approaches are related to the pioneering works by [50] and [3], within which a material system tends to evolve towards a minimum state of free energy.

Chapter 2 introduces a thermodynamically consistent framework for combining nonlinear elastoplasticity and multivariant phase-field approach at large strains [252]. In accordance with the Clausius-Duhem inequality in section 2.2, the Helmholtz free energy and time-dependent constitutive relations give rise to displacive driving forces for pressure-induced martensitic phase transitions in materials. Inelastic forces are obtained by using a representation of the energy landscape that involves the concept of reaction pathways with respect to the point group symmetry operations of crystal lattices [76]. Using the element-free Galerkin method with high-performance computing resources, the finite deformation framework is used to analyze the polymorphic $\alpha$- into $\epsilon$-Fe iron phase transitions under high hydrostatic compression [252] and shock-wave [253] loadings, as detailed in sections 2.3 and 2.4, respectively, while a recent application to twinning and retwinning in tantalum can be found in Ref. [44]. The three-dimensional nonlinear simulations accurately reproduce observable characteristics reported by the experimental literature, for which
the crucial role played by the plastic deformation is analyzed with respect to the peculiar formation of interface-dominated microstructure with a specific selection of high-pressure variants.

## 2.2 A phase-field model coupled with finite elastoplasticity

This section is concerned with a thermodynamically consistent phase-field formalism for solid-state transitions. The model is formulated in a Lagrangian framework for finite strains, motivated by obtaining isothermal driving forces and constitutive relations at a material point.

### 2.2.1 Kinematics

An arbitrary material point X is defined in a homogeneous reference configuration $\Omega_0 \subset \mathbb{R}^3$, for which the motion of $\Omega_0$ is given by the mapping $x = \chi(X, t) : \Omega_0 \to \Omega \subset \mathbb{R}^3$ with respect to time t. The total deformation gradient F is related to the following multiplicative decomposition [158, 159, 148, 161], i.e.,

$$F = \frac{\partial \chi}{\partial X}|_t = \nabla \chi = \mathbf{Fe} \cdot \mathbf{Ft} \cdot \mathbf{Fp}, \qquad (2.1)$$

with $\nabla$ the material gradient with respect to X. Here, the reference configuration is associated with the initial single-crystal bcc iron, and, the total deformation gradient is decomposed into elastic Fe, plastic Fp, and, transformational Ft distortions, leading to the pressure-induced phase transformation from the bcc to hcp phases.

Similarly to classical crystal elastoplasticity theories [151, 156], the decomposition eq. (2.1) is not uniquely defined and different ordering relations have been taken into account in the literature [241]. Because the local irreversible plastic deformation Fp of the neighborhood of X, e.g. caused by dislocation glides, does not alter the crystal orientation and structure of the lattice vectors, the transformational component Ft occurs between Fp and Fe, where the elastic contribution accounts for the lattice stretching Ue and rotation Re. The polar decomposition to Fe reads: $Fe = Re \cdot Ue$, with $Ue^2 = Fe^t \cdot Fe$, and, $det Fe = det Ue = j_e$. The superscript t denotes the transpose operation. Although the controversy regarding the decomposition is beyond the scope of this paper, both tensors Fp and Ft describe here two intermediate configurations, $\Omega p$ and $\Omega t$, as shown in Fig. (2.1). For more justifications regarding the three-term multiplication decomposition eq. (2.1) for nonlinear elasticity coupled to martensitic phase transformations and plasticity, the reader is referred to the recent analysis on combined kinematics in Ref. [161]. It follows from eq. (2.1) that the total spatial velocity gradient tensor L is given by

$$L = \dot{F} \cdot F^{-1} = Le + Fe \cdot Lt \cdot Fe^{-1} + Fet \cdot Lp \cdot Fet^{-1}, \qquad (2.2)$$

with $Fet = Fe \cdot Ft$. The superposed dot in eq. (2.2) denotes the time derivative. The elastic Le, transformational Lt, and, plastic Lp velocity strain tensors are similarly defined by

$$Le = \dot{Fe} \cdot Fe^{-1}, \quad Lt = \dot{Ft} \cdot Ft^{-1}, \quad \text{and} \quad Lp = \dot{Fp} \cdot Fp^{-1}, \qquad (2.3)$$

which are related to the current and the intermediate configurations, i.e., $\Omega$, $\Omega t$ and $\Omega p$, respectively. Furthermore, two basic kinematic assumptions are considered in the present theory:

1. The measures of volume changes after each deformation processes satisfy:

$$j_e = \det Fe > 0, \quad j_t = \det Ft > 0, \quad \text{and} \quad \det Fp = 1, \qquad (2.4)$$

so that, Fe, Ft and Fp are invertible, and, the plastic flow preserves the volume.
<figure><img src="image_3.png" /><figcaption>FIGURE 2.1: Schematics of the reference $\Omega_0$, intermediate, $\Omega_p$ and $\Omega_t$, and, current $\Omega$ configurations, for which the total deformation gradient tensor $F$ is decomposed multiplicatively into plastic $F_p$, transformational $F_t$ and elastic $F_e$ distortions.</figcaption></figure>

2. The model is restricted to isotropic plastic theories with irrotational plastic flows. Therefore,

$$ \dot{F}_p = D_p \cdot F_p, \text{ with } D_p = \text{sym } L_p = L_p, \qquad (2.5) $$

where $\text{sym } L_p$ denotes the symmetric part of $L_p$.

Figure (2.1) illustrates the multiplicative split of the total deformation gradient tensor $F$. In agreement with the conservation law of mass, the determinant of $F$ gives the volume change between the current (with a volume $V$) and the reference ($V_0$) configurations, i.e., $j = \det F = \rho_0/\rho = V/V_0$, where $\rho$ ($\rho_0$ with $\dot{\rho}_0 = 0$) is the current (reference) mass density.

### 2.2.2 Balance laws

During the different deformation processes, the equilibrium equations of force must be fulfilled. In the Lagrangian description, the local form of the linear momentum balance is given by

$$ \nabla \cdot \mathbf{P} + \rho_0 \mathbf{b} = \rho_0 \ddot{\mathbf{u}} \quad \text{in } \Omega_0, \qquad (2.6) $$

where $\mathbf{P}$ is the first (non-symmetric) Piola-Kirchhoff stress tensor, $\mathbf{b}$ are external body forces per unit mass, and, $\ddot{\mathbf{u}} = \ddot{\mathbf{x}}(\mathbf{X}, t)$ is the acceleration of the material point $\mathbf{X}$, with $\mathbf{u}$ the corresponding displacement field, defined by $\mathbf{u} = \mathbf{x}(\mathbf{X}, t) - \mathbf{X}$.

An appropriate formulation of the constitutive relations for isothermal and irreversible processes of deformation requires a thermodynamically consistent formalism, within which the balance law in eq. (2.6) holds at all points $\mathbf{X}$ in the domain of $\Omega_0$.

### 2.2.3 The Clausius-Duhem inequality

The martensitic phase-field approach coupled with large elastoplastic deformations is derived within a thermodynamic framework in which the second law of thermodynamics plays a crucial role. When the thermal effects are ignored, the fundamental Clausius-Duhem inequality is expressed in terms of stress power per unit reference volume [67] as

$$ \int_{\Omega_0} (\mathbf{P} : \dot{\mathbf{F}} - \rho_0 \dot{\psi}) \, d\Omega_0 \geq 0, \qquad (2.7) $$

where $:$ denotes the double inner tensor product, and, $\psi$ the specific Helmholtz free energy. Equation (2.7) shows that the first Piola-Kirchhoff stress tensor $\mathbf{P}$ and the deformation gradient $\mathbf{F}$ are work-conjugate variables, while $\mathbf{P}:\dot{\mathbf{F}}$ defines the mechanical stress power per unit volume in the Lagrangian formulation.
Within the model of the multiplicative decomposition in finite strains, it is conveniently pos-
tulated that the Helmholtz free energy can be written in the following form:

$$
\psi \doteq \psi (\mathbf{Fe}, \mathbf{Ft}, \nabla \mathbf{Ft}), \qquad (2.8)
$$

where **∇**Ft is a phenomenological third-order gradient term that acts as a penalty for spatial nonuniformity to produce diffuse interfaces. Because the elastic response is not affected by the plastic activities, the elastic part of the Helmholtz free energy is supposed to depend on the elastic and transformational distortions only. Moreover, it is assumed that both transformational and plastic works are not dependent on each other, so that the free energy may be additively decomposed into elastic ψₑ, transformational ψₜ, and, purely empirical gradient penalty ψ_∇ contributions. With the aforementioned considerations, the Helmholtz free energy can thus be written as

$$
\psi \doteq \psi_e(\mathbf{Fe}, \mathbf{Ft}) + \psi_t(\mathbf{Ft}) + \psi_\nabla(\nabla \mathbf{Ft}), \quad (2.9)
$$

which, in contrast with ab-initio electronic structure calculations, is not uniquely defined. How-
ever, such elastic/inelastic splitting, comparable to the classical phase-field models with elastic
and chemical potentials [278, 10], is fundamental for applications that exhibit a strong coupling
between acoustic waves and phase transformations, e.g. wave propagation influencing the early
stages of the phase transitions induced by shock loadings. Thus, eqs. (2.1) and (2.9) yield to the
rates of the total deformation and free energy, i.e.,

$$
\begin{align}
\dot{\mathbf{F}} &= \dot{\mathbf{Fe}} \cdot \mathbf{Ft} \cdot \mathbf{Fp} + \mathbf{Fe} \cdot \dot{\mathbf{Ft}} \cdot \mathbf{Fp} + \mathbf{Fe} \cdot \mathbf{Ft} \cdot \dot{\mathbf{Fp}} \nonumber \\
\dot{\psi} &= \left. \frac{\partial \psi_e}{\partial \mathbf{Fe}} \right|_{\mathbf{Ft}} : \dot{\mathbf{Fe}} + \left. \frac{\partial \psi_e}{\partial \mathbf{Ft}} \right|_{\mathbf{Fe}} : \dot{\mathbf{Ft}} + \frac{\partial \psi_t}{\partial \mathbf{Ft}} : \dot{\mathbf{Ft}} + \frac{\partial \psi_\nabla}{\partial \nabla \mathbf{Ft}} : \nabla \dot{\mathbf{Ft}}, \tag{2.10}
\end{align}
$$

where . denotes the triple inner tensor product. Inserting eqs. (2.10) into the global form of the
Clausius-Duhem inequality (2.7) and applying the chain rule, the non-negative requirement leads
therefore to

$$
\int_{\Omega_0} \left\{ \left( \mathbf{P} \cdot \mathbf{Fp}^t \cdot \mathbf{Ft}^t - \rho_0 \frac{\partial \psi_e}{\partial \mathbf{Fe}} \Big|_{\mathbf{Ft}} \right) \cdot \dot{\mathbf{Fe}} + \left( \mathbf{Fe}^t \cdot \mathbf{P} \cdot \mathbf{Fp}^t - \rho_0 \frac{\partial \psi_e}{\partial \mathbf{Ft}} \Big|_{\mathbf{Fe}} - \rho_0 \frac{\partial \psi_t}{\partial \mathbf{Ft}} \right) \cdot \dot{\mathbf{Ft}} + \boldsymbol{\Sigma}_*: \mathbf{Dp} \right. \\
\left. - \rho_0 \frac{\partial \psi_\nabla}{\partial \nabla \mathbf{Ft}} : \nabla \dot{\mathbf{Ft}} \right\} d\Omega_0 \ge 0,
\quad (2.11)
$$

where Σ\* is a work-conjugate stress measure related to the first Piola-Kirchhoff P, as follows

$$
\Sigma_* = \mathrm{Fet}^t \cdot \mathbf{P} \cdot \mathbf{Fp}^t . \qquad (2.12)
$$

Using the permutability of time and space differentiation in the reference configuration and
the Gauss theorem, the last right-hand side term in eq. (2.11) can be rewritten, i.e.,

$$
\int_{\Omega_0} \left( \frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{Ft}} : \cdot \nabla \mathbf{Ft} \right) d\Omega_0 = - \int_{\Omega_0} \left( \nabla \cdot \frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{Ft}} : \mathbf{Ft} \right) d\Omega_0 + \int_{\Sigma_0} \underbrace{\left( \mathbf{Ft} : \frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{Ft}} \cdot \mathbf{n} \right)}_{\text{surface dissipation}} d\Sigma_0, \quad (2.13)
$$

where Σ₀ is a boundary of Ω₀ with unit outward normal n. Assuming that the surface dissipa-
tion is absent during the transformational process, additional boundary conditions as set of nine
equations for phase transitions may also be derived by

$$
\frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{Ft}} \cdot \mathbf{n} = 0, \text{ with } \dot{\mathbf{Ft}} \neq 0 \text{ at } \Sigma_0, \qquad (2.14)
$$
corresponding to the orthogonality relations between $\nabla Ft$ and the external surfaces $\Sigma_0$. Thus, eqs. (2.11–2.14) yield to a local formulation of the free energy imbalance in terms of dissipation per unit reference volume of mechanical energy $D$, as follows

$$D = \left( \mathbf{P} \cdot \mathbf{Fp}^{\mathrm{t}} \cdot \mathbf{Ft}^{\mathrm{t}} - \rho_0 \frac{\partial \psi_e}{\partial \mathbf{Ft}} \bigg|_{\mathbf{Ft}} \right) : \dot{\mathbf{F}}e + \mathbf{Xt} : \dot{\mathbf{F}}t + \mathbf{\Sigma}_* : \mathbf{Dp} \geq 0, \quad (2.15)$$

where the dissipative forces $\mathbf{Xt}$, conjugated to dissipative rate $\dot{\mathbf{F}}t$, are given by

$$\mathbf{Xt} = \mathbf{Fe}^{\mathrm{t}} \cdot \mathbf{P} \cdot \mathbf{Fp}^{\mathrm{t}} - \rho_0 \frac{\partial \psi_e}{\partial \mathbf{Ft}} \bigg|_{\mathbf{Fe}} - \rho_0 \frac{\partial \psi_t}{\partial \mathbf{Ft}} + \rho_0 \nabla \cdot \frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{Ft}}. \quad (2.16)$$

The relation (2.16) defines the thermodynamic displacive driving forces for change in $\mathbf{Ft}$, acting on a material point $X$ under isothermal conditions. Although the plastic deformation is not integrated as an internal state variable, e.g. via a defect-energy term as in Refs. [107, 2], but rather as a kinematic variable, the plastic contribution may significantly alter the state of residual stress and also play an important role in dictating the morphology of the microstructural changes and in modeling the irreversibility of phase transitions.

### 2.2.4 Constitutive equations

Constitutive equations for reversible elastic deformations and irreversible processes of deformable material bodies undergoing phase and plastic deformations are required to be consistent with the Clausius-Duhem inequality.

#### Hyperelasticity

The standard assumption that the rate of dissipation is independent of $\dot{\mathbf{F}}e$ in eq. (2.15), i.e., elasticity is a non-dissipative process, results in the hyperelasticity constitutive relation in terms of the first Piola-Kirchhoff stress field, as follows

$$\mathbf{P} = \rho_0 \frac{\partial \psi_e}{\partial \mathbf{Fe}} \bigg|_{\mathbf{Ft}} \cdot \mathbf{Ft}^{-\mathrm{t}} \cdot \mathbf{Fp}^{-\mathrm{t}}. \quad (2.17)$$

A quadratic form for the strain energy density per unit reference volume is assumed, for which a dependence of $\psi_e$ on $\mathbf{Fe}$ and $\mathbf{Ft}$ manifests explicitly via the anisotropic elastic components:

$$\rho_0 \psi_e = \frac{1}{2} \mathbf{Ee} : \mathbf{D}(\mathbf{Cet}) : \mathbf{Ee}, \quad (2.18)$$

where $\mathbf{Ee}$ is the elastic Green-Lagrange strain tensor, defined by

$$\mathbf{Ee} = \frac{1}{2} (\mathbf{Ce} - \mathbf{I}), \quad (2.19)$$

with $\mathbf{Ce} = \mathbf{Fe}^{\mathrm{t}} \cdot \mathbf{Fe}$ the right elastic Cauchy-Green deformation tensor, so that $\mathbf{Cet} = \mathbf{Ft}^{\mathrm{t}} \cdot \mathbf{Ce} \cdot \mathbf{Ft}$. Inserting eq. (2.18) into the hyperelasticity condition (2.17), the nonlinear stress-elastic strain constitutive relation is rewritten as follows

$$\mathbf{P} = \mathbf{Fe} \cdot \mathbf{Se} \cdot \mathbf{Ft}^{-\mathrm{t}} \cdot \mathbf{Fp}^{-\mathrm{t}} + \mathbf{Fet} \cdot \left( \mathbf{Ee} : \frac{\partial \mathbf{D}(\mathbf{Cet})}{\partial \mathbf{Cet}} : \mathbf{Ee} \right) \cdot \mathbf{Fp}^{-\mathrm{t}}, \quad (2.20)$$

where $\mathbf{Se} = \mathbf{D}(\mathbf{Cet}) : \mathbf{Ee}$ is an elastic stress measure associated with $\mathbf{Ee}$, and, $\partial_{\mathbf{Cet}} \mathbf{D}$ is a sixth-order tensor, i.e., the derivative of $\mathbf{D}$ with respect of $\mathbf{Cet}$. It is worth pointing out that the anisotropic pressure-dependent elastic stiffness tensors of both bcc and hcp phases are explicitly taken into account in the present formalism.
With use of the non-dissipative properties of hyperelasticity, the local dissipation considered in the Clausius-Duhem inequality (2.15) can also be conceptually divided into transformational $D_t$ and plastic $D_p$ dissipative rates per unit reference volume, i.e.,

$$D \doteq D_t + D_p \geq 0, \qquad (2.21)$$

due to the onset of the phase transitions or the movements of interface during phase transitions, and, to the plastic deformation in materials, respectively. For simplicity, it is assumed that both transformational and plastic dissipative processes are thermodynamically uncoupled such that the inequality (2.21) splits into two stronger non-negative inequalities, as follows

$$D_t = \mathbf{X}t: \dot{\mathbf{F}}t \geq 0 \text{ and, } D_p = \mathbf{\Sigma}_*: \mathbf{D}p \geq 0. \qquad (2.22)$$

Kinetic constitutive relations that relate the rates $\dot{\mathbf{F}}t$ and $\mathbf{D}p$ to the associated driving forces for both dissipative processes in hyperelastic materials must also be defined such that the inequalities in eqs. (2.22) are satisfied. These steps are carried out in the two subsequent sections.

## Kinetics of phase transitions

For solid-state structural transformations, a linear kinetic equation that relates the rate of the transformational distortion $\dot{\mathbf{F}}t$ to the displacive driving forces $\mathbf{X}t$ is suggested, i.e.,

$$v \dot{\mathbf{F}}t = \mathbf{X}t, \qquad (2.23)$$

where $v > 0$ is a viscosity-like parameter. For example, the case with $v \to 0$ represents an instantaneous relaxation. The evaluation of the kinetic equations for martensitic phase transitions is still a subject of intense debates, within which the average transformational kinetics may be influenced by the nucleation processes, interface mobilities, collective dislocation behaviors, as well as inertial effects. In the context of the time-dependent Ginzburg-Landau formalism, a detailed modeling of the kinetics of phase transitions in iron is not the purpose of the present analysis. However, the linear form of the driving forces $\mathbf{X}t$ gives rise to thermodynamic consistency conditions for phase transformations, so that the dissipation inequality in eq. (2.22) is unequivocally satisfied, as follows

$$D_t = v |\mathbf{X}t|^2 \geq 0, \qquad (2.24)$$

with $|\mathbf{X}t|$ the Frobenius norm of $\mathbf{X}t$. A nonequilibrium thermodynamic system is also characterized when $D_t > 0$, e.g. corresponding to mobile solid-solid interfaces when $\mathbf{X}t > 0$. Using eqs. (2.16) and (2.20), eq. (2.23) yields

$$v \dot{\mathbf{F}}t = \mathbf{X}t = \underbrace{\mathbf{C}\mathbf{e} \cdot (\mathbf{D}(\mathbf{C}\mathbf{e}):\mathbf{E}\mathbf{e}) \cdot \mathbf{F}t^{-t}}_{\text{forces due to elastic energy}} - \underbrace{\rho_0 \frac{\partial \psi_t}{\partial \mathbf{F}t} + \rho_0 \nabla \cdot \frac{\partial \psi_{\nabla}}{\partial \nabla \mathbf{F}t}}_{\text{transformational forces}}, \qquad (2.25)$$

including mechanical elastically and transformational inelastically induced driving forces, with a gradient-related term for interface energy. Equation (2.25) shows competition between driving forces due to elastic energy and the inelastic transformational forces related to microstructure evolution processes in materials. In particular, the (meta)stable equilibrium configurations are achieved when $\mathbf{X}t = 0$, exhibiting a force balance between the elastic and inelastic contributions.

A general quadratic form for the gradient energy penalty that is localized at the diffuse interfaces between two phases may be defined by

$$\rho_0 \psi_{\nabla} = \frac{1}{2} \nabla \mathbf{F}t \cdot \mathbf{\Lambda} \cdot \nabla \mathbf{F}t, \qquad (2.26)$$
where $^6\Lambda$ is a positive definite symmetric (major symmetry) sixth-order tensor that takes into account the gradient-energy interaction between different phases. Assuming an isotropic description of the interface energy and neglecting the interactions between all phases [160] such that $^6\Lambda = \lambda$ $^6\mathbf{I}$, with $^6\mathbf{I}$ the sixth-rank identity tensor, eq. (2.26) reduces to

$$ \rho_0 \psi_{\nabla} = \frac{1}{2} \lambda \nabla \mathbf{F} t \cdot \dots \mathbf{I} \cdot \nabla \mathbf{F} t = \frac{1}{2} \lambda \nabla \mathbf{F} t \cdot \nabla \mathbf{F} t, \quad (2.27) $$

where the positive scalar $\lambda$ controls phenomenologically the finite width of interfaces. The latter distance may be correlated to the short-range elastic fields produced by discrete intrinsic dislocation arrays between bcc/hcp semicoherent heterophase interfaces and also computed by using a recent formalism linking the Frank-Bilby equation and anisotropic elasticity theory, as investigated in chapter 3. Finally, the driving forces expressed in the Ginzburg-Landau formalism are given by

$$ v \dot{\mathbf{F}} t = \mathbf{X} t = \mathbf{C} e \cdot (\mathbf{D}(\mathbf{C} e) : \mathbf{E} e) \cdot \mathbf{F} t^{-t} - \rho_0 \frac{\partial \psi_t}{\partial \mathbf{F} t} + \lambda \nabla^2 \mathbf{F} t, \quad (2.28) $$

with $\nabla^2$ the Laplacian operator.

### Plastic flow rule

Macroscopic quasi-perfectly plastic regimes have been observed in polycrystalline bcc iron samples under high-strain rate compressions [132]. To go beyond the elastic limit, the large strain perfectly plastic $J_2$ flow theory has also been incorporated in the present model. Accordingly, the evolution of the plastic distortion $\mathbf{F}p$, given in terms of the plastic deformation rate $\mathbf{D}p$, is determined by considering the postulate of maximum dissipation [118]. The space of admissible stresses $\mathcal{E}_\sigma$ is written as

$$ \mathcal{E}_\sigma = \{ \sigma | \phi(\sigma) < 0 \}, \quad (2.29) $$

where the yield function $\phi$ is expressed in terms of the Cauchy stress $\sigma$, defined by

$$ \sigma = j^{-1} \mathbf{P} \cdot \mathbf{F}^t = j^{-1} \mathbf{F} e \cdot \mathbf{S} e \cdot \mathbf{F} e^t + j^{-1} \mathbf{F} \mathrm{f} \cdot \left( \mathbf{E} e : \frac{\partial \mathbf{D}(\mathbf{C} e)}{\partial \mathbf{C} e} : \mathbf{E} e \right) \cdot \mathbf{F} \mathrm{f}^t, \quad (2.30) $$

according to eq. (2.20). The work-conjugate stress $\Sigma_*$ in eq. (2.12) may also be related to the Cauchy stress tensor $\sigma$ by

$$ \Sigma_* = j \mathbf{F} \mathrm{f}^t \cdot \sigma \cdot \mathbf{F} \mathrm{f}^{-t} = \mathbf{F} \mathrm{f}^t \cdot \Sigma \cdot \mathbf{F} \mathrm{f}^{-t}, \quad (2.31) $$

where $\Sigma = j \mathbf{F} \mathrm{f}^t \cdot \sigma \cdot \mathbf{F} \mathrm{f}^{-t}$. Thus, the rate of plastic deformation $\mathbf{D}p$ is given by the associated flow rule, as follows

$$ \mathbf{D}p = \dot{\eta} \mathbf{F} \mathrm{f}^{-1} \cdot \frac{\partial \phi}{\partial \sigma} \cdot \mathbf{F} \mathrm{f} = \dot{\eta} \mathbf{H}, \quad (2.32) $$

with $\dot{\eta} \ge 0$ a non-negative scalar-valued factor, so-called the plastic multiplier, that is required to satisfy the consistency relation: $\dot{\eta} \phi = 0$. The outward normal to the yield surface is given by $\mathbf{H}$ in the stress space, for which the yield function $\phi$ in eqs. (2.29) and (2.32) is described with the von Mises yield criterion, i.e.,

$$ \phi(\sigma) = \sqrt{3 J_2(\sigma)} - \sigma_0 \quad \text{with,} \quad J_2 = \frac{1}{2} \mathrm{dev} \sigma : \mathrm{dev} \sigma, \quad (2.33) $$

where $\sigma_0 > 0$ is the yield stress measure, and, $\mathrm{dev} \sigma$ denotes the deviatoric part of $\sigma$. Finally, including the direction of the plastic flow into the rate $\mathbf{D}p$, eq. (2.32) yields

$$ \mathbf{D}p = \frac{3}{2} \dot{\eta} \mathbf{F} \mathrm{f}^{-1} \cdot \frac{\mathrm{dev} \sigma}{\sigma_0} \cdot \mathbf{F} \mathrm{f}, \quad (2.34) $$
for which the dissipation inequality for plastic flow in eq. (2.22) with (2.31) is satisfied, i.e.,

$$ \mathcal{D}_p = \frac{3}{2} j \dot{j} \frac{|\mathrm{dev}\sigma|^2}{\sigma_0} \geq 0. \qquad (2.35) $$

According to eqs. (2.24) and (2.35), the present formalism is also thermodynamically consistent since the Clausius-Duhem inequality (2.21) is fulfilled.

### 2.2.5 Multiple reaction pathways and energy landscape

In what follows in section 2.2.5, focus is on the $\alpha \leftrightarrow \epsilon$ phase transitions in iron, where the energy landscape is defined by reaction pathways for multivariants with respect to the point group symmetry properties of the bcc and hcp lattices.

#### The bcc-to-hcp transition mechanism

As illustrated in Fig. (2.2a), the considered crystallographic relations in the bcc-to-hcp martensitic phase transition are given by the Mao-Bassett-Takahashi mechanism [179], as follows

$$ [001]_{\text{bcc}} \parallel [2\bar{1}\bar{1}0]_{\text{hcp}} \text{ and, } (110)_{\text{bcc}} \parallel (0001)_{\text{hcp}}, \qquad (2.36) $$

which differs from the transformation path proposed in Ref. [49] by a rotation of $\sim \pm 5.2^\circ$ around the $[0001]_{\text{hcp}}$ axis [271]. The structural relations in eq. (2.36) are achieved by considering two transformation operations, as shown in Fig. (2.2b). The hcp phase may be obtained by applying a shear to a $(110)_{\text{bcc}}$ plane, which consists of an elongation and a compression along the $[1\bar{1}0]_{\text{bcc}}$ and the $[001]_{\text{bcc}}$ directions, respectively. This transformation is required to form a regular hexagon (in red in Fig. 2.2b) and may be related to a homogeneous linear mapping $\mathbf{U}$, i.e.,

$$ \mathbf{U} = \begin{bmatrix} \frac{3}{4\sqrt{2}} + \frac{1}{4}\sqrt{\frac{3}{2}}c/a & -\frac{3}{4\sqrt{2}} + \frac{1}{4}\sqrt{\frac{3}{2}}c/a & 0 \\ -\frac{3}{4\sqrt{2}} + \frac{1}{4}\sqrt{\frac{3}{2}}c/a & \frac{3}{4\sqrt{2}} + \frac{1}{4}\sqrt{\frac{3}{2}}c/a & 0 \\ 0 & 0 & \frac{\sqrt{3}}{2} \end{bmatrix}, \qquad (2.37) $$

where $c/a = c/a$ is the lattice ratio for the pure $\epsilon$-Fe phase [55], while the volume change accompanying the phase transition is determined by $\det \mathbf{U} = 9c/a/16$. Then, the mechanism involves a shuffle $t$, which corresponds to atomic displacements of every other deformed $(110)_{\text{bcc}}$ plane in one of the two possible opposite $[1\bar{1}0]_{\text{bcc}}$ directions. The close-packed structure of hcp is also obtained, where a ratio $c/a$ of $1.603 \pm 0.001$ has been experimentally determined along this bcc-to-hcp path in iron [179, 78], reflecting a $\sim 10\%$ volume reduction.

In the described case, the transformations $\mathbf{U}$ and $t$ are illustrated separately but can occur simultaneously, as discussed by using first-principles simulations [81]. For both scenarios, the shuffle does not induce any lattice-distortion transformations and has therefore no direct coupling with the overall stress in the deforming materials. Although not visible for a given deformation state at the macroscopic scale, the shuffling modes, however, may have important implications on the free energy along the reaction pathways as well as the kinetics of phase transitions, which are not taken into account in the present formalism. Assuming to take place at a smaller time scale compared to the lattice-distortion transformations, additional variables of state (also, additional associated kinetic equations) should therefore be introduced to characterize such atomistic displacements. With the aforementioned considerations, and because of the required number of finite element meshes for three-dimensional calculations, the example applications to high-pressure compression in sections 2.3 and 2.4 focus on the first cycle of forward and reverse martensitic transitions only, for which the shuffle does not modify the point group symmetries. For higher-order
cycles, this mechanism may be responsible for the generation of an unbounded set of variants. The notion of transformation cycles has been addressed in Ref. [76], where a two-dimensional simulation has shown that variants could hierarchically nucleate into previously created ones over up to five levels of transformations for the square to hexagonal martensitic phase transitions.

When $c/a$ is experimentally chosen to determine eq. (2.37), the corresponding homogeneous mapping $U$ contains obviously and inseparably both elastic and irreversible part of the deformation in samples. A homogeneous distortion $Ut$ is therefore introduced to identify the pure transformational component of the total deformation provided by experimental data under high hydrostatic pressure, i.e.,

$$ Ut = \kappa U, \qquad (2.38) $$

where $\kappa$ is an elastic correction factor, as discussed in Ref. [252].

<figure><img src="image_4.png" /><figcaption>FIGURE 2.2: Crystallographic relations in the bcc-to-hcp martensitic phase transition established in Refs. [179, 23]. (a) Red atoms in a bcc atomic-side unit cell are located at a (110)$_{bcc}$ layer and the blue atoms at the adjacent layers. (b) The transition path consists of two transformations. First, a shear deformation $U$ leads to an elongation and a compression along the [110]$_{bcc}$ and the [001]$_{bcc}$ directions, respectively. The deformation transforms a polygon in blue into a regular hexagon in red, corresponding to the (0001)$_{hcp}$ hcp basal plane. Then, a shuffle $t$ is applied to the entire plane that contains the blue atoms, e.g. by shifting all these atoms in the [110]$_{bcc}$ direction.</figcaption></figure>

### Multiple symmetry-related variants

During the forward $\alpha \to \epsilon$ and the reverse $\epsilon \to \alpha'$ martensitic transformations, significant differences in orientation from the initial $\alpha$-Fe phase may exist. To make the clear distinction in phase orientation between variant formation and selection, $\alpha'$ denotes here the reversed $\alpha$ phase, as depicted by the two-dimensional schematic network in Fig. (2.3a).

A rigorous link between the standard crystallographic concepts of holohedry with group-subgroup relationships, crystal system and Bravais lattice type (cubic and hexagonal), is explicitly included into the phase-field formalism. For the forward $\alpha \to \epsilon$ transition, the generation of all hcp variants $\epsilon_n^\alpha Ut$ from the linear mapping $Ut$ is described by

$$ \epsilon_n^\alpha Ut = R_{\text{bcc}}^t \cdot Ut \cdot R_{\text{bcc}}, \qquad (2.39) $$
where $R_{bcc}$ is a rotation matrix in the point group of cubic lattice $n\mathcal{H}_{bcc}$ and $n$ the number of hcp variants [206]. Because of the high symmetry of the considered phase, a total number of 6 hcp variants are generated, i.e., $n = 1, \dots, 6$, within which 18 operations in the basic group of 24 rotations for cubic lattices are redundant. To complete the phase transformations with the reverse $\epsilon \to \alpha'$ transitions, the bcc variants $\alpha'_{cn} Ut$ are deduced by performing the following operation:

$$ \alpha'_{cn} Ut = R_{bcc}^t \cdot R_{hcp}^t \cdot Ut^{-1} \cdot R_{hcp} \cdot Ut \cdot R_{bcc}, \qquad (2.40) $$

where $R_{hcp}$ is a rotation matrix in the point group of hexagonal lattice $m\mathcal{H}_{hcp}$ and $m$ the number of bcc variants [206]. Equation (2.40) consists in generating 12 bcc variants, i.e., $m = 1, \dots, 12$, so that a total of 19 variants (including the identity as the 19th variant) are identified to describe the complete bcc-hcp-bcc transition in terms of multiple symmetry-related variant structures. Figure (2.3a) depicts the forward transition of the initial bcc phase, leading to six equivalent hcp phases, and, the reverse transition from each hcp phase that leads to three bcc phases.

All tabulated hcp and bcc variants with the corresponding holohedral subgroups $m\mathcal{H}_{bcc}$ and $n\mathcal{H}_{hcp}$ are given in Tab. 1 from Ref. [252], where the rotation axes are expressed in the hcp and bcc lattice basis, respectively. For clarity, the matrices defined by eqs. (2.39) and (2.40) are written in the following as $kUt$ with $k = 1, \dots, 18$, i.e.,

$$ kUt = \begin{cases} \epsilon_n Ut & 1 \le k \le 6 \\ \alpha'_{cn} Ut & 7 \le k \le 18, \end{cases} \qquad (2.41) $$

which are associated with the variant of interest $V_k$ for the forward ($1 \le k \le 6$) and the reverse ($7 \le k \le 18$) transformations.

## Reaction pathways in strain spaces

Instead of introducing the Landau thermodynamic potential [162], where the classical Landau-type approach with polynomials is not convenient to apply for reconstructive transitions due to the large numbers of potential energy wells [28], the concept of reaction pathways [77, 76] is used to describe the phase transitions in iron. In particular, the minimum inelastic energy density profile between two different pure phases is represented by a single reaction pathway, along which the associated function $\hat{\psi}_{tk}$ is assumed to possess the same symmetries as all symmetry-related variants $V_k$, and, to satisfy the principle of material objectivity [27], e.g.,

$$ \hat{\psi}_{tk} \doteq \hat{\psi}_{tk} (k\hat{C}t), \qquad (2.42) $$

where $k\hat{C}t$ are the transformational Cauchy-Green strain measures for all pure phases, given by

$$ k\hat{C}t = kUt^t \cdot kUt, \qquad (2.43) $$

as listed in Appendix A from Ref. [252], with the aid of eqs. (2.39–2.41). Here and in the following, the superimposed caret will be used to indicate quantities strictly defined along the pathways. To model continuous forward and the reverse transformations, each transition pathway $k$ is represented by linear interpolation between starting $k\hat{C}t_{start}$ and ending $k\hat{C}t_{end}$ strain states, as follows

$$ k\hat{C}t(s_k) = (1 - s_k)^k \hat{C}t_{start} + s_k^k \hat{C}t_{end}, \qquad (2.44) $$

with $s_k \in [0, 1]$ the curvilinear coordinates along $k$. For instance, hcp variants $V_k$ are parameterized by: $k\hat{C}t_{start} = I$ and $k\hat{C}t_{end} = Ut^2$, with $1 \le k \le 6$. Generating the reaction pathways with eqs. (2.37–2.44) and using projection matrices $C_1, C_2$ and $C_3$, an example of three-dimensional
<figure><img src="image_5.png" /><figcaption>FIGURE 2.3: (a) Schematic illustration of the multiple symmetry-related variants for the forward $\alpha \to \epsilon$ (in red) and the reverse $\epsilon \to \alpha'$ (blue) phase transitions in iron. (b) The corresponding reaction pathway network in a specific $\{C_1, C_2, C_3\}$ strain space, within which the transformational Cauchy-Green tensor $\mathbf{C} = \mathbf{F}^t \cdot \mathbf{F}^t$ as well as some quantities described in the text, are defined.</figcaption></figure>

representation of the network is shown in Fig. (2.3b), within which each pathway connects continuously and linearly with two pure bcc/hcp variants $V_k$ in the $\{C_1, C_2, C_3\}$ strain space. The projection is not unique and the specific strain space in Fig. (2.3b) is characterized by using the following matrices:

$$ \mathbf{C}_1 = \begin{bmatrix} 1 & 1 & 0 \\ 1 & -3 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{C}_2 = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & -3 \end{bmatrix}, \quad \mathbf{C}_3 = \begin{bmatrix} -3 & 0 & 1 \\ 0 & 0 & 0 \\ 1 & 0 & 1 \end{bmatrix}. \qquad (2.45) $$

The reaction pathway network describes also a six-dimensional energy landscape, for which each straight segment represents a minimum-energy reaction pathway that connects two stable/(meta)stable states with possible (if any) saddle points [28].

### Inelastic energy landscape

In order to define the total inelastic energy landscape $\psi_t$ in a whole strain space, e.g. not only restricted along the pathways as $\hat{\psi}_{t_k}$, the partition of unity approach is used as a weighted sum of the contribution $\psi_{t_k}$ of each individual pathway $k$. Thus, the overall inelastic energy density $\psi_t$ is formally defined by introducing the weighting functions $\omega_k(\mathbf{C}t)$, i.e.,

$$ \psi_t(\mathbf{C}t) = \sum_{k=1}^{18} \omega_k(\mathbf{C}t) \psi_{t_k}(\mathbf{C}t), \qquad (2.46) $$

for any transformational Cauchy-Green tensor $\mathbf{C} = \mathbf{F}^t \cdot \mathbf{F}^t$. Without loss of generality, for any given tensor $\mathbf{A}$, e.g. $\mathbf{C}t$ and $\mathbf{C}et$, these functions $\omega_k(\mathbf{A})$ satisfy the partition of unity condition,
namely:

$$ \sum_{k=1}^{18} \omega_k(\mathbf{A}) = 1 \quad \text{with,} \quad \omega_k(\mathbf{A}) = \frac{d_k^{-h}(\mathbf{A})}{\sum_{i=1}^{18} d_i^{-h}(\mathbf{A})}, \qquad (2.47) $$

where h is a positive parameter that controls the weighted average of all pathways. The quantities $d_k(\mathbf{A})$ correspond to the minimum Euclidean distances in the strain space between $\mathbf{A}$ and the pathways k, defined by

$$ d_k(\mathbf{A}) = |^k \Pi(\mathbf{A})| = \min_{\zeta_k \in [0,1]} |\mathbf{A} - {}^k \hat{\mathbf{A}}(\zeta_k)|, \qquad (2.48) $$

where ${}^k \hat{\mathbf{A}}(\zeta_k)$ are also mapped onto the reaction pathways with $\zeta_k(\mathbf{A})$ the corresponding reaction coordinates. For example, when $\mathbf{A} = \mathbf{Ct}$: Fig. (2.3b) shows the projected tensor $^1\hat{\mathbf{C}}t(\zeta_1)$ onto the forward pathway 1, between the initial single-crystal bcc phase and the hcp variant $V_1$. Introducing the convenient curvilinear coordinates $\zeta_k^\infty(\mathbf{A})$ for fictitious unbounded pathways, as follows

$$ \zeta_k^\infty(\mathbf{A}) = {}^k \hat{\mathbf{D}} : (\mathbf{A} - {}^k \hat{\mathbf{C}}t_{\text{start}}) = \frac{{}^k \hat{\mathbf{C}}t_{\text{end}} - {}^k \hat{\mathbf{C}}t_{\text{start}}}{|{}^k \hat{\mathbf{C}}t_{\text{end}} - {}^k \hat{\mathbf{C}}t_{\text{start}}|} : (\mathbf{A} - {}^k \hat{\mathbf{C}}t_{\text{start}}), \qquad (2.49) $$

where ${}^k \hat{\mathbf{D}}$ defines the normalized direction of the pathway k, the argmin $\zeta_k$ in eq. (2.48) is also determined by solving $\partial_{\zeta_k} d_k(\mathbf{A}) = 0$ for a given $\mathbf{Ct}$, leading to

$$ \zeta_k(\mathbf{A}) = \begin{cases} \zeta_k^\infty(\mathbf{A}) & \text{if: } \zeta_k^\infty(\mathbf{A}) \in [0, 1] \\ 0 & \text{if: } \zeta_k^\infty < 0 \\ 1 & \text{if: } \zeta_k^\infty > 1, \end{cases} \qquad (2.50) $$

so that the distance measure $d_k(\mathbf{A})$ in eq. (2.48) with (2.50) represents the minimum distance from $\mathbf{A}$ to a given segment in $\mathbb{R}^6$.

On the other hand, it is assumed that each potential $\psi_{t_k}$ in eq. (2.46) is related to the minimum energy density $\hat{\psi}_{t_k}$ combining with an additional out-of-path component, i.e.,

$$ \psi_{t_k}(\mathbf{Ct}) = \hat{\psi}_{t_k}(\zeta_k(\mathbf{Ct})) + \underbrace{\sigma d_k(\mathbf{Ct}) + \pi |\text{tr}^k \Pi(\mathbf{Ct})|}_{\text{out-of-path component}}, \qquad (2.51) $$

such that $\partial_{\mathbf{Ct}} \text{tr}^k \Pi(\mathbf{Ct})$ and $^k \hat{\mathbf{D}}$ are orthogonal to each other, i.e., $\partial_{\mathbf{Ct}} \text{tr}^k \Pi(\mathbf{Ct}) : {}^k \hat{\mathbf{D}} = 0$. Here, $\text{tr} \mathbf{A}$ denotes the trace of $\mathbf{A}$. The parameters $\sigma$ and $\pi$ in eq. (2.51) scale two different out-of-path energy barriers: one component is linearly proportional to the Euclidean distance from the pathways with $\sigma$, while the second coefficient $\pi$ is used to distinguish different force magnitudes for isochoric and volumetric transformational deformations, when $\pi \neq 0$.

Figure (2.4) illustrates the construction of the overall inelastic energy landscape $\psi_t$ defined by eq. (2.46) with eq. (2.51), for all $\mathbf{Ct}$ of the neighborhood of the associated reaction pathway network in Fig. (2.3b). In accordance with the model parameters discussed in section 2.3.1, Fig. (2.4a) shows the given (invariant) minimum energy density $\hat{\psi}_{t_k}$ along all reaction coordinates $\zeta_k(\mathbf{Ct})$ of the individual pathways k. Then, the weighting functions $\omega_k(\mathbf{Ct})$ are used to extrapolate each contribution into the whole space: Fig. (2.4b) depicts a $5 \times 10^8$ J.m$^{-3}$-iso-surface of the extended inelastic energy part $\omega_k(\mathbf{Ct}) \hat{\psi}_{t_k}$ in the $\{\mathbf{C}_1, \mathbf{C}_2, \mathbf{C}_3\}$ strain space. As illustrated by arrows, the iso-surface is perpendicular to the reaction pathways and the energy profile is "sombrero-shaped" along the axis $\mathbf{C}_1 + \mathbf{C}_2 + \mathbf{C}_3$. Figure (2.4c) shows a $10^9$ J.m$^{-3}$-iso-volume related to the out-of-path contribution $\sigma d_k(\mathbf{Ct})$ only, i.e., with $\pi = 0$ in eq. (2.51). For sake of clarity, this additional energy potential is depicted in Fig. (2.4d) onto two planes passing by variants $V_1$ and $V_3$ (upper plane) and variants $V_5$ and $V_6$ (lower plane). It is also shown that the energy profile is exclusively controlled by the iso-distances around the paths, as illustrated by the cylinders around the paths and by the half-spheres at their ends. Finally, Fig. (2.4e) gives the same $10^9$ J.m$^{-3}$-iso-volume of
<figure><img src="image_6.png" /><figcaption>FIGURE 2.4: Construction of the total inelastic energy landscape $\psi_t$ associated with the multiple reaction pathways $V_k$ in iron. (a) Invariant and minimum energy profile along the individual reaction pathways $k$ from 0 (in dark red, for the pure bcc phases) to $\sim 8 \times 10^8 \text{ J.m}^{-3}$ (in white, for hcp phases). (b) Extrapolation of the minimum energy potential in the whole $\{C_1, C_2, C_3\}$ strain space, e.g. $5 \times 10^8 \text{ J.m}^{-3}$-iso-surface. (c) shows a $10^9 \text{ J.m}^{-3}$-iso-volume of the out-of-path contribution $\sigma d_k$ with $\pi = 0$, whereas (d) illustrates the energy profile onto two planes passing by variants $V_1$ and $V_3$ (upper plane) and variants $V_5$ and $V_6$ (lower plane). (e) and (f) are similar to (d) and (e) for the total inelastic energy $\psi_t$ landscape, respectively.</figcaption></figure>

the total inelastic energy $\psi_t$ landscape, within which the volume in (c) is plotted with transparency as well, for comparison. In contrast with Figs. (2.4c) and (d), it is shown that the total energy has a "cone-shaped" profile, exhibiting the directional character of the transformations toward the pure hcp phases, as distinctly depicted onto both planes in Fig. (2.4f).

### Transformational inelastic forces

The calculation of the inelastic driving forces for phase transformations in eq. (2.28) is deduced by computing the derivative of $\psi_t$ with respect to $\text{Ft}$, which can be expressed as follows

$$ \frac{\partial \psi_t}{\partial \text{Ft}} = 2 \text{Ft} \cdot \frac{\partial \psi_t (\text{Ct})}{\partial \text{Ct}}. \qquad (2.52) $$
According to eq. (2.46), the derivative of the energy function in the right-hand side of eq. (2.52) yields

$$ \frac{\partial \psi_t(\mathbf{Ct})}{\partial \mathbf{Ct}} = \sum_{k=1}^{18} \psi_{t_k}(\mathbf{Ct}) \frac{\partial \omega_k(\mathbf{Ct})}{\partial \mathbf{Ct}} + \omega_k(\mathbf{Ct}) \frac{\partial \psi_{t_k}(\mathbf{Ct})}{\partial \mathbf{Ct}}, \quad (2.53) $$

where the derivative of the weighting functions $\omega_k(\mathbf{A})$ with respect to $\mathbf{A}$ is given, without loss of generality, by

$$ \frac{\partial \omega_k(\mathbf{A})}{\partial \mathbf{A}} = h \sum_{i=1}^{18} \frac{\omega_i(\mathbf{A})}{d_i(\mathbf{A})} (\omega_k(\mathbf{A}) - \delta_{ik}) {}^i\mathbf{N}(\mathbf{A}), \quad (2.54) $$

with $\delta_{ik}$ the Kronecker delta, i.e., $\delta_{ik} = 1$ if $i = k$ and $= 0$, otherwise, and, ${}^i\mathbf{N}(\mathbf{A})$ represents the normal tensor to the pathway $i$ in the direction of $\mathbf{A}$, obtained in the following form:

$$ {}^i\mathbf{N}(\mathbf{A}) = \frac{\partial d_i(\mathbf{A})}{\partial \mathbf{A}} = \frac{i \mathbf{\Pi}(\mathbf{A})}{d_i(\mathbf{A})}, \quad (2.55) $$

such that $|{}^i\mathbf{N}(\mathbf{A})| = 1$, and, ${}^i\mathbf{N}(\mathbf{A}) : {}^i\hat{\mathbf{D}} = 0$ when $\zeta_k^\infty(\mathbf{A}) \in [0, 1]$. Moreover, the derivative of $\psi_{t_k}$ with respect to $\mathbf{Ct}$ in eq. (2.53) leads to

$$ \frac{\partial \psi_{t_k}(\mathbf{Ct})}{\partial \mathbf{Ct}} = \frac{\partial \hat{\psi}_{t_k}(\zeta_k(\mathbf{Ct}))}{\partial \zeta_k} {}^k\hat{\mathbf{D}} + \sigma {}^k\mathbf{N}(\mathbf{Ct}) + \pi \operatorname{sgn}(\operatorname{tr} {}^k\mathbf{\Pi}(\mathbf{Ct})) (\mathbf{I} - {}^k\hat{\mathbf{D}} \operatorname{tr} {}^k\hat{\mathbf{D}}). \quad (2.56) $$

Substituting eqs. (2.54) and (2.56) into eq. (2.53), and, then into eq. (2.52), it is also shown that two directions are included in the transformational inelastic forces: one component is related to the longitudinal directions ${}^k\hat{\mathbf{D}}$ along the reaction pathways, while the second component is associated with the normal directions ${}^k\mathbf{N}(\mathbf{Ct})$ towards $\mathbf{Ct}$.

## Mechanical elastic forces

Since the phase-field model aims at modeling high-pressure phase transitions in iron, particular attention is paid to the configuration within which the nonlinear elastic stiffness tensor is expressed. The out-of-path elasticity tensor $\mathcal{D}(\mathbf{Cet})$ in eq. (2.28), which depends on the elastic and transformational deformation states, is given in the whole strain space by

$$ \mathcal{D}(\mathbf{Cet}) = \sum_{k=1}^{18} \omega_k(\mathbf{Cet}) \, {}^k\mathcal{D}(\zeta_k(\mathbf{Cet})), \quad (2.57) $$

where ${}^k\mathcal{D}$ are the elasticity tensors associated with the reaction pathways $k$, and, $\zeta_k(\mathbf{Cet})$ are the reaction coordinates that minimize the Euclidean distance $d_k(\mathbf{Cet})$ between $\mathbf{Cet}$ and the individual paths $k$. The weighting functions $\omega_k(\mathbf{Cet})$ are also defined by eq. (2.47), where the partition of unity is written as a function of $\mathbf{Cet}$. The projected tensors ${}^k\hat{\mathbf{Cet}}$ are also mapped onto the reaction pathways (as well as ${}^k\hat{\mathbf{Ct}}$) and the corresponding reaction coordinates are consistently determined by solving $\partial_{\zeta_k} d_k(\mathbf{Cet}) = 0$. Imposing $\partial_{\mathbf{Cet}}\mathcal{D}(\mathbf{Cet}) = \mathbf{0}$ for all pure (meta)stable phases (i.e., at the ends of all reaction pathways $k$, when $\zeta_k = 0$ and $\zeta_k = 1$) with $\mathbf{Cet} = {}^k\hat{\mathbf{Cet}}$, the elasticity tensors ${}^k\mathcal{D}$ in eq. (2.57) may be represented by a cubic interpolation function to ensure numerical stability, i.e.,

$$ {}^k\mathcal{D}(\zeta_k(\mathbf{Cet})) = (1 - 3\zeta_k^2 + 2\zeta_k^3) \mathcal{D}^\alpha + (3\zeta_k^2 - 2\zeta_k^3) \mathcal{D}^\epsilon, \quad (2.58) $$

with $\mathcal{D}^\alpha$ and $\mathcal{D}^\epsilon$ the elastic stiffness tensors of the pure bcc and hcp iron phases, respectively. In particular, if $\zeta_k^\infty(\mathbf{Cet}) < 0$ ($> 1$), also ${}^k\mathcal{D}(\zeta_k(\mathbf{Cet})) = \mathcal{D}^\alpha (= \mathcal{D}^\epsilon)$. For instance, for such pure hcp
$\epsilon$-Fe phases, the finite hyperelasticity condition from eq. (2.18) leads therefore to

$$ \mathbb{D}^{\epsilon} = \rho_0 \frac{\partial^2 \psi_e}{\partial \mathbf{Ee} \partial \mathbf{Ee}} \Big|_{\epsilon} = \frac{\partial \mathbf{Se}}{\partial \mathbf{Ee}} \Big|_{\epsilon}, \qquad (2.59) $$

where $\mathbb{D}^{\epsilon}$ is defined in the reference configuration $\Omega_0$ and obeys the left and right minor symmetries. However, the elastic tensor differs from experimental or computed (e.g. using atomistic simulations) elasticity tensors $\mathbb{B}^{\epsilon}$, expressed in the current and deformed $\Omega$ by

$$ \mathbb{B}^{\epsilon} = \left. \frac{\partial \sigma}{\partial \epsilon} \right|_{\epsilon}, \qquad (2.60) $$

with $\epsilon$ the Eulerian strain tensor [269, 232, 66]. In the present work, the relevant tensor $\mathbb{B}^{\epsilon}$ for a pure hydrostatic compression is obtained by considering the following two-step deformation state along the reaction pathways: first, a material is subjected to a volumetric deformation $\mathbf{F}_{\text{vol}}$ from initial volume $V_0$ to the final volume $V = jV_0$ with $\mathbf{F}_{\text{vol}} = j^{1/3}\mathbf{I}$, and then, to a small and symmetric shear isochoric deformation $\mathbf{F}_{\text{iso}} = \mathbf{I} + \epsilon$ with $|\epsilon| \ll 1$, such as it is commonly performed using density functional theory calculations [163, 143]. Without plasticity, the total deformation gradient is also given by $\mathbf{F} = j^{1/3}(\mathbf{I} + \epsilon)$. Using eq. (2.60) with the aid of eq. (2.30) and considering $\sigma = p\mathbf{I}$ with the Cauchy pressure $p \le 0$, the relation between the elasticity tensors $\mathbb{D}^{\epsilon}$ and the incremental tangent modulus $\mathbb{B}^{\epsilon}$ is reduced to

$$ \mathbb{D}_{ijkl}^{\epsilon} = j^{-1/3} (\mathbb{B}_{ijkl}^{\epsilon} + p^{\epsilon} (\delta_{ij}\delta_{kl} - \delta_{ik}\delta_{jl} - \delta_{il}\delta_{jk})), \qquad (2.61) $$

exhibiting the same symmetries as in eq. (2.59). Thus, the elasticity tensors $\mathbb{D}^{\epsilon}$ is obtained by identifying the values of $\mathbb{B}^{\epsilon}$ as well as the corresponding equilibrium pressures $p^{\epsilon}$ from experiments or atomistic calculations. Inserting eq. (2.61) into eq. (2.58) with $\mathbb{D}^{\alpha} = \mathbb{B}^{\alpha}$ (here, the bcc $\alpha$-Fe phase is thermodynamically stable at zero pressure and zero temperature), and, then into eq. (2.57), the mechanical elastic driving forces in eq. (2.25) may therefore be determined in a computational Lagrangian framework.

### 2.2.6 Computational framework

The present model is implemented in a three-dimensional Lagrangian code using an element-free Galerkin least-squares formulation [25] with explicit time integration that handles acoustic wave propagation and rapid phase changes. The objective is to obtain solutions of the 12 unknown primary solution variables (i.e., degrees of freedom (DoFs) at integration nodes) namely, the displacement field $\mathbf{u}$ (3 DoFs) and the non-symmetric transformational distortion $\mathbf{Ft}$ (9 DoFs) at each reference point $\mathbf{X}$ in $\Omega_0$, by solving the system of partial differential equations, as follows

$$ \left\{
\begin{aligned}
\rho_0 \ddot{\mathbf{u}} &= \sum_{k=1}^{18} \nabla \cdot \left\{ \omega_k (\mathbf{Cet}) \mathbf{Fe} \cdot ({}^k\mathbb{D}(\zeta_k (\mathbf{Cet})) : \mathbf{Ee}) \mathbf{Ft}^{-t} \mathbf{Fp}^{-t} \right. \\
&\qquad \left. + \omega_k (\mathbf{Cet}) (\mathbf{Ee} : {}^k\mathbb{D}'(\zeta_k (\mathbf{Cet})) : \mathbf{Ee}) \mathbf{Fet} \cdot {}^k\mathbf{D} \cdot \mathbf{Fp}^{-t} + (\mathbf{Ee} : {}^k\mathbb{D}(\zeta_k (\mathbf{Cet})) : \mathbf{Ee}) \mathbf{Fet} \cdot \frac{\partial \omega_k (\mathbf{Cet})}{\partial \mathbf{Cet}} \cdot \mathbf{Fp}^{-t} \right\} \\
\nu \dot{\mathbf{Ft}} &= \sum_{k=1}^{18} \omega_k (\mathbf{Cet}) \mathbf{Ce} \cdot ({}^k\mathbb{D}(\zeta_k (\mathbf{Cet})) : \mathbf{Ee}) \mathbf{Ft}^{-t} + \lambda \nabla^2 \mathbf{Ft} \\
&\qquad - 2\rho_0 (\hat{\psi}_{tk}(\zeta_k (\mathbf{Ct})) + \sigma d_k (\mathbf{Ct}) + \pi |\text{tr} \, {}^k\mathbf{\Pi} (\mathbf{Ct})|) \mathbf{Ft} \cdot \frac{\partial \omega_k (\mathbf{Ct})}{\partial \mathbf{Ct}} \\
&\qquad - 2\rho_0 \omega_k (\mathbf{Ct}) (\hat{\psi}_{tk}'(\zeta_k (\mathbf{Ct}))) \mathbf{Ft} \cdot {}^k\mathbf{D} + \sigma \mathbf{Ft} \cdot {}^k\mathbf{N} (\mathbf{Ct}) + \pi \text{sgn}(\text{tr} \, {}^k\mathbf{\Pi} (\mathbf{Ct})) (\mathbf{Ft} - \mathbf{Ft} \cdot {}^k\mathbf{D} \text{tr} \, {}^k\mathbf{D})),
\end{aligned}
\right.
\qquad (2.62) $$
where ′ denotes the derivative with respect to the reaction coordinates $\zeta_k$, while the derivatives of the weighting functions $\omega_k$ with respect to Ct and Cet are determined by using eq. (2.54). The calculation of the first Piola-Kirchhoff stress tensor in the linear momentum balance in eqs. (2.62) is given in Appendix B from Ref. [252].

## 2.3 Pure hydrostatic compression

The phase-field formalism coupled with finite elastoplastic deformations is applied to analyze the $\alpha$-Fe into $\epsilon$-Fe phase transitions under high hydrostatic compression. The simulations exhibit the major role played by the plastic deformation in the morphological and microstructure evolution processes.

### 2.3.1 Material and model inputs

Tables (2) and (3) in Ref. [252] list the values for the material and model parameters for iron under high pressure compression, respectively, which have been collected from a variety of sources.

In the present phase-field model, the elastic pressure-dependent properties of iron are defined by four pressures: {$p^\alpha, p^\epsilon$}, for which the crystalline phases are fully bcc, and, fully converted to hcp, respectively; and: {$p^{\alpha \to \epsilon}, p^{\epsilon \to \alpha}$}, which characterize the transition states where the forward and reverse transformations start, respectively. Here, the equilibrium pressures {$p^\alpha = 0, p^\epsilon = -20$} GPa, with the corresponding atomic volumes {$v^\alpha = 11.75, v^\epsilon = 10.20$} Å³/at, are selected from Ref. [78]. In accordance with these experimental measures, the associated elastic components $B^\alpha$ and $B^\epsilon$ for both pure bcc and hcp phases are given in Ref. [163], while the stiffness tensor $D^\epsilon$ is expressed in the current configuration by using eq. (2.61), and, $D^\alpha = B^\alpha$ at zero pressure.

The ratio $c/a = 1.603$ of the hcp close-packed structure has been experimentally determined in Ref. [179], so that $\det U = 9c/a/16 = 0.902$. However, U corresponds to the complete phase transformation into the hcp iron sample at $p^\epsilon = -20$ GPa, for which the experimental measurements contain indistinctly elastic and transformational distortions. According to eq. (2.38) and following the procedure in Appendix C from Ref. [252], the transformational part Ut is related to U as follows

$$ Ut = \kappa U = \sqrt{2} \left( 1 + \sqrt{1 + \frac{8 j_{\text{exp}} p^{\epsilon}}{3 D^{\epsilon}}} \right)^{-1/2} U, \quad (2.63) $$

where $D^\epsilon$ is the hcp bulk modulus, and, $j_{\text{exp}} = v^\epsilon/v^\alpha$ is the experimental volume change from the initial pure bcc sample, at $p^\alpha = 0$ GPa, to the final pure polycrystalline hcp iron, at $p^\epsilon = -20$ GPa.

In the present perfect plasticity theory, a constant yield stress is chosen to analyze the crucial role of plasticity on nucleation and selection of variants during phase transformations, i.e., $\sigma_0 = 0.25$ GPa, which is fairly of the same order of magnitude with Hugoniot elastic limits in Ref. [214].

The positive parameter $h$ of the weighting functions controls the energetic part of the phase transition during a possible jump from one reaction pathway to the neighboring branches. The energy variation for such transition may be determined using molecular dynamics simulations [77], for which the exponent can be tuned to reproduce the atomistic results. However, without relevant information about the bcc-bcc and hcp-hcp phase transitions in iron, it is therefore assumed that all reaction pathways are mainly controlled by their immediate surroundings. This consideration may be achieved by imposing large magnitudes for $h$, e.g. $h = 10$, as well as large values for the energy barrier parameters $\sigma$ and $\pi$. The relation $\pi = 10\sigma$ (in GPa) is used in the energy penalty part of eq. (2.51) to consider higher pull-back forces onto the pathways for the volumetric than the isochoric phase transformations, which are conveniently applied to non-zero strain states that are out of the transition pathways, i.e., for any Ct with $^k\Pi(Ct) \neq 0$.

The onset of a new crystalline phase can be viewed as the product of a morphological instability involving elastic energy, interfacial energy, inelastic energy, transformational dissipation,
plastic dissipation, additional energies due to the long-range elastic interactions between variants, etc. Because of the complexity in modeling such phase instability, a phenomenological form is adopted to define the minimum energy density $\hat{\psi}_{tk}$ as a function of the reaction coordinate $\zeta_k$ along each individual pathway $k$, i.e.,

$$ \rho_0 \hat{\psi}_{tk} (\zeta_k (\mathbf{C}t)) = \frac{1}{2} c_1 \zeta_k^2 + c_2 \zeta_k, \qquad (2.64) $$

with $c_1$ and $c_2$ (in J.m$^{-3}$) two parameters that may be calibrated to experimental data. As described in Appendix D from Ref. [252], these parameters are given by

$$ c_2 = \frac{1}{2} j^{\alpha \to \epsilon} p^{\alpha \to \epsilon} \text{tr}^k \hat{\mathbf{D}}, \quad \text{and,} \quad c_1 = \frac{1}{2} j^{\epsilon \to \alpha} p^{\epsilon \to \alpha} ({}^k \mathbf{U} t^{-2} : {}^k \hat{\mathbf{D}}) - c_2, \qquad (2.65) $$

with $j^{\alpha \to \epsilon} = v^{\alpha \to \epsilon} / v^{\alpha}$ and $j^{\epsilon \to \alpha} = v^{\epsilon \to \alpha} / v^{\alpha}$ the experimental volume changes from the initial pure bcc sample to the Hugoniot states where the forward and reverse transitions occur, at $p^{\alpha \to \epsilon}$ and $p^{\epsilon \to \alpha}$, respectively. According to the recent experimental results from Ref. [78], the forward transition starts at $p^{\alpha \to \epsilon} = -14.9$ GPa, with the corresponding volume $v^{\alpha \to \epsilon} = 11.0$ Å$^3$/at, and, the reverse at $p^{\epsilon \to \alpha} = -12.0$ GPa, with $v^{\epsilon \to \alpha} = 10.6$ Å$^3$/at. The minimum energy density profile along the individual reaction pathways from eq. (2.64) with eq. (2.65), for which the values of $c_1$ and $c_2$ are provided in Tab. 3 from Ref. [252], is depicted in Fig. (2.4a).

The parameter $v$ in the relaxation eq. (2.23) is akin to viscosity in classical viscoplastic approaches. For the face-centered cubic (fcc) to bcc phase transitions in Fe$_3$Ni, an attempt to fit the magnitude $v = 14$ mPa.s, comparable to the viscosity of liquid metals, has been investigated by using molecular dynamics simulations [77]. Such quantitative data analysis is not available for the bcc-hcp transformations in iron, but it is assumed that the amount of stress state due to the viscosity is lower than the yield stress, i.e., $v\dot{\epsilon}_t < \sigma_0$, where $\dot{\epsilon}_t$ is a measure of the transformational strain rate. This measure can be estimated by $\dot{\epsilon}_t = \epsilon_t / \Delta t = \frac{1}{2} |\mathbf{C}t - \mathbf{I}| / \Delta t$ during a time interval $\Delta t$ awaited for the transformation, with $\epsilon_t$ the norm of the transformational Green-Lagrange deformation tensor. Thus, it follows that $v < \sigma_0 t_f / \epsilon_t$, with $t_f$ the final simulation time. According to the mentioned material inputs and time characteristics discussed in the following section, it is also considered that $v \approx \sigma_0 t_f / \epsilon_t \approx 2.6$ kPa.s.

Finally, the Laplacian operator in eq. (2.28) can be approximated using the mesh discretization in the finite element framework, such that $\lambda = \lambda^*/\ell^2$, where $\lambda^* = 0.5$ GPa is a mesh-size parameter and $\ell$ is an average element size of the simulation grid.

### 2.3.2 Analysis of the pressure-volume responses

The simulated material is a cube containing 1 million finite elements with full periodic boundary conditions, within which each element volume is $V_e = \ell^3 = 1$ µm$^3$. In the present dynamic continuum mechanics framework, the final simulation time $t_f$ is related to the physical time $t_c$, needed for acoustic waves to travel through the samples. Assuming that $t_f = 100 t_c$, the latter relation also means that the acoustic waves run over 100 times during the entire simulations for each sample, which ensures the quasi-static loading conditions. Thus, $t_c = L/c_L$, with $L = 100\ell = 0.1$ mm, the initial box length, and $c_L$, the longitudinal wave celerity in iron, i.e., $c_L = \sqrt{b_{11}^\alpha / \rho_0}$. It therefore follows that: $c_L = 5850$ m.s$^{-1}$, and, $t_f \approx 1.7$ µs, corresponding to the duration of the all performed calculations. Here and in the following, the subscript $f$ will denote the final state.

The initial single-crystal bcc iron is subjected to a three-step loading, as follows. First, all edges are continuously and proportionally decreased to a global volume reduction imposed by $j = V/V_0 = 0.86$, for which the volume change is achieved within a time step from $t_0$ to $t = 0.4$ $t_f$. Then, a constant volume is maintained from $t = 0.4$ $t_f$ to $0.6$ $t_f$, and, finally the volume is released back to the initial volume, so that $j_f = V_f/V_0 = 1$, at $t = t_f$.

Figure (2.5) illustrates the volume change $j$ as a function of the pressure $p$ in GPa. Although different in shape and magnitude, both hysteresis loops characterize martensitic transitions over
a wide range of pressure, involving an important stored elastic strain energy caused by the coexistence of numerous solid-state phases. The difference in both phase transformation hysteresis is due to plastic deformation in samples, which exhibits a larger width for the case with plasticity than without. When increased pressure, the appearance of the high pressure hcp phase is reached at -25.6 GPa, followed by a sudden drop to -23.1 GPa (without) and -19.7 GPa (with plasticity), due to dissipative effects during the forward $\alpha \rightarrow \epsilon$ transitions. However, the reverse $\epsilon \rightarrow \alpha$ transition without plasticity is characterized by a slow martensitic transformation, compared to an instantaneous volume change that occurs between -7.4 and -2.1 GPa with plasticity. Significantly, the forward transformation pressures predicted by the present model are higher than the experimental values for bcc samples that have been fully converted to hcp phases, within the range of -18.4 GPa [78] and -23.0 GPa [239] at room temperature. The experimental measurements from Refs. [100, 78] have been plotted in Fig. (2.5) with symbols, where the more recent data in Ref. [78] for high-purity Fe single crystals in helium pressure medium (shown by the oriented blue arrows) can be compared to the simulated hysteresis widths. Within the pressure range of coexistence of both phases, the experimental bcc (open symbols) and hcp (solid symbols) atomic volumes are separately deduced from X-ray diffraction measurements of lattice parameters at each applied pressure step. On the other hand, the computed results (solid lines) are obtained using the average pressure and volume states over the simulation samples. In addition, the pressure discrepancies are possibly due to the approximations/presumptions in the present coupled formalism and, more precisely, to the absence of free boundaries in the prescribed simulation setups. For instance, simulations in a helium pressure media, which is a fluid with a very low viscosity, together with a dislocation density-based crystal plasticity model, should give rise to a better description of the nonhydrostatic effects and anisotropic stresses in the transition pressures, and also of the hysteresis widths of iron. In accordance with the present calculations with periodic boundary conditions, classical molecular dynamics simulations using an embedded atom method potential have shown that the simulated transition pressure of the hcp and face-centered cubic (fcc) phases is significantly higher for uniform (31 – 33 GPa) than uniaxial (14 GPa) compression [270]. Although the simulated coexistence domain is larger than the experimental domain under quasi-hydrostatic conditions, the present P-V equation-of-state curves behave in good agreement with experimental responses when increasing (from 0 to -18 GPa) and decreasing (from -23 to -7 GPa) pressures [100, 78].

Figure (2.6) illustrates the partitioning of the total energy $\psi$ in terms of elastic $\psi_e/\psi$ (in blue) and inelastic ($\psi_t + \psi_\nabla)/\psi$ (green) energy ratios as a function of the dimensionless simulation time $t^* = t/t_f$, for calculations without and with plasticity. It is also shown that the total energy is mainly composed by the elastic strain energy until the nucleation of the first hcp phases in iron occurs at $t^* \approx 0.28$, as depicted by the two vertical arrows. When the volume is maintained constant, Fig. (2.6a) shows that the dissipative transformational process leads to 38% decrease in the amount of elastic energy, while the latter represents 54% of the total energy. During the early stages of the pressure release (as shown by a double-headed arrow), the stress state decreases, but the pressure remains sufficiently high to maintain the newly formed phases, as depicted by $*$ in Fig. (2.6a) when the internal elastic stored energy increases then to $t^* = 0.90$, before completely releasing back to zero. However, plastic deformation allows for a considerably higher stress relaxation between variants when phase transformations occur at large volume change states, as shown in Fig. (2.6b), where the upper thin curve for the elastic energy ratio without plasticity has been included for comparison. It also emphasizes the reduction of the stored elastic energy due to the plastic dissipation, for which the elastic strain energy falls down to 42% (compared to 54% without plasticity) of the total energy and remains constant during the second loading step. When the volume increases back to the initial volume, the elastic energy is then dramatically reduced to zero, significantly dissipated by plastic deformation.
<figure><img src="image_7.png" /><figcaption>FIGURE 2.5: Volume change *j* as a function of the pressure *p* in GPa, for calculations without (black dotted line) and with (red full line) plasticity, with $\sigma_0 = 0.25$ GPa. The experimental bcc (open symbols) and hcp (solid symbols) atomic volumes are separately deduced from X-ray diffraction measurements of lattice parameters at each applied pressure step, while the computed results (solid lines) are obtained using the average pressure and volume states over the simulation samples.</figcaption></figure>

### 2.3.3 Microstructure and variant selection

Figure (2.7a) illustrates the microstructure texture variation of transition-induced volume change *j* versus the dimensionless time *t*\* in the form of histograms. These histograms are obtained by splitting the simulation volume change (ranging from *j* = 0.80 to 1) into 100 bins of constant width, within which the phase fraction of materials is computed for all time steps. Coexistence of α-Fe and ε-Fe phases with different equilibrium volumes therefore leads to a multimodal histogram in the large range of pressure, where the grayscale represents the volume fractions of phases. For both simulations, the single-crystal volume is homogeneously decreased with respect to the prescribed hydrostatic conditions, as depicted by the points A. Without plasticity, Fig. (2.7a) shows a single-mode histogram: the volume change is slightly spread out over a large time interval, starting from the first forward phase transitions at *t*\* = 0.28 (point B). This spreading regime is spatially correlated to the strong elastic interactions between numerous variants that have partially been reversed into hcp phases only, from point B to D. However, continued pressure release results in a decrease in the proportion of the hcp phase compensated by an increase of the bcc phase between C and D. When the simulated iron is transformed back to the initial single-crystal material (point D), the volume exhibits no spatial variation, corresponding to a sharp single-mode histogram. With plasticity, the volume spreading is dramatically reduced after a brief fluctuation (point B) and remains a single mode until the first reverse phase transitions occur. Between *t*\* ≈ 0.75 and 0.90, a mixed-mode regime can be pointed out, which exhibits the structural texture formation of heterogeneous microstructure. The higher volume (point C') is greater in magnitude than the average prescribed volume, until all reversions are achieved (point D'). The second mode (point C) corresponds to a volume that remains constant and slightly increases during the reversions (point D). According to these different modes, a particular microstructure texture evolution in iron associated with preferential variant selection during the phase transitions is also expected.

Figure (2.7b) shows the volume fractions of each variant $V_k$ as a function of the simulation time *t*\*. Without plasticity, Fig. (2.7b) illustrates that the initial phase is partially transformed into the 6 possible hcp variants with comparable phase fractions, within which a residual amount of bcc phase persists in the microstructure, even for a large pressure range up to -25 GPa. When the compression is released to the original volume, all hcp variants are transformed back to the initial
<figure><img src="image_8.png" /><figcaption>FIGURE 2.6: Partitioning of the total energy into the elastic and inelastic components as a function of the dimensionless simulation time $t^*$, for calculations (a) without and (b) with plasticity.</figcaption></figure>

single-crystal bcc iron, behaving partially as a shape-memory alloy. For this case, most of transformations to $\epsilon$-Fe phases are partial only. These pseudo-hcp structures break the symmetries of fully formed hcp lattice, and, cannot lead to the formation of reversed $\alpha'$-Fe phases. Because the mismatch between bcc and hcp phases is not taken into account in the present formalism, the elastic strain state due to the interaction between variants is mainly responsible for the incomplete polymorphic phase transformations without plasticity. Therefore, when numerous hcp nucleus are considered, the long-range elastic interactions between variants dramatically increase the overall elastic energy, which in turn hinder the forward $\alpha \rightarrow \epsilon$ phase transitions. Because plasticity dissipates considerably the stored elastic strain energy, the onset of plasticity screens the elastic interactions between variants and thus decreases the energy cost to form the hcp variants. It also appears as an essential mechanism to enhance phase transformations by relaxing stresses due to elastic interactions, so that the complete formation of a polycrystalline iron formed by the 6 hcp variants is energetically favorable, as shown in Fig. (2.7b). In addition, a sudden burst of reversed $\alpha'$-Fe nucleation of variants occurs at $t^* \approx 0.90$, with ~2% volume fraction for each $\{V_{12}, V_{13}, V_{15}\}$, ~1% for each $\{V_{11}, V_{14}, V_{16}\}$, and, ~0.5% for each of the 6 other bcc variants. Thus, both initial $\alpha$-Fe and reversed $\alpha'$-Fe phases coexist at $t^* = 1.0$, without any retained hcp phases. However, the initial $\alpha$-Fe phase orientation largely dominates the forward and reverse transitions, while the volume fraction of $\alpha'$ inclusions is ~12.3% in the final microstructure.

To summarize, Fig. (2.8) illustrates the microstructure evolution under hydrostatic pressure at $t^* = 0.6$ and $t^* = 1.0$, defined in both strain and current mesh spaces. As shown in Fig. (2.8a), the non-flat sample surfaces capture the signature of the local unconstrained deviatoric stress component of the externally applied hydrostatic conditions. For the simulation without plasticity, the initial bcc $\alpha$-Fe phase (in gray) is not completely converted into hcp $\epsilon$-Fe phases, with a retained ~26.6% volume fraction of bcc phase at $t^* = 0.6$. However, the calculation with plasticity exhibits a polycrystalline iron that has been entirely transformed into 6 hcp $\epsilon$-Fe grain variants (red gradient). Such close-packed grains have been observed by performing large-scale molecular dynamics simulations under shock loading [137]. It is worth mentioning that various morphologies of hcp phases have been observed for structural phase transformations in iron, e.g. needle-like $\epsilon$-Fe phases [276], lath-like $\epsilon$-Fe regions [55], and, ellipsoidal $\epsilon$-Fe particles [202], for which the $\alpha \leftrightarrow \omega$ martensitic transitions exhibit similar discrepancies in zirconium [18]. On the release of hydrostatic pressure, the calculation without plasticity transforms back to the initial single-crystal bcc iron at $t^* = 1.0$, while the calculation with plasticity leads to 12 reversed bcc $\alpha'$-Fe, heterogeneously nucleated in pairs (e.g. $\{V_{11}, V_{12}\}$, in light and dark green) from one single $\epsilon$-Fe variant.
<figure><img src="image_9.png" /><figcaption>FIGURE 2.7: Evolution of (a) the volume change $j$ and (b) the phase fractions of variants $V_k$ as a function of the dimensionless simulation time $t^*$, for calculations without and with plasticity.</figcaption></figure>

## 2.4 Shock wave propagation

The numerical shock wave calculations accurately describe some important features reported by the experimental literature, and strongly complement our understanding of the phase-change dynamics in iron at larger time and length scales than hitherto explored by molecular dynamics simulations in the last two decades. The numerical model is able to reproduce unstable shock waves (which break up into elastic, plastic and phase-transition waves), providing new stress-informed insights into the coupling between the high strain-rate plasticity and microstructure evolution during the displacive phase transitions.

### 2.4.1 The internal structure of shock waves

In the following dynamical analyses, the three-dimensional iron samples are oriented along the [100] directions, and the shock waves are generated along the $z \parallel [001]_{bcc}$ direction, using 80 × 80 × 1280 element-free Galerkin nodes (~8.2 millions), with periodic boundary conditions transverse to the direction of shock front propagation, i.e., to $x \parallel [100]_{bcc}$ and $y \parallel [010]_{bcc}$. The initial shock compression is induced by imposing a velocity of 850 m.s$^{-1}$ on the rear face along $z \parallel [001]_{bcc}$, while the free surface is located at the extremity of the rectangular parallelepiped-shaped samples, as depicted in Fig. (2.9a). The unshocked material is at rest at $t = t_0 = 0$, while the final simulation time $t_f$ is related to the physical time $t_c$ for acoustic waves to travel through the sample. The dynamical loading conditions are controlled by assuming that $t_f = 2.5 t_c$, such that the acoustic waves run over 2.5 times the samples during the entire simulations. Thus, $t_c = L_z/c_L$, where $L_z$ is the initial box length in the [001]$_{bcc}$ shock direction, with $L_z = 16 L_x = 16 L_y = 1.28$ mm, and $c_L$ is the longitudinal wave celerity in iron, defined by $c_L = \sqrt{b_{11}^\alpha/\rho_0}$, with $b_{11}^\alpha = 271$ GPa the corresponding low-pressure elastic component of the pure bcc iron [163]. It therefore follows that: $c_L = 5850$ m.s$^{-1}$, so that $t_f \approx 0.55$ µs, which corresponds to the duration of all calculations. For convenience, a dimensionless time $t^*$ is defined as $t^* = t/t_f$, while the dimensionless length
<figure><img src="image_10.png" /></figure>

FIGURE 2.8: Transformational states defined in both strain and current mesh spaces at (a) $t^* = 0.6$ and (b) $t^* = 1.0$, for calculations without and with plasticity. Each black dot in the strain space represents the current transformational strain Ct for all mesh elements, while the colors along the pathways are associated with the corresponding phases and variants in the 3D simulated microstructures. Without plasticity, the initial bcc α-Fe phase remains in a large fraction (~26.6%, in gray) at $t^* = 0.6$, whereas the calculation with plasticity exhibits a polycrystalline iron formed by the 6 hcp e-Fe variants only (red gradient). On the release of hydrostatic pressure, the former is transformed back to the initial single-crystal bcc iron at $t^* = 1.0$, while the latter shows the presence of 12 reversed bcc α'-Fe with ~12.3% volume fraction.

$L^*$ along $z$ is given by $L^* = z/L_z$, so that both quantities $t^*$ and $L^*$ are ranged between 0 and 1. Moreover, the classical sign convention in continuum mechanics is used, so that compressive (extensive) volumetric stresses have negative (positive) signs.

The capability of the continuum element-free Galerkin model to reproduce the experimental multiple split-wave structure is illustrated in Fig. (2.9) by displaying the spatial heterogeneous distribution of the pressure behind the incident compressive wave. Figures (2.9c) and (2.9d) show the corresponding two- and three-wave structures for representative simulations without and with plasticity at $t^* = 0.35$, respectively. Different regions, namely, the initial unshocked, the elastically compressed bcc iron, and the transformed regions with high-pressure hcp Fe multivariants are also depicted. Furthermore, the plastically deformed bcc iron can be displayed for the calculation with plasticity in Fig. (2.9d). A sharp PT wave front is exhibited without plasticity, while a more complex rough PT front (see inset in Fig. (2.9d)) is shown to generate multiple planar pulses (as depicted by the vertical double-headed arrows) that propagate toward the leading plastic front. These localized traveling-wave fronts are suddenly produced by the dynamical phase transitions with high velocity in the compressed bcc region with high-pressure elastic properties. The consequences of the complex three-wave structure and competing wave interactions in the evolving
<figure><img src="image_11.png" /></figure>

FIGURE 2.9: (a) Schematics of the finite deformation framework that combines nonlinear elasto-viscoplasticity and multivariant phase-field theory to model the shock-induced response of single-crystal iron along the [001]$_{bcc}$ direction. (b) Distribution of the pressure resulting from three-dimensional simulation without plasticity. The unstable shock wave breaks up into the elastic precursor and the phase-transition wave, which leads to different internal deformation states at material points. (c) Similar calculation with plasticity, within which an intermediate plastic wave front propagates between the elastic and phase-transition wave fronts. The inset shows a rough phase-transition front, leaving behind a complex high-pressure microstructure with preferred selection and evolution of hcp variants.

deformation microstructure are elucidated in the following sections.

The shocked-induced microstructure during the martensitic phase transitions (also, the PT front) is analyzed in the six-dimensional Cauchy-Green strain space, as illustrated in Fig. (2.10). Thus, the deformation states that are mapped and visualized by colored points correspond to the local transformational distortions experienced by the iron samples. Each color is associated with the index of the nearest first-rank variant $V_k$. Figures (2.10a) and (2.10b) depict the corresponding states that are captured when the elastic fronts reach the free surfaces for calculations without plasticity and with plasticity, respectively. The former shows that two hcp variants are nucleated without plasticity, denoted by $V_1$ and $V_2$. These two preferential $\epsilon$-Fe variants are formed with different volume fractions, i.e., 62% for $V_1$ and 35% for $V_2$, and are thoroughly promoted by the [001]$_{bcc}$ direction of the shock. On the other hand, although the calculations with plasticity initiate the early formation of the same two variants, the four companion hcp variants are rapidly nucleated behind the PT wave front with comparable volume fractions. This microstructural fingerprint exhibits a crucial role played by the plastic deformation in nucleating and selecting all six energetically equivalent hcp variants in Fig. (2.10b). According to the previous simulations under high-pressure hydrostatic compression, the single-crystal iron has been transformed at high pressure into a polycrystalline microstructure that consists of the same six hcp variants, without any retained initial bcc phase. It is therefore suggested that the present high strain-rate plastic
<figure><img src="image_12.png" /></figure>

FIGURE 2.10: (a) The three-rank network of reaction pathways is projected in a {C₁, C₂, C₃} strain space, within which the local transformational Cauchy-Green Ct strain states at all material points are displayed with different colors (each color is associated with a specific hcp variant from the first-rank group symmetry operation). The results are related to the simulation without plasticity, captured at the instant when the elastic front reaches the free surface, revealing the nucleation of two (from amongst six possible variants) preferred hcp variants. (b) Similar simulation with plasticity at the same time instant as in (a), where the other four energetically equivalent hcp variants are activated in the transformed polymorphic microstructure. Such structural features indicate that the high strain-rate plastic deformation is locally capable of producing a nearly relaxed hydrostatic state from the uniaxial strain state produced by the shock-wave compression.

deformation can locally achieve a similar nearly relaxed three-dimensional hydrostatic state from
the uniaxial strain state produced by the shock-wave compression. The nucleation of all (also, six)
high-pressure hcp variants have never been described by atomistic calculations of shock-loaded
iron, certainly because of the small dimensions that hinder plastic relaxation needed to nucleate
these four companion hcp variants. For instance, two twinned hcp variants, separated by nonco-
herent grain boundaries (GBs), are observed in Refs. [137, 138].

**2.4.2 Effect of plasticity in shock-loaded iron**

Because the deformation processes act as distinctive signatures in shock-compressed samples, reflecting the history the solid experienced (in terms of velocity, shock pressure, etc.), three averaged quantities over the computational samples are plotted in Fig. (2.11). Slice-averaged quantities within spatial planar bins (of one element width) are also used to quantify the role of plasticity in tailoring the complex microstructure from the uniaxial strain deformation, namely the free-surface velocity $v_z$ in Fig. (2.11a), the pressure $p = -(\sigma_{xx} + \sigma_{yy} + \sigma_{zz})/3$ in Fig. (2.11b), and the von Mises stress $\sigma_{\text{VM}}$ in Fig. (2.11c) with respect to $t^*$, obtained without (gray curves) and with (black curves) plasticity. Both averaged quantities $p$ and $\sigma_{\text{VM}}$ are displayed with respect to $L^*$ along the $z \parallel [001]_{\text{bcc}}$ loading direction of the samples.

Figure (2.11a) shows the presence of two distinct plateaus for the free-surface velocity profile without plasticity (gray curve), supporting by the split two-wave structure into the noticeable fastest elastic and the phase-transition (denoted by PT, see arrow) waves. The elastic wave is characterized by the elastic precursor Ep with $v_z = 255 \text{ m.s}^{-1}$, while the phase-transition front produces a considerable increase of the velocity at free surface, up to $v_z = 1660 \text{ m.s}^{-1}$. On the other hand, the simulation with plasticity shows a much more complex velocity profile, where the multiple-wave structure consists inter alia of the elastic precursor Ep with the same velocity as the case without plasticity, the plastic (P wave) front, and the elastic wave reverberation with the on-going PT wave, i.e., the rEp wave. This wave profile is comparable to those reported in experimental works with distinct three-wave structures [19, 130]. The instants when both P and
<figure><img src="image_13.png" /></figure>

FIGURE 2.11: (a) Free surface velocity histories from shock-loaded iron samples without (gray curve) and with (black curve) plasticity. The former is caused by the arrival of the elastic precursor (denoted by Ep) and of the phase-transition (PT) wave. The latter is decomposed by Ep, the plastic (P) wave front, and rEp that results from the interaction between the reflected Ep front at the free surface and the on-coming PT wave. (b) The representative profiles of pressure in GPa along the [001]bcc direction for both calculations without and with plasticity. The slice-averaged values within spatial planar bins of one finite element width correspond to the three-dimensional microstructures in Figs. (2.9c) and (2.9d). (c) The von Mises stress in GPa for both calculations without and with plasticity.

rEp waves reach the free surface are displayed by the double-headed arrows in Fig. (2.11a), corresponding to $v_z = 880 \text{ m.s}^{-1}$ and $v_z = 1170 \text{ m.s}^{-1}$, respectively. It is worth noting that both reflected Ep and P waves that propagate back in the elastically compressed and plastically deformed bcc iron (thus, along the [001]bcc direction) produce a residual stress state that does not favor the mandatory forward $\alpha \rightarrow \epsilon$ phase transitions. The interaction in releasing the stresses between the reflected Ep and P waves with the PT wave encourages therefore the reverse $\epsilon \rightarrow \alpha$ phase transitions, without retaining any $\epsilon$-Fe hcp phase nor without forming any $\alpha'$-Fe bcc variants. Interestingly, this feature differs from the pure hydrostatic compression loading, for which a significant residual volume fraction ($\sim 12\%$) of $\alpha'$ bcc inclusions has been obtained in the microstructure after the reverse phase transformations. Consequently, the incident PT wave cannot reach the free surface for calculations with plasticity, in contrast to the simulation case without plasticity. Additionally, it is worth mentioning that the amplitude of the steady-state free-surface velocity with plasticity is close to the one without plasticity, i.e., $v_z = 1707 \text{ m.s}^{-1}$, which is roughly twice the particle velocity of $850 \text{ m.s}^{-1}$ imposed on the rear face behind the incident shock as a loading condition, consistently with the traction-free conditions at free surfaces.

Both calculations without and with plasticity in Fig. (2.11b) exhibit a similar elastic state where
compression remains uniaxial in the [001]$_{bcc}$ direction, characterized by a pressure $p_E = 3.9$ GPa in the elastically compressed bcc phase. By considering this threshold pressure as the Hugoniot elastic limit for the plasticity-free case, the value of 3.9 GPa is defined between two reference experimental data in polycrystalline iron samples, i.e., ~2.1 GPa [292] and ~5.5 GPa [230]. It is worth mentioning that the similar computed values for both uniaxial elastic limits without and with plasticity are fortuitous since the former corresponds to the transformational front (accompanied by both hydrostatic and deviatoric stresses), while the latter is related to the plastic front (mainly controlled by deviatoric stresses). In practice, once the phase transformation operated by one specific variant is initiated, the excess free energy between both bcc and hcp iron phase promotes a partially-to-complete shock-induced transition that behaves differently than pure pressure, as quantified by eq. (2.37). The corresponding released stress state after this phase transformation is much more complex than the stress state behind the deviatoric stress-driven plastic front. The changes from the uniaxial shock compression to a complex stress state after phase transitions in the plastically deformed iron cannot therefore be captured by a usual pressure-shock velocity (e.g. represented by a Rayleigh line), yielding an important distinction between the shock physics described at the macroscopic scale and ones described at the grain scale. Behind the traveling Ep wave front, the pressure profile depicts the presence of one (two) plateaus for calculations without (with) plasticity. The former exhibits the presence of the PT wave front as the pressure dramatically increases up to $p_{PT} = 37.7$ GPa. The latter profile shows an intermediate pressure plateau that characterizes the plastically deformed bcc region, within which the forward $\alpha \rightarrow \epsilon$ phase transitions start roughly at the onset pressure $p_{PT} = 18.2$ GPa, as indicated by the dotted line in Fig. (2.11b). This value is on the range of experimental values for single-crystal iron under hydrostatic pressure [78], and in excellent agreement with large-scale molecular dynamics simulations in single-crystal iron as well, i.e., 18 GPa along the same [001]$_{bcc}$ shock direction [274]. Here, the value deviates from the conventional macroscopic threshold from experiments on polycrystalline Fe samples (occurring at 13 GPa [17, 19]), for which the GBs with pre-existing intrinsic defects reduce the amplitude of the forward transition pressure [105, 275, 297]. Achieved after the complete phase transformation of the bcc into hcp variants, the upper plateau is governed by the load intensity and is reached at $p = 44.1$ GPa, slightly higher than the pressure without plasticity. This value is in very good agreement with recent results from molecular dynamics simulations in shocked iron [4], where a maximum mean pressure of ~40 GPa has been measured by applying a comparable piston velocity of $800 \text{ m.s}^{-1}$.

Figure (2.11c) shows the corresponding values for the von Mises stress, with $\sigma_{vM} = 2.7$ GPa for both simulations in the elastically compressed bcc iron. Then, the large von Mises stress profile increases inhomogeneously in the sample without plasticity, which is due to a heterogeneous distribution of both hcp variants $V_1$ and $V_2$ in the microstructure with lamellar arrangements along the shock direction (not shown here). The maximum value is $\sigma_{vM} = 18.1$ GPa. With plasticity, however, the volume-preserving plastic deformation relaxes significantly the internal von Mises stress to reach an averaged von Mises stress of $\sigma_{vM} = 1.1$ GPa ($< 3.9$ GPa, at the peak Hugoniot elastic state) in the shocked-induced hcp multivariant region. This difference asserts the role played by plasticity to release the shear stress state produced by the uniaxial strain compression to obtain a roughly hydrostatic state with 6 high-pressure hcp variants (instead of 2 variants without plasticity) in the transformed heterogeneous microstructure.

Figures (2.12a) and (2.12b) capture the evolution of the longitudinal stress component in the shock direction $\sigma_{zz}$ in the Lagrangian adimensional position-time ($L^*, t^*$) diagrams, without and with plasticity, respectively. The non-steady-state regimes of the present elastic precursor (Ep, solid lines), plastic (P, dashed), and phase-transition (PT, dotted) waves —moving with different average speed so that net distances between the respective fronts increase with time— exhibit a more complicated picture for the three-wave structure with the high strain-rate plasticity than the corresponding diagram without plasticity. The reflection of the incident fronts from the free surfaces are depicted as well.
<figure><img src="image_14.png" /></figure>

FIGURE 2.12: (a) Slice-averaged maps of the longitudinal stress component $\sigma_{zz}$ in the Lagrangian adimensional position-time ($L^*, t^*$) diagram from simulation without plasticity. The two-wave (composed of the elastic Ep and phase-transition PT waves) structure with the reflection of both waves at the free surface are shown using different line types. (b) The three-wave structure with the presence of the intermediate plastic wave (P wave) front illustrates a considerably more complicated scenario of nonlinear wave interaction. As depicted by the arrows, this calculation with plasticity reveals two nucleation events at $t^* = 0.10$ and $t^* = 0.27$, which result in the inhomogeneous propagation of the trailing PT wave and in the presence of a stress-release envelope. The latter travels faster than the leading shock and is characterized by a lower longitudinal stress in magnitude.

The leading E wave front, traveling at $5412 \text{ m.s}^{-1}$ ($5541 \text{ m.s}^{-1}$) for calculation with (without) plasticity, leaves the iron system in an elastically compressed state with high-pressure properties. The former value is in excellent agreement with the computed shock velocity of $5409 \text{ m.s}^{-1}$ using atomistics simulations in single-crystal iron without pre-existing defects [137], which is consistent with the present calculations. Without plasticity, the trailing PT front travels homogeneously in the sample at $4655 \text{ m.s}^{-1}$. For the three-wave structure, the nearly over-driven P front (but not over-run, i.e., characterized by a finite separation between the E and P waves) propagates at $5059 \text{ m.s}^{-1}$, while the slower heterogeneous PT front travels with intermittent regimes at $3002 \pm 99 \text{ m.s}^{-1}$, which is much lower than the homogeneous PT front without plasticity. In contrast to the case without plasticity, the intermittent propagation of the PT front with plasticity reveals the presence of i) sudden nucleation events of hcp variants (as depicted by the arrows in Fig. (2.12b)), and consequently of ii) a so-called traveling release-stress envelope. This envelope propagates by reflection between the rear surface on the left-hand side of the sample and the PT wave before interacting the (unloading) reflected Ep wave with the free surface, as displayed by the asterisk \* in Fig. (2.12b). It precedes always the slower wave, i.e., the PT wave, but travels faster than the elastic wave at $8312 \text{ m.s}^{-1}$ in the transformed high-pressure regions of iron (i.e., with high pressure-induced stiffness and density). These distinct nucleation sites of hcp variants are not experienced for calculations without plasticity, exhibiting again the specific role played by the plastic deformation in governing such microstructural features. Analogous distinct nucleation events in position-time diagrams have been observed in shocked crystalline 1,3,5-triamino-2,4,6-trinitrobenzene using large values for the input parameter $\sigma$ in molecular dynamics simulations [152].

### 2.4.3 Residual stresses in the plastically-deformed microstructure

Figures (2.13a) display three shock-induced microstructures $M_1$, $M_2$, and $M_3$ in Fig. (2.12b) that are associated with $t^* = 0.21$, $t^* = 0.27$, and $t^* = 0.44$, respectively, for the calculation with plasticity only. For these microstructures, various stress-related quantities, i.e., the longitudinal Cauchy
<figure><img src="image_15.png" /><figcaption>FIGURE 2.13: (a) Three-dimensional time snapshots of shock-induced microstructures, designated M₁ at t* = 0.21, (b) M₂ at t* = 0.27, and (c) M₃ at t* = 0.44, from the position-time diagram in Fig. (2.12b) for simulation with plasticity. From top to the bottom, each panel captures the heterogeneous distribution of various stress quantities, namely, the longitudinal Cauchy stress σ<sub>zz</sub>, the shear stress τ, the stress-related quantities s<sub>n</sub> and s<sub>s</sub> to the second invariant J<sub>2</sub> using eq. (2.66), as well as the polycrystalline high-pressure domains composed of six hcp variants. These variants are colored using the same code as in Fig. (2.10b), while the transparent zones are associated with the initial unshocked bcc iron. As displayed by *, a dynamic instability in the polymorphic phase transitions is observed in M₂, leading to the nucleation of a large monovariant with columnar growth in the microstructure that is still visible (**) after the propagation of the incident phase-transition wave front. (b) The color legends associated with the stress-related quantities.</figcaption></figure>

stress tensor component in the shock direction σ<sub>zz</sub>, the shear stress τ = (σ<sub>zz</sub> - (σ<sub>xx</sub> + σ<sub>yy</sub>) / 2) / 2,
s<sub>n</sub>, and s<sub>s</sub>, as well as the corresponding hcp variant selection, are displayed. Both stress quantities
s<sub>n</sub> and s<sub>s</sub> are related to the second invariant of the stress deviator J<sub>2</sub> and the von Mises stress σ<sub>vM</sub>
by

$$
3J_2 = \sigma_{\text{vM}}^2 = \frac{3}{2} \operatorname{dev} \sigma : \operatorname{dev} \sigma = \frac{1}{2}s_n + 6s_s, \quad (2.66)
$$

where dev σ is the deviatoric part of σ, so that s_n and s_s are defined by

$$
\begin{align}
s_n &= (\sigma_{xx} - \sigma_{yy})^2 + (\sigma_{yy} - \sigma_{zz})^2 + (\sigma_{xx} - \sigma_{zz})^2 \nonumber \\
s_s &= \sigma_{xy}^2 + \sigma_{yz}^2 + \sigma_{xz}^2, \tag{2.67}
\end{align}
$$

with σ<sub>xy</sub>, σ<sub>xz</sub>, and σ<sub>yz</sub> being the orthogonal shear stresses. As a signed quantity, the shear stress τ, which equals to the von Mises stress if the off-diagonal terms are neglected, can also have positive (in red) or negative (green) values depending on the magnitude of σ<sub>zz</sub> with respect to (σ<sub>xx</sub> + σ<sub>yy</sub>) / 2. All color legends for the stress-related quantities are displayed in Fig. (2.13b).

At instant t* = 0.21, the split three-wave structure into the Ep, P, and PT wave fronts is clearly
distinguishable by the change in magnitude of σ<sub>zz</sub> in Fig. (2.13a). Close to the phase-transition
front, the transformed region with 6 high-pressure hcp variants is characterized by positive val-
ues of the shear stress τ (values in red). Between the PT and P wave fronts, the shear stress τ
is negative (green), the stress field s<sub>s</sub> is zero, while the quantity s<sub>n</sub> exhibits the presence planar
surfaces as pulses generated by the PT front that dynamically nucleates the hcp variants. These
six variants are pictured with the same colors as in Fig. (2.10b). Behind the complex rough PT
front, some hcp grains grow preferentially into flaky morphology with (110)<sub>bcc</sub> and (11̄0)<sub>bcc</sub> habit
<figure><img src="image_16.png" /><figcaption>FIGURE 2.14: (a) Slice-averaged magnitude of the plastic Green-Lagrange strain tensor Ep in the Lagrangian adimensional position-time ($L^*, t^*$) diagram from simulation with plasticity. The elastic and plastic wave fronts with constant velocities are shown using different line types as well as two (primary and secondary) phase-transition zones that are associated with specific nucleation of release and reload variants (see text for details). (b) The corresponding selection of hcp variants, categorized into two groups, so-called $G_1 = \{V_1 + V_2\}$ (in blue) and $G_2 = \{V_3 + V_4 + V_5 + V_6\}$ (red).</figcaption></figure>

planes of the bcc iron, which are transformed into the $(0001)_{hcp}$ close-packed planes after the phase transition.

At $t^* = 0.27$, the presence of a dynamical instability in the compressed and plastically deformed microstructure is shown. This occurs under a complex stress state that is responsible to an extremely rapid nucleation of a large single-crystal hcp variant $V_1$ (in orange, as depicted by * in $M_2$ in Fig. (2.13a)) with columnar growth in the direction of the shock loading. This spontaneous nucleation is characterized by a notable change in sign of the shear stress $\tau$ from negative (green) to positive (red) values. The ideal volume-reducing transition path of the strain-free monovalent variant $V_1$ requires a compression of ~12.5% along the $z \parallel [001]_{bcc}$ direction, as defined by the $zz$ component in eq. (2.37). This sudden nucleation event gives rise to the aforementioned traveling release-stress envelop in Fig. (2.12b), which is also characterized by a finite domain with positive shear stress values, as depicted by white double-sided arrows in $M_2$ and $M_3$. Surrounded by the initial bcc phase, the variant $V_1$ is able to grow in the shock direction, whereas the confined region between the PT front and $V_1$ in $M_2$ becomes an unstable zone for nucleation of high-pressure variants. Similar shock-driven regions of instabilities, within which local nucleation of hcp embryos occur, have been observed by Wang et al. using atomistic simulations [273].

Although $V_1$ is still visible at instant $t^* = 0.44$, the phase-transition wave front continues to propagate in the shock direction, exhibiting the coalescence of the hcp variants and also a specific morphological fingerprint of shock-induced hcp variants with large transformed bands (due to the periodic boundary conditions) at high pressure. A thickness of ~77 µm for $V_1$ is found in the $z \parallel [001]_{bcc}$ direction, which also depends on the shock velocity (results not shown here). Overall, $s_n$ exhibits large values in the elastically compressed zones, which significantly decrease as soon as the nucleation of growth of hcp variants take place during the polymorphic phase transitions. In turn, because both quantities $s_n$ and $s_s$ quantities play a complementary role in the present $J_2$ plasticity theory, $s_s$ gives rise to large values in the phase-transformed hcp regions. The aforementioned transformed bands are therefore considered here as an important mechanism of stress relaxation under shock compression at high strain rate, thus providing novel guidelines for future experimental diagnostics of shock wave propagation in iron.
<figure><img src="image_17.png" /><figcaption>FIGURE 2.15: Schematics of the presence of internal dislocation structures at solid-solid interfaces.</figcaption></figure>

### 2.4.4 Dynamical instability in structural phase transitions

Figure (2.14a) illustrates the shock-induced instability in the structural phase transitions by means of the magnitude of the plastic Green-Lagrange deformation Ep, defined by $|Ep| = |Fp^t \cdot Fp - I|/2$. This quantity is plotted in the Lagrangian adimensional position-time ($L^*, t^*$) diagram, where $t^*$ is restricted between 0 and 0.5 for clarity, so that the multiple reflections of incident waves from the free surface are conveniently omitted in the following discussion. It is shown that the propagation of the PT front gives rise to a spatially (not temporally) heterogeneous distribution of $|Ep|$ with local values up to 0.25. This localization of plastic deformation is therefore strongly correlated with the specific selection of shock-induced hcp variants, which can be separated into two pertinent groups, so-called $G_1 = \{V_1 + V_2\}$ and $G_2 = \{V_3 + V_4 + V_5 + V_6\}$, each set involving different features of the microstructural fingerprints in shock-loaded iron. Thus, Fig. (2.14b) displays the variant selection during the shock wave propagation using a linear interpolation of color to distinguish the presence of both groups $G_1$ (blue) and $G_2$ (red) in the microstructure. As already mentioned, both $V_1$ and $V_2$ variants (from amongst six available variants) of $G_1$ are promoted by the shock direction in the first instants of the shock wave propagation. Since the two-phase mixture induces a large contraction along the loading $z \parallel [001]_{bcc}$ direction, the corresponding group $G_1$ is composed of variants designated by "release variants". However, the second group $G_2$, which consists of a mixture of the complementary four variants with identical volume fractions, experiences an expansive reaction in the shock direction. In contrast to $G_1$, these newly-formed variants of $G_2$ are also expected to generate an expansion (or reloading) wave, which are therefore not promoted by the initial compressive (or loading) wave. In the following, the four variants of $G_2$ are denoted by "reload variants", for which the nucleation is accompanied by severe plastic deformation with large values of $|Ep|$, as indicated by Fig. (2.14a).

## 2.5 Limitations

While the present phase-field approach is capable of considering the elastic mismatch between low- and high-pressure variants during the pure pressure- and shock-induced phase transformations in iron, the coexistence of both solid-state phases with different crystal structures (e.g. lattice parameters) yields to the loss of lattice coherence at the interfaces. This also means that the
perfect lattice correspondence across the bcc/hcp interfaces as well as the misoriented hcp/hcp
grain boundaries obtained in Fig. (2.8) becomes an implausible model assumption, and that the
current description of the crystalline interfaces during solid-solid phase transitions remains obvi-
ously incomplete. In fact, experimental observations of such interfaces show that coherent inter-
faces break down through the formations of misfit dislocation structures, as sketched in Fig. (2.15)
with internal hexagonal dislocation patterns. The resulting "semc coherent interfaces" consist of
coherent regions separated by these interfacial dislocation structures. Since the earliest observa-
tions of dislocation arrangements into periodic patterns along solid-state interfaces in a variety
of conditions [6, 54, 5, 65], the advantages/inconveniences introduced by the presence of such
crystal defects in high-technology applications have been addressed in interdisciplinary materials
science and engineering [237, 95], involving chemistry, physics, electronics, metallurgy, mechan-
ics, etc. Extensive investigations have indicated that the interfacial dislocation patterns at grain
and interphase boundaries may, however, be designed to increase the unprecedented levels of
high strength [7], ductility [300], and radiation-induced damage tolerances [26] in nanocrystalline
polycrystals, nanolayered laminated composites, precipitation-strengthened alloys, and epitaxial
free-standing thin films. In part, the fundamental problem of characterizing the dislocation struc-
tures and energetics at heterophase interfaces is treated in the following chapter 3.
Chapter 3

Dislocation structures and energetics at
heterophase interfaces

Selected peer-reviewed articles

[P28] A. Vattré, E. Pan. *Dislocation singularities in layered magneto-electro-elastic plates*. International Journal of Engineering Science, 181, 103765, 2022.

[P21] A. Vattré, E. Pan. *Semicoherent heterophase interfaces with core-spreading dislocation structures in magneto-electro-elastic multilayers under external surface loads*. Journal of the Mechanics and Physics of Solids, 124, 929-956, 2019.

[P20] A. Vattré, N. Abdolrahim, S. Navale, M. Demkowicz. *The relaxed structure of intrinsic dislocation networks in semicoherent interfaces: predictions from anisotropic elasticity theory and comparison with atomistic simulations*. Extreme Mechanics Letters, 28, 50-57, 2019.

[P19] T. Jourdan, **A. Vattré.** *A continuous model including elastodiffusion for sink strength calculation of interfaces*. Computational Materials Science, 153, 473-478, 2018.

[P18] **A. Vattré**, E. Pan. *Three-dimensional interaction and movements of various dislocations in anisotropic bicrystals with semicoherent interfaces*. Journal of the Mechanics and Physics of Solids, 116, 185-216, 2018.

[P17] **A. Vattré**, E. Pan. *Interaction between semicoherent interfaces and Volterra-type dislocations in dissimilar anisotropic materials*. Journal of Materials Research, 32, 3947-3957, 2017.

[P16] **A. Vattré**. *Elastic strain relaxation in interfacial dislocation patterns: II. From long- and short-range interactions to local reactions*. Journal of the Mechanics and Physics of Solids, 105, 283-305, 2017.

[P15] **A. Vattré**. *Elastic strain relaxation in interfacial dislocation patterns: I. A parametric energy-based framework*. Journal of the Mechanics and Physics of Solids, 105, 254-282, 2017.

[P14] **A. Vattré**, T. Jourdan, H. Ding, C. Marinica, M. Demkowicz. *Non-random walk diffusion en-hances the sink strength of semicoherent interfaces*. Nature communications, 7, 1-10, 2016.

[P13] **A. Vattré**. *Elastic interactions between interface dislocations and internal stresses in finite-thickness nanolayered materials*. Acta Materialia, 114, 184-197, 2016.

[P10] **A. Vattré**. *Mechanical interactions between semicoherent heterophase interfaces and free surfaces in crystalline bilayers*. Acta Materialia, 93, 46-59, 2015.

[P9] **A. Vattré**, M. Demkowicz. *Partitioning of elastic distortions at a semicoherent heterophase interface between anisotropic crystals*. Acta Materialia, 82, 234-243, 2015.
[P8] **A. Vattré**, N. Abdolrahim, K. Kolluri, M. Demkowicz. *Computational design of patterned interfaces using reduced order models*. Scientific report, Nature, 4, 2014.

[P7] **A. Vattré**, M. Demkowicz. *Effect of interface dislocation Burgers vectors on elastic fields in anisotropic bicrystals*. Computational Materials Science, 88, 110-115, 2014.

[P5] **A. Vattré**, M. Demkowicz. *Determining the Burgers vectors and elastic strain energies of interface dislocation arrays using anisotropic elasticity theory*. Acta Materialia, 61, 5172-5187, 2013.

### 3.1 Motivation

Far from being featureless dividing surfaces between neighboring crystals, interfaces in homo- and hetero-phase solids have internal structures of their own. These structures depend on in- terface crystallographic character (misorientation and interface plane orientation) and affect the physical properties of interfaces, such as interface energy [73], resistivity [42], permeability [125], mechanical properties [141], morphology and variant selection of precipitates [210], point defect sink efficiencies [225], and mobilities [147]. To better understand and control the properties of interfaces, it is desirable to be able to predict their internal structures. The first part of this chap- ter 3 presents a method for predicting a specific interface structural feature: the Burgers vectors of intrinsic dislocations in semicoherent homophase and heterophase interfaces. This informa- tion is then used to compute the interface elastic strain energies in standard tilt and twist GBs as well as the partition of elastic distortions at complex heterophase interfaces. An application to the sink strength of semicoherent interfaces is described in section 3.5, for which the non-random walk diffusion of radiation-induced defects is highly sensitive to the detailed character of inter- facial stresses. The follow-up extensions to the elastic strain relaxation in interfacial dislocation patterns and to the elastic interaction with extrinsic dislocation arrays and loops are investigated in sections 3.6 and 3.7, respectively.

One way of studying interface structure is through atomistic simulations, which explicitly ac- count for all the atoms that make up an interface. However, this approach is not always practical or efficient: it can be very resource-intensive because it requires a separate simulation for each in- dividual interface. Thus, it does not lend itself to rapidly scanning over many different interfaces, for example if one were searching for trends in interface structures or for tailored interfaces with a specific structure. Low-cost, analytical techniques for predicting interface structure would be preferable in such situations.

One widely used analytical approach applies to semicoherent interfaces and describes inter- face structures in terms of intrinsic dislocations using the closely related Frank-Bilby [93, 30] and O-lattice [32, 296, 237] techniques. Both procedures require the selection of a reference state, within which the Burgers vectors of individual interface dislocations are defined. Because this choice does not affect the calculated spacing and line directions of interface dislocations, it has some- times been viewed as if it were arbitrary. In practice, one of the adjacent crystals [145, 110, 289] or a "median lattice" [91] have often been used as the reference state.

However, the choice of reference state does influence the values of far-field stresses, strains, and rotations associated with interface dislocations. These, in turn, are usually subject to con- straints, namely that the far-field stresses be zero and that the far-field rotations be consistent with a prescribed misorientation. Thus, the choice of reference state is in fact not arbitrary. As discussed by Hirth and co-workers [119, 122, 123], the importance of selecting proper reference states has often been overlooked in part because the best-known applications of interface dislo- cation models are to interfaces of relatively high symmetry, such as symmetric tilt or twist GBs, for which correct reference states are easy to guess. Furthermore, many analyses assume uniform isotropic elasticity, which leads to equal partitioning of interface dislocation elastic fields between the neighboring crystals. In general, however, interfaces need not have high symmetry and the
neighboring crystals may have unlike, anisotropic elastic constants. By use of heterogeneous and anisotropic elasticity theory, the correct selection of reference states in such general cases is far more challenging.

Elasticity theory for analyzing semicoherent interfaces and determining the field solutions pro-
duced by interface dislocations has been initiated by van der Merwe [244]. The concept of misfit
dislocations, which act as stress annihilators to free the total stress fields far from the interfaces,
has been introduced using the Peierls-Nabarro model to formulate a misfit dislocation theory for
critical thicknesses of strained and layer systems during epitaxial growth of structures with two
isotropic semi-infinite solids [245, 243]. The problem of single straight screw and edge disloca-
tions and dislocation arrays situated at the interface between two anisotropic elastic half-spaces
has received special attention in the literature [282, 60, 21, 34, 283, 277, 260, 146], for which the
dislocation-based calculations and also mechanisms may be significantly altered when isotropic
elastic approximation is considered.

By means of the Stroh sextic formalism [234, 235] with a Fourier series-based technique, the geometry of interface dislocation patterns as well as the corresponding Burgers vectors have been solved using anisotropic elasticity theory in bicrystals with two sets of dislocations [260, 258, 261]. This computational method for structural and energetic properties of individual heterophase interfaces has been extended by taking into account the presence of free surfaces in bi- and tri-layered materials [247, 248] and the local reactions between planar and crossing dislocation arrays to form new dislocation arrangements [249, 250]. Application examples have revealed the significant influence played by elastic anisotropy in the interactions between the semicoherent interfaces and radiation-induced point defects [255] as well as extrinsinc dislocation loops [256].

**3.2 Determining the Burgers vectors of interface dislocation arrays**

The notion of introducing Volterra dislocations into a reference state for constrained interfaces is
consistently defined with the Frank-Bilby equation that are free of far-field stresses.

**3.2.1 Planar interfaces in linear elastic bicrystals**

In the present analysis, planar interfaces are considered formed by joining two semi-infinite linear elastic crystals, for which the crystallography of the interfaces has been specified completely. For a GB, this requires five parameters: three to describe the relative misorientation between neighbor-ing crystals and two to describe the orientation of the GB plane [237]. For a heterophase interface, the number of crystallographic DoFs may be higher. For example, an interface between two fcc crystals such as Al and Ni would require the lattice parameters of the two neighboring metals to be given in addition to the five parameters needed for a GB. Interfaces between materials with differing crystal structures may require further parameters.

To describe completely the crystallography of a heterophase interface between elements A and B, the notion of a "reference" state for the interface is adopted: in the reference state, the interface is coherent, i.e. the two separate crystals that meet at the interface are rotated and strained [131, 237] such that they are in perfect registry with each other across the interface plane after bonding. Thus, the reference state has the interface structure of a single perfect crystal.

Starting from the reference state, materials A and B are mapped separately into new configura-
tions that yield an interface with the required crystallographic character and zero far-field stresses,
as shown in Fig. (3.1). Following Hirth, Pond, and co-workers [123], the state of the interface after
this mapping is referred as the "natural" state. For a GB, the maps applied to materials A and B are
proper rotations while for a pure misfit interface they are pure strains. To account for both cases
as well as for heterophase interfaces between misoriented crystals, the maps are described as uni-
form displacement gradients $_{A}$F and $_{B}$F. In the reference state, the neighboring crystals might not
<figure><img src="image_18.png" /><figcaption>FIGURE 3.1: Mapping from a coherent reference state to the natural state using displacement gradients $_{A}$F and $_{B}$F. Volterra dislocations introduced into the reference state remove coherency stresses and may change the misorientation of the neighboring crystals.</figcaption></figure>

be stress free, but the interface is coherent. In the natural state, the interface is not coherent, but
the neighboring crystals are both free of far-field stresses.

This framework is sufficiently general to describe the crystallography of many commonly
studied heterophase interfaces, e.g. ones formed by fcc and bcc metals [73, 75], but not all. For
example, mapping from a common reference state to an interface between a cubic and hcp crystal
cannot directly be accomplished by a displacement gradient alone and requires an internal shuffle
rearrangement, as mentioned in section 2.2.5. The present chapter 3 is also focused on materials
that may be mapped to a common reference state using displacement gradients alone.

The crystallographic considerations described above do not require a single, unique reference state. On the contrary, an infinite number of new reference states may be generated from an original one by applying to it any uniform displacement gradient $_{R}$F. If the original reference state may be mapped to the natural state with $_{A}$F and $_{B}$F, then the new reference state may be mapped to the same natural state using $_{A}$F$_{R}$ F$^{-1}$ and $_{B}$F$_{R}$ F$^{-1}$. However, a consistent description of the elastic fields of a discrete dislocation network in an interface of specified crystallography and free of far-field stresses does require a single specific reference state.

3.2.2 Volterra dislocations in the reference state

The atomic structures of real interfaces are not like those generated by the linear mappings from a
reference state. Instead, for any given interface crystallography, the atomic structure may undergo
a variety of local relaxations or reconstructions that lower its energy. In many low-misorientation
GBs and low-misfit heterophase interfaces, these changes lead to formation of regions of co-
herency (which generally have low energies) separated by networks of intrinsic dislocations.
Many such interface dislocation networks have been imaged using transmission electron mi-
croscopy [5].

There are two common ways of describing interface dislocations. In one, they are viewed not as conventional Volterra dislocations, but rather as special kinds of interface defects with short-range elastic fields that are formed when the interface atomic structure in the natural state relaxes [120, 37]. The superimposed elastic fields of all such defects residing within an interface decay away to zero at long range and therefore do not alter the far-field stress state or the crystallography of the natural interface state.

Another description—the one adopted here—views interface dislocations as genuine Volterra dislocations with resultant elastic stress fields that need not decay to zero at long range. For example, the structure of some pure misfit heterophase interfaces may be described as an array of equally spaced edge dislocations residing on the same glide plane [181]. Such an array of Volterra
dislocations has a non-zero far-field stress [7]. Certain symmetric tilt GBs may be described as arrays of edge dislocations lying directly one above the other on separate glide planes. These Volterra dislocation arrays have zero far-field strains (hence, also zero stresses [7]), but possess non-zero rotations at long range [213, 164]. In general, arrays of Volterra dislocations may have non-zero far-field strains, rotations, or both.

In the work described here, interface dislocations are viewed as Volterra dislocations that have been introduced into the reference state, as shown in Fig. (3.1). Therefore, the far-field stresses due to these dislocations $A\sigma_{dis}^{\infty}$ and $B\sigma_{dis}^{\infty}$ are equal and opposite to the coherency stresses $A\sigma_{c}$ and $B\sigma_{c}$ in the reference state respectively, leading to the removal of all far-field stresses in the natural state:

$$A\sigma_{c} + A\sigma_{dis}^{\infty} = 0, \quad \text{and,} \quad B\sigma_{c} + B\sigma_{dis}^{\infty} = 0. \qquad (3.1)$$

Although free of long-range stresses, interface dislocation networks in the natural state have non-zero short-range elastic fields as a result of the superposition of the non-uniform stress fields of the Volterra dislocation networks and the uniform coherency stresses in the reference state. Additionally, the far-field rotations due to the Volterra dislocations are required to conform to the given interface crystallographic character. These requirements restrict the choice of reference states to a single specific one.

The notion of introducing Volterra dislocations into the reference state primarily is treated as a hypothetical operation. However, this operation may be a physically meaningful analog of processes occurring at some real interfaces. For example, the transformation of certain coherent heterophase interfaces into ones that are not coherent, but free of far-field stresses, occurs by the deposition on the interface of Volterra dislocations that glide through the neighboring crystalline layers [181, 182]. Similarly, subgrain boundaries are thought to assemble from glide dislocations formed during plastic deformation of polycrystals [8].

### 3.2.3 Crystallographic constraints on interface dislocations

A variety of shapes of interface dislocation networks have been observed [5], such that the ones that may be represented by $j \le 2$ arrays of parallel dislocations with Burgers vectors $b_j$, line directions $\xi_j$, and inter-dislocation spacings $d_j$. Following previous investigators [93, 30, 237], these quantities are related to the density of admissible Volterra dislocations in the reference state and interface crystallography as

$$B = \sum_{i=1}^{j} \left( \frac{n \times \xi_i}{d_i} \cdot p \right) b_i = (A F^{-1} - B F^{-1}) p = T p, \qquad (3.2)$$

where $n$ is a unit vector normal to the interface and the so-called probe vector $p$ is any vector contained within the interface plane. Equation (3.2) is known as the quantized Frank-Bilby equation [237, 289], where $T$ corresponds to an average operation that maps $p$ to the resultant Burgers vector $B$ of interface dislocations intersected by $p$.

The individual Burgers vectors $b_i$ of interface dislocations are assumed to be related to the crystal structure of the reference state. For example, if the reference state is an fcc crystal of lattice parameter $a$, values of $b_i$ may be drawn from a set of $a/2 \langle 110 \rangle$-type glide or $a/6 \langle 112 \rangle$-type Shockley partial dislocation Burgers vectors. Once the set of admissible Burgers vectors is known, well-studied methods stemming from Bollmann's O-lattice theory [32] may be used to compute $n$, $\xi_i$, and $d_i$ [145, 289] from the O-lattice vectors $p_i^o$, defined by

$$b_i = T p_i^o. \qquad (3.3)$$

The O-lattice vectors $p_i^o$ and therefore both $\xi_i$ and $d_i$ do not depend on the choice of reference state. If an original reference state is mapped to a new one using displacement gradient $\mathbf{r}_F$, then $b_i$
is mapped to $\tilde{b}_i = {}_R F b_i$. Here and in the following, the superimposed inverse caret will be used to indicate trial values of variables. The new reference state may also be mapped to the natural state using ${}_A \tilde{F} = {}_A F {}_R F^{-1}$ and ${}_B \tilde{F} = {}_B F {}_R F^{-1}$, as discussed in section 3.2.1. Assuming that rank $\mathbf{T} = 3$, the O-lattice vectors computed from the original and new reference states are identical:

$$ \mathbf{p}_i^o = \mathbf{T}^{-1} \mathbf{b}_i = ({}_A \tilde{\mathbf{F}}^{-1} - {}_B \tilde{\mathbf{F}}^{-1})^{-1} \tilde{\mathbf{b}}_i = \tilde{\mathbf{p}}_i^o. \quad (3.4) $$

This conclusion may also be shown for matrix $\mathbf{T}$ of rank 2 or 1. Thus, for a given set of Burgers vectors $b_i$, interface crystallography uniquely determines interface dislocation line directions $\zeta_i$ and spacings $d_i$, but not the reference state. Based on this result, some authors have argued that the choice of reference state is truly arbitrary [32]. However, in different reference states, $b_i$ will clearly have different magnitudes and directions, both of which influence the magnitudes of the elastic fields generated by interface dislocations (the latter by altering their characters).

### 3.2.4 Solution strategy

Determining the elastic fields of semicoherent interfaces requires finding the correct interface dislocation Burgers vectors, which are defined in the coherent reference state. The following five-step strategy is applied to determine the specific reference state that meets the constraints of interface crystallographic character and zero far-field stresses.

#### Step 1: Solving for geometry of dislocation networks

As shown in section 3.2.3, the geometry of interface dislocations (their line directions and spacings) is independent of the choice of reference state. Thus, a reference state is chosen identical to one of the crystals adjacent to the interface in its natural state. This choice provides an initial guess of the interface dislocation Burgers vectors. Then, the interface dislocation geometry is determined by using standard methods [31, 145, 110]. Multiple dislocation geometries are possible in some interfaces, but attention is restricted in this section to interfaces with unique geometries.

#### Step 2: Solving for interface dislocation elastic fields

The complete elastic fields, produced by the arrays of dislocations found in step 1, are determined using anisotropic linear elasticity theory in bicrystals. The elastic fields are assumed to follow the periodicity of the two-dimensional dislocation structures predicted in step 1 and must also satisfy specific boundary conditions at the interfaces.

#### Step 3: Solving for far-field distortions

The far-field distortions associated with each set of parallel dislocations are computed separately and then superimposed to obtain the resultant far-field distortions of the interface dislocation network as a whole. These elastic distortions are key for determining the correct reference state for the interfaces of interest. Far-field strains, stresses, and rotations may also be deduced.

#### Step 4: Solving for the reference state

The correct reference state is the one in which the superposition of the strains produced by interface dislocation arrays eliminate the coherency strains, giving a bicrystal that is free of far-field stresses and has far-field rotations that agree with the given interface crystallographic character. This condition is met by continuously adjusting the reference state along a specified transformation pathway, starting with the initial guess selected in step 1.
<figure><img src="image_19.png" /><figcaption>FIGURE 3.2: (a) Schematic illustration of a planar interface dislocation network formed by bonding materials A and B. (b) The geometry of an interface containing two sets of dislocations described by O-lattice vectors $\mathbf{p}_1^0$ and $\mathbf{p}_2^0$. Open circles represent O-lattice points and filled circles illustrate atoms with nearly matching positions in materials A and B.</figcaption></figure>

## Step 5: Solving for the interface elastic strain energy

Incomplete cancellation of the coherency and Volterra fields near the interface gives rise to short-range stresses and strains. These stresses and strains are used to compute the elastic energies of semicoherent interfaces.

### 3.2.5 Elastic fields of interface dislocation arrays

This section is focused on interfaces containing up to two arrays of infinitely long straight, and uniformly spaced parallel dislocations at equilibrium, as illustrated in Fig. (3.2a). The Stroh formalism of anisotropic linear elasticity [234, 235, 59] and a Fourier series-based solution technique are used to compute the elastic fields outside the cores of interface dislocations [22, 68, 34]. For clarity in this section, the pre-subscripts A and B in the field expressions will be omitted if no distinction between materials is required.

#### Problem formulation

The geometry of a dislocation network consisting of two arrays of straight parallel dislocations may be described by two O-lattice vectors $\mathbf{p}_1^0 \neq \mathbf{p}_2^0$ in the interface of interest using a Cartesian coordinate system with basis vectors ($x_1, x_2, x_3$), as shown in Fig. (3.2b). An interface containing only one array of straight parallel dislocations is a special case of this more general geometrical description. The unit vector normal to the interface is $\mathbf{n} \parallel x_2$, with the interface located at $x_2 = 0$: $x_2 > 0$ for material A, and $x_2 < 0$ for material B. The dislocation line direction $\xi_1$ is parallel to $\mathbf{p}_2^0$ and $\xi_2 \parallel \mathbf{p}_1^0$, as illustrated in previous studies [110, 237, 289].

A representative interface unit cell of the dislocation pattern is illustrated in Fig. (3.2b). Translations of the unit cell by the basis vectors $\mathbf{p}_1^0$ and $\mathbf{p}_2^0$ tessellate the interface plane. It is also convenient to identify a non-orthogonal (oblique) frame with basis vectors ($x'_1, x'_2, x'_3$), where $x'_1 \parallel \mathbf{p}_1^0 \parallel \xi_2$ and $x'_3 \parallel \mathbf{p}_2^0 \parallel \xi_1$. The oriented angle between $\xi_2$ and $\xi_1$ is denoted by $\phi$, so that $x'_1 = x_1 \csc \phi$ and $x'_3 = x_3 - x_1 \operatorname{ctg} \phi$. Thus, any position vector in this non-orthogonal frame may be expressed as $\mathbf{r} = x'_1 \mathbf{p}_1^0 + x'_3 \mathbf{p}_2^0$.
Due to the periodicity of the interface dislocation structure, it is useful to seek a complete set of wavevectors $k$ such that the elastic fields in the interface may be analyzed using plane waves $e^{i2\pi k \cdot r}$. The set of all $k$ is conveniently written as $k = n p_1^\times + m p_2^\times$ with respect to the reciprocal vectors $p_1^\times$ and $p_2^\times$, defined by the orthogonality conditions $p_\alpha^\times \cdot p_\beta^\times = \delta_{\alpha\beta}$, where $n, m$ are integers.

The complete elastic distortion field $\mathbf{D}$ is the superposition of the uniform coherency and the Volterra dislocation distortions, $\mathbf{D}_c$ and $\mathbf{D}_{\text{dis}}$, as discussed in section 3.2.2. Following the seminal work of Bonnet [34, 35, 36], outside of dislocation cores, $\mathbf{D}$ may be expressed as the biperiodic Fourier series, i.e.

$$ \mathbf{D}(\mathbf{x}) = \mathbf{D}_c + \mathbf{D}_{\text{dis}}(\mathbf{x}) = \mathbf{D}_c + \text{Re} \sum_{k \neq 0} e^{i2\pi k \cdot r} \mathbf{D}^k(x_2), \quad (3.5) $$

with $i = \sqrt{-1}$, while Re stands for the real part of a complex quantity and the sum spans over all non-zero wavevectors $k$. The Fourier amplitudes of the complete distortion waves $\mathbf{D}^k(x_2)$ are required to converge (not necessary to zero) in the far-field, i.e. $x_2 \to \pm\infty$. The components $k_1$ and $k_3$ of the wavevector $k$ satisfy

$$ k \cdot r = k_1 x_1 + k_3 x_3 = \left( \frac{n \csc \phi}{|\mathbf{p}_1^0|} - \frac{m \operatorname{ctg} \phi}{|\mathbf{p}_2^0|} \right) x_1 + \frac{m}{|\mathbf{p}_2^0|} x_3. \quad (3.6) $$

The complete displacement field $u$ may be found by integrating eq. (3.5) as

$$ u(\mathbf{x}) = \underbrace{u_0 + \mathbf{D}_c \mathbf{x}}_{\text{affine part}} + \text{Re} \sum_{k \neq 0} e^{i2\pi k \cdot r} u^k(x_2) = u_{\text{aff}}(\mathbf{x}) + u_{\text{dis}}(\mathbf{x}), \quad (3.7) $$

where $u_0$ is an arbitrary constant displacement. The complete displacement field $u$ may be decomposed into an affine part $u_{\text{aff}}$ corresponding to $\mathbf{D}_c$ and a biperiodic Fourier series representation of displacement fields $u_{\text{dis}}$ generated by the Volterra dislocations.

The Fourier amplitudes in eqs. (3.5) and (3.7) are determined from linear elasticity in the absence of body forces and subject to boundary conditions associated with interface dislocations. The complete displacement gradients $\mathbf{D}(\mathbf{x}) = \text{grad } u(\mathbf{x})$ in crystals A and B must fulfill the partial differential equations of mechanical equilibrium

$$ \operatorname{div}(\mathbf{C} : \operatorname{grad} u(\mathbf{x})) = 0, \quad (3.8) $$

where $:$ denotes the double inner product and $\mathbf{C}$ is a fourth-order anisotropic elasticity tensor.

## Complete field solutions

Substituting the displacement field eq. (3.7) into eq. (3.8), the second-order differential equation applied to both half-spaces is obtained as follows

$$ w_1 \mathbf{W}_1 u^k(x_2) + w_2 (\mathbf{W}_2 + \mathbf{W}_2^t) \frac{\partial u^k(x_2)}{\partial x_2} + \mathbf{W}_3 \frac{\partial^2 u^k(x_2)}{\partial x_2^2} = 0. \quad (3.9) $$

with $w_1 = -4\pi^2$ and $w_2 = i2\pi$. Here, $^t$ denotes the matrix transpose and $\mathbf{W}_1, \mathbf{W}_2$, and $\mathbf{W}_3$ are 3 × 3 real matrices related to the wavevectors (i.e. interface geometry) and the stiffness constants
(i.e. elasticity) indexed in Voigt notation:

$$
\begin{align}
\mathbf{W}_1 = \mathbf{W}_1^\mathrm{t} &= \begin{bmatrix}
k_1^2 c_{11} + 2k_1 k_3 c_{15} + k_3^2 c_{55} & k_1^2 c_{16} + k_1 k_3 (c_{14} + c_{56}) + k_3^2 c_{45} & k_1^2 c_{15} + k_1 k_3 (c_{13} + c_{55}) + k_3^2 c_{35} \\
& \ddots & \ddots \\
& & k_1^2 c_{56} + k_1 k_3 (c_{36} + c_{45}) + k_3^2 c_{34} & k_1^2 c_{66} + 2k_1 k_3 c_{46} + k_3^2 c_{44} \\
& & & k_1^2 c_{55} + 2k_1 k_3 c_{35} + k_3^2 c_{33}
\end{bmatrix} \tag{3.10} \\
\mathbf{W}_2 &= \begin{bmatrix}
k_1 c_{16} + k_3 c_{56} & k_1 c_{12} + k_3 c_{25} & k_1 c_{14} + k_3 c_{45} \\
k_1 c_{66} + k_3 c_{46} & k_1 c_{26} + k_3 c_{24} & k_1 c_{46} + k_3 c_{44} \\
k_1 c_{56} + k_3 c_{36} & k_1 c_{25} + k_3 c_{23} & k_1 c_{45} + k_3 c_{34}
\end{bmatrix} \nonumber \\
\mathbf{W}_3 &= \mathbf{W}_3^\mathrm{t} = \begin{bmatrix}
c_{66} & c_{26} & c_{46} \\
& c_{22} & c_{24} \\
& & c_{44}
\end{bmatrix} . \nonumber
\end{align}
$$

As demonstrated in Appendix A from Ref. [260], the complete displacement field (3.7) is written
as follows

$$
u(\mathbf{x}) = u_0 + \mathbf{D}_c \mathbf{x} + \mathrm{Re} \frac{1}{i2\pi} \sum_{k \neq 0} e^{i2\pi k \cdot \mathbf{r}} \sum_{\alpha=1}^{3} \lambda^{\alpha} e^{i2\pi p^{\alpha} x_{2}} \mathbf{a}^{\alpha} + \zeta^{\alpha} e^{i2\pi p_{*}^{\alpha} x_{2}} \mathbf{a}_{*}^{\alpha}, \quad (3.11)
$$

where the eigenvalues $p^\alpha$ and eigenvectors $a^\alpha$ are calculated by solving the sextic algebraic equa-
tion of the Stroh formalism [234, 235] for each material A and B. The asterisk indicates complex
conjugates of solutions with positive imaginary parts, i.e. $p_*^{\alpha+3} = p_*^{\alpha}$ and $a_*^{\alpha+3} = a_*^{\alpha}$, indexed by
$\alpha = 1, 2, 3$. The complete elastic strains and stresses are also deduced from eq. (3.11) by

$$
\begin{align}
\mathbf{E}(\mathbf{x}) &= \{\mathbf{D}(\mathbf{x})\} = \frac{1}{2} (\operatorname{grad} \mathbf{u}(\mathbf{x}) + \operatorname{grad} \mathbf{u}^{\mathrm{t}}(\mathbf{x})) \nonumber \\
\sigma(\mathbf{x}) &= \mathbf{C} : \mathbf{E}(\mathbf{x}), \tag{3.12}
\end{align}
$$

respectively. Equation (3.12a) gives the strain-displacement relationship, where {$\mathbf{D}(\mathbf{x})$} denotes
the symmetric component of the distortion field, while eq. (3.12b) is the Hooke’s law for small
strains that determines the stress field. The general solutions of elastic fields of eqs. (3.11–3.12)
are expressed as linear combinations of the eigenfunctions given by eq. (3.76), and include $\lambda^\alpha$
and $\zeta^\alpha$ as complex unknown quantities that are to be determined by the boundary conditions, as
follows.

**Boundary condition 1: Convergence of far-field solutions**

In accordance with Saint Venant's principle, the convergence of the Fourier amplitudes $u^k(x_2)$ when $x_2 \to \pm\infty$ leads to the requirement that $\zeta_5^\alpha = 0$ and $\lambda^\alpha = 0$. This condition applies to infinite bicrystals and would not be appropriate for bicrystals terminated with free-surfaces.

**Boundary condition 2: Absence of far-field strains**

The elimination of the coherency strains $E_c$ by the far-field strains of the interface Volterra dislocations $E_{\text{dis}}^{\infty}$ is taken into account by requiring the total elastic strain field $\mathbf{E}$ to decay to zero when $x_2 \to \pm\infty$, i.e.

$$
\lim_{x_2 \to \pm\infty} \mathbf{E}(\mathbf{x}) = \mathbf{E}^{\infty} = \mathbf{E}_c + \mathbf{E}_{\text{dis}}^{\infty} = \mathbf{0}, \quad (3.13)
$$

where $\mathbf{E}_c = \{\mathbf{D}_c\}$ and $\mathbf{E}_{\text{dis}}^{\infty} = \{\mathbf{D}_{\text{dis}}^{\infty}\}$ is the far-field strain produced by the interface dislocations.
Equation (3.13) is equivalent to eqs. (3.1) expressed using strains rather than stresses. As detailed
in Appendix B from Ref. [260], the far-field distortions, calculated individually for each set of
dislocations, $i = 1$ and 2, and then superposed, are given as follows

$$
\mathbf{D}_{\text{dis}}^{\infty} = -\operatorname{sgn}(x_2) \operatorname{Re} \sum_{i=1}^{2} d_i^{-1} \sum_{\alpha=1}^{3} \bar{\lambda}_i^{\alpha} \mathbf{G}_i^{\alpha} + \bar{\zeta}_i^{\alpha} \mathbf{G}_{i*}^{\alpha}. \quad (3.14)
$$
Here, $A\bar{\zeta}_1^\alpha = A\bar{\zeta}_2^\alpha = 0$ and $B\bar{\lambda}_1^\alpha = B\bar{\lambda}_2^\alpha = 0$ for the reasons described in boundary condition 1. Superimposed bars are used to indicate quantities related to the far-field boundary conditions, while the complex constants $A\bar{\lambda}_i^\alpha$ and $B\bar{\zeta}_i^\alpha$ are determined by solving a specific system of equations, as described in Ref. [260].

### Boundary condition 3: Disregistry due to interface Volterra dislocations

Disregistry is the discontinuity of displacements across an interface [7], expressed in terms of the relative displacements between neighboring atomic planes. Each dislocation produces a stepwise change in disregistry at its core with magnitude equals its Burgers vector. The disregistry at $x_2 = 0$ of a network of two sets of dislocations may be represented by the staircase functions

$$ \Delta \mathbf{u}(x_1, x_3) = A \mathbf{u}(x_1, x_3) - B \mathbf{u}(x_1, x_3) = -b_1 \left[ \frac{\csc \phi x_1}{|\mathbf{p}_1^0|} \right] - b_2 \left[ \frac{x_3 - \operatorname{ctg} \phi x_1}{|\mathbf{p}_2^0|} \right], \quad (3.15) $$

as illustrated in Fig. (3.3), where only one set has been displayed for clarity, for which the complete displacement discontinuity at the interface can therefore be expressed as

$$ \Delta \mathbf{u}(x_1, x_3) = \Delta \mathbf{u}_{\text{aff}}(x_1, x_3) + \Delta \mathbf{u}_{\text{dis}}(x_1, x_3). \quad (3.16) $$

The left-hand side of eq. (3.16) gives the relative displacement field $\Delta \mathbf{u}_{\text{aff}}$ at the interface generated by the uniform macroscopic distortions $A\mathbf{D}_c$ and $B\mathbf{D}_c$ in the affine form

$$ \Delta \mathbf{u}_{\text{aff}}(x_1, x_3) = \Delta \mathbf{u}_0 + \left[ (A\mathbf{D}_c - B\mathbf{D}_c) \mathbf{x} \right]_{x_2=0}, \quad (3.17) $$

where $\Delta \mathbf{u}_0 = -\frac{1}{2}(\mathbf{b}_1 + \mathbf{b}_2)$ is chosen, without loss of generality. As shown in Fig. (3.3), eq. (3.17) may be interpreted as a continuous distribution of (fictitious) Volterra dislocations with infinitesimal Burgers vectors and spacing [30, 195].

The right-hand side of eq. (3.16) is the displacement discontinuity $\Delta \mathbf{u}_{\text{dis}}$ produced by equilibrium interface dislocations in the natural state, shown as $\Delta$ in Fig. (3.1). According to eqs. (3.7) and (3.11), the quantity $\Delta \mathbf{u}_{\text{dis}}$ is given in Ref. [260]

$$ \Delta \mathbf{u}_{\text{dis}}(x_1, x_3) = \frac{1}{i2\pi} \sum_{k \neq 0} e^{i2\pi k \cdot r} \sum_{\alpha=1}^{3} A \lambda_{A}^{\alpha} \mathbf{a}^{\alpha} - B \zeta_{B}^{\alpha} \mathbf{a}_{*}^{\alpha}, \quad (3.18) $$

which may be represented by sawtooth functions [87, 34, 90], as illustrated in Fig. (3.3). Using the Fourier sine series analysis and superposing the sawtooth-shaped functions associated with the two sets of dislocations, eq. (3.18) can be expressed as

$$ \Delta \mathbf{u}_{\text{dis}}(x_1, x_3) = \underbrace{\sum_{n=1}^{\infty} -\frac{b_1}{n\pi} \sin 2\pi n \frac{\csc \phi x_1}{|\mathbf{p}_1^0|}}_{\text{set 1}} + \underbrace{\sum_{m=1}^{\infty} -\frac{b_2}{m\pi} \sin 2\pi m \frac{x_3 - \operatorname{ctg} \phi x_1}{|\mathbf{p}_2^0|}}_{\text{set 2}}. \quad (3.19) $$

Thus, the boundary condition in eq. (3.19) for equilibrium interface dislocations, combined with eq. (3.18), leads a set of 6 linear equations:

$$ \Sigma_1: \left\{ \begin{array}{l} \operatorname{Re} \sum_{\alpha=1}^{3} A \lambda_{A}^{\alpha} \mathbf{a}^{\alpha} - B \zeta_{B}^{\alpha} \mathbf{a}_{*}^{\alpha} = \mathbf{v} \\ \\ \operatorname{Im} \sum_{\alpha=1}^{3} A \lambda_{A}^{\alpha} \mathbf{a}^{\alpha} - B \zeta_{B}^{\alpha} \mathbf{a}_{*}^{\alpha} = \mathbf{0}, \end{array} \right. \quad (3.20) $$
where Im stands for the imaginary part of a complex quantity and $\vartheta$ is given by

$$ \vartheta = \begin{cases} -\frac{b_1}{n} & \text{if } m = 0 \quad (n \ge 1) \\ -\frac{b_2}{m} & \text{if } n = 0 \quad (m \ge 1) \\ 0 & \text{if } nm \neq 0 \quad (n, m \ge 1). \end{cases} \qquad (3.21) $$

### Boundary condition 4: No net tractions along the interfaces

The solution must satisfy the traction-free boundary condition along the interfaces:

$$ A \sigma (x_1, 0, x_3) n = B \sigma (x_1, 0, x_3) n, \qquad (3.22) $$

where $\sigma(x_1, 0, x_3)$ is reduced to the short-range stress field produced by the interface equilibrium dislocations when eqs. (3.1) are satisfied. In that case, the tractions at the interface read

$$ \sigma(x_1, 0, x_3) n = \operatorname{sgn}(x_2) \sum_{k \neq 0} e^{i2\pi k \cdot r} \sum_{\alpha=1}^{3} \lambda^{\alpha} h^{\alpha} + \zeta^{\alpha} h_{*}^{\alpha}, \qquad (3.23) $$

where the subsidiary complex vectors $h^\alpha$ are related to the vectors $a^\alpha$ by

$$ h^{\alpha} = (\mathbf{W}_{2}^{\mathrm{t}} + p^{\alpha} \mathbf{W}_{3}) a^{\alpha} = -p^{\alpha}^{-1} (\mathbf{W}_{1} + p^{\alpha} \mathbf{W}_{2}) a^{\alpha}, \qquad (3.24) $$

with $h_k^\alpha = H_{k2}^\alpha$. Boundary condition in eq. (3.22) together with eq. (3.23) leads the additional system of 6 linear equations:

$$ \Sigma_2 : \left\{ \begin{array}{l} \operatorname{Re} \sum_{\alpha=1}^{3} A \lambda^{\alpha} A h^{\alpha} - B \zeta^{\alpha} B h_{*}^{\alpha} = 0 \\ \\ \operatorname{Im} \sum_{\alpha=1}^{3} A \lambda^{\alpha} A h^{\alpha} - B \zeta^{\alpha} B h_{*}^{\alpha} = 0. \end{array} \right. \qquad (3.25) $$

The two latter conditions 3. and 4. may be rewritten in a eigenvalue problem for equilibrium interface dislocation arrays. Indeed, the elastic fields of these dislocations in an anisotropic bicrystal free of far-field strains are given in terms of the 12 eigenvalues Eval and 12 corresponding eigenvectors Evec with $\alpha = 1, 2, 3$, i.e.

$$ \begin{aligned} \operatorname{Eval} &= \{ \operatorname{Re}_A p^\alpha, \operatorname{Im}_A p^\alpha, \operatorname{Re}_B p^\alpha, \operatorname{Im}_B p^\alpha \} \\ \operatorname{Evec} &= \{ A a^\alpha, B a^\alpha, A h^\alpha, B h^\alpha \}. \end{aligned} \qquad (3.26) $$

All these quantities are determined by solving a 6-dimensional eigenvalue problem that may be recast with the aid of eqs. (3.24) into the form

$$ \mathbb{N} \begin{bmatrix} a^{\alpha} \\ h^{\alpha} \end{bmatrix} = p^{\alpha} \begin{bmatrix} a^{\alpha} \\ h^{\alpha} \end{bmatrix} \qquad (3.27) $$

where the real nonsymmetric 6 × 6 matrices N depend on the wavevectors and the stiffness constants for crystals A and B through the W matrices given by eqs. (3.10), i.e.

$$ \mathbb{N} = \begin{bmatrix} -\mathbf{W}_3^{-1} \mathbf{W}_2^{\mathrm{t}} & \mathbf{W}_3^{-1} \\ -\mathbf{W}_1 + \mathbf{W}_2 \mathbf{W}_3^{-1} \mathbf{W}_2^{\mathrm{t}} & -\mathbf{W}_2 \mathbf{W}_3^{-1} \end{bmatrix}. \qquad (3.28) $$
<figure><img src="image_20.png" /><figcaption>FIGURE 3.3: The disregistry $\Delta u$ due to interface Volterra dislocations is a staircase function. It may be decomposed into an affine part $\Delta u_{aff}$ generated by a uniform distortion (represented by a continuous distribution of fictitious infinitesimal dislocations) and a sawtooth function $\Delta u_{dis}$ associated with the equilibrium interface dislocations in the natural state.</figcaption></figure>

<table><caption>TABLE 3.1: Material properties for copper, niobium, iron, aluminium, and nickel. The values of lattice parameters a for all materials are those listed by Gray [104] and elastic components $c_{11}$, $c_{12}$, and $c_{44}$ by Hirth and Lothe [7].</caption><thead><tr><th colspan="2">Properties</th><th colspan="5">Materials</th></tr><tr><th>Symbol</th><th>Unit</th><th>Cu</th><th>Nb</th><th>Fe</th><th>Al</th><th>Ni</th></tr></thead><tbody><tr><td>a</td><td>Å</td><td>3.615</td><td>3.301</td><td>2.866</td><td>4.050</td><td>3.524</td></tr><tr><td>c<sub>11</sub></td><td>GPa</td><td>168.4</td><td>246.0</td><td>242.0</td><td>108.2</td><td>246.5</td></tr><tr><td>c<sub>12</sub></td><td>GPa</td><td>121.4</td><td>134.0</td><td>146.5</td><td>61.3</td><td>147.3</td></tr><tr><td>c<sub>44</sub></td><td>GPa</td><td>75.4</td><td>28.7</td><td>112.0</td><td>28.5</td><td>124.7</td></tr></tbody></table>

Finally, the linear systems $\Sigma_1$ and $\Sigma_2$ are solved numerically to determine the 12 real constants $Ecst$, i.e.

$$ Ecst = \{ \mathrm{Re}_{\mathrm{A}}\lambda^{\alpha}, \mathrm{Im}_{\mathrm{A}}\lambda^{\alpha}, \mathrm{Re}_{\mathrm{B}}\zeta^{\alpha}, \mathrm{Im}_{\mathrm{B}}\zeta^{\alpha} \}, \quad (3.29) $$

completing the solutions of the elastic fields.

### 3.2.6 Interface elastic strain energy

Using the divergence theorem, the elastic strain energy $\gamma_e$ of equilibrium interface dislocation arrays may be expressed as a surface integral over a unit cell A of the interface dislocation network, i.e.

$$ \gamma_e(r_0) = \frac{1}{2A} \iint_{A(r_0)} \sigma(x_1, 0, x_3) \mathbf{n} \cdot \Delta \mathbf{u}_{\text{dis}}(x_1, x_3) \, dS, \quad (3.30) $$

where $\sigma(x_1, 0, x_3) \mathbf{n}$ is the total traction vector produced at the interface of interest. Stress fields at dislocation cores diverge, so regions near the cores must be excluded from the integral in eq. (3.30). Following standard practice [7], the domain of integration is limited to parts of the interface unit cell that are not within a pre-determined cutoff distance $r_0$ of the dislocation cores.

## 3.3 Symmetric example applications

The model described in the foregoing sections is applied to simple example interfaces: symmetric tilt and twist GBs as well as a pure misfit heterophase interface. The materials properties used in these examples are listed in Table 3.1.
### 3.3.1 Pure tilt grain boundary

Tilt boundaries that contain one set of interfacial dislocations have been discussed extensively [237]. To illustrate and validate the present method, a symmetrical tilt boundary with [001] tilt axis and tilt angle $\theta = 2^\circ$ is analyzed in detail. The calculations are carried out for Cu, which has a moderately high anisotropy ratio, $A_{Cu} = 2c_{44}/(c_{11} - c_{12}) = 3.21$. The boundary consists of one set of straight parallel dislocations with Burgers vector content $B$, expressed as

$$B = \left( \frac{n \times \xi}{d} \cdot p \right) b = \underbrace{(\mathbf{R}_+^{-1} - \mathbf{R}_-^{-1})}_{\mathbf{T}} p = 2 \sin \theta/2 \ p \times \omega. \qquad (3.31)$$

Here, the "median lattice" is used as the obvious reference state: the mapping matrices $F$ have been replaced by rotation matrices $R$, with $R_+$ representing a rotation of the upper crystal by angle $\theta_+ = \theta/2$ about the tilt axis and $R_-$ the rotation $\theta_- = -\theta/2$ of the adjacent lower crystal. Equation (3.31) is known as Frank's formula [91, 48], which gives the density of interface dislocations needed to create the tilt boundary. Selecting $b = a_{Cu}$ [010] $\parallel n$, eq. (3.31) shows that $\xi = [001]$ and $d = 10.3567$ nm.

As expected, the far-field stresses vanish for this choice of reference lattice, and only non-zero stresses are short-ranged. Figure (3.4) plots interface stresses as a function of $x_1$ and $x_2$ (the stresses are invariant along the dislocation line direction, $x_3$). The red contour illustrates where the stresses fall to zero when $|x_2| \ge 7 - 10$ nm (depending on the stress components), showing that their range is comparable to the dislocation spacing. The far-field rotations may be calculated from the antisymmetric part of the far-field distortions, i.e. $\Omega^\infty = \mathbf{D}_{\text{dis}}^\infty$. They satisfy $\Omega_+^\infty - \Omega_-^\infty = \mathbf{T}$ and yield a net non-vanishing rotation about the tilt axis, as excepted [180, 121]:

$$\omega = \omega_+^\infty - \omega_-^\infty = - \begin{pmatrix} 0 \\ 0 \\ 0.03490 \end{pmatrix} = - \frac{x_1 \times b}{d}. \qquad (3.32)$$

The disregistry $\Delta u_2$ and the displacement discontinuity $\Delta u_{2\text{dis}}$ associated with the Volterra and equilibrium tilt boundary dislocations are plotted in Fig. (3.5a). They are in good quantitative agreement with the applied boundary conditions, represented by staircase and sawtooth curves.

The average elastic energy per unit interface area $\gamma_e$ is determined for several values of the core cutoff parameter $r_0$. Following eq. (3.30), $\gamma_e$ may be written as

$$\gamma_e (r_0) = \frac{1}{2d} \int_{r_0}^{d-r_0} \underbrace{\sigma_{22} (x_1, 0, 0) \Delta u_{2\text{dis}} (x_1, 0)}_{W} dx_1. \qquad (3.33)$$

The variation of stress component $\sigma_{22}$ at $x_2 = 0$ with $x_1$ is plotted as a black line in Fig. (3.5b). The core region is shaded in grey. Local contributions to the interface elastic energy $W$ (values of the integrand in eq. (3.33)) are plotted in red. The average elastic energy per unit interface area will depend on the choice of $r_0$. For example, $\gamma_e = 142.8 \text{ mJ.m}^{-2}$ with $r_0 = b/2$ and $\gamma_e = 167.8 \text{ mJ.m}^{-2}$ with $r_0 = b/3$, where $b$ is the magnitude of $b$. An appropriate $r_0$ value is selected by comparing the interface elastic energies computed with the present dislocation-based method to experimentally measured energies of small angle [001] tilt boundaries [102], plotted as solid triangles in Fig. (3.6). The calculations using $r_0 = b/2$ are in good agreement with the experiments up to $\sim 5^\circ$ while $r_0 = b/3$ fits better in the range of $\sim 5 - 12^\circ$. The classical energy per unit area given by Read and Shockley [213], $E_{\text{RS}}(\theta) = 1450 \theta (-3 - \ln \theta) \text{ mJ.m}^{-2}$, is also shown in Fig. (3.6). It compares well with the calculations for $r_0 = b/3$.
<figure><img src="image_21.png" /><figcaption>FIGURE 3.4: Contour plots of stress components (a) $\sigma_{11}$ and (b) $\sigma_{22}$, for the 2° symmetric tilt boundary described in the text. The negative values (compression) are plotted in light grey, and the positive values (extension) in dark grey. The stresses decay away over distances comparable to the interface dislocation spacing. In red, the stress field values are equal to zero.</figcaption></figure>

### 3.3.2 Twist grain boundary

As shown in Fig. (3.7a), small-angle (010) twist GBs contain two sets of dislocations, so their dislocation content $B$ is expressed as

$$B = \left( \frac{n \times \xi_1}{d_1} \cdot p \right) b_1 + \left( \frac{n \times \xi_2}{d_2} \cdot p \right) b_2 = (\mathbf{R}_+^{-1} - \mathbf{R}_-^{-1}) p. \quad (3.34)$$

The twist boundaries of angle $\theta = 2°$ is considered in Cu, where the rotation axis is perpendicular to the boundary, $\omega = x_2 = [010]$. As in the case of the tilt boundary, the obvious reference state for twist boundaries is the "median lattice" suggested by Frank [92]. In this state, the total rotation across the boundary is equally partitioned between the two grains. However, to illustrate the importance of selecting the correct reference state, other possible reference states are considered. A common choice is to use of the adjacent crystal grains as the reference state. There is a continuum of other possible reference states between these two extremes, and the angle $\theta_c = -\kappa\theta$ is introduced to define the rotation of the reference state from the case where the upper crystal above the boundary has been chosen as the reference lattice. Here, $\kappa$ is a dimensionless parameter that varies from 0 to 1. Equipartitioning of rotations between the adjacent crystals (i.e. the "median lattice") occurs when $\kappa = 1/2$.

Section 3.2.3 demonstrated that interface dislocation geometry is independent of reference state. In this example, the twist boundary contains an orthogonal grid of dislocations with line directions $\xi_1 = 1/\sqrt{2}[\bar{1}01]$ and $\xi_2 = 1/\sqrt{2}[101]$. The spacings between successive parallel dislocations are $d_1 = d_2 = d = 7.3233$ nm. Because of the pure twist misorientations, the coherency stress fields are zero for all possible reference states. Figure (3.7b) plots the dependence of non-vanishing far-field stress components on $\kappa$. If a reference state with $\kappa = 0$ is chosen, then the interface dislocations deviate by $1°$ from pure screw character and possess non-zero far-field stress components $\sigma_{11+}^{\infty} = \sigma_{33+}^{\infty}$ and $\sigma_{11-}^{\infty} = \sigma_{33-}^{\infty}$. This demonstrates that $\kappa = 0$ does not represent the correct reference state since eqs. (3.1) (and eqs. (3.13) via eq. (3.12b)) are not satisfied. Furthermore, the far-field rotation with $\kappa = 0$ does not equal $2°$, but discrepancies on the order of $0.001°$ between the rotation vector component and the prescribed misorientation are found. As $\kappa$ increases, the
<figure><img src="image_22.png" /><figcaption>FIGURE 3.5: (a) Disregistries $\Delta u_2$ (staircase function) and $\Delta u_{2dis}$ (sawtooth function) computed using 100 harmonics for the 2° symmetric tilt boundary described in the text. (b) Stress distribution $\sigma_{22}$ and local elastic energy density $\gamma_e$ at the GB.</figcaption></figure>

far-field stresses decrease and eventually reach zero at $\kappa = 1/2$, as expected. The interface dislocations have perfect screw characters for this reference state, where non-zero far-field stresses are again obtained when $\kappa$ is increased beyond $\kappa = 1/2$.

Taking $\kappa = 1/2$, the elastic strain energy per unit area $\gamma_e$ is calculated for the twist GB using the expression:

$$ \gamma_e (r_0) = \frac{1}{2A} \iint_{r_0}^{d-r_0} (W_{(1)} + W_{(2)} + W_{(1-2)}) dx_1 dx_3, \quad (3.35) $$

with $A = |\mathbf{p}_1^0 \times \mathbf{p}_2^0|$ the area of the interface unit cell. Equation (3.35) is decomposed into self-energy densities $W_{(1)}$ and $W_{(2)}$ for each set of parallel dislocations and the interaction energy density $W_{(1-2)}$ between the two sets. These energies are obtained from the separate elasticity solutions for each set of dislocations:

$$ \begin{aligned} W_{(1)} + W_{(2)} &= \sigma_{23(1)}(x_1, 0, 0) \Delta u_{3dis(1)}(x_1, 0) + \sigma_{12(2)}(0, 0, x_3) \Delta u_{1dis(2)}(0, x_3) \\ W_{(1-2)} &= \sigma_{23(1)}(x_1, 0, 0) \Delta u_{1dis(2)}(0, x_3) + \sigma_{12(2)}(0, 0, x_3) \Delta u_{3dis(1)}(x_1, 0). \end{aligned} \quad (3.36) $$

The local self- and interaction energies are shown in Figs. (3.8a) and (b), respectively. The integral of the interaction energy $W_{(1-2)}$ over area $A$ is zero for any value $r_0$, in agreement with the classical dislocation theory result that orthogonal screw dislocations do not exert any forces on each other [7]. The total elastic energy is plotted in Fig. (3.9) as a function of the twist angle up to 12° for three core cutoff parameters: $r_0 = b_1/2$, $r_0 = b_1/3$, and $r_0 = b_1/4$.

### 3.3.3 Pure misfit interface

Lastly, the model is illustrated on an Al/Ni heterophase interface. The terminal planes of both adjacent crystals are fcc (010) planes. The [100] and [001] directions of both crystals are parallel in the interface plane. Thus, the interface is in the cube-on-cube orientation and contains two sets of
<figure><img src="image_23.png" /><figcaption>FIGURE 3.6: Interface elastic energies $\gamma_e$ computed using two different core cutoff parameters $r_0$ for a [001] tilt GB in Cu as a function of the tilt angle $\theta$. The gray line shows the Read-Shockley solution. Experimental values are shown with solid triangles [102].</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 3.7: (a) Small-angle twist GB on a (010) plane containing two sets of orthogonal dislocations. (b) Dependence of far-field stresses on $\kappa$ for the 2° twist boundary described in the text.</figcaption></figure>

parallel dislocations. Following eq. (3.2), the Burgers vector content $\mathbf{B}$ is written as

$$ \mathbf{B} = \left( \frac{\mathbf{n} \times \boldsymbol{\zeta}_1}{d_1} \cdot \mathbf{p} \right) \mathbf{b}_1 + \left( \frac{\mathbf{n} \times \boldsymbol{\zeta}_2}{d_2} \cdot \mathbf{p} \right) \mathbf{b}_2 = \underbrace{(\mathbf{A}_{\text{Al}}\mathbf{S}^{-1}(r_{\text{Al}}) - \mathbf{N}_{\text{Ni}}\mathbf{S}^{-1}(r_{\text{Ni}}))}_{\mathbf{T}} \mathbf{p}. \quad (3.37) $$

The reference state for this interface is a crystal oriented identically to the Al and Ni in their natural state, but strained such that its lattice constant in the interface plane is $a_c$, with $a_{Ni} \le a_c \le a_{Al}$. Only strains within the interface are necessary to ensure coherency: normal strains are not required. Thus, the matrix $\mathbf{T}$ in eq. (3.37) is composed of two equibiaxial stretch matrices (no rotations), $\mathbf{A}_{\text{Al}}\mathbf{S}^{-1} = \mathbf{A}_{\text{Al}}\mathbf{E}_c + \mathbf{I}$ and $\mathbf{N}_{\text{Ni}}\mathbf{S}^{-1} = \mathbf{N}_{\text{Ni}}\mathbf{E}_c + \mathbf{I}$, where $\mathbf{I}$ represents the identity matrix. These mapping matrices depend on the ratios of lattice parameters between Al and Ni in their natural and reference states, $r_{\text{Al}} = a_{\text{Al}}/a_c \ge 1$ and $r_{\text{Ni}} = a_{\text{Ni}}/a_c \le 1$. The matrix $\mathbf{T}$ in eq. (3.37) may also be rewritten as the difference between the coherency strains prescribed in Al and Ni:

$$ \mathbf{A}_{\text{Al}}\mathbf{E}_c - \mathbf{N}_{\text{Ni}}\mathbf{E}_c = \mathbf{T}. \quad (3.38) $$

Following the procedure described in section 3.2.4, Ni is initially chosen as the reference lattice, so that $r_{\text{Al}} = a_{\text{Al}}/a_{\text{Ni}}$ and $r_{\text{Ni}} = 1$, and identify $\tilde{\mathbf{b}}_1 = a_{\text{Ni}}/\sqrt{2} [101]$ and $\tilde{\mathbf{b}}_2 = a_{\text{Ni}}/\sqrt{2} [10\bar{1}]$. Then, using eq. (3.3), an interface that consists of an orthogonal grid of edge dislocations with $\boldsymbol{\zeta}_1 = 1/\sqrt{2} [101]$
<figure><img src="image_25.png" /><figcaption>FIGURE 3.8: Local (a) self- {$W_{(1)} + W_{(2)}$} and (b) interaction $W_{(1-2)}$ elastic energies arising from two sets of orthogonal screw dislocations in a 2° twist boundary on a (010) plane in Cu.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIGURE 3.9: Elastic energies per unit area $\gamma_e$ as a function of the rotation angle $\theta$ of twist GBs along (010) planes in Cu for three core cutoff parameters $r_0$.</figcaption></figure>

and $\xi_2 = 1/\sqrt{2}[101]$ is found, and the corresponding dislocation spacings $d_1 = d_2 = 1.902$ nm. Using this choice of reference state, the far-field strains produced by the interface dislocations are:

$$ \text{Al}\mathbf{E}_{\text{dis}}^{\infty} = \begin{bmatrix} 0.10133 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0.10133 \end{bmatrix}, \quad \text{and,} \quad \text{Ni}\mathbf{E}_{\text{dis}}^{\infty} = \begin{bmatrix} -0.03243 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -0.03243 \end{bmatrix}, \quad (3.39) $$

such that the matrices in eqs. (3.39) satisfy

$$ -(\text{Al}\mathbf{E}_{\text{dis}}^{\infty} - \text{Ni}\mathbf{E}_{\text{dis}}^{\infty}) = \mathbf{T}. \quad (3.40) $$

Combining eqs. (3.38) and (3.40), it follows

$$ \text{Al}\mathbf{E}_c + \text{Al}\mathbf{E}_{\text{dis}}^{\infty} = \underbrace{\text{Ni}\mathbf{E}_c}_{=0} + \text{Ni}\mathbf{E}_{\text{dis}}^{\infty} = \begin{bmatrix} -0.03243 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -0.03243 \end{bmatrix} \neq \mathbf{0} \quad (\Leftrightarrow \text{Al}\mathbf{E}^{\infty} = \text{Ni}\mathbf{E}^{\infty}), \quad (3.41) $$
<figure><img src="image_27.png" /><figcaption>FIGURE 3.10: Dependence of the total far-field strain component $a_{Al}E_{11}^{\infty}$ in Al on $\delta$ for a Al/Ni heterophase interface. The red dotted line gives the unique reference state, for which the far-field decay to zero and the coherent parameter $a_c$ is defined. The lattice parameter $\bar{a} = 2a_{Al}a_{Ni} / (a_{Al} + a_{Ni})$, which is a good approximation for an interface between crystals of different lattice parameters but identical elastic constants [93, 131], is marked by a grey cross symbol.</figcaption></figure>

with $a_{Ni}E_c = 0$ here, because Ni has been chosen as the reference lattice. However, according to eq. (3.41b), condition 2. given by eq. (3.13) is not satisfied since the total far-field strains in each individual material do not decay to zero when $x_2 \to \pm\infty$. This demonstrates that the initial choice of reference state is not correct.

To find the correct reference state, a variable $\delta$, with $0 \le \delta \le 1$ that interpolates $a_c$ between $a_{Al}$ and $a_{Ni}$ is introduced as follows

$$a_c = \delta a_{Al} + (1 - \delta) a_{Ni}. \qquad (3.42)$$

It is shown that the far-field strains in Al and Ni are equal for all $\delta$, so that eq. (3.41a) is always satisfied, i.e. $a_{Al}E_c^{\infty} = a_{Ni}E_c^{\infty}$ with $a_{Ni}E_c = 0$ if $\delta = 0$, and $a_{Al}E_c = 0$ if $\delta = 1$. However, only one unique reference state (corresponding to an unique value of $\delta$) gives vanishing far-field strains in the bicrystal in its natural state by satisfying eq. (3.13) as well. The pure misfit interface example serves to show that eq. (3.41a) is a necessary, but not sufficient condition for determining the reference state.

The total far-field strain component $a_{Al}E_{11}^{\infty}$ in Al is plotted in Fig. (3.10) as a function of $\delta$ and is identical to the component $a_{Ni}E_{33}^{\infty}$, according to the interface symmetry (all other strain components are zero). Because eq. (3.41a) is verified for all $\delta$, the same components in Ni give the same plot as in Fig. (3.10). The far-field strains vary linearly with $\delta$ and become zero when $\delta = 0.21787$, so that $a_c = 0.36386$ nm. This value of $a_c$ is the unique coherent reference state for which the pure misfit Al/Ni interface of interest is consistent with the Frank-Bilby equation. It is closer to $a_{Ni}$ than to $a_{Al}$ because Ni is the stiffer of these two materials and so carries a lower coherency strain in the reference state. The far-field rotations are zero for all values of $\delta$, as excepted.

To demonstrate the errors that come about from ignoring the unequal partitioning of elastic fields and to validate the current calculation, $a_c$ is recomputed under the assumption that both sides of the interface have the same stiffness (equal to that of Al or Ni), but different natural lattice parameters ($a_{Al}$ and $a_{Ni}$, as the original calculation). For this case, the calculated value for $a_c$ is in very good agreement with the well-known approximate result $\bar{a} = 2a_{Al}a_{Ni} / (a_{Al} + a_{Ni}) = 0.37687$ nm [93, 131], corresponding to $\delta = 0.46521$. This value, however, is far from the correct lattice parameter of the reference state when the differing stiffnesses of Al and Ni are taken into account, as illustrated by cross symbols in Fig. (3.10). It is also shown that $\bar{a}$ deviates from the prediction and is not consistent with the Frank-Bilby equation when the heterogeneous distortions of bicrystals are explicitly described at equilibrium.
<figure><img src="image_28.png" /><figcaption>FIGURE 3.11: (a) The reference and natural states of an interface are related by transformation matrices $_{fcc}$F and $_{bcc}$F. (b) The correspondence between a closed right-handed circuit enclosing the probe vector $\mathbf{p}$ in the natural state and its corresponding path with closure failure B in the reference state.</figcaption></figure>

3.4 Partitioning of elastic distortions at fcc/bcc interfaces

In this section, the study is focused on semicoherent heterophase interfaces comprised of two
sets of dislocations and formed along closest-packed planes in fcc/bcc bimetals, especially for
fcc{111}/bcc{110} (Cu/Nb, Ag/V, and Cu/Mo) interfaces in the Nishiyama-Wassermann (NW)
orientation relations (OR) [279, 192] as well as in ORs that differ from the NW by an in-plane
twist rotation. It is showed that elastic distortions, i.e. strains as well as tilt and twist rotations,
are in general unequally partitioned at such interfaces. The correct partitioning of these fields
determines the coherent reference state for which the bicrystal of interest is free of far-field strains.
Using these results, the stress fields generated by misfit dislocation patterns are computed and
analyzed for the Cu/Nb system in the NW and Kurdjumov-Sachs (KS) [154] ORs. The dislocation
structure (i.e. the Burgers vectors, spacings, and line directions) is also determined in lowest strain
energy solutions of the Frank-Bilby equation along a specific transformation pathway between the
NW and KS ORs.

Similarly to Fig. (3.1), the concept of reference and natural states of an interface is depicted
in Fig. (3.11). The natural state contains an interface formed by joining two crystals with pre-
scribed misorientation and interface planes as well as vanishing far-field strains. This state is also
related to a single crystal, coherent reference state by uniform displacement gradients $_{A}$**F = $_{fcc}$**F
and $_{B}$**F = $_{bcc}$**F, which map the reference state to the natural state, as shown in Fig. (3.11a). In the
reference state, the two adjacent materials that meet at the interface are rotated and strained such
that they are in perfect registry with each other across the **X** − **Z** interface plane after bonding. In
general, these displacement gradients entail interface misorientations that have both tilt and twist
components [7, 237, 123]. Again, the interface along the **x** − **z** plane is not coherent in the natural
state, but rather semicoherent due to the presence of misfit dislocations.

The atomically sharp fcc{111}/bcc{110} interfaces in NW and in-plane twisted-NW ORs con-
tain two periodic arrays of infinitely long, straight, and uniformly spaced dislocations. In the NW
OR, one of the <110> directions in a fcc {111} plane lies parallel to the <100> direction in a bcc
{110} plane [279, 192]. The in-plane twisted-NW ORs considered here differ from the NW OR
only by a twist rotation of one crystal (here, the bcc material) with respect to the adjacent (fcc)
crystal about the axis normal to the interface. The procedure described in section 3.2.4 is adopted
<table><caption>TABLE 3.2: Material properties for copper (Cu), niobium (Nb), silver (Ag), vanadium (V), and molybdenum (Mo). The values of stiffness constants ĉ₁₁, ĉ₁₂, ĉ₄₄, and lattice parameters a for all materials are those listed in Ref. [258].</caption><thead><tr><th>Systems</th><th>ĉ<sub>11</sub> (GPa)</th><th>ĉ<sub>12</sub> (GPa)</th><th>ĉ<sub>44</sub> (GPa)</th><th>a (Å)</th></tr></thead><tbody><tr><td>Cu</td><td>178.8</td><td>122.6</td><td>81.03</td><td>3.615</td></tr><tr><td>Nb</td><td>245.6</td><td>133.7</td><td>28.8</td><td>3.3008</td></tr><tr><td>Ag</td><td>124.2</td><td>93.9</td><td>46.1</td><td>4.090</td></tr><tr><td>V</td><td>220.15</td><td>130.7</td><td>42.8</td><td>3.039</td></tr><tr><td>Cu</td><td>187.8</td><td>125.7</td><td>70.6</td><td>3.615</td></tr><tr><td>Mo</td><td>545.9</td><td>219.3</td><td>108.8</td><td>3.147</td></tr></tbody></table>

to determine the unique reference states that meet the condition of vanishing far-field strains and prescribed misorientation for such interfaces. Thus, the dislocation content B of an interface, intersected by a probe vector p contained within the interface plane as illustrated in Fig. (3.11b), is described by the Frank-Bilby equation in eq. (3.2). For interfaces in the NW OR, a transformation pathway is defined by continuously adjusting the reference state from the strain-free state of the fcc crystal present at the interface to that of the adjacent bcc crystal. For all reference states along this path, the method described in section 3.2 is used to compute the superposition of the uniform coherency strains, E_c, needed to maintain perfect registry and the far-field strain fields produced by the Volterra dislocation arrays, E_dis. In the correct reference state, these quantities cancel and the total far-field strain field E vanishes in both upper fcc (ŷ > 0) and lower bcc (ŷ < 0) materials, as defined by eqs. (3.13), as

$$ \lim_{\hat{y} \to \pm\infty} \mathbf{E}(\hat{x}, \hat{y}, \hat{z}) = \mathbf{0} \Leftrightarrow \begin{cases} \mathbf{f}_{\text{fcc}} \mathbf{E}^{\infty} = \mathbf{f}_{\text{fcc}} \mathbf{E}_{\text{c}} + \mathbf{f}_{\text{fcc}} \mathbf{E}_{\text{dis}}^{\infty} = \mathbf{0} \\ \mathbf{b}_{\text{bcc}} \mathbf{E}^{\infty} = \mathbf{b}_{\text{bcc}} \mathbf{E}_{\text{c}} + \mathbf{b}_{\text{bcc}} \mathbf{E}_{\text{dis}}^{\infty} = \mathbf{0}, \end{cases} \quad (3.43) $$

for which the far-field rotation state in the NW OR is consistent with the given crystallographic character (interface plane and misorientation).

To find the reference state for interfaces differing from those in the NW ORs by an in-plane twist angle θ, a second pathway is defined by rotating the previously determined reference state in the NW OR from 0 to θ. Along this second path, the rotated reference state, for which eqs. (3.43) are satisfied, also yields far-field rotations that must be consistent with the in-plane prescribed twist misorientations. Using the correct reference states for all ORs, the short-range interface strains and stresses that arise from the incomplete cancellation of the coherency and Volterra dislocation fields near the interfaces are also computed as well as the interface elastic energy γ_e from eq. (3.30) as a surface integral over a unit cell. The domain of integration is related to a pre-determined cutoff distance r₀ of the dislocation cores to determine the likeliest interface misfit dislocation configurations whenever the Frank-Bilby equation (eq. (3.2)) admits multiple solutions.

In this section, a detailed discussion of partitioning of distortions at Cu/Nb interfaces is presented, while analogous results for Ag/V and Cu/Mo interfaces are shown, albeit without going into detail. The material properties (elastic constants and lattice parameters) used in all calculations for these three interface types are listed in Table 3.2.

### 3.4.1 Mapping between states in the Nishiyama-Wassermann orientations

Without loss of generality, the following specific relation is used among 12 possible equivalent variants of the NW OR [106] to construct the mapping from the fcc to the bcc crystal:

$$ \text{NW}: \left\{ \begin{array}{l} \hat{x} \parallel \mathbf{x}_{\text{fcc}} = [11\bar{2}]_{\text{fcc}} \parallel \mathbf{x}_{\text{bcc}} = [01\bar{1}]_{\text{bcc}} \\ \mathbf{n} \parallel \hat{y} \parallel \mathbf{y}_{\text{fcc}} = [111]_{\text{fcc}} \parallel \mathbf{y}_{\text{bcc}} = [011]_{\text{bcc}} \\ \hat{z} \parallel \mathbf{z}_{\text{fcc}} = [1\bar{1}0]_{\text{fcc}} \parallel \mathbf{z}_{\text{bcc}} = [100]_{\text{bcc}}. \end{array} \right. \quad (3.44) $$