equation

(1) $\BODY$

Habilitation à Diriger des Recherches
Interfaces in crystalline materials
Aurélien Vattré
Maître de Recherche à l’ONERA
<EMAIL_ADDRESS>

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

Ioan Ionescu |  | Président
---|---|---
Brigitte Bacroix |  | Rapporteure
Stéphane Berbenni |  | Rapporteur
Marc Fivel |  | Rapporteur
Renald Brenner |  | Examinateur

## Abstract

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

Far from being featureless dividing surfaces between neighboring crystals,
elucidating features 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
polymorphic 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 combining nonlinear
elastoplasticity and multivariant phase-field approach at large strains. The
calculations investigate the crucial role played by the plastic deformation in
the morphological and microstructure evolution processes under high
hydrostatic compression and shock-wave conditions. 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 sextic
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 defects and
semicoherent interfaces.

###### Contents

  1. 1 Introduction
  2. 2 Crystalline interfaces during solid-solid phase transitions in iron
    1. 2.1 Motivation
    2. 2.2 A phase-field model coupled with finite elastoplasticity
      1. 2.2.1 Kinematics
      2. 2.2.2 Balance laws
      3. 2.2.3 The Clausius-Duhem inequality
      4. 2.2.4 Constitutive equations
      5. 2.2.5 Multiple reaction pathways and energy landscape
      6. 2.2.6 Computational framework
    3. 2.3 Pure hydrostatic compression
      1. 2.3.1 Material and model inputs
      2. 2.3.2 Analysis of the pressure-volume responses
      3. 2.3.3 Microstructure and variant selection
    4. 2.4 Shock wave propagation
      1. 2.4.1 The internal structure of shock waves
      2. 2.4.2 Effect of plasticity in shock-loaded iron
      3. 2.4.3 Residual stresses in the plastically-deformed microstructure
      4. 2.4.4 Dynamical instability in structural phase transitions
    5. 2.5 Limitations
  3. 3 Dislocation structures and energetics at heterophase interfaces
    1. 3.1 Motivation
    2. 3.2 Determining the Burgers vectors of interface dislocation arrays
      1. 3.2.1 Planar interfaces in linear elastic bicrystals
      2. 3.2.2 Volterra dislocations in the reference state
      3. 3.2.3 Crystallographic constraints on interface dislocations
      4. 3.2.4 Solution strategy
      5. 3.2.5 Elastic fields of interface dislocation arrays
      6. 3.2.6 Interface elastic strain energy
    3. 3.3 Symmetric example applications
      1. 3.3.1 Pure tilt grain boundary
      2. 3.3.2 Twist grain boundary
      3. 3.3.3 Pure misfit interface
    4. 3.4 Partitioning of elastic distortions at fcc/bcc interfaces
      1. 3.4.1 Mapping between states in the Nishiyama-Wassermann orientations
      2. 3.4.2 Far-field strains and rotations
      3. 3.4.3 Spurious fields from incorrect reference states
      4. 3.4.4 Orientations differing from the Nishiyama-Wassermann relations
      5. 3.4.5 Short-range elastic fields
      6. 3.4.6 Comparison with atomistic simulations
    5. 3.5 Application to the sink strength of semicoherent interfaces
      1. 3.5.1 Computational multi-model strategy
      2. 3.5.2 Kinetic Monte Carlo simulations with elastic interactions
      3. 3.5.3 Effect of elastic interactions on interface sink strength
    6. 3.6 Elastic strain relaxation in interfacial dislocation patterns
      1. 3.6.1 General considerations on hexagonal-shaped dislocation patterns
      2. 3.6.2 Solution methodology for strain-relaxed rearrangements
      3. 3.6.3 Parametric energy-based framework
      4. 3.6.4 Boundary conditions with surface/interface constitutive relations
      5. 3.6.5 Application to Au/Cu heterosystems
      6. 3.6.6 Comparison with atomistic simulations
    7. 3.7 Interaction with extrinsic dislocations in bimaterials
      1. 3.7.1 Extrinsic dislocation arrays and loops
      2. 3.7.2 Internal forces on intrinsic and extrinsic dislocations
      3. 3.7.3 On the piled-up dislocations in the (111)Cu/(011)Nb bimaterial
      4. 3.7.4 Limitations
    8. 3.8 Extension to non-singular fields in multilayered magneto-electro-elastic plates
      1. 3.8.1 Boundary-value problem and singularity-free field solutions
      2. 3.8.2 A primary case: 2D bilayered composites
      3. 3.8.3 Energy-based criterion for interlayers in A/B/A trilayers
      4. 3.8.4 Dislocation-induced response under applied external loading
  4. 4 Conclusion and future works
    1. 4.1 Concluding remarks
    2. 4.2 Perspectives
      1. 4.2.1 Thermoelasticity of semicoherent interfaces
      2. 4.2.2 Distributed dislocations for periodic networks of cracks
      3. 4.2.3 Towards a general treatment for {interfaces, dislocations, cracks}

## 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}=\\{$ Achille (4 mois), Anton (2 ans),
Aurélie $\\}^{\,\star}$ 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 1: ${}^{\star}\,$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 »

##  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 broder 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 introductive 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 this 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 materials. 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 dislocation 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 interface 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
corresponding 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 bilayers 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 elasticity 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

  * [P1]

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.

  * [P2]

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.

  * [P3]

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.

  * [P4]

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.

### 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 [236]. 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 (diffusionless) 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
$\sim 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 $\alpha$-Fe) has shown a magnetic and structural transition
to the nonmagnetic hexagonal close-packed phase (hcp $\epsilon$-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 [221], 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 [144,
292]. 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 [48, 183, 23], a major problem
deals with the accuracy in determining the energy landscape for the bcc-to-hcp
transition [83, 172]. 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
[141, 142, 143]. 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 [133, 219] and large [152, 181] 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
[161, 9, 146, 74, 296]. Treatments of thermodynamics and kinetic relations in
phase-field approaches are related to the pioneering works by [49] and [7],
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 [259, 44]. 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 [75]. 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 [259] and shock-wave [266] 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.

### 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_{\scalebox{0.65}{0}}\subset\mathbb{R}^{3}$, for which
the motion of $\Omega_{\scalebox{0.65}{0}}$ is given by the mapping
$\textbf{{x}}=\boldsymbol{\chi}\left(\textbf{{X}},t\right)\colon\\!\Omega_{\scalebox{0.65}{0}}\rightarrow\Omega\subset\mathbb{R}^{3}$
with respect to time $t$. The total deformation gradient $\mathbb{F}$ is
related to the following multiplicative decomposition [162, 163, 152, 166],
i.e.,

(2.1)
$\mathbb{F}=\left.\frac{\partial\boldsymbol{\chi}}{\partial\textbf{{X}}}\right|_{\scriptsize
t}=\boldsymbol{\nabla}\boldsymbol{\chi}=\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}\\!\cdot\\!\mathbb{F}\mathrm{p}\,,$

with $\boldsymbol{\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
$\mathbb{F}\mathrm{e}$, plastic $\mathbb{F}\mathrm{p}$, and, transformational
$\mathbb{F}\mathrm{t}$ distortions, leading to the pressure-induced phase
transformation from the bcc to hcp phases.

Similarly to classical crystal elastoplasticity theories [155, 160], the
decomposition eq. (2.1) is not uniquely defined and different ordering
relations have been taken into account in the literature [245]. Because the
local irreversible plastic deformation $\mathbb{F}\mathrm{p}$ 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 $\mathbb{F}\mathrm{t}$ occurs between $\mathbb{F}\mathrm{p}$ and
$\mathbb{F}\mathrm{e}$, where the elastic contribution accounts for the
lattice stretching $\mathbb{U}\mathrm{e}$ and rotation $\mathbb{R}\mathrm{e}$.
The polar decomposition to $\mathbb{F}\mathrm{e}$ reads:
$\mathbb{F}\mathrm{e}=\mathbb{R}\mathrm{e}\cdot\mathbb{U}\mathrm{e}$, with
$\mathbb{U}\mathrm{e}^{2}=\mathbb{F}\mathrm{e}^{\mathrm{t}}\cdot\mathbb{F}\mathrm{e}$,
and,
$\mathop{\mathrm{det}}\mathbb{F}\mathrm{e}=\mathop{\mathrm{det}}\mathbb{U}\mathrm{e}=j_{\mathrm{e}}$.
The superscript t denotes the transpose operation. Although the controversy
regarding the decomposition is beyond the scope of this paper, both tensors
$\mathbb{F}\mathrm{p}$ and $\mathbb{F}\mathrm{t}$ describe here two
intermediate configurations, $\Omega\mathrm{p}$ and $\Omega\mathrm{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. [166]. It follows from eq.
(2.1) that the total spatial velocity gradient tensor $\mathbb{L}$ is given by

(2.2)
$\mathbb{L}=\dot{\mathbb{F}}\\!\cdot\\!\mathbb{F}^{-1}=\mathbb{L}\mathrm{e}+\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{L}\mathrm{t}\\!\cdot\\!\mathbb{F}\mathrm{e}^{-1}+\mathbb{F}\mathrm{et}\\!\cdot\\!\mathbb{L}\mathrm{p}\\!\cdot\\!\mathbb{F}\mathrm{et}^{-1}\,,$

with
$\mathbb{F}\mathrm{et}\\!=\\!\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}$.
The superposed dot in eq. (2.2) denotes the time derivative. The elastic
$\mathbb{L}\mathrm{e}$, transformational $\mathbb{L}\mathrm{t}$, and, plastic
$\mathbb{L}\mathrm{p}$ velocity strain tensors are similarly defined by

(2.3)
$\mathbb{L}\mathrm{e}=\dot{\mathbb{F}}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{e}^{-1},~{}\mathbb{L}\mathrm{t}=\dot{\mathbb{F}}\mathrm{t}\\!\cdot\\!\mathbb{F}\mathrm{t}^{-1},\mbox{~{}and}~{}~{}\mathbb{L}\mathrm{p}=\dot{\mathbb{F}}\mathrm{p}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-1}\,,$

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

  1. 1.

