# A model for the interaction of dislocations with planar defects based on
Allen-Cahn type microstructure evolution coupled to strain gradient elasticity

M. Budnitzki<EMAIL_ADDRESS>S. Sandfeld<EMAIL_ADDRESS>Institute for Advanced Simulation (IAS-9: Materials Data Science and
Informatics), Forschungszentrum Jülich GmbH, 52428 Jülich, Germany  TU
Bergakademie Freiberg, Institute of Mechanics and Fluid Dynamics,
Lampadiusstr. 4, 09599 Freiberg

###### Abstract

In classical elasticity theory the stress-field of a dislocation is
characterized by a $1/r$-type singularity. When such a dislocation is
considered together with an Allen-Cahn-type phase-field description for
microstructure evolution this leads to singular driving forces for the order
parameter, resulting in non-physical (and discretization-dependent)
predictions for the interaction between dislocations and phase-, twin- or
grain-boundaries. We introduce a framework based on first strain gradient
elasticity to regularize the dislocation core. It is shown that the use of
strain energy density that is quadratic in the gradient of elastic deformation
results in non-singular stresses but may result in singular driving forces,
whereas a strain energy, which is quadratic in the gradient of the full
deformation tensor, regularizes both stresses and driving forces for the order
parameter and is therefore a suitable choice. The applicability of the
framework is demonstrated using a comprehensive example.

###### keywords:

strain gradient elasticity , phase field , dislocation

††journal: JMPS

capbtabboxtable[][]

##  1 Introduction

Phase field approaches have proven to be very powerful for the investigation
of the formation and evolution of microstructures due to solid-solid phase
transformations and twinning. This appears to be the natural framework for the
investigation of the interaction of planar crystal defects such as phase- or
twin-boundaries with line defects (dislocations, disclinations). A typical
phase field model for diffusionless (martensitic) transformations comprises of
evolution equations of Allen-Cahn-type for the order parameters $\phi_{\beta}$

$M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}-\rho\partial_{\phi_{\beta}}\psi\,,$
(1)

where $M$ and $\alpha$ are constants, $\rho$ denotes the mass density, and
$\psi$ is a bulk specific free energy. The subscript $\beta$ indicates the
number of the phase, grain or twin variant. Assuming a small perturbation
setting, the linear strain tensor111Nomenclature: We denote vectors by bold
lower case latin $\boldsymbol{a}$ and greek $\boldsymbol{\alpha}$ letters. The
dot operator “$\cdot$” denotes the scalar product. Second order tensors are
denoted by bold uppercase latin letters $\boldsymbol{A}$. We introduce a
scalar product between second order tensors denoted by “:” as
$\boldsymbol{A}:\boldsymbol{B}:=\text{tr}\boldsymbol{A}\cdot\boldsymbol{B}^{\top}$,
where $\boldsymbol{B}^{\top}$ is the transpose of $\boldsymbol{B}$ and
$\text{tr}(\cdot)$ denotes the trace operator. Similarly, we denote third
order tensors $\boldsymbol{\mathcal{A}}$ by bold calligraphic capital letters
and “$\smash{\,\smash{\vdots}\,}$” is the corresponding scalar product. We use
black-board capital letters $\boldsymbol{\mathbb{C}}$ for fourth-order
tensors. Whenever index-notation is used, summation over latin indices
appearing twice is implied and spatial derivatives are denoted using the comma
operator, e.g. $\partial_{x_{i}}y\equiv y^{,i}$. $\boldsymbol{E}$ can be
additively decomposed into elastic $\boldsymbol{E}^{\text{e}}$ and inelastic
(i.e., eigenstrain) $\boldsymbol{E}^{\text{in}}(\phi_{\beta})$ contributions,
such that
$\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)=\boldsymbol{E}-\boldsymbol{E}^{\text{in}}(\phi_{\beta})$.
Assuming linear elasticity, the stress $\boldsymbol{S}$ is given by
$\boldsymbol{S}=\mathbb{C}:\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)$,
and the specific free energy takes the form

$\psi\bigl{(}\boldsymbol{E}\,,\,\phi_{\beta}\,,\,\theta\bigr{)}=\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\mathbb{C}:\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)+\psi_{\text{b}}\bigl{(}\phi_{\beta}\,,\,\theta\bigr{)}\,.$
(2)

As a consequence, the evolution equation (1) can be rewritten as

$M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\rho\partial_{\phi_{\beta}}\psi_{\text{b}}\,.$
(3)

In linear elastic Volterra theory, the stresses diverge as the dislocation
line is approached. In particular for dislocations the singularity is of
$1/r$-type. As per Eq. (3), this results in singular driving forces for the
evolution of the order parameters, effectively negating the concepts such as a
nucleation barrier or a pile-up stress. Different approaches to regularize the
stress in the core region exist in literature based either on the concept of a
distributed Burger’s vector (Lothe, 1992; Cai et al., 2006), which are
inspired by richer microscopic models for dislocations (Peierls, 1940;
Nabarro, 1947), or generalized continuum theories (Lazar et al., 2005, 2006;
Lazar and Po, 2015; Po et al., 2018). However, the first strain gradient
approach advocated by Po et al. (2018) has the advantage that the obtained
regularization is independent of the type of defect in question and therefore
does not require any defect-specific information for the determination of
model parameters. In principle, these parameters can directly be obtained from
atomistic interaction potentials (Admal et al., 2017).

The purpose of this work is to follow a micromorphic approach and to derive a
framework which consistently couples first strain gradient elasticity to
Allen-Cahn-type microstructure evolution ensuring non-singular driving forces
on the order parameters in the presence of line defects.

##  2 Balance equations and boundary conditions

The principle of virtual power (PVP) provides a systematic way of deriving
field equations and boundary conditions for arbitrary mechanical and coupled
problems (cf. Maugin, 1980; Germain, 1973; Del Piero, 2009). In the present
work it is used in the following form: The virtual power of the inertia forces
$\mathscr{P}^{*}_{\text{a}}$ balances the virtual power
$\mathscr{P}^{*}_{\text{int}}$ of the internal and
$\mathscr{P}^{*}_{\text{ext}}$ of the external forces acting on any sub-domain
$\mathscr{S}$ of the material body $\mathscr{B}$ for any admissible virtual
velocity field $\boldsymbol{v}^{*}$ and virtual rate of order parameter field
${\dot{\phi\mkern 5.0mu}\mkern-5.0mu}{}^{*}_{\beta}$, i.e.,

$\mathscr{P}^{*}_{\text{a}}=\mathscr{P}^{*}_{\text{int}}+\mathscr{P}^{*}_{\text{ext}}\,.$
(4)

For the sake of simplicity we disregard any higher order inertia terms Mindlin
(1964) as well as inertial forces acting on the order parameter, resulting in

$\mathscr{P}^{*}_{\text{a}}=\int_{\mathscr{S}}\rho\dot{\boldsymbol{v}}\cdot\boldsymbol{v}^{*}\,\text{d}V\,.$
(5)

The power of internal forces is given by

$\mathscr{P}^{*}_{\text{int}}=-\int_{\mathscr{S}}\left(\boldsymbol{S}^{\top}:\boldsymbol{L}^{*}+\boldsymbol{\mathcal{T}}\,\smash{\vdots}\,\operatorname*{grad}{\boldsymbol{L}^{*}}-\pi_{\beta}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}+\boldsymbol{\xi}_{\beta}\cdot\operatorname*{grad}{{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}}\right)\,\text{d}V\,,$ (6)

with $\boldsymbol{L}^{*}:=\operatorname*{grad}{\boldsymbol{v}^{*}}$. Here
$\boldsymbol{S}$ and $\boldsymbol{\mathcal{T}}$ are the Cauchy and higher
order stresses, respectively, while $\pi_{\beta}$ and
$\boldsymbol{\xi}_{\beta}$ are thermodynamic forces that directly correspond
to the internal microforce and microstress introduced by Gurtin (1996). We
note that the invariance requirement of $\mathscr{P}^{*}_{\text{int}}$ with
respect to superimposed rigid body motions is satisfied sufficiently by
assuming $\boldsymbol{S}=\boldsymbol{S}^{\top}$ and
$\boldsymbol{\mathcal{T}}\cdot\boldsymbol{a}=(\boldsymbol{\mathcal{T}}\cdot\boldsymbol{a})^{\top}$
for arbitrary vectors $\boldsymbol{a}$. For the power of external forces we
consider the very simple case of no body or contact forces acting on
$\boldsymbol{L}^{*}$ and $\operatorname*{grad}{{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}}$, and only a contact (micro)force
$\zeta_{\beta}$ acting ${\dot{\phi\mkern 5.0mu}\mkern-5.0mu}{}^{*}_{\beta}$

$\mathscr{P}^{*}_{\text{ext}}=\int_{\mathscr{S}}\boldsymbol{f}\cdot\boldsymbol{v}^{*}\rho\,\text{d}V+\int_{\partial\mathscr{S}}\left(\boldsymbol{t}\cdot\boldsymbol{v}^{*}+\zeta_{\beta}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}\right)\,\text{d}a\,.$ (7)

In order to obtain the consequences of the PVP, the integrals in Eq. (6) are
transformed using the following identities

$\displaystyle\operatorname*{div}{(\boldsymbol{S}\cdot\boldsymbol{v}^{*})}=(\operatorname*{div}{\boldsymbol{S}})\cdot\boldsymbol{v}^{*}+\boldsymbol{S}:\boldsymbol{L}^{*}\,,$
(8)
$\displaystyle\operatorname*{div}{(\boldsymbol{\mathcal{T}}:\boldsymbol{L}^{*})}=(\operatorname*{div}{\boldsymbol{\mathcal{T}}}):\boldsymbol{L}^{*}+\boldsymbol{\mathcal{T}}\,\smash{\vdots}\,\operatorname*{grad}{\boldsymbol{L}^{*}}\,,$
(9)
$\displaystyle\operatorname*{div}{\bigl{(}(\operatorname*{div}{\boldsymbol{\mathcal{T}}})\cdot\boldsymbol{v}^{*}\bigr{)}}=(\operatorname*{div}{\operatorname*{div}{\boldsymbol{\mathcal{T}}}})\cdot\boldsymbol{v}^{*}+(\operatorname*{div}{\boldsymbol{\mathcal{T}}}):\boldsymbol{L}^{*}\,,$
(10)
$\displaystyle\operatorname*{div}{(\boldsymbol{\xi}_{\beta}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta})}=(\operatorname*{div}{\boldsymbol{\xi}_{\beta}})\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}+\boldsymbol{\xi}_{\beta}\cdot\operatorname*{grad}{{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}}\,,$ (11)

and the divergence theorem, resulting in

