#  Simulations of multivariant Si I to Si II phase transformation in
polycrystalline silicon with finite-strain scale-free phase-field approach

Hamed Babaei  Raghunandan Pratoori  Valery I. Levitas  Department of Aerospace
Engineering, Iowa State University, Ames, IA 50011, USA  Departments of
Mechanical Engineering, Iowa State University, Ames, IA 50011, USA  Ames
Laboratory, Division of Materials Science and Engineering, Ames, IA 50011, USA

###### Abstract

Scale-free phase-field approach (PFA) at large strains and corresponding
finite element method (FEM) simulations for multivariant martensitic phase
transformation (PT) from cubic Si I to tetragonal Si II in a polycrystalline
aggregate are presented. Important features of the model are large and very
anisotropic transformation strain tensor
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\\{0.1753;0.1753;-0.447\\}$ and stress-
tensor dependent athermal dissipative threshold for PT, which produce
essential challenges for computations. 3D polycrystals with 55 and 910
stochastically oriented grains are subjected to uniaxial strain- and stress-
controlled loadings under periodic boundary conditions and zero averaged
lateral strains. Coupled evolution of discrete martensitic microstructure,
volume fractions of martensitic variants and Si II, stress and transformation
strain tensors, and texture are presented and analyzed. Macroscopic variables
effectively representing multivariant transformational behavior are
introduced. Macroscopic stress-strain and transformational behavior for 55 and
910 grains are close (less than 10% difference). This allows the determination
of macroscopic constitutive equations by treating aggregate with a small
number of grains. Large transformation strains and grain boundaries lead to
huge internal stresses of tens GPa, which affect microstructure evolution and
macroscopic behavior. In contrast to a single crystal, the local mechanical
instabilities due to PT and negative local tangent modulus are stabilized at
the macroscale by arresting/slowing the growth of Si II regions by the grain
boundaries and generating the internal back stresses. This leads to increasing
stress during PT. The developed methodology can be used for studying similar
PTs with large transformation strains and for further development by including
plastic strain and strain-induced PTs.

###### keywords:

Multivariant martensitic phase transformation , Silicon polycrystal , Scale-
free phase-field approach , Finite element simulations , Stress-dependent
effective threshold , Finite strain

##  1 Introduction

Silicon is the second most abundant material in the earth’s crust. The
semiconducting Si I phase (cubic diamond lattice, $Fd3m$ space group) is
extensively used in microelectronics, integrated circuits, photovoltaics, and
MEMS/ NEMS technologies. Single-crystal Si is also used in high-power lasers.
Polycrystalline Si is widely used in solar panels [1], thin transistors [2],
and very large-scale integration (VLSI) manufacturing. It also has low
toxicity and high stability. Due to high demands, the recent CHIPS and Science
Act will provide new funding to boost the research and manufacturing of
semiconductors in the US. Under the pressure of 10-16 GPa, semiconducting Si I
transforms to metallic phase Si II ($\beta$-tin structure, $I4_{1}/amd$ space
group). Si I is very strong and brittle, and hence its bulk hardness is 12 GPa
and is determined by Si I$\rightarrow$Si II PT rather then dislocational
plasticity [3, 4]. Stresses exceed 10 GPa for machining (turning, polishing,
scratching, etc.) of single and polycrystalline Si [5]; such loads cause
plastic flow, Si I$\rightarrow$Si II and some other PTs, e.g., amorphization.
High-pressure torsion of Si at 24 GPa is used to produce nanostructured
metastable phases [6, 7]. Machining of strong brittle semiconducting Si-I is
accompanied by microcrack propagation inside the bulk. PT from Si I to ductile
and weaker Si II is utilized to develop ductile machining regimes [8], which
reduces forces, energy, and damage. This also may eliminate the necessity of
using chemical additives during machining, which brings definite environmental
benefits by reducing pollution.

Since transformation strain tensor that describes the transformation of cubic
to tetragonal Si I $\rightarrow$ Si II PT has large and very anisotropic
principal components
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\\{0.1753;0.1753;-0.447\\}$ [9, 10,
11], it is clear from thermodynamics that deviatoric part of the stress tensor
should strongly affect this PT. Both pressure- and stress-induced PTs start at
pre-existing defects (different dislocation configurations, grain boundaries),
which represent stress concentrators. However, in many of the applications,
like turning, polishing, scratching, friction, high-pressure torsion and ball
milling, PTs occur during plastic deformation. According to classification
[12, 13, 14], such PTs are called plastic strain-induced PT under high
pressure, and they occur at defects permanently generated during plastic flow.
Strain-induced PTs require completely different thermodynamic, kinetic, and
experimental treatments than pressure- and stress-induced PTs.

There are numerous very strong effects of plastic deformations on PTs,
summarized in [15, 16, 17, 18, 12, 13, 14]; one of the most important is a
drastic reduction in PT pressure. Thus, plastic strain-induced PTs from
graphite to hexagonal and cubic diamonds were obtained at 0.4 and 0.7 GPa, 50
and 100 times lower than under hydrostatic loading, respectively, and well
below the phase equilibrium pressure of 2.45 GPa [18]. About an order of
magnitude reduction in PT pressure was reported for PT from rhombohedral to
cubic BN [19], hexagonal to wurtzitic BN [20], and from $\alpha$ to $\omega$
Zr [21, 22].

The effect of plastic straining on PTs in Si is also very strong but more
sophisticated. Thus, under compression and shear in rotational diamond anvil
cell [23], Si II was obtained at 5.2 GPa, but not directly from Si I, but via
Si III. However, in these experiments, only optical, pressure, and electric
resistivity measurements were utilized without in situ x-ray diffraction. In
our recent in situ x-ray diffraction experiments [24], the PTs pressure for
direct Si I$\rightarrow$Si II PT was reduced from 13.5 (hydrostatic loading)
to 2.5 GPa (plastic straining) for micron size Si particles and from 16.2 GPa
to 0.4 GPa for 100 nm Si nanoparticles, i.e., by a factor of 40.5 (and 26.3
below the phase equilibrium pressure of 10.5 GPa [25]).

To understand the effect of stress tensor and plastic strain on Si
I$\rightarrow$Si II PT, various techniques at multiple scales are used. With
first principle simulations, the lattice instability for Si I under two-
parametric loadings was studied in [26, 27, 28, 29]. Although it is not
specified, it is related to Si I$\rightarrow$Si II PT. The stress-strain
behavior and elastic instabilities leading to Si I$\rightarrow$Si II PT under
all six components of the stress tensor were determined in [30, 31]. While
this seems to be impossible due to a large number of combinations, the
following solution was found. First, the analytical expression for crystal
lattice instability criterion was formulated within the PFA [32, 10, 11]. Then
it was confirmed and quantified with first principle simulations in [30].

A similar approach was realized earlier with molecular dynamics simulations
[10, 11] under the action of three normal stresses. Molecular dynamics
simulations of PTs in Si during various loadings, nanoindentation, scratching,
and surface processing were performed in [33, 34, 35, 36, 4, 37].

Nanoscale PFA for PT Si I$\rightarrow$Si II in a single crystal with
corresponding simulations was developed in [38, 39, 40, 41, 42]. It was
calibrated by results of molecular dynamics simulations in [10, 11]. The
effect of a single dislocation on this PT under uniaxial compression is
modeled in [40]. The main general problem with nanoscale PFA is that the width
of the martensitic interface is 1 nm and one needs at least 3-5 finite
elements within the interface [43], making problem computationally expensive
for large samples. Hence this model can be used to treat nano-sized samples
only.

Here, we will consider scale-free PFA for modeling discrete martensitic
microstructure. It was developed for small strains in [44, 45] and updated and
applied to NiTi shape memory alloy in [46]. The only finite-strain
generalization of the scale-free model for finite strains is presented in
[47]. This model was applied for simulations of Si I$\rightarrow$Si II PT in a
single crystal.

The main differences between scale-free and nanoscale PFAs are:

  1. 1.

The energy term that includes gradients of the order parameters and determines
the phase interface widths and the interface energies is excluded. This makes
the model scale-free. Interface width is getting equal to a single finite
element, which is much more computationally economical than in the nanoscale
approach, where 3-5 elements are usually required to reproduce analytical
solution for an interface [43]. However, this leads to two problems. (a) The
solution is mesh dependent. However, detailed computational experiments in
[44, 45, 47] show that the solution becomes practically mesh-independent after
the mesh size is 80 times smaller than the sample size. (b) Since the volume
fraction of martensite varies from 0 to 1 within the one-element thick
interface, for large transformation strains, there are large strain gradients
within the interface, which leads often to divergence in the FEM solution.

  2. 2.

The interfaces between individual martensitic variants are not resolved. Each
of these variant has a width of $d\simeq 10\;\text{nm}$, i.e., thousands of
interfaces at the microscale, making the problem computationally impractical.
Here, martensite is considered a mixture of martensitic variants with
corresponding volume fractions.

  3. 3.

Since there is no need to reproduce atomic level energy landscape versus order
parameters, the linear mixture rule is applied for all material properties.
This is in contrast to higher-order polynomials in order parameters for
nanoscale PFA. Also, the athermal dissipative threshold for interface
propagation (interface friction) can be easily introduced in the scale-free
model, while this is a problem for the nanoscale model.

  4. 4.

The volume fraction of the martensite is the order parameter, i.e., it is
responsible for the material instability and transformation strain
localization at the interface between austenite and martensite. The volume
fractions of the individual martensitic variants are only internal variables
and do not produce any instabilities. That is how we eliminate interfaces
between martensitic variants.

All the above features allow us to economically model multivariant martensitic
PTs in a sample of arbitrary size.