The measures of volume changes after each deformation processes satisfy:

(2.4)
$j_{\mathrm{e}}=\mathop{\mathrm{det}}\mathbb{F}\mathrm{e}>0,~{}j_{\mathrm{t}}=\mathop{\mathrm{det}}\mathbb{F}\mathrm{t}>0,\mbox{~{}and}~{}~{}\mathop{\mathrm{det}}\mathbb{F}\mathrm{p}=1\,,$

so that, $\mathbb{F}\mathrm{e}$, $\mathbb{F}\mathrm{t}$ and
$\mathbb{F}\mathrm{p}$ are invertible, and, the plastic flow preserves the
volume.

  2. 2.

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

(2.5)
$\dot{\mathbb{F}}\mathrm{p}=\mathbb{D}\mathrm{p}\\!\cdot\\!\mathbb{F}\mathrm{p},\mbox{~{}with}~{}\mathbb{D}\mathrm{p}=\mathrm{sym}\,\mathbb{L}\mathrm{p}=\mathbb{L}\mathrm{p}\,,$

where $\mathrm{sym}\,\mathbb{L}\mathrm{p}$ denotes the symmetric part of
$\mathbb{L}\mathrm{p}$.

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

Figure 2.1: Schematics of the reference $\Omega_{0}$, intermediate,
$\Omega\mathrm{p}$ and $\Omega\mathrm{t}$, and, current $\Omega$
configurations, for which the total deformation gradient tensor $\mathbb{F}$
is decomposed multiplicatively into plastic $\mathbb{F}\mathrm{p}$,
transformational $\mathbb{F}\mathrm{t}$ and elastic $\mathbb{F}\mathrm{e}$
distortions.

####  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

(2.6) $\boldsymbol{\nabla}\hskip 1.13791pt\\!\cdot\\!\hskip
1.28032pt\mathbb{P}+\rho_{\scalebox{0.65}{0}}\,\textbf{{b}}=\rho_{\scalebox{0.65}{0}}\,\ddot{\textbf{{u}}}\mbox{~{}~{}in}~{}\Omega_{\scalebox{0.65}{0}}\,,$

where $\mathbb{P}$ is the first (non-symmetric) Piola-Kirchhoff stress tensor,
b are external body forces per unit mass, and,
$\ddot{\textbf{{u}}}=\ddot{\textbf{{x}}}\left(\textbf{{X}},t\right)$ is the
acceleration of the material point X, with u the corresponding displacement
field, defined by
$\textbf{{u}}=\textbf{{x}}\left(\textbf{{X}},t\right)-\textbf{{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 X in
the domain of $\Omega_{\scalebox{0.65}{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 [64] as

(2.7)
$\int_{{}_{\\!\\!\Omega_{\mbox{0}}}}\left(\mathbb{P}\colon\\!\dot{\mathbb{F}}-\rho_{\scalebox{0.65}{0}}\,\dot{\psi}\right)\,\mathrm{d}\Omega_{\scalebox{0.65}{0}}\geq
0\,,$

where $\colon\\!$ denotes the double inner tensor product, and, $\psi$ the
specific Helmholtz free energy. Equation (2.7) shows that the first Piola-
Kirchhoff stress tensor $\mathbb{P}$ and the deformation gradient $\mathbb{F}$
are work-conjugate variables, while $\mathbb{P}\colon\\!\dot{\mathbb{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 postulated that the Helmholtz free energy can be written in the
following form:

(2.8)
$\psi\doteq\psi\left(\mathbb{F}\mathrm{e},\mathbb{F}\mathrm{t},\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\right)\,,$

where $\boldsymbol{\nabla}\mathbb{F}\mathrm{t}$ 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 $\psi_{\mathrm{e}}$, transformational $\psi_{\mathrm{t}}$, and,
purely empirical gradient penalty $\psi_{{}_{\nabla}}$ contributions. With the
aforementioned considerations, the Helmholtz free energy can thus be written
as

(2.9)
$\psi\doteq\psi_{\mathrm{e}}\left(\mathbb{F}\mathrm{e},\mathbb{F}\mathrm{t}\right)+\psi_{\mathrm{t}}\left(\mathbb{F}\mathrm{t}\right)+\psi_{{}_{\nabla}}\left(\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\right)\,,$

which, in contrast with ab-initio electronic structure calculations, is not
uniquely defined. However, such elastic/inelastic splitting, comparable to the
classical phase-field models with elastic and chemical potentials [275, 9], 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.,

(2.10) $\displaystyle\dot{\mathbb{F}}$
$\displaystyle=\dot{\mathbb{F}}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}\\!\cdot\\!\mathbb{F}\mathrm{p}+\mathbb{F}\mathrm{e}\\!\cdot\\!\dot{\mathbb{F}}\mathrm{t}\\!\cdot\\!\mathbb{F}\mathrm{p}+\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}\\!\cdot\\!\dot{\mathbb{F}}\mathrm{p}$
$\displaystyle\dot{\psi}$
$\displaystyle=\left.\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{e}}\right|_{\scriptsize\mathbb{F}\mathrm{t}}\\!\\!\colon\\!\dot{\mathbb{F}}\mathrm{e}+\left.\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{t}}\right|_{\scriptsize\mathbb{F}\mathrm{e}}\\!\\!\colon\\!\dot{\mathbb{F}}\mathrm{t}+\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}\colon\\!\dot{\mathbb{F}}\mathrm{t}+\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\therefore\boldsymbol{\nabla}\dot{\mathbb{F}}\mathrm{t}\,,$

where $\therefore$ 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

(2.11) $\displaystyle\int_{{}_{\\!\\!\Omega_{\mbox{0}}}}\bigg{\\{}$
$\displaystyle\bigg{(}\mathbb{P}\\!\cdot\\!\mathbb{F}\mathrm{p}^{\mathrm{t}}\\!\cdot\\!\mathbb{F}\mathrm{t}^{\mathrm{t}}-\rho_{\scalebox{0.65}{0}}\left.\\!\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{e}}\right|_{\scriptsize\mathbb{F}\mathrm{t}}\,\bigg{)}\\!\colon\\!\dot{\mathbb{F}}\mathrm{e}+\bigg{(}\mathbb{F}\mathrm{e}^{\mathrm{t}}\\!\cdot\\!\mathbb{P}\\!\cdot\\!\mathbb{F}\mathrm{p}^{\mathrm{t}}-\rho_{\scalebox{0.65}{0}}\left.\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{t}}\right|_{\scriptsize\mathbb{F}\mathrm{e}}-\rho_{\scalebox{0.65}{0}}\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}\,\bigg{)}\\!\colon\\!\dot{\mathbb{F}}\mathrm{t}+\boldsymbol{\Sigma}_{\ast}\colon\\!\mathbb{D}\mathrm{p}$
$\displaystyle-$
$\displaystyle\rho_{\scalebox{0.65}{0}}\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\therefore\boldsymbol{\nabla}\dot{\mathbb{F}}\mathrm{t}\bigg{\\}}\,\mathrm{d}\Omega_{\scalebox{0.65}{0}}\geq
0\,,$

where $\boldsymbol{\Sigma}_{\ast}$ is a work-conjugate stress measure related
to the first Piola-Kirchhoff $\mathbb{P}$, as follows

(2.12)
$\displaystyle\boldsymbol{\Sigma}_{\ast}=\mathbb{F}\mathrm{et}^{\mathrm{t}}\\!\cdot\\!\mathbb{P}\\!\cdot\\!\mathbb{F}\mathrm{p}^{\mathrm{t}}\,.$

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.,

(2.13)
$\displaystyle\int_{{}_{\\!\\!\Omega_{\mbox{0}}}}\bigg{(}\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\therefore\boldsymbol{\nabla}\dot{\mathbb{F}}\mathrm{t}\bigg{)}\,\mathrm{d}\Omega_{\scalebox{0.65}{0}}=-$
$\displaystyle\int_{{}_{\\!\\!\Omega_{\mbox{0}}}}\bigg{(}\boldsymbol{\nabla}\\!\cdot\\!\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\colon\\!\dot{\mathbb{F}}\mathrm{t}\bigg{)}\,\mathrm{d}\Omega_{\scalebox{0.65}{0}}+\int_{{}_{\\!\\!\Sigma_{\mbox{0}}}}\underbrace{\bigg{(}\dot{\mathbb{F}}\mathrm{t}\colon\\!\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\\!\cdot\textbf{{n}}\bigg{)}}_{\mbox{\footnotesize
surface dissipation}}\,\mathrm{d}\Sigma_{\scalebox{0.65}{0}}\,,$