$\mathscr{P}^{*}_{\text{int}}=\int_{\mathscr{S}}\bigl{(}\operatorname*{div}{\boldsymbol{S}}-\operatorname*{div}{\operatorname*{div}{\boldsymbol{\mathcal{T}}}}\bigr{)}\cdot\boldsymbol{v}^{*}\,\text{d}V-\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\bigl{(}\boldsymbol{S}^{\top}-\operatorname*{div}{\boldsymbol{\mathcal{T}}}\bigr{)}\cdot\boldsymbol{v}^{*}\,\text{d}a-\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}:\boldsymbol{L}^{*}\,\text{d}a\\\
+\int_{\mathscr{S}}\bigl{(}\pi_{\beta}+\operatorname*{div}{\boldsymbol{\xi}_{\beta}}\bigr{)}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}\,\text{d}V-\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\xi}_{\beta}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}^{*}_{\beta}\,\text{d}a\,.$ (12)

Introducing the surface gradient operator

$\operatorname*{grad_{S}}(\cdot)=\operatorname*{grad}(\cdot)-\partial_{\boldsymbol{n}}(\cdot)\otimes\boldsymbol{n}\,,$
(13)

where $\partial_{\boldsymbol{n}}$ is the directional derivative in the
direction of the outward normal $\boldsymbol{n}$, the third integral in
expression (12) can be rewritten as

$\displaystyle\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}:\boldsymbol{L}^{*}\,\text{d}a$
$\displaystyle=\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}:\operatorname*{grad_{S}}{\boldsymbol{v}^{*}}\,\text{d}a+\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}:\partial_{\boldsymbol{n}}\boldsymbol{v}^{*}\otimes\boldsymbol{n}\,\text{d}a$
(14)
$\displaystyle=\int_{\partial\mathscr{S}}\operatorname*{div_{S}}{\bigl{(}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}\cdot\boldsymbol{v}^{*}\bigr{)}}\,\text{d}a-\int_{\partial\mathscr{S}}\operatorname*{div_{S}}{\bigl{(}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}\bigr{)}}\cdot\boldsymbol{v}\,\text{d}a+\int_{\partial\mathscr{S}}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}:\partial_{\boldsymbol{n}}\boldsymbol{v}^{*}\otimes\boldsymbol{n}\,\text{d}a,.$

Finally, applying the surface divergence theorem and, for the sake of
simplicity, neglecting any wedge line and corner contributions, we find

$\int_{\partial\mathscr{S}}\operatorname*{div_{S}}{\bigl{(}\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}\cdot\boldsymbol{v}^{*}\bigr{)}}\,\text{d}a=\int_{\partial\mathscr{S}}\bigl{(}\operatorname*{div_{S}}{\boldsymbol{n}}\bigr{)}\boldsymbol{n}\otimes\boldsymbol{n}:\boldsymbol{\mathcal{T}}\cdot\boldsymbol{v}^{*}\,\text{d}a\,.$
(15)

Enforcing Eq. (4) we arrive after a number of straight forward algebraic
manipulations at the following field equations on $\mathscr{B}$

$\displaystyle\rho\dot{\boldsymbol{v}}=\operatorname*{div}{\left(\boldsymbol{S}-\operatorname*{div}\boldsymbol{\mathcal{T}}\right)}+\rho\boldsymbol{f}\,,$
(16a) $\displaystyle
0=\operatorname*{div}{\boldsymbol{\xi}_{\beta}}+\pi_{\beta}\,,$ (16b) and
boundary conditions on $\partial\mathscr{B}$
$\displaystyle\boldsymbol{t}=\left(\boldsymbol{S}-\operatorname*{div}\boldsymbol{\mathcal{T}}\right)\cdot\boldsymbol{n}-\operatorname*{div_{S}}{\left(\boldsymbol{n}\cdot\boldsymbol{\mathcal{T}}\right)}\,,$
(16c)
$\displaystyle\zeta_{\beta}=\boldsymbol{\xi}_{\beta}\cdot\boldsymbol{n}\,.$
(16d)

We note that, introducing the total stress

$\boldsymbol{S}_{\text{t}}:=\boldsymbol{S}-\operatorname*{div}{\boldsymbol{\mathcal{T}}}\,,$
(17)

the balance of linear momentum (16a) regains its standard form for simple
materials

$\displaystyle\rho\dot{\boldsymbol{v}}=\operatorname*{div}{\boldsymbol{S}_{\text{t}}}+\rho\boldsymbol{f}\,,$
(18)

which is convenient for the numerical implementation.

##  3 Constitutive equations

The following equations are formulated assuming a geometrically linear
setting, i.e., the displacement gradient is considered to be small
$||\operatorname*{grad}{\boldsymbol{u}}||\ll 1$. In this case the deformation
is characterized by the linear strain tensor
$\boldsymbol{E}=\frac{1}{2}\left(\operatorname*{grad}{\boldsymbol{u}}+(\operatorname*{grad}{\boldsymbol{u}})^{\top}\right)$.
Its gradient will be denoted by
$\boldsymbol{\mathcal{Y}}:=\operatorname*{grad}\boldsymbol{E}$.

###  3.1 Laws of state

We choose the following ansatz for the specific free energy and thermodynamic
forces

$\displaystyle\psi=\psi\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)\,,$
$\displaystyle\boldsymbol{S}=\boldsymbol{S}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)\,,$
$\displaystyle\boldsymbol{\mathcal{T}}=\boldsymbol{\mathcal{T}}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)\,,$
$\displaystyle\pi_{\beta}=\pi_{\beta}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\,,\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}\right)\,,$
$\displaystyle\boldsymbol{\xi}_{\beta}=\boldsymbol{\xi}_{\beta}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)\,.$

The second law of the thermodynamics in the form of the Clausius-Duhem
inequality given for the isothermal case by

$\left(\boldsymbol{S}-\rho\partial_{\boldsymbol{E}}\psi\right):{\dot{\boldsymbol{E}\mkern
5.0mu}\mkern-5.0mu}{}+\left(\boldsymbol{\mathcal{T}}-\rho\partial_{\boldsymbol{\mathcal{Y}}}\psi\right)\,\smash{\vdots}\,{\dot{\boldsymbol{\mathcal{Y}}\mkern
5.0mu}\mkern-5.0mu}{}-\left(\pi_{\beta}+\rho\partial_{\phi}\psi\right){\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}+\left(\boldsymbol{\xi}_{\beta}-\rho\partial_{\operatorname*{grad}\phi_{\beta}}\psi\right)\cdot\operatorname*{grad}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}\geqslant 0$ (19)

can be exploited using the classical Coleman-Noll procedure to arrive at the
laws of state

$\displaystyle\boldsymbol{S}=\rho\partial_{\boldsymbol{E}}\psi\,,$
$\displaystyle\boldsymbol{\mathcal{T}}=\rho\partial_{\boldsymbol{\mathcal{Y}}}\psi\,,$
$\displaystyle\boldsymbol{\xi}_{\beta}=\rho\partial_{\operatorname*{grad}\phi_{\beta}}\psi$
(20)

and the residual dissipation inequality

$\displaystyle-\pi^{\text{d}}_{\beta}\,{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}\geqslant
0\,,\quad\text{with}\quad\pi^{\text{d}}_{\beta}:=\pi_{\beta}+\rho\partial_{\phi_{\beta}}\psi\,.$
(21)

###  3.2 Free energy and dissipation potential

As customary in phase field models for solid-solid transformations, the
specific free energy can be split into an elastic, a bulk chemical and an
interface contribution

$\psi=\psi_{\text{e}}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\theta\right)+\psi_{\text{b}}\left(\phi_{\beta}\,,\,\theta\right)+\psi_{\text{i}}\left(\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)\,,$
(22)

as indicated by the subscripts “e” (elastic), “b” (bulk chemical) and “i”
(interface). In our formulation, the elastic free energy is of Helmholtz-type,
i.e.,

$\displaystyle\rho\psi_{\text{e}}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\theta\right)=\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)+\frac{1}{2}\bigl{(}\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\bigr{)}\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}\,,\text{
or}$ (23)
$\displaystyle\rho\psi_{\text{e}}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\theta\right)=\frac{1}{2}\mathbb{C}^{ijkl}(\phi_{\beta})E_{ij}^{\text{e}}(\boldsymbol{E}\,,\,\phi_{\beta})E_{kl}^{\text{e}}(\boldsymbol{E}\,,\,\phi_{\beta})+\frac{1}{2}\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})\boldsymbol{\mathcal{Y}}_{ij}^{n}\boldsymbol{\mathcal{Y}}_{kl}^{m}\,$
(24)

where $\boldsymbol{E}^{\text{in}}(\phi_{\beta})$ is the inelastic strain,
$\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):=\boldsymbol{E}-\boldsymbol{E}^{\text{in}}(\phi_{\beta})$
is the elastic strain, $\mathbb{C}(\phi_{\beta})$ the stiffness tensor and
$\boldsymbol{\Lambda{}}(\phi_{\beta})$ a gradient length scale tensor (cf. Po
et al., 2018). The specific choice of functional dependence of
$\boldsymbol{E}^{\text{in}}(\phi_{\beta})$,
$\psi_{\text{b}}\left(\phi_{\beta}\,,\,\theta\right)$ and
$\psi_{\text{i}}\left(\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right)$
on the order parameter $\phi_{\beta}$ is of no relevance at this point;
however, we will assume that the interface energy is of the form

$\displaystyle\rho\psi_{\text{i}}\left(\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right):=\frac{\alpha}{2}\,||\operatorname*{grad}\phi_{\beta}||^{2}+g(\phi_{\beta}\,,\,\theta)\,,$
$\displaystyle\rho\psi_{\text{i}}\left(\phi_{\beta}\,,\,\operatorname*{grad}\phi_{\beta}\,,\,\theta\right):=\frac{\alpha}{2}\,\phi_{\beta}^{,i}\phi_{\beta}^{,i}+g(\phi_{\beta}\,,\,\theta)\,.$
(25)

Using the laws of state (20) we immediately find

$\displaystyle\boldsymbol{S}=\mathbb{C}(\phi_{\beta}):\bigl{(}\boldsymbol{E}-\boldsymbol{E}^{\text{in}}(\phi_{\beta})\bigr{)}\,,$
$\displaystyle
S^{ij}=\mathbb{C}^{ijkl}(\phi_{\beta})\left(E_{kl}-E_{kl}^{\text{in}}(\phi_{\beta})\right)\,,$
(26a)
$\displaystyle\boldsymbol{\mathcal{T}}=\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\,,$
$\displaystyle\boldsymbol{\mathcal{T}}^{ij}_{n}=\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})\boldsymbol{\mathcal{Y}}_{kl}^{m}\,,$
(26b)
$\displaystyle\boldsymbol{\xi}_{\beta}=\alpha\operatorname*{grad}\phi_{\beta}\,,$
$\displaystyle\xi_{\beta}^{i}=\alpha\phi_{\beta}^{,i}\,,$ (26c)

and combining the first two equations