The next natural step is to simulate PTs in Si and discrete martensitic
microstructure in a polycrystalline sample, which was not done yet. This is
also an important step in studying plastic strain-induced PTs. The main
hypothesis [12] is that they initiate at the tip of dislocation pileups
against grain boundaries, which causes strong stress concentration for all
stress components proportional to the number of dislocations in a pileup.
These stresses drastically reduce the required external pressure needed for
nucleation and further growth. Nanoscale PFA allows us to treat a bicrystal
and to qualitatively prove this hypothesis [48, 49, 50], however, for 2D and
small strain formulation. Our scale-free PFA [51, 52] has been applied to 2D
polycrystalline samples with several dozen grains and also for small strain
formulation.

In this work, we extend the modeling presented in [47] to Si I$\rightarrow$Si
II PT in 3D polycrystalline samples with up to 1000 grains. The main challenge
is to reach convergence of the computational procedure due to strong
nonlinearities, large transformation strain localization with large gradients
within a single-element diffuse interface, and potential elastic
instabilities. For a single crystal, just two loadings, uniaxial and
hydrostatic, have been considered in [47]. However, each grain is subjected to
different complex and heterogeneous loading in a polycrystal. While we used
the simplest quadratic in elastic Lagrangian strain expression (7) for the
elastic energy, it is shown in [53] that even simple uniaxial compression
leads to elastic instability. This instability can cause additional strain
localization and divergence of solution. Due to the variety of heterogeneous
complex loadings, different for different grains, chances for various elastic
instabilities and divergence are very high. That is why several computational
parameters have been varied to reach convergence of the solution. Periodic
boundary conditions along the lateral surfaces played important part in the
avoiding divergence. Another problem is to adequately quantitatively present
the evolution of the local and average volume fraction of martensitic
variants; a straightforward approach leads to contradictory results.

The paper is organized as follows. In Section 2 complete system of coupled PFA
and nonlinear mechanics equations is presented. Materials parameters for the
model are given in Section 3. Problem formulations, results of simulations,
and their analyses are presented in Section 4. Concluding remarks are
summarized in Section 5. Vectors and tensors are designated with boldface
symbols. We designate contractions of tensors
${\mbox{\boldsymbol{$A$}}}=\\{A_{ij}\\}$ and
${\mbox{\boldsymbol{$B$}}}=\\{B_{ji}\\}$ over one and two indices as
${\mbox{\boldsymbol{$A$}}}{\mbox{\boldsymbol{$\cdot$}}}{\mbox{\boldsymbol{$B$}}}=\\{A_{ij}\,B_{jk}\\}$
and
${\mbox{\boldsymbol{$A$}}}{\mbox{\boldsymbol{$:$}}}{\mbox{\boldsymbol{$B$}}}=A_{ij}\,B_{ji}$.
The transpose of $A$ is ${\mbox{\boldsymbol{$A$}}}^{T}$; $I$ is the unit
tensor; $\mbox{\boldsymbol{$\nabla$}}_{0}$ is the gradient operator in the
undeformed state.

##  2 Model description

Here, we expand the microscale model for multivariant martensitic PTs
developed in [47] for a polycrystalline elastic materials.

###  2.1 Kinematics:

Let us consider polycrystalline Si I aggregate in the undeformed stress-free
configuration $\Omega_{0}$. The orientation of each grain is characterized by
the rotation tensor ${\mbox{\boldsymbol{$R$}}_{g}}$ in the undeformed
configuration, which rotates local cubic crystallographic axes of the grain to
the global coordinate system. Tensors ${\mbox{\boldsymbol{$R$}}_{g}}$ do not
evolve during loading and PTs. The deformation gradient $F$ is
multiplicatively split into the elastic $\mbox{\boldsymbol{$F$}}_{e}$ and the
transformational $\mbox{\boldsymbol{$F$}}_{t}$ parts:

${\mbox{\boldsymbol{$F$}}}:=\mbox{\boldsymbol{$\nabla$}}_{0}{\mbox{\boldsymbol{$r$}}}={\mbox{\boldsymbol{$F$}}}_{e}{\mbox{\boldsymbol{$\cdot$}}}{\mbox{\boldsymbol{$F$}}}_{t};\qquad{\mbox{\boldsymbol{$F$}}}_{e}={\mbox{\boldsymbol{$R$}}}_{e}\cdot{\mbox{\boldsymbol{$U$}}}_{e};\qquad{\mbox{\boldsymbol{$F$}}}_{t}={\mbox{\boldsymbol{$F$}}}_{t}^{T}.$
(1)

Here, $r$ is the position vector in the current deformed configuration
$\Omega$; ${\mbox{\boldsymbol{$F$}}}_{t}$ is defined as $F$ after complete
local stress release (producing the intermediate configuration $\Omega_{t}$)
and is considered to be rotation-free; ${\mbox{\boldsymbol{$U$}}}_{e}$ is the
symmetric elastic right stretch tensor and ${\mbox{\boldsymbol{$R$}}}_{e}$ is
the orthogonal lattice rotation tensor during loading and PT, which determines
texture evolution. The transformation deformation gradient is defined using
the mixture rule for $m$ variants as

$\mbox{\boldsymbol{$F$}}_{t}=\mbox{\boldsymbol{$I$}}+\mbox{\boldsymbol{$\varepsilon$}}_{t}=\mbox{\boldsymbol{$I$}}+\sum_{i=1}^{m}{\mbox{\boldsymbol{$R$}}}_{g}\cdot\tilde{\mbox{\boldsymbol{$\varepsilon$}}}_{ti}\cdot{\mbox{\boldsymbol{$R$}}}_{g}^{T}c_{i},$
(2)

where $\mbox{\boldsymbol{$\varepsilon$}}_{t}$ is the transformation strain,
$\tilde{\mbox{\boldsymbol{$\varepsilon$}}}_{ti}$ is the transformation strain
of the $i^{th}$ martensitic variant in the local crystallographic basis of the
grain, and the $c_{i}$ is the volume fraction of the $i^{th}$ variant in terms
of volumes in the reference configuration.

The total, elastic, and transformational Lagrangian strains are defined as

$\mbox{\boldsymbol{$E$}}=\frac{1}{2}(\mbox{\boldsymbol{$F$}}^{T}{\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$F$}}-\mbox{\boldsymbol{$I$}});\quad\mbox{\boldsymbol{$E$}}_{e}=\frac{1}{2}(\mbox{\boldsymbol{$F$}}_{e}^{T}{\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$F$}}_{e}-\mbox{\boldsymbol{$I$}});\quad\mbox{\boldsymbol{$E$}}_{t}=\frac{1}{2}(\mbox{\boldsymbol{$F$}}_{t}^{T}{\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$F$}}_{t}-\mbox{\boldsymbol{$I$}}).\qquad$
(3)

Utilization of multiplicative decomposition eq. 1 results in the following
relationship:

$\mbox{\boldsymbol{$E$}}_{e}=\mbox{\boldsymbol{$F$}}_{t}^{-1}{\mbox{\boldsymbol{$\cdot$}}}(\mbox{\boldsymbol{$E$}}-\mbox{\boldsymbol{$E$}}_{t}){\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$F$}}_{t}^{-1}.$
(4)

We will also need the ratios of elemental volumes $dV$ and mass densities
$\rho$ in the different configurations, which are described by the Jacobian
determinants:

$J=\frac{dV}{dV_{0}}=\frac{\rho_{0}}{\rho}=\det{\mbox{\boldsymbol{$F$}}};\quad
J_{t}=\frac{dV_{t}}{dV_{0}}=\frac{\rho_{0}}{\rho_{t}}=\det{\mbox{\boldsymbol{$F$}}_{t}};\quad
J_{e}=\frac{dV}{dV_{t}}=\frac{\rho_{t}}{\rho}=\det{\mbox{\boldsymbol{$F$}}_{e}}.$
(5)

###  2.2 Helmholtz Free energy

In defining the Helmholtz free energy $\psi$ per unit undeformed volume in
$\Omega_{0}$ of the mixture of austenite and $m$ martensitic variants,
contributions from the elastic $\psi^{e}$, thermal $\psi_{i}^{\theta}$, and
interaction $\psi^{in}$ energy are given by

$\psi(\mbox{\boldsymbol{$F$}}_{e},c_{i},\theta)=J_{t}\psi^{e}(\mbox{\boldsymbol{$F$}}_{e},c_{i})+\psi^{\theta}(\theta,c_{i})+\psi^{in}(c_{i}).$
(6)

Here, $\psi^{e}$ is defined in $\Omega_{t}$, and the Jacobian $J_{t}$ maps it
into $\Omega_{0}$; $\psi_{i}^{\theta}$ includes the thermal driving force for
the PT and depends on the temperature $\theta$;
$\psi^{in}=\mathcal{A}cc_{0}\geq 0$ includes the interactions between
austenite and martensite, the energy of internal stresses, as well as the
austenite-martensite phase interface energy; interaction between martensitic
variants is neglected to avoid formation of variant-variant interfaces.
Positive $\mathcal{A}$ results in a negative tangent modulus of the
equilibrium stress-strain curve during PT, which results in local mechanical
instability and the formation of the localized transformation bands/regions of
the product phase. In such a way, a discrete martensitic structure is
reproduced, similar to the nanoscale PFA.

The elastic energy is expressed as

$\psi^{e}=\frac{1}{2}\mbox{\boldsymbol{$E$}}_{e}\mbox{\boldsymbol{$:$}}\tilde{\mbox{\boldsymbol{$C$}}}\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$E$}}_{e}=\frac{1}{2}\tilde{C}^{ijkl}{E}_{e}^{ij}{E}_{e}^{kl};\qquad\tilde{C}^{ijkl}=\sum_{p=0}^{m}\tilde{C}^{ijkl}_{p}c_{p},$
(7)

where the components of the forth-rank elastic moduli tensor
$\tilde{C}^{ijkl}_{p}$ in the global coordinate system are defined in terms of
components in the local for each grain crystallographic system
${C}^{rqnm}_{p}$ and grain rotations $R_{g}^{ml}$ as

$\tilde{C}^{ijkl}_{p}=R_{g}^{lm}R_{g}^{kn}R_{g}^{jq}R_{g}^{ir}{C}^{rqnm}_{p}.$
(8)

Since the thermal energy of all martensitic variants is the same,
$\psi_{i}^{\theta}=\psi_{j}^{\theta}=\psi_{M}^{\theta}$, the thermal energy of
the mixture is

$\psi^{\theta}=\sum_{i=0}^{m}c_{i}\psi_{i}^{\theta}(\theta)=c_{0}\psi_{A}^{\theta}(\theta)+c\psi_{M}^{\theta}(\theta);\qquad
c=\sum_{i=1}^{m}c_{i};\qquad c_{0}=1-c,$ (9)

where $\psi_{A}^{\theta}$ is the thermal energy of austenite, $c$ and $c_{0}$
are the volume fractions of the martensite and austenite.

###  2.3 Dissipation inequality

The Plank’s inequality for isothermal processes is

$D=\mbox{\boldsymbol{$P$}}^{T}\mbox{\boldsymbol{$:$}}\dot{\mbox{\boldsymbol{$F$}}}-\dot{\psi}\geq
0,$ (10)

where $D$ is the dissipation rate per unit undeformed volume; $P$ is the first
Piola-Kirchhoff stress. After traditional thermodynamic manipulations, $D$ can
be expressed as the product of the thermodynamic driving forces for
$A\rightarrow M_{i}$, $X_{i0}$, and $M_{j}\rightarrow M_{i}$, $X_{i0}$,
transformations and conjugate rates:

$\displaystyle
D=\sum_{i=1}^{m}X_{i0}\dot{c}_{i0}+\sum_{j=1}^{m-1}\sum_{i=j+1}^{m}X_{ij}\dot{c}_{ij}\geq
0;$ $\displaystyle
X_{i0}=W_{i0}-\frac{J_{t}}{2}\mbox{\boldsymbol{$E$}}_{e}\mbox{\boldsymbol{$:$}}(\bar{\mbox{\boldsymbol{$C$}}}_{i}-\bar{\mbox{\boldsymbol{$C$}}}_{0})\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$E$}}_{e}-\frac{J_{t}}{2}\left(\mbox{\boldsymbol{$E$}}_{e}\mbox{\boldsymbol{$:$}}\bar{\mbox{\boldsymbol{$C$}}}(c)\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$E$}}_{e}\right){\mbox{\boldsymbol{$F$}}}_{t}^{-1}\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$\varepsilon$}}_{ti}-\Delta\psi^{\theta}-\mathcal{A}(1-2c);$
$\displaystyle
X_{ij}=W_{ij}-\frac{J_{t}}{2}\mbox{\boldsymbol{$E$}}_{e}\mbox{\boldsymbol{$:$}}(\bar{\mbox{\boldsymbol{$C$}}}_{i}-\bar{\mbox{\boldsymbol{$C$}}}_{j})\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$E$}}_{e}-\frac{J_{t}}{2}\left(\mbox{\boldsymbol{$E$}}_{e}\mbox{\boldsymbol{$:$}}\bar{\mbox{\boldsymbol{$C$}}}(c)\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$E$}}_{e}\right){\mbox{\boldsymbol{$F$}}}_{t}^{-1}\mbox{\boldsymbol{$:$}}(\mbox{\boldsymbol{$\varepsilon$}}_{ti}-\mbox{\boldsymbol{$\varepsilon$}}_{tj});$
$\displaystyle
W_{i0}=\mbox{\boldsymbol{$P$}}^{T}\mbox{\boldsymbol{$\cdot$}}\mbox{\boldsymbol{$F$}}_{e}\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$\varepsilon$}}_{ti}=J{\mbox{\boldsymbol{$F$}}}^{T}_{e}{\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$\sigma$}}\mbox{\boldsymbol{$\cdot$}}{\mbox{\boldsymbol{$F$}}}_{e}^{t-1}\cdot{\mbox{\boldsymbol{$F$}}}_{t}^{-1}\mbox{\boldsymbol{$:$}}\mbox{\boldsymbol{$\varepsilon$}}_{ti};\quad
W_{ij}=\mbox{\boldsymbol{$P$}}^{T}\mbox{\boldsymbol{$\cdot$}}\mbox{\boldsymbol{$F$}}_{e}\mbox{\boldsymbol{$:$}}(\mbox{\boldsymbol{$\varepsilon$}}_{ti}-\mbox{\boldsymbol{$\varepsilon$}}_{tj})=J{\mbox{\boldsymbol{$F$}}}^{T}_{e}{\mbox{\boldsymbol{$\cdot$}}}\mbox{\boldsymbol{$\sigma$}}\mbox{\boldsymbol{$\cdot$}}{\mbox{\boldsymbol{$F$}}}_{e}^{t-1}\cdot{\mbox{\boldsymbol{$F$}}}_{t}^{-1}\mbox{\boldsymbol{$:$}}(\mbox{\boldsymbol{$\varepsilon$}}_{ti}-\mbox{\boldsymbol{$\varepsilon$}}_{tj}).$
(11)

Here, $\dot{c}_{i0}$ and $\dot{c}_{ij}$ are the rate of change of volume
fraction of variant $i$ due to transformation to the austenite and variant
$j$, respectively; $W_{i0}$ is the transformation work for austenite to
martensite PT, $W_{ij}$ is the transformation work for variant $j$ to variant
$i$ transformation, $\Delta\psi^{\theta}$ is the jump in thermal energy during
transformation, and
$\mbox{\boldsymbol{$\sigma$}}=J^{-1}\mbox{\boldsymbol{$P$}}\mbox{\boldsymbol{$\cdot$}}\mbox{\boldsymbol{$F$}}^{T}$
is the Cauchy (true) stress tensor.

###  2.4 Kinetic equations

The kinetic equations are formulated as follows for the $A\leftrightarrow
M_{i}$ PTs

$\begin{cases}\begin{split}\dot{c}_{i0}=\lambda_{i0}(X_{i0}-k_{i-0})\quad\text{if}\;&\\{X_{i0}-k_{i-0}(c_{i},\sigma_{i})>0\;\&\;c_{i}<1\;\&\;c_{0}>0\\}\;\qquad
A\rightarrow M_{i}\\\
\dot{c}_{i0}=\lambda_{i0}(X_{i0}+k_{i-0})\quad\text{if}\;&\\{X_{i0}+k_{i-0}(c_{i},\sigma_{i})<0\;\&\;c_{i}>0\;\&\;c_{0}<1\\}\qquad\;M_{i}\rightarrow
A\end{split}\\\ \dot{c}_{i0}=0\qquad\qquad\qquad\;\text{otherwise;}\qquad
i=1,2,...,m,\end{cases}$ (12)

and for $Mj\leftrightarrow M_{i}$ PTs

$\begin{cases}\begin{split}\dot{c}_{ij}=\lambda_{ij}X_{ij}\quad\text{if}\;&\\{X_{ij}>0\;\&\;c_{i}<1\;\&\;c_{j}>0\\}\qquad
j\rightarrow i\\\ \text{or}\;&\\{X_{ij}<0\;\&\;c_{i}>0\;\&\;c_{j}<1\\}\qquad
i\rightarrow j\end{split}\\\
\dot{c}_{ij}=0\qquad\qquad\qquad\qquad\qquad\quad\text{otherwise;}\qquad
i,j=1,2,...,m,\end{cases}$ (13)

where $k_{i-0}$ is the athermal threshold and $\lambda_{i0}$ and
$\lambda_{ij}$ are the kinetic coefficients. We neglect the athermal threshold
for $Mj\leftrightarrow M_{i}$ PTs. The non-strict inequalities for the volume
fraction of phases in Eqs. (12)-(13) imply that the PT from any phase does not
occur if the parent phase does not exist or if the product phase is complete.

###  2.5 Macroscopic parameters

Macroscopic Cauchy stress and the first Piola-Kirchhoff stress, the
deformation gradient, transformation strain, and volume fraction of Si II and
each martensitic variant, averaged over the sample, are defined as [54, 55,
56]

$\displaystyle\bar{\mbox{\boldsymbol{$\sigma$}}}=\frac{1}{V}\int_{V}\mbox{\boldsymbol{$\sigma$}}dV;\qquad\bar{\mbox{\boldsymbol{$P$}}}=\frac{1}{V_{0}}\int_{V_{0}}{\mbox{\boldsymbol{$P$}}}dV_{0};$
(14)
$\displaystyle\bar{\mbox{\boldsymbol{$F$}}}=\frac{1}{V_{0}}\int_{V_{0}}\mbox{\boldsymbol{$F$}}dV_{0};\qquad\bar{\mbox{\boldsymbol{$F$}}}_{t}\simeq\frac{1}{V_{0}}\int_{V_{0}}\mbox{\boldsymbol{$F$}}_{t}dV_{0};\qquad\bar{\mbox{\boldsymbol{$\varepsilon$}}}_{t}\simeq\frac{1}{V_{0}}\int_{V_{0}}{\mbox{\boldsymbol{$\varepsilon$}}}_{t}dV_{0};$
(15)
$\displaystyle\bar{c}=\frac{1}{V_{0}}\int_{V_{0}}cdV_{0};\qquad\bar{c}_{i}=\frac{1}{V_{0}}\int_{V_{0}}c_{i}dV_{0}.$
(16)

For the Cauchy stress and the first Piola-Kirchhoff stress, the deformation
gradient, and volume fractions $c$ and $c_{i}$, averaging equations are
strict; for the transformation deformation gradient and strain, they are
approximate because, in the unloaded stress-free state, they generally are not
compatible due to residual elastic strain. We use Eq. (15) because the exact
equation is quite bulky and is for $\dot{\mbox{\boldsymbol{$F$}}}_{t}$ instead
of ${\mbox{\boldsymbol{$F$}}}_{t}$. Eq. (16) for $c_{i}$, while formally
correct, is misleading (Section 4.1). Other macroscopic parameters are defined
via $\bar{\mbox{\boldsymbol{$F$}}}$ and $\bar{\mbox{\boldsymbol{$P$}}}$ by
equations similar to the corresponding local equations:

$\displaystyle\bar{\mbox{\boldsymbol{$\sigma$}}}=(\det\bar{\mbox{\boldsymbol{$F$}}})^{-1}\bar{\mbox{\boldsymbol{$P$}}}\mbox{\boldsymbol{$\cdot$}}\bar{\mbox{\boldsymbol{$F$}}}^{T};\qquad\bar{\mbox{\boldsymbol{$E$}}}=\frac{1}{2}(\bar{\mbox{\boldsymbol{$F$}}}^{T}{\mbox{\boldsymbol{$\cdot$}}}\bar{\mbox{\boldsymbol{$F$}}}-\mbox{\boldsymbol{$I$}}).$
(17)

Eq. (17) for the Cauchy stress is used instead of Eq. (14) because integration
over the fixed parallelepiped with regular cubic mesh is much simpler and
faster than the integration over the deformed volume.

##  3 Model parameters

The material parameters required for the implementation of the model, the same
as in [47], are presented in section 2 are listed in table 1. The
transformation strain tensors are taken from the MD simulations in [10, 11],
and in the local crystallographic axes for all three variants are given by

(18)

The elastic constants for both phases are collected from [57, 31]. The
constants $C_{0}^{ij}$ in table 1 denote the independent elastic constants of
the austenite, and ${C}_{1}^{ij}$ denote those of the first variant of the
martensite, both in the local crystallographic axes. The components of the
tensor of elastic moduli $C^{ijkl}$ can be calculated using

$\displaystyle
C^{ijkl}=\sum_{n=1}^{3}[\lambda^{n}\delta^{in}\delta^{jn}\delta^{kn}\delta^{ln}+\mu^{n}(\delta^{in}\delta^{jn}\delta^{kl}+\delta^{ij}\delta^{kn}\delta^{ln})$
$\displaystyle+\nu^{n}(\delta^{in}\delta^{jk}\delta^{ln}+\delta^{jn}\delta^{ik}\delta^{ln}+\delta^{in}\delta^{jl}\delta^{kn}+\delta^{jn}\delta^{il}\delta^{kn})],$
(19)

where $\lambda^{n}$, $\mu^{n}$ and $\nu^{n}$ for a cubic and tetragonal
crystal lattice are given by sections 3 and 3, respectively:

$\displaystyle\lambda^{1}=\lambda^{2}=\lambda^{3}=C^{11}-C^{12}-2C^{44},$
$\displaystyle 2\mu^{1}=2\mu^{2}=2\mu^{3}=C^{12},$ $\displaystyle
2\nu^{1}=2\nu^{2}=2\nu^{3}=C^{44}.$ (20)
$\displaystyle\lambda^{1}=\lambda^{2}=C^{11}-(C^{12}+2C^{66}),$
$\displaystyle\lambda^{3}=C^{33}+C^{12}+2C^{66}-2(C^{13}+2C^{44}),$
$\displaystyle 2\mu^{1}=2\mu^{2}=C^{12},\qquad\;\;2\mu^{3}=2C^{13}-C^{12},$
$\displaystyle 2\nu^{1}=2\nu^{2}=C^{66},\qquad\;\;2\nu^{3}=2C^{44}-C^{66}.$
(21)

Under hydrostatic conditions, the phase equilibrium pressure
$p^{eq}_{0}=10.5GPa$ [25] at $J_{t}=0.764$. The jump in the thermal energy
$\Delta\psi^{\theta}=-p^{eq}_{0}({J}_{t}-1)=2.47GPa$, where elastic strain and
change in elastic moduli are neglected. Localization of strain is required to
reproduce discrete microstructure. As noted in [47], to obtain strain
localization, the strain rate should be commensurate with the rate of
transformation. The kinetic coefficient $\lambda$ and interaction parameter
$A$ are selected to ensure that this condition is satisfied.

It is found with the first principle and molecular dynamics simulations [30,
10, 11] that the criteria for for Si I$\leftrightarrow$Si II PTs are linearly
dependent on the Cauchy stress components normal to the cubic faces. We assume
the same for the microscale experiments, where the role of defects is
effectively included. For cubic to tetragonal PTs, the PT criteria are given
by [40, 47]

$\displaystyle A\rightarrow M_{i}:\qquad
a^{d}(\sigma_{1}+\sigma_{2})+b^{d}\sigma_{3}>c^{d};$ $\displaystyle
M_{i}\rightarrow A:\qquad a^{r}(\sigma_{1}+\sigma_{2})+b^{r}\sigma_{3}<c^{r},$
(22)

where $a^{d},b^{d},c^{d},a^{r},b^{r}$ and $c^{r}$ are constants determined
empirically. To make the thermodynamic PT conditions consistent with
experimental conditions, the athermal threshold $k_{i-0}$ is considered to be
stress and volume fraction dependent, as described in detail in [47], and
calculated based on the relations

$k_{i-0}=J[a_{1}(c_{i})(\sigma_{1}+\sigma_{2})+a_{3}(c_{i})\sigma_{3}];\qquad
a_{k}(c_{i})=d_{k}+(r_{k}-d_{k})c_{i};\;\;k=1\;\rm{and}\;3,$ (23)

where $d_{k}$ and $r_{k}$ are the fitting parameters, given in table 1. They
are determined by substituting expressions for $k_{i-0}$ from eq. 23 in the
transformation criteria in eq. 12 with the driving forces from section 2.3 and
some empirical data.

Table 1:  Material parameters including kinetic coefficient $\lambda\,\rm{(Pa\cdot s)}^{-1}$, dimensionless constants in the expression for effective thresholds, as well as interaction coefficient $A$, jump in the thermal energy $\Delta\psi^{\theta}$, and elastic constants, all in GPa.  $\lambda$ | $A$ | $\Delta\psi^{\theta}$ | $d_{1}$ | $d_{3}$ | $r_{1}$ | $r_{3}$ |  |
---|---|---|---|---|---|---|---|---
0.02 | 2 | 2.47 | 0.082 | 0.111 | -0.90 | 0.338 |  |
$C_{0}^{11}$ | $C_{0}^{44}$ | $C_{0}^{12}$ | $C_{1}^{11}$ | $C_{1}^{33}$ | $C_{1}^{44}$ | $C_{1}^{66}$ | $C_{1}^{12}$ | $C_{1}^{13}$
167.5 | 80.1 | 65.0 | 174.76 | 136.68 | 60.24 | 42.22 | 102.0 | 68.0

Figure 1: Grain distribution generated from DREAM.3D for (a) 55 grains showing
the local orientations of individual grains (red=1, yellow=2, green=3) and (b)
910 grains. Each grain has a different orientation that is randomly assigned
to make sure the sample, on the whole, is texture-free. Files with complete
information about orientation of each grain for both samples are given in
supplementary material.

##  4 Study of multivariant microstructure evolution

A finite element implementation of the scale-free model is developed in the
open-source FEM code deal.II [58] using 8-noded 3D cubic linear elements with
first-order interpolation and full integration. For such elements,
calculations of volume averaged of any parameter $a$ in the undeformed
configuration is the sum of values $a$ in each quadrature point divided by the
number of quadrature points. Microstructure evolution is studied for
polycrystalline samples containing two different numbers of grains, 55 and
910, as shown in fig. 1. The total number of the finite elements and
integration points were $\simeq$2.1 million and $\simeq$16.7 million,
respectively, for both cases. If we consider position vector $r$ and volume
fractions $c_{1}$, $c_{2}$, and $c_{3}$ as the primary independent variables,
then the total number of degrees of freedom is six times the number of the
integration points.

A 3D unit cube sample is constructed in DREAM.3D [59] with the required number
of grains. Each grain in the sample is assigned an orientation randomly so
that the overall sample remains texture free. Two types of compression are
applied, strain-controlled and stress-controlled. For both cases, periodic
boundary conditions are applied in directions 1 and 2 with averaged zero
averaged strains in these directions. Application of the periodic boundary
conditions significantly simplified elimination of divergence of the
computational procedure in comparison with other conditions. Such conditions
are realized when Si is clamped in the 1-2 plane; a thin Si layer is attached
to the rigid substrate, or within a shock wave. For the strain-controlled
case, periodic boundary conditions are applied in the direction 3 as well, and
a sample is subjected to a uniaxial compressive averaged strain in direction
3. For the stress-controlled loading, on the external planes orthogonal to
axis 3, the sample is subjected to homogeneous compressive normal Cauchy
stress in the deformed configuration with zero shear stresses, i.e., like
under the action of the liquid. Stress-controlled loading ends with constant
”pressure in liquid ” of 11 GPa. In total, four simulations are run with two
different numbers of grains and two different loading conditions.

###  4.1 Strain-controlled loading

For strain-controlled loading, the initial strain rate used, in the elastic
regime, is $1\times 10^{-2}s^{-1}$ and reduced to $5\times 10^{-5}s^{-1}$ once
the material starts to transform. While we can choose any strain rate, such
low strain rates are chosen because they are not achievable in atomistic
simulations. Figures 2 and 3 show the evolution of martensite at three
intermediate stages (25%, 50%, and 75% of the simulation) and the final stage
of the simulation. The first row shows the total volume fraction of martensite
$c$ followed by those for individual martensitic variants. The black lines in
figs. 2 and 3 outline the grain boundaries for all the grains. The videos
demonstrating the evolution of the volume fractions are provided in the
supplementary material. Volume fractions of each martensitic variant
$\bar{c}_{i}$ and martensite $\bar{c}$ averaged over the sample based on Eq.
16 vs. strain are shown in Fig. 6. Nucleation starts mostly at the grain
boundaries and triple junctions, which represent stress concentrators.
Nucleation starts with the dominant variant, which produces maximum
transformation work; two other variants appear in the same regions to
accommodate the transformation strain and reduce internal stresses. A similar
proportion between variants remains during further growth. In grains with
different orientations, growth occurs either along the grain boundaries or
inside the grain, which is arrested either at grain boundaries or other
martensitic units. In each unit full transformation to Si II occurs quickly,
thus forming discrete Si II microstructure, as desired. With increasing
strain, both Si II broadening of existing Si II regions and the appearance of
new nuclei occur. At the end of the loading, PT is completed almost
everywhere, with small residual Si I pockets. They are caused by internal
stresses due to large transformation strains. For both numbers of grains,
$\bar{c}$ saturates at 96%.

It can be clearly seen from figs. 2, 3 and 6 that the second and the third
variants are dominating over the first one, which looks contradictory. To
better understand the reasons for these results, let, for simplicity and
illustration, consider 5 grains: grain 1 oriented with $[100]$ direction along
axis 3 (variant 1); grains 2 and 3, rotated by $\pm 90^{o}$ about axis 1
(variant 1), and grains 4 and 5, rotated by $\pm 90^{o}$ about axis 2 (variant
3). Due to symmetry, this is still the same single crystal. If we treat it
like a single crystal and transform it homogeneously till completion, we
obtain $c_{1}=1$, $c_{2}=c_{3}=0$, and
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\mbox{\boldsymbol{$\varepsilon$}}_{t1}$.
However, if we treat it as a polycrystal, in grain 1, we will have the same
$c_{1}=1$, $c_{2}=c_{3}=0$ and
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\mbox{\boldsymbol{$\varepsilon$}}_{t1}$.
In the local coordinate system of grains 2 and 3, this transformation strain
looks like
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\mbox{\boldsymbol{$\varepsilon$}}_{t2}$,
i.e., $c_{2}=1$, $c_{1}=c_{3}=0$. Similarly, in the local coordinate system of
grains 4 and 5, the transformation strain is
$\mbox{\boldsymbol{$\varepsilon$}}_{t3}$, i.e., $c_{3}=1$, $c_{1}=c_{2}=0$.
After averaging these volume fractions over the entire sample, we obtain
$c_{1}=1/5$, $c_{2}=c_{3}=2/5$, what we approximately observe in Fig. 6. The
main point is that ${c}_{i}$ for each grain does not have a lot of sense
unless the orientation of the grain is shown. That is why we give complete
information about orientation of each grain in supplementary material and show
orientations in fig. 1a. The average over the sample $\bar{c}_{i}$ does not
have any physical sense because the orientation of grains is not taken into
account; they cannot be used in the macroscopic theories. The average over the
sample $\bar{c}$ has clear physical sense.

We suggest the following ways to present a multivariant structure in the
polycrystal. For local presentation, one can show a triad of local
crystallographic axes in each grain along with fields $c_{i}$ (fig. 1a). One
can present fields of six components of
$\mbox{\boldsymbol{$\varepsilon$}}_{ti}c_{i}$ for each variant $M_{i}$ in the
global coordinate system, which is 18 fields. For the above example with 5
grains, in each grain, we will have components of
$\mbox{\boldsymbol{$\varepsilon$}}_{t3}$ in the global coordinate system,
which is consistent with the treatment of the aggregate as a single crystal.
More compact is to present six fields of six components of the total
transformation strain
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\sum_{i=1}^{m}\mbox{\boldsymbol{$\varepsilon$}}_{ti}c_{i}$
in the global coordinate system (Figs. 4 and 5), which takes into account in
more averaged sense both an orientation of grain and fields $c_{i}$. For the
averaged description, one can use a plot of six components of
$\bar{\mbox{\boldsymbol{$\varepsilon$}}}_{t}$.

As it follows from Figs. 4 and 5, normal components of
$\mbox{\boldsymbol{$\varepsilon$}}_{t}$ vary between extremes -0.447 and
0.1753 corresponding to the full local transformation in single variant
oriented along the global coordinate axes, despite zero averaged total lateral
strains. Shear components vary between extremes $~{}\pm 0.31$, also despite
zero averaged. Colors corresponding to zero value of each component of
transformation strain is clear from colors of large Si I regions at 25%
transformation progress corresponding $c=0$ in Figs. 2 and 3. For both grain
sizes one sees strong heterogeneity of transformation strains from grain to
grain and within grains. Note that fields at the surface do not completely
represent fields in the entire volume, that is why some results may look
counterintuitive. For example, while averaged $\varepsilon_{t}^{33}$ is larger
than $\varepsilon_{t}^{11}$ and $\varepsilon_{t}^{22}$, this is not evident
from the surface fields.

Despite the fact that for the smaller number of grains size of Si II units is
larger, difference in $\bar{c}(\bar{E}_{33})$ (and even
$\bar{c}_{i}(\bar{E}_{33})$) for different number of grains is small (fig. 6).
Initiation of PT occurs at the same stress $\sigma_{zz}=-10.33$ GPa for both
grain sizes, and stress-strain curves in fig. 9 also differ insignificantly.
That means that the current model does not describe experimentally observed
effect of the grain size on the Si I to Si II PT pressure or stress [24, 60,
61, 62]. The reason is that the current model does not include dislocations as
local stress concentrators in bulk and at grain boundaries. This will be the
next step in developing the current model. We can use the same approach for
introducing discrete dislocations via the solution of the contact problem, as
it was done in [51, 52] for small strains and 2D formulations. Of course, it
is much more challenging to do this for 3D and large strains.

Figures 7 and 8 shows the evolution of all components of Cauchy stresses
throughout the simulation. It can be noticed from figs. 7 and 8 that the grain
boundaries and triple junctions have the highest stresses of both signs, which
cause nucleation of the dominating variant accompanied by two other variants
and growth, often along the grain boundaries. Next, a large stress
concentration of both signs appears at the phase interface, causing further
nucleation in bulk (so-called autocatalytic effect). The peak stresses are
huge for both normal and even shear stresses. Thus, for small grains,
$\sigma_{22}$ varies from $-45$ to $35$ GPa and $\sigma_{33}$ varies from
$-40$ to $10$ GPa. Similar, shear stress $\sigma_{12}$ varies from $-20$ to
$20$ GPa and $\sigma_{23}$ varies from $-15$ to $15$ GPa. For large grains,
the magnitude of extremes in stresses is smaller by $5$ to $10$ GPa. Despite
this difference for different grain sizes, the macroscopic PT initiation
stress $\sigma_{zz}=-10.33$ GPa is the same, and the entire
$\sigma_{zz}-\varepsilon_{zz}$ curve do not differ significantly. The
initiation stress is determined by the transformation work, which depends on
all components of the stress tensor. While stresses are different, the
transformation work may be approximately the same for both grain sizes. This
is similar to results in [40], where a strong stress concentrator due to a
dislocation in Si produced a relatively small contribution to the
transformation work.

Figure 2: Evolution of volume fractions of phases for strain-controlled
loading of a sample with 55 grains. The figure shows snapshots at different
stages of completion of simulation (25%, 50%, 75%, and 100%) in different
columns for Si II, $c$, and each martensitic variant, $c_{i}$.

Figure 3: Evolution of volume fractions for strain-controlled loading of a
sample with 910 grains.

Figure 4: Evolution of components of the transformation strain tensor
$\mbox{\boldsymbol{$\varepsilon$}}_{t}$ for strain-controlled loading of a
sample with 55 grains.

Figure 5: Evolution of the components of the transformation strain tensor
$\mbox{\boldsymbol{$\varepsilon$}}_{t}$ for strain-controlled loading of a
sample with 910 grains.

Figure 6: Volume fraction of each martensitic variant ($\bar{c}_{i}$) and the
total martensite ($\bar{c}$) averaged over the sample based on Eq. (16) vs.
strain for strain-controlled loading. The lack of physical sense for
$\bar{c}_{i}$ is described in the text.

Figure 7: Evolution of components of the Cauchy stress tensor for strain-
controlled loading of a sample with 55 grains.

Figure 8: Evolution of the components of the Cauchy stress tensor for strain-
controlled loading of a sample with 910 grains.

The stress-strain plots for strain-controlled loading are given in fig. 9.
After PT starts at $\bar{E}_{33}=-0.06$ and $\bar{\sigma}_{33}=-10.33$ GPa for
both small and large grains, the stress-strain plots continue along the
elastic curve of the austenite with small nonlinearities. This behavior is
because, initially, the transformation rate is very low, as can be noticed
from fig. 6. Si II nuclei are localized near stress concentrators without
essential growth. Only at $\bar{E}_{33}=-0.07$ does intense growth start,
which leads to a strong reduction in tangent modulus. Note that for a single
crystal, the tangent modulus is getting negative at the onset of the PT,
causing macroscopic instability [47]. In contrast, for polycrystals, the local
mechanical instabilities due to PT and negative local tangent modulus are
stabilized at the macroscale by arresting/slowing the growth of Si II regions
by the grain boundaries and generating the internal back stresses. This is
reflected by the positive tangent moduli in the stress-strain plots in fig. 9.
While intuitively, the more grain boundaries we have, the higher the tangent
moduli should be, in fact, the response for 910 grains is slightly softer than
that of 55 grains fig. 9. The reasons are: (a) more triple junctions and
nucleation sites for smaller grains leading to Si II regions; (b) smaller
misorientation between neighboring grains leading to easier transfer of PT
growth from grain to grain, and (c) a larger number of surrounding grains
giving more chances to find proper orientation for new nucleation caused by
internal stresses in the transforming grains. Although there is a noticeable
difference in the response for both cases, the maximum difference is $<1$ GPa
or $<10\%$. Such a small difference implies that it is unnecessary to treat a
sample with such a large number of grains to estimate the macroscopic behavior
of a polycrystal. Note that the Lagrangian elastic strain at the end of
simulation, at $\bar{E}_{33}=-0.22$, is $~{}-0.1$, i.e., comparable to the
transformation strain.

Figure 9: Averaged Cauchy stress - Lagrangian strain plot for strain-
controlled loading for 55 grains and 910 grains. The dot marker represents the
onset of phase transformation.

Figures 10 and 11 shows the pole figures for 55 and 910 grains, respectively.
Figure 10(a) and Figure 11(a) for initial austenite demonstrate quite an even
spread in all the directions because the random texture of the stress-free
initial configuration was chosen. Figure 10(a) shows few low-density regions
because of the lower number of grains. Comparing the initial (Figure 10(a))
and final (Figure 10(b)) austenite for 55 grains, we can clearly notice
depletion of the austenite grains due to their transformation to martensite.
Similar depletion is not noticed in Figure 11(b) for 910 grains because many
grains do not completely transform to martensite. The increase in density in
Figure 10(b) or Figure 11(b) is because of the rotation of part of or entire
austenitic grains, which was observed for single crystal as well [47].

Figure 10(c)-(e) and Figure 11(c)-(e) present the pole figures for the three
martensitic variants at the last stage of simulations. As expected, there is a
90∘ rotation relation between the first and the second or the third variant
about the $(001)$ axis. The trends in the volume fractions of the individual
martensitic variant in Figures 10 and 11 match that of Figure 6. The first
variant is the lowest in both cases, as only the grains that have their
$c$-axis aligned with the loading direction can transform to the first
variant. As for the second or the third variants, they primarily need to be
oriented such that their $a$ or $b$-axis needs to be aligned with the loading
direction. But because of the boundary conditions used, there are lateral
stresses which can contribute to giving a resultant load in the preferred axes
assisting in the transformation. This phenomenon is not possible for the first
variant as the lateral stresses are only a fraction of the applied load and
cannot be the primary contributor to fulfilling the phase transformation
criteria.

Figure 10: Pole figures of (a) initial austenite, (b) residual austenite, (c)
final martensitic variant $M_{1}$, (d) $M_{2}$, and (e) $M_{3}$ for strain-
controlled loading of 55 grains. The minimum and maximum intensities for each
direction are shown in the legend corresponding to each figure.

Figure 11: Pole figures of (a) initial austenite, (b) residual austenite, (c)
final martensitic variant $M_{1}$, (d) $M_{2}$, and (e) $M_{3}$ for strain-
controlled loading of 55 grains. The minimum and maximum intensities for each
direction are shown in the legend corresponding to each figure.

###  4.2 Stress-controlled loading

In stress-controlled, the stress is applied at a rate of -4 MPa/s till it
reaches -11 GPa and held constant thereafter. Figures 14 and 15 and figs. 18
and 19 show the evolution of the volume fraction of Si II and the individual
martensitic variants as well as all stress components at different simulation
stages for 55 and 910 grins. The videos showing the evolution of the volume
fractions are provided in the supplementary material. Volume fraction of each
martensitic variant $\bar{c}_{i}$ averaged over the sample based on Eq. 16 and
$\bar{c}$ vs. strain is shown in Fig. 13. There is no significant difference,
in comparison with the strain-controlled case, in nucleation at grain
boundaries and triple junctions and character of growth of martensitic units,
stress concentrations, and that $\bar{c}_{2}\simeq\bar{c}_{3}\simeq
2\bar{c}_{1}$. The same discussion on the luck of the physical sense in
$\bar{c}_{i}$ is valid; $\bar{c}$ has a physical meaning only.

Figure 12: True stress-strain plot for stress-controlled loading for 55 and
910 grains. The dot marker represents the onset of phase transformation.

Figure 13: Volume fraction of each martensitic variant $\bar{c}_{i}$ and total
martensite $\bar{c}$) averaged over the sample based on Eq. 16 vs. strain for
stress-controlled loading. The lack of physical sense for $\bar{c}_{i}$ is
described in the text.

Figure 14: Evolution of volume fractions of phases for stress-controlled
loading of a sample with 55 grains.

Figure 15: Evolution of volume fractions of phases for strain-controlled
loading of a sample with 910 grains.

Figure 16: Evolution of the components of the transformation strain tensor for
stress-controlled loading of a sample with 55 grains.

Figure 17: Evolution of the components of the transformation strain tensor for
strain-controlled loading of a sample with 910 grains.

Figure 18: Evolution of the components of the Cauchy stress tensor for stress-
controlled loading of a sample with 55 grains.

Figure 19: Evolution of the components of the Cauchy stress tensor for strain-
controlled loading of a sample with 910 grains.

Stress-controlled compression for both numbers of grains is up to 42%
transformation to Si II, after which process diverges. The stress-strain plots
for strain-controlled loading are given in Figure 12. The PT starts
immediately at -11 GPa and continues to the end of simulations at
$E_{33}=-0.116$. These values are larger than those for strain-controlled
loading: at -11 GPa, we have $E_{33}=-0.063$ and $c=2.4\times 10^{-3}$ only;
however, at $E_{33}=-0.116$ one gets $c=0.362$ for strain-controlled loading.
The possibility of larger strains and transformation progress at -11 GPa is
related to less constraint deformation at the horizontal external surfaces.
Periodic conditions for displacements along the axis 3 for strain-controlled
loading lead to more homogenous PT near both horizontal surfaces and small
deviations from the flat surfaces. In contrast, for stress-controlled loading,
lack of periodic conditions along the axis 3 leads to much pronounced PT near
the upper horizontal surface, which spreads in the upper part of the sample.
This leads to the loss of stability of Si II nuclei near stress concentrators,
their fast growth, interaction, coalescence, and more pronounced auto-
catalytic effect. Thus, specific boundary conditions are very influential,
which is necessary to take into account in the problem formulation. The
difference in volume fractions of $M_{i}$ for 55 and 910 grains in fig. 6 and
$\bar{c}$ is much smaller than for the strain-controlled loading in Fig. 6,
again due to less restrictive boundary conditions.

The peak stresses are large but smaller than for strain-controlled loading
(figs. 18 and 19). Thus, for small grains, $\sigma_{22}$ varies from $-33$ to
$30$ GPa and $\sigma_{33}$ varies from $-25$ to $5$ GPa. Similar, shear stress
$\sigma_{12}$ and $\sigma_{23}$ vary from $-15$ to $15$ GPa. For large grains,
the magnitude of extremes in stresses are smaller by $5$ to $10$ GPa, like for
strain-controlled loading. Lower peak stresses are partially caused by a
smaller volume fraction of Si II.

Figures 20 and 21 show the pole figures for 55 and 910 grains, respectively,
for the stress-controlled loading. Figure 20(a) and Figure 21(a) show initial
austenite texture for both cases, which are the same as the ones chosen for
strain-controlled loading. For both 55 and 910 grains, there is no significant
depletion of the austenite, unlike for the strain-controlled loading, because
the $\bar{c}=0.42$ only. This is reflected in the small difference between
fig. 20(a) and fig. 20(b), and in sparsely distributed fig. 20(c)-(e).
Contrary to this, fig. 21(c)-(e) show the pole figures with more uniform
distributions because grains in many different orientations start transforming
to martensite like in the case of fig. 11.

Figure 20: Pole figures of (a) initial austenite, (b) residual austenite, (c)
final martensitic variant $M_{1}$, (d) $M_{2}$, and (e) $M_{3}$ for stress-
controlled loading of 55 grains.

Figure 21: Pole figures of (a) initial austenite, (b) residual austenite, (c)
final martensitic variant $M_{1}$, (d) $M_{2}$, and (e) $M_{3}$ for stress-
controlled loading of 910 grains.

##  5 Concluding remarks

In the paper, the first scale-free PFA modeling of the multivariant
martensitic PT from cubic Si I to tetragonal Si II in a polycrystalline
aggregate with up to 1000 grains is presented. All computational challenges
related to large and very anisotropic transformation strain tensors, the
stress-tensor dependent athermal dissipative threshold for the PT, and
potential elastic instabilities due to a variety of complex loadings in each
grain, are overcome. The importance of the simulations should also be stressed
by the fact that since Si II does not exist under ambient conditions, its
microstructure cannot be studied using traditional post-mortem methods (SEM,
TEM, etc.). For a single crystal, positions of Si I-Si II interfaces can be
determined using in-situ high-pressure Laue diffraction, but still, in
combination with molecular dynamics [37]. For polycrystals, this is currently
impossible.

Coupled evolution of discrete martensitic microstructure, volume fractions of
martensitic variants, and Si II, stress and transformation strain tensors, and
texture are presented and analyzed. It is demonstrated that the volume
fraction of each martensitic variant ${c}_{i}$ in each grain does not have a
lot of sense unless the orientation of the grain is explicitly shown. The
average over the sample $\bar{c}_{i}$ does not have any physical sense because
the orientation of grains is not taken into account; they are misleading and
cannot be used in the macroscopic theories. Macroscopic variables effectively
representing multivariant transformational behavior are introduced. One can
present fields of six components of
$\mbox{\boldsymbol{$\varepsilon$}}_{ti}c_{i}$ for each variant $M_{i}$ in the
global coordinate system. More compact is to present six components of the
total transformation strain
$\mbox{\boldsymbol{$\varepsilon$}}_{t}=\sum_{i=1}^{m}\mbox{\boldsymbol{$\varepsilon$}}_{ti}c_{i}$
in the global coordinate system. For the averaged description, one can utilize
six components of $\bar{\mbox{\boldsymbol{$\varepsilon$}}}_{t}$.

For strain-controlled uniaxial compression with periodic conditions in all
directions, almost complete (96%) PT was reached with small pockets of
residual Si I. For stress-controlled uniaxial compression without periodic
conditions in the loading direction, 42% of completion of was achieved at -11
GPa, much larger than 0.24% for the same axial stress for strain-controlled
loading, but relatively close volume fraction of Si II was reached for the
same strain. Lack of periodic conditions in the loading direction results in
less constraint deformation at the horizontal external surfaces and more
localized transformation near one of them, which leads to loss of stability of
Si II nuclei near stress concentrators, their fast growth, interaction, and
coalescence. Thus, tiny detail in the boundary conditions is very influential,
which is necessary to take into account in the problem formulation.

In contrast to a single crystal, the local mechanical instabilities due to PT
and negative local tangent modulus are stabilized at the macroscale by
arresting/slowing the growth of Si II regions by the grain boundaries and
generating the internal back stresses. This leads to increasing magnitude of
stress during PT.

Large transformation strains and grain boundaries lead to huge internal
stresses, which affect the microstructure evolution and macroscopic behavior.
The peak stresses reach 45 GPa in compression, 35 GPa in tension, and 20 GPa
in shear for strain-controlled loading of 910 grains; for 55 grains, the
magnitude of extremes in stresses are smaller by $5$ to $10$ GPa. For stress-
controlled loading and small grains, the peak stresses reach 35 GPa in
compression, 30 GPa in tension, and 15 GPa in shear; for large grains, they
are smaller by 5-10 GPa, like for strain-controlled loading. Lower peak
stresses for stress-controlled loading are partially caused by a smaller
volume fraction of Si II. Despite these differences, the macroscopic (overall)
stress-strain and transformational behavior for 55 and 910 grains are quite
close and differ less than by 10% for strain-controlled loading and even less
for stress-controlled one. On the good side, this allows the determination of
the macroscopic constitutive equations by treating aggregate with a small
number of grains. On the bad side, that means that the current model does not
describe experimentally observed effect of the grain size on the Si I to Si II
PT pressure or stress [24, 60, 61, 62]. The reason is that the current model
does not include dislocations as local stress concentrators in bulk and at
grain boundaries. This will be the next step in developing the current model.
We can use the same approach for introducing discrete dislocations via the
solution of the contact problem, as it was done in [51, 52] for small strains
and 2D formulations. Of course, it is much more challenging to do this for 3D
and large strains. A more realistic model of the grain boundaries is required
as well. The developed methodology can be used for studying various PTs with
large transformation strains (e.g., hexagonal and rhombohedral graphite to
hexagonal and cubic diamond, similar PTs from graphite-like BN to superhard
diamond-like BN, PTs in semiconducting Ge and GaSb, etc.) and for further
development for plastic strain-induced PTs.

Acknowledgement

The support of NSF (CMMI-1943710) and Iowa State University (Vance Coffman
Faculty Chair Professorship) is gratefully acknowledged. The simulations were
performed at Extreme Science and Engineering Discovery Environment (XSEDE),
allocation TG-MSS170015.

## References

  * [1] G. A. Heath, T. J. Silverman, M. Kempe, M. Deceglie, D. Ravikumar, T. Remo, H. Cui, P. Sinha, C. Libby, S. Shaw, et al., Research and development priorities for silicon photovoltaic module recycling to support a circular economy, Nature Energy 5 (7) (2020) 502–510.
  * [2] J. Chelikowsky, Introduction: Silicon in all its forms, in: P. Siffert, E. F. Krimmel (Eds.), Silicon: Evolution and Future of a Technology, Springer Berlin Heidelberg, Berlin, Heidelberg, 2004, pp. 1–22.
  * [3] V. Domnich, D. Ge, Y. Gogotsi, Indentation-induced phase transformations in semiconductors, in: V. Domnich, Y. Gogotsi (Eds.), High Pressure Surface Science and Engineering; Section 5.1, Institute of Physics, Bristol, 2004, pp. 381–442.
  * [4] M. S. Kiran, B. Haberl, J. E. Bradby, J. S. Williams, Chapter Five - Nanoindentation of Silicon and Germanium, in: L. Romano, V. Privitera, C. Jagadish (Eds.), Defects in Semiconductors, Vol. 91 of Semiconductors and Semimetals, Elsevier, 2015, pp. 165–203.  doi:https://doi.org/10.1016/bs.semsem.2014.12.002.
  * [5] S. Goel, X. Luo, A. Agrawal, R. L. Reuben, Diamond machining of silicon: a review of advances in molecular dynamics simulation, International Journal of Machine Tools and Manufacture 88 (2015) 131–164.
  * [6] Y. Ikoma, K. Hayano, K. Edalati, K. Saito, Q. Guo, Z. Horita, Phase transformation and nanograin refinement of silicon by processing through high-pressure torsion, Applied Physics Letters 101 (12) (2012) 121908.
  * [7] Y. Ikoma, K. Hayano, K. Edalati, K. Saito, Q. Guo, Z. Horita, T. Aoki, D. J. Smith, Fabrication of nanograined silicon by high-pressure torsion, Journal of Materials Science 49 (2014) 6565–6569.
  * [8] J. A. Patten, H. Cherukuri, J. Yan, Ductile-regime machining of semiconductors and ceramics, in: V. Domnich, Y. Gogotsi (Eds.), High-Pressure Surface Science and Engineering, Institute of Physics, Bristol, 2004, pp. 543–632.
  * [9] Z. Malyushitskaya, Mechanisms responsible for the strain-induced formation of metastable high-pressure Si, Ge, and GaSb phases with distorted tetrahedral coordination, Inorganic Materials 35 (5) (1999) 425–430.
  * [10] V. I. Levitas, H. Chen, L. Xiong, Lattice instability during phase transformations under multiaxial stress: Modified transformation work criterion, Phys. Rev. B 96 (2017) 054118.  doi:10.1103/PhysRevB.96.054118.
  * [11] V. I. Levitas, H. Chen, L. Xiong, Triaxial-Stress-Induced Homogeneous Hysteresis-Free First-Order Phase Transformations with Stable Intermediate Phases, Phys. Rev. Lett. 118 (2017) 025701.  doi:10.1103/PhysRevLett.118.025701.
  * [12] V. I. Levitas, High-pressure mechanochemistry: Conceptual multiscale theory and interpretation of experiments, Phys. Rev. B 70 (2004) 184118.  doi:10.1103/PhysRevB.70.184118.
  * [13] V. I. Levitas, High pressure phase transformations revisited, Journal of Physics: Condensed Matter 30 (16) (2018) 163001.  doi:10.1088/1361-648X/aab4b0.
  * [14] V. I. Levitas, High-pressure phase transformations under severe plastic deformation by torsion in rotational anvils, Materials Transactions 60 (7) (2019) 1294–1301.
  * [15] V. D. Blank, E. I. Estrin, Phase transitions in solids under high pressure, CRC Press, 2013.
  * [16] P. W. Bridgman, Effects of High Shearing Stress Combined with High Hydrostatic Pressure, Phys. Rev. 48 (1935) 825–847.  doi:10.1103/PhysRev.48.825.
  * [17] K. Edalati, Z. Horita, A review on high-pressure torsion (HPT) from 1935 to 1988, Materials Science and Engineering: A 652 (2016) 325–352.  doi:https://doi.org/10.1016/j.msea.2015.11.074.
  * [18] Y. Gao, Y. Ma, Q. An, V. Levitas, Y. Zhang, B. Feng, J. Chaudhuri, W. A. Goddard, Shear driven formation of nano-diamonds at sub-gigapascals and 300 K, Carbon 146 (2019) 364–368.  doi:https://doi.org/10.1016/j.carbon.2019.02.012.
  * [19] V. I. Levitas, L. K. Shvedov, Low-pressure phase transformation from rhombohedral to cubic BN: Experiment and theory, Phys. Rev. B 65 (2002) 104109\.  doi:10.1103/PhysRevB.65.104109.
  * [20] C. Ji, V. I. Levitas, H. Zhu, J. Chaudhuri, A. Marathe, Y. Ma, Shear-induced phase transition of nanocrystalline hexagonal boron nitride to wurtzitic structure at room temperature and lower pressure, Proceedings of the National Academy of Sciences 109 (47) (2012) 19108–19112.
  * [21] K. Pandey, V. I. Levitas, In situ quantitative study of plastic strain-induced phase transformations under high pressure: Example for ultra-pure Zr, Acta Materialia 196 (2020) 338–346.  doi:https://doi.org/10.1016/j.actamat.2020.06.015.
  * [22] V. Levitas, F. Lin, K. Pandey, S. Yesudhas, C. Park, Laws of high-pressure phase and nanostructure evolution and severe plastic flow, September 9, Research Square (2022) 29doi:https://doi.org/10.21203/rs.3.rs-1998605/v1.
  * [23] M. Aleksandrova, V. Blank, S. Buga, Phase transitions in Ge and Si under shear deformation at pressure up to 12 GPa conditions and P-T- gamma[shear] diagrams of these elements, Physics of the Solid State 35 (5) (1993) 1308–17.
  * [24] S. Yesudhas, V. Levitas, F. Lin, K. Pandey, J. Smith, Plastic strain induced phase transformations in micron and nano silicon, in preparation (2023).
  * [25] G. Voronin, C. Pantea, T. Zerda, L. Wang, Y. Zhao, In situ x-ray diffraction study of silicon at pressures up to 15.5 GPa and temperatures up to 1073 K, Physical Review B 68 (2) (2003) 020102.
  * [26] J. Pokluda, M. Černỳ, M. Šob, Y. Umeno, Ab initio calculations of mechanical properties: Methods and applications, Progress in Materials Science 73 (2015) 127–158.
  * [27] Y. Umeno, M. Černỳ, Effect of normal stress on the ideal shear strength in covalent crystals, Physical Review B 77 (10) (2008) 100101.
  * [28] R. Telyatnik, A. Osipov, S. Kukushkin, Ab initio modelling of nonlinear elastoplastic properties of diamond-like C, SiC, Si, Ge crystals upon large strains., Materials Physics & Mechanics 29 (1) (2016) 1–16.
  * [29] M. Černỳ, P. Řehák, Y. Umeno, J. Pokluda, Stability and strength of covalent crystals under uniaxial and triaxial loading from first principles, Journal of Physics: Condensed Matter 25 (3) (2012) 035401.
  * [30] N. A. Zarkevich, H. Chen, V. I. Levitas, D. D. Johnson, Lattice instability during solid-solid structural transformations under a general applied stress tensor: Example of Si I $\rightarrow$ Si II with metallization, Physical Review Letters 121 (16) (2018) 165701.
  * [31] H. Chen, N. A. Zarkevich, V. I. Levitas, D. D. Johnson, X. Zhang, Fifth-degree elastic energy for predictive continuum stress–strain relations and elastic instabilities under large strain and complex loading in silicon, NPJ Computational Materials 6 (1) (2020) 115.
  * [32] V. I. Levitas, Phase-field theory for martensitic phase transformations at large strains, International Journal of Plasticity 49 (2013) 85–118.
  * [33] P. Valentini, W. W. Gerberich, T. Dumitrică, Phase-Transition Plasticity Response in Uniaxially Compressed Silicon Nanospheres, Phys. Rev. Lett. 99 (2007) 175701.  doi:10.1103/PhysRevLett.99.175701.
  * [34] D. Chrobak, N. Tymiak, A. Beaber, O. Ugurlu, W. W. Gerberich, R. Nowak, Deconfinement leads to changes in the nanoscale plasticity of silicon, Nature Nanotechnology 6 (8) (2011) 480–484.
  * [35] H. Chen, V. I. Levitas, L. Xiong, Amorphization induced by 60∘ shuffle dislocation pileup against different grain boundaries in silicon bicrystal under shear, Acta Materialia 179 (2019) 287–295.
  * [36] L. C. Zhang, W. C. D. Cheong, Molecular dynamics simulation of phase transformations in monocrystalline silicon, in: V. Domnich, Y. Gogotsi (Eds.), High Pressure Surface Science and Engineering, Institute of Physics, Bristol, 2004, pp. 57–119.
  * [37] H. Chen, V. I. Levitas, D. Popov, N. Velisavljevic, Nontrivial nanostructure, stress relaxation mechanisms, and crystallography for pressure-induced Si-I $\rightarrow$ Si-II phase transformation, Nature Communications 13 (1) (2022) 982.  doi:10.1038/S41467-022-28604-1.
  * [38] V. I. Levitas, Phase field approach for stress- and temperature-induced phase transformations that satisfies lattice instability conditions. Part I. General theory, International Journal of Plasticity 106 (2018) 164–185.  doi:10.1016/J.IJPLAS.2018.03.007.
  * [39] H. Babaei, V. I. Levitas, Phase-field approach for stress- and temperature-induced phase transformations that satisfies lattice instability conditions. Part 2. Simulations of phase transformations Si I $\leftrightarrow$ Si II, International Journal of Plasticity 107 (2018) 223–245.  doi:10.1016/j.ijplas.2018.04.006.
  * [40] H. Babaei, V. I. Levitas, Effect of 60∘ dislocation on transformation stresses, nucleation, and growth for phase transformations between silicon I and silicon II under triaxial loading: Phase-field study, Acta Materialia 177 (2019) 178–186.
  * [41] H. Babaei, V. I. Levitas, Stress-Measure Dependence of Phase Transformation Criterion under Finite Strains: Hierarchy of Crystal Lattice Instabilities for Homogeneous and Heterogeneous Transformations, Phys. Rev. Lett. 124 (2020) 075701.  doi:10.1103/PhysRevLett.124.075701.
  * [42] H. Babaei, A. Basak, V. I. Levitas, Algorithmic aspects and finite element solutions for advanced phase field approach to martensitic phase transformation under large strains, Computational Mechanics 64 (4) (2019) 1177–1197.  doi:10.1007/S00466-019-01699-Y.
  * [43] V. I. Levitas, M. Javanbakht, Phase-field approach to martensitic phase transformations: Effect of martensite-martensite interface energy, International Journal of Materials Research 102 (6) (2011) 652–665.  doi:10.3139/146.110529.
  * [44] V. I. Levitas, A. V. Idesman, D. L. Preston, Microscale simulation of martensitic microstructure evolution, Physical Review Letters 93 (10) (2004) 105701\.
  * [45] A. V. Idesman, V. I. Levitas, D. L. Preston, J.-Y. Cho, Finite element simulations of martensitic phase transitions and microstructures based on a strain softening model, Journal of the Mechanics and Physics of Solids 53 (3) (2005) 495–523.
  * [46] S. E. Esfahani, I. Ghamarian, V. I. Levitas, P. C. Collins, Microscale phase field modeling of the martensitic transformation during cyclic loading of NiTi single crystal, International Journal of Solids and Structures 146 (2018) 80–96.
  * [47] H. Babaei, V. I. Levitas, Finite-strain scale-free phase-field approach to multivariant martensitic phase transformations with stress-dependent effective thresholds, Journal of the Mechanics and Physics of Solids 144 (2020) 104114.
  * [48] V. I. Levitas, M. Javanbakht, Phase transformations in nanograin materials under high pressure and plastic shear: nanoscale mechanisms, Nanoscale 6 (1) (2014) 162–166.
  * [49] M. Javanbakht, V. I. Levitas, Nanoscale mechanisms for high-pressure mechanochemistry: A phase field study, Journal of Materials Science 53 (19) (2018) 13343–13363.  doi:10.1007/S10853-018-2175-X.
  * [50] M. Javanbakht, V. I. Levitas, Phase field simulations of plastic strain-induced phase transformations under high pressure and large shear, Phys. Rev. B 94 (2016) 214104.  doi:10.1103/PhysRevB.94.214104.
  * [51] V. I. Levitas, S. E. Esfahani, I. Ghamarian, Scale-Free Modeling of Coupled Evolution of Discrete Dislocation Bands and Multivariant Martensitic Microstructure, Physical Review Letters 121 (20) (2018) 205701.  doi:10.1103/PHYSREVLETT.121.205701/FIGURES/4/MEDIUM.
  * [52] S. E. Esfahani, I. Ghamarian, V. I. Levitas, Strain-induced multivariant martensitic transformations: A scale-independent simulation of interaction between localized shear bands and microstructure, Acta Materialia (2020) 430–443.
  * [53] V. I. Levitas, Elastic model for stress tensor-induced martensitic transformation and lattice instability in silicon under large strains, Materials Research Letters 5 (8) (2017) 554–561.  doi:10.1080/21663831.2017.1362054.
  * [54] R. Hill, On Macroscopic Effects Of Heterogeneity In Elastoplastic Media At Finite Strain, Mathematical Proceedings of the Cambridge Philosophical Society 95 (3) (1984) 481–494.  doi:10.1017/S0305004100061818.
  * [55] V. I. Levitas, Some relations for finite inelastic deformation of microheterogeneous materials with moving discontinuity surfaces, in: A. Pineau, A. Zaoui (Eds.), IUTAM Symposium on Micromechanics of Plasticity and Damage of Multiphase Materials, Springer Netherlands, Dordrecht, 1996, pp. 313–320.
  * [56] H. Petryk, Macroscopic rate-variables in solids undergoing phase transformation, Journal of the Mechanics and Physics of Solids 46 (5) (1998) 873–894.  doi:10.1016/S0022-5096(97)00099-9.
  * [57] J. D. Schall, G. Gao, J. A. Harrison, Elastic constants of silicon materials calculated as a function of temperature using a parametrization of the second-generation reactive empirical bond-order potential, Physical Review B 77 (11) (2008) 115209.
  * [58] W. Bangerth, R. Hartmann, G. Kanschat, deal.IIA general-purpose object-oriented finite element library, ACM Transactions on Mathematical Software (TOMS) 33 (4) (8 2007).  doi:10.1145/1268776.1268779.
  * [59] M. A. Groeber, M. A. Jackson, DREAM. 3D: A digital representation environment for the analysis of microstructure in 3D, Integrating Materials and Manufacturing Innovation 3 (1) (2014) 56–72.
  * [60] Z. Zeng, Q. Zeng, M. Ge, B. Chen, H. Lou, X. Chen, J. Yan, W. Yang, H.-k. Mao, D. Yang, et al., Origin of plasticity in nanostructured silicon, Physical Review Letters 124 (18) (2020) 185701.
  * [61] Y. Xuan, L. Tan, B. Cheng, F. Zhang, X. Chen, M. Ge, Q. Zeng, Z. Zeng, Pressure-induced phase transitions in nanostructured silicon, The Journal of Physical Chemistry C 124 (49) (2020) 27089–27096.
  * [62] S. H. Tolbert, A. B. Herhold, L. E. Brus, A. Alivisatos, Pressure-induced structural transformations in Si nanocrystals: surface and shape effects, Physical Review Letters 76 (23) (1996) 4384.