This manuscript has been authored by UT-Battelle, LLC, under Contract No.
DE6AC05-00OR22725 with the U.S. Department of Energy. The U.S. Government is
authorized to reproduce and distribute reprints for Government purposes
notwithstanding any copyright notation hereon. The Department of Energy will
provide public access to these results of federally sponsored research in
accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-
publicaccess-plan).

# Modeling Kinetic Effects of Charged Vacancies on Electromechanical Responses
of Ferroelectrics: Rayleighian Approach

Rajeev Kumar<EMAIL_ADDRESS>Center for Nanophase Materials Sciences, Oak
Ridge National Laboratory, Oak Ridge, TN-37831  Shuaifang Zhang  Center for
Nanophase Materials Sciences, Oak Ridge National Laboratory, Oak Ridge,
TN-37831  P. Ganesh<EMAIL_ADDRESS>Center for Nanophase Materials Sciences,
Oak Ridge National Laboratory, Oak Ridge, TN-37831

###### Abstract

Understanding time-dependent effects of charged vacancies on electromechanical
responses of materials is at the forefront of research for designing materials
exhibiting metal-insulator transition, and memresistive behavior. A
Rayleighian approach is used to develop a model for studying the non-linear
kinetics of reaction leading to generation of vacancies and electrons by the
dissociation of vacancy-electron pairs. Also, diffusion and elastic effects of
charged vacancies are considered to model polarization-electric potential and
strain-electric potential hysteresis loops. The model captures multi-physics
phenomena by introducing couplings among polarization, electric potential,
stress, strain, and concentrations of charged (multivalent) vacancies and
electrons (treated as classical negatively charged particles), where the
concentrations can vary due to association-dissociation reactions. Derivation
of coupled time-dependent equations based on the Rayleighian approach is
presented. Three limiting cases of the governing equations are considered
highlighting effects of 1) non-linear reaction kinetics on the generation of
charged vacancies and electrons, 2) the Vegard’s law (i.e., the concentration-
dependent local strain) on asymmetric strain-electric potential relations, and
3) coupling between a fast component and the slow component of the net
polarization on the polarization-electric field relations. The Rayleighian
approach discussed in this work should pave the way for developing a multi-
scale modeling framework in a thermodynamically consistent manner while
capturing multi-physics phenomena in ferroelectric materials.

##  I Introduction

Fundamental principles, which can help in design of materials exhibiting
metal-insulator transitionmit_review  and memory effects, are under extensive
scrutiny by a number of researchers in the fields of
spintronicsawschalom2007challenges , ferroelectricsrabe2007modern , and
neuromorphic computingfe_neuro . Experimental discoveriesparkin_science
connecting oxygen vacancies to unconventional ferroelectricity in films and
metal-insulator transitions have led to research activities focused on
understanding effects of the vacancies on material properties. These
activities include experimental measurements, which can decouple contributions
to switching of polarization from polar crystal phases and vacancy migration
under applied electric fields, transient effects of vacancy
orderingkelley2022oxygen , and electrochemical effects due to
surfacesyang2017mixed . In parallel, models capturing effects of oxygen
vacancies at different levels of approximations have been developed, which
include molecular models based on density functional theorykelley2022oxygen ,
reactive force fieldsakbarian2019understanding , and continuum
modelsxiao2008continuum  capturing long-wavelength physics. These research
activities, experimental and theoretical, have unequivocally established that
electrostatic effects due to vacancies are of paramount importance in
controlling properties of materials at length and time scales relevant to
application of materials in microelectronics. For example, Lin et
al.lin2020oxygen  proposed a hypothesis that enhanced remnant polarization of
a composite containing ferroelectric barium titante (BaTiO3) and metallic non-
ferroelectric oxide of SrRuO3 may be attributed to the accumulation of oxygen
vacancies at the BaTiO3/SrRuO3 interface. Similarly, thin films ($\sim 80$ nm)
of $\mbox{BaTiO}_{3}$ grown on $\mbox{SrRuO}_{3}//\mbox{SrTiO}_{3}$ have been
studied recentlykelley2022oxygen  to understand effects of oxygen vacancy
injection on electromechanical responses of the films. It has been
establishedkelley2022oxygen  that enhanced electromechanical responses in
these films can be sustained by injection of oxygen vacancies and kinetic
effects due to the vacancies were postulated to be responsible for the
responses.

Modeling kinetic effects due to vacancies in thin films of ferroelectric
materials requires use of a simulation method, which can capture effects of
multiple phases (varying in polarization and crystal symmetry) and multi-
physics phenomena due to couplings among polarization, electrostatic
potential, strain, and vacancies. For long length and time scales, phase field
method has been used to model ferroelectricli2001phase ; saha2019phase ;
wang2019understanding  materials, using parameters, which can be either
inferred from experimental data or estimated using density functional theory
(DFT). The phase field method has been used to model phenomena occurring at
long time scales (of the order of seconds), such as domain wall motion
li2001phase  and domain nucleationsaha2019phase . The method incorporates a
thermodynamic free energy similar to the originally proposed by Landau,
Ginzburg, and Devonshire for ferroelectricity in the bulk (i.e., without
surfaces), generalized to inhomogeneous systems and time domain. The free
energy is typically written in terms of a space and time dependent
polarization vector, and its gradients. Coefficients of various terms in the
free energy can be either kept phenomenological or estimated from first
principles calculations such as DFT after identifying underlying origins of
polarization in terms of atomic displacements . For example, the free energies
for some ferroelectrics materials are summarized by Chen chen2007appendix .

In the last two decades, a number of phase field models have been developed
for understanding the effects of vacancies in ferroelectric materials. For
example, Zhang et al. zhang2010oxygen  developed a phase field model for
ferroelectrics with oxygen vacancies by treating them as defect dipoles and
investigated the oxygen-vacancy-induced memory effect and large recoverable
strain in BaTiO3. Cao et al. cao2014effect  treated the vacancies as charged
defects and developed a phase field method for BaTiO3 by coupling time
dependent equations for polarization with the Poisson-Nernst-Planck equations
for density of charges. Shindo et al.shindo2015phase  used a phase field model
to study the electromechanical response of polycrystal BaTiO3 with oxygen
vacancies by treating the vacancies as defect dipoles. Recently, Fedeli et
al.fedeli2019phase  presented a phase field model for ferroelectrics by
treating defects in the single crystal and polycrystal structures as voids,
charged point defect and polarization pinned objects, such that the
polarization pinned objects retain their polarization during the cycling of
electric field. Lastly, Kelley et al. kelley2022oxygen  proposed a new free
energy by including the impact of oxygen vacancies on not only polarization,
but also on strains, which they used to understand asymmetry in the
polarization-electrostatic potential ($P-\psi$) loops with different
concentrations of oxygen vacancies. In order to simply the semi-analytical
analysis, they invoked an assumption that the vacancies move much faster than
the polarization so that a steady state approximation for relaxation of the
densities of the vacancies can be justified. Despite these modeling efforts
spanning two decades related to the modeling of vacancies and ferroelectrics,
a unified theoretical framework, which can be used to construct
thermodynamically consistent kinetic models for thin films of ferroelectric
materials with vacancies is still lacking. Development of such a theoretical
framework is the main goal of this paper. Limiting cases are discussed to
highlight key features of this framework.

In this manuscript, a Rayleighian approach is used to develop time and space
dependent equations for polarization, densities of vacancies, electrons, and
their pairs, electrostatic potential, and strains. The approach has its basis
in principles of linear irreversible thermodynamics and ensures that the
second law of thermodynamics is obeyed. Numerical solutions of coupled set of
equations will be presented in another manuscript for the purpose of
constructing average polarization-electric potential and strain-electric
potential hysteresis loops. The manuscript is organized as follows: details of
energy functional, the Rayleighian, and the coupled equations are presented in
the section II. Different limits of the coupled set of equations are analyzed
in section III. Conclusions are presented in section IV.

##  II Model building for thin films of ferroelectrics using the Rayleighian
approach

A set of time-dependent equations can be derived for simulating ferroelectric
materials like $\mbox{BaTiO}_{3}$ with vacancies. These equations can be
derived using principles of irreversible thermodynamics. For the derivation,
we decompose local polarization into fast and slow components resulting from
electronic and atomic/orientational motions, respectively. The fast component
of the polarization ($=\mathbf{P}_{e}(\mathbf{r},t)$) at a location
$\mathbf{r}$ at time $t$ is assumed to be in equilibrium with the electric
field ($\mathbf{E}(\mathbf{r},t)=-\nabla\psi(\mathbf{r},t)$) i.e.,
$\mathbf{P}_{e}(\mathbf{r},t)=\epsilon_{0}\epsilon_{\infty}\mathbf{E}(\mathbf{r},t)$
so that $\epsilon_{0}$ is the permittivity of vacuum, $\epsilon_{\infty}$ is
the infinite frequency dielectric constant of the material, and $\psi$ is the
electrostatic potential. In contrast, the slow component
($=\mathbf{P}(\mathbf{r},t)$) can be out-of-equilibrium and time-dependent
equations are derived for this component. Such a decomposition of the
polarization is similar to the work by Marcusmarcus_1  focused on electron
transfer processes. Coupling of the slow component of polarization with the
fast component is considered by introducing lattice strain. Diffusion of
charged vacancies and electrons are considered by using their local densities
as additional order parameters, and following our previous work related to ion
transport in polymerized ionic liquidskumar2017 . Furthermore, reaction
kinetics leading to electrochemical effects, which result from charge
generation and recombination of positively charged vacancies with electrons,
are considered. Strain and concentration of the vacancies are coupled by the
Vegard’s lawashcroftPRA , which is an empirical rule leading to a linear
relation between the strain and the concentration of the vacancies.
Mathematical framework leading to equations, which can introduce couplings
among polarization, strain, diffusion, and reaction kinetics, is presented
below.

Any material perturbed from equilibrium by the application of a small external
force will try to approach equilibrium through various dissipative processes.
Dissipative processes and the path must lead to the positive rate of entropy
production as per the second law of thermodynamicsonsager1 ; groot1962non ;
kumar2017 . These statements were cast into a rigorous mathematical framework
by Onsageronsager1  after defining a functional called the Rayleighian, which
consists of the rate of change of “free” energy (or entropy for isothermal
processes) and a dissipation function. According to Onsager’s variational
principle, the true dynamics of a system is one that demonstrates the least
dissipation. Linear relations among various fluxes and forces can hence be
derived by minimizing the Rayleighian with respect to the velocities leading
to the maximum rate of entropy production and the resulting equations define
the most probable path towards an equilibriummachulp1953 ; doi2019 . Details
of this so-called Onsager variational principle are presented
elsewhereonsager1 ; machulp1953 ; Doi_2011 ; doi2019 . Construction of the
Rayleighian requires identification of independent fluxes, reaction rates, and
constraints, which are discussed in the next subsection.

###  II.1 Independent Fluxes, Rates, and Constraints

A model for understanding effects of vacancies can be developed using the
principles of irreversible thermodynamicsonsager1 ; groot1962non ; kumar2017 .
In particular, we use Onsager’s variational principleDoi_2011  to derive a set
of equations for the polarization ($\boldsymbol{P}(\mathbf{r},t)$), strain
tensor ($\boldsymbol{\varepsilon}(\mathbf{r},t)$), number densities of
electroactive vacancies ($\rho_{+}(\mathbf{r},t)$), electrons
($\rho_{-}(\mathbf{r},t)$), and their pairs ($\rho_{\pm}(\mathbf{r},t)$),
which are coupled via the electric field ($\boldsymbol{E}(\mathbf{r},t)$). The
number densities, $\rho_{+}(\mathbf{r},t)$, $\rho_{-}(\mathbf{r},t)$, and
$\rho_{\pm}(\mathbf{r},t)$ are assumed to satisfy

$\displaystyle\frac{\partial\rho_{i}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$
$\displaystyle-\nabla\cdot\mathbf{j}_{i}(\mathbf{r},t)+S_{i}(\mathbf{r},t)\quad\mbox{for}\quad
i=+,-,\pm$ (1)

where
$\mathbf{j}_{i}(\mathbf{r},t)=\rho_{i}(\mathbf{r},t)\mathbf{v}_{i}(\mathbf{r},t)$
is the diffusive flux so that $\mathbf{v}_{i}(\mathbf{r},t)$ is the collective
velocity of $i=+,-,\pm$. In Eq. 1, $S_{i}(\mathbf{r},t)$ is the rate of change
of the number density of $i$ resulting from the dissociation and the formation
of vacancy-electron pairs. We assume that the number densities satisfy the no-
void condition at all locations at all times so that the sum of the volume
fractions defined as $\phi_{i}(\mathbf{r},t)=\rho_{i}(\mathbf{r},t)/\rho_{io}$
is unity i.e.,

$\displaystyle\sum_{i=+,-,\pm}\phi_{i}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle 1$ (2)

where $1/\rho_{io}$ is the molar volume of $i$. Using Eq. 1 and 2, a
constraint on the fluxes and the reaction rates is obtained, which is written
as

$\displaystyle\sum_{i=+,-,\pm}\left[-\nabla\cdot\left\\{\frac{\mathbf{j}_{i}(\mathbf{r},t)}{\rho_{io}}\right\\}+\frac{S_{i}(\mathbf{r},t)}{\rho_{io}}\right]$
$\displaystyle=$ $\displaystyle 0$ (3)

Now, we need to construct relations between $\mathbf{j}_{i}(\mathbf{r},t)$,
$S_{i}(\mathbf{r},t)$ and thermodynamic forces. In the following, this is
accomplished by considering a functional called the Rayleighian ($R$). Before
defining the Rayleighian, we need to identify relevant indepedent fluxes
($\mathbf{j}_{i}$), rates ($S_{i}$), and thermodynamic forces.

For identifying independent fluxes ($\mathbf{j}_{i}$), rates ($S_{i}$), and
other constraints, we consider the rate of change of the total number
($=M(t)$) of particles in the vacancies, the electrons, and the pairs, defined
as $M(t)=\int
d\mathbf{r}\,\left[\rho_{+}(\mathbf{r},t)+\rho_{-}(\mathbf{r},t)+(z_{+}+1)\rho_{\pm}(\mathbf{r},t)\right]$.
Due to the fixed number of particles at all times, the rate of change of
$M(t)$ must be zero and the rate is given by

$\displaystyle\frac{dM(t)}{dt}$ $\displaystyle=$ $\displaystyle\int
d\mathbf{r}\,\left[-\nabla\cdot\left\\{\mathbf{j}_{+}(\mathbf{r},t)+\mathbf{j}_{-}(\mathbf{r},t)+(z_{+}+1)\mathbf{j}_{\pm}(\mathbf{r},t)\right\\}\right.$
(4)
$\displaystyle\left.+\left\\{S_{+}(\mathbf{r},t)+S_{-}(\mathbf{r},t)+(z_{+}+1)S_{\pm}(\mathbf{r},t)\right\\}\right]$
$\displaystyle=$ $\displaystyle\int
d\mathbf{\Gamma}\,\mathbf{\hat{n}}\cdot\left\\{\mathbf{j}_{+}(\mathbf{r},t)+\mathbf{j}_{-}(\mathbf{r},t)+(z_{+}+1)\mathbf{j}_{\pm}(\mathbf{r},t)\right\\}$

where we have used Eqs. 1 and assumed
$S_{+}(\mathbf{r},t)+S_{-}(\mathbf{r},t)+(z_{+}+1)S_{\pm}(\mathbf{r},t)=0$ so
that the total number of particles remains the same due to the dissociation-
association reactions. Here, $\Gamma$ represents the surface enclosing the
volume under consideration so that $\mathbf{\hat{n}}$ is an outward normal at
the surface and we have used the divergence theorem. Now, equating
$dM(t)/dt=0$, Eq. 4 is satisfied if

$\displaystyle\mathbf{\hat{n}}\cdot\left\\{\mathbf{j}_{+}(\mathbf{r},t)+\mathbf{j}_{-}(\mathbf{r},t)+(z_{+}+1)\mathbf{j}_{\pm}(\mathbf{r},t)\right\\}$
$\displaystyle=$ $\displaystyle 0$ (5)

at the surface. Another constraint for the boundary fluxes is obtained by
considering the rate of change of the total charge ($=C(t)$) of the vacancies,
and the electrons, defined as $C(t)=e\int
d\mathbf{r}\left[z_{+}\rho_{+}(\mathbf{r},t)-\rho_{-}(\mathbf{r},t)\right]$,
which must be zero due to the global electroneutrality at all times. The rate
of change of $C(t)$ is given by

$\displaystyle\frac{dC(t)}{dt}$ $\displaystyle=$ $\displaystyle e\int
d\mathbf{r}\,\left[-\nabla\cdot\left\\{z_{+}\mathbf{j}_{+}(\mathbf{r},t)-\mathbf{j}_{-}(\mathbf{r},t)\right\\}+\left\\{z_{+}S_{+}(\mathbf{r},t)-S_{-}(\mathbf{r},t)\right\\}\right]$
(6) $\displaystyle=$ $\displaystyle\int
d\mathbf{\Gamma}\,\mathbf{\hat{n}}\cdot\left\\{z_{+}\mathbf{j}_{+}(\mathbf{r},t)-\mathbf{j}_{-}(\mathbf{r},t)\right\\}=0$

where we have assumed that $z_{+}S_{+}(\mathbf{r},t)-S_{-}(\mathbf{r},t)=0$
i.e., the rate of charge generation due to the dissociation-association
reactions is taken to be zero. Combining
$S_{+}(\mathbf{r},t)+S_{-}(\mathbf{r},t)+(z_{+}+1)S_{\pm}(\mathbf{r},t)=0$
with $z_{+}S_{+}(\mathbf{r},t)-S_{-}(\mathbf{r},t)=0$ leads to

$\displaystyle S_{+}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{z_{+}}S_{-}(\mathbf{r},t)=-S_{\pm}(\mathbf{r},t)\equiv-S(\mathbf{r},t)$
(7)

for the association-dissociation reactions involving electroneutral vacancy-
electron pairs, positively charged vacancies with charge $z_{+}e$ (so that $e$
is the charge of an electron) and neutralizing electrons. For example,
$z_{+}=2$ in the case of oxygen vacancies. Eq. 7 implies that there is only
one independent reaction rate, which is taken to be $S(\mathbf{r},t)$ and
defined by Eq. 7. Now, Eqs. 3 and 7 imply that there are only two independent
fluxes out of the three fluxes, $\mathbf{j}_{i}$. Based on a molecular
description of the diffusion resulting from frictional forces and relative
motion of molecules, we work with the relative fluxes defined as
$\mathbf{\hat{j}}_{i}(\mathbf{r},t)=\phi_{i}(\mathbf{r},t)\left[\mathbf{v}_{i}(\mathbf{r},t)-\mathbf{v}(\mathbf{r},t)\right]$,
where

$\displaystyle\mathbf{v}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\sum_{i=+,-,\pm}\phi_{i}(\mathbf{r},t)\mathbf{v}_{i}(\mathbf{r},t)$
(8)

These relative fluxes are related by the relation

$\displaystyle\sum_{i=+,-,\pm}\mathbf{\hat{j}}_{i}(\mathbf{r},t)$
$\displaystyle=$ $\displaystyle 0$ (9)

and the Rayleighian can be written in terms of any two relative fluxes. In
here, we choose $\mathbf{\hat{j}}_{+}$ and $\mathbf{\hat{j}}_{-}$ to study the
effects of diffusion of the vacancies and the electrons, respectively. Eqs. 1
are rewritten in terms of these independent fluxes and rates in the form

$\displaystyle\frac{\partial\phi_{k}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$
$\displaystyle-\nabla\cdot\mathbf{\hat{j}}_{k}(\mathbf{r},t)+\frac{1}{\rho_{ko}}S_{k}(\mathbf{r},t)-\nabla\cdot\left[\mathbf{\phi}_{k}(\mathbf{r},t)\mathbf{v}(\mathbf{r},t)\right]\quad\mbox{for}\,k=+,-.$
(10)

and Eq. 3 can be written as

$\displaystyle\nabla\cdot\mathbf{v}(\mathbf{r},t)+S(\mathbf{r},t)\left\\{\frac{1}{\rho_{+o}}+\frac{z_{+}}{\rho_{-o}}-\frac{1}{\rho_{\pm
o}}\right\\}$ $\displaystyle=$ $\displaystyle 0$ (11)

Eqs. 4 and 6 are rewritten in terms of
$\mathbf{\hat{j}}_{i=+,-,\pm}\cdot\mathbf{\hat{n}}$ as

$\displaystyle\mathbf{\hat{n}}\cdot\left\\{\rho_{+o}\mathbf{\hat{j}_{+}}(\mathbf{r},t)+\rho_{-o}\mathbf{\hat{j}_{-}}(\mathbf{r},t)+(z_{+}+1)\rho_{\pm
o}\mathbf{\hat{j}_{\pm}}(\mathbf{r},t)\right\\}$
$\displaystyle+\mathbf{\hat{n}}\cdot\left\\{\left[\rho_{+o}\phi_{+}(\mathbf{r},t)+\rho_{-o}\phi_{-}(\mathbf{r},t)+(z_{+}+1)\rho_{\pm
o}\phi_{\pm}(\mathbf{r},t)\right]\mathbf{v}(\mathbf{r},t)\right\\}$
$\displaystyle=$ $\displaystyle 0$ (12)
$\displaystyle\mathbf{\hat{n}}\cdot\left\\{z_{+}\rho_{+o}\mathbf{\hat{j}_{+}}(\mathbf{r},t)-\rho_{-o}\mathbf{\hat{j}_{-}}(\mathbf{r},t)+\left[z_{+}\rho_{+o}\phi_{+}(\mathbf{r},t)-\rho_{-o}\phi_{-}(\mathbf{r},t)\right]\mathbf{v}(\mathbf{r},t)\right\\}$
$\displaystyle=$ $\displaystyle 0$ (13)

These equations show that there are two boundary fluxes, which are
independent. These are chosen to be
$\mathbf{\hat{j}}_{k=+,-}\cdot\mathbf{\hat{n}}$. In summary, independent
fluxes and rates are
$\mathbf{\hat{j}}_{k=+,-},\mathbf{\hat{j}}_{k=+,-}\cdot\mathbf{\hat{n}}$, and
$S$, respectively. In addition, we need to consider the constraint written as
Eq. 11 and discussed above.

###  II.2 Rayleighian

For applying the Onsager’s variational principle, a Rayleighian ($R(t)$) for
the thin films of ferroelectric materials can be defined askumar2017 ;
cummings2018

$\displaystyle R(t)$ $\displaystyle=$ $\displaystyle\frac{dH(t)}{dt}+W(t)-\int
d\mathbf{r}\,p(\mathbf{r},t)\left[\mathbf{\nabla}\cdot\mathbf{v}(\mathbf{r},t)+S(\mathbf{r},t)\left\\{\frac{1}{\rho_{+o}}+\frac{z_{+}}{\rho_{-o}}-\frac{1}{\rho_{\pm
o}}\right\\}\right]$ (14)

where $H(t)$ is the time-dependent energy functional (presented in the next
subsection) and $p(\mathbf{r},t)$ is a Lagrange’s multiplier to enforce the
constraint (cf. Eq. 11), and $W$ is the dissipation function. We should point
out that the functional form for the dissipation function is assumed to be
known in order to use the variational principle and there is no prescription
for deriving it. In this paper, we present a functional form for the
dissipation function by including various multi-physics phenomena and ensuring
that the derived equations lead to known relations in various limiting cases.
For example, $W$ has contributions from the coupling of the polarization with
the lattice velocity, derived by Hubbard and Onsagerhubbard1977dielectric
using the approximation of fast rotational relaxationgroot1962non ;
felderhof1999hydrodynamics . Also, contributions from friction, which result
from relative motions of charged vacancies and electrons with respect to the
velocity $\mathbf{v}(\mathbf{r},t)$ (cf. Eq. 8), are includedkumar2017 ;
Doi_2011 . Adding these contributions, $W$ can be written as

$\displaystyle W(t)$ $\displaystyle=$ $\displaystyle\frac{\tau_{p}}{2}\int
d\mathbf{r}\,\left[\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}+\left[\mathbf{v}_{l}(\mathbf{r},t)\cdot\nabla\right]\mathbf{P}(\mathbf{r},t)+\frac{1}{2}\mathbf{P}(\mathbf{r},t)\times\left[\nabla\times\mathbf{v}_{l}(\mathbf{r},t)\right]\right]^{2}$
(15) $\displaystyle+\frac{1}{2}\int
d\mathbf{r}\sum_{k=+,-}\sum_{k^{\prime}=+,-}L_{kk^{\prime}}(\mathbf{r},t)\rho_{ko}\rho_{k^{\prime}o}\mathbf{\hat{j}}_{k}(\mathbf{r},t)\cdot\mathbf{\hat{j}}_{k^{\prime}}(\mathbf{r},t)$
$\displaystyle+\frac{1}{2}\int
d\mathbf{r}\sum_{k=+,-}\sum_{k^{\prime}=+,-}M_{kk^{\prime}}(\mathbf{r},t)\rho_{ko}\rho_{k^{\prime}o}\left[\mathbf{\hat{j}}_{k}(\mathbf{r},t)\cdot\mathbf{\hat{n}}\right]\left[\mathbf{\hat{j}}_{k^{\prime}}(\mathbf{r},t)\cdot\mathbf{\hat{n}}\right]$
$\displaystyle+\frac{1}{2}\int
d\mathbf{r}\,\omega(\mathbf{r},t)S^{2}(\mathbf{r},t)$

Here, $\tau_{p}$ and $L_{kk^{\prime}}$ are parameters characterizing time-
scale for change in polarization $\mathbf{P}$, and friction coefficient for
motion of the vacancies, and electrons relative to $\mathbf{v}(\mathbf{r},t)$
(defined by Eq. 8). Similarly, $M_{kk^{\prime}}$ are the parameters
characterizing the dissipation due to the relative fluxes at the boundaries.
Also, $\mathbf{v}_{l}=\sum_{l=1,2,3}\left[\frac{\partial
u_{l}(\mathbf{r},t)}{\partial t}\right]\hat{i}_{l}$ is net velocity of an
underlying lattice so that $\hat{i}_{l}$ are unit vectors and
$u_{l}(\mathbf{r},t)$ is the displacement of underlying atoms at location
$\mathbf{r}$ at time t. In general, $L_{kk^{\prime}}$ and $M_{kk^{\prime}}$
can be concentration dependent but each matrix with either $L_{kk^{\prime}}$
or $M_{kk^{\prime}}$ as its elements must be positive definite for the
positive entropy production. The last term in Eq. 15 is the dissipation due to
the reaction with a prefactor $\omega$, which will be related to the rate of
vacancy-electron pair dissociation and recombination.

In the following, we present explicit expression for $H(t)$ by using
thermodynamic free energy and by generalizing it after considering additional
effects of the charged vacancies.

###  II.3 Energy Functional

For deriving a set of time-dependent equations, we use a time ($t$)-dependent
energy functional, written as

$\displaystyle H(t)$ $\displaystyle=$ $\displaystyle\int
d\mathbf{r}\,\left[H_{LGD}\left\\{\boldsymbol{P}\right\\}+H_{grad}\left\\{\boldsymbol{\nabla
P}\right\\}+H_{mech}\left\\{\boldsymbol{P},\boldsymbol{\varepsilon},\rho_{+}\right\\}+H_{self}\left\\{\rho_{+},\rho_{-},\rho_{\pm}\right\\}\right.$
(16)
$\displaystyle\left.+H_{elec}\left\\{\boldsymbol{P},\boldsymbol{E}=-\nabla\psi,\rho_{+},\rho_{-}\right\\}+H_{mix}\left\\{\rho_{+},\rho_{-},\rho_{\pm},\nabla\rho_{+},\nabla\rho_{-},\nabla\rho_{\pm}\right\\}\right],$
$\displaystyle\equiv$ $\displaystyle\int
d\mathbf{r}\,h\left\\{\boldsymbol{P},\boldsymbol{\nabla
P},\rho_{+},\rho_{-},\rho_{\pm},\nabla\rho_{+},\nabla\rho_{-},\nabla\rho_{\pm},\boldsymbol{E},\boldsymbol{\varepsilon}\right\\}$
(17)

where $H_{LGD}$ is the Landau-Ginzburg-Devonshire (LGD) energy
densityrabe2007modern  and written in terms of time-dependent polarization by
invoking local-equilibrium approximationgroot1962non . Similarly, $H_{grad}$
is the gradient/interfacial energy density capturing the effects of
inhomogeneous polarization in the long-wavelength limit. Explicit expressions
for these contributions are presented in Appendix A for $\mbox{BaTiO}_{3}$ in
a seminal work by Chen and co-workerschen2007appendix . Coupling between the
polarization and the strain is encoded in $H_{mech}$, which is the mechanical
strain energy density so that $\boldsymbol{\varepsilon}$ is the strain tensor.
Considering limit of small deformation, the mechanical strain energy
densityLandau1986  can be defined as

$H_{mech}=\frac{1}{2}\boldsymbol{\sigma}(\mathbf{r},t):\left[\boldsymbol{\varepsilon}(\mathbf{r},t)-\boldsymbol{\varepsilon^{0}}(\mathbf{r},t)\right]\equiv\frac{C_{ijkl}}{2}\left[\varepsilon_{ij}(\mathbf{r},t)-\varepsilon_{ij}^{0}(\mathbf{r},t)\right]\left[\varepsilon_{kl}(\mathbf{r},t)-\varepsilon_{kl}^{0}(\mathbf{r},t)\right],$
(18)

where $\varepsilon_{ij}(\mathbf{r},t)=\left[\partial u_{i}/\partial
x_{j}+\partial u_{j}/\partial x_{i}\right]/2$ is the $ij$ element of the total
lattice-strain tensormura1987general  so that $x_{i}$ are components of the
spatial vector $\mathbf{r}$, $u_{i}(\mathbf{r},t)$ is $i^{th}$ component of
the displacement vector of lattice, $\boldsymbol{\sigma}$ is the stress
tensor, and $C_{ijkl}$ is the rank four elasticity tensor. Ferroelectric
materials can have spontaneous strain even in stress-free conditions and such
a strain tensor (eigenstrainmura1987general ) is denoted as
$\varepsilon_{ij}^{0}$, results from electrostriction and Vegard effects, and
can be defined aschen2002phase ; chen2007appendix ; kelley2022oxygen ,

$\varepsilon_{ij}^{0}(\mathbf{r},t)=Q_{ijkl}P_{k}(\mathbf{r},t)P_{l}(\mathbf{r},t)+w_{ij}^{v}\rho_{+}(\mathbf{r},t),$
(19)

where $Q_{ijkl}$ is a rank four order tensor. In this notation,
$\varepsilon_{ij}(\mathbf{r},t)$ is the tensor containing both, elastic and
eigenstrain, components of the strain. Furthermore, effects of the vacancies
on the strain is included by the last term in Eq. 19, where $w_{ij}^{v}$ are
phenomenological parameters. Here, Einstein’s notation of sum over repeated
indices is used.

$H_{self}$ is the self energy density for creating the vacancies, electrons
and their pairs, written as

$H_{self}=\sum_{i=+,-,\pm}G_{io}\rho_{i}(\mathbf{r},t)$ (20)

where $G_{io}$ is the self-energyisraelachvili92a  for creating $i$.
$H_{elec}$ is the excess electrical energy density written asmarcus_1

$H_{elec}=\left[z_{+}\rho_{+}(\mathbf{r},t)-\rho_{-}(\mathbf{r},t)\right]e\psi(\mathbf{r},t)-\frac{\epsilon_{0}\epsilon_{\infty}}{2}\mathbf{E}^{2}(\mathbf{r},t)-\left[\mathbf{P}(\mathbf{r},t)\cdot\mathbf{E}(\mathbf{r},t)\right],$
(21)

where $\mathbf{E}(\mathbf{r},t)=-\nabla\psi(\mathbf{r},t)$ is the electric
field and $\psi$ is the electrostatic potential. $z_{+}$ is valency of oxygen
vacancy and $e$ is charge of an electron. Here, the pairs of the vacancies and
the electrons are assumed to carry no charge. Furthermore,
$-\frac{\delta\left[\int d\mathbf{r}H_{elec}\right]}{\delta
E(\mathbf{r}^{\prime},t)}=\mathbf{D}(\mathbf{r}^{\prime},t)=\epsilon_{0}\epsilon_{\infty}\mathbf{E}(\mathbf{r}^{\prime},t)+\mathbf{P}(\mathbf{r}^{\prime},t)$
can be readily identified as the dielectric displacement vector.

Excess entropy of mixing vacancies, electrons, and their pairs is defined as
$H_{mix}$ along with the entropic cost of generating their inhomogeneous
density profiles, written askumar2017

$\displaystyle H_{mix}$ $\displaystyle=$ $\displaystyle
k_{B}T\sum_{i=+,-,\pm}\left[\rho_{i}(\mathbf{r},t)\ln\left[\frac{\rho_{i}(\mathbf{r},t)}{\rho_{io}}\right]+\frac{1}{2}\kappa_{i}|\nabla\rho_{i}(\mathbf{r},t)|^{2}\right]$
(22)

Logarithmic terms in Eq. 22 can be derived by considering the number of ways
in which the vacancies, the electrons, and the pairs can be distributed in
space, such that their total number remain fixed. $\kappa_{i}$ is the
coefficient of the square-gradient termkumar2014 , which penalizes
inhomogeneous density profiles of $i$.

With the dissipation and energy functional given by Eqs. 15 and 17,
respectively, the Rayleighian has been specified (cf. Eq. 14) completely with
the quantities like $\tau_{p},L_{kk^{\prime}},M_{kk^{\prime}},$ and $\omega$
taken as inputs. After constructing the Rayleighian, a set of equations can be
systematically derived by optimizing $R$ with respect to
$\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t},\mathbf{\hat{j}}_{i=+,-},\mathbf{\hat{j}}_{i=+,-}\cdot\mathbf{\hat{n}}$,
$S$, and $\mathbf{v}_{l}(\mathbf{r},t)$. The set is complemented by two
additional equations: one for the Lagrange’s multiplier, $p$, and other one
for the electrostatic potential, $\psi$. In total, nine coupled equations are
derived in the next subsection.

###  II.4 Governing Equations: Linear Irreversible Thermodynamics

We assume that the electrostatic potential adjust itself so fast that
stationary condition $\frac{\delta H}{\delta\psi(\mathbf{r},t)}=0$ is
satisfied at all times and at all locations. Explicitly, this leads to

$\epsilon_{0}\epsilon_{\infty}\nabla^{2}\psi(\mathbf{r},t)-\nabla\cdot\mathbf{P}(\mathbf{r},t)+e\left[z_{+}\rho_{+}(\mathbf{r},t)-\rho_{-}(\mathbf{r},t)\right]=0,$
(23)

or equivalently,

$\nabla\cdot\mathbf{D}(\mathbf{r},t)=e\left[z_{+}\rho_{+}(\mathbf{r},t)-\rho_{-}(\mathbf{r},t)\right],$
(24)

where the right hand side is the local charge density. Evaluating $\delta
R(t)/\delta\left\\{\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}\right\\}=0$, we get governing equations for the three components of
polarization $P_{1},P_{2},P_{3}$ such that

$\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}+\left[\mathbf{v}_{l}(\mathbf{r},t)\cdot\nabla\right]\mathbf{P}(\mathbf{r},t)+\frac{1}{2}\mathbf{P}(\mathbf{r},t)\times\left[\nabla\times\mathbf{v}_{l}(\mathbf{r},t)\right]=-\frac{1}{\tau_{p}}\frac{\delta\left[\int
d\mathbf{r}^{\prime}h\left\\{\mathbf{P}(\mathbf{r}^{\prime},t)\right\\}\right]}{\delta\mathbf{P}(\mathbf{r},t)}\equiv\mathbf{P}^{\star}(\mathbf{r},t),$
(25)

Functional $h$ is defined in Eq. 17 and its dependencies on variables other
than the polarization are suppressed here. Similarly, $\delta
R(t)/\delta\mathbf{v}_{l}(\mathbf{r},t)=0$ leads to

$\displaystyle\nabla\cdot\boldsymbol{\sigma}(\mathbf{r},t)=\frac{\tau_{p}}{2}\left\\{\nabla\left[\mathbf{P}^{\star}(\mathbf{r},t)\cdot\mathbf{P}(\mathbf{r},t)\right]-\mathbf{P}(\mathbf{r},t)\times\left[\nabla\times\mathbf{P}^{\star}(\mathbf{r},t)\right]-\mathbf{P}^{\star}(\mathbf{r},t)\times\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]\right.$
$\displaystyle\left.-\mathbf{P}(\mathbf{r},t)\left[\nabla\cdot\mathbf{P}^{\star}(\mathbf{r},t)\right]+\mathbf{P}^{\star}(\mathbf{r},t)\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]\right\\}$
(26)

where we have used
$dH_{mech}/dt=\boldsymbol{\sigma}(\mathbf{r},t):\left[\nabla\mathbf{v}_{l}(\mathbf{r},t)-\frac{\partial\varepsilon_{ij}^{0}(\mathbf{r},t)}{\partial
t}\right]$. With the constitutive relation,

$\sigma_{ij}(\mathbf{r},t)=C_{ijkl}\left[\varepsilon_{kl}(\mathbf{r},t)-\varepsilon_{kl}^{0}(\mathbf{r},t)\right]$
(27)

Eq. 26 can be used to study effects of strains. We should point out that in
general, $\sigma_{ij}(\mathbf{r},t)$ should be computed from gradients of
lattice velocity. Here, we choose a simpler and intuitive linear constitutive
relation (cf. Eq. 27) between stress and strain. Optimizations of $R(t)$ with
respect to $\mathbf{\hat{j}}_{k=+,-}$ and
$\mathbf{\hat{j}}_{k=+,-}\cdot\mathbf{\hat{n}}$ give

$\displaystyle\nabla\left[\mu_{k}(\mathbf{r},t)-\mu_{\pm}(\mathbf{r},t)\right]$
$\displaystyle=$
$\displaystyle-\sum_{k^{\prime}=+,-}L_{kk^{\prime}}(\mathbf{r},t)\rho_{ko}\rho_{k^{\prime}o}\mathbf{\hat{j}}_{k^{\prime}}(\mathbf{r},t)\quad\mbox{in
the volume}$ (28) $\displaystyle\mu_{k}(\mathbf{r},t)-\mu_{\pm}(\mathbf{r},t)$
$\displaystyle=$
$\displaystyle\sum_{k=+,-}M_{kk^{\prime}}(\mathbf{r},t)\rho_{ko}\rho_{k^{\prime}o}\left[\mathbf{\hat{j}}_{k^{\prime}}(\mathbf{r},t)\cdot\mathbf{\hat{n}}\right]\quad\mbox{at
the surface,}$ (29)

respectively. Here, $\mu_{k}(\mathbf{r},t)=\delta\left[\int
d\mathbf{r}^{\prime}h\left\\{\rho_{k}(\mathbf{r}^{\prime},t)\right\\}\right]/\delta\phi_{k}(\mathbf{r},t)$
and becomes the local chemical potential of $k$ in the steady state.

Evaluating $\delta R(t)/\delta S(\mathbf{r},t)=0$ gives

$\displaystyle S(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{1}{\omega(\mathbf{r},t)}\left[\frac{\mu_{+}(\mathbf{r},t)}{\rho_{+o}}+z_{+}\frac{\mu_{-}(\mathbf{r},t)}{\rho_{-o}}-\frac{\mu_{\pm}(\mathbf{r},t)}{\rho_{\pm
o}}+p(\mathbf{r},t)\left\\{\frac{1}{\rho_{+o}}+\frac{z_{+}}{\rho_{-o}}-\frac{1}{\rho_{\pm
o}}\right\\}\right]$ (30)

Plugging the expression for $S$ from Eq. 30 in Eq. 11, an equation for $p$ is
obtained, which closes this set of equations along with the boundary
conditions shown in Eqs. 12 and 13.

###  II.5 Non-linear Reaction Kinetics

Eq. 30 shows that the reaction rate depends linearly on the chemical
potentials, which limits the validity of the model presented here. We can
improve this by introducing non-linear relations between the reaction rate and
the chemical potentials e.g., similar to those forming the basis of the
Eyring’s rate of reactionsGLA41 . In this subsection, we present such a model,
which can lead to the non-linear reaction kinetics such as in autocatalytic
reactionskumar2021harnessing . For such a purpose, we rewrite Eq. 30 as a
limiting case of the non-linear relation

$\displaystyle S_{NL}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle\frac{k_{B}T}{\omega(\mathbf{r},t)}\left[\exp\left[\frac{\mu_{+}(\mathbf{r},t)+p(\mathbf{r},t)}{\rho_{+o}k_{B}T}+z_{+}\frac{\mu_{-}(\mathbf{r},t)+p(\mathbf{r},t)}{\rho_{-o}k_{B}T}\right]\right.$
(31)
$\displaystyle\left.-\exp\left[\frac{\mu_{\pm}(\mathbf{r},t)+p(\mathbf{r},t)}{\rho_{\pm
o}k_{B}T}\right]\right]$

Now, the Rayleighian, $R\equiv R_{NL}$, based on $S_{NL}$ can be constructed
using Eq. 14 and it can be shown that

$\displaystyle R_{NL}(t)$ $\displaystyle=$ $\displaystyle R(t)+\frac{1}{2}\int
d\mathbf{r}\,\omega(\mathbf{r},t)\left[S_{NL}(\mathbf{r},t)-S(\mathbf{r},t)\right]^{2}$
(32)

Note that $R(t)=-W(t)<0$ as per the Onsager’s variational principle, where the
inequality is valid for a positive dissipation function and the governing
equations derived in the last subsection. In contrast, $R_{NL}(t)$ can have
either sign for $\omega(\mathbf{r},t)>0$. We should point out that
$1/\omega(\mathbf{r},t)=\exp\left[\frac{\mu^{\star}(\mathbf{r},t)}{\rho^{\star
o}k_{B}T}\right]\sim\rho^{\star o}>0$ is expected on the basis of the Eyring’s
rate of reactions so that $\mu^{\star}(\mathbf{r},t)$ and $1/\rho^{\star o}$
are the chemical potential and volume of an activated complexGLA41 ,
respectively. Physically, this relation between the prefactor
$\omega(\mathbf{r},t)$ and $\mu^{\star}(\mathbf{r},t)$ implies that the rate
of reaction is linearly proportional to the concentration of the activated
complexesGLA41 . More importantly, probability of realizing a kinetic path
with the non-linear reaction rates is given by the Onsager-Machulp
integralmachulp1953 ; doi2019 , based on time-integral of the difference
$R_{NL}(t)-R(t)$. In particular, the probability of realizing the non-linear
reaction rates should be
$\exp\left[-\int_{0}^{t}dt^{\prime}\,\left[R_{NL}(t^{\prime})-R(t^{\prime})\right]/2k_{B}T\right]$
as per the theoretical works of Onsager and Machulpmachulp1953 . For this
probability to be non-zero and significant enough, $\omega(\mathbf{r},t)$
needs to be chosen in such a manner so that the exponential of the negative of
the Onsager-Machulp integral remain close to unity. In the following, we make
such a choice for $\omega(\mathbf{r},t)$ and keep all of the governing
equations the same except that we replace $S$ by $S_{NL}$ to capture effects
of non-linear reaction kinetics in the model for the ferroelectrics developed
here.

Now, using Eqs. 10, 28 and replacing $S$ with $S_{NL}$ (cf. Eq. 31),

$\displaystyle\frac{\partial\phi_{+}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$
$\displaystyle\nabla\cdot\left[\sum_{k^{\prime}=+,-}\tilde{L}_{+k^{\prime}}^{-1}(\mathbf{r},t)\nabla\left[\mu_{k^{\prime}}(\mathbf{r},t)-\mu_{\pm}(\mathbf{r},t)\right]\right]-\frac{1}{\rho_{+o}}S_{NL}(\mathbf{r},t)-\nabla\cdot\left[\mathbf{\phi}_{+}(\mathbf{r},t)\mathbf{v}(\mathbf{r},t)\right]$
$\displaystyle\frac{\partial\phi_{-}(\mathbf{r},t)}{\partial t}$
$\displaystyle=$
$\displaystyle\nabla\cdot\left[\sum_{k^{\prime}=+,-}\tilde{L}_{-k^{\prime}}^{-1}(\mathbf{r},t)\nabla\left[\mu_{k^{\prime}}(\mathbf{r},t)-\mu_{\pm}(\mathbf{r},t)\right]\right]-\frac{z_{+}}{\rho_{-o}}S_{NL}(\mathbf{r},t)-\nabla\cdot\left[\mathbf{\phi}_{-}(\mathbf{r},t)\mathbf{v}(\mathbf{r},t)\right]$

where $\tilde{L}_{kk^{\prime}}^{-1}$ is the $kk^{\prime}$ element of the
inverse matrix of $L_{kk^{\prime}}\rho_{ko}\rho_{k^{\prime}o}$. Three
independent elements $\tilde{L}_{kk^{\prime}}^{-1}$ can be
interpretednewman2021electrochemical  in terms of the ionic conductivity,
transference number of the vacancies related to their partial ionic currents,
and diffusion constants of the vacancy-electron pairs.

##  III Results

In here, we present analysis of some limiting cases to highlight key effects
of the vacancies and novel aspects of the model. Numerical results obtained by
solving the coupled equations will be presented in a separate publication. In
the model presented here, we can capture non-linear effects of barriers on
reaction rates and study characteristic time for reactions. In the following,
we consider three limiting cases: 1) reaction dominated regime leading to
identification of a characteristic time, 2) steady state analysis highlighting
coupling between the strain, electric potential, and vacancies, 3) vacancy-
free regime, where a coupling between the fast and the slow component of the
polarization can lead to effects of geometry manifesting in the stabilization
of new topological configurations.

###  III.1 Reaction Dominated Regime: Characteristic time

We consider a limiting case, where vacancies, electrons, and their pairs are
homogeneously distributed so that diffusive flux is minimal. In this limit, we
need to consider $\phi_{i=+,-,\pm}(\mathbf{r},t)\equiv\phi_{i}^{h}(t)$, which
satisfy (cf. Eqs. LABEL:eq:phiplusfinal\- LABEL:eq:phinegfinal)

$\displaystyle\frac{\partial\phi_{+}^{h}(t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{1}{\rho_{+o}}S_{NL}^{h}(t)$ (35)
$\displaystyle\frac{\partial\phi_{-}^{h}(t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{z_{+}}{\rho_{-o}}S_{NL}^{h}(t)$ (36)

Now, consider a case, when $\frac{1}{\rho_{\pm
o}}=\frac{1}{\rho_{+o}}+\frac{z_{+}}{\rho_{-o}}$ so that
$\phi_{+}^{h}(t)+\phi_{-}^{h}(t)+\phi_{\pm}^{h}(t)=1$ is satisfied for non-
zero reaction rate, $S_{NL}(\mathbf{r},t)\equiv S_{NL}^{h}(t)$, given by (cf.
Eq. 31)

$\displaystyle S_{NL}^{h}(t)$ $\displaystyle=$ $\displaystyle
K_{0}(t)\left[\exp\left[\frac{\mu_{+}^{h}(t)}{\rho_{+o}k_{B}T}+z_{+}\frac{\mu_{-}^{h}(t)}{\rho_{-o}k_{B}T}\right]-\exp\left[\frac{\mu_{\pm}^{h}(t)}{\rho_{\pm
o}k_{B}T}\right]\right]$ (37)

where we have defined
$K_{0}(t)=\frac{k_{B}T}{\omega^{h}(t)}\exp\left[\frac{p^{h}(t)}{\rho_{\pm
o}k_{B}T}\right]$ and used the notation
$\omega(\mathbf{r},t)\equiv\omega^{h}(t),p(\mathbf{r},t)\equiv
p^{h}(t),\mu_{+}(\mathbf{r},t)\equiv\mu_{+}^{h}(t),\mu_{-}(\mathbf{r},t)\equiv\mu_{-}^{h}(t),\mu_{\pm}(\mathbf{r},t)\equiv\mu_{\pm}^{h}(t)$.
From Eq. 17, we get

$\displaystyle\frac{\mu_{+}^{h}(t)}{\rho_{+0}k_{B}T}$ $\displaystyle=$
$\displaystyle\ln\phi_{+}^{h}(t)+\frac{G_{+0}+z_{+}e\psi^{h}(t)}{k_{B}T}-\frac{w_{ij}^{v}\sigma_{ij}^{h}(t)}{k_{B}T}$
(38) $\displaystyle\frac{\mu_{-}^{h}(t)}{\rho_{-0}k_{B}T}$ $\displaystyle=$
$\displaystyle\ln\phi_{-}^{h}(t)+\frac{G_{-0}-e\psi^{h}(t)}{k_{B}T}$ (39)
$\displaystyle\frac{\mu_{\pm}^{h}(t)}{\rho_{\pm 0}k_{B}T}$ $\displaystyle=$
$\displaystyle\ln\phi_{\pm}^{h}(t)+\frac{G_{\pm 0}}{k_{B}T}$ (40)

where
$\sigma_{ij}(\mathbf{r},t)\equiv\sigma_{ij}^{h}(t),\psi(\mathbf{r},t)\equiv\psi^{h}(t)$.
Using Eqs. 38\- 40, we can write Eq. 37 as

$\displaystyle S_{NL}^{h}(t)$ $\displaystyle=$ $\displaystyle
K_{A}(t)\phi_{+}^{h}(t)\left[\phi_{-}^{h}(t)\right]^{z_{+}}-K_{D}(t)\phi_{\pm}^{h}(t)$
(41) $\displaystyle K_{A}(t)$ $\displaystyle=$ $\displaystyle
K_{0}(t)\exp\left[\frac{G_{+0}+z_{+}G_{-0}}{k_{B}T}-\frac{w_{ij}^{v}\sigma_{ij}^{h}(t)}{k_{B}T}\right]$
(42) $\displaystyle K_{D}(t)$ $\displaystyle=$ $\displaystyle
K_{0}(t)\exp\left[\frac{G_{\pm 0}}{k_{B}T}\right]$ (43)

where $K_{A}(t)$ and $K_{D}(t)$ are defined as the time-dependent association
and dissociation constants, respectively. It should be noted that the reaction
rate $S_{NL}^{h}(t)$ is independent of the electrostatic potential due to the
assumption of a uniform potential. In general, small variations of the
electrostatic potential will lead to the dependence of the reaction rate on
the potential. Furthermore, the reaction rate depends on the stress due to the
Vegard’s law.

Now, we consider the case of oxygen vacancies so that $z_{+}=2$, and assume
that $K_{A}(t)\equiv K_{A0}$ and $K_{D}(t)\equiv K_{D0}$. In this case, we can
solve Eqs. 35, 36, and 41 in terms of a time-dependent parameter, $\alpha(t)$
so that
$\phi_{+}^{h}(t)=\alpha(t)/\rho_{+0},\phi_{-}^{h}(t)=z_{+}\alpha(t)/\rho_{-0}$,
and $\phi_{\pm}^{h}(t)=1-(1/\rho_{+0}+z_{+}/\rho_{-0})\alpha(t)$ and
$\alpha(t)$ satisfies

$\displaystyle\frac{\partial\alpha(t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{4K_{A0}}{\rho_{+o}\rho_{-o}^{2}}\left[\alpha^{3}(t)+\left(\frac{1}{\rho_{+0}}+\frac{2}{\rho_{-0}}\right)\frac{K_{D0}}{K_{A0}}\frac{\rho_{+o}\rho_{-o}^{2}}{4}\alpha(t)-\frac{K_{D0}}{K_{A0}}\frac{\rho_{+o}\rho_{-o}^{2}}{4}\right]$
(44)

Although Eq. 44 can be solved exactly but the exact solution obscure
identification of a characteristic time. In here, we consider a situation so
that Eq. 44 can be approximated as

$\displaystyle\frac{\partial\alpha(t)}{\partial t}$ $\displaystyle=$
$\displaystyle-\frac{4K_{A0}}{\rho_{+o}\rho_{-o}^{2}}\left[\alpha^{3}(t)-\frac{K_{D0}}{K_{A0}}\frac{\rho_{+o}\rho_{-o}^{2}}{4}\right]$
(45)

Integrating Eq. 45 leads to

$\displaystyle\ln\left[\frac{\frac{\alpha(t)}{\alpha(\infty)}-1}{\frac{\alpha(0)}{\alpha(\infty)}-1}\right]^{2}-\ln\left[\frac{\left(\frac{\alpha(t)}{\alpha(\infty)}\right)^{2}+\frac{\alpha(t)}{\alpha(\infty)}+1}{\left(\frac{\alpha(0)}{\alpha(\infty)}\right)^{2}+\frac{\alpha(0)}{\alpha(\infty)}+1}\right]=-\frac{t}{\tau_{0}}$
$\displaystyle+2\sqrt{3}\arctan\left[\frac{1}{\sqrt{3}}+\frac{2}{\sqrt{3}}\frac{\alpha(t)}{\alpha(\infty)}\right]-2\sqrt{3}\arctan\left[\frac{1}{\sqrt{3}}+\frac{2}{\sqrt{3}}\frac{\alpha(0)}{\alpha(\infty)}\right]$
(46)

where $\tau_{0}$ is the characteristic time for the dissociation of a divalent
oxygen vacancy-electron pair and is given by

$\displaystyle\tau_{0}$ $\displaystyle=$
$\displaystyle\frac{\rho_{+o}\rho_{-o}^{2}}{24K_{A0}\alpha^{2}(\infty)}=\frac{1}{6}\left[\frac{\rho_{+o}\rho_{-o}^{2}}{4}\right]^{1/3}\frac{1}{K_{A0}^{1/3}K_{D0}^{2/3}}$
(47)

Solution of Eq. 46 is plotted in Fig. 1, which shows that for $t/\tau_{0}>10$,
the dissociation of vacancy-electron pair is almost complete i.e.,
$\alpha(t)/\alpha(\infty)\rightarrow 1$. This implies that the model developed
here can be simplified to some extent by considering only dissociated
vacancies and electrons for capturing kinetic effects at times much greater
than $\tau_{0}$. For example, construction of polarization-electrostatic
potential loop can be constructed without considering the vacancy-electron
pairs at time scales much greater than $\tau_{0}$, typical of experimental
rates at which voltages are sweeped across ferroelectric filmskelley2022oxygen
. But our complete formalism also allows us to describe ultra-fast timescales
within the same framework with relative ease.

Figure 1: Solution of Eq. 46 is shown here for different initial conditions.

###  III.2 Steady State Analysis: Implications of the Vegard’s law

Spatiotemporal responses of oxygen vacancies and electrons to the
electrostatic potential and the strain encoded in Eqs. LABEL:eq:phiplusfinal
and LABEL:eq:phinegfinal are asymmetric due to the fact that $z_{+}=2$ and the
use of the Vegard’s law so that only the vacancies affect the strain. In order
to make this clearer, the functional derivatives (exchange chemical
potentials) appearing in Eqs. LABEL:eq:phiplusfinal and LABEL:eq:phinegfinal
are evaluated as

$\displaystyle\left[\frac{\mu_{+}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]$
$\displaystyle=$
$\displaystyle\ln\left(\frac{\phi_{+}(\mathbf{r},t)}{1-\phi_{+}(\mathbf{r},t)-\phi_{-}(\mathbf{r},t)}\right)+\frac{z_{+}e\psi(\mathbf{r},t)}{k_{B}T}+\frac{G_{+o}-G_{\pm
o}}{k_{B}T}$
$\displaystyle-\frac{w_{ij}^{v}\sigma_{ij}(\mathbf{r},t)}{k_{B}T}-\frac{\kappa_{+}\rho_{+o}}{k_{B}T}\nabla^{2}\phi_{+}(\mathbf{r},t)-\frac{\kappa_{\pm
o}\rho_{\pm
o}}{k_{B}T}\nabla^{2}(\phi_{+}(\mathbf{r},t)+\phi_{-}(\mathbf{r},t))$
$\displaystyle\left[\frac{\mu_{-}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]$
$\displaystyle=$
$\displaystyle\ln\left(\frac{\phi_{-}(\mathbf{r},t)}{1-\phi_{+}(\mathbf{r},t)-\phi_{-}(\mathbf{r},t)}\right)-\frac{e\psi(\mathbf{r},t)}{k_{B}T}+\frac{G_{-o}-G_{\pm
o}}{k_{B}T}$ (49)
$\displaystyle-\frac{\kappa_{+}\rho_{+o}}{k_{B}T}\nabla^{2}\phi_{+}(\mathbf{r},t)-\frac{\kappa_{\pm
o}\rho_{\pm
o}}{k_{B}T}\nabla^{2}(\phi_{+}(\mathbf{r},t)+\phi_{-}(\mathbf{r},t))$

Note here that local stress-dependent term only appears in Eq.
LABEL:eq:vacancyflux and results from our assumption that the vacancies affect
the local strain, which appear in Eq. 19. In here, we show that the local
polarization depends on spatial distribution of the vacancies and electrons
via the electrostatic potential and doesn’t depend solely on the electric
field as in the classical ferroelectrics. First, we consider local
equilibrium, which correspond to a steady state of the time-dependent
equations. The local equilibrium is defined as the conditions
$\left[\frac{\mu_{+}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]=0$
and
$\left[\frac{\mu_{-}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]=0$.

At the steady state and representing local equilibrium,
$\mathbf{P}^{\star}(\mathbf{r},t)=0,p=\mbox{constant}$ and hence, leads to
$\nabla\cdot\mathbf{\sigma}(\mathbf{r},t)=0$. Representing all of the
variables (independent of time, $t$ at $t\rightarrow\infty$) in the steady
state by subscript $s$, we get
$\sigma_{ij}(\mathbf{r},t)\equiv\sigma_{ij,s}(\mathbf{r})=0$, where the latter
equality represents equilibrium (i.e., a stress-free state). This, in turn,
implies that
$\epsilon_{ij}(\mathbf{r},t)\equiv\epsilon_{ij,s}(\mathbf{r})=\epsilon_{ij,s}^{0}(\mathbf{r})=Q_{ijkl}P_{k,s}(\mathbf{r})P_{l,s}(\mathbf{r})+w_{ij}^{v}\rho_{+o}\phi_{+,s}(\mathbf{r})$.
For a weakly inhomogeneous distribution of the vacancies and the electrons so
that the derivative terms in Eqs. LABEL:eq:vacancyflux and 49 can be ignored,
$\left[\frac{\mu_{+}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]=0$
and
$\left[\frac{\mu_{-}(\mathbf{r},t)}{k_{B}T}-\frac{\mu_{\pm}(\mathbf{r},t)}{k_{B}T}\right]=0$
lead to

$\displaystyle\phi_{+,s}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp\left[\frac{ez_{+}\psi_{s}(\mathbf{r})-w_{ij}^{v}\sigma_{ij,s}(\mathbf{r})}{k_{B}T}+\frac{G_{+o}-G_{\pm
o}}{k_{B}T}\right]\left[1+\exp\left[\frac{e\left[\psi_{s}(\mathbf{r})-\left\\{G_{-o}-G_{\pm
o}\right\\}\right]}{k_{B}T}\right]\right]}$ (50)
$\displaystyle\phi_{-,s}^{-}(\mathbf{r})$ $\displaystyle=$
$\displaystyle\frac{1}{1+\exp\left[-\frac{e\left[\psi_{s}(\mathbf{r})-\left\\{G_{-o}-G_{\pm
o}\right\\}\right]}{k_{B}T}\right]\left[1+\exp\left[-\frac{\left[ez_{+}\psi_{s}(\mathbf{r})-w_{ij}^{v}\sigma_{ij,s}(\mathbf{r})\right]}{k_{B}T}-\frac{G_{+o}-G_{\pm
o}}{k_{B}T}\right]\right]}$ (51)

where we have used the notation
$\psi(\mathbf{r},t\rightarrow\infty)=\psi_{s}(\mathbf{r})$ and
$\sigma_{ij}(\mathbf{r},t\rightarrow\infty)=\sigma_{ij,s}(\mathbf{r})$. Using
these equations, it is clear that the effects of vacancies on the total strain
appear via the Vegard strain and leads to a result (valid at the steady state)

$\displaystyle\epsilon_{ij,s}(\mathbf{r})=Q_{ijkl}P_{k,s}(\mathbf{r})P_{l,s}(\mathbf{r})+\frac{w_{ij}^{v}\rho_{+o}}{1+\exp\left[\frac{ez_{+}\psi_{s}(\mathbf{r})-w_{ij}^{v}\sigma_{ij,s}(\mathbf{r})}{k_{B}T}+\frac{G_{+o}-G_{\pm
o}}{k_{B}T}\right]\left[1+\exp\left[\frac{e\left[\psi_{s}(\mathbf{r})-\left\\{G_{-o}-G_{\pm
o}\right\\}\right]}{k_{B}T}\right]\right]}$

A similar result without any consideration of the self-energy terms ( i.e.,
without $G_{io}$) was derived in Ref. kelley2022oxygen . An additional effect
of the vacancies is to affect the local electric field. At the local
equilibrium,
$\mathbf{P}_{s}(\mathbf{r})=\alpha(\mathbf{r})\mathbf{E}(\mathbf{r})$, where
the prefactor $\alpha(\mathbf{r})$ depends on the specific form of $H$ and can
be determined numerically for any functional form of $H$. As
$\mathbf{E}(\mathbf{r})=\mathbf{E}_{0}(\mathbf{r})+\mathbf{E}_{1}(\mathbf{r})$,
where
$\nabla\cdot\left[\left\\{\epsilon_{0}\epsilon_{\infty}+\alpha(\mathbf{r})\right\\}\mathbf{E}_{0}(\mathbf{r})\right]=0$
and
$\nabla\cdot\left[\left\\{\epsilon_{0}\epsilon_{\infty}+\alpha(\mathbf{r})\right\\}\mathbf{E}_{1}(\mathbf{r})\right]=ez_{+}\rho_{+o}\phi_{+,s}(\mathbf{r})-e\rho_{-o}\phi_{-,s}(\mathbf{r})$.
In other words, local electric field has contribution, $\mathbf{E}_{1}$,
resulting from inhomogeneous distribution of vacancies and oppositely charged
carriers. It should be noted that $\mathbf{E}_{1}=\mathbf{E}_{0}$ in the
absence of the vacancies but $\mathbf{E}_{1}\neq\mathbf{E}_{0}$ in the
presence of the charged vacancies. This, in turn, implies that the local
polarization
$\mathbf{P}_{s}(\mathbf{r})=\alpha(\mathbf{r})\left[\mathbf{E}_{0}(\mathbf{r})+\mathbf{E}_{1}(\mathbf{r})\right]$
is intimately connected with spatial distribution of vacancies and electrons.
In summary, the strain-electrostatic potential loop will be asymmetric with
respect to the sign of the electrostatic potential, in qualitative agreements
with recent experimentskelley2022oxygen .

###  III.3 Coupling between the fast and the slow components of the
polarization: topological effects in vacancy-free regime

The model developed here is based on the decomposition of the net local
polarization into a slowly-varying component, $\mathbf{P}$ and a fast
component, $\mathbf{P}_{e}$. A coupling between these two components appear in
the form of $\epsilon_{\infty}$ in the model, which affects the electrostatic
potential, $\psi(\mathbf{r},t)$. In addition, another coupling appears in Eq.
25 in the form of the lattice velocity, $\mathbf{v}_{l}$, which can be
interpreted in terms of the rate of the change of the net local polarization.
Molecular origin of this interpretation is the fact that local displacement of
electrons and ions of ferroelectric crystals contribute towards the net
polarization, which are treated in the model as $\mathbf{P}_{e}$ and
$\mathbf{P}$, respectively. This, in turn, leads to the relation
$\mathbf{v}_{l}(\mathbf{r},t)\sim\partial(\mathbf{P}_{e}(\mathbf{r},t)+\mathbf{P}(\mathbf{r},t))/\partial
t$. In the most of the phase field models, $\mathbf{v}_{l}(\mathbf{r},t)$ is
taken to be zero, which implies that
$\mathbf{v}_{l}(\mathbf{r},t)\simeq\partial(\mathbf{P}_{e}(\mathbf{r},t))/\partial
t\rightarrow 0$ at the time-scales relevant to the models. This implies that
Eq. 25 can be written as

$\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}=\mathbf{P}^{\star}(\mathbf{r},t),$ (53)

where Eqs. 17 and 25 allow us to identify

$\displaystyle\mathbf{P}^{\star}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau_{p}}\left[\frac{\delta\left[\int
d\mathbf{r}^{\prime}\,H_{LGD}\left\\{\mathbf{P}(\mathbf{r}^{\prime},t)\right\\}+H_{grad}\left\\{\nabla\mathbf{P}(\mathbf{r}^{\prime},t)\right\\}\right]}{\delta\mathbf{P}(\mathbf{r},t)}-\mathbf{E}(\mathbf{r},t))\right]$
(54)

In general, $H_{LGD}+H_{grad}$ can be written in powers of $\mathbf{P}$ so
that

$\displaystyle
H_{LGD}\left\\{\mathbf{P}\right\\}+H_{grad}\left\\{\nabla\mathbf{P}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{2\chi}\mathbf{P}^{2}(\mathbf{r},t)+\frac{\kappa_{p,1}}{2}\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]^{2}+\frac{\kappa_{p,2}}{2}\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]^{2}$
(55)
$\displaystyle+\kappa_{p,3}\mathbf{P}(\mathbf{r},t)\cdot\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]$

Eqs. 54 and 55 lead to

$\displaystyle\mathbf{P}^{\star}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau_{p}}\left[\frac{1}{\chi}\mathbf{P}(\mathbf{r},t)-\kappa_{p,1}\nabla\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]+\kappa_{p,2}\nabla\times\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]+2\kappa_{p,3}\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]\right.$
(56) $\displaystyle\left.-\mathbf{E}(\mathbf{r},t)\right],$

Using Eq. 23 in the absence of the vacancies and the electrons; and
$\nabla\times\mathbf{E}(\mathbf{r},t)=0$,

$\displaystyle\nabla\cdot\mathbf{P}^{\star}(\mathbf{r},t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau_{p}}\left[\left\\{\frac{1}{\chi}+\frac{1}{\epsilon_{0}\epsilon_{\infty}}\right\\}\nabla\cdot\mathbf{P}(\mathbf{r},t)-\kappa_{p,1}\nabla^{2}\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]\right]$
(57) $\displaystyle\nabla\times\mathbf{P}^{\star}(\mathbf{r},t)$
$\displaystyle=$
$\displaystyle-\frac{1}{\tau_{p}}\left[\frac{1}{\chi}\nabla\times\mathbf{P}(\mathbf{r},t)+\kappa_{p,2}\nabla\times\left(\nabla\times\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]\right)+2\kappa_{p,3}\nabla\times\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]\right]$

Operating with divergence and curl on Eq. 53, and using Eqs. 57,
LABEL:eq:curlPstar, we get

$\displaystyle\frac{\partial\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]}{\partial
t}$ $\displaystyle=$
$\displaystyle\left[\frac{\kappa_{p,1}}{\tau_{p}}\nabla^{2}-\frac{1}{\tau_{L}}\right]\left[\nabla\cdot\mathbf{P}(\mathbf{r},t)\right]$
(59)
$\displaystyle\frac{\partial\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]}{\partial
t}$ $\displaystyle=$
$\displaystyle\left[\frac{\kappa_{p,2}}{\tau_{p}}\nabla^{2}-\frac{2\kappa_{p,3}}{\tau_{p}}\left(\nabla\times\right)-\frac{1}{\tau_{p}\chi}\right]\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]$
(60)

where
$\tau_{L}=\tau_{p}/\left\\{\frac{1}{\chi}+\frac{1}{\epsilon_{0}\epsilon_{\infty}}\right\\}$
is the characteristic time for the change of
$\nabla\cdot\mathbf{P}(\mathbf{r},t)$ and $\tau_{p}\chi$ is the chracteristic
time for the change of $\nabla\times\mathbf{P}(\mathbf{r},t)$. Note that
$\tau_{L}/(\tau_{p}\chi)=1/(1+\chi/(\epsilon_{0}\epsilon_{\infty}))\ll 1$ for
$\chi\gg 1$. This means that for $t\gg\tau_{L}$,
$\nabla\cdot\mathbf{P}(\mathbf{r},t)=0$. This allows us to identify relation
between the local polarization and the electric field by using Eq. 53. For
example, in a steady state, $\mathbf{P}^{\star}(\mathbf{r},t)=0$, which leads
to a relation between the local polarization and the local electric field as

$\displaystyle\frac{1}{\chi}\mathbf{P}_{s}(\mathbf{r})-\kappa_{p,1}\nabla\left[\nabla\cdot\mathbf{P}_{s}(\mathbf{r})\right]+\kappa_{p,2}\nabla\times\left[\nabla\times\mathbf{P}_{s}(\mathbf{r})\right]+2\kappa_{p,3}\left[\nabla\times\mathbf{P}_{s}(\mathbf{r})\right]$
$\displaystyle=$ $\displaystyle\mathbf{E}_{s}(\mathbf{r}),$ (61)

where we have used the notation
$\mathbf{P}_{s}(\mathbf{r})=\mathbf{P}(\mathbf{r},\infty)$ and
$\mathbf{E}_{s}(\mathbf{r})=\mathbf{E}(\mathbf{r},\infty)$. For a given
volume, solution of Eq. 61 depends on the geometry and boundary conditions. A
particular set of solutions, which enforces
$\nabla\cdot\mathbf{P}_{s}(\mathbf{r})=0$ everywhere in space including the
boundaries, will be discussed here. In particular,
$\nabla\times\mathbf{P}_{s}(\mathbf{r})=\lambda\mathbf{P}_{s}(\mathbf{r})$,
which enforces $\nabla\cdot\mathbf{P}_{s}(\mathbf{r})=0$ for a constant
$\lambda$ will be discussed here. For such a divergence-free polarization
vector, Eq. 61 demands

$\displaystyle\left[\frac{1}{\chi}+\kappa_{p,2}\lambda^{2}+2\kappa_{p,3}\lambda\right]\mathbf{P}_{s}(\mathbf{r})$
$\displaystyle=$ $\displaystyle\mathbf{E}_{s}(\mathbf{r}),$ (62)

$\lambda$ can be obtained by solving for an eigenvalue of the equation
$\nabla\times\mathbf{P}_{s}(\mathbf{r})=\lambda\mathbf{P}_{s}(\mathbf{r})$
with boundary conditions. Nevertheless, Eq. 62 shows that
$\mathbf{P}_{s}(\mathbf{r})$ and $\mathbf{E}_{s}(\mathbf{r})$ are parallel to
each other at the steady state, which may not be true in general. In fact, Eq.
53 can be used to determine the time-dependent electric field required for
setting up known polarization vector in space and time.

In contrast to the above analysis, if we consider
$\mathbf{v}_{l}(\mathbf{r},t)=\beta\partial\mathbf{P}(\mathbf{r},t)/\partial
t,\beta$ being a constant, then Eq. 25 can be written as

$\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}+\beta\left[\frac{\partial\mathbf{P}(\mathbf{r},t)}{\partial
t}\cdot\nabla\right]\mathbf{P}(\mathbf{r},t)+\frac{\beta}{2}\mathbf{P}(\mathbf{r},t)\times\frac{\partial\left[\nabla\times\mathbf{P}(\mathbf{r},t)\right]}{\partial
t}=\mathbf{P}^{\star}(\mathbf{r},t),$ (63)

Eq. 63 shows that the coupling between the slow and the fast component of the
net polarization will affect spatiotemporal distribution of
$\mathbf{P}(\mathbf{r},t)$, while the steady state behavior of
$\mathbf{P}(\mathbf{r},t)\equiv\mathbf{P}_{s}(\mathbf{r})$ remains the same.
Eq. 63 opens up a way to study effects of confinement on the spatiotemporal
distribution of topologically non-trivial polarization, which may not be
realized at an equilibrium.

##  IV Conclusions

We developed a thermodynamically consistent time-dependent model for
understanding the effects of multivalent vacancies on relations among
polarization-electric potential and strain-electric potential in thin films of
ferroelectrics. In contrast to the most of the phase field models, non-linear
effects of the reaction kinetics leading to generation of charged vacancies
and electrons from their pairs are introduced in the model. In addition,
diffusion and elastic effects of the charged vacancies are considered, which
are shown to exhibit asymmetric reponses of the strain to the electric
potential. Furthermore, the model introduces coupling between the slow
component and a fast component of the net polarization, which is expected to
affect time-dependent relation between the polarization and the electric
field. Impedance response from the thin films of ferroelectrics with vacancies
should exhibit defect dipole behaviording_influence_2018 ; kelley2022oxygen ;
akbarian2019understanding ; DHAKANE2023100264  which will dynamically couple
with the polar matrix giving rise to frequency-dependent behavior. In the
presence of electrodes, one can also expect electrode polarization
phenomenonkumar2017  due to localization of the vacancies near an oppositely
charged electrode. These localized vacancy induced polarization phenomenon can
be used to extract diffusion constant of the vacancies. A simplified
modelkumar2017  developed using the Rayleighian approach has been used
previously to fit the impedance spectra as a function of frequency and
temperature for thin films of ionic polymers. We envision that in future,
experimental results for impedance spectra from thin films of ferroelectrics
with vacancies can be fitted using the model developed here and extract the
diffusion constant of the vacancies. The Rayleighian approach to build the
model is shown to be general enough for constructing non-linear reaction
kinetics. Application of the model to simulate motion of domain walls in
ferroelectrics in the presence of vacancies will be presented in a forthcoming
publication.

##  V Acknowledgements

This work was supported by the Center for Nanophase Materials Sciences, which
is a US DOE, Office of Science User Facility at Oak Ridge National Laboratory.

##  VI Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * [1] Masatoshi Imada, Atsushi Fujimori, and Yoshinori Tokura.  Metal-insulator transitions.  Rev. Mod. Phys., 70:1039–1263, Oct 1998.
  * [2] David D Awschalom and Michael E Flatté.  Challenges for semiconductor spintronics.  Nature physics, 3(3):153–159, 2007.
  * [3] Karin M Rabe, Matthew Dawber, Céline Lichtensteiger, Charles H Ahn, and Jean-Marc Triscone.  Modern physics of ferroelectrics: Essential background.  Physics of Ferroelectrics, pages 1–30, 2007.
  * [4] S. Oh, H. Hwang, and I. K. Yoo.  Ferroelectric materials for neuromorphic computing.  APL Materials, 7(9), 09 2019.  091109\.
  * [5] Jaewoo Jeong, Nagaphani Aetukuri, Tanja Graf, Thomas D. Schladt, Mahesh G. Samant, and Stuart S. P. Parkin.  Suppression of metal-insulator transition in $\mbox{VO}_{2}$ by electric field induced oxygen vacancy formation.  Science, 339(6126):1402–1405, 2013.
  * [6] Kyle P Kelley, Anna N Morozovska, Eugene A Eliseev, Vinit Sharma, Dundar E Yilmaz, Adri CT van Duin, Panchapakesan Ganesh, Albina Borisevich, Stephen Jesse, Peter Maksymovych, et al.  Oxygen vacancy injection as a pathway to enhancing electromechanical response in ferroelectrics.  Advanced Materials, 34(2):2106426, 2022.
  * [7] Sang Mo Yang, Anna N Morozovska, Rajeev Kumar, Eugene A Eliseev, Ye Cao, Lucie Mazet, Nina Balke, Stephen Jesse, Rama K Vasudevan, Catherine Dubourdieu, et al.  Mixed electrochemical–ferroelectric states in nanoscale ferroelectrics.  Nature Physics, 13(8):812–818, 2017.
  * [8] Dooman Akbarian, Dundar E Yilmaz, Ye Cao, Panchapakesan Ganesh, Ismaila Dabo, Jason Munro, Renee Van Ginhoven, and Adri CT Van Duin.  Understanding the influence of defects and surface chemistry on ferroelectric switching: a reaxff investigation of batio 3.  Physical Chemistry Chemical Physics, 21(33):18240–18249, 2019.
  * [9] Yu Xiao and Kaushik Bhattacharya.  A continuum theory of deformable, semiconducting ferroelectrics.  Archive for Rational Mechanics and Analysis, 189:59–95, 2008.
  * [10] Jun Liang Lin, Ri He, Zengxing Lu, Yi Lu, Zhiming Wang, Zhicheng Zhong, Xiang Zhao, Run-Wei Li, Zhi Dong Zhang, and Zhan Jie Wang.  Oxygen vacancy enhanced ferroelectricity in bto: Sro nanocomposite films.  Acta Materialia, 199:9–18, 2020.
  * [11] YL Li, SY Hu, ZK Liu, and LQ Chen.  Phase-field model of domain structures in ferroelectric thin films.  Applied Physics Letters, 78(24):3878–3880, 2001.
  * [12] Atanu K Saha, Kai Ni, Sourav Dutta, Suman Datta, and Sumeet Gupta.  Phase field modeling of domain dynamics and polarization accumulation in ferroelectric hzo.  Applied Physics Letters, 114(20):202903, 2019.
  * [13] Jian-Jun Wang, Bo Wang, and Long-Qing Chen.  Understanding, predicting, and designing ferroelectric domain structures and switching guided by the phase-field method.  Annual Review of Materials Research, 49:127–152, 2019.
  * [14] Long-Qing Chen.  Appendix a–landau free-energy coefficients.  Physics of ferroelectrics: a modern perspective, pages 363–372, 2007.
  * [15] Yihui Zhang, Jiangyu Li, and Daining Fang.  Oxygen-vacancy-induced memory effect and large recoverable strain in a barium titanate single crystal.  Physical Review B, 82(6):064103, 2010.
  * [16] Ye Cao, Jie Shen, Clive Randall, and Long-Qing Chen.  Effect of ferroelectric polarization on ionic transport and resistance degradation in batio3 by phase-field approach.  Journal of the American Ceramic Society, 97(11):3568–3575, 2014\.
  * [17] Yasuhide Shindo, Fumio Narita, and Takuya Kobayashi.  Phase field simulation on the electromechanical response of poled barium titanate polycrystals with oxygen vacancies.  Journal of Applied Physics, 117(23):234103, 2015.
  * [18] Patrick Fedeli, Marc Kamlah, and Attilio Frangi.  Phase-field modeling of domain evolution in ferroelectric materials in the presence of defects.  Smart Materials and Structures, 28(3):035021, 2019.
  * [19] R. A. Marcus.  Electrostatic Free Energy and Other Properties of States Having Nonequilibrium Polarization. I.  The Journal of Chemical Physics, 24(5):979–989, 12 2004.
  * [20] Rajeev Kumar, Jyoti P. Mahalik, Vera Bocharova, Eric W. Stacy, Catalin Gainaru, Tomonori Saito, Mallory P. Gobet, Steve Greenbaum, Bobby G. Sumpter, and Alexei P. Sokolov.  A Rayleighian approach for modeling kinetics of ionic transport in polymeric media.  The Journal of Chemical Physics, 146(6), 02 2017.  064902\.
  * [21] A. R. Denton and N. W. Ashcroft.  Vegard’s law.  Phys. Rev. A, 43:3161–3164, Mar 1991.
  * [22] Lars Onsager.  Reciprocal relations in irreversible processes. i.  Phys. Rev., 37:405–426, Feb 1931.
  * [23] S.R. Groot and P. Mazur.  Non-equilibrium Thermodynamics.  Dover books on physics and chemistry. North-Holland Publishing Company, 1962.
  * [24] L. Onsager and S. Machlup.  Fluctuations and irreversible processes.  Phys. Rev., 91:1505–1512, Sep 1953.
  * [25] Masao Doi, Jiajia Zhou, Yana Di, and Xianmin Xu.  Application of the onsager-machlup integral in solving dynamic equations in nonequilibrium systems.  Phys. Rev. E, 99:063303, Jun 2019.
  * [26] Masao Doi.  Onsager’s variational principle in soft matter.  Journal of Physics: Condensed Matter, 23(28):284118, jun 2011.
  * [27] John Cummings, John S. Lowengrub, Bobby G. Sumpter, Steven M. Wise, and Rajeev Kumar.  Modeling solvent evaporation during thin film formation in phase separating polymer mixtures.  Soft Matter, 14:1833–1846, 2018.
  * [28] J Hubbard and L Onsager.  Dielectric dispersion and dielectric friction in electrolyte solutions. i.  The Journal of chemical physics, 67(11):4850–4857, 1977.
  * [29] BU Felderhof and HJ Kroh.  Hydrodynamics of magnetic and dielectric fluids in interaction with the electromagnetic field.  The Journal of chemical physics, 110(15):7403–7411, 1999.
  * [30] L. D. Landau, L. P. Pitaevskii, E. M. Lifshitz, and A. M. Kosevich.  Theory of Elasticity.  Butterworth-Heinemann, 3 edition, 1986.
  * [31] Toshio Mura and Toshio Mura.  General theory of eigenstrains.  Micromechanics of defects in solids, pages 1–73, 1987.
  * [32] Long-Qing Chen.  Phase-field models for microstructure evolution.  Annual review of materials research, 32(1):113–140, 2002.
  * [33] J. Israelachvili.  Intermolecular and Surface Forces.  Academic, London, 1992.
  * [34] Rajeev Kumar, Bobby G. Sumpter, and M. Muthukumar.  Enhanced phase segregation induced by dipolar interactions in polymer blends.  Macromolecules, 47(18):6491–6502, 2014.
  * [35] S. Glasstone, K. Laidler, and H. Eyring.  The theory of rate processes.  McGraw Hill, New York, 1941.
  * [36] Rajeev Kumar, Zening Liu, Brad Lokitz, Jihua Chen, Jan-Michael Carrillo, Jacek Jakowski, C Patrick Collier, Scott Retterer, and Rigoberto Advincula.  Harnessing autocatalytic reactions in polymerization and depolymerization.  MRS Communications, 11:377–390, 2021.
  * [37] John Newman and Nitash P Balsara.  Electrochemical systems.  John Wiley & Sons, 2021.
  * [38] Jilai Ding, Janakiraman Balachandran, Xiahan Sang, Wei Guo, Jonathan S. Anchell, Gabriel M. Veith, Craig A. Bridges, Yongqiang Cheng, Christopher M. Rouleau, Jonathan D. Poplawsky, Nazanin Bassiri-Gharb, Raymond R. Unocic, and P. Ganesh.  The Influence of Local Distortions on Proton Mobility in Acceptor Doped Perovskites.  Chemistry of Materials, 30(15):4919–4925, August 2018.  Publisher: American Chemical Society.
  * [39] Abhijeet Dhakane, Tian Xie, Dundar E. Yilmaz, Adri C.T. van Duin, Bobby G. Sumpter, and P. Ganesh.  A graph dynamical neural network approach for decoding dynamical states in ferroelectrics.  Carbon Trends, 11:100264, 2023.


Web-based telluric correction made in Spain: spectral fitting of Vega-type
telluric standards
de la Fuente, D.1,2, Marco, A.1, Patrick, L. R.1,3,4, Rübke, K.1, López, I.1,
Fernández, A.1, Conejero, S.1, Navarro, J.1, Palazón, M.1, and Negueruela, I.4
1Departamento de Física, Ingeniería de Sistemas y Teoría de la Señal,
Universidad de Alicante, San Vicente del Raspeig, Spain

2CREOL, The College of Optics and Photonics, University of Central Florida,
Orlando, FL, USA

3Departamento de Astrofísica, Centro de Astrobiología (CSIC-INTA), Torrejón de
Ardoz, Spain

4Departamento de Física Aplicada, Universidad de Alicante, San Vicente del
Raspeig, Spain

## Abstract

Infrared spectroscopic observations from the ground must be corrected from
telluric contamination to make them ready for scientific analyses. However,
telluric correction is often a tedious process that requires significant
expertise to yield accurate results in a reasonable time frame. To solve these
inconveniences, we present a new method for telluric correction that employs a
roughly simultaneous observation of a Vega analog to measure atmospheric
transmission. After continuum reconstruction and spectral fitting, the stellar
features are removed from the observed Vega-type spectrum and the result is
used for cancelling telluric absorption features on science spectra. This
method is implemented as TelCorAl (Telluric Correction from Alicante), a
Python-based web application with a user-friendly interface, whose beta
version will be released soon.

##  1 Introduction

Infrared spectra that are observed from ground-based facilities are heavily
contaminated by absorption features from the Earth’s atmosphere. The accuracy
of the decontamination process, known as telluric correction, is crucial for
producing reliable results from subsequent spectroscopic analyses.

Telluric correction basically consists in dividing the science spectra by
atmospheric transmission. However, the latter cannot be easily determined
because the atmosphere behaves in a complex and unpredictable manner. The
intensities and shapes of telluric absorptions depend not only on
observational settings (airmass, instrumental profile, the observer’s
altitude) but also on atmospheric properties that can change dramatically
(abundances for the involved molecules, precipitable water vapor, seeing).

The currently available methods mainly differ on how this transmission
spectrum is obtained. On the one hand, theoretical methods (e.g. [11]) fit a
transmission model to the observations. On the other hand, observational
methods (e.g. [12]) extract the transmission spectrum from a roughly
simultaneous observation of a standard star whose intrinsic features are well
known and easily separated from telluric absorptions. There also are some
hybrid approaches (e.g. [2]). The user often has to make a complex choice,
based on instrumental and operational features (spectral resolution, exposure
times, observation schedule, etc.), and other requirements related to the
planned treatment (intended signal-to-noise ratio, integration into a
pipeline, etc.).

The majority of publicly available algorithms are designed for specific
instruments or tailored to the scientific goals of some particular team. A few
others have a more general purpose and are targeted to a wider community,
although they have practical disadvantages in common. For example, users often
complain about problems for installation and usability of these software, or
about the amount of experience needed to obtain accurate results in a
reasonable amount of time.

This situation has encouraged us to create a cross-platform, user-friendly
application for telluric correction. To avoid the complications of fitting sky
models or empirical libraries to the telluric component of science spectra, we
have chosen the approach of measuring atmospheric transmission on an Observed
Standard Star Spectrum (OSSS). We refer to [10] and [12] for the theoretical
background for this kind of methods.

##  2 Methodology

Vega-type standards (i.e. spectral types around A0V) are preferred, as their
intrinsic spectra only display wide hydrogen absorptions in most cases (see
[9] for a notable counterexample). These lines are usually well described by
Voigt functions (see [8] and references therein), which facilitates the
modeling and removal of the OSSS stellar features. Model fitting, however,
cannot be directly performed within wavelength regions where the continuum
disappears from the observation (owing to broad telluric bands or to
overlapping high-order hydrogen lines).

The extraction of the atmospheric transmission spectrum in such a way is the
first phase of telluric correction, and is carried out in five steps:

Figure 1:  Spectra involved in the process of fitting a spectrum of a telluric
standard star, as well as the normalized version of the OSSS for comparison
purposes.

  1. 1.

A “default telluric correction” of the standard star is made by dividing the
OSSS by a synthetic trasmission spectrum (respectively shown in black and
light blue in Fig. 1). To generate the latter, the [4] telluric models for
average atmospheric conditions are interpolated at the airmass of the OSSS,
convolved with a Gaussian profile that matches the spectral resolution of the
instrument, and resampled into the OSSS wavelength axis. The [4] models also
provide variability information that is used for computing the uncertainty of
the default correction.

  2. 2.

A simple model of hydrogen absorption spectrum is created as a combination of
Voigt profiles whose widths and relative intensities are preset for the
spectral type of the OSSS. The profile heights are then multiplied by a
scaling constant to match the line heights of the default telluric correction
(brown line in Fig. 1), thus obtaining a first model of the standard star
(dashed orange line).

  3. 3.

The ratio between the results of previous steps (i.e. brown line divided by
orange line) is used as an approximation for the continuum level (shown in
green in Fig. 1). Weighted least-squares polynomial fits with boundary
conditions are performed for increasing degrees of the polynomial, and the
optimal result (magenta dashed line) is chosen on the basis of a goodness-of-
fit test.

  4. 4.

The spectrum that results from the default correction is normalized (i.e.
brown/magenta = red in Fig. 1) and used for fitting a model for the
combination of all H line profiles. The result is multiplied by the continuum
fit (magenta dashed line) to produce a full model for the intrinsic features
of the OSSS (See upper graph in Fig. 2).

  5. 5.

The telluric transmission spectrum (Fig. 2, lower plot) is obtained as the
ratio of the OSSS to its full model.

Figure 2:  TelCorAl interface, showing an example of the fitting process for
intrinsic features of a telluric standard star.

If one or both of the initial synthetic spectra (like those in Fig. 1, top
pannel) diverges significantly from the OSSS features, the fitting procedures
will produce inaccurate results. In this case, the method will iterate on
steps 1–5 so that the model and transmission results of an iteration are used
as initial spectra for the next one.

Once an accurate result for telluric transmission is achieved, the second
phase consists in using the resulting spectrum for correcting the science
spectra. In an ideal case, the telluric features imprinted on all these
spectra would be identical, and could be perfectly removed by dividing each
science spectrum by the transmission. In practice, however, there are
intensity and shape variations, caused not only by variability of the sky (for
non-simultaneous observations), but also by factors related with the reduction
process (e.g. wavelength calibrations of co-added frames) or with implicit
approximations (e.g. the definition of airmass). To deal with these
differences, the following tasks are carried out for each science spectrum:

  * •

The transmission spectrum is slightly red- or blue-shifted according to a
factor that corresponds to the maximum cross-correlation with the science
spectrum. As far as possible, the cross-correlation function should be
computed in a wavelength range where telluric features dominate over the
intrinsic features of the science spectrum.

  * •

The ratio between airmass factors for the science and transmission spectra is
used for scaling the strength of telluric features through an exponential
relation ([7]). The result can be manually fine-tuned by the user until
optimal telluric cancellation is achieved.

##  3 Implementation and current status

We have developed TelCorAl (Telluric Correction from Alicante), a user-
friendly web application that performs telluric correction in a semiautomatic
way. Apart from carrying out all the telluric correction process as explained
in Sect. 2, the users may choose to perform just one phase, or both phases
independently, in order to test the method or make the correction in a more
customized way. In all cases, the users are required to upload their input
spectra in comma-separated value format, and the results are delivered in the
same format together with information on the parameters used throughout the
process. All calculations are internally made by a Python 3 code, and users
interact with it through a web interface that also displays helpful plots
(Fig. 2). Everything is processed in the RAM of the server, and no data or
results are ever stored outside it, thus keeping them private at all times.

The combined H line models are built and fitted through lmfit 1.0.3 ([6, 5]),
by employing a line profile function of the user’s choice. Based on our
experience, the pseudo-Voigt approximation ([14]) is usually the option that
works better in general, although Gaussian or Lorentzian profiles are also
offered, as they are more suitable under certain circunstances (e.g. extreme
values for resolution and $v\sin i$ of the standard star). Other critical
components of fitting and processing the spectra for telluric correction make
use of the SciPy ([13]), PyAstronomy ([1]), and specutils ([3]) libraries.

A preliminary version of TelCorAl is currently running on the server of our
group, although being only accessible from the University of Alicante network
for the moment. Once the final telluric correction code is fine-tuned and
comprehensive testing is carried out, a beta version will be released to the
public.

##  4 Concluding remarks

Our forthcoming web application solves issues that make other telluric
correction software impractical for widespread use. The user-friendly webpage
format, which does not require installation, along with the semiautomatic
approach, makes our software suitable even for inexperienced users that aim at
promptly getting a decent correction of their infrared data. TelCorAl is also
capable of dealing with a wide variety of instrumental and observational
features of the input data.

Compared to other telluric correction methods that are based on spectral
fitting of Vega-type stars ([12, 10]), our code uses a novel approach that
consists in reconstructing the continuum level (green line in Fig. 1). This
allows us to achieve accurate fits on regions where hydrogen line profiles
overlap with each other or with wide telluric bands (e.g. between 1.54 and
$1.62~{}\mu\mathrm{m}$, see Fig. 1).

## Acknowledgments

This work has been financially supported by Generalitat Valenciana through
grants APOSTD/2020/228, APOSTD/2020/247, and PROMETEO/2019/041. The Spanish
Ministerio de Ciencia e Innovación and Agencia Estatal de Investigación
(MCIN/AEI/10.13039/501100011033/FEDER, UE) has also provided financial support
through grants PGC2018-093741-B-C21 and PID2021-122397NB-C22.

## References

  * [1] Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, ascl:1906.010
  * [2] Cook, N. J., Artigau, É., Doyon, R., et al. 2022, arXiv:2211.01358
  * [3] Earl, N., Tollerud, E., Jones, C., et al. 2022, Zenodo
  * [4] Moehler, S., Modigliani, A., Freudling, W., et al. 2014, A&A 568, A9
  * [5] Newville, M., Otten, R., Nelson, A., et al. 2021, Zenodo
  * [6] Newville, M., Stensitzki, T., Allen, D. B., et al. 2014, Zenodo
  * [7] Perliski, L. M. & Solomon, S. 1993, JGR, 98, 10363
  * [8] Posener, D. W. 1959, AuJPh., 12, 184
  * [9] Sameshima, H., Ikeda, Y., Matsunaga, N., et al. 2018, ApJS, 239, 19
  * [10] Sameshima, H., Matsunaga, N., Kobayashi, N., et al. 2018, PASP, 130, 074502
  * [11] Smette, A., Sana, H. Noll, S. et al. 2015, A&A 576, A77
  * [12] Vacca, W. D., Cushing, M. C., & Rayner, J. T. 2003, PASP, 115, 389
  * [13] Virtanen P., Gommers R., Oliphant T. E., et al. 2020, NatMe, 17, 261
  * [14] Whiting, E. E. 1968, JQSRT, 8, 1379


# Using Machine Learning to predict Characteristics of Microstrip Line and
Microstrip Patch Antenna

Bharath Balaji
Department of Electronics and Communication Engineering
National Institute of Technology
Trichy
<EMAIL_ADDRESS>
&Dr. S. Raghavan
Department of Electronics and Communication Engineering
National Institute of Technology
Trichy
<EMAIL_ADDRESS>

###### Abstract

This study, conducted in 2017, explores the use of Machine learning algorithms
to predict Characteristics of Transmission Lines such as Impedance or
resonance frequency using design parameters of Transmission Lines. Using
formulas and equations that define the characteristics of Transmission lines,
training data was generated. We trained different models for this dataset. The
extent of deviation of predicted output from the actual output was measured in
terms of maximum error and average error. This helped determine how well an
algorithm worked for a particular transmission line. Further, the best-suited
algorithm for each transmission line under consideration was found based on
the error. 111This research was conducted in the year 2017.

##  1 Introduction

Machine learning has significantly advanced the field of microwave research,
particularly in the areas of modeling, simulation, and optimization. Notably,
the work by Naser-Moghaddasi et al. [5] demonstrated the application of
heuristic artificial neural networks to analyze and synthesize the performance
characteristics of rectangular microstrip antennas. Their approach highlighted
the potential of neural networks to predict antenna behavior effectively,
thereby streamlining the design process. Building upon such foundational
studies, this paper extends the use of machine learning to model a broader
range of transmission line characteristics. Specifically, we expand the
application of machine learning models to include various types of
transmission lines such as Microstrip, Slotline, Stripline, Co-Planar
Waveguide (CPW), Co-Planar Strip (CPS), and Microstrip Patch Antenna,
employing both Linear Regression and advanced Neural Networks. In this
publication, we explore the results for the Microstrip line and Microstrip
patch Antenna.

This research leverages standardized historical data to train models that can
predict the electrical properties of planar transmission lines, based on
physical dimension parameters. This method simplifies the practical design
process, enabling the prediction of output parameters for all combinations of
inputs before fabrication. Such predictive modeling saves considerable time
and resources in antenna design and simulation efforts. Furthermore, the
trained models can be utilized to predict the characteristics of new
transmission lines, eliminating the need for extensive empirical simulation.
This study aims not only to validate the effectiveness of machine learning
models in replicating known transmission line behaviors but also to enhance
their accuracy and generalizability compared to earlier works. Thus, this
research not only confirms the utility of machine learning in this domain but
also advances its capability to accommodate a wider array of transmission line
configurations and complexities.

##  2 Transmission Line Models

According to a Microstrip design proposed by K.C.Gupta and Ramesh Garg[1]: For
w/h > 1

$\epsilon_{eff}=\frac{\epsilon_{r}+1}{2}+\frac{\epsilon_{r}-1}{2}\left(\frac{1}{\sqrt{1+12\frac{H}{W}}}\right)$

$Z_{0}=\frac{120\pi}{\sqrt{\epsilon_{eff}}}\left(\frac{W}{H}+1.393+\frac{2}{3}\ln\left(\frac{W}{H}+1.444\right)\right)\Omega$

A proposed design of Microstrip Patch Antenna by Bablu Kumar Singh[3]:

$\epsilon_{\text{eff}}=\frac{\epsilon_{r}+1}{2}+\frac{\epsilon_{r}-1}{2}\left(1+12\frac{h}{W}\right)^{-\frac{1}{2}}$

$\Delta
L=0.412h\left(\frac{\epsilon_{\text{eff}}+0.3}{\epsilon_{\text{eff}}-0.258}\left(\frac{W}{h}+0.264\right)\right)$

$f_{r}=\frac{c}{2\sqrt{\epsilon_{\text{eff}}}(L+2\Delta L)}$

###  2.1 Microstrip Lines

Microstrip is a planar transmission line used to carry Electro-magnetic waves
(EM waves) or microwave frequency signals. It consists of 3 layers, conducting
strip, dielectric , and Ground plane. It is used to design and fabricate RF
and microwave components such as directional couplers, power
dividers/combiners, filters, antennae, MMIC, etc.

$\varepsilon_{r}$ for all the plots using ANN The dashed line implies
Predicted Output The dotted line implies Actual Output

###  2.2 Results for Microstrip

Figure 1: Impedance Vs. w/h for Microstrip Line using neural networks Figure 2:  Absolute Error Vs. w/h for Microstrip Line using neural networks Figure 3: Impedance Vs. w/h for Microstrip Line using linear regression Figure 4:  Absolute Error Vs. w/h for Microstrip Line using linear regression Table 1: Predicted Impedance values for Microstrip Line using neural networks and Linear Regression w/h | Impedance (actual) | neural networks | Linear Regression | % Error (neural networks) | % Error (Linear Regression)
---|---|---|---|---|---
1.000 | 98.525 | 97.323 | 71.658 | 1.220 | 27.269
1.500 | 80.819 | 81.450 | 68.265 | 0.781 | 15.533
2.000 | 68.774 | 68.488 | 64.871 | 0.416 | 5.675
2.500 | 59.999 | 60.280 | 61.477 | 0.468 | 2.463
3.000 | 53.296 | 53.129 | 58.084 | 0.313 | 8.984
3.500 | 47.993 | 48.268 | 54.690 | 0.573 | 13.954
4.000 | 43.686 | 43.408 | 51.297 | 0.636 | 17.422
4.500 | 40.113 | 39.950 | 47.903 | 0.406 | 19.420
5.000 | 37.098 | 37.376 | 44.509 | 0.749 | 19.977
5.500 | 34.517 | 34.802 | 41.116 | 0.826 | 19.118
6.000 | 32.281 | 32.228 | 37.722 | 0.164 | 16.855
7.000 | 30.325 | 29.801 | 34.329 | 1.728 | 13.204
7.500 | 28.598 | 28.517 | 30.935 | 0.283 | 8.172
8.000 | 27.061 | 27.233 | 27.541 | 0.636 | 1.774
8.500 | 25.685 | 25.949 | 24.148 | 1.028 | 5.984
9.000 | 24.445 | 24.665 | 20.754 | 0.900 | 15.099

###  2.3 Microstrip Patch Antenna

In this study, the rectangular microstrip antennas consist of a rectangular
patch of varying dimensions, namely width (W) and length (L), positioned over
a ground plane separated by a substrate with thickness (h). The substrate
material’s dielectric constant $\varepsilon_{r}$ plays a critical role in
determining the antenna’s properties. For the purpose of the analyses
presented, a dielectric constant $\varepsilon_{r}$ of 6 is used, which falls
within the typical range for microstrip antenna substrates
$2.2<\varepsilon_{r}<12$. All plots and results are derived using this
specified dielectric constant.

###  2.4 Results for Patch Antenna

Figure 5: Impedance Vs. w/h for Patch using neural networks Figure 6:  Absolute Error Vs. w/h for Patch using neural networks Figure 7: Impedance Vs. w/h for Patch using linear regression Figure 8:  Absolute Error Vs. w/h for Patch using linear regression Table 2: Frequency vs. w/h for Microstrip Line with Error Percentage w/h | Fr in MHz (actual) | neural networks | Linear Regression | % Error (neural networks) | % Error (Linear Regression)
---|---|---|---|---|---
1.0 | 7710.557 | 7692.914 | 7529.473 | 0.229 | 2.349
1.5 | 7585.017 | 7593.029 | 7491.309 | 0.106 | 1.235
2.0 | 7489.211 | 7493.143 | 7453.145 | 0.053 | 0.482
2.5 | 7411.944 | 7397.386 | 7414.982 | 0.196 | 0.041
3.0 | 7347.491 | 7348.007 | 7376.818 | 0.007 | 0.399
3.5 | 7292.474 | 7298.629 | 7338.654 | 0.084 | 0.633
4.0 | 7244.707 | 7251.310 | 7300.490 | 0.091 | 0.770
4.5 | 7202.688 | 7204.505 | 7262.326 | 0.025 | 0.828
5.0 | 7165.336 | 7157.707 | 7224.162 | 0.106 | 0.821
5.5 | 7131.843 | 7123.447 | 7185.999 | 0.118 | 0.759
6.0 | 7101.594 | 7100.492 | 7147.835 | 0.016 | 0.651
6.5 | 7074.102 | 7077.537 | 7109.671 | 0.049 | 0.503
7.0 | 7048.982 | 7054.582 | 7071.507 | 0.079 | 0.320
7.5 | 7025.921 | 7031.627 | 7033.343 | 0.081 | 0.106
8.0 | 7004.661 | 7008.672 | 6995.180 | 0.057 | 0.135
8.5 | 6984.987 | 6985.717 | 6957.016 | 0.010 | 0.400
9.0 | 6966.719 | 6962.762 | 6918.852 | 0.057 | 0.687
9.5 | 6949.706 | 6939.807 | 6880.688 | 0.142 | 0.993

##  3 Analysis

This study has systematically investigated various neural modeling techniques
for predicting the characteristics of transmission lines, specifically
examining their applicability and accuracy compared to traditional
computational methods. Our analysis utilized actual impedance and frequency
measurements from microstrip lines to evaluate the performance of neural
networks and linear regression algorithms. The results demonstrated that while
no single method uniformly excels across all scenarios, neural networks
generally provided closer approximations to actual measurements than linear
regression. For instance, in impedance predictions, neural networks showed a
maximum error of approximately 0.229% for a w/h ratio of 1.000, whereas linear
regression exhibited a significantly higher error of 2.349% under the same
conditions. Similarly, in frequency estimations, neural networks consistently
maintained lower error margins, indicating its superior ability to adapt to
the nonlinear dynamics of microstrip line characteristics. The traditional
approaches, often implemented through MATLAB toolboxes or similar simulation
software, while effective, are somewhat rigid, limiting their adaptability to
different antenna architectures. The models discussed in this paper offer a
more flexible and scalable alternative. By employing neural networks, we have
described synthesis and analysis methods that not only handle varying
transmission line parameters but also demonstrate scalability and adaptability
to diverse antenna designs.

##  4 Conclusion

Our findings underscore the importance of choosing the right model based on
the specific characteristics, requirements and architecture of the
transmission line. The errors observed in our simulations were within
acceptable limits, affirming the viability of using advanced machine learning
techniques for the design and analysis of transmission lines. These insights
pave the way for future research to further refine these models, potentially
integrating more complex[2][4] transmission lines and diverse datasets to
enhance predictive accuracy and reliability across a broader range of
transmission line applications.

## References

  * [1] K.C. Gupta and Ramesh Garg, _Microstrip Lines and Slotlines_ , Artech Inc., pp. 43-50, Jun. 1979.
  * [2] Cam Nguyen, _Analysis Methods for RF, Microwave, and Millimeter-Wave Planar Transmission Line Structures_ , Texas A and M University, John Wiley and sons, Inc., pp. 68-80, Aug. 2000.
  * [3] Bablu Kumar Singh, “Design of Rectangular Microstrip Patch Antenna based on Artificial Neural Network Algorithm,” Jodhpur Institute of Engineering and Technology, 2nd International Conference on SPIN, 2015.
  * [4] S. Gevorgian, H. Berg, H. Jacobsson, and T. Lewin, “Application Notes Of Basic Parameters Of Coplanar-Strip Waveguides On Multilayer Dielectric/Semiconductor Substrates,” _IEEE Microwave Magazine_ , vol. 4, issue 2, Jun. 2003.
  * [5] M. Naser-Moghaddasi, Pouya Derakhshan Barjoei, Alireza Naghsh, “A Heuristic Artificial Neural Network for Analyzing and synthesizing Rectangular Microstrip Antenna,” _IJCSNS International Journal of Computer Science and Network Security_ , vol.7, no.12, Dec. 2007.