where $\Sigma_{\scalebox{0.65}{0}}$ is a boundary of
$\Omega_{\scalebox{0.65}{0}}$ with unit outward normal n. Assuming that the
surface dissipation is absent during the transformational process, additional
boundary conditions as set of nine equations for phase transitions may also be
derived by

(2.14)
$\displaystyle\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\\!\cdot\textbf{{n}}=\mathbb{0},\mbox{~{}with}~{}~{}\dot{\mathbb{F}}\mathrm{t}\neq\mathbb{0}\mbox{~{}at}~{}\Sigma_{\scalebox{0.65}{0}}\,,$

corresponding to the orthogonality relations between
$\boldsymbol{\nabla}\mathbb{F}\mathrm{t}$ and the external surfaces
$\Sigma_{\scalebox{0.65}{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 $\pazocal{D}$, as follows

(2.15)
$\displaystyle\pazocal{D}=\bigg{(}\mathbb{P}\\!\cdot\\!\mathbb{F}\mathrm{p}^{\mathrm{t}}\\!\cdot\\!\mathbb{F}\mathrm{t}^{\mathrm{t}}-\rho_{\scalebox{0.65}{0}}\left.\\!\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{e}}\right|_{\scriptsize\mathbb{F}\mathrm{t}}\,\bigg{)}\\!\colon\\!\dot{\mathbb{F}}\mathrm{e}+\mathbb{X}\mathrm{t}\colon\\!\dot{\mathbb{F}}\mathrm{t}+\boldsymbol{\Sigma}_{\ast}\colon\\!\mathbb{D}\mathrm{p}\geq
0\,,$

where the dissipative forces $\mathbb{X}\mathrm{t}$, conjugated to dissipative
rate $\dot{\mathbb{F}}\mathrm{t}$, are given by

(2.16)
$\displaystyle\mathbb{X}\mathrm{t}=\mathbb{F}\mathrm{e}^{\mathrm{t}}\\!\cdot\\!\mathbb{P}\\!\cdot\\!\mathbb{F}\mathrm{p}^{\mathrm{t}}-\rho_{\scalebox{0.65}{0}}\left.\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{t}}\right|_{\scriptsize\mathbb{F}\mathrm{e}}-\rho_{\scalebox{0.65}{0}}\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}+\rho_{\scalebox{0.65}{0}}\,\boldsymbol{\nabla}\\!\cdot\\!\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}\,.$

The relation (2.16) defines the thermodynamic displacive driving forces for
change in $\mathbb{F}\mathrm{t}$, 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. [110,
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{\mathbb{F}}\mathrm{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

(2.17)
$\mathbb{P}=\rho_{\scalebox{0.65}{0}}\left.\frac{\partial\psi_{\mathrm{e}}}{\partial\mathbb{F}\mathrm{e}}\right|_{\scriptsize\mathbb{F}\mathrm{t}}\\!\\!\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}\,.$

A quadratic form for the strain energy density per unit reference volume is
assumed, for which a dependence of $\psi_{\mathrm{e}}$ on
$\mathbb{F}\mathrm{e}$ and $\mathbb{F}\mathrm{t}$ manifests explicitly via the
anisotropic elastic components:

(2.18)
$\rho_{\scalebox{0.65}{0}}\psi_{\mathrm{e}}=\tfrac{1}{2}\mathbb{E}\mathrm{e}\colon\\!\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)\colon\\!\mathbb{E}\mathrm{e}\,,$

where $\mathbb{E}\mathrm{e}$ is the elastic Green-Lagrange strain tensor,
defined by

(2.19)
$\mathbb{E}\mathrm{e}=\tfrac{1}{2}\left(\mathbb{C}\mathrm{e}-\mathbb{I}\right)\,,$

with
$\mathbb{C}\mathrm{e}=\mathbb{F}\mathrm{e}^{\mathrm{t}}\cdot\mathbb{F}\mathrm{e}$
the right elastic Cauchy-Green deformation tensor, so that
$\mathbb{C}\mathrm{et}=\mathbb{F}\mathrm{t}^{\mathrm{t}}\\!\cdot\\!\mathbb{C}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}$.
Inserting eq. (2.18) into the hyperelasticity condition (2.17), the nonlinear
stress-elastic strain constitutive relation is rewritten as follows

(2.20)
$\displaystyle\mathbb{P}=\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{S}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}+\mathbb{F}\mathrm{et}\\!\cdot\\!\left(\mathbb{E}\mathrm{e}\colon\\!\frac{\partial\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)}{\partial\mathbb{C}\mathrm{et}}\colon\\!\mathbb{E}\mathrm{e}\right)\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}\,,$

where
$\mathbb{S}\mathrm{e}=\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)\colon\\!\mathbb{E}\mathrm{e}$
is an elastic stress measure associated with $\mathbb{E}\mathrm{e}$, and,
$\partial_{\mathbb{C}\mathrm{et}}\mathbb{D}$ is a sixth-order tensor, i.e.,
the derivative of $\mathbb{D}$ with respect of $\mathbb{C}\mathrm{et}$. 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 $\pazocal{D}_{\mathrm{t}}$ and
plastic $\pazocal{D}_{\mathrm{p}}$ dissipative rates per unit reference
volume, i.e.,

(2.21) $\pazocal{D}\doteq\pazocal{D}_{\mathrm{t}}+\pazocal{D}_{\mathrm{p}}\geq
0\,,$

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

(2.22)
$\pazocal{D}_{\mathrm{t}}=\mathbb{X}\mathrm{t}\colon\\!\dot{\mathbb{F}}\mathrm{t}\geq
0\mbox{~{}~{}and},~{}~{}\pazocal{D}_{\mathrm{p}}=\boldsymbol{\Sigma}_{\ast}\colon\\!\mathbb{D}\mathrm{p}\geq
0\,.$

Kinetic constitutive relations that relate the rates
$\dot{\mathbb{F}}\mathrm{t}$ and $\mathbb{D}\mathrm{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{\mathbb{F}}\mathrm{t}$ to the displacive driving forces
$\mathbb{X}\mathrm{t}$ is suggested, i.e.,

(2.23)
$\displaystyle\upsilon\,\dot{\mathbb{F}}\mathrm{t}=\mathbb{X}\mathrm{t}\,,$

where $\upsilon>0$ is a viscosity-like parameter. For example, the case with
$\nu\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 $\mathbb{X}\mathrm{t}$ gives
rise to thermodynamic consistency conditions for phase transformations, so
that the dissipation inequality in eq. (2.22) is unequivocally satisfied, as
follows

(2.24)
$\displaystyle\pazocal{D}_{\mathrm{t}}=\upsilon\,\lvert\mathbb{X}\mathrm{t}\lvert^{2}\geq
0\,,$

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

(2.25)
$\displaystyle\upsilon\,\dot{\mathbb{F}}\mathrm{t}=\mathbb{X}\mathrm{t}=\underbrace{\mathbb{C}\mathrm{e}\\!\cdot\\!\left(\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)\colon\\!\mathbb{E}\mathrm{e}\right)\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}}_{\mbox{\footnotesize
forces due to elastic
energy}}-\underbrace{~{}\rho_{\scalebox{0.65}{0}}\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}+\rho_{\scalebox{0.65}{0}}\,\boldsymbol{\nabla}\\!\cdot\\!\frac{\partial\psi_{{}_{\nabla}}}{\partial\boldsymbol{\nabla}\mathbb{F}\mathrm{t}}~{}}_{\mbox{\footnotesize
transformational forces}}\,,$

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 $\mathbb{X}\mathrm{t}=\mathbb{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

(2.26)
$\rho_{\scalebox{0.65}{0}}\,\psi_{{}_{\nabla}}=\tfrac{1}{2}\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\,\therefore{{}^{\scalebox{0.65}{6}}\boldsymbol{\Lambda}}\,\therefore\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\,,$

where ${{}^{\scalebox{0.65}{6}}\boldsymbol{\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 [165] such that
${{}^{\scalebox{0.65}{6}}\boldsymbol{\Lambda}}={\lambda\;^{\scalebox{0.65}{6}}\mathbb{I}}$,
with ${{}^{\scalebox{0.65}{6}}\mathbb{I}}$ the sixth-rank identity tensor, eq.
(2.26) reduces to

(2.27)
$\rho_{\scalebox{0.65}{0}}\,\psi_{{}_{\nabla}}=\tfrac{1}{2}\lambda\,\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\,\therefore{{}^{\scalebox{0.65}{6}}\mathbb{I}}\,\therefore\boldsymbol{\nabla}\mathbb{F}\mathrm{t}=\tfrac{1}{2}\lambda\,\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\,\therefore\boldsymbol{\nabla}\mathbb{F}\mathrm{t}\,,$

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

(2.28)
$\displaystyle\upsilon\,\dot{\mathbb{F}}\mathrm{t}=\mathbb{X}\mathrm{t}=\mathbb{C}\mathrm{e}\\!\cdot\\!\left(\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)\colon\\!\mathbb{E}\mathrm{e}\right)\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}-\rho_{\scalebox{0.65}{0}}\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}+\lambda\,\boldsymbol{\nabla}^{2}\mathbb{F}\mathrm{t}\,,$

with $\boldsymbol{\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 [136]. 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 $\mathbb{F}\mathrm{p}$, given in terms of
the plastic deformation rate $\mathbb{D}\mathrm{p}$, is determined by
considering the postulate of maximum dissipation [119]. The space of
admissible stresses $\pazocal{E}_{\boldsymbol{\sigma}}$ is written as

(2.29)
$\pazocal{E}_{\boldsymbol{\sigma}}=\\{\boldsymbol{\sigma}\mid\phi\left(\boldsymbol{\sigma}\right)<0\\}\,,$

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

(2.30)
$\boldsymbol{\sigma}=j^{-1}\,\mathbb{P}\\!\cdot\\!\mathbb{F}^{\mathrm{t}}=j^{-1}\,\mathbb{F}\mathrm{e}\\!\cdot\\!\mathbb{S}\mathrm{e}\\!\cdot\\!\mathbb{F}\mathrm{e}^{\mathrm{t}}+j^{-1}\,\mathbb{F}\mathrm{et}\\!\cdot\\!\left(\mathbb{E}\mathrm{e}\colon\\!\frac{\partial\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)}{\partial\mathbb{C}\mathrm{et}}\colon\\!\mathbb{E}\mathrm{e}\right)\\!\\!\cdot\\!\mathbb{F}\mathrm{et}^{\mathrm{t}}\,,$

according to eq. (2.20). The work-conjugate stress
$\boldsymbol{\Sigma}_{\ast}$ in eq. (2.12) may also be related to the Cauchy
stress tensor $\boldsymbol{\sigma}$ by

(2.31)
$\displaystyle\boldsymbol{\Sigma}_{\ast}=j\,\mathbb{F}\mathrm{et}^{\mathrm{t}}\\!\cdot\\!\boldsymbol{\sigma}\\!\cdot\\!\mathbb{F}\mathrm{et}^{-\mathrm{t}}=\mathbb{F}\mathrm{t}^{\mathrm{t}}\\!\cdot\\!\boldsymbol{\Sigma}\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}\,,$

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

(2.32)
$\mathbb{D}\mathrm{p}=\dot{\eta}\;\mathbb{F}\mathrm{et}^{-1}\\!\cdot\\!\frac{\partial\phi}{\partial\boldsymbol{\sigma}}\\!\cdot\\!\mathbb{F}\mathrm{et}=\dot{\eta}\;\mathbb{H}\,,$

with $\dot{\eta}\geq 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
$\mathbb{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.,

(2.33)
$\phi\left(\boldsymbol{\sigma}\right)=\sqrt{3\,J_{2}\left(\boldsymbol{\sigma}\right)}-\sigma_{\scalebox{0.65}{0}}\mbox{~{}~{}with},~{}~{}J_{2}=\tfrac{1}{2}\mathop{\mathrm{dev}}\boldsymbol{\sigma}\colon\\!\mathop{\mathrm{dev}}\boldsymbol{\sigma}\,,$

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

(2.34)
$\mathbb{D}\mathrm{p}=\tfrac{3}{2}\,\dot{\eta}\;\mathbb{F}\mathrm{et}^{-1}\\!\cdot\\!\frac{\mathop{\mathrm{dev}}\boldsymbol{\sigma}}{\sigma_{\scalebox{0.65}{0}}}\\!\cdot\\!\mathbb{F}\mathrm{et}\,,$

for which the dissipation inequality for plastic flow in eq. (2.22) with
(2.31) is satisfied, i.e.,

(2.35)
$\pazocal{D}_{\mathrm{p}}=\tfrac{3}{2}\,j\,\dot{\eta}\;\frac{\left\lvert\mathrm{dev}\boldsymbol{\sigma}\right\rvert^{2}}{\sigma_{\scalebox{0.65}{0}}}\geq
0\,.$

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 [183], as follows

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

which differs from the transformation path proposed in Ref. [48] by a rotation
of $\sim\pm 5.2^{\circ}$ around the $[0001]_{\textrm{hcp}}$ axis [276]. 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)_{\textrm{bcc}}$ plane, which consists
of an elongation and a compression along the $[1\bar{1}0]_{\textrm{bcc}}$ and
the $[001]_{\textrm{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 $\mathbb{U}$, i.e.,

(2.37)
$\mathbb{U}=\small\begin{bmatrix}{\color[rgb]{1,1,1}-}\dfrac{3}{4\sqrt{2}}+\dfrac{1}{4}\sqrt{\dfrac{3}{2}}\;c_{\mbox{\scriptsize/{a}}}&-\dfrac{3}{4\sqrt{2}}+\dfrac{1}{4}\sqrt{\dfrac{3}{2}}\;c_{\mbox{\scriptsize/{a}}}&\footnotesize
0\\\
-\dfrac{3}{4\sqrt{2}}+\dfrac{1}{4}\sqrt{\dfrac{3}{2}}\;c_{\mbox{\scriptsize/{a}}}&{\color[rgb]{1,1,1}-}\dfrac{3}{4\sqrt{2}}+\dfrac{1}{4}\sqrt{\dfrac{3}{2}}\;c_{\mbox{\scriptsize/{a}}}&0\\\
0&0&\dfrac{\sqrt{3}}{2}\end{bmatrix}\,,$

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

In the described case, the transformations $\mathbb{U}$ and t are illustrated
separately but can occur simultaneously, as discussed by using first-
principles simulations [79]. 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. [75], 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_{\mbox{\scriptsize/{a}}}$ is experimentally chosen to determine eq.
(2.37), the corresponding homogeneous mapping $\mathbb{U}$ contains obviously
and inseparably both elastic and irreversible part of the deformation in
samples. A homogeneous distortion $\mathbb{U}\mathrm{t}$ is therefore
introduced to identify the pure transformational component of the total
deformation provided by experimental data under high hydrostatic pressure,
i.e.,

(2.38) $\mathbb{U}\mathrm{t}=\kappa\mathbb{U}\,,$

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

Figure 2.2: Crystallographic relations in the bcc-to-hcp martensitic phase
transition established in Refs. [183, 23]. (a) Red atoms in a bcc atomic-side
unit cell are located at a $(110)_{\textrm{bcc}}$ layer and the blue atoms at
the adjacent layers. (b) The transition path consists of two transformations.
First, a shear deformation $\mathbb{U}$ leads to an elongation and a
compression along the $[1\bar{1}0]_{\textrm{bcc}}$ and the
$[001]_{\textrm{bcc}}$ directions, respectively. The deformation transforms a
polygon in blue into a regular hexagon in red, corresponding to the
$(0001)_{\textrm{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 $[1\bar{1}0]_{\textrm{bcc}}$ direction.

##### Multiple symmetry-related variants

During the forward $\alpha\to\epsilon$ and the reverse
$\epsilon\to\alpha^{\prime}$ 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^{\prime}$ 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 $\prescript{\epsilon n}{\alpha}{\mathbb{U}\mathrm{t}}$ from the
linear mapping $\mathbb{U}\mathrm{t}$ is described by

(2.39) $\prescript{\epsilon
n}{\alpha}{\mathbb{U}\mathrm{t}}=\mathbb{R}_{\textrm{bcc}}^{\mathrm{t}}\\!\cdot\\!\mathbb{U}\mathrm{t}\\!\cdot\\!\mathbb{R}_{\textrm{bcc}}\,,$

where $\mathbb{R}_{\textrm{bcc}}$ is a rotation matrix in the point group of
cubic lattice $\prescript{n\\!}{}{\pazocal{H}_{\textrm{bcc}}}$ and $n$ the
number of hcp variants [210]. 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^{\prime}$ transitions, the bcc variants
$\prescript{\alpha^{\prime}\\!m}{\epsilon n}{\mathbb{U}\mathrm{t}}$ are
deduced by performing the following operation:

(2.40) $\prescript{\alpha^{\prime}\\!m}{\epsilon
n}{\mathbb{U}\mathrm{t}}=\mathbb{R}_{\textrm{bcc}}^{\mathrm{t}}\\!\cdot\\!\mathbb{R}_{\textrm{hcp}}^{\mathrm{t}}\\!\cdot\\!\mathbb{U}\mathrm{t}^{-1}\\!\cdot\\!\mathbb{R}_{\textrm{hcp}}\\!\cdot\\!\mathbb{U}\mathrm{t}\\!\cdot\\!\mathbb{R}_{\textrm{bcc}}\,,$

where $\mathbb{R}_{\textrm{hcp}}$ is a rotation matrix in the point group of
hexagonal lattice $\prescript{m\\!}{}{\pazocal{H}_{\textrm{hcp}}}$ and $m$ the
number of bcc variants [210]. 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 19${}^{\textrm{th}}$ 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
$\prescript{m\\!}{}{\pazocal{H}_{\textrm{bcc}}}$ and
$\prescript{n\\!}{}{\pazocal{H}_{\textrm{hcp}}}$ are given in Tab. 1 from Ref.
[259], 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 ${{}^{k}\mathbb{U}\mathrm{t}}$ with
$k=1,\dots,18$, i.e.,

(2.41) ${{}^{k}\mathbb{U}\mathrm{t}}=\left\\{\begin{matrix}\begin{aligned}
\prescript{\epsilon n}{\alpha}{\mathbb{U}\mathrm{t}}&~{}~{}~{}~{}~{}~{}1\leq
k\leq 6\\\ \prescript{\alpha^{\prime}\\!m}{\epsilon
n}{\mathbb{U}\mathrm{t}}&~{}~{}~{}~{}~{}~{}7\leq k\leq
18\,,\end{aligned}\end{matrix}\right.$

which are associated with the variant of interest $V_{k}$ for the forward
$(1\leq k\leq 6)$ and the reverse $(7\leq k\leq 18)$ transformations.

Figure 2.3: (a) Schematic illustration of the multiple symmetry-related
variants for the forward $\alpha\rightarrow\epsilon$ (in red) and the reverse
$\epsilon\rightarrow\alpha^{\prime}$ (blue) phase transitions in iron. (b) The
corresponding reaction pathway network in a specific
$\\{\mathbb{C}_{1},\mathbb{C}_{2},\mathbb{C}_{3}\\}$ strain space, within
which the transformational Cauchy-Green tensor
$\mathbb{C}\mathrm{t}=\mathbb{F}\mathrm{t}^{\mathrm{t}}\cdot\mathbb{F}\mathrm{t}$
as well as some quantities described in the text, are defined.

##### Reaction pathways in strain spaces

Instead of introducing the Landau thermodynamic potential [164], 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 [74, 75] 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_{\mathrm{t}}}_{k}}$
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.,

(2.42)
$\displaystyle{{\hat{\psi_{\mathrm{t}}}_{k}}}\doteq{{\hat{\psi_{\mathrm{t}}}_{k}}}\big{(}{{}^{k}\hat{\mathbb{C}}\mathrm{t}}\big{)}\,,$

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

(2.43)
${{}^{k}\hat{\mathbb{C}}\mathrm{t}}={{}^{k}\mathbb{U}\mathrm{t}^{\mathrm{t}}}\cdot{{}^{k}\mathbb{U}\mathrm{t}}\,,$

as listed in Appendix A from Ref. [259], 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{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}$ and ending
${{}^{k}\mathbb{C}\mathrm{t}}_{\scriptsize\mbox{end}}$ strain states, as
follows

(2.44)
${{}^{k}\hat{\mathbb{C}}\mathrm{t}}\left({s_{k}}\right)=\left({1-s_{k}}\right){{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}+{s_{k}}\,{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{end}}\,,$

with ${s_{k}}\in[0,1]$ the curvilinear coordinates along $k$. For instance,
hcp variants $V_{k}$ are parameterized by:
${{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}=\mathbb{I}$ and
${{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{end}}={{}^{k}\mathbb{U}\mathrm{t}^{2}}$,
with $1\leq k\leq 6$. Generating the reaction pathways with eqs. (2.37$-$2.44)
and using projection matrices $\mathbb{C}_{1}$, $\mathbb{C}_{2}$ and
$\mathbb{C}_{3}$, an example of three-dimensional 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
$\\{\mathbb{C}_{1},\mathbb{C}_{2},\mathbb{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:

(2.45) $\displaystyle\mathbb{C}_{1}$
$\displaystyle=\small\begin{bmatrix}1&{\color[rgb]{1,1,1}-}1&0\\\ 1&-3&0\\\
0&{\color[rgb]{1,1,1}-}0&0\end{bmatrix}\,,~{}~{}\mathbb{C}_{2}=\small\begin{bmatrix}0&0&{\color[rgb]{1,1,1}-}0\\\
0&1&{\color[rgb]{1,1,1}-}1\\\ 0&1&-3\end{bmatrix}\,,~{}~{}\mathbb{C}_{3}$
$\displaystyle=\small\begin{bmatrix}-3&0&1\\\ {\color[rgb]{1,1,1}-}0&0&0\\\
{\color[rgb]{1,1,1}-}1&0&1\end{bmatrix}\,.$

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_{\mathrm{t}}$ in
a whole strain space, e.g. not only restricted along the pathways as
${{\hat{\psi_{\mathrm{t}}}_{k}}}$, the partition of unity approach is used as
a weighted sum of the contribution ${\psi_{\mathrm{t}}}_{k}$ of each
individual pathway $k$. Thus, the overall inelastic energy density
$\psi_{\mathrm{t}}$ is formally defined by introducing the weighting functions
$\omega_{k}\left(\mathbb{C}\mathrm{t}\right)$, i.e.,

(2.46)
$\psi_{\mathrm{t}}\left(\mathbb{C}\mathrm{t}\right)=\sum_{k=1}^{18}\omega_{k}\left(\mathbb{C}\mathrm{t}\right)\,{\psi_{\mathrm{t}}}_{k}\left(\mathbb{C}\mathrm{t}\right)\,,$

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

(2.47)
$\sum_{k=1}^{18}\omega_{k}\left(\mathbb{A}\right)=1\mbox{~{}~{}with},~{}~{}\omega_{k}\left(\mathbb{A}\right)=\frac{d_{k}^{\,-h}\left(\mathbb{A}\right)}{\displaystyle\sum\nolimits_{i=1}^{18}d_{i}^{\,-h}\left(\mathbb{A}\right)}\,,$

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

(2.48)
$d_{k}\left(\mathbb{A}\right)=\lvert{{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{A}\right)}\lvert=\min_{{\zeta_{k}}\in[0,1]}\;\lvert\mathbb{A}-{{}^{k}\hat{\mathbb{A}}}\left({\zeta_{k}}\right)\lvert\,,$

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

(2.49)
${\zeta^{\infty}_{k}}\left(\mathbb{A}\right)={{}^{k}{\hat{\mathbb{D}}}}\,\colon\\!\big{(}\mathbb{A}-{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}\big{)}=\frac{{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{end}}-{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}}{\lvert{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{end}}-{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}\lvert}\,\colon\\!\big{(}\mathbb{A}-{{}^{k}\hat{\mathbb{C}}\mathrm{t}}_{\scriptsize\mbox{start}}\big{)}\,,$

where ${{}^{k}{\hat{\mathbb{D}}}}$ defines the normalized direction of the
pathway $k$, the $\mathrm{argmin}~{}{\zeta_{k}}$ in eq. (2.48) is also
determined by solving $\partial_{\zeta_{k}}\,d_{k}\left(\mathbb{A}\right)=0$
for a given $\mathbb{C}\mathrm{t}$, leading to

(2.50) $\zeta_{k}\left(\mathbb{A}\right)=\left\\{\begin{matrix}\begin{aligned}
&{\zeta^{\infty}_{k}}\left(\mathbb{A}\right)&&\mbox{if}:~{}~{}~{}{\zeta^{\infty}_{k}}\left(\mathbb{A}\right)\in[0,1]\\\
&0&&\mbox{if}:~{}~{}~{}{\zeta^{\infty}_{k}}<0\\\
&1&&\mbox{if}:~{}~{}~{}{\zeta^{\infty}_{k}}>1\,,\end{aligned}\end{matrix}\right.$

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

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

(2.51)
${\psi_{\mathrm{t}}}_{k}\left(\mathbb{C}\mathrm{t}\right)={{\hat{\psi_{\mathrm{t}}}_{k}}}\left({\zeta_{k}}\left(\mathbb{C}\mathrm{t}\right)\right)+\underbrace{\sigma\,d_{k}(\mathbb{C}\mathrm{t})+\pi\,|\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{C}\mathrm{t}\right)}|\;}_{\mbox{\footnotesize
out-of-path component}}\,,$

such that
$\partial_{\mathbb{C}\mathrm{t}}\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}}\left(\mathbb{C}\mathrm{t}\right)$
and ${{}^{k}{\hat{\mathbb{D}}}}$ are orthogonal to each other, i.e.,
$\partial_{\mathbb{C}\mathrm{t}}\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}}\left(\mathbb{C}\mathrm{t}\right)\colon{{}^{k}{\hat{\mathbb{D}}}}=0$.
Here, $\mathrm{tr}\,\mathbb{A}$ denotes the trace of $\mathbb{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_{\mathrm{t}}$ defined by eq. (2.46) with eq. (2.51), for all
$\mathbb{C}\mathrm{t}$ 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_{\mathrm{t}}}_{k}}}$ along all reaction coordinates
${\zeta_{k}}\left(\mathbb{C}\mathrm{t}\right)$ of the individual pathways $k$.
Then, the weighting functions $\omega_{k}\left(\mathbb{C}\mathrm{t}\right)$
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}\left(\mathbb{C}\mathrm{t}\right){{\hat{\psi_{\mathrm{t}}}_{k}}}$
in the $\\{\mathbb{C}_{1},\mathbb{C}_{2},\mathbb{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
$\mathbb{C}_{1}+\mathbb{C}_{2}+\mathbb{C}_{3}$. Figure (2.4c) shows a
$10^{9}\;$J.m-3-iso-volume related to the out-of-path contribution $\sigma
d_{k}\left(\mathbb{C}\mathrm{t}\right)$ 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 the total inelastic energy $\psi_{\mathrm{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).

Figure 2.4: Construction of the total inelastic energy landscape
$\psi_{\mathrm{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}\;$J.m-3 (in white, for hcp phases). (b) Extrapolation of the minimum
energy potential in the whole
$\\{\mathbb{C}_{1},\mathbb{C}_{2},\mathbb{C}_{3}\\}$ strain space, e.g.
$5\times 10^{8}\;$J.m-3-iso-surface. (c) shows a $10^{9}\;$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_{\mathrm{t}}$ landscape, respectively.

##### 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_{\mathrm{t}}$ with
respect to $\mathbb{F}\mathrm{t}$, which can be expressed as follows

(2.52)
$\displaystyle\frac{\partial\psi_{\mathrm{t}}}{\partial\mathbb{F}\mathrm{t}}=2\,\mathbb{F}\mathrm{t}\\!\cdot\\!\frac{\partial\psi_{\mathrm{t}}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}\,.$

According to eq. (2.46), the derivative of the energy function in the right-
hand side of eq. (2.52) yields

(2.53)
$\displaystyle\frac{\partial\psi_{\mathrm{t}}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}=\sum_{k=1}^{18}{\psi_{\mathrm{t}}}_{k}\left(\mathbb{C}\mathrm{t}\right)\frac{\partial\omega_{k}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}+\omega_{k}\left(\mathbb{C}\mathrm{t}\right)\frac{\partial{\psi_{\mathrm{t}}}_{k}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}\,,$

where the derivative of the weighting functions
$\omega_{k}\left(\mathbb{A}\right)$ with respect to $\mathbb{A}$ is given,
without loss of generality, by

(2.54)
$\displaystyle\frac{\partial\omega_{k}\left(\mathbb{A}\right)}{\partial\mathbb{A}}=h\,\sum_{i=1}^{18}\frac{\omega_{i}\left(\mathbb{A}\right)}{d_{i}\left(\mathbb{A}\right)}\left(\omega_{k}\left(\mathbb{A}\right)-\delta_{ik}\right)\,{{}^{i}\mathbb{N}}\left(\mathbb{A}\right)\,,$

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

(2.55) $\displaystyle{{}^{i}\mathbb{N}}\left(\mathbb{A}\right)=\dfrac{\partial
d_{i}\left(\mathbb{A}\right)}{\partial\mathbb{A}}=\dfrac{{{}^{i}{\boldsymbol{\Pi}}}\left(\mathbb{A}\right)}{d_{i}\left(\mathbb{A}\right)}\,,$

such that $\lvert{{}^{i}\mathbb{N}}\left(\mathbb{A}\right)\lvert\;=1$, and,
${{}^{i}\mathbb{N}}\left(\mathbb{A}\right)\colon\\!\\!{{}^{i}\hat{\mathbb{D}}}=0$
when ${\zeta^{\infty}_{k}}\left(\mathbb{A}\right)\in[0,1]$. Moreover, the
derivative of ${\psi_{\mathrm{t}}}_{k}$ with respect to $\mathbb{C}\mathrm{t}$
in eq. (2.53) leads to

(2.56)
$\displaystyle\frac{\partial{\psi_{\mathrm{t}}}_{k}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}=$
$\displaystyle\frac{\partial{{\hat{\psi_{\mathrm{t}}}_{k}}}\left({\zeta_{k}}\left(\mathbb{C}\mathrm{t}\right)\right)}{\partial{\zeta_{k}}}{{}^{k}{\hat{\mathbb{D}}}}+\sigma\;{{}^{k}{\mathbb{N}}\left(\mathbb{C}\mathrm{t}\right)}+\pi\,\mathrm{sgn}\big{(}\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{C}\mathrm{t}\right)}\big{)}\,\big{(}\mathbb{I}-{{}^{k}{\hat{\mathbb{D}}}}~{}\mathrm{tr}\,{{}^{k}{\hat{\mathbb{D}}}}\big{)}\,.$

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{\mathbb{D}}}}$ along the reaction
pathways, while the second component is associated with the normal directions
${{}^{k}{\mathbb{N}}\left(\mathbb{C}\mathrm{t}\right)}$ towards
$\mathbb{C}\mathrm{t}$.

##### 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 $\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)$ in eq. (2.28), which
depends on the elastic and transformational deformation states, is given in
the whole strain space by

(2.57)
$\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)=\sum_{k=1}^{18}\omega_{k}\left(\mathbb{C}\mathrm{et}\right)\;{{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)\,,$

where ${{}^{k}\mathbb{D}}$ are the elasticity tensors associated with the
reaction pathways $k$, and, $\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)$ are
the reaction coordinates that minimize the Euclidean distance
$d_{k}\left(\mathbb{C}\mathrm{et}\right)$ between $\mathbb{C}\mathrm{et}$ and
the individual paths $k$. The weighting functions
$\omega_{k}\left(\mathbb{C}\mathrm{et}\right)$ are also defined by eq. (2.47),
where the partition of unity is written as a function of
$\mathbb{C}\mathrm{et}$. The projected tensors
${{}^{k}\hat{\mathbb{C}}\mathrm{et}}$ are also mapped onto the reaction
pathways (as well as ${{}^{k}\hat{\mathbb{C}}\mathrm{t}}$) and the
corresponding reaction coordinates are consistently determined by solving
$\partial_{\zeta_{k}}\,d_{k}\left(\mathbb{C}\mathrm{et}\right)=0$. Imposing
$\partial_{\mathbb{C}\mathrm{et}}\mathbb{D}\left(\mathbb{C}\mathrm{et}\right)=\mathbb{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
$\mathbb{C}\mathrm{et}={{}^{k}\hat{\mathbb{C}}\mathrm{et}}$, the elasticity
tensors ${{}^{k}\mathbb{D}}$ in eq. (2.57) may be represented by a cubic
interpolation function to ensure numerical stability, i.e.,

(2.58)
${{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)=(1-3{\zeta_{k}}^{2}+2{\zeta_{k}}^{3})\;{\mathbb{D}^{\alpha}}+(3{\zeta_{k}}^{2}-2{\zeta_{k}}^{3})\;{\mathbb{D}^{\epsilon}}\,,$

with ${\mathbb{D}^{\alpha}}$ and ${\mathbb{D}^{\epsilon}}$ the elastic
stiffness tensors of the pure bcc and hcp iron phases, respectively. In
particular, if ${\zeta^{\infty}_{k}}\left(\mathbb{C}\mathrm{et}\right)<0$
($>1$), also
${{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)={\mathbb{D}^{\alpha}}$
($={\mathbb{D}^{\epsilon}}$). For instance, for such pure hcp $\epsilon$-Fe
phases, the finite hyperelasticity condition from eq. (2.18) leads therefore
to

(2.59)
${\mathbb{D}^{\epsilon}}=\rho_{\scalebox{0.65}{0}}\left.\\!\dfrac{\partial^{2}\psi_{\mathrm{e}}}{\partial\mathbb{E}\mathrm{e}\;\partial\mathbb{E}\mathrm{e}}\right|_{\epsilon}=\left.\\!\frac{\partial\mathbb{S}\mathrm{e}}{\partial\mathbb{E}\mathrm{e}}\right|_{\epsilon}\,,$

where ${\mathbb{D}^{\epsilon}}$ is defined in the reference configuration
$\Omega_{\scalebox{0.65}{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 ${\text{{b}}^{\epsilon}}$, expressed
in the current and deformed $\Omega$ by

(2.60)
$\displaystyle{\text{{b}}^{\epsilon}}=\left.\\!\frac{\partial\boldsymbol{\sigma}}{\partial\boldsymbol{\varepsilon}}\right|_{\epsilon}\,,$

with $\boldsymbol{\varepsilon}$ the Eulerian strain tensor [274, 237, 63]. In
the present work, the relevant tensor $\text{{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 $\mathbb{F}_{\scriptsize\mbox{vol}}$ from initial
volume $V_{\scalebox{0.65}{0}}$ to the final volume
$V=j\,V_{\scalebox{0.65}{0}}$ with
$\mathbb{F}_{\scriptsize\mbox{vol}}=j^{1/3}\mathbb{I}$, and then, to a small
and symmetric shear isochoric deformation
$\mathbb{F}_{\scriptsize\mbox{iso}}=\mathbb{I}+\boldsymbol{\varepsilon}$ with
$\lvert\boldsymbol{\varepsilon}\lvert\ll 1$, such as it is commonly performed
using density functional theory calculations [167, 147]. Without plasticity,
the total deformation gradient is also given by
$\mathbb{F}=j^{1/3}\left(\mathbb{I}+\boldsymbol{\varepsilon}\right)$. Using
eq. (2.60) with the aid of eq. (2.30) and considering
$\boldsymbol{\sigma}=p\,\mathbb{I}$ with the Cauchy pressure $p\leq 0$, the
relation between the elasticity tensors $\mathbb{D}^{\epsilon}$ and the
incremental tangent modulus $\text{{b}}^{\epsilon}$ is reduced to

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

exhibiting the same symmetries as in eq. (2.59). Thus, the elasticity tensors
${\mathbb{D}^{\epsilon}}$ is obtained by identifying the values of
${\text{{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}={\text{{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 u (3 DoFs) and the non-symmetric transformational
distortion $\mathbb{F}\mathrm{t}$ (9 DoFs) at each reference point X in
$\Omega_{\scalebox{0.65}{0}}$, by solving the system of partial differential
equations, as follows

(2.62) $\left\\{\begin{matrix}\begin{aligned}
\rho_{\scalebox{0.65}{0}}\,\ddot{\textbf{{u}}}=&\sum_{k=1}^{18}\boldsymbol{\nabla}\\!\cdot\\!\hskip
1.28032pt\Big{\\{}\omega_{k}\left(\mathbb{C}\mathrm{et}\right)\;\mathbb{F}\mathrm{e}\\!\cdot\\!\big{(}{{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)\colon\\!\mathbb{E}\mathrm{e}\big{)}\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}\\\
+&\,\omega_{k}\left(\mathbb{C}\mathrm{et}\right)\big{(}\mathbb{E}\mathrm{e}\colon\\!{{}^{k}\mathbb{D}^{\prime}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)\colon\\!\mathbb{E}\mathrm{e}\big{)}\;\mathbb{F}\mathrm{et}\\!\cdot\\!{{}^{k}{\hat{\mathbb{D}}}}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}+\big{(}\mathbb{E}\mathrm{e}\colon\\!{{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)\colon\\!\mathbb{E}\mathrm{e}\big{)}\;\mathbb{F}\mathrm{et}\\!\cdot\\!\frac{\partial\omega_{k}\left(\mathbb{C}\mathrm{et}\right)}{\partial\mathbb{C}\mathrm{et}}\\!\cdot\\!\mathbb{F}\mathrm{p}^{-\mathrm{t}}\Big{\\}}\\\
\upsilon\dot{\mathbb{F}}\mathrm{t}=&\sum_{k=1}^{18}\omega_{k}\left(\mathbb{C}\mathrm{et}\right)\;\mathbb{C}\mathrm{e}\\!\cdot\\!\big{(}{{}^{k}\mathbb{D}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{et}\right)\right)\colon\\!\mathbb{E}\mathrm{e}\big{)}\\!\cdot\\!\mathbb{F}\mathrm{t}^{-\mathrm{t}}+\lambda\,\boldsymbol{\nabla}^{2}\mathbb{F}\mathrm{t}\\\
-&\,2\rho_{\scalebox{0.65}{0}}\big{(}{{\hat{\psi_{\mathrm{t}}}_{k}}}\left({\zeta_{k}}\left(\mathbb{C}\mathrm{t}\right)\right)+\sigma\,d_{k}(\mathbb{C}\mathrm{t})+\pi\,|\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{C}\mathrm{t}\right)}|\big{)}\,\mathbb{F}\mathrm{t}\\!\cdot\\!\frac{\partial\omega_{k}\left(\mathbb{C}\mathrm{t}\right)}{\partial\mathbb{C}\mathrm{t}}\\\
-&\,2\rho_{\scalebox{0.65}{0}}\,\omega_{k}\left(\mathbb{C}\mathrm{t}\right)\Big{(}{\hat{\psi_{\mathrm{t}}}_{k}}^{\prime}\left({\zeta_{k}}\left(\mathbb{C}\mathrm{t}\right)\right)\;\mathbb{F}\mathrm{t}\\!\cdot\\!{{}^{k}{\hat{\mathbb{D}}}}+\sigma\;\mathbb{F}\mathrm{t}\\!\cdot\\!{{}^{k}{\mathbb{N}}\left(\mathbb{C}\mathrm{t}\right)}+\pi\,\mathrm{sgn}\big{(}\mathop{\mathrm{tr}}\,{{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{C}\mathrm{t}\right)}\big{)}\,\big{(}\mathbb{F}\mathrm{t}-\mathbb{F}\mathrm{t}\\!\cdot\\!{{}^{k}{\hat{\mathbb{D}}}}~{}\mathrm{tr}\,{{}^{k}{\hat{\mathbb{D}}}}\big{)}\Big{)}\,,\end{aligned}\end{matrix}\right.$

where ′ denotes the derivative with respect to the reaction coordinates
${\zeta_{k}}$, while the derivatives of the weighting functions $\omega_{k}$
with respect to $\mathbb{C}\mathrm{t}$ and $\mathbb{C}\mathrm{et}$ 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. [259].

### 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. [259] 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\\}$
$\mathring{\mbox{A}}$3/at, are selected from Ref. [76]. In accordance with
these experimental measures, the associated elastic components
$\text{{b}}^{\alpha}$ and $\text{{b}}^{\epsilon}$ for both pure bcc and hcp
phases are given in Ref. [167], while the stiffness tensor
$\mathbb{D}^{\epsilon}$ is expressed in the current configuration by using eq.
(2.61), and, $\mathbb{D}^{\alpha}=\text{{b}}^{\alpha}$ at zero pressure.

The ratio $c_{\mbox{\scriptsize/{a}}}=1.603$ of the hcp close-packed structure
has been experimentally determined in Ref. [183], so that
$\mathop{\mathrm{det}}\mathbb{U}=9c_{\mbox{\scriptsize/{a}}}\,/16=0.902$.
However, $\mathbb{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.
[259], the transformational part $\mathbb{U}\mathrm{t}$ is related to
$\mathbb{U}$ as follows

(2.63)
$\mathbb{U}\mathrm{t}=\kappa\mathbb{U}=\sqrt{2}\left(1+\sqrt{1+\frac{8}{3}\frac{j_{\mbox{\scriptsize
exp}}\,p^{\epsilon}}{D^{\epsilon}}}\right)^{\\!\\!\\!-1/2}\mathbb{U}\,,$

where $D^{\epsilon}$ is the hcp bulk modulus, and, $j_{\mbox{\scriptsize
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_{\scalebox{0.65}{0}}=0.25$ GPa,
which is fairly of the same order of magnitude with Hugoniot elastic limits in
Ref. [218].

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 [74], 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 $\mathbb{C}\mathrm{t}$ with
${{}^{k}{\boldsymbol{\Pi}}\left(\mathbb{C}\mathrm{t}\right)}\neq\mathbb{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_{\mathrm{t}}}_{k}}}$ as a function of the reaction coordinate
$\zeta_{k}$ along each individual pathway $k$, i.e.,

(2.64)
$\rho_{\scalebox{0.65}{0}}{{\hat{\psi_{\mathrm{t}}}_{k}}}\left(\zeta_{k}\left(\mathbb{C}\mathrm{t}\right)\right)=\tfrac{1}{2}c_{1}\,\zeta_{k}^{2}+c_{2}\,\zeta_{k}\,,$

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. [259], these
parameters are given by

(2.65) $\displaystyle c_{2}$
$\displaystyle=\tfrac{1}{2}j^{\alpha\to\epsilon}p^{\alpha\to\epsilon}\,\mathrm{tr}\,{{}^{k}{\hat{\mathbb{D}}}}\,,~{}~{}\mbox{and},~{}~{}c_{1}$
$\displaystyle=\tfrac{1}{2}j^{\epsilon\to\alpha}p^{\epsilon\to\alpha}\,({{}^{k}\mathbb{U}\mathrm{t}^{-2}}\colon\\!{{}^{k}{\hat{\mathbb{D}}}})-c_{2}\,,$

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. [76], the forward transition starts at
$p^{\alpha\to\epsilon}=-14.9$ GPa, with the corresponding volume
$v^{\alpha\to\epsilon}=11.0$ $\mathring{\mbox{A}}$3/at, and, the reverse at
$p^{\epsilon\to\alpha}=-12.0$ GPa, with $v^{\epsilon\to\alpha}=10.6$
$\mathring{\mbox{A}}$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. [259], is
depicted in Fig. (2.4a).

The parameter $\upsilon$ 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 Fe3Ni, an attempt to fit the magnitude $\upsilon=14$
mPa.s, comparable to the viscosity of liquid metals, has been investigated by
using molecular dynamics simulations [74]. 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.,
$\upsilon\dot{\varepsilon}_{\mathrm{t}}<\sigma_{\scalebox{0.65}{0}}$, where
$\dot{\varepsilon}_{\mathrm{t}}$ is a measure of the transformational strain
rate. This measure can be estimated by
$\dot{\varepsilon}_{\mathrm{t}}=\varepsilon_{\mathrm{t}}/\Delta
t=\frac{1}{2}\lvert\mathbb{C}\mathrm{t}-\mathbb{I}\lvert/\Delta t$ during a
time interval $\Delta t$ awaited for the transformation, with
$\varepsilon_{\mathrm{t}}$ the norm of the transformational Green-Lagrange
deformation tensor. Thus, it follows that
$\upsilon<\sigma_{\scalebox{0.65}{0}}t_{\\!f}/\varepsilon_{\mathrm{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
$\upsilon\approx\sigma_{\scalebox{0.65}{0}}t_{\\!f}/\varepsilon_{\mathrm{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^{\ast}/\ell^{2}$, where $\lambda^{\ast}=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~{}{\mu}$m3. 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^{\alpha}_{11}/\rho_{\scalebox{0.65}{0}}}$. It therefore follows
that: $c_{L}=5850~{}$m.s-1, and, $t_{\\!f}\approx 1.7~{}\mu$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_{\scalebox{0.65}{0}}=0.86$, for
which the volume change is achieved within a time step from
$t_{\scalebox{0.65}{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_{\scalebox{0.65}{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\to\epsilon$ transitions.
However, the reverse $\epsilon\to\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 [76]
and $-23.0$ GPa [243] at room temperature. The experimental measurements from
Refs. [101, 76] have been plotted in Fig. (2.5) with symbols, where the more
recent data in Ref. [76] 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 [278]. 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 [101, 76].

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.

Figure (2.6) illustrates the partitioning of the total energy $\psi$ in terms
of elastic $\psi_{\mathrm{e}}/\psi$ (in blue) and inelastic
$(\psi_{\mathrm{t}}+\psi_{{}_{\nabla}})/\psi$ (green) energy ratios as a
function of the dimensionless simulation time $t^{\ast}=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^{\ast}\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 $\ast$
in Fig. (2.6a) when the internal elastic stored energy increases then to
$t^{\ast}=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 2.6: Partitioning of the total energy into the elastic and inelastic
components as a function of the dimensionless simulation time $t^{\ast}$, for
calculations (a) without and (b) with plasticity.

####  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^{\ast}$ 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 $\alpha$-Fe and $\epsilon$-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^{\ast}=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^{\ast}\approx 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.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^{\ast}$, for calculations without and with plasticity.

Figure (2.7b) shows the volume fractions of each variant $V_{k}$ as a function
of the simulation time $t^{\ast}$. 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 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^{\prime}$-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\to\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^{\prime}$-Fe
nucleation of variants occurs at $t^{\ast}\approx 0.90$, with $\sim$2% volume
fraction for each $\\{V_{12},V_{13},V_{15}\\}$, $\sim$1% for each
$\\{V_{11},V_{14},V_{16}\\}$, and, $\sim$0.5% for each of the 6 other bcc
variants. Thus, both initial $\alpha$-Fe and reversed $\alpha^{\prime}$-Fe
phases coexist at $t^{\ast}=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^{\prime}$ inclusions
is $\sim 12.3\%$ in the final microstructure.

To summarize, Fig. (2.8) illustrates the microstructure evolution under
hydrostatic pressure at $t^{\ast}=0.6$ and $t^{\ast}=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
$\sim$26.6% volume fraction of bcc phase at $t^{\ast}=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 [141]. 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 [279], lath-like $\epsilon$-Fe regions [54], and, ellipsoidal
$\epsilon$-Fe particles [205], 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^{\ast}=1.0$,
while the calculation with plasticity leads to 12 reversed bcc
$\alpha^{\prime}$-Fe, heterogeneously nucleated in pairs (e.g.
$\\{V_{11},V_{12}\\}$, in light and dark green) from one single $\epsilon$-Fe
variant.

Figure 2.8: Transformational states defined in both strain and current mesh
spaces at (a) $t^{\ast}=0.6$ and (b) $t^{\ast}=1.0$, for calculations without
and with plasticity. Each black dot in the strain space represents the current
transformational strain $\mathbb{C}\mathrm{t}$ 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 $\alpha$-Fe phase remains in a large fraction ($\sim$26.6%, in gray) at
$t^{\ast}=0.6$, whereas the calculation with plasticity exhibits a
polycrystalline iron formed by the 6 hcp $\epsilon$-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^{\ast}=1.0$, while the
latter shows the presence of 12 reversed bcc $\alpha^{\prime}$-Fe with
$\sim$12.3% volume fraction.

### 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 $\textbf{{z}}\parallel[001]_{\textrm{bcc}}$ direction, using $80\times
80\times 1280$ element-free Galerkin nodes ($\sim 8.2$ millions), with
periodic boundary conditions transverse to the direction of shock front
propagation, i.e., to $\textbf{{x}}\parallel[100]_{\textrm{bcc}}$ and
$\textbf{{y}}\parallel[010]_{\textrm{bcc}}$. The initial shock compression is
induced by imposing a velocity of $850$ m.s-1 on the rear face along
$\textbf{{z}}\parallel[001]_{\textrm{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_{\scalebox{0.65}{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]_{\textrm{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^{\alpha}_{11}/\rho_{\scalebox{0.65}{0}}}$, with
$b^{\alpha}_{11}=271~{}$GPa the corresponding low-pressure elastic component
of the pure bcc iron [167]. It therefore follows that: $c_{L}=5850~{}$m.s-1,
so that $t_{\\!f}\approx 0.55~{}\mu$s, which corresponds to the duration of
all calculations. For convenience, a dimensionless time $t^{\ast}$ is defined
as $t^{\ast}=t/t_{\\!f}$, while the dimensionless length $L^{\ast}$ along z is
given by $L^{\ast}=z/L_{z}$, so that both quantities $t^{\ast}$ and $L^{\ast}$
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.

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]_{\textrm{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.

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^{\ast}=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 deformation
microstructure are elucidated in the following sections.

Figure 2.10:  (a) The three-rank network of reaction pathways is projected in
a $\\{\mathbb{C}_{1},\mathbb{C}_{2},\mathbb{C}_{3}\\}$ strain space, within
which the local transformational Cauchy-Green $\mathbb{C}\mathrm{t}$ 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.

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]_{\textrm{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 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
noncoherent grain boundaries (GBs), are observed in Refs. [141, 142].

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

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]_{\textrm{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.

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_{\textrm{vM}}$ in Fig. (2.11c) with
respect to $t^{\ast}$, obtained without (gray curves) and with (black curves)
plasticity. Both averaged quantities $p$ and $\sigma_{\textrm{vM}}$ are
displayed with respect to $L^{\ast}$ along the
$\textbf{{z}}\parallel[001]_{\textrm{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$ m.s-1, while the phase-transition
front produces a considerable increase of the velocity at free surface, up to
$v_{z}=1660$ 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, 134]. The instants when both P and rEp waves reach the free
surface are displayed by the double-headed arrows in Fig. (2.11a),
corresponding to $v_{z}=880$ m.s-1 and $v_{z}=1170$ 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
$[00\bar{1}]_{\textrm{bcc}}$ direction) produce a residual stress state that
does not favor the mandatory forward $\alpha\to\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\to\alpha$ phase
transitions, without retaining any $\epsilon$-Fe hcp phase nor without forming
any $\alpha^{\prime}$-Fe bcc variants. Interestingly, this feature differs
from the pure hydrostatic compression loading, for which a significant
residual volume fraction ($\sim 12\%$) of $\alpha^{\prime}$ 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$ m.s-1, which is roughly twice the particle
velocity of $850$ 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]_{\textrm{bcc}}$ direction, characterized by a pressure
$p_{\textrm{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., $\sim
2.1~{}$GPa [297] and $\sim 5.5~{}$GPa [234]. 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_{\textrm{PT}}=37.7~{}$GPa. The latter profile shows an intermediate
pressure plateau that characterizes the plastically deformed bcc region,
within which the forward $\alpha\to\epsilon$ phase transitions start roughly
at the onset pressure $p_{\,\textrm{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 [76], and in excellent
agreement with large-scale molecular dynamics simulations in single-crystal
iron as well, i.e., $18~{}$GPa along the same $[001]_{\textrm{bcc}}$ shock
direction [280]. 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 [106, 282, 302].
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 $\sim
40~{}$GPa has been measured by applying a comparable piston velocity of 800
m.s-1.

Figure (2.11c) shows the corresponding values for the von Mises stress, with
$\sigma_{\textrm{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_{\textrm{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_{\textrm{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.

Figure 2.12:  (a) Slice-averaged maps of the longitudinal stress component
$\sigma_{zz}$ in the Lagrangian adimensional position-time
$(L^{\ast},t^{\ast})$ 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^{\ast}=0.10$
and $t^{\ast}=0.27$, which result in the inhomogeneous propagation of the 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.

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^{\ast},t^{\ast})$ 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.

The leading E wave front, traveling at $5412$ m.s-1 ($5541$ 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$ m.s-1
using atomistics simulations in single-crystal iron without pre-existing
defects [141], which is consistent with the present calculations. Without
plasticity, the trailing PT front travels homogeneously in the sample at
$4655$ 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$ m.s-1, while the slower heterogeneous PT
front travels with intermittent regimes at $3002\pm 99$ 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$ 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 [156].

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

Figure 2.13:  (a) Three-dimensional time snapshots of shock-induced
microstructures, designated M1 at $t^{\ast}=0.21$, (b) M2 at $t^{\ast}=0.27$,
and (c) M3 at $t^{\ast}=0.44$, from the position-time diagram in Fig. (2.12b)
for simulation with plasticity. From top to the bottom, each panel captures