$\displaystyle\boldsymbol{\mathcal{T}}=\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\left(\mathbb{C}^{-1}(\phi_{\beta}):\boldsymbol{S}\right)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})+\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\boldsymbol{E}^{\text{in}}(\phi_{\beta})}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\,,\text{
or}$ (27)
$\displaystyle\boldsymbol{\mathcal{T}}^{ij}_{n}=\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})\left(\mathbb{C}^{-1}_{klpq}(\phi_{\beta})S^{pq}\right)^{,m}+\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})E_{kl}^{\text{in},m}(\phi_{\beta})\,.$
(28)

Equation (17) can now be used in two ways: In conjunction with the laws of
state (26a) and (26b) it is a constitutive equation for the total stress
$\boldsymbol{S}_{\text{t}}$, which enters the balance of linear momentum (18)

$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right]}\,,\text{
or}$ (29) $\displaystyle
S_{\text{t}}^{ij}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}^{ijkl}(\phi_{\beta}):E^{\text{e}}_{kl}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\left(\mathbb{C}^{ijkl}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}_{kl}^{m}\Lambda_{mn}\right)^{,n}\,.$
(30)

When combined with Eq. (27), Eq. (17) can be used to determine the true stress
$\boldsymbol{S}$ from the total stress $\boldsymbol{S}_{\text{t}}$

$\displaystyle\boldsymbol{S}-\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\left(\mathbb{C}^{-1}(\phi_{\beta}):\boldsymbol{S}\right)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right]}=\boldsymbol{S}_{\text{t}}+\operatorname*{div}{\bigl{(}\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\boldsymbol{E}^{\text{in}}(\phi_{\beta})}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\bigr{)}}\,,\text{
or}$ (31) $\displaystyle
S^{ij}-\left[\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})\left(\mathbb{C}^{-1}_{klpq}(\phi_{\beta})S^{pq}\right)^{,m}\right]^{,n}=S_{\text{t}}^{ij}+\left[\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})E_{kl}^{\text{in},m}(\phi_{\beta})\right]^{,n}\,.$
(32)

In order to complete the phase field formulation we require a constitutive
equation for $\pi^{\text{d}}_{\beta}$, which is obtained in the spirit of
classical irreversible thermodynamics as

${\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=-\partial_{\pi^{\text{d}}_{\beta}}\Omega\left(\pi^{\text{d}}_{\beta}\right)$
(33)

from a dissipation potential $\Omega\left(\pi^{\text{d}}_{\beta}\right)$ that
is homogeneous of degree two

$\Omega\left(\pi^{\text{d}}_{\beta}\right):=\frac{1}{2}M\left(\pi^{\text{d}}_{\beta}\right)^{2}\,,$
(34)

where $M$ is the so called mobility constant. Combining equations (16b), (21),
(26c), (33) and (34) we find the classical Allen-Cahn equation

$\displaystyle M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}-\rho\partial_{\phi_{\beta}}\psi\,,$
$\displaystyle M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\phi_{\beta}^{,ii}-\rho\partial_{\phi_{\beta}}\psi\,,$
(35)

or, explicitely writing down the partial derivatives of $\psi$,

$M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\frac{1}{2}\left(\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\cdot\partial_{\phi_{\beta}}\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\\\
-\frac{1}{2}\left(\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\rho\partial_{\phi_{\beta}}\psi_{\text{b}}(\phi_{\beta}\,,\,\theta)-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,,$
(36)

or

$M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\phi_{\beta}^{,ii}+S^{ij}\partial_{\phi_{\beta}}E^{\text{in}}_{kl}(\phi_{\beta})-\frac{1}{2}\partial_{\phi_{\beta}}\mathbb{C}_{ijkl}(\phi_{\beta})E^{\text{e}}_{ij}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)E^{\text{e}}_{kl}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\frac{1}{2}\partial_{\phi_{\beta}}\mathbb{C}^{ijkl}(\phi_{\beta})\Lambda_{mn}(\phi_{\beta})\boldsymbol{\mathcal{Y}}_{ij}^{n}\boldsymbol{\mathcal{Y}}_{kl}^{m}-\\\
-\frac{1}{2}\mathbb{C}^{ijkl}(\phi_{\beta})\partial_{\phi_{\beta}}\Lambda_{mn}(\phi_{\beta})\boldsymbol{\mathcal{Y}}_{ij}^{n}\boldsymbol{\mathcal{Y}}_{kl}^{m}-\rho\partial_{\phi_{\beta}}\psi_{\text{b}}(\phi_{\beta}\,,\,\theta)-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,,$

Note that all terms that appear in the driving force, and as per Lazar et al.
(2005) the Cauchy stress $\boldsymbol{S}$ in particular, are non-singular even
in the presence of dislocations. Interestingly, this is not true for an
elastic specific free energy that is quadratic in
$\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right):=\operatorname*{grad}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)$
rather that $\boldsymbol{\mathcal{Y}}$ (cf. A).

###  3.3 Formulation for specific cases

For phase transformations the crystal lattice on both sides of the interface
will, in general, be different leading to different elastic properties and a
different shape of the dislocation core. In this case the equations (18),
(29), (31) and (36) retain their full complexity. However, the strength of the
general formulation is that it also covers simplified special cases. In the
following, we consider scenarios for which these equations can be strongly
reduced and which therefore elucidates the structure of the whole formalism.

####  3.3.1 Homogeneous bulk material

In the bulk phase the order parameter does not vary in space, i.e.,
$\operatorname*{grad}\phi_{\beta}=\boldsymbol{0}$,
$\mathbb{C}(\phi_{\beta})=\mathbb{C}$,
$\boldsymbol{\Lambda{}}(\phi_{\beta})=\boldsymbol{\Lambda{}}$,
$\boldsymbol{E}^{\text{in}}(\phi_{\beta})=\boldsymbol{0}$. The Allen-Cahn
equation is fulfilled automatically and Eqs. (31) and (29) recover the form
derived by Po et al. (2018)

$\displaystyle\boldsymbol{S}-\operatorname*{div}{\bigl{(}(\operatorname*{grad}{\boldsymbol{S}})\cdot\boldsymbol{\Lambda{}}}\bigr{)}=\boldsymbol{S}_{\text{t}}\,,$
with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\bigr{)}=\mathbb{C}:\left[\boldsymbol{E}-\operatorname*{div}{\bigl{(}\boldsymbol{\mathcal{Y}}\cdot\boldsymbol{\Lambda{}}}\bigr{)}\right]\,.$
(37a) For materials with cubic symmetry the gradient length scale tensor
$\boldsymbol{\Lambda{}}$ is isotropic, i.e.,
$\boldsymbol{\Lambda{}}=\mathcal{l}^{2}\boldsymbol{I}$, and the above
expressions can be further simplified to the form derived by Lazar et al.
(2005)
$\displaystyle\boldsymbol{S}-\mathcal{l}^{2}\Delta\boldsymbol{S}=\boldsymbol{S}_{\text{t}}\,,$
with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\bigr{)}=\mathbb{C}:\bigl{(}\boldsymbol{E}-\mathcal{l}^{2}\operatorname*{div}{\boldsymbol{\mathcal{Y}}}\bigr{)}=\mathbb{C}:\bigl{(}\boldsymbol{E}-\mathcal{l}^{2}\Delta\boldsymbol{E}\bigr{)}\,.$
(37b)

####  3.3.2 Boundaries between grains without inelastic strain

The crystal lattices on both sides of a grain boundary differ only by a
rotation $\boldsymbol{Q}(\phi_{\beta})$. Hence, we assume that the chemical
bulk energy is independent of the order parameter, i.e.,
$\psi_{\text{b}}\left(\phi_{\beta}\,,\,\theta\right)=\psi_{\text{b}}\left(\theta\right)$.
Then the elastic stiffness $\mathbb{C}(\phi_{\beta})$ and the gradient length
scale tensor $\boldsymbol{\Lambda{}}(\phi_{\beta})$ can be expressed as
$\mathbb{C}(\phi_{\beta})=\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}$ and
$\boldsymbol{\Lambda{}}(\phi_{\beta})=\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}$,
respectively. In the absence of inelastic strain, we have
$\boldsymbol{E}^{\text{in}}(\phi_{\beta})=\boldsymbol{0}$. For this case Eqs.
(31), (29) and (36) take the form

$\displaystyle\boldsymbol{S}-\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\operatorname*{grad}{\Bigl{(}\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}^{-1}\bigr{)}:\boldsymbol{S}\Bigr{)}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right]}=\boldsymbol{S}_{\text{t}}\,,$
(38a) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}-\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right]}\,,$
(38b) and $M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}-\frac{1}{2}\boldsymbol{E}:\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{E}-\frac{1}{2}\left(\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}*\boldsymbol{\Lambda{}}\bigr{)}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\\\
-\frac{1}{2}\left(\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(38c)

The isotropy of the gradient length scale tensor $\boldsymbol{\Lambda{}}$ for
cubic crystals implies that
$\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}=\boldsymbol{\Lambda{}}=\mathcal{l}^{2}\boldsymbol{I}$,
which simplifies Eqs. (38) to the following form

$\displaystyle\boldsymbol{S}-\mathcal{l}^{2}\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\operatorname*{grad}{\Bigl{(}\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}^{-1}\bigr{)}:\boldsymbol{S}\Bigr{)}}\right]}=\boldsymbol{S}_{\text{t}}\,,$
(39a) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}-\mathcal{l}^{2}\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\right]}\,,$
(39b) and $\displaystyle M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}-\frac{1}{2}\boldsymbol{E}:\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{E}-\frac{1}{2}\mathcal{l}^{2}\left(\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(39c)

####  3.3.3 Twin boundaries and boundaries between grains with inelastic
strain

Since the twin variants on both sides of the boundary are related by mirror
and/or rotational symmetry transformations between the unit cells, we can - as
in the case of grain boundaries - assume that the bulk chemical energy remains
unchanged, i.e.,
$\psi_{\text{b}}\left(\phi_{\beta}\,,\,\theta\right)=\psi_{\text{b}}\left(\theta\right)$,
and the elastic stiffness $\mathbb{C}(\phi_{\beta})$ and the gradient length
scale tensor $\boldsymbol{\Lambda{}}(\phi_{\beta})$ can be expressed using an
orthogonal tensor $\boldsymbol{Q}(\phi_{\beta})$ as
$\mathbb{C}(\phi_{\beta})=\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}$ and
$\boldsymbol{\Lambda{}}(\phi_{\beta})=\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}$,
respectively. Under these assumptions we find

$\displaystyle\boldsymbol{S}-\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\operatorname*{grad}{\Bigl{(}\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}^{-1}\bigr{)}:\boldsymbol{S}\Bigr{)}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right]}=\boldsymbol{S}_{\text{t}}+\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\bigl{(}\boldsymbol{E}^{\text{in}}(\phi_{\beta})\bigr{)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})}\right]}\,,$
(40a) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}_{\text{e}}(\boldsymbol{E})-\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right]}\,,$
(40b) and $M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\\\
-\frac{1}{2}\left(\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}*\boldsymbol{\Lambda{}}\bigr{)}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\frac{1}{2}\left(\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\cdot\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\boldsymbol{\Lambda{}}\bigr{)}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(40c)

For cubic lattices these expressions simplify to

$\displaystyle\boldsymbol{S}-\mathcal{l}^{2}\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\operatorname*{grad}{\Bigl{(}\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}^{-1}\bigr{)}:\boldsymbol{S}\Bigr{)}}\right]}=\boldsymbol{S}_{\text{t}}+\mathcal{l}^{2}\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\boldsymbol{E}^{\text{in}}(\phi_{\beta})}\right]}\,,$
(41a) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\mathcal{l}^{2}\operatorname*{div}{\left[\bigl{(}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\right]}\,,$
(41b) and $M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\\\
-\frac{1}{2}\mathcal{l}^{2}\left(\bigl{(}\partial_{\phi_{\beta}}\boldsymbol{Q}(\phi_{\beta})*\mathbb{C}\bigr{)}:\boldsymbol{\mathcal{Y}}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(41c)

####  3.3.4 Phase boundaries between cubic phases

In the case of phase boundaries between different cubic phases the gradient
length scale tensor $\boldsymbol{\Lambda{}}$ is isotropic on both sides of the
interface, even though not necessarily constant across the interface, i.e.,
$\boldsymbol{\Lambda{}}=\mathcal{l}(\phi_{\beta})^{2}\boldsymbol{I}$. This
allows us to reduce Eqs. (31), (29) and (36) to the following form

$\displaystyle\boldsymbol{S}-\operatorname*{div}{\left[\mathcal{l}(\phi_{\beta})^{2}\,\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\bigl{(}\mathbb{C}^{-1}(\phi_{\beta}):\boldsymbol{S}\bigr{)}}\right]}=\boldsymbol{S}_{\text{t}}+\operatorname*{div}{\bigl{(}\mathcal{l}(\phi_{\beta})^{2}\,\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\boldsymbol{E}^{\text{in}}(\phi_{\beta})}\bigr{)}}\,,$
(42a) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\operatorname*{div}{\left(\mathcal{l}(\phi_{\beta})^{2}\,\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\right)}\,,$
(42b) and $M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})+\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\mathcal{l}(\phi_{\beta})\partial_{\phi_{\beta}}\mathcal{l}(\phi_{\beta})\left(\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\\\
-\frac{\mathcal{l}(\phi_{\beta})^{2}}{2}\left(\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}-\rho\partial_{\phi_{\beta}}\psi_{\text{b}}(\phi_{\beta}\,,\,\theta)-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(42c)

##  4 Examples

To demonstrate the key properties of the above model numerical simulations
using the finite element method are performed using the commercial software
“COMSOL Multiphysics”222https://www.comsol.com/. A uniform mesh with
quadratic333Independent of the chosen shape functions, Comsol does not provide
third spatial derivatives of the degrees of freedom. Therefore, in order to
obtain the second derivative of strain (third spatial derivative of the
displacement), the “Distributed ODE” feature is used in order to introduce
additional degrees of freedom, corresponding to the second spatial derivatives
of the displacement. For this “Distributed ODE” linear shape functions are
employed., quadrilateral elements is used for the domain discretization. The
element size is 0.2 nm. Time stepping is performed using the BDF method. Based
on the assumptions of the small perturbation hypothesis444Both the
displacement $\boldsymbol{u}$ as well as the displacement gradient are
considered to be small, i.e., $|\boldsymbol{u}|\ll L$ and
$||\operatorname*{grad}{\boldsymbol{u}}||\ll 1$. (Maugin, 1992), we apply
traction boundary conditions to the undeformed geometry whenever required. We
assume elastostatics with an isotropic stiffness tensor
$\boldsymbol{\mathbb{C}}$. Material parameters have been chosen to represent
$\alpha$-iron with the values of the elastic constants, $E=200$ GPa,
$\nu=0.29$, and the Burger’s vector $b=0.285$ nm.

parameter name | symbol | value
---|---|---
Young’s modulus | $E$ | 200 GPa
Poisson’s ratio | $\nu$ | 0.29
Burger’s vector | $b$ | 0.285 nm
coefficient | $a$ | 2.98
coefficient | $A$ | $1.155\times 10^{8}\text{ J/m}^{3}$
coefficient | $B$ | $-3.43\times 10^{7}\text{ J/m}^{3}$
coefficient | $C$ | $-2.78\times 10^{8}\text{ J/m}^{3}$
mobility | $M$ | $2\text{ m}^{3}/\text{Js}$
gradient coefficient | $\alpha$ | $5\times 10^{-11}\text{ N}$
Table 1: Model parameters used for the numerical example in Sec. 4.2.

###  4.1 Regularization in the dislocation core

As shown in Sec. 3.3.1, the present model reduces to the set of equations
proposed by Po et al. (2018) in the homogeneous bulk phase. Here, we apply
this formulation to a single edge dislocation in an infinite elastic medium:
Fig. 1 shows the shear stress component $S_{12}$ in the plane perpendicular of
this dislocation with and without regularization ($\mathcal{l}=2$ Å). In the
“classical” case without regularization, the stress in the dislocation core is
singular, whereas it is well defined and finite for the regularized solution,
in analogy to what one would expect from a real atomistic configuration.

\begin{overpic}[percent,scale={0.95}]{disl_stress}
\put(62.0,15.0){\includegraphics[scale={0.6}]{sample_1dislocation.pdf}}
\end{overpic} (a) Shear stress component $S_{12}$ in the glide plane.

(b) Density plot the of stress component $S_{12}$.

Figure 1: Shear stress component $S_{12}$ for a single edge dislocation. The
inset in (a) shows the simulation setup.

###  4.2 Effect of the regularization on the interaction of dislocations with
a moving interface

This examples demonstrates the interaction of dislocations with a moving
interface between phase variant 1 (indicated by the superscript “V1”) and
variant 2 (indicated by the superscript “V2”). The phase mesostructure is
described by one single order parameter $\phi$. The only difference between
the two variants is with respect to the eigenstrain induced by the phase
transformation. This inelastic strain is given as a function of the order
parameter by

$\boldsymbol{E}^{\text{in}}(\phi_{\beta})(\phi)=\begin{cases}\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V1}}\,\varphi(\phi)&\textrm{if}\;\;\phi\geqslant
0\,,\\\
\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V2}}\,\varphi(\phi)&\textrm{if}\;\;\phi<0\\\
\end{cases}\,,$ (43)

where $\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V1}}$ and
$\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V2}}$ are the eigenstrains of
the phases V1 and V2, respectively,

$\displaystyle\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V1}}=\left(\begin{matrix}0&0.076\\\
0.076&0\end{matrix}\right)\,,$
$\displaystyle\boldsymbol{E}^{\text{in}}(\phi_{\beta})^{\text{V2}}=\left(\begin{matrix}0&-0.076\\\
-0.076&0\end{matrix}\right)\,.$ (44)

$\varphi(\phi)$ is a polynomial chosen in accordance with Levitas and Preston
(2002)

$\varphi(\phi)=\frac{a}{2}\phi^{2}+(3-a)\phi^{4}+\frac{1}{2}(a-4)\phi^{6}\,.$
(45)

The symmetric bulk chemical free energy takes the following form

$\rho\psi_{\text{b}}(\phi)=A\phi^{6}+B\phi^{4}+C\phi^{2}\,,$ (46)

and the interface energy density is assumed as

$\rho\psi_{\text{i}}\left(\operatorname*{grad}\phi\right)=\frac{\alpha}{2}\,||\operatorname*{grad}\phi||^{2}.$
(47)

For this specific case the resulting set of partial differential equations
(18) and (42a) can be further simplified to

$\displaystyle\operatorname*{div}{\boldsymbol{S}_{\text{t}}}=\boldsymbol{0}\,,$
(48a)
$\displaystyle\boldsymbol{S}-\mathcal{l}^{2}\Delta\boldsymbol{S}=\boldsymbol{S}_{\text{t}}+\mathcal{l}^{2}\,\mathbb{C}:\Delta\boldsymbol{E}^{\text{in}}(\phi_{\beta})\,,$
(48b) with
$\displaystyle\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi\bigr{)}=\mathbb{C}:\left(\boldsymbol{E}-\boldsymbol{E}^{\text{in}}(\phi_{\beta})\right)-\mathcal{l}^{2}\,\mathbb{C}:\operatorname*{div}{\boldsymbol{\mathcal{Y}}}\,,$
(48c) and $\displaystyle M^{-1}{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}=\alpha\Delta\phi+\boldsymbol{S}:\partial_{\phi}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\rho\partial_{\phi}\psi_{\text{b}}(\phi)\,.$
(48d)

These equations are solved for the displacement field $\boldsymbol{u}$, the
order parameter $\phi$ and the true stress $\boldsymbol{S}$. All parameters
and coefficients occurring in the above equations are summarized in Tab. 1.
The resulting interface energy, computed for a stationary flat interface, is
$\gamma=0.22\text{ J/m}^{2}$. The timescale in the simulation is controlled by
the mobility constant $M$. Since the simulation time can be arbitrarily re-
scaled using the mobility, in our simulations we treat it as dimensionless
pseudo time.

variant 2
---

variant 1
---

1
---
2
---
3
---
4
---
5
---
6
---
(a) Schematic representation of interface and dislocation arrangement. The
system is assumed to be periodic in vertical direction. Only the domain
indicated by the dashed box is used in the simulations. The false-color plot
indicated domains with positive (red) and negative (blue) in-plane shear
stress $S_{12}$.

(b) In-plane shear stress $S_{12}$ due to the dislocation for the two
different regularization lengths $\mathcal{l}$ used in this example.

Figure 2: Model problem with initially flat phase boundary driven by pure
shear loading towards a periodic arrangement of dislocations.  2
---
2.1
---
2.3
---
4.5
---
5
---
5.2
---
4
---
3.5
---
3.3
---
2.95
---
3.05
---
0
---
1
---
1.5
---
1.7
---
1.85
---
2
---
2.3 - 6
---
a)
---
b)
---
Figure 3: Propagation of the interface. The labeled lined correspond to the
center of the interface at the pseudo-times (in $\upmu$s) denoted by the
corresponding labels. a) $\mathcal{l}=0.6$ Å: The interface is arrested at the
dislocation array. b) $\mathcal{l}=1.9$ Å: The interface sweeps over the
dislocation array.

(a) Evolution of the V1 phase content for different regularization lengths
$\mathcal{l}$. For $\mathcal{l}=0.6$ Å the interface is arrested, whereas for
$\mathcal{l}=1.9$ Å it moves past the dislocation array.

(b) Rate of the evolution of the V1 phase content for different regularization
lengths $\mathcal{l}$. The large rate at time 5.2 is an artifact of
approaching the boundary of the simulation domain.

(c) Phase-boundary positions at equidistant time intervals for regularization
lengths $\mathcal{l}=0.6$ Å (left) and $\mathcal{l}=1.9$ Å (right). An
increasing distance between the contours indicates an acceleration of the
interface, while a decreasing distance indicates deceleration. The interface
positions at times 1, 2, 3, 4, 5, 6 $\upmu$s are shown in red (from left to
right).

(d) Evolution of the V1 phase content for $\mathcal{l}=0.6$ Å. A comparison
between the overall phase content, the phase content along the centerline of
the simulation box and the top of the simulation box.

(e) Evolution of the V1 phase content for $\mathcal{l}=1.9$ Å. A comparison
between the overall phase content, the phase content along the centerline of
the simulation box and the top of the simulation box.

Figure 4: Evolution of the phase content of variant 1 over pseudo-time.

The following scenario considers an initially flat interface between variants
V1 and V2, and a periodic, immobile dislocation structure with a dislocation
spacing of 10 nm within variant 2. For a pictorial representation, see Fig.
2a. The structure is assumed to be infinite in vertical direction, allowing us
to reduce the simulation domain to the dashed 40 nm wide and 10 nm high box in
Fig. 2a with periodic boundary conditions in vertical direction. The domain is
loaded under pure shear conditions with an in-plain shear stress of 85 MPa,
under which V1 is energetically more favorable, i.e., the interface will move
to the right.

Simulations are carried out using different regularization lengths
$\mathcal{l}$ ($\mathcal{l}=0.6$ Å and $\mathcal{l}=1.9$ Å) resulting in
different peak stresses in the dislocation core (see Fig. 2b). Fig. 3 shows
the positions of the V1-V2 interface for different points in simulation time.
As the interface approaches the dislocations it bows out due to the
interaction with the stress field of the dislocation core. The smaller
regularization length results in a larger stress magnitude in the dislocation
core region, leading to an arrest of the interface (see Fig. 3a). For the
larger regularization length the stress in the vicinity of the dislocation is
low enough in order to allow the interface to pass over the dislocation as
shown in Fig. 3b.

To analyze the temporal evolution in more detail Fig. 4 visualizes a number of
different aspects of the investigated system. The overall V1 phase content,
i.e., the area containing phase variant V1 divided by the whole area, as a
function of (pseudo) time is shown in Fig. 4a for the two different
regularization lengths. There, the most obvious characteristic is the arrest
of the interface for both values of $\mathcal{l}$ happening simultaneously
shortly after $t=2$ µs. While the system with the smaller regularization
length has already reached a stationary state, the other system shows that the
interface “detaches” from the dislocation and swipes the same area per time as
before, which shows in the same inclination of the respective line in Fig. 4a.

How does the _rate_ of the V1 phase content evolution change shortly before
and after the arresting of the interface took place? In Fig. 4a the peaks at
$t\approx 2$ and $\approx 3$ µs indicate that the phase boundary is
accelerating towards the dislocation until its velocity is significantly
reduced in the vicinity of the dislocation. The second peak shows that the
interface effectively accelerates again after passing the dislocation. This
stage is followed by another dip (in between $\approx 3$ and $\approx 3.5$ µs)
where the interface motion right of the dislocation is again decelerated. This
behaviors is also visualized in Fig. 4c, which shows the interface position at
equidistant points in time.

Two different phenomenons operate here, which can be understood from the
change of sign of the shear stress field of an edge dislocation as shown in
Fig. 2a. Recall that the inelastic strain of the interface is governed only by
the shear components of the strain tensor. Once the phase boundary is getting
close enough to interact with the dislocation, the upper and lower sections of
phase V2 are in the regions 3/6 of the dislocation (compare Fig. 2a). The
driving force is effectively directed in positive x-direction and causes the
acceleration. The central regions of the phase boundary, is located in region
4 of the dislocation stress field and therefore experiences a net driving
force that is directed in _opposite_ direction. This interplay between the
directions of the two driving forces is also responsible for the curvature of
the interface. Once the interface has passed the dislocation, the central
regions of the interface experiences a very large driving force in positive
direction (region 1 of the dislocation). The top and bottom sections of the
interface, however, are located in regions with negative driving force 2/6.
When the interface moves further towards the right, the magnitude of the
driving force from the dislocation acting in section 1 decreases as $1/r$ and
this region of the interface decelerates. At the same time, the driving force
acting on the top and bottom of the interface increases only slightly,
explaining the second dip in Fig 4b around $t=3...3.75$µs. At a sufficient
distance from the dislocation the top and bottom parts of the interface
accelerate, which leads to a decrease in curvature of the interface.

This is further shown in Figs. 4d and 4e, which relate the motion of the
curved interface to the motion of a flat interface for the case without a
dislocation structure (indicated by the dotted line).

##  5 Summary

In this paper we developed a framework for coupling a phase-field description
of planar defects such as phase or twin-boundaries with a discrete
representation of dislocations within (anisotropic) first-strain-gradient
elasticity. Its main features and advantages in contrast to phase-field within
classical elasticity are:

  * 1.

Non-singular stresses at the dislocation core that can be easily calibrated to
match molecular statics predictions using the approach of Admal et al. (2017)

  * 2.

Non-singular driving forces for the evolution of the phase-field evolution in
the presence of dislocation. This ensures a mesh-independent numerical
solution and is a necessary condition for modeling the _interaction_ of
dislocations with interfaces such as phase-, grain- or twin-boundaries.

We have shown that in order to ensure regularized driving forces in the
dislocation core, a Helmholtz-type elastic free energy that is quadratic in
the gradient of the total rather than the elastic strain must be used.

We implemented the proposed framework in the Comsol Multiphics Modeling
Software and demonstrated its feasibility and basic properties based on a
number of examples. Coupled to a dislocation-dynamics code, we expect this
phase-field framework to be a valuable tool for understanding microstructure-
evolution on a small scale.

Acknowledgements
The authors gratefully acknowledge the Deutsche Forschungsgemeinschaft (DFG)
for supporting this work carried out within the framework of Collaborative
Research Center SFB 799. SS acknowledges financial support from the European
Research Council through the ERC Grant Agreement No. 759419 MuDiLingo (“A
Multiscale Dislocation Language for Data-Driven Materials Science”).

##  Appendix A Energy that is quadratic in
$\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)$

Starting with a setup identical to Sec. 3.2 but for the Helmholtz-type elastic
free energy

$\rho\psi_{\text{e}}\left(\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\,,\,\theta\right)=\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)+\frac{1}{2}\left(\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\,,$
(49)

we find

$\displaystyle\boldsymbol{S}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)=\mathbb{C}(\phi_{\beta}):\left(\boldsymbol{E}-\boldsymbol{E}^{\text{in}}(\phi_{\beta})\right)\,,$
(50)
$\displaystyle\boldsymbol{\mathcal{T}}=\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\,,$
(51)
$\displaystyle\boldsymbol{\xi}_{\beta}=\alpha\operatorname*{grad}\phi_{\beta}-\boldsymbol{\mathcal{T}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})\,,$
(52)

and once again combining the first two equations

$\displaystyle\boldsymbol{\mathcal{T}}=\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\left(\mathbb{C}^{-1}(\phi_{\beta}):\boldsymbol{S}\right)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\,.$
(53)

From Eq. (17) we find the constitutive equation for the total stress
$\boldsymbol{S}_{\text{t}}$

$\boldsymbol{S}_{\text{t}}\bigl{(}\boldsymbol{E}\,,\,\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\bigr{)}=\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\bigl{(}\boldsymbol{\mathcal{Y}}-\operatorname*{grad}{\boldsymbol{E}^{\text{in}}(\phi_{\beta})}\bigr{)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right]}\,.$
(54)

the equation to determine true stress $\boldsymbol{S}$ from the total stress
$\boldsymbol{S}_{\text{t}}$

$\boldsymbol{S}-\operatorname*{div}{\left[\mathbb{C}(\phi_{\beta}):\operatorname*{grad}{\left(\mathbb{C}^{-1}(\phi_{\beta}):\boldsymbol{S}\right)}\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right]}=\boldsymbol{S}_{\text{t}}\,.$
(55)

The evolution equation for the order parameter obtained using the same
procedure as in Sec. 3.2 is

$M{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\operatorname*{div}{\left[\boldsymbol{\mathcal{T}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})\right]}-\rho\partial_{\phi_{\beta}}\psi\,.$
(56)

The divergence on the right hand side of (56) is easily evaluated:

$\displaystyle\operatorname*{div}{\left[\boldsymbol{\mathcal{T}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})\right]}$
$\displaystyle=-\operatorname*{div}{\left[\boldsymbol{\mathcal{T}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)\right]}$
$\displaystyle=-\operatorname*{div}{\boldsymbol{\mathcal{T}}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\boldsymbol{\mathcal{T}}\,\smash{\vdots}\,\partial_{\phi_{\beta}}\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)$
$\displaystyle=\boldsymbol{S}_{\text{t}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\boldsymbol{S}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\boldsymbol{\mathcal{T}}\,\smash{\vdots}\,\partial_{\phi_{\beta}}\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)$
$\displaystyle=\boldsymbol{S}_{\text{t}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})+\rho\partial_{\phi_{\beta}}\psi-\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-$
$\displaystyle\hskip
71.13188pt\frac{1}{2}\left(\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\partial_{\phi_{\beta}}\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)-$
$\displaystyle\hskip
99.58464pt\frac{1}{2}\left(\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)$

Finally, we find the expression

$M{\dot{\phi\mkern
5.0mu}\mkern-5.0mu}{}_{\beta}=\alpha\Delta\phi_{\beta}+\boldsymbol{S}_{\text{t}}:\partial_{\phi_{\beta}}\boldsymbol{E}^{\text{in}}(\phi_{\beta})-\frac{1}{2}\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right):\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{E}^{\text{e}}\left(\boldsymbol{E}\,,\,\phi_{\beta}\right)-\\\
\quad\frac{1}{2}\left(\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\partial_{\phi_{\beta}}\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)-\frac{1}{2}\left(\partial_{\phi_{\beta}}\mathbb{C}(\phi_{\beta}):\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)\cdot\boldsymbol{\Lambda{}}(\phi_{\beta})\right)\,\smash{\vdots}\,\boldsymbol{\mathcal{Y}}^{\text{e}}\left(\boldsymbol{\mathcal{Y}}\,,\,\phi_{\beta}\right)-\\\
\rho\partial_{\phi_{\beta}}\psi_{\text{b}}(\phi_{\beta}\,,\,\theta)-\partial_{\phi_{\beta}}g(\phi_{\beta}\,,\,\theta)\,.$
(57)

where the total stress $\boldsymbol{S}_{\text{t}}$ appears in the driving
force. In general, this stress cannot be assumed to be bounded in the
dislocation-core. This is illustrated in Fig. 5 that shows the maximum shear
stress in the dislocation core for different “thicknesses” of the dislocation,
i.e., different discretizations. While the true stress $S_{12}$ does not
change noticeably once the discretization is sufficiently fine, the total
stress $S_{\text{t}12}$ keeps increasing with decreasing thickness of the
dislocation.

Figure 5: The maximum shear stress in the dislocation core as a function of
the dislocation “thickness”.

## References

  * Admal et al. (2017) Admal, N.C., Marian, J., Po, G., 2017.  The atomistic representation of first strain-gradient elastic tensors.  Journal of the Mechanics and Physics of Solids 99, 93–115.  doi:10.1016/j.jmps.2016.11.005.
  * Cai et al. (2006) Cai, W., Arsenlis, A., Weinberger, C.R., Bulatov, V.V., 2006\.  A non-singular continuum theory of dislocations.  Journal of the Mechanics and Physics of Solids 54, 561–587.  doi:10.1016/j.jmps.2005.09.005.
  * Del Piero (2009) Del Piero, G., 2009.  On the method of virtual power in continuum mechanics.  Journal of Mechanics of Materials and Structures 4, 281–292.  doi:10.2140/jomms.2009.4.281.
  * Germain (1973) Germain, P., 1973.  The Method of Virtual Power in Continuum Mechanics. Part 2: Microstructure.  SIAM Journal on Applied Mathematics 25, 556–575.  doi:10.1137/0125053.
  * Gurtin (1996) Gurtin, M., 1996.  Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance.  Physica D: Nonlinear Phenomena 92, 178–192.
  * Lazar et al. (2005) Lazar, M., Maugin, G.A., Aifantis, E.C., 2005.  On dislocations in a special class of generalized elasticity.  physica status solidi (b) 242, 2365–2390.  doi:10.1002/pssb.200540078.
  * Lazar et al. (2006) Lazar, M., Maugin, G.A., Aifantis, E.C., 2006.  Dislocations in second strain gradient elasticity.  International Journal of Solids and Structures 43, 1787–1817.  doi:10.1016/j.ijsolstr.2005.07.005.
  * Lazar and Po (2015) Lazar, M., Po, G., 2015.  The non-singular Green tensor of Mindlin’s anisotropic gradient elasticity with separable weak non-locality.  Physics Letters A 379, 1538–1543.  doi:10.1016/j.physleta.2015.03.027.
  * Levitas and Preston (2002) Levitas, V.I., Preston, D.L., 2002\.  Three-dimensional Landau theory for multivariant stress-induced martensitic phase transformations. II. Multivariant phase transformations and stress space analysis.  Physical Review B 66, 134207\.  doi:10.1103/PhysRevB.66.134207.
  * Lothe (1992) Lothe, J., 1992.  Dislocations in Continuous Elastic Media, in: Indenbom, V.L., Lothe, J. (Eds.), Elastic Strain Fields and Dislocation Mobility. Elsevier. volume 31 of Modern Problems in Condensed Matter Sciences, pp. 175–235.  doi:10.1016/B978-0-444-88773-3.50008-X.
  * Maugin (1980) Maugin, G., 1980.  The method of virtual power in continuum mechanics: Application to coupled fields.  Acta Mechanica 35, 1–70.
  * Maugin (1992) Maugin, G., 1992.  The Thermomechanics of Plasticity and Fracture.  Cambridge University Press.
  * Mindlin (1964) Mindlin, R., 1964.  Micro-structure in linear elasticity.  Archive for Rational Mechanics and Analysis 16, 51–78.
  * Nabarro (1947) Nabarro, F.R.N., 1947.  Dislocations in a simple cubic lattice.  Proceedings of the Physical Society 59, 256.  doi:10.1088/0959-5309/59/2/309.
  * Peierls (1940) Peierls, R., 1940.  The size of a dislocation.  Proceedings of the Physical Society 52, 34.  doi:10.1088/0959-5309/52/1/305.
  * Po et al. (2018) Po, G., Lazar, M., Admal, N.C., Ghoniem, N., 2018.  A non-singular theory of dislocations in anisotropic crystals.  International Journal of Plasticity 103, 1–22.  doi:10.1016/j.ijplas.2017.10.003.


# Generator coordinate method for transition-state dynamics in nuclear fission

G.F. Bertsch  Department of Physics and Institute for Nuclear Theory, Box
351560, University of Washington, Seattle, Washington 98915, USA  K. Hagino
Department of Physics, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

The existence of transition channels across the fission barrier has been
central to the theory of induced fission, but there has been no microscopic
theory applicable to energies at the barrier top. We propose a microscopic
treatment motivated by the Generator Coordinate Method (GCM) and Gaussian
Overlap Approximation (GOA) to parameterize both the dynamics within the
channels and their incoherent couplings to states outside the barrier. The
physical characteristics of the channels (often called ”transition states”)
examined here are their effective bandwidths for crossing the barrier and the
quality of the coupling to compound-nucleus states as measured by the
transmission factor $T$. We also investigate the spacing of GCM states with
respect to their degree of overlap. We find that a rather coarse mesh provides
an acceptable accuracy for estimating the bandwidths and transmission factors.
The common numerical stability problem in using the GCM is avoided due to the
choice of meshes and the finite bandwidths of the channels.

##  I Introduction

Transition-state theory111 The term “channel” describes its role in reaction
theory better than “state” and we shall use that designation for the models
presented here. has been at the foundation of the theory of induced fission
since the original paper by Bohr and Wheeler in 1939 bo39  and continuing up
to the present era bo13 ; capote2009 ; chadwick2006 ; chadwick2011 ; cap2011 ;
lu2016 ; schmidt1991 . It is encapsulated in the formula for the decay rate
$\Gamma$

$\Gamma_{\rm BW}=\frac{1}{2\pi\rho}\sum_{i}T_{i}$ (1)

where $i$ labels channels, $\rho$ is the level density of the compound
nucleus, and $T_{i}$ is a transmission coefficient or conductance. It is also
identical to the penetration factor in subbarrier conductance. It satisfies
the bounds $0\leq T\leq 1$. Typically $T$ is assumed to depend on energy as a
particle traversing an inverted parabolic barrier, but that is a pure guess
absent a microscopic understanding of the Hamiltonian dynamics.

The goal of this paper is to carry out the first steps of building a
microscopic theory of the barrier-crossing dynamics applicable to heavy
nuclei. In the theory of large nuclei, the starting point is the wave
functions of self-consistent mean-field theory, such as those given by the
energy density functionals of Skyrme, Gogny, or relativistic formulations
bender03 . Besides the self-consistent solutions of the Hartree-Fock (HF) or
Hartree-Fock-Bogoliubov (HFB) equations, an adequate basis of states for
studying transport properties can be constructed using the Generator
Coordinate Method (GCM). This requires the calculation of mean-field
configurations that are constrained by one or more single-particle fields. The
GCM has been used previously for modeling fission dynamics near the barrier
top go05 ; ta17 . In those works, the authors used GCM with two constraining
fields and the Gaussian Overlap Approximation (GOA) to map the Hamiltonian
onto a two-dimensional Schrödinger equation. However, the steps needed to
arrive at a Schrödinger equation ignore the statistical aspects of the decay
and gives no hint of a connection to Eq. (1). In an earlier paper be21  we
showed how one could derive the transition-state formula in a highly
simplified configuration-interaction approach. Here we shall use the same
reaction theory formalism to calculate transmission coefficients, but with a
much more realistic description of the channels. An important advantage of the
reaction theory is that statistical aspects of the theory can be easily
included in the formalism ha21 .

A technical obstacle in the GCM approach is the nonorthogonality of the basis
configurations. As will be shown, any formal difficulties are avoided in
present-day theory based on the many-body Green’s function in Eq. (4). A
related problem is the danger of numerical instabilities when overlaps between
configurations are large. We will show that for the reaction theory that
problem does not arise, because one can use coarse bases without much loss of
accuracy.

For investigating transition channels in fermionic systems, the general
characteristics can be derived independently of the details of the
constraining field. A configuration is labeled by the expectation value of the
field; we shall call the expectation value $q$ and further specify it as
$q_{i}$ for a configuration $i$ in a finite-dimensional basis. We will start
the analysis with a model that assumes the validity of the Gaussian Overlap
Approximation.

Besides the internal properties of the channel, one needs specific information
about the coupling of the reservoirs of states on either side of the channel.
The situation is very similar to treating couplings to electrical cables. The
cable has a characteristic impedance, and conductance depends on impedance
matching. An optimally matched coupling yields a transition conductance $T=1$.
Mismatches on either side decrease it, except for resonances within the cable.
To model the complete conductance we define reservoirs of states outside the
channel and some details of the interaction connecting them to channel states.
In this paper we treat the outside states and their couplings schematically.

##  II GCM methodology for transmission channels

The usual procedure for applying the GCM to nuclear spectroscopy consists of
the following steps.

1) Define a set of configurations calculated in mean-field theory and
constrained by some physical one-body field such as the mass quadrupole moment
$Q$. The set of expectation values of the field $(q_{1},q_{2},...q_{N})$
defines an $N$-dimensional basis for the configuration space.

2) Calculate the matrix $\boldsymbol{N}$ of overlaps between configurations
and the matrix $\boldsymbol{H}$ of the Hamiltonian or the energy functional
that plays the role of the Hamiltonian in the mean-field theory. Here and
below we use boldface symbols for matrices.

3) Solve the nonhermitian eigenvalue problem (i.e., the Hill-Wheeler equation)

$\boldsymbol{H}\psi=E\boldsymbol{N}\psi$ (2)

for energies $E$ and corresponding $N$-dimensional wave functions $\psi$.

4) Check for convergence by varying the number of configurations $N$ in the
calculation. The effect on the properties in the low-energy part of the
spectrum should be small.

Steps 1) and 2) are the same for calculating reaction rates in the GCM, but
the remaining steps are completely different. Namely, the new steps are:

3’.) In a new step 3), the Hamiltonian is made complex by adding imaginary
terms $-\boldsymbol{\Gamma}_{j}/2$

$\boldsymbol{H}^{\prime}=\boldsymbol{H}-i\sum_{j}\boldsymbol{\Gamma}_{j}/2.$
(3)

Each $\boldsymbol{\Gamma}_{j}$ is a matrix of decay rates to states $j$
outside of the model space. For the present problem, there are two decay
modes, one corresponding to the set of compound nucleus states and the other
to states in the second well. We label $a$ and $b$, respectively, in the
equation below.

4’.) The next step is to calculate the Green’s function for the Hamiltonian
$\boldsymbol{H}^{\prime}$ by the matrix inversion,

$\boldsymbol{G}(E)=(\boldsymbol{H}-i\boldsymbol{\Gamma}_{a}/2-i\boldsymbol{\Gamma}_{b}/2-\boldsymbol{N}E)^{-1}.$
(4)

This replaces the matrix diagonalization in the old step 3).

5’.) As a final step, the transmission factor $T_{ab}$ between reservoir $a$
and $b$ is computed in $S$-matrix theory as encapsulated by the Datta formula
da95 ; al21

$T_{ab}=\sum_{ijkl}(\boldsymbol{\Gamma}_{a})_{ij}\boldsymbol{G}_{jk}(\boldsymbol{\Gamma}_{b})_{kl}\boldsymbol{G}^{*}_{li}={\rm
Tr}\,\left(\boldsymbol{\Gamma}_{a}\boldsymbol{G}\boldsymbol{\Gamma}_{b}\boldsymbol{G}^{*}\right).$
(5)

Calculated this way, $T_{ab}$ is a continuous real function of $E$ in the
range $[0,1]$. As in the procedure for spectroscopic studies, one gains
confidence by varying the dimension of the configuration spaces.

##  III Decay widths

Constructing the matrices $\boldsymbol{H}$ and $\boldsymbol{N}$ is
straightforward with available tools for calculating properties of large
nuclei in the HF or HFB approximations. The calculation of the internal widths
$\boldsymbol{\Gamma}_{j}$ is more subtle. The guiding principle is Fermi’s
Golden Rule for the decay of configuration $i$ into a set of states $j$:

$\Gamma_{j}=2\pi\overline{\langle i|H|j\rangle^{2}}\rho_{j}.$ (6)

Here $\rho_{j}$ is the density of final states, and the overline indicates an
average over them. The state $i$ is part of the basis for the transmission
channel, while the states $j$ might include quasiparticle excitations of the
GCM configurations beyond the channel entry and exit points. The Golden Rule
assumes that the states are defined in an orthogonal basis, even though it can
be extended to a non-orthogonal basis bertsch19 . It is easy to keep $i$
nearly orthogonal to the set $j$ by demanding that the orbital fillings be
different.

The transmission channel includes multiple configurations $i$ and $i^{\prime}$
and the coherence between them needs to be taken into account. This requires a
matrix form for $\boldsymbol{\Gamma}$,

$(\boldsymbol{\Gamma}_{j})_{ii^{\prime}}=2\pi\overline{\langle
i|H|j\rangle\langle i^{\prime}|H|j\rangle}\rho_{j}.$ (7)

We follow this structure and use the separable form

$(\boldsymbol{\Gamma}_{j})_{ii^{\prime}}=\gamma_{j}g_{i_{j}}(i)g_{i_{j}}(i^{\prime}).$
(8)

to parameterize them. Here $i_{j}$ labels e.g., the cm positions of the final
state configurations. We assume that the decays take place near the positions
of the end states, so that $i_{e}=1$ or $N$.

##  IV Examples of GOA Hamiltonians

Before specializing to the internal transmission channels that are the subject
of this paper, we show how the GCM/GOA methodology can be applied to calculate
the free propagation of a composite particle. A basic assumption is that the
mean-field wave function $\Psi_{q_{i}}$ of the many-body system can be
factorized into a product of an internal wave function $\Psi(\xi)$ and a
Gaussian wave packet in the center-of-mass coordinate $x$ be19c ,

$\Psi_{z}(x)\sim\Psi(\xi)\exp\left(-(x-z)^{2}/2s^{2}\right).$ (9)

Here $z$ locates the center of the wave packet and $s$ is its r.m.s. width.
The factorization is exact for the ground states in the harmonic oscillator
model of the single-particle Hamiltonian, but is an uncontrolled approximation
for more realistic mean-field potentials. On the other hand, the Gaussian form
of the overlap function in Eq. (10) below can be checked and is quite well
satisfied in practice be19c . With this approximations the overlap matrix
$\boldsymbol{N}=\langle\Psi_{z}|\Psi_{x^{\prime}}\rangle$ is given by

$\boldsymbol{N}_{ij}=\exp(-(z-z^{\prime})^{2}/4s^{2}).$ (10)

Under the GOA the Hamiltonian matrix $\boldsymbol{H}$ can be approximated as a
quadratic polynomial times the overlap function,

$\boldsymbol{H}_{ij}=E_{K}(1-(z-z^{\prime})^{2}/2s^{2})\boldsymbol{N}_{ij}.$
(11)

The prefactor $E_{K}$ is the kinetic energy associated with the center-of-mass
wave function, $E_{K}=\hbar^{2}/4Ms^{2}$ with $M$ the mass of the composite
particle. When the formula is applied to other GCM collective coordinates,
$E_{K}$ is related to the zero-point energy associated with that coordinate.
Eq. (11) depends only on the coordinate difference $z-z^{\prime}$ when the
Hamiltonian is translationally invariant.

Translational invariance also implies that the eigenfunctions of the center-
of-mass coordinate are plane waves. A plane wave of momentum $k$ is
constructed by the integral

$\Psi_{k}=\int^{\infty}_{-\infty}dz\,\Psi_{z}e^{ikz}.$ (12)

The corresponding eigenvalues can be computed as the ratio of matrix elements

$E_{k}=\frac{\langle\Psi_{z}|H|\Psi_{k}\rangle}{\langle\Psi_{z}|\Psi_{k}\rangle}=\frac{\hbar^{2}}{2M}k^{2}.$
(13)

Here $z$ is arbitrary.

Turning now to the computational treatment of internal channels, we consider a
chain of $N$ states spanning an interval $[q_{1},q_{N}]$ with spacing $\Delta
q=(q_{N}-q_{1})/N$ between configurations. As discussed above, imaginary decay
matrices centered at end points $q_{1}$ and $q_{N}$ are added to
$\boldsymbol{H}$. If the cm wave functions of the reservoir states have the
same Gaussian form, $g$ in Eq. (8) may be taken as
$g_{i_{e}}(i)=\boldsymbol{N}_{i_{e}i}$. The resulting Hamiltonian is then

$\boldsymbol{H}^{\prime}_{ij}=\boldsymbol{H}_{ij}-i\gamma\boldsymbol{N}_{i1}\boldsymbol{N}_{j1}-i\gamma\boldsymbol{N}_{iN}\boldsymbol{N}_{jN}.$
(14)

###  IV.1 A single flat channel

The first model we investigate is a flat chain composed of $N=4$
configurations with overlaps between them determined by the parameter
$s=1/\sqrt{5}$. This choice of $s$ was shown in Ref. be19c  to give a good
compromise between accuracy and computational effort. The channel is the
topmost one depicted in Fig. 1. The states indicated by black circles are the
ones included in the $N=4$ model. We will also examine the same model with 7
states; the added states are shown as the red circles.

Figure 1: Relationship between states in the models described in Sections
IV.A,B, and C. The states in the 4-state and 7-state channels are shown as
black and black+red circles respectively. The real part of the Hamiltonian
couples the states in the channel or channels; the couplings to the reservoirs
are parameterized by the imaginary part of the Hamiltonian.

The diagonal energies of the GCM states $E_{K}$ are taken to be $E_{K}=5/4$
and the strengths of the absorption at the ends are $\gamma=1$. With these
parameters the overlap between neighboring states is fairly small,
$\boldsymbol{N}_{i,i+1}=0.28$. The resulting transmission factor $T(E)$
calculated by Eq. (5) is shown in Fig. 2 as the black solid line.

Figure 2: Transmission factor for a chain of length $q=3\sqrt{5}/2$ comparing
GCM calculations for 4 and 7 states in the chain (solid black and dashed red
lines, respectively). The parameters of the Hamiltonian are
$(s,E_{0},\gamma)=(1/\sqrt{5},5/4,1)$. See the Supplementary Material for the
computer scripts used to calculate the data presented here and later in the
Figures.

One sees a structure of 4 peaks, each close to an eigenvalue of the Hill-
Wheeler Eq. (2). Physically, the peaked structure arises from the wave
reflection at the ends of the channel. Note also that the range of $T(E)$
satisfies $0<T<1$ as required by the unitarity of the $S$ matrix. Note also
that the channel starts conducting near $E\approx 0$, as would be the case for
a classical channel. The adequacy of the mesh spacing can be assessed by
shrinking it. Decreasing it by a factor of 2, the same interval contains 7
seven states instead of 4. The resulting transmission factor is shown as the
dashed red line in Fig. 2. One sees that in the low-energy region it is quite
similar to the 4-state approximation. However, it has 3 additional peaks at
higher energy, corresponding to the high-energy eigenfunctions of the 7-d
model. These peaks are much narrower than the lower ones and can be neglected
in calculating integrated transmission rates. The same behavior would continue
with finer mesh spacings; there would remain 4 peaks in the energy region
$[0,2]$ and the additional narrow peaks would appear at higher and higher
energies. The qualitative aspects of this behavior can be easily understood.
With a finite mesh spacing of Gaussian wave packets one can approximate plane
wave with a good fidelity for low momentum, but there is a momentum cutoff
controlled by the mesh spacing. In the transmission channel as parameterized,
the momentum at the injection and exit point is controlled by the Gaussian
width parameter $s$. The momentum match to the channel parameters suppresses
the transmission to the high-momentum modes in the channel. We conclude that
fairly sparse meshes are adequate for representing the overall conductivity of
flat transmission channels.

As mentioned earlier, very fine mesh spacings often lead to numerical
instabilities in the spectroscopic applications of the GCM. The usual fix is
to make a singular value decomposition of the overlap matrix, throwing out
eigenfunctions that have small norms. It is instructive to see what happens
when the same procedure is applied here.

Figure 3: Transmission factor for a chain of length $q=4$ comparing GCM
calculations for 4 and 7 states in the chain (solid black and dashed red
lines, respectively). The difference with Fig. 2 is that the 7-d space was
truncated to 4 dimensions by the singular value decomposition of the overlap
matrix. The Hamiltonians are the same as in Fig. 2.

Fig. 3 compares the $4$-state model with the 7-state model truncated to 4
states by eliminate the small-norm vectors. One sees that the resonance
positions are rather close and the widths are also very similar. There is no
obvious benefit from starting out with a larger space. Since there is no need
to truncate the space for reasons of numerical stability, this aspect of the
usual methodology can be dropped.

We next examine how $T(E)$ depends on the strength of the absorption at the
ends of the channel. Fig. 4 shows $T(E)$ for a range of absorption strengths
$\gamma$.

Figure 4: Transmission factors in the 4-d model for several values of
absorption strength: $\gamma=0.5$ (solid black line); $\gamma=1.0$ (red dashed
line); $\gamma=2.0$ (blue dotted line).

Obviously, for small $\gamma$ the channel acts as a resonant cavity with
sharply defined resonances and the overall conductance is low. For the
intermediate value the reflection amplitude is small and the individual peak
structure disappears. Very large $\gamma$’s are probably unphysical, as the
couplings are strong enough to concentrate the decay strength into a small
energy region.

###  IV.2 A parabolic channel

In this section we extend the model to include a peaked barrier. We take the
shape of the barrier as an inverted parabola, as is often assumed in
phenomenological treatments.

Under the factorization Ansatz Eq. (9) the GCM matrix elements of a potential
depending only on the cm position are given by

$\displaystyle\langle i|V|j\rangle=$
$\displaystyle\frac{1}{s\pi^{1/2}}\int^{\infty}_{-\infty}dx\,V(x)e^{-(x-z_{i})^{2}/2s^{2}-(x-z_{j})^{2}/2s^{2}}.$
(15)

Here $V(x)$ is taken as the parabolic form $V(x)=V_{2}(x-x_{b})^{2}$ where
$x_{b}$ is at the center of the barrier. The resulting GCM matrix elements are

$\langle
i|V|j\rangle=V_{2}\left[\left(\frac{z_{i}+z_{j}}{2}-x_{b}\right)^{2}+\frac{s^{2}}{2}\right]\boldsymbol{N}_{ij}.$
(16)

The matrix $\boldsymbol{V}$ of these elements are added to the Hamiltonian
defined in Eq. (11) and (14). Note that the diagonal potential matrix elements
are slightly below the defining potential due to the second term in Eq. (16).
The diagonal energies are indicated in the channel marked “B” in Fig. 1. For a
numerical example we take $V_{2}=-1/2$. The channel Hamiltonian then has
eigenenergies ranging from $-$0.4 to 3.3. Fig. 5 shows the transmission factor
as a function of energy taking $\gamma=1.0$.

Figure 5: Transmission factor in the 4-d and 7-d models with a parabolic
barrier. Solid black line: 4x4 model; red dashed line: 7x7 model.

It may be seen that the coupling to outside is strongest at the energy of the
end channel states. However, at that energy the barrier suppresses the
conductance and $T$ is a fraction of the maximum. At higher energies the $T$
can approach its maximum value of one, but the coupling is weaker and the
peaks are narrow.

We believe this behavior is generic for channels that follow the topography of
the potential energy surface. This is the case when they are constructed using
an adiabatic approximation. To see that the results are not an artifact of the
GCM mesh spacing, we also show the transmission factor taking a finer mesh
with 7 GCM configurations instead of 4. One sees that the low-energy
conductance is almost the same. At higher energies, the narrowing of the peaks
is also similar, although the peaks are somewhat shifted in energy.

###  IV.3 Two crossing channels

To understand better the adiabatic treatment, we consider a model in which the
adiabatic channels arise from coupling between diabatic ones. We start with
two diabatic channels that cross as depicted in the lowest diagram in Fig. 1.
The dashed black lines link configurations that have large matrix elements in
HF mean field theory; resulting chains are the diabatic paths in the dynamics.
Adiabatic dynamics arise when one first diagonalizes the Hamiltonian within
the subspace at fixed $z_{i}$. These are indicated by the curved red dotted
lines in the Figure. The picture of adiabatic channels peaking at the barrier
top is unavoidable in transition-channel theory as implemented in Eq. (1).

For the Hamiltonian model we add linear potentials to generate the diabatic
paths together with a constant interaction between intrinsic states at the
same positions $q$. The matrix elements for a potential having a constant
slope $V_{a}(x)=v_{a}x$ are given by

$\boldsymbol{V}_{ai;a^{\prime}j}=v_{a}\left(\frac{z_{i}+z_{j}}{2}-x_{b}\right)\boldsymbol{N}_{ij}\delta_{aa^{\prime}}.$
(17)

Here $a$ and $a^{\prime}$ label the two diabatic channels. The other term to
be added to the Hamiltonian is the coupling $h_{c}$ between intrinsic states
of the two diabatic channels. We take the coupling as a constant independent
of $z$. Again invoking factorization hypothesis, the matrix elements are

$\boldsymbol{H}_{ai,a^{\prime}j}=h_{c}(1-\delta_{a,a^{\prime}})N_{ij}.$ (18)

For the numerical example, we take $v_{a}=\pm 0.5$ in the channels and
$h_{c}=0.8$ for the coupling.

As depicted in Fig. 1 there are now 4 decay matrices to be added to the
Hamiltonian. We may assume that final states are all orthogonal, so we can
apply the transmission formula Eq. (5) with an incoherent sum over all four
combinations $(a,a^{\prime})\times(b,b^{\prime})$.

In the adiabatic approximation, only the transmission factor from the two
lowest states at the ends are included,

$T_{\rm adiabatic}\approx T_{ab^{\prime}}.$ (19)

It is shown as the red dashed line in Fig. 6. The dotted blue line shows the
combined transmission factor that includes including the upper adiabatic
channel as well,

$T^{\prime}_{\rm adiabatic}\approx T_{a^{\prime}b}.$ (20)

These are to be compared to the full transmission factor (solid black line)
including all contributions,

$T=T_{ab}+T_{a^{\prime}b^{\prime}}+T_{ab^{\prime}}+T_{b^{\prime}a}.$ (21)

Figure 6: Transmission factor for the Hamiltonian “C8” depicted in Fig. 1.
Solid black line: all contributions to $T$; red line: lowest adiabatic
contribution; blue line: both adiabatic contributions

One sees that the adiabatic approximation works well overall when both
channels are included. The second channel adds hardly anything where the lower
channel is open, but fills out the higher region $1.0<E<2.5$. Another
interesting finding, not very visible in the figure, is that the adiabatic
approximation significantly under predicts the transmission factor at the
lowest energies. This inadequacy of the approximation was noted earlier in
Refs. ha20 ; gi14 .

We have also calculated $T$ without any coupling between the diabatic
channels. As expected, that treatment seriously underpredicts the transmission
coefficient.

##  V General conclusions

.

A few tentative conclusions may be drawn from the simple models presented
here. First of all, one does not need fine collective-coordinate meshes in the
GCM configuration space. A mesh spacing giving overlaps of 0.3 between a
configuration and its diabatic neighbor seems adequate; smaller mesh spacings
will place the resonances at different positions but the coarse properties of
the channel will remain the same. The second conclusion is that momentum
matching is an important consideration in the channel coupling to the
reservoir states. It produces an effective energy cutoff in the conductance of
the channel. The energy scale for this effect is given by the zero-point
energy of the collective coordinate in the mean-field wave function. To give a
sense of that, we present in Table I some characteristics of the transmission
function $T(E)$ for the models discussed in the previous section. The first
characteristic is the integrated transmission factor. This is reported in the
Table in units of $E_{K}$,

$I_{T}=\int_{-\infty}^{\infty}dE\,T(E)/E_{K}.$ (22)

Another important characteristic is the energy difference between the nominal
barrier height and the energies of the states in the reservoirs that are
strongly coupled by the channel. A measure of that is the median energy
$E_{m}$ in the integrated $T(E)$,

$\int_{-\infty}^{E_{m}}dE\,T(E)=\frac{1}{2}\int_{-\infty}^{\infty}dE\,T(E),$
(23)

measured with respect to the GCM energy at the barrier, $\langle
x_{b}|H|x_{b}\rangle$ in the single-channel examples and the lower eigenvalue
of the Hamiltonian for the case of two coupled channels.

Model | $I_{T}$ | $E_{m}/E_{K}$
---|---|---
A4 | 1.7 | 0.0
A7 | 1.7 | 0.0
B4 | 0.6 | -0.1
B7 | 0.6 | -0.3
C8 | 3.1 | 0.9
C${}^{a}_{8}$ | 2.5 | 0.9
Table 1: Integrated channel properties of the models discussed in the text.
The model labels refer to the subsection in Sect. IV where they were
discussed. The subscript refers to the dimension of the GCM space. $E_{m}$ is
computed with respect to the diagonal GCM energy at the barrier peak. In case
C, the energy is the adiabatic one computed by diagonalizing the 2x2 matrix
mixing the two GCM states. The row marked C${}^{a}_{8}$ ignores the coupling
between the adiabatic channels.

The energy difference scaled to $E_{K}$ is given in the third column of Table
I.

Both quantities exhibited in the Table support our conclusion that one can
safely use coarse meshes to define the channels. Comparing the 4-state models
with the 7-state models one sees little or no change in the integrated
transmission factor and in its median energy. Model C with two interacting
diabatic channels has about twice the integrated transmission strength as
model A, which is hardly surprising.

Another point to note about the phenomenological transition-state theory is
that the coupling between adiabatic channels is completely neglected. We can
test this assumption by the entry marked “C${}^{a}_{8}$”. That model ignores
the coupling between the adiabatic channels; its transmission factor is lower
by 20%. This is not beyond the expected uncertainties in a fully microscopic
theory, but a computational framework taking all interactions into account is
obviously to be preferred.

## References

  * (1) N. Bohr and J.A. Wheeler, Phys. Rev. 56, 426 (1939).
  * (2) O. Bouland, J.E. Lynn, and P. Talou, Phys. Rev. C 88 054612 (2013).
  * (3) R. Capote et al., Nucl. Data Sheets 110, 3107 (2009).
  * (4) M.B. Chadwick et al., Nucl. Data Sheets 107, 2931 (2006).
  * (5) M.B. Chadwick et al., Nucl. Data Sheets 112, 2887 (2011).
  * (6) T. Cap, K. Siwek-Wilczynska, and J. Wilczynski, Phys. Rev. C83, 054602 (2011).
  * (7) H. Lu, A. Marchix, Y. Abe and D. Boilley, Comp. Phys. Comm. 200, 381 (2016).
  * (8) K.-H. Schmidt and W. Morawek, Rep. Prog. Phys. 54, 949 (1991).
  * (9) M. Bender, P.H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75, 121 (2003).
  * (10) H. Goutte, et al., Phys. Rev. C 71 024316 (2005).
  * (11) H. Tao, et al., Phys. Rev. C 96 024319 (2017).
  * (12) G.F. Bertsch and K. Hagino, J. Phys. Soc. Jpn. 90 114005 (2012).
  * (13) K. Hagino and G.F. Bertsch, Phys. Rev. E104, L052104 (2021).
  * (14) S. Datta, Electronic Transport in Mesoscopic Systems, (Cambridge University Press, Cambridge, 1995), Eq. (3.5.20).
  * (15) Y. Alhassid, G.F. Bertsch, and P. Fanto, Ann. Phys. 424 168381 (2021).
  * (16) G.F. Bertsch and L.M. Robledo, Phys. Rev. C100, 044606 (2019).
  * (17) G.F. Bertsch and W.Younes, Ann. Phys. 403 68 (2019).
  * (18) K. Hagino and G.F. Bertsch, Phys. Rev. C 101 064317 (2020).
  * (19) S.A. Guilini, L.M. Robledo, and R. Rodriguez-Guzmán, Phys. Rev. C. 90 054311 (2014).