A multiphase phase-field study of three-dimensional martensitic twinned
microstructures at large strains

Anup Basak1 and Valery I. Levitas2,3
1 Department of Mechanical Engineering, Indian Institute of Technology
Tirupati, Tirupati, A.P. 517506, India.
2 Departments of Aerospace Engineering, Mechanical Engineering, and Material
Science and Engineering, Iowa State University, Ames, IA 50011, USA.
3Division of Materials Science and Engineering, Ames Laboratory, Ames, IA
50011, USA.

A thermodynamically consistent multiphase phase-field approach for stress and
temperature-induced martensitic phase transformation at the nanoscale and
under large strains is developed. A total of $N$ independent order parameters
are considered for materials with $N$ variants, where one of the order
parameters describes $\sf A\leftrightarrow\sf M$ transformations and the
remaining $N-1$ independent order parameters describe the transformations
between the variants. A non-contradictory gradient energy is used within the
free energy of the system to account for the energies of the interfaces. In
addition, a non-contradictory kinetic relationships for the rate of the order
parameters versus thermodynamic driving forces is suggested. As a result, a
system of consistent coupled Ginzburg-Landau equations for the order
parameters are derived. The crystallographic solution for twins within twins
is presented for the cubic to tetragonal transformations. A 3D complex twins
within twins microstructure is simulated using the developed phase-field
approach and a large-strain-based nonlinear finite element method. A
comparative study between the crystallographic solution and the simulation
result is presented.

Keywords: Multiphase phase-field approach; Martensitic transformations; Twins
within twins; Crystallographic solution; Large strains; Finite element method.

##  1 Introduction

Martensitic transformations and microstructures. Martensitic transformations
(MTs) are diffusionless solid-solid phase transformations observed in many
metallic and nonmetallic crystalline solids, minerals, and various compounds,
where a parent phase called austenite (high temperature phase) transforms into
the product phase called martensite (low temperature phase) [2, 3]. The
martensitic phase (here denoted by $\sf M$) has lower crystallographic
symmetry than the austenite phase (here denoted by $\sf A$), and generally has
multiple variants. Very complex microstructures such as austenite-twinned
martensite, twins within twins, twins within twins within twins, wedge,
X-interfaces, etc., are observed within the materials undergoing MTs [4, 2, 3,
5, 6]. The evolution of such microstructures play a central role in, for
example, strengthening of steel, shape memory effect in various alloys,
ferromagnetic effect, caloric effects, etc [7, 2].

In continuum theories for MTs, such phase-changing materials are modeled as
nonlinear elastic materials having multiple wells in the free energy density
function [5, 2, 3, 8, 9]. Fine twinned microstructures associated with
austenite-martensite interfaces, which were observed under the microscopes
[10, 11, 12, 13], are usually obtained as the minimizers of such non-convex
energies within the continuum theories [5]. The analytical crystallographic
solutions for twins between a pair of variants and the austenite-twinned
martensite interfaces are well-known within the small as well as finite
deformation theories [5, 2, 3, 14]. These solutions have been further used for
more complex wedge microstructures [6, 15] and X-interfaces [16, 17]. Another
important complex microstructure is twins within twins [4, 2] for which the
general crystallographic equations were established in [2, 11, 12]. Though the
governing equations are well-known for twins within twins, the analytical
solutions to such microstructures are still missing to the best of our
knowledge.

Phase-field approach to MTs. The phase-field approach based on the Ginzburg-
Landau equations [18] (similar to Allen-Cahn’s approach [19]), which provides
an ideal framework for studying the MTs, have been widely used for studying
nucleation, growth of the phases, and evolution of complex microstructures
[20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 32, 31, 33, 34, 35, 36, 37, 38,
39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50]. A set of sufficiently smooth
scalar internal variables called the order parameters, are used to describe
the phases. The volume fraction based (e.g. [51, 52, 53, 54, 55, 56, 57, 58,
59]) or the transformation strains based (e.g. [20, 21, 24, 25, 27]) order
parameters have been used. Important requirements for the interpolation
functions are formulated in [46, 47, 48, 60, 49]. Within the multiphase phase-
field approaches, the order parameters should be constrained to some specified
surfaces in order to control the transformation paths. To control the
transformation paths, various constraint hypersurfaces such as hypersphere
[61], planar surfaces [62], and straight lines [63, 21] have been used; see
[21] for a review. The double-well [20, 21, 27, 63, 61] or double-obstacle
[23] based thermal energies are usually used. The free energies are considered
to be smooth functions of the order parameters, and the transformation strains
are accepted as linear [51, 53, 54, 56, 55, 56, 57, 58, 59] or nonlinear
functions [20, 21, 22, 64, 24, 25] of the order parameters smoothly varying
between all the phases, in particular, satisfying some additional requirements
[46, 47]. First large-strain phase-field theory and computational approaches
were presented in [38, 48, 22]. They utilized the methods of repetitive
superposition of large strains, developed by Levin in [67, 65, 66] for
viscoelastic materials, extended to materials with phase transformations.

A gradient (of the order parameters) based nonlocal energy is considered which
introduces a finite interface width between the phases. The time evolution of
these order parameters describing the kinetics of the PTs is derived using the
laws of thermodynamics yielding a system of coupled Ginzburg-Landau equations.
The interfacial stresses, consisting of the elastic and structural components
and which play an important role in the nucleation of the phases and also in
their kinetics and growth, have been considered [36, 68, 69]. A detailed
comparison between these various multiphase phase-field approaches to MTs is
presented in [21].

The multiphase phase-field model for studying the multivariant MTs developed
by the authors in [21] yields non-contradictory results for a two-variant
system. However, for a system with more than two variants, some contradictions
have been observed in relation to the gradient energy and the system of
kinetic equations. One of the aims of this work is to discuss those issues
from that model and present a non-contradictory multiphase phase-field model
for MTs. A gradient energy proposed therein simplifies consistently for a two-
variant system and matches with the well-established result (see [21] for the
discussion), but a contradiction is observed when the system contains more
than two variants as discussed in Sec. 2.3. An alternative form of the
gradient energy has been used here which has similarities with the gradient
energy used in [62, 24, 23, 54] and yields non-contradictory results for any
number of variants. In the present model, we however multiply this gradient
term with the determinant of the total deformation gradient while determining
the total system energy to ensure an appropriate form of the structural stress
tensor, which was however not considered in [62, 24, 23, 54]. Furthermore, we
point out in Sec. 2.5 that the coupled kinetic equations for the order
parameters are non-contradictory for a two-variant system, but leads to
contradictions for an $N$-variant system for $N>2$, if the kinetic
coefficients are assumed to be constants. We thus introduce a system of
kinetic equations with kinetic coefficients which are step-wise functions of
the order parameters and the driving forces, which is motivated from Ref.
[56].

Contribution of the paper. The contributions of this paper are mainly two-
fold:

i) We present a thermodynamically consistent nanoscale phase-field approach
for multivariant MTs considering non-contradictory gradient energies and the
local energies including the barrier, chemical, and elastic energy, as well as
energies penalizing the multiphase junctions, while the deviations of the
transformation paths for $\sf A\leftrightarrow\sf M$ and $\sf
M_{i}\leftrightarrow\sf M_{j}$ PTs from the specified paths are appropriately
controlled. The issues with the existing gradient energy models are discussed.
Furthermore, a consistent kinetic model for coupled Ginzburg-Landau equations
is derived, and the issues with the existing kinetic models are discussed. The
present model can be used for MTs with any number of variants.

ii) The crystallographic solution for the twins within twins microstructure
arising in cubic to tetragonal MTs are derived. The evolution and formation of
3D twins within twins microstructures in a single grain is studied using the
present phase-field approach. The simulation results are in good agreement
with the crystallographic solution and the experimental results.

Notations. The multiplication and the inner product between two arbitrary
second order tensors ${\boldsymbol{A}}$ and ${\boldsymbol{D}}$ are denoted by
$({\boldsymbol{A}}\cdot{\boldsymbol{D}})_{ab}=A_{ac}D_{cb}$ and
${\boldsymbol{A}}:{\boldsymbol{D}}=A_{ab}D_{ba}$, respectively, where $A_{ab}$
and $D_{ab}$ are the components of the tensors in a right-handed orthonormal
Cartesian basis
$\\{\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\\}$. The repeated
indices imply Einstein’s summation. The Euclidean norm of $\boldsymbol{A}$ is
defined by $|\boldsymbol{A}|=\sqrt{{\boldsymbol{A}}:{\boldsymbol{A}}^{T}}$.
The second-order identity tensor is denoted by ${\boldsymbol{I}}$.
$\boldsymbol{A}^{T}$, $tr\,\boldsymbol{A}$, $det\,\boldsymbol{A}$,
$sym(\boldsymbol{A})$, and $skew(\boldsymbol{A})$ denote the transpose, trace,
determinant, symmetric part, and skew part of $\boldsymbol{A}$, respectively.
For an invertible tensor $\boldsymbol{A}$, its inverse is denoted by
$\boldsymbol{A}^{-1}$. The tensor or dyadic product between two arbitrary
vectors $\boldsymbol{a}$ and $\boldsymbol{b}$ is denoted by
$\boldsymbol{a}\otimes\boldsymbol{b}$. The reference, stress-free
intermediate, and deformed or current configurations are denoted by
$\Omega_{0}$, $\Omega_{t}$, and $\Omega$, respectively. The volumes in the
reference and current configurations are denoted by $V_{0}$ and $V$, and their
external boundaries are denoted by $S_{0}$ and $S$, respectively. The symbols
$\nabla_{0}(\cdot)$ and $\nabla(\cdot)$ denote the gradient operators in
$\Omega_{0}$ and $\Omega$, respectively. The Laplacian operators in
$\Omega_{0}$ and $\Omega$ are designated by
$\nabla_{0}^{2}:=\nabla_{0}\cdot\nabla_{0}$ and
$\nabla^{2}:=\nabla\cdot\nabla$, respectively. The symbol $:=$ implies
equality by definition.

##  2 Coupled mechanics and phase-field model

In this section, we describe our multiphase phase-field model, mostly based on
[21] but with further development in terms of non-contradictory gradient
energy and kinetic relationships between the rate of the order parameters and
conjugated thermodynamical driving forces. The free energy used in the model
is presented. The coupled elasticity equations and a new system of coupled
Ginzburg-Landau equations are derived. A comparison of the present model with
the previous models from the literature is also presented.

###  2.1 Order parameters

For the MTs in a system with austenite and $N$ martensitic variants we
consider $N+1$ order parameters
$\eta_{0},\eta_{1},\ldots,\eta_{i},\eta_{j},\ldots,\eta_{N}$, where $\eta_{0}$
describes ${\sf A}\leftrightarrow{\sf M}$ transformations such that
$\eta_{0}=0$ in $\sf A$ and $\eta_{0}=1$ in $\sf M$, and $\eta_{i}$ (for
$i=1,\ldots,N$) describes the variant ${\sf M}_{i}$ such that $\eta_{i}=1$ in
${\sf M}_{i}$ and $\eta_{i}=0$ in ${\sf M}_{j}$ for all $j\neq i$. Such
descriptions for the order parameters were introduced by the authors in
earlier work [21]. The order parameters $\eta_{1},\,\eta_{2},\ldots,\eta_{N}$
are constrained to lie on a plane by satisfying (see [21] for details)

$\sum_{i=1}^{N}\eta_{i}=1.$ (2.1)

###  2.2 Kinematics

The position vector of a particle in the deformed configuration $\Omega$ at
time instance $t$ is given by
$\boldsymbol{r}(\boldsymbol{r}_{0},t)=\boldsymbol{r}_{0}+\boldsymbol{u}(\boldsymbol{r}_{0},t)$,
where $\boldsymbol{r}_{0}$ is the position vector of that particle in
$\Omega_{0}$, and $\boldsymbol{u}$ is the displacement vector. We consider the
following multiplicative decomposition for the total deformation gradient
$\boldsymbol{F}$ [68]:

${\boldsymbol{F}}:=\nabla_{0}\boldsymbol{r}=\boldsymbol{F}_{e}\cdot\boldsymbol{F}_{t}=\boldsymbol{V}_{e}\cdot\boldsymbol{R}\cdot\boldsymbol{U}_{t},$
(2.2)

where the subscripts $e$ and $t$ denote elastic and transformational parts,
respectively, and $\boldsymbol{F}_{e}$ and $\boldsymbol{F}_{t}$ are the
elastic and transformational parts of $\boldsymbol{F}$. The tensors
$\boldsymbol{F}_{e}$ and $\boldsymbol{F}_{t}$ are also decomposed into
$\boldsymbol{F}_{e}=\boldsymbol{V}_{e}\cdot\boldsymbol{R}_{e}$ and
$\boldsymbol{F}_{t}=\boldsymbol{R}_{t}\cdot\boldsymbol{U}_{t}$, respectively,
where $\boldsymbol{V}_{e}$ is the left elastic stretch tensor (symmetric),
$\boldsymbol{U}_{t}$ is the right transformation stretch tensor (symmetric),
and $\boldsymbol{R}_{e}$ and $\boldsymbol{R}_{t}$ are rotations. In Eq. (2.2)
we have used $\boldsymbol{R}=\boldsymbol{R}_{e}\cdot\boldsymbol{R}_{t}$. We
denote $J=det\,\boldsymbol{F}$, $J_{t}=det\,\boldsymbol{F}_{t}$, and
$J_{e}=det\,\boldsymbol{F}_{e}$. Hence, by Eq. (2.2), $J=J_{t}J_{e}$. The
Lagrangian total and elastic strain tensors are defined as

${\boldsymbol{E}}:=0.5(\boldsymbol{C}-{\boldsymbol{I}}),\quad\text{and}\quad{\boldsymbol{E}}_{e}:=0.5(\boldsymbol{C}_{e}-{\boldsymbol{I}}),$
(2.3)

respectively, where
$\boldsymbol{C}={\boldsymbol{F}}^{T}\cdot{\boldsymbol{F}}$, and
$\boldsymbol{C}_{e}={\boldsymbol{F}}_{e}^{T}\cdot{\boldsymbol{F}}_{e}$. We
define another measure of the total and elastic strain tensors as

${\boldsymbol{b}}=0.5(\boldsymbol{B}-{\boldsymbol{I}}),\quad\text{and
}\quad{\boldsymbol{b}}_{e}=0.5(\boldsymbol{B}_{e}-{\boldsymbol{I}}),$ (2.4)

respectively, where
$\boldsymbol{B}={\boldsymbol{F}}\cdot{\boldsymbol{F}}^{T}=\boldsymbol{V}^{2}$,
$\boldsymbol{B}_{e}={\boldsymbol{F}}_{e}\cdot{\boldsymbol{F}}_{e}^{T}=\boldsymbol{V}_{e}^{2}$,
and $\boldsymbol{V}=\sqrt{{\boldsymbol{B}}}$ is the left total stretch tensor.

Kinematic model for $\boldsymbol{F}_{t}$. We consider $\boldsymbol{F}_{t}$ as
a linear combination of the Bain strains multiplied by the interpolation
functions related to the order parameters [21]:

${\boldsymbol{F}}_{t}={\boldsymbol{U}}_{t}=\boldsymbol{I}+\sum_{i=1}^{N}(\boldsymbol{U}_{ti}-\boldsymbol{I})\varphi(a_{\varepsilon},\eta_{0})\phi_{i}(\eta_{i}),$
(2.5)

where $\boldsymbol{U}_{ti}$ is the Bain stretch tensor for $\mathsf{M}_{i}$.
We take the interpolation functions $\varphi(a_{\varepsilon},\eta_{0})$ and
$\phi_{i}(\eta_{i})$ as [21]

$\displaystyle\varphi(a_{\varepsilon},\eta_{0})$ $\displaystyle=$
$\displaystyle
a_{\varepsilon}\eta_{0}^{2}+(4-2a_{\varepsilon})\eta_{0}^{3}+(a_{\varepsilon}-3)\eta_{0}^{4},\quad\text{and}$
$\displaystyle\phi_{i}(\eta_{i})$ $\displaystyle=$
$\displaystyle\eta_{i}^{2}(3-2\eta_{i})\quad\text{for all }i=1,2,\ldots,N,$
(2.6)

respectively, which satisfy the following conditions derived from the
requirement of thermodynamic equilibrium of the homogeneous phases:

$\displaystyle\varphi(a_{\varepsilon},0)=0,\quad\varphi(a_{\varepsilon},1)=1,\quad\text{and}\quad\frac{\partial\varphi(a_{\varepsilon},0)}{\partial\eta_{0}}=\frac{\partial\varphi(a_{\varepsilon},1)}{\partial\eta_{0}}=0;$
$\displaystyle\phi_{i}(0)=0,\quad\phi_{i}(1)=1,\quad\text{and}\quad\frac{\partial\phi_{i}(\eta_{i}=0)}{\partial\eta_{i}}=\frac{\partial\phi_{i}(\eta_{i}=1)}{\partial\eta_{i}}=0\quad\text{for
all }i=1,2,\ldots,N.$ (2.7)

The constant $a_{\varepsilon}$ in Eq. (2.6) must be within the range $0\leq
a_{\varepsilon}\leq 6$ [68].

###  2.3 Free energy

We assume the Helmholtz free energy per unit mass of the body as [21, 68]:

$\displaystyle\psi(\boldsymbol{F},\theta,\eta_{0},\eta_{i},\nabla\eta_{0},\nabla\eta_{i})=\psi^{l}(\boldsymbol{F}_{e},\theta,\eta_{0},\eta_{i})+J\psi^{\nabla}(\eta_{0},\nabla\eta_{0},\nabla\eta_{i}),$
(2.8)

where $\psi^{l}(\boldsymbol{F}_{e},\theta,\eta_{0},\eta_{i})$ is the local
part of the free energy density and
$\psi^{\nabla}(\eta_{0},\nabla\eta_{0},\nabla\eta_{i})$ is the gradient based
nonlocal energy accounting for the energies of all the interfaces. We have
taken $\psi^{l}$ as

$\displaystyle\psi^{l}(\boldsymbol{F}_{e},\theta,\eta_{0},\eta_{i})$
$\displaystyle=$
$\displaystyle\frac{J_{t}}{\rho_{0}}\psi_{e}(\boldsymbol{F}_{e},\theta,\eta_{0},\eta_{i})+J\breve{\psi}^{\theta}(\theta,\eta_{0},\eta_{i})+\tilde{\psi}^{\theta}(\theta,\eta_{0},\eta_{i})+\psi_{p}(\eta_{0},\eta_{i}),$
(2.9)

where $\psi_{e}$ is the strain energy per unit volume in $\Omega_{t}$,
$\breve{\psi}^{\theta}$ is the barrier energy related to ${\sf
A}\leftrightarrow{\sf M}$ PT and all the variant$\leftrightarrow$variant
transformations, $\tilde{\psi}^{\theta}$ is the thermal/chemical energy for
${\sf A}\leftrightarrow{\sf M}$ transformations, $\psi_{p}$ penalizes various
triple and higher junctions between all the phases and also accounts for the
penalization in energies for the deviation of the transformation paths from
the assigned ones, $\theta>0$ is the absolute temperature, and $\rho_{0}$ is
the density of the solid in $\Omega_{0}$. In Eqs. (2.8) and (2.9), the barrier
energy and the gradient energy are multiplied by $J$ following [68] in order
to obtain the desired expression for the structural stresses; see Eq. (2.29).
Any material property $B$ for each particle of the body are determined using
[21]

$B(\eta_{0},\eta_{i},\theta,{\boldsymbol{F}})=B_{0}(1-\varphi(a,\eta_{0}))+\varphi(a,\eta_{0})\sum_{i=1}^{N}B_{i}\phi_{i}(\eta_{i}),$
(2.10)

where $B_{0}$ and $B_{i}$ are the material properties of the homogeneous
phases $\mathsf{A}$ and $\mathsf{M}_{i}$, respectively, $\varphi(a,\eta_{0})$
has the same expression as $\varphi(a_{\varepsilon},\eta_{0})$ given by Eq.
(2.6)1 when $a_{\varepsilon}$ is replaced by $a$ therein.

The expressions for all the energies introduced in Eqs. (2.8) and (2.9) are
given below:

(i) Strain energy: We consider $\psi_{e}$ as [21]

$\psi_{e}=0.5\boldsymbol{E}_{e}:\hat{\boldsymbol{\mathcal{C}}}_{e}(\eta_{0},\eta_{i}):\boldsymbol{E}_{e},$
(2.11)

where the fourth-order elastic modulus tensor at any particle is taken
following Eq. (2.10) as [21]

$\hat{\boldsymbol{\mathcal{C}}}_{e}(\eta_{0},\eta_{i})=(1-\varphi(a,\eta_{0}))\hat{\boldsymbol{\mathcal{C}}}_{(e)0}+\varphi(a,\eta_{0})\sum_{i=1}^{N}\phi_{i}(\eta_{i})\hat{\boldsymbol{\mathcal{C}}}_{(e)i},$
(2.12)

and $\hat{\boldsymbol{\mathcal{C}}}_{(e)0}$ and
$\hat{\boldsymbol{\mathcal{C}}}_{(e)i}$ are the fourth order elastic modulus
tensors of $\sf A$ and ${\mathsf{M}}_{i}$, respectively.

(ii) Barrier energy: The total energy of the barriers between $\sf A$ and $\sf
M$ and between all the variants is [21]

$\breve{\psi}^{\theta}=A_{0M}\,\eta_{0}^{2}(1-\eta_{0})^{2}+\varphi(a_{b},\eta_{0})\tilde{A}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\eta_{i}^{2}\,\eta_{j}^{2},$
(2.13)

where $A_{0M}$ and $\tilde{A}$ are the coefficients for the barrier energies
between ${\sf A}$ and ${\sf M}$, and ${\sf M}_{i}$ and ${\sf M}_{j}$ (for all
$i\neq j$), respectively.

(iii) Thermal energy: The thermal energy of a particle undergoing ${\sf
A}\leftrightarrow{\sf M}$ PTs is taken as [46, 47, 48, 21]

$\tilde{\psi}^{\theta}=\psi_{0}^{\theta}(\theta)+\eta_{0}^{2}(3-2\eta_{0})\,\Delta\psi^{\theta}(\theta),\quad\text{where
}\Delta\psi^{\theta}=-\Delta s_{0M}(\theta-\theta_{e}),$ (2.14)

$\psi_{0}^{\theta}$ is the thermal energy of $\sf A$,
$\Delta\psi^{\theta}=\psi^{\theta}_{M}-\psi^{\theta}_{0}$ is the thermal
energy difference between $\sf A$ and $\sf M$ phases, $\Delta
s_{0M}=s_{M}-s_{0}$, $s_{0}$ and $s_{M}$ are the entropies of $\sf A$ and $\sf
M$, respectively, per unit volume, and $\theta_{e}$ is the thermodynamic
equilibrium temperature between $\sf A$ and $\sf M$ phases.

(iv) Gradient energy: We consider the gradient energy, taking all the
interfacial energies into consideration, as

$\psi^{\nabla}=\frac{\beta_{0M}}{2\rho_{0}}{\left|{\nabla\eta_{0}}\right|^{2}}+\frac{1}{2\rho_{0}}\tilde{\varphi}(\eta_{0},a_{\beta},a_{c})\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}{\beta_{ij}}\nabla\eta_{i}\cdot\nabla\eta_{j},$
(2.15)

where the interpolation function $\tilde{\varphi}$ is taken as [21]

$\tilde{\varphi}(a_{\beta},a_{c},\eta_{0})=a_{c}+{a_{\beta}}{\eta_{0}^{2}}-2[{a_{\beta}}-2(1-a_{c})]{\eta_{0}^{3}}+[{a_{\beta}}-3(1-a_{c})]{\eta_{0}^{4}},$
(2.16)

where the constant is taken as $0<a_{c}\ll 1$, and the purpose of considering
it in Eq. (2.16) is discussed in [21]. In Eq. (2.15), $\beta_{0M}$ and
$\beta_{ij}$ are the gradient energy coefficients for ${\sf A}$-${\sf M}$ and
${\sf M}_{i}$-${\sf M}_{j}$ interfaces, respectively. In Eq. (2.16)
$a_{\beta}$ and $a_{c}$ are constant parameters. The gradient energy similar
to Eq. (2.15) was earlier used in [54, 62]. Notably, the authors earlier
introduced another form of the gradient energy in [21] given by

$\displaystyle\psi^{\nabla}=\frac{1}{2\rho_{0}}\left[{\beta_{0M}}{\left|{\nabla\eta_{0}}\right|^{2}}+\sum_{i=1}^{N}\sum_{j=1,\neq
i}^{N}\frac{\beta_{ij}}{8}|\nabla\eta_{i}-\nabla\eta_{j}|^{2}\tilde{\varphi}(\eta_{0},a_{\beta},a_{0})\right].$
(2.17)

This energy given by Eq. (2.17) simplifies to

$\psi^{\nabla}=\frac{1}{2\rho_{0}}{\beta_{0M}}{\left|{\nabla\eta_{0}}\right|^{2}}+\frac{\tilde{\varphi}(\eta_{0},a_{\beta},a_{0})}{8\rho_{0}}\left(\beta_{12}{\left|\nabla\eta_{1}-\nabla\eta_{2}\right|^{2}}+\beta_{21}{\left|\nabla\eta_{2}-\nabla\eta_{1}\right|^{2}}\right),$
(2.18)

for a system with two variants, where applying the constraint
$\eta_{1}+\eta_{2}=1$ and $\beta_{12}=\beta_{21}$ due to the symmetry [68] we
further simplify it to

$\psi^{\nabla}=\frac{1}{2\rho_{0}}{\beta_{0M}}{\left|{\nabla\eta_{0}}\right|^{2}}+\frac{\tilde{\varphi}(\eta_{0},a_{\beta},a_{0})}{2\rho_{0}}\beta_{12}{\left|\nabla\eta_{1}\right|^{2}},$
(2.19)

which is consistent with the results of earlier models; see e.g. [68] and the
references therein. We would like to mention that the coefficients
$\beta_{ij}$ in Eq. (2.17) for the variant pairs which are in twin
relationships would be much smaller than that for the variant pairs which are
not in twin relationships. For a system with three variants, Eq. (2.17)
reduces to

$\displaystyle\psi^{\nabla}$ $\displaystyle=$
$\displaystyle\frac{{\beta_{0M}}}{2\rho_{0}}{\left|{\nabla\eta_{0}}\right|^{2}}+\frac{1}{8\rho_{0}}\left[(\beta_{12}+\beta_{13})\left|{\nabla\eta_{1}}\right|^{2}+(\beta_{12}+\beta_{23})\left|{\nabla\eta_{2}}\right|^{2}+(\beta_{13}+\beta_{23})\left|{\nabla\eta_{3}}\right|^{2}-\right.$
(2.20)
$\displaystyle\left.2\,(\beta_{12}\nabla\eta_{1}\cdot\nabla\eta_{2}+\beta_{23}\nabla\eta_{2}\cdot\nabla\eta_{3}+\beta_{13}\nabla\eta_{1}\cdot\nabla\eta_{3})\right]\tilde{\varphi}(\eta_{0},a_{\beta},a_{0}).$

Let us consider a region for such a three-variant system where only $\sf
M_{1}$ and $\sf M_{2}$ coexist and $\sf M_{3}$ is absent, i.e. $\eta_{3}=0$
and $\eta_{1}+\eta_{2}=1$. The gradient energy given by Eq. (2.20) in that
region is rewritten by applying these conditions as

$\displaystyle\psi^{\nabla}$ $\displaystyle=$
$\displaystyle\frac{{\beta_{0M}}}{2\rho_{0}}{\left|{\nabla\eta_{0}}\right|^{2}}+\frac{1}{8\rho_{0}}(\beta_{23}+\beta_{13}+4\beta_{12})\left|{\nabla\eta_{1}}\right|^{2}\tilde{\varphi}(\eta_{0},a_{\beta},a_{0}).$
(2.21)

Obviously, the energy parameters $\beta_{23}$ and $\beta_{13}$ are going to
influence the interfacial energy $\sf M_{1}$ and $\sf M_{2}$ variants, which
is nonphysical; see also [70] for an analysis. We have shown that the gradient
energy given by Eq. (2.17) yields a nonphysical contribution for an interface
between two variants from the gradient coefficients which are not related to
that interface, and hence this form of gradient energy is not acceptable.
However, the energy given by Eq. (2.15) is non-contradictory for any number of
variants, and hence it is accepted in this paper.

(v) Penalization energy for junctions: We penalize the triple and higher
junctions between all the phases and the deviations of the transformation
paths between the variants using [21]

$\displaystyle\psi_{p}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}K_{ij}(\eta_{i}+\eta_{j}-1)^{2}\eta_{i}^{2}\eta_{j}^{2}+[1-\varphi(a_{K},\eta_{0})]\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}K_{0ij}\eta_{0}^{2}\eta_{i}^{2}\eta_{j}^{2}+\sum_{i=1}^{N-2}\sum_{j=i+1}^{N-1}\sum_{k=j+1}^{N}K_{ijk}\eta_{i}^{2}\eta_{j}^{2}\eta_{k}^{2}+$
$\displaystyle[1-\varphi(a_{K},\eta_{0})]\sum_{i=1}^{N-2}\sum_{j=i+1}^{N-1}\sum_{k=j+1}^{N}K_{0ijk}\eta_{0}^{2}\eta_{i}^{2}\eta_{j}^{2}\eta_{k}^{2}+\sum_{i=1}^{N-3}\sum_{j=i+1}^{N-2}\sum_{k=j+1}^{N-1}\sum_{l=k+1}^{N}K_{ijkl}\eta_{i}^{2}\eta_{j}^{2}\eta_{k}^{2}\eta_{l}^{2},\quad\text{where}$
$\displaystyle K_{ii}$ $\displaystyle=$ $\displaystyle
K_{0ii}=K_{iji}=K_{iik}=K_{0iji}=K_{0iik}=K_{ijil}=K_{iikl}=K_{ijjl}=K_{ijki}=K_{ijkk}=0.$
(2.22)

In Eq. (2.3), the parameter $K_{ij}\geq 0$ controls the penalization of the
deviation of the ${\sf M}_{j}\leftrightarrow{\sf M}_{i}$ transformation path
from the straight line $\eta_{j}+\eta_{i}=1$ for all $\eta_{k}=0$ and $k\neq
j,i$; the constant coefficients $K_{0ij}\geq 0$, $K_{ijk}\geq 0$,
$K_{0ijk}\geq 0$, and $K_{ijkl}\geq 0$ are related to the penalization of the
junctions between ${\sf A}$-${\sf M}_{i}$-${\sf M}_{j}$, ${\sf M}_{i}$-${\sf
M}_{j}$-${\sf M}_{k}$, ${\sf A}$-${\sf M}_{i}$-${\sf M}_{j}$-${\sf M}_{k}$,
and ${\sf M}_{i}$-${\sf M}_{j}$-${\sf M}_{k}$-${\sf M}_{l}$, respectively.

In the absence of all the penalty terms, i.e. when
$K_{ij}=K_{0ij}=K_{0ik}=K_{0jk}=K_{ijk}=K_{0ijk}=0$, we can show that for a
martensitic region ($\eta_{0}=1$) with three variants, say, $\sf M_{1}$, $\sf
M_{2}$ and $\sf M_{3}$, the barrier energy (see Eq. (2.13)) in the middle of
the triple junction region, i.e., at the point with
$\eta_{1}=\eta_{2}=\eta_{3}=1/3$ is $3\times\bar{A}/81=\bar{A}/27$ which is
less than the barrier energy $\bar{A}/16$ at the middle line of any variant-
variant interface (a line with, say, $\eta_{1}=\eta_{2}=1/2$ and
$\eta_{3}=0$). Hence when $K_{123}\neq 0$ the total energy at a martensitic
particle with $\eta_{1}=\eta_{2}=\eta_{3}=1/3$ is

$\left.E_{TJ}\right|_{\eta_{1}=\eta_{2}=\eta_{3}=1/3}=\frac{\bar{A}}{27}+\frac{K_{123}}{729}.$
(2.23)

It is to be noted that Tóth et al. [71] and Bollada et al. [70] considered
barrier energy similar to ours given by Eq. (2.13) but with a common
multiplication factor to incorporate higher energy at the junction region as
compared to the respective interface regions. In this paper, we have, however,
followed a different and simpler approach for that purpose, where we have
introduced the penalty terms in the free energy, and by varying the
coefficients $K_{0ij},\,K_{0ik},\,K_{0jk},\,K_{ijk},\,K_{0ijk}$ we can control
the energy of all the junction regions. For example, by tuning the parameter
$K_{123}$ in Eq. (2.23), we can make the barrier energy height at the junction
region higher than the barrier energy in the interfacial region. A
quantitative comparison between our formulation and the approach in [71, 70]
is however not given here.

###  2.4 Governing equations

We now present the governing equations. Applying the principle of balance of
linear and angular momentum, and the first and second law of thermodynamics,
and using an approach similar to [68, 21], we derive the mechanical
equilibrium equation, and the dissipation inequalities listed below.

####  2.4.1 Mechanical equilibrium equations and stresses

The mechanical equilibrium equation is given by [68, 21]

$\nabla_{0}\cdot\boldsymbol{P}=\boldsymbol{0}\quad\text{in
}\Omega_{0},\quad\text{or}\quad\nabla\cdot\boldsymbol{\sigma}=\boldsymbol{0}\quad\text{in
}\Omega,$ (2.24)

where the body forces and inertia are neglected, $\boldsymbol{P}$ is the total
first Piola-Kirchhoff stress tensor, and $\boldsymbol{\sigma}$ is the total
Cauchy stress tensor which is symmetric. The total stresses are composed of
their respective elastic and structural parts [21]:

$\boldsymbol{P}=\boldsymbol{P}_{e}+\boldsymbol{P}_{st},\quad\text{and}\quad\boldsymbol{\sigma}=\boldsymbol{\sigma}_{e}+\boldsymbol{\sigma}_{st}.$
(2.25)

The elastic stresses are given by [21, 68]

$\boldsymbol{P}_{e}=J_{t}\boldsymbol{F}_{e}\cdot\hat{\boldsymbol{S}}_{e}\cdot\boldsymbol{F}_{t}^{-T},\quad\text{and}\quad\boldsymbol{\sigma}_{e}=J_{e}^{-1}\boldsymbol{F}_{e}\cdot\hat{\boldsymbol{S}}_{e}\cdot\boldsymbol{F}_{e}^{T},$
(2.26)

where
$\hat{\boldsymbol{S}}_{e}=\frac{\partial\psi_{e}(\boldsymbol{E}_{e})}{\partial\boldsymbol{E}_{e}}$.
For an isotropic elastic response, $\boldsymbol{P}_{e}$ and
$\boldsymbol{\sigma}_{e}$ can alternatively be expressed as

$\boldsymbol{P}_{e}=J_{t}\boldsymbol{V}_{e}^{2}\cdot\frac{\partial\psi_{e}(\boldsymbol{b}_{e})}{\partial\boldsymbol{b}_{e}}\cdot\boldsymbol{F}^{-T},\quad\text{and}\quad\boldsymbol{\sigma}_{e}=J^{-1}_{e}\boldsymbol{V}_{e}^{2}\cdot\frac{\partial\psi_{e}(\boldsymbol{b}_{e})}{\partial\boldsymbol{b}_{e}},$
(2.27)

respectively. The general forms of structural stress tensors are given by [21,
68]

$\displaystyle{\boldsymbol{P}}_{st}$ $\displaystyle=$ $\displaystyle
J\rho_{0}(\breve{\psi}^{\theta}+\psi^{\nabla}){\boldsymbol{F}}^{-T}-J\rho_{0}\left(\nabla\eta_{0}\otimes\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{0}}+\sum_{i=1}^{N}\nabla\eta_{i}\otimes\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{i}}\right)\cdot\boldsymbol{F}^{-T},\quad\text{and}$
$\displaystyle\boldsymbol{\sigma}_{st}$ $\displaystyle=$
$\displaystyle\rho_{0}(\breve{\psi}^{\theta}+\psi^{\nabla})\boldsymbol{I}-\rho_{0}\left(\nabla\eta_{0}\otimes\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{0}}+\sum_{i=1}^{N}\nabla\eta_{i}\otimes\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{i}}\right).$
(2.28)

Using the gradient energy given by Eq. (2.15) in Eq. (2.28)1,2, the exact form
of the structural stresses are obtained as

$\displaystyle{\boldsymbol{P}}_{st}$ $\displaystyle=$ $\displaystyle
J\rho_{0}(\breve{\psi}^{\theta}+\psi^{\nabla}){\boldsymbol{F}}^{-T}-J\beta_{0M}\nabla\eta_{0}\otimes\nabla\eta_{0}\cdot{\boldsymbol{F}}^{-T}-J\tilde{\varphi}\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\beta_{ij}\nabla\eta_{i}\otimes\nabla\eta_{j}\right)\cdot{\boldsymbol{F}}^{-T},\quad\text{and}$
$\displaystyle\boldsymbol{\sigma}_{st}$ $\displaystyle=$
$\displaystyle\rho_{0}(\breve{\psi}^{\theta}+\psi^{\nabla})\boldsymbol{I}-\beta_{0M}\nabla\eta_{0}\otimes\nabla\eta_{0}-\tilde{\varphi}\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\beta_{ij}\nabla\eta_{i}\otimes\nabla\eta_{j}\right).$
(2.29)

###  2.5 Dissipation inequality and Ginzburg-Landau equations

The dissipation inequalities for the order parameter $\eta_{0}$ and also for
the order parameter $\eta_{1},\ldots,\eta_{N}$ are obtained as (see [21] for
the details)

$D_{0}=\dot{\eta}_{0}X_{0}\geq 0\quad\text{and}\quad
D_{M}=\sum_{i=1}^{N}\dot{\eta}_{i}X_{i}\geq 0,$ (2.30)

respectively, where the conjugate ‘forces’ $X_{0}$ and $X_{i}$ correspond to
the ‘fluxes’ $\dot{\eta}_{0}$ and $\dot{\eta}_{i}$, respectively, are given by

$\displaystyle
X_{k}=-\rho_{0}\frac{\partial\psi}{\partial\eta_{k}}+\nabla_{0}\cdot\left(\rho_{0}J\frac{\partial\psi^{\nabla}}{\partial\nabla_{0}\eta_{k}}\right)\quad\text{for
all }k=0,1,\ldots,N.$ (2.31)

Using Eqs. (2.8) and (2.9) in (2.31) and also applying the following
identities (for all $k=0,1,\ldots,N$)

$\nabla_{0}\eta_{k}=\boldsymbol{F}^{t}\cdot\nabla\eta_{k},$ (2.32)
$\displaystyle\frac{\partial\psi^{\nabla}}{\partial\nabla_{0}\eta_{k}}=\boldsymbol{F}^{-1}\cdot\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{k}},$
(2.33)

which can be easily proved using the indicial notations, we get the conjugate
forces $X_{0}$ and $X_{i}$ (for all $i=1,\ldots,N$) as

$\displaystyle X_{0}$ $\displaystyle=$
$\displaystyle\left(\boldsymbol{P}_{e}^{T}\cdot\boldsymbol{F}-J_{t}\psi_{e}\boldsymbol{I}\right):\boldsymbol{F}_{t}^{-1}\cdot\frac{\partial\boldsymbol{F}_{t}}{\partial\eta_{0}}-J_{t}\left.\frac{\partial\psi_{e}}{\partial\eta_{0}}\right|_{\boldsymbol{F}_{e}}-\rho_{0}(6\eta_{0}-6\eta_{0}^{2})\Delta\psi^{\theta}-J\rho_{0}\tilde{A}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\eta_{i}^{2}\eta_{j}^{2}\frac{\partial\varphi(a_{b},\eta_{0})}{\partial\eta_{0}}-$
(2.34) $\displaystyle
J\rho_{0}[A_{0M}(\theta)+(a_{\theta}-3)\Delta\psi^{\theta}(\theta)](2\eta_{0}-6\eta_{0}^{2}+4\eta_{0}^{3})-\frac{J}{2}\frac{\partial\tilde{\varphi}(a_{\beta},a_{c},\eta_{0})}{\partial\eta_{0}}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\beta_{ij}(\boldsymbol{C}^{-1}\cdot\nabla_{0}\eta_{i})\cdot\nabla_{0}\eta_{j}-$
$\displaystyle\rho_{0}\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}K_{0ij}\eta_{i}^{2}\eta_{j}^{2}+\sum_{i=1}^{N-2}\sum_{j=i+1}^{N-1}\sum_{k=j+1}^{N}K_{0ijk}\eta_{i}^{2}\eta_{j}^{2}\eta_{k}^{2}\right)\left[2(1-\varphi(a_{K},\eta_{0}))\eta_{0}-\frac{\partial\varphi(a_{K},\eta_{0})}{\partial\eta_{0}}\eta_{0}^{2}\right]+$
$\displaystyle\nabla_{0}\cdot\left(J\beta_{0M}\nabla_{0}\eta_{0}\right);$
$\displaystyle X_{i}$ $\displaystyle=$
$\displaystyle\left(\boldsymbol{P}_{e}^{T}\cdot\boldsymbol{F}-J_{t}\psi_{e}\boldsymbol{I}\right):\boldsymbol{F}_{t}^{-1}\cdot\frac{\partial\boldsymbol{F}_{t}}{\partial\eta_{i}}-J_{t}\left.\frac{\partial\psi_{e}}{\partial\eta_{i}}\right|_{\boldsymbol{F}_{e}}-2J\rho_{0}\tilde{A}\sum_{j=1,\neq
i}^{N}\eta_{i}\eta_{j}^{2}\varphi(a_{b},\eta_{0})-2\rho_{0}\sum_{j=1}^{N}K_{ij}(\eta_{i}+\eta_{j}-1)\times$
(2.35)
$\displaystyle(2\eta_{i}+\eta_{j}-1)\eta_{j}^{2}\eta_{i}-2\rho_{0}\left(\sum_{j=1}^{N}K_{0ij}\eta_{j}^{2}+\sum_{j=1}^{N-1}\sum_{k=j+1}^{N}K_{0ijk}\eta_{j}^{2}\eta_{k}^{2}\right)\eta_{0}^{2}\eta_{i}(1-\varphi(a_{K},\eta_{0}))-$
$\displaystyle
2\rho_{0}\sum_{j=1}^{N-1}\sum_{k=j+1}^{N}K_{ijk}\eta_{i}\eta_{j}^{2}\eta_{k}^{2}-2\rho_{0}\sum_{j=1}^{N-2}\sum_{k=j+1}^{N-1}\sum_{l=k+1}^{N}K_{ijkl}\eta_{i}\eta_{j}^{2}\eta_{k}^{2}\eta_{l}^{2}+$
$\displaystyle\nabla_{0}\cdot\left(\tilde{\varphi}(a_{\beta},a_{c},\eta_{0})J\sum_{j=1}^{N}\beta_{ij}\nabla_{0}\eta_{j}\right)\quad\text{for
all }i=1,2,3,\ldots,N.$

In Eqs. (2.34) and (2.35) the conjugate forces are expressed with respect to
the field variables in the reference configuration $\Omega_{0}$.

Alternatively, these forces can be expressed with respect to the field
variables in $\Omega$ as follows. Using the following identity (see e.g.
Chapter 2 of [72])

$\displaystyle\nabla_{0}\cdot(Cof\,\boldsymbol{F})=\nabla_{0}\cdot(J\boldsymbol{F}^{-t})=\boldsymbol{0},$
(2.36)

we can rewrite the conjugate force $X_{l}$ given by Eq. (2.31) in terms of the
derivative with respect to $\Omega(t)$ as

$\displaystyle
X_{k}=-J\rho\frac{\partial\psi}{\partial\eta_{k}}+J\nabla\cdot\left(\rho_{0}\frac{\partial\psi^{\nabla}}{\partial\nabla\eta_{k}}\right)=J\left[-\frac{\partial(\rho\psi)}{\partial\eta_{k}}+\nabla\cdot\left(\frac{\partial(\rho
J\psi^{\nabla})}{\partial\nabla\eta_{k}}\right)\right]\quad\text{for all
}k=0,1,\ldots,N.$ (2.37)

The conjugate forces given by Eqs. (2.34) and (2.35) can be rewritten in terms
of the Cauchy stress and the gradient operator $\nabla(\cdot)$ as

$\displaystyle X_{0}$ $\displaystyle=$
$\displaystyle\left(J\boldsymbol{F}^{-1}\cdot\boldsymbol{\sigma}_{e}\cdot\boldsymbol{F}-J_{t}\psi_{e}\boldsymbol{I}\right):\boldsymbol{F}_{t}^{-1}\cdot\frac{\partial\boldsymbol{F}_{t}}{\partial\eta_{0}}-J_{t}\left.\frac{\partial\psi_{e}}{\partial\eta_{0}}\right|_{\boldsymbol{F}_{e}}-\rho_{0}(6\eta_{0}-6\eta_{0}^{2})\Delta\psi^{\theta}-J\rho_{0}\tilde{A}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\eta_{i}^{2}\eta_{j}^{2}\frac{\partial\varphi(a_{b},\eta_{0})}{\partial\eta_{0}}-$
(2.38) $\displaystyle
J\rho_{0}[A_{0M}(\theta)+(a_{\theta}-3)\Delta\psi^{\theta}(\theta)](2\eta_{0}-6\eta_{0}^{2}+4\eta_{0}^{3})-\frac{J}{2}\frac{\partial\tilde{\varphi}(a_{\beta},a_{c},\eta_{0})}{\partial\eta_{0}}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\beta_{ij}\nabla\eta_{i}\cdot\nabla\eta_{j}-$
$\displaystyle\rho_{0}\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}K_{0ij}\eta_{i}^{2}\eta_{j}^{2}+\sum_{i=1}^{N-2}\sum_{j=i+1}^{N-1}\sum_{k=j+1}^{N}K_{0ijk}\eta_{i}^{2}\eta_{j}^{2}\eta_{k}^{2}\right)\left[2(1-\varphi(a_{K},\eta_{0}))\eta_{0}-\frac{\partial\varphi(a_{K},\eta_{0})}{\partial\eta_{0}}\eta_{0}^{2}\right]+$
$\displaystyle J\nabla\cdot\left(\beta_{0M}\nabla\eta_{0}\right);$
$\displaystyle X_{i}$ $\displaystyle=$
$\displaystyle\left(J\boldsymbol{F}^{-1}\cdot\boldsymbol{\sigma}_{e}\cdot\boldsymbol{F}-J_{t}\psi_{e}\boldsymbol{I}\right):\boldsymbol{F}_{t}^{-1}\cdot\frac{\partial\boldsymbol{F}_{t}}{\partial\eta_{i}}-J_{t}\left.\frac{\partial\psi_{e}}{\partial\eta_{i}}\right|_{\boldsymbol{F}_{e}}-2J\rho_{0}\tilde{A}\sum_{j=1,\neq
i}^{N}\eta_{i}\eta_{j}^{2}\varphi(a_{b},\eta_{0})-2\rho_{0}\sum_{j=1}^{N}K_{ij}(\eta_{i}+\eta_{j}-1)\times$
(2.39)
$\displaystyle(2\eta_{i}+\eta_{j}-1)\eta_{j}^{2}\eta_{i}-2\rho_{0}\left(\sum_{j=1}^{N}K_{0ij}\eta_{j}^{2}+\sum_{j=1}^{N-1}\sum_{k=j+1}^{N}K_{0ijk}\eta_{j}^{2}\eta_{k}^{2}\right)\eta_{0}^{2}\eta_{i}(1-\varphi(a_{K},\eta_{0}))-$
$\displaystyle
2\rho_{0}\sum_{j=1}^{N-1}\sum_{k=j+1}^{N}K_{ijk}\eta_{i}\eta_{j}^{2}\eta_{k}^{2}-2\rho_{0}\sum_{j=1}^{N-2}\sum_{k=j+1}^{N-1}\sum_{l=k+1}^{N}K_{ijkl}\eta_{i}\eta_{j}^{2}\eta_{k}^{2}\eta_{l}^{2}+$
$\displaystyle
J\nabla\cdot\left(\tilde{\varphi}(a_{\beta},a_{c},\eta_{0})\sum_{j=1}^{N}\beta_{ij}\nabla\eta_{j}\right)\quad\text{for
all }i=1,2,3,\ldots,N.$

From the dissipation inequality (2.30)1 we derive the kinetic law for
$\eta_{0}$ as

$\dot{\eta}_{0}=L_{0M}X_{0},$ (2.40)

where $L_{0M}>0$ is the kinetic coefficient for $\sf A\leftrightarrow\sf M$
PTs. In order to derive the kinetic laws for the order parameters
$\eta_{1},\ldots,\eta_{N}$ using the inequality (2.30)2, we introduce

$\dot{\eta}_{i}=\sum_{j=1}^{N}\dot{\eta}_{ij},\qquad\text{where
}\dot{\eta}_{ij}=-\dot{\eta}_{ji}\text{ and }\dot{\eta}_{ii}=0\quad\text{for
all }i,j=1,\ldots,N.$ (2.41)

Using Eq. (2.41), the dissipation rate due to the evolution of the martensitic
variants given by Eq. (2.30)2 is rewritten as

$\displaystyle D_{M}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\dot{\eta}_{i}X_{i}=\sum_{i=1}^{N}\sum_{j=1}^{N}X_{i}\dot{\eta}_{ij}\qquad\text{(using
Eq. \eqref{vari_toj}${}_{1}$)}$ (2.42) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}X_{ij}\dot{\eta}_{ij}+\sum_{i=1}^{N}\sum_{j=1}^{N}X_{j}\dot{\eta}_{ij}\qquad\text{(using
$X_{ij}=X_{i}-X_{j}$)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}X_{ij}\dot{\eta}_{ij}-\sum_{i=1}^{N}\sum_{j=1}^{N}X_{j}\dot{\eta}_{ji}\qquad\text{(using
Eq. \eqref{vari_toj}${}_{2}$)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}X_{ij}\dot{\eta}_{ij}-\sum_{i=1}^{N}\sum_{j=1}^{N}X_{i}\dot{\eta}_{ij}\qquad\text{(swapping
the indices in the second term)}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}X_{ij}\dot{\eta}_{ij}-\sum_{i=1}^{N}X_{i}\dot{\eta}_{i}.$

Noticing that the second term on the right-hand side of Eq. (2.42) is equal to
$D_{M}$ (compare with Eq. (2.30)2) we obtain

$\displaystyle D_{M}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}X_{ij}\dot{\eta}_{ij}=\sum_{j=1}^{N-1}\sum_{i=j+1}^{N}X_{ij}\dot{\eta}_{ij}\geq
0.$ (2.43)

Using the inequality (2.43) we derive the Ginzburg-Landau equations for the
evolution of the variants as

$\displaystyle\dot{\eta}_{ij}=L_{ij}(X_{i}-X_{j}),$ (2.44)

where $L_{ij}\geq 0$ is the kinetic coefficient for transformations between
$\sf M_{i}$ and $\sf M_{j}$ and is taken as in [56]

$\displaystyle L_{ij}\left\\{\begin{array}[]{@{}lr@{}}\neq 0&\text{if
}(X_{i}-X_{j})\geq
0\quad\text{and}\quad\\{0\leq\eta_{i}<1\,\,\&\,\,0<\eta_{j}\leq 1\\}\\\ \neq
0&\text{if }(X_{i}-X_{j})\leq 0\quad\text{and}\quad\\{0<\eta_{i}\leq
1\,\,\&\,\,0\leq\eta_{j}<1\\}\\\ =0&\text{if }(X_{i}-X_{j})\geq
0\quad\text{and}\quad\\{\eta_{i}=1\,\,\text{or}\,\,\eta_{j}=0\\}\\\
=0&\text{if }(X_{i}-X_{j})\leq
0\quad\text{and}\quad\\{\eta_{i}=0\,\,\text{or}\,\,\eta_{j}=1\\}.\end{array}\right.$
(2.49)

Substituting Eq. (2.44) in Eq. (2.41) the Ginzburg-Landau equations for all
$N$ order parameters $\eta_{1},\ldots,\eta_{N}$ are obtained as

$\displaystyle\dot{\eta}_{i}=\sum_{j=1,j\neq i}^{N}L_{ij}(X_{i}-X_{j}).$
(2.50)

We note that $L_{ij}$ in Eq. (2.50) and defined in Eq. (2.49) are piece-wise
constant, jumping between their finite values and zero depending on the
driving forces and the order parameters. If $L_{ij}$ is simply assumed to be
constants similar to our earlier work in [21], an issue would arise. To
understand it clearly, let us consider, without any loss of generality, a
three-variant martensitic system (where $\eta_{0}=1$) with $\sf M_{1}$, $\sf
M_{2}$, and $\sf M_{3}$. Using the constraint $\eta_{1}+\eta_{2}+\eta_{3}=1$,
the two independent Ginzburg-Landau equations from Eq. (2.50), when expressed
for $\dot{\eta}_{1}$ and $\dot{\eta}_{2}$, are given by

$\displaystyle\dot{\eta}_{1}=L_{12}(X_{1}-X_{2})+L_{13}(X_{1}-X_{3})\quad\text{and}\quad\dot{\eta}_{2}=L_{12}(X_{2}-X_{1})+L_{23}(X_{2}-X_{3}).$
(2.51)

We now consider a martensitic region in the domain where $\sf M_{3}$ is absent
and only the variants $\sf M_{1}$ and $\sf M_{2}$ are evolving within an
arbitrary time interval. The order parameters $\eta_{1}$ and $\eta_{2}$ hence
must be determined using the equation
$\dot{\eta}_{1}=\dot{\eta}_{2}=L_{12}(X_{1}-X_{2})$ as $\eta_{1}+\eta_{2}=1$
therein, and this is possible if and only if $L_{13}=L_{23}=0$ therein within
that time interval. However, if the coefficients $L_{13}$ and $L_{23}$ are
taken as nonzero constants, the contributions from the terms
$L_{13}(X_{1}-X_{3})$ and $L_{23}(X_{2}-X_{3})$ would be there unwantedly,
since the driving forces $X_{1}-X_{3}$ and $X_{2}-X_{3}$ might be nonzero
there. The desired condition can obviously be fulfilled by $L_{ij}$ given by
Eq. (2.49), but not by constant $L_{ij}$ considered in [21]. The essence of
the third and fourth conditions in Eq. (2.49) is that if variant $i$ is
absent, it cannot be transformed into other variants [56].

###  2.6 Boundary conditions

The boundary conditions for the phase-field equations and the mechanics
problem namely the Dirichlet, Neumann, and periodic BCs are listed here.

Phase-field problem. We have applied the periodic BC for all the order
parameters. If two boundaries $S_{p1\eta_{k}}\subset S_{0}$ and
$S_{p2\eta_{k}}\subset S_{0}$ (where $S_{p1\eta_{k}}\cap S_{p2\eta_{k}}$ is
empty), having opposite unit normals in $\Omega_{0}$, i.e.
$({\boldsymbol{n}_{0}})_{S_{p1\eta_{k}}}=-(\boldsymbol{n}_{0})_{S_{p2\eta_{k}}}$,
are subjected to a periodic BC related to the order parameters $\eta_{k}$ (for
$k=0,1,2,\ldots,N$), then the order parameters and their gradients must
satisfy

$\eta_{k}|_{S_{p1\eta_{k}}}=\eta_{k}|_{S_{p2\eta_{k}}}\quad\text{and}\quad(\nabla_{0}\eta_{k}\cdot\boldsymbol{n}_{0})_{S_{p1\eta_{k}}}=(\nabla_{0}\eta_{k}\cdot\boldsymbol{n}_{0})_{S_{p2\eta_{k}}}\qquad\text{for
all }k=0,1,2,\ldots,N.$ (2.52)

Mechanics problem. On the traction boundary $S_{0T}\subset S_{0}$, the
traction vector is specified (denoted by $\boldsymbol{p}^{sp}$), and on the
displacement boundary $S_{0u}\subset S_{0}$, the displacements are specified
(denoted by $\boldsymbol{u}^{sp}$), i.e.

$\boldsymbol{P}\cdot\boldsymbol{n}_{0}=\boldsymbol{p}^{sp}\qquad\text{on
}S_{0T},$ (2.53) $\boldsymbol{u}=\boldsymbol{u}^{sp}\qquad\text{on }S_{0u}.$
(2.54)

If two boundaries $S_{pu1}\subset S_{0}$ and $S_{pu2}\subset S_{0}$ (where
$S_{pu1}\cap S_{pu2}$ is empty) are subjected to a periodic BC on the
displacements $\boldsymbol{u}$, then the displacements on these boundaries are
related by

${\boldsymbol{u}}|_{S_{pu1}}={\boldsymbol{u}}|_{S_{pu2}}+({\boldsymbol{F}}_{h}-\boldsymbol{I})\cdot{\boldsymbol{r}}_{0},$
(2.55)

where ${\boldsymbol{F}}_{h}$ is a specified homogeneous deformation gradient.
The mixed boundary conditions where, on a single surface, some components of
displacements are specified and some components of the traction are specified,
are also used. Furthermore, one or more displacement components on a boundary
may be related by the periodic BC given by Eq. (2.55) with another appropriate
boundary.

##  3 Crystallographic solutions for twins within twins microstrucrure

Figure 1: A schematic of twins within twins.

In this section, we obtain the analytical solutions for the twins within twins
microstructures for the cubic to tetragonal MTs using the crystallographic
theory (see e.g. [2]). A schematic of twins within twins is shown in Fig. 1
where the twins formed by a variants pair $\sf M_{i}$ and $\sf M_{j}$, and the
twins formed by another variants pair $\sf M_{k}$ and $\sf M_{l}$ form an
interface of finite thickness $\delta_{k}$ shown by a shaded region. The
volume fractions of $\sf M_{i}$ and $\sf M_{k}$ in the respective twins are
$\kappa_{1}$ and $\kappa_{2}$. The variants $\sf M_{i}$ and $\sf M_{j}$ need
not be in twin relationship with the other two variants $\sf M_{k}$ or $\sf
M_{l}$ [2].

###  3.1 Crystallographic equations and general solutions

Crystallographic equations for twin-twin: The equations for the twins between
the pairs $\sf M_{i}$-$\sf M_{j}$, and $\sf M_{k}$-$\sf M_{l}$ are (see e.g.
Chapter 7 of [2] and [11, 12])

$\boldsymbol{R}_{2}\cdot\boldsymbol{U}_{tj}-\boldsymbol{R}_{1}\cdot\boldsymbol{U}_{ti}={\boldsymbol{a}}_{1}^{\prime}\otimes\boldsymbol{n}_{1},\quad\text{and}\quad\boldsymbol{R}_{4}\cdot\boldsymbol{U}_{tl}-\boldsymbol{R}_{3}\cdot\boldsymbol{U}_{tk}={\boldsymbol{a}}_{2}^{\prime}\otimes\boldsymbol{n}_{2},$
(3.1)

respectively, where $\boldsymbol{R}_{1}$, $\boldsymbol{R}_{2}$,
$\boldsymbol{R}_{3}$, and $\boldsymbol{R}_{4}$ are the rotation tensors;
$\boldsymbol{n}_{1}$ and $\boldsymbol{n}_{2}$ are the unit normals to the
respective twin boundaries such that $\boldsymbol{n}_{1}$ points into $\sf
M_{j}$ and $\boldsymbol{n}_{2}$ points into $\sf M_{l}$ (see Fig. 1); the
vectors ${\boldsymbol{a}}_{1}^{\prime}$ and ${\boldsymbol{a}}_{2}^{\prime}$
are related to the simple shear deformations. The governing equation for the
twins within twins shown in Fig. 1 is (Chapter 7 of [2] and [11, 12])

$(\kappa_{1}\boldsymbol{R}_{2}\cdot\boldsymbol{U}_{tj}+(1-\kappa_{1})\boldsymbol{R}_{1}\cdot\boldsymbol{U}_{ti})-(\kappa_{2}\boldsymbol{R}_{4}\cdot\boldsymbol{U}_{tl}+(1-\kappa_{2})\boldsymbol{R}_{3}\cdot\boldsymbol{U}_{tk})={\boldsymbol{b}}^{\prime}\otimes{\boldsymbol{m}},$
(3.2)

where $\kappa_{1}$ and $\kappa_{2}$ are the volume fractions of $\sf M_{i}$
and $\sf M_{k}$ in the respective twins, ${\boldsymbol{m}}$ is the unit normal
to the twin-twin interface shown in the figure, and
${\boldsymbol{b}}^{\prime}$ is a vector related to the deformation.

Equations (3.1)1 and (3.1)2 can be rewritten as

$\boldsymbol{Q}_{1}\cdot\boldsymbol{U}_{tj}-\boldsymbol{U}_{ti}=\boldsymbol{a}_{1}\otimes\boldsymbol{n}_{1},\quad\text{and}\quad\qquad\boldsymbol{Q}_{2}\cdot\boldsymbol{U}_{tl}-\boldsymbol{U}_{tk}=\boldsymbol{a}_{2}\otimes\boldsymbol{n}_{2},$
(3.3)

where $\boldsymbol{Q}_{1}=\boldsymbol{R}_{1}^{t}\cdot\boldsymbol{R}_{2},$
$\boldsymbol{a}_{1}=\boldsymbol{R}_{1}^{t}\cdot\boldsymbol{a}_{1}^{\prime}$,
$\boldsymbol{Q}_{2}=\boldsymbol{R}_{3}^{t}\cdot\boldsymbol{R}_{4}$, and
$\boldsymbol{a}_{2}=\boldsymbol{R}_{3}^{t}\cdot\boldsymbol{a}_{2}^{\prime}$.
Similarly, using Eq. (3.3)1,2, Eq. (3.2) is rewritten as

$\boldsymbol{Q}_{3}\cdot(\boldsymbol{U}_{ti}+\kappa_{1}{\boldsymbol{a}}_{1}\otimes\boldsymbol{n}_{1})-(\boldsymbol{U}_{tk}+\kappa_{2}{\boldsymbol{a}}_{2}\otimes\boldsymbol{n}_{2})={\boldsymbol{b}}\otimes\boldsymbol{m},$
(3.4)

where $\boldsymbol{Q}_{3}=\boldsymbol{R}_{3}^{t}\cdot\boldsymbol{R}_{1}$ and
$\boldsymbol{b}=\boldsymbol{R}_{3}^{t}\cdot{\boldsymbol{b}}^{\prime}$. In
order to solve Eq. (3.4), we post-multiply it with
$(\boldsymbol{U}_{tk}+\kappa_{2}{\boldsymbol{a}}_{2}\otimes\boldsymbol{n}_{2})^{-1}$
and rearrange the terms to rewrite the equation as [11]

$\boldsymbol{Q}_{3}\cdot\tilde{\boldsymbol{A}}=\boldsymbol{I}+{\boldsymbol{b}}\otimes\tilde{\boldsymbol{m}},$
(3.5)

where

$\tilde{\boldsymbol{A}}=(\boldsymbol{U}_{ti}+\kappa_{1}{\boldsymbol{a}}_{1}\otimes\boldsymbol{n}_{1})\cdot(\boldsymbol{U}_{tk}+\kappa_{2}{\boldsymbol{a}}_{2}\otimes\boldsymbol{n}_{2})^{-1},\quad\text{
and
}\quad\tilde{\boldsymbol{m}}=(\boldsymbol{U}_{tk}+\kappa_{2}{\boldsymbol{a}}_{2}\otimes\boldsymbol{n}_{2})^{-t}\cdot{\boldsymbol{m}}.$
(3.6)

The variables to be determined from the above equations (3.3) to (3.6) are
$\kappa_{1}$, $\kappa_{2}$, $\delta_{\kappa}$, $\boldsymbol{a}_{1}$,
$\boldsymbol{a}_{2}$, $\boldsymbol{n}_{1}$, $\boldsymbol{n}_{2}$,
$\boldsymbol{b}$, $\boldsymbol{m}$, $\boldsymbol{Q}_{1}$,
$\boldsymbol{Q}_{2}$, and $\boldsymbol{Q}_{3}$.

Twins within twins solution: The solution for Eq. (3.3)1 and (3.3)2 are well
known (see e.g. Chapter 5 of [2]) and is enlisted here for completeness. To do
so, we define a symmetric tensor
$\boldsymbol{G}_{1}=\boldsymbol{U}_{tj}^{-1}\cdot\boldsymbol{U}_{ti}^{2}\cdot\boldsymbol{U}_{tj}^{-1}$
corresponding to Eq. (3.3)1. The eigenvalues of $\boldsymbol{G}_{1}$ are
denoted by $\lambda_{1}$, $\lambda_{2}$, and $\lambda_{3}$, which are all
positive, and the corresponding normalized eigenvectors are denoted by
$\boldsymbol{i}_{1}$, $\boldsymbol{i}_{2}$, and $\boldsymbol{i}_{3}$,
respectively. Equation (3.3)1 has a solution if and only if $\lambda_{1}\leq
1$, $\lambda_{2}=1$, and $\lambda_{3}\geq 1$ (assuming
$\lambda_{1}\leq\lambda_{2}\leq\lambda_{3}$). The expressions for
$\boldsymbol{a}_{1}$ and $\boldsymbol{n}_{1}$ are given by

$\displaystyle\boldsymbol{a}_{1}$ $\displaystyle=$
$\displaystyle\zeta_{1}\left(\sqrt{\frac{\lambda_{3}(1-\lambda_{1})}{\lambda_{3}-\lambda_{1}}}\,\,\boldsymbol{i}_{1}+\xi\sqrt{\frac{\lambda_{1}(\lambda_{3}-1)}{\lambda_{3}-\lambda_{1}}}\,\,\boldsymbol{i}_{3}\right),$
$\displaystyle\boldsymbol{n}_{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\lambda_{3}}-\sqrt{\lambda_{1}}}{\zeta_{1}\sqrt{\lambda_{3}-\lambda_{1}}}\left(-\sqrt{{1-\lambda_{1}}}\,\,\boldsymbol{U}_{tj}\boldsymbol{i}_{1}+\xi\sqrt{\lambda_{3}-1}\,\,\boldsymbol{U}_{tj}\boldsymbol{i}_{3}\right)\,$
(3.7)

where $\xi=\pm 1$, and $\zeta_{1}$ is such that $|\boldsymbol{n}_{1}|=1$. The
solutions $\boldsymbol{a}_{2}$ and $\boldsymbol{n}_{2}$ for the twins between
$\sf M_{k}$ and $\sf M_{l}$ are similarly obtained using the eigenpairs of
$\boldsymbol{U}_{tl}^{-1}\cdot\boldsymbol{U}_{tk}^{2}\cdot\boldsymbol{U}_{tl}^{-1}$
in Eq. (3.7). The rotations $\boldsymbol{Q}_{1}$ and $\boldsymbol{Q}_{2}$ can
then be obtained using Eqs. (3.3)1,2.

We now solve the twins within twins equation (3.5) which has got a form
similar to the $\sf A$-twinned martensite interface equation using the
procedure of Ball and James [5] (also see Chapter 7 of [11]). Noticing that
the Bain stretches $\boldsymbol{U}_{ti}$, $\boldsymbol{U}_{tj}$,
$\boldsymbol{U}_{tk}$, and $\boldsymbol{U}_{tl}$ are given for a material, we
first obtain $\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$, $\boldsymbol{n}_{1}$,
and $\boldsymbol{n}_{2}$ using, e.g., Eq. (3.7). The procedure for obtaining
the remaining unknowns $\kappa_{1}$, $\kappa_{2}$, $\boldsymbol{Q}_{3}$,
$\boldsymbol{b}$, and $\tilde{\boldsymbol{m}}$ is derived here. In obtaining
these unknowns, let us first assume that the parameters $\kappa_{1}$ and
$\kappa_{2}$ are given, and solve for the remaining unknown variables. We
introduce

$\boldsymbol{G}_{2}=\tilde{\boldsymbol{A}}^{t}\cdot\tilde{\boldsymbol{A}},$
(3.8)

which is obviously symmetric and positive-definite, and all the eigenvalues
are positive numbers which we denote by $\Lambda_{1}$, $\Lambda_{2}$, and
$\Lambda_{3}$. The corresponding normalized eigenvectors are denoted by
$\boldsymbol{j}_{1}$, $\boldsymbol{j}_{2}$, and $\boldsymbol{j}_{3}$,
respectively. Equation (3.5) has a solution if and only if $\Lambda_{1}\leq
1$, $\Lambda_{2}=1$, and $\Lambda_{3}\geq 1$ assuming
$\Lambda_{1}\leq\Lambda_{2}\leq\Lambda_{3}$. The solutions for the vectors
$\boldsymbol{b}$ and $\tilde{\boldsymbol{m}}$ are obtained as (see e.g.
Chapter 6 of [2])

$\displaystyle\boldsymbol{b}$ $\displaystyle=$
$\displaystyle\frac{\zeta_{2}}{\sqrt{\Lambda_{3}-\Lambda_{1}}}\left(\sqrt{\Lambda_{3}(1-\Lambda_{1})}\,\,\boldsymbol{j}_{1}+\xi\sqrt{\Lambda_{1}(\Lambda_{3}-1)}\,\,\boldsymbol{j}_{3}\right),\quad\text{and}$
$\displaystyle\tilde{\boldsymbol{m}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\Lambda_{3}}-\sqrt{\Lambda_{1}}}{\zeta_{2}\sqrt{\Lambda_{3}-\Lambda_{1}}}\left(-\sqrt{{1-\Lambda_{1}}}\,\,\boldsymbol{j}_{1}+\xi\sqrt{\Lambda_{3}-1}\,\,\boldsymbol{j}_{3}\right).$
(3.9)

The unit normal to the twin-twin boundary $\boldsymbol{m}$ is finally obtained
using Eq. (3.9)2 in Eq. (3.6)2

$\displaystyle{\boldsymbol{m}}=\frac{\sqrt{\Lambda_{3}}-\sqrt{\Lambda_{1}}}{\zeta_{2}\sqrt{\Lambda_{3}-\Lambda_{1}}}(\boldsymbol{U}_{tk}+\kappa_{2}\boldsymbol{n}_{2}\otimes{\boldsymbol{a}}_{2})\cdot\left(-\sqrt{{1-\Lambda_{1}}}\,\,\boldsymbol{j}_{1}+\xi\sqrt{\Lambda_{3}-1}\,\,\boldsymbol{j}_{3}\right),$
(3.10)

where $\zeta_{2}$ in Eqs. (3.9) and (3.10) is such that $|\boldsymbol{m}|=1$.
The rotation $\boldsymbol{Q}_{3}$ is then determined using Eqs. (3.9)1,2 in
Eq. (3.5). Note that the middle eigenvalues $\Lambda_{2}$ obtained would be an
expression as a function of the volume fractions $\kappa_{1}$ and
$\kappa_{2}$. Setting $\Lambda_{2}=1$, which is required for the existence of
the twins within twins solution [11], would give a relation between
$\kappa_{1}$ and $\kappa_{2}$. However, it is not possible to obtain the exact
solutions for $\kappa_{1}$ and $\kappa_{2}$ from the limited governing
equations at hand. The thickness of the transition layer $\kappa_{t}$ is also
indeterminate.

###  3.2 Twins within twins solutions for cubic to tetragonal MTs

The solutions for twins within twins for cubic to tetragonal MTs are now
obtained. The three Bain stretch tensors for such transformations are [2]

$\displaystyle\boldsymbol{U}_{t1}$ $\displaystyle=$
$\displaystyle\chi\,\boldsymbol{c}_{1}\otimes\boldsymbol{c}_{1}+\alpha\,\boldsymbol{c}_{2}\otimes\boldsymbol{c}_{2}+\alpha\,\boldsymbol{c}_{3}\otimes\boldsymbol{c}_{3},$
$\displaystyle\boldsymbol{U}_{t2}$ $\displaystyle=$
$\displaystyle\alpha\,\boldsymbol{c}_{1}\otimes\boldsymbol{c}_{1}+\chi\,\boldsymbol{c}_{2}\otimes\boldsymbol{c}_{2}+\alpha\,\boldsymbol{c}_{3}\otimes\boldsymbol{c}_{3},$
$\displaystyle\boldsymbol{U}_{t3}$ $\displaystyle=$
$\displaystyle\alpha\,\boldsymbol{c}_{1}\otimes\boldsymbol{c}_{1}+\alpha\,\boldsymbol{c}_{2}\otimes\boldsymbol{c}_{2}+\chi\,\,\boldsymbol{c}_{3}\otimes\boldsymbol{c}_{3},$
(3.11)

where $\alpha<1$ and $\chi>1$ are the material constants and
$\\{\boldsymbol{c}_{1},\,\boldsymbol{c}_{2},\,\boldsymbol{c}_{3}\\}$ is a
right-handed standard Cartesian basis in the cubic unit cell of $\sf A$ such
that the basis vectors are parallel to three mutually orthogonal sides of that
unit cell. The solutions for twins between $\sf M_{1}$-$\sf M_{2}$ are
($\boldsymbol{n}$ pointing into $\sf M_{2}$), $\sf M_{1}$-$\sf M_{3}$
($\boldsymbol{n}$ pointing into $\sf M_{3}$), and $\sf M_{2}$-$\sf M_{3}$
($\boldsymbol{n}$ pointing into $\sf M_{3}$) are (see e.g. Chapter 5 of [2])

$\displaystyle\boldsymbol{a}$ $\displaystyle=$
$\displaystyle\sqrt{2}\upsilon(-\chi\boldsymbol{c}_{1}\pm\alpha\boldsymbol{c}_{2}),\qquad\boldsymbol{n}=(\boldsymbol{c}_{1}\pm\boldsymbol{c}_{2})/\sqrt{2};$
$\displaystyle\boldsymbol{a}$ $\displaystyle=$
$\displaystyle\sqrt{2}\upsilon(-\chi\boldsymbol{c}_{1}\pm\alpha\boldsymbol{c}_{3}),\qquad\boldsymbol{n}=(\boldsymbol{c}_{1}\pm\boldsymbol{c}_{3})/\sqrt{2};$
$\displaystyle\boldsymbol{a}$ $\displaystyle=$
$\displaystyle\sqrt{2}\upsilon(-\chi\boldsymbol{c}_{2}\pm\alpha\boldsymbol{c}_{3}),\qquad\boldsymbol{n}=(\boldsymbol{c}_{2}\pm\boldsymbol{c}_{3})/\sqrt{2},$
(3.12)

respectively, where
$\upsilon={(\chi^{2}-\alpha^{2})}/({\chi^{2}+\alpha^{2}})$. Since all the
variants for tetragonal $\sf M$ phase are in a twin relationship [2], the
combinations of possible twins within twins solutions are $\\{\sf M_{1},\sf
M_{2}\\}$-$\\{\sf M_{1},\sf M_{3}\\}$, $\\{\sf M_{2},\sf M_{1}\\}$-$\\{\sf
M_{2},\sf M_{3}\\}$, and $\\{\sf M_{2},\sf M_{3}\\}$-$\\{\sf M_{2},\sf
M_{1}\\}$, where the corresponding twin solutions
$\\{\boldsymbol{a}_{1},\boldsymbol{n}_{1}\\}$ and
$\\{\boldsymbol{a}_{2},\boldsymbol{n}_{2}\\}$ are to be considered from Eq.
(3.12).

Using Eq. (3.6)1 and the solutions to the corresponding twin pairs from Eq.
(3.12), we get ${\boldsymbol{G}}_{2}$ as the following diagonal tensor for all
the possible combinations of twins within twins listed above:

${\boldsymbol{G}}_{2}=\tilde{\boldsymbol{A}}^{t}\cdot\tilde{\boldsymbol{A}}=\Lambda_{1}\boldsymbol{j}_{1}\otimes\boldsymbol{j}_{1}+\Lambda_{2}\boldsymbol{j}_{2}\otimes\boldsymbol{j}_{2}+\Lambda_{3}\boldsymbol{j}_{3}\otimes\boldsymbol{j}_{3},\quad\text{where}$
(3.13) $\displaystyle\Lambda_{1}$ $\displaystyle=$
$\displaystyle\left[1-\kappa_{2}\upsilon\right]^{2}<1,$
$\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle[(1+\kappa_{1})\alpha^{2}+(1-\kappa_{1})\chi^{2}]^{2}[(1+\kappa_{2})\chi^{2}+(1-\kappa_{2})\alpha^{2}]^{2}/(\chi^{2}+\alpha^{2})^{4}=1,\quad\text{and}\quad$
$\displaystyle\Lambda_{3}$ $\displaystyle=$
$\displaystyle\left[1+\kappa_{1}\upsilon\right]^{2}>1,$ (3.14)

and we have used the facts $\chi>1$, $\alpha<1$, $0<\upsilon<1$, and
$0<\kappa_{1},\kappa_{2}<1$. The eigenvectors $\boldsymbol{j}_{1}$,
$\boldsymbol{j}_{2}$, and $\boldsymbol{j}_{3}$ are functions of the vectors
$\boldsymbol{c}_{1}$, $\boldsymbol{c}_{2}$, and $\boldsymbol{c}_{3}$ and
depend on the combinations of the variants in the twins as obtained below. We
have imposed the condition $\Lambda_{2}=1$ in Eq. (3.14)2, which is a
requirement for $\boldsymbol{G}_{2}$ given by Eq. (3.13) to represent a twins
within twins as discussed above, and from that the following two conditions on
the parameters $\chi,\,\alpha,\,\kappa_{1},\,\kappa_{2}$ are obtained:

$\displaystyle\frac{1}{\kappa_{1}}-\frac{1}{\kappa_{2}}$ $\displaystyle=$
$\displaystyle\upsilon,\quad\text{or}$
$\displaystyle(\kappa_{1}-\kappa_{2})\frac{1}{\upsilon}+\kappa_{1}\kappa_{2}$
$\displaystyle=$ $\displaystyle\frac{2}{\upsilon^{2}}.$ (3.15)

Since $0<\upsilon<1$, Eq. (3.15)1 is satisfied if and only if
$\kappa_{1}<\kappa_{2}$. For $\kappa_{1}=\kappa_{2}$ Eq. (3.15)1 yields the
trivial condition $\chi=\alpha$ which does not yield twins within twins. It is
easy to verify that for all $\chi>1$ and $\alpha<1$ no
$0<\kappa_{1},\kappa_{2}<1$ satisfy Eq. (3.15)2, and hence we disregard this
relation. Finally, considering $\kappa_{1}$ and $\kappa_{2}$ are related by
Eq. (3.15)1, and using the expressions for $\Lambda_{1}$ and $\Lambda_{3}$
given by Eqs. (3.14)1,3 into Eqs. (3.7) and (3.10), we get the solutions for
$\boldsymbol{b}$ and $\boldsymbol{m}$ for different twins within twins as
listed below.

##### Case I

For $\\{\sf M_{1},\sf M_{2}\\}$-$\\{\sf M_{1},\sf M_{3}\\}$ twin pairs

For $\\{\sf M_{1},\sf M_{2}\\}$-$\\{\sf M_{1},\sf M_{3}\\}$ twin pairs the
indices are $i=k=1$, $j=2$, and $l=3$, and $\kappa_{1}$ and $\kappa_{2}$ are
the volume fractions of $\sf M_{1}$ in the respective twins. The eigenvectors
of $\boldsymbol{G}_{2}$ tensor are obtained as
$\boldsymbol{j}_{1}=\boldsymbol{c}_{3}$,
$\boldsymbol{j}_{2}=\boldsymbol{c}_{1}$, and
$\boldsymbol{j}_{3}=\boldsymbol{c}_{2}$. The vectors $\boldsymbol{b}$ and
$\tilde{\boldsymbol{m}}$ are obtained using Eqs. (3.14)1,3 in Eq. (3.9) as

$\displaystyle\boldsymbol{b}$ $\displaystyle=$
$\displaystyle\frac{\zeta_{2}}{\sqrt{(\kappa_{1}+\kappa_{2})[2+(\kappa_{1}-\kappa_{2})\upsilon]}}\left[\xi(1-\kappa_{2}\upsilon)\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\boldsymbol{c}_{2}+(1+\kappa_{1}\upsilon)\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\boldsymbol{c}_{3}\right],\quad\text{and}$
$\displaystyle\tilde{\boldsymbol{m}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{2}-\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\,\boldsymbol{c}_{3}\right],$
(3.16)

respectively. Using Eq. (3.16)2 and (3.7) in Eq. (3.10) we finally get
$\boldsymbol{m}$ as

$\displaystyle\boldsymbol{m}=\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon\alpha}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[-\upsilon\kappa_{2}^{1.5}\sqrt{2-\kappa_{2}\upsilon}\,\,\boldsymbol{c}_{1}+\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{2}+\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}(1+\kappa_{2}\upsilon)\,\,\boldsymbol{c}_{3}\right],$
(3.17)

where the condition $|\boldsymbol{m}|=1$ yields

$\displaystyle\zeta_{2}=\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon\alpha}{\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[\kappa_{2}(2-\kappa_{2}\upsilon)(1+2\kappa_{2}\upsilon+2\kappa_{2}^{2}\upsilon^{2})+\kappa_{1}(2+\kappa_{1}\upsilon)\right].$
(3.18)

##### Case II

For $\\{\sf M_{2},\sf M_{1}\\}$-$\\{\sf M_{2},\sf M_{3}\\}$ twin pairs

In this case the indices are $i=k=2$, $j=1$, and $l=3$, and $\kappa_{1}$ and
$\kappa_{2}$ are the volume fractions of $\sf M_{2}$ in the respective twins.
The eigenvectors for $\boldsymbol{G}_{2}$ tensor are given by
$\boldsymbol{j}_{1}=\boldsymbol{c}_{3}$,
$\boldsymbol{j}_{2}=\boldsymbol{c}_{2}$, and
$\boldsymbol{j}_{3}=\boldsymbol{c}_{1}$. The vectors $\boldsymbol{b}$ and
$\tilde{\boldsymbol{m}}$ are obtained using Eqs. (3.7) and (3.9) as

$\displaystyle\boldsymbol{b}$ $\displaystyle=$
$\displaystyle\frac{\zeta_{2}}{\sqrt{(\kappa_{1}+\kappa_{2})[2+(\kappa_{1}-\kappa_{2})\upsilon]}}\left[(1+\kappa_{1}\upsilon)\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\boldsymbol{c}_{3}+\xi(1-\kappa_{2}\upsilon)\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\boldsymbol{c}_{1}\right],\quad\text{and}$
$\displaystyle\tilde{\boldsymbol{m}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[-\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\,\boldsymbol{c}_{3}+\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{1}\right],$
(3.19)

respectively. Using Eq. (3.16)2 in Eq. (3.10) we finally get $\boldsymbol{m}$
as

$\displaystyle\boldsymbol{m}=\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon\alpha}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[-\xi\upsilon\kappa_{2}^{1.5}\sqrt{2-\kappa_{2}\upsilon}\,\,\boldsymbol{c}_{2}+\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{1}+\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}(1+\kappa_{2}\upsilon)\,\,\boldsymbol{c}_{3}\right],$
(3.20)

where $\zeta_{2}$ is given by Eq. (3.18).

##### Case III

For $\\{\sf M_{3},\sf M_{1}\\}$-$\\{\sf M_{3},\sf M_{2}\\}$ twin pairs

In this case the indices are $i=k=3$, $j=1$, and $l=2$, and $\kappa_{1}$ and
$\kappa_{2}$ are the volume fractions of $\sf M_{3}$ in the respective twins.
The eigenvectors for $\boldsymbol{G}_{2}$ tensor are given by
$\boldsymbol{j}_{1}=\boldsymbol{c}_{2}$,
$\boldsymbol{j}_{2}=\boldsymbol{c}_{3}$, and
$\boldsymbol{j}_{3}=\boldsymbol{c}_{1}$. The vectors $\boldsymbol{b}$ and
$\tilde{\boldsymbol{m}}$ are obtained using Eqs. (3.7) and (3.10) as

$\displaystyle\boldsymbol{b}$ $\displaystyle=$
$\displaystyle\frac{\zeta_{2}}{\sqrt{(\kappa_{1}+\kappa_{2})[2+(\kappa_{1}-\kappa_{2})\upsilon]}}\left[(1+\kappa_{1}\upsilon)\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\boldsymbol{c}_{2}+\xi(1-\kappa_{2}\upsilon)\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\boldsymbol{c}_{1}\right],\quad\text{and}$
$\displaystyle\tilde{\boldsymbol{m}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[-\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}\,\,\boldsymbol{c}_{2}+\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{1}\right],$
(3.21)

respectively. Using Eq. (3.16)2 in Eq. (3.10) we finally get $\boldsymbol{m}$
as

$\displaystyle\boldsymbol{m}=\frac{\sqrt{\kappa_{1}+\kappa_{2}}\,\,\upsilon\alpha}{\zeta_{2}\sqrt{2+(\kappa_{1}-\kappa_{2})\upsilon}}\left[-\xi\upsilon\kappa_{2}^{1.5}\sqrt{2-\kappa_{2}\upsilon}\,\,\boldsymbol{c}_{3}+\xi\sqrt{\kappa_{1}(2+\kappa_{1}\upsilon)}\,\,\boldsymbol{c}_{1}+\sqrt{\kappa_{2}(2-\kappa_{2}\upsilon)}(1+\kappa_{2}\upsilon)\,\,\boldsymbol{c}_{2}\right],$
(3.22)

where $\zeta_{2}$ is given by Eq. (3.18).

In summary, we have obtained the general analytical solution for
$\boldsymbol{a}_{1}$, $\boldsymbol{a}_{2}$, $\boldsymbol{n}_{1}$,
$\boldsymbol{n}_{2}$, $\boldsymbol{b}$, and $\boldsymbol{m}$ listed in Eqs.
(3.7), (3.9), and (3.10). The rotation tensors $\boldsymbol{Q}_{1}$,
$\boldsymbol{Q}_{2}$, and $\boldsymbol{Q}_{3}$ can be finally obtained using
Eqs. (3.4)1,2 and (3.3). The twins within twins solutions for cubic to
tetragonal MTs are listed in Eqs. (3.12) and (3.16) to (3.22). Since the
volume fractions $\kappa_{1}$ and $\kappa_{2}$ satisfy the relation given by
Eq. (3.15)1, there are obviously many solutions possible for each of the twins
within twins listed in Cases (I), (II) and (III) (also see Chapter 7 of [2]).
The width of the twins within twins interface (shaded region in Fig. 1)
$\delta_{\kappa}$ is indeterminate within the set of governing equations at
hand.

##  4 Results and discussions

We now present the simulation results for twins within twins microstructure
obtained using our phase-field approach. The materials properties for NiAl
alloy, which exhibits cubic to tetragonal MTs, are considered and listed in
Section 4.1. In Section 4.2 the phase-field results are compared with the
crystallographic solution obtained in Section 3.

###  4.1  Material parameters

Figure 2: Unit cells for cubic austenite and three tetragonal martensitic
variants.

The material parameters for NiAl alloy are enlisted here. We consider the
interfacial widths and energies as $\delta_{0M}=1$ nm, $\gamma_{0M}=0.2$ N/m,
$\delta_{12}=\delta_{13}=\delta_{23}=0.75$ nm, and
$\gamma_{12}=\gamma_{13}=\gamma_{23}=0.1$ N/m. Using the following analytical
relations between the interfacial thickness and energy, and the phase-field
parameters [68, 73]

$\delta_{ij}=\sqrt{18\beta_{ij}/A_{ij}};\qquad\beta_{ij}=\gamma_{ij}\delta_{ij}\quad\text{for
}\sf A\text{-}\sf M\text{ and all the variant-variant interfaces},$ (4.1)

which were obtained by solving an 1D Ginzburg-Landau equation neglecting
mechanics, we obtain $\rho_{0}A_{0M}=3600$ MPa, $\rho_{0}\bar{A}=2400$ MPa,
$\beta_{0M}=2\times 10^{-10}$ N, and
$\beta_{12}=\beta_{13}=\beta_{23}=7.5\times 10^{-11}$ N. We take
$\theta_{e}=215$ K and $\Delta s=-1.47$ MPa K-1, using which we calculate the
critical temperatures for ${\sf A}\to{\sf M}$ and ${\sf M}\to{\sf A}$
transformations as (see [60]) $\theta^{c}_{{\sf A}\to{\sf
M}}=\theta_{e}+A_{0M}/(3\Delta s)=0$ K and $\theta^{c}_{{\sf M}\to{\sf
A}}=\theta_{e}-A_{0M}/(3\Delta s)=430$ K, respectively. The Lamé constants
assuming isotropic elasticity are taken to be identical for all the phases
$\sf A$, ${\sf M}_{1}$, ${\sf M}_{2}$, and ${\sf M}_{3}$:
$\lambda_{0}=\lambda_{1}=\lambda_{2}=\lambda_{3}=74.62$ GPa,
$\mu_{0}=\mu_{1}=\mu_{2}=\mu_{3}=72$ GPa . The other material constants are
taken as $a_{\beta}=a_{\varepsilon}=a_{K}=3$, $a_{0}=10^{-3}$,
$K_{12}=K_{23}=K_{13}=50$ GPa, $K_{012}=K_{023}=K_{013}=5$ GPa,
$K_{0123}=K_{123}=50$ GPa, $L_{0M}=2600$ (Pa-s)-1 and
$L_{12}=L_{13}=L_{23}=12600$ (Pa-s)-1. The transformation stretches are
$\alpha=0.922$ and $\chi=1.215$ [46].

###  4.2 Numerical result for twins within twins

Figure 3: Evolution of twin within twin microstructure in a $22$ nm$\times 22$
nm$\times 22$ nm cube shown by color plot of
$\eta_{eq}=\eta_{0}(1-0.67\eta_{1}-0.33\eta_{2})$: $\eta_{eq}=0.33$ denotes
$\sf M_{1}$; $\eta_{eq}=0.67$ denotes $\sf M_{2}$; $\eta_{eq}=1$ denotes $\sf
M_{3}$.

For the simulation we consider a $22$ nm$\times 22$ nm$\times 22$ nm cube as
the reference configuration $\Omega_{0}$ as shown in Figs. 3(a). The periodic
BCs given by Eq. (2.52) are used for all the order parameters on the
respective opposite faces of the cube domain. We have used the periodic BC for
the normal component of the displacement vector on the opposite faces of the
cube given by Eq. (2.55) with the homogeneous deformation gradient as
$\boldsymbol{F}_{h}=\boldsymbol{I}+0.98\,\boldsymbol{n}_{0}\otimes\boldsymbol{n}_{0}$,
where $\boldsymbol{n}_{0}$ is the unit normal to the opposite faces of the
cube in $\Omega_{0}$. On each face of the cube domain, we have used the
traction-free BC for the tangential components of the first Piola-Kirchhoff
traction vector. The temperature of the sample is taken as $\theta=0$ K. The
material properties listed in Section 4.1 are used. The Bain tensors listed in
Eq. (3.11) are used in Eq. (2.5), where the basis vectors
$\boldsymbol{c}_{1}$, $\boldsymbol{c}_{2}$, and $\boldsymbol{c}_{3}$ of $\sf
A$ unit cell are parallel to basis vectors $\boldsymbol{e}_{1}$,
$\boldsymbol{e}_{2}$, and $\boldsymbol{e}_{3}$, respectively, attached to the
sample (see Fig. 3 (a)). The initial distribution of the order parameters are
taken between $0\leq\eta_{0}\leq 1$, $0\leq\eta_{1}\leq 0.8$, and
$0\leq\eta_{2}\leq 0.8$, all distributed randomly shown in Figs. 3(a). The
other order parameter $\eta_{3}$ is calculated using Eq. (2.1) for all the
times $t\geq 0$. We have used a finite element procedure similar to that of
[74] developed by the authors. A finite element code has been developed using
an open-source code deal.ii [75]. The domain is discretized spatially with
quadratic brick element and it is ensured that at least three grid points lie
across all the interfaces. The mesh density in the 3D domain is shown in Fig.
4(a), and the mesh density on one of the boundaries is shown in Fig. 4(b). The
time derivatives of the order parameters are discretized using the backward
difference scheme of order two described in [74]. A constant time step size of
$2\times 10^{-16}$ s is used for the simulation.

Figure 4: (a) Mesh density in the 3D computational domain. (b) Mesh density on
one of the external boundaries.

The evolution of the microstructure is shown at different time instances in
Figs. 3. We have in fact shown a colour plot for an equivalent order parameter
defined as $\eta_{eq}=\eta_{0}(1-0.67\eta_{1}-0.33\eta_{2})$ which obviously
takes the following values in different phases: $\eta_{eq}=0$ in $\sf A$,
$\eta_{eq}=0.33$ in $\sf M_{1}$, $\eta_{eq}=0.67$ in $\sf M_{2}$, and
$\eta_{eq}=1$ in $\sf M_{3}$. Figure 3(a) shows the initial distribution of
$\eta_{eq}$. Figures 3(b) and (c) show the intermediate microstructures at
different time instances which are approaching a twinned microstructure. We
finally obtain twins within twins microstructures between the twin pairs $\sf
M_{1}$-$\sf M_{2}$ and $\sf M_{1}$-$\sf M_{3}$ as shown in Fig. 3(d). The
microstructure shown in Fig. 3(d) is a little far from being a stationary one.
However, the twinned microstructure obtained can be compared with the
analytical solution as there is no further significant change in the
orientations of the twin boundaries and twin-twin boundaries with time as
observed in Figs. 3 (c) and (d). The microstructures on three faces of the
domain with unit normals parallel to $\boldsymbol{e}_{3}$,
$\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{1}$ are shown in Fig. 5(a), (b), and
(c), respectively.

Figure 5: The twin microstructures on three mutually perpendicular faces of
the sample shown in Fig. 3. The colour plots of $\eta_{eq}$ are shown:
$\eta_{eq}=0.33$ denotes $\sf M_{1}$; $\eta_{eq}=0.67$ denotes $\sf M_{2}$;
$\eta_{eq}=1$ denotes $\sf M_{3}$.

The plots for the components of the Cauchy elastic stress tensor (in GPa) at
$t=2.25$ ps (corresponding to the microstructure shown in Fig. 3 (d)) are
shown in Fig. 6. The internal stresses are concentrated mainly across the
twin-twin boundaries and also across the twin boundaries (see e.g. [76] for
experimental results). The stresses vary from compressive to tensile between
two adjacent variant plates within and near the twin-twin interfaces. For a
better understanding of the elastic stresses across the interfaces, we have
shown a plot for the components across the lines joined by the points O and A,
and C and B (see Fig. 3 (d)) in Fig. 7(a) and (b), respectively. Obviously, in
these two figures, the elastic stresses within $\sf M_{1}$-$\sf M_{3}$ and
$\sf M_{1}$-$\sf M_{2}$ twin boundaries are plotted. Figure 7(b) also shows
the stresses across a twin-twin boundary. All the normal stresses
$\sigma_{(e)11}$, $\sigma_{(e)22}$, and $\sigma_{(e)33}$ are significantly
higher across the twin boundaries compared to the adjacent phases. The reason
for such large elastic stresses within the twin boundaries is studied in
detail by the authors in [77]. The shear stresses $\sigma_{(e)12}$ and
$\sigma_{(e)13}$ on the corresponding external boundary (having unit normal
$\boldsymbol{e}_{1}$ in $\Omega_{0}$) are much lower due to the traction-free
BC applied in the tangential plane. From the elastic stress plots shown in
Fig. 6 it is clear that stresses are mainly concentrated within the twin
boundaries and the twins-twins boundaries. In fact, the stresses within twins
within twins boundaries are usually much higher than that across the twin
boundaries (also see Fig. 7). This can be justified by noticing that the twin-
twin boundaries are compatible in an average sense, whereas the twin
boundaries are compatible in Hadamard’s sense according to the
crystallographic theory (see Chapter 5 of [2]). Understanding the stress
distributions across these interfaces is important from the materials design
perspective [76].

Figure 6: Components of the Cauchy elastic stress tensor (in GPa) in the
domain shown in Fig. 3(d).

Figure 7:  Components of the Cauchy elastic stress tensor across the lines
drawn between points (a) O and A, and (b) C and B shown in Fig. 3(d).
Obviously, figures (a) and (b) show the elastic stress distribution across the
twin boundaries between variants $\sf M_{1}$-$\sf M_{3}$ and $\sf M_{1}$-$\sf
M_{2}$, respectively.

Table-1: Crystallographic solutions for twins between the variants for NiAl
alloy in $\\{\boldsymbol{c}_{1}\,,\boldsymbol{c}_{2}\,,\boldsymbol{c}_{3}\\}$
basis.  Variant pair $\boldsymbol{a}$ $\boldsymbol{n}$ $\sf M_{1}$-$\sf M_{2}$
$-0.4625\,\boldsymbol{c}_{1}\pm 0.3510\,\boldsymbol{c}_{2}$
$0.7071\,\boldsymbol{c}_{1}\pm 0.7071\,\boldsymbol{c}_{2}$ $\sf M_{2}$-$\sf
M_{3}$ $0.4625\,\boldsymbol{c}_{1}\pm 0.3510\,\boldsymbol{c}_{3}$
$-0.7071\,\boldsymbol{c}_{1}\pm 0.7071\,\boldsymbol{c}_{3}$ $\sf M_{3}$-$\sf
M_{1}$ $0.4625\,\boldsymbol{c}_{2}\pm 0.3510\,\boldsymbol{c}_{3}$
$-0.7071\,\boldsymbol{c}_{2}\pm 0.7071\,\boldsymbol{c}_{3}$

Figure 8: Plot for $\kappa_{2}$ versus $\kappa_{1}$ given by Eq. (3.15)1 for
NiAl alloy. Figure 9: A plane within a twins within twins interface (within
the domain shown in Fig. 3(d)) is shown by the rectangle. The unit normal to
that plane is
$\boldsymbol{m}=0.1020\,\boldsymbol{c}_{1}-0.7204\,\boldsymbol{c}_{2}+0.6859\,\boldsymbol{c}_{3}$.

##### Comparison of crystallographic and numerical solutions.

We now present a comparative study of the microstructure obtained numerically
with the crystallographic solution given in Section 3.2. We have shown in Fig.
3 the twins within twins microstructures obtained using the present phase-
field approach, where a twin with a variant pair ${\sf M}_{1}$-${\sf M}_{2}$
is forming interfaces with another twin made of variants pair ${\sf
M}_{1}$-${\sf M}_{3}$, i.e. the indices of Fig. 1 are $i=k=1$, $j=2$, and
$l=3$. The normal to the interfaces between ${\sf M}_{1}$ and ${\sf M}_{2}$
plates are making approximately $45^{\circ}$ with both $\boldsymbol{e}_{1}$
and $\boldsymbol{e}_{2}$ axes, and the normal to the interfaces between ${\sf
M}_{1}$ and ${\sf M}_{3}$ plates are making approximately $45^{\circ}$ with
both $\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{3}$ axes. The normals to these
twins are obviously in good agreement with the crystallographic solution
listed in Table-1. The volume fractions of ${\sf M}_{1}$ in the respective
twin pairs are calculated as $\kappa_{1}=0.46$ and $\kappa_{2}=0.56$. The unit
normal to one of the twins within twins interface $\boldsymbol{m}$ is obtained
as

$\boldsymbol{m}=0.1020\,\boldsymbol{c}_{1}-0.7204\,\boldsymbol{c}_{2}+0.6859\,\boldsymbol{c}_{3},$
(4.2)

where a rectangular plane (on red lines) lying on the twin-twin interface is
shown in Fig. 9. We now calculate the analytical expression for
$\boldsymbol{m}$. As mentioned earlier, the volume fractions $\kappa_{1}$ and
$\kappa_{2}$, satisfying the relation given by Eq. (3.15)1 and plotted in Fig.
8, cannot be uniquely determined from the limited governing equations at hand.
In view of that, we simply assume $\kappa_{1}=0.46$ for the analytical
solution and obtain $\kappa_{2}=0.5250$ using Eq. (3.15)1 which differs from
the numerical result by $7.7$%. Finally, using Eqs. (3.16), (3.17), and (3.18)
we obtain $\zeta_{2}=0.2634$ and the vectors $\boldsymbol{b}$,
$\tilde{\boldsymbol{m}}$, and $\boldsymbol{m}$ as

$\displaystyle\boldsymbol{b}$ $\displaystyle=$ $\displaystyle\pm
0.1362\,\boldsymbol{c}_{2}+0.2093\,\boldsymbol{c}_{3},$
$\displaystyle\tilde{\boldsymbol{m}}$ $\displaystyle=$ $\displaystyle\pm
0.7119\,\boldsymbol{c}_{2}-0.7115\,\boldsymbol{c}_{3},\quad\text{and}$
$\displaystyle\boldsymbol{m}$ $\displaystyle=$
$\displaystyle-0.0928\,\boldsymbol{c}_{1}\pm
0.6560\,\boldsymbol{c}_{2}+0.7492\,\boldsymbol{c}_{3}.$ (4.3)

The maximum difference in the components of analytical and numerically
obtained unit normals $\boldsymbol{m}$ is within 10%, and the orientation of
the analytical $\boldsymbol{m}$ differs by $12.3^{\circ}$ from the numerical
one, which is acceptable. One of the main sources of difference is that in the
numerical solution local stress fields and their relaxation by incomplete
martensitic variants at the twin within twin interface is automatically taken
into account.

##  5 Concluding remarks

A thermodynamically consistent Ginzburg-Landau type nanoscale multiphase
phase-field approach to multivariant martensitic transformations is formulated
at large strain and taking the interfacial stresses into consideration. A
total of $N$ independent order parameters are assumed to describe the
austenite and $N$ martensitic variants. A non-contradictory form for the
gradient energy is considered, and the gradient energies with the previous
models are compared. Furthermore, a non-contradictory kinetic model for the
coupled Ginzburg-Landau equations is derived for all the order parameters. The
kinetic models from the previous models are compared and their shortcomings
are discussed. A general crystallographic solution for the twins within twins
microstructure is obtained, and the solutions for cubic to tetragonal
transformations are presented. Using the present phase-field approach 3D twins
within twins microstructures evolution in a single grain is studied, and the
numerical results are compared with the crystallographic solution. The present
phase-field model can be used for studying more complex microstructures with
more variants such as cubic$\leftrightarrow$orthorhombic and
cubic$\leftrightarrow$monoclinic MTs, and also in polycrystalline solids. Note
that for large transformation strains for the MTs Si I to Si II, a new
martensitic microstructure, which does not obey the mathematical theory of
martensite [4, 2, 3, 5, 6], was obtained with molecular dynamic simulations in
[78]. It will be a challenge for the current (and any other) large-strain
theory to reproduce such a microstructure.

Acknowledgments

AB acknowledges the support from SERB, Government of India (Project number
SRG/2020/001194), and IIT Tirupati. VL acknowledges the support from NSF
(CMMI-1943710 and DMR-1904830) and Iowa State University (Vance Coffman
Faculty Chair Professorship). The simulations were performed at Extreme
Science and Engineering Discovery Environment (XSEDE), allocations TG-
MSS140033 and MSS170015.

## References

  * [1]
  * [2] K. Bhattacharya, Microstructure of Martensite: Why It Forms and How It Gives Rise to the Shape-Memory Effect, Oxford University Press, Oxford, 2004.
  * [3] M. Pitteri, G. Zanzotto, Continuum Models for Phase Transitions and Twinning in Crystals, Chapman & Hall/CRC, Boca Raton, 2003.
  * [4] K. Adachi, C.M. Wayman, Transformation behavior of nearly stoichiometric Ni-Mn alloy, Metall. Trans. A 16 (1985) 1567-1579.
  * [5] J.M. Ball, R.D. James, Fine Phase Mixtures as Minimizers of Energy, Arch. Ration. Mech. Anal. 100 (1987) 13-52.
  * [6] K.Bhattacharya, Wedge-like microstructure in martensites, Acta Metall. Mater. 39 (1991) 2431-2444.
  * [7] D.A. Porter, K.E. Easterling, M.Y. Sherif, Phase Transformations in Metals and Alloys. CRC Press, Boca Raton, (2000).
  * [8] J. Wang, P. Steinmann, On the modeling of equilibrium twin interfaces in a single-crystalline magnetic shape memory alloy sample. I: theoretical formulation, Continuum Mech. Thermodyn. 26 (2014) 563-592.
  * [9] S. Govindjee, K. Hackl, R. Heinen, An upper bound to the free energy of mixing by twin-compatible lamination for $n$-variant martensitic phase transformations, Continuum Mech. Thermodyn. 18 (2007) 443-453.
  * [10] D. Schryvers, P. Boullaya, P. L. Potapov, R.V. Kohn, J.M. Ball, Microstructures and interfaces in Ni–Al martensite: comparing HRTEM observations with continuum theories, Int. J. Solids Struct. 39 (2002) 3543-3554.
  * [11] C. Chu, Hysteresis and microstructures: a study of biaxial loading on compound twins of copper-aluminum-nickel single crystals. PhD Dissertation, University of Minnesota (1993).
  * [12] R. Abeyaratne, C. Chu, R.D. James, Kinetics of materials with wiggly energies: Theory and application to the evolution of twinning microstructures in a Cu-Al-Ni shape memory alloy. Philos. Mag. A 73 (1996) 457-497.
  * [13] D. Schryvers, D. Holland-Moritz, Austenite and martensite microstructures in splat-cooled Ni-Al, Intermetallics, 6 (1998) 427-436.
  * [14] R.D. James, K.F. Hane. Martensitic transformations and shape memory materials, Acta Mater. 48 (2000) 197-222.
  * [15] K.F. Hane, T.W. Shield, Microstructure in a copper–aluminium–nickel shape-memory alloy, Proc. R. Soc. Lond. A 455 (1999) 3901-3915.
  * [16] G. Ruddock, A microstructure of martensite which is not a minimiser of energy: the X-interface, Arch. Rat. Mech. Anal. 127 (1994) 1-39.
  * [17] S. Stupkiewicz, A.G. Lengiewicz, Almost compatible X-microstructures in CuAlNi shape memory alloy, Continuum Mech. Thermodyn. 24 (2012) 149-164.
  * [18] A. Umantsev, Field Theoretic Method in Phase Transformations, Springer, Heidelberg, 2012.
  * [19] S.M. Allen, J.W. Cahn, A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening, Acta Mater. 27 (1979) 1085-1095.
  * [20] Y.M. Jin, A. Artemev, A.G. Khachaturyan, Three-dimensional phase field model of low-symmetry martensitic transformation in polycrystal: Simulation of $\zeta_{2}$ martensite in AuCd alloys, Acta Mater. 49 (2001) 2309-2320.
  * [21] A. Basak, V.I. Levitas, Nanoscale multiphase phase field approach for stress- and temperature-induced martensitic phase transformations with interfacial stresses at finite strains, J. Mech. Phys. Solids 113 (2018) 162-196.
  * [22] V.A. Levin, V.I. Levitas, K.M. Zingerman, E.I. Freiman, Phase-field simulation of stress-induced martensitic phase transformations at large strains, Int. J. Solids Struct. 50 (2013) 2914-2928.
  * [23] K. Tůma,, M. R. Hajidehib, J. Hrona, P.E. Farrellc, S. Stupkiewicz, Phase field modeling of multivariant martensitic transformation at finite-strain: Computational aspects and large-scale finite-element simulations, Comput. Methods Appl. Mech. Engrg. 377 (2021) 113705.
  * [24] A. Artemev, Y. Wang, A.G. Khachaturyan, Three-dimensional phase field model and simulation of martensitic transformation in multilayer systems under applied stresses, Acta. Mater. 48 (2000) 2503-2518.
  * [25] A. Artemev, Y. Jin, A.G. Khachaturyan, Three-dimensional phase field model of proper martensitic transformation, Acta. Mater. 49 (2001) 1165-1177.
  * [26] O.U. Salman, A. Finel, R. Delville, D. Schryvers, The role of phase compatibility in martensite, J. Appl. Phys. 111 (2012) 103517.
  * [27] A. Artemev, J. Slutsker, A.L. Roytburd, Phase field modeling of self-assembling nanostructures in constrained films, Acta. Mater. 53 (2005) 3425–3432.
  * [28] D.J. Seol, S.Y. Hu, Y.L. Li, L.Q. Chen, K.H. Oh, Computer simulation of martensitic transformation in constrained films, Mater. Sci. Forum 408-412 (2002) 1645-1650.
  * [29] T.W. Heo, L.Q. Chen, Phase-field modeling of displacive phase transformations in elastically anisotropic and inhomogeneous polycrystals, Acta. Mater. 76 (2014) 68-81.
  * [30] D. Fan, L.Q. Chen, Computer Simulation of Twin Formation during the Displacive c$\to$t’ Phase Transformation in the Zirconia‐Yttria System, J. Amer. Ceramic Soc. 78 (1995) 769-773.
  * [31] O. Shchyglo, G. Du, J.K. Engels, I. Steinbach, Phase-field simulation of martensite microstructure in low-carbon steel, Acta Mater. 175 (2019) 415-425.
  * [32] O. Shchyglo, U. Salman, A. Finel, Martensitic phase transformations in Ni–Ti-based shape memory alloys: The Landau theory, Acta Mater. 60 (2012) 6784-6792.
  * [33] A. Planes, P. Lloveras, T. Castán, A. Saxena, M. Porta, Ginzburg–Landau modelling of precursor nanoscale textures in ferroelastic materials, Continuum Mech. Thermodyn. 24 (2012) 619-627.
  * [34] J. Kundin, D. Raabe, H. Emmerich, A phase-field model for incoherent martensitic transformations including plastic accommodation processes in the austenite, J. Mech. Phys. Solids 59 (2011) 2082-2102.
  * [35] M. Mamivand, M.A. Zaeem, H.E. Kadiri, A review on phase field modeling of martensitic phase transformation, Comput. Mater. Sci. 77 (2013) 304-311.
  * [36] V.I. Levitas, M. Javanbakht, Surface tension and energy in multivariant martensitic transformations: Phase-field theory, simulations, and model of coherent interface, Phys. Rev. Lett. 105 (2010) 165701.
  * [37] A.Y. Woldman, C.M. Landis, Phase-field modeling of ferroelectric to paraelectric phase boundary structures in single-crystal barium titanate, Smart Mater. Struct., 25 (2016) 035033.
  * [38] V.I. Levitas, V.A. Levin, K.M. Zingerman, E.I. Freiman, Displacive phase transitions at large strains: Phase-field theory and simulations, Phys. Rev. Lett. 103 (2009) 025702.
  * [39] V.I. Levitas, M. Javanbakht, Phase-field approach to martensitic phase transformations: Effect of martensite-martensite interface energy, Int. J. Mater. Res. 102 (2011) 652-665.
  * [40] J.C.H. Lei, L.J. Li, Y.C. Shu, J.Y. Li, Austenite-martensite interface in shape memory alloys, Appl. Phys. Lett. 96 (2010) 141910.
  * [41] F.E. Hildebrand, C. Miehe, A phase field model for the formation and evolution of martensitic laminate microstructure at finite strains, Phil. Mag. 92 (2012) 1-41.
  * [42] J.D. Clayton, J. Knap, A phase field model of deformation twinning: Nonlinear theory and numerical simulations, Physica D 240 (2011a) 841-858.
  * [43] J.D. Clayton, Nonlinear thermodynamic phase field theory with application to fracture and dynamic inelastic phenomena in ceramic polycrystals. J. Mech. Phys. Solids 157 (2021) 104633.
  * [44] J.D. Clayton, J. Knap, Continuum modeling of twinning, amorphization, and fracture: theory and numerical simulations, Continuum Mech. Thermodyn. 30 (2018) 421-455.
  * [45] H. She, Y. Liu, B. Wang, D. Ma Finite element simulation of phase field model for nanoscale martensitic transformation, Comput. Mech. 52 (2013) 949-958.
  * [46] V.I. Levitas, D.L. Preston, Three-dimensional Landau theory for multivariant stress-induced martensitic phase transformations. I. Austenite $\leftrightarrow$ Martensite, Phys. Rev. B 66 (2002) 134206.
  * [47] V.I. Levitas, D.L. Preston, Three-dimensional Landau theory for multivariant stress-induced martensitic phase transformations. II. Multivariant phase transformations and stress-space analysis, Phys. Rev. B 66 (2002) 134207.
  * [48] V.I. Levitas, Phase-field theory for martensitic phase transformations at large strains, Int. J. Plast. 49 (2013) 85-118.
  * [49] V.I. Levitas, Phase field approach for stress- and temperature-induced phase transformations that satisfies lattice instability conditions. Part 1. General theory, Int. J. Plast. 106 (2018) 164-185.
  * [50] 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, Int. J. Plast. 107 (2018) 223-245.
  * [51] K. Tůma, S. Stupkiewicz, Phase-field study of size-dependent morphology of austenite-twinned martensite interface in CuAlNi, Int. J. Solids Struct. 97-98 (2016) 89-100.
  * [52] K. Tůma, S. Stupkiewicz, H. Petryk, Size effects in martensitic microstructures: Finite-strain phase field model versus sharp-interface approach, J. Mech. Phys. Solids 95 (2016) 284-307.
  * [53] K. Tůma, S. Stupkiewicz, H. Petryk, Rate-independent dissipation in phase-field modelling of displacive transformations, J. Mech. Phys. Solids 114 (2018) 117-142.
  * [54] M.R. Hajidehi, S. Stupkiewicz, Phase-field modeling of multivariant martensitic microstructures and size effects in nano-indentation, Mech. Mater. 141 (2020) 103267.
  * [55] V.I. Levitas, A.V. Idesman, D.L. Preston, Microscale simulation of martensitic microstructure evolution, Phys. Rev. Lett. 93 (2004) 105701.
  * [56] A.V. Idesman, V.I. Levitas, D.L. Preston, J.-Y. Cho, Finite element simulations of martensitic phase transitions and microstructure based on strain softening model, J. Mech. Phys. Solids 53 (2005) 495-523.
  * [57] Levitas V.I., Esfahani S.E., and Ghamarian I. Scale-free modeling of coupled evolution of discrete dislocation bands and multivariant martensitic microstructure. Phys. Rev. Lett. 121 (2018) 205701.
  * [58] Esfahani S.E., Ghamarian I., and Levitas V.I. Strain-induced multivariant martensitic transformations: A scale-independent simulation of interaction between localized shear bands and microstructure. Acta Mate. 196 (2020) 430-443.
  * [59] Babaei H. and Levitas V.I. Finite-strain scale-free phase-field approach to multivariant martensitic phase transformations with stress-dependent effective thresholds. J. Mech. Phys. Solids 144 (2020) 104114.
  * [60] V.I. Levitas, A.M. Roy, Multiphase phase field theory for temperature-induced phase transformations: Formulation and application to interfacial phases, Acta Mater. 105 (2016) 244-257.
  * [61] V.I. Levitas, A.M. Roy, D.L. Preston, Multiple twinning and variant-variant transformations in martensite: Phase-field approach, Phys. Rev. B 88 (2013) 054113.
  * [62] I. Steinbach, F. Pezzolla, B. Nestler, M.Seeßelberg, R.Prieler, G.J.Schmitz, J.L.L.Rezende. A phase field concept for multiphase systems, Physica D. 94 (1996) 135-147.
  * [63] V.I. Levitas, A.M. Roy, Multiphase phase field theory for temperature- and stress-induced phase transformations, Phys. Rev. B 91 (2015) 174109.
  * [64] J.Y. Cho, A.V. Idesman, V.I. Levitas, T. Park, Finite element simulations of dynamics of multivariant martensitic phase transitions based on Ginzburg-Landau theory, Int. J. Solids Struct. 49 (2012) 1973-1992.
  * [65] Levin, V.A., Zingerman, K.M., 1998. Interaction and Microfracturing Pattern for Successive Origination (Introduction) of Pores in Elastic Bodies: Finite Deformation. Trans. ASME. J. Appl. Mech. 65, 431-435.
  * [66] Levin, V.A., 1999. Multiple superposition of large deformation in elastic and viscoelastic bodies. Nauka, Moscow. (in Russian).
  * [67] Levin, V.A., 1998. Theory of Repeated Superposition of Large Deformations. Elastic and Viscoelastic Bodies. Int. J. Solids Struct. 35, 2585-2600.
  * [68] V.I. Levitas, Phase field approach to martensitic phase transformations with large strains and interface stresses, J. Mech. Phys. Solids 70 (2014) 154-189.
  * [69] Levitas V.I. and Warren J. A. Phase field approach with anisotropic interface energy and interface stresses: large strain formulation. J. Mech. Phys. Solids 91 (2016) 94-125.
  * [70] P.C. Bollada, P.K. Jimack, A.M. Mullis, Multiphase field modelling of alloy solidification, Computat. Mater. Sci. 171 (2020) 109085.
  * [71] G.I. Tóth, T. Pusztai, L. Granasi, Consistent multiphase-field theory for interface driven multidomain dynamics, Phys. Rev. B 92 (2015) 184105.
  * [72] C.S. Jog, Continuum Mechanics: Foundations and Applications of Mechanics (Volume I), Cambridge Univ. Press, New Delhi, 2015.
  * [73] I. Steinbach, Phase-field models in materials science, Model. Simul. Mater. Sci. Eng. 17 (2009) 073001.
  * [74] A. Basak, V.I. Levitas, Finite element procedure and simulations for a multiphase phase field approach to martensitic phase transformations at large strains and with interfacial stresses. Comput. Methods Appl. Mech. Engrg. 343 (2019) 368–406.
  * [75] W. Bangerth, D. Davydov, T. Heister, L. Heltai, G. Kanschat, M. Kronbichler, M. Maier, B. Turcksin, D. Wells, The deal.II library, version 8.4, J. Numer. Math. 24 (2016).
  * [76] S.M. Ueland, C.A. Schuh, Grain boundary and triple junction constraints during martensitic transformation in shape memory alloys, J. Appl. Phys. 114(39) (2013) 053503.
  * [77] A. Basak, V.I. Levitas, Interfacial stresses within boundary between martensitic variants: analytical and numerical finite strain solutions for three phase field models, Acta Mater. 139 (2017) 174-187.
  * [78] Chen H., Levitas V.I., Popov D., and Velisavljevic N. Nontrivial nanostructure, stress relaxation mechanisms, and crystallography for pressure-induced Si-I $\rightarrow$ Si-II phase transformation. Nat. Commun. 2022, DoI: 10.1038/s41467-022-28604-1.