\section{Introduction}

Distributed parameter systems described by partial differential equations
arise in systems theory from a modeling and a control theoretic point
of view and are without doubt a challenging research problem, where
lot of progress has been achieved in the last years. Also the port-Hamiltonian
setting, originally developed in the finite dimensional scenario has
been transfered to infinite-dimensional systems, where e.g. the well-known
approach based on (Stokes-)Dirac structures (also known from the lumped
parameter scenario) is available, see e.g. \cite{Maschke2005,SchaftMaschke,Machelli2004I,Machelli2004II,MachelliSIAM,Machelli2005CDC}
and references therein.

Also in mathematical physics, systems described by partial differential
equations (pdes) are interpreted in a Hamiltonian setting, e.g. in
\cite{Giachetta,Gotay,Olver} and references therein, but in most
cases systems with trivial boundary conditions are considered, which
is not the case in many engineerings applications. Therefore, the
approach based on Stokes-Dirac structures has been setup to overcome
the problem of non-zero energy flow through the boundary.

A different port-Hamiltonian approach is based on a bundle structure
with respect to independent and dependent coordinates (not necessarily
relying on an underlying Stokes-Dirac structure), see also \cite{Enns,SchoeberlMCMDSPiezo2008,SchlacherHam2008,SchoeberlMCMDS2010,SiukaActa,SchoeberlCDC,SchoeberlLHMNL}
which are all based on \cite{Olver} but adapted to control purposes,
i.e. modified in a sense, such that non-zero energy flow through the
boundary can be considered (such that boundary ports are included)
and furthermore control inputs on the domain and/or the boundary can
be included.

The main difference of the approach relying on Stokes-Dirac structures
and the approach using bundles, is that the Stokes-Dirac scenario
is based on the choice of proper energy variables (flows and efforts)
for which the power balance is formulated, whereas the second approach
is based on a given Hamiltonian density (the total energy density)
and the evaluation of the power balance is performed based on the
underlying bundle formalism in order to restructure the pdes such
that the energy flows are linked to the physics. This will have the
consequence, that the variational derivative is interpreted differently
and the choice of state variables is different, in the two mentioned
approaches.

The purpose of this paper is, that based on the well-known Lagrangian
setting for first order field theories i) the partial differential
equations and the boundary conditions derived using a variational
principle are reinterpreted using two different port-Hamiltonian settings,
which describe the same physical phenomenon but using a completely
different port-Hamiltonian representation, ii) by using the example
of a Mindlin plate all these concepts are visualized and compared
in great detail.

\section{Notation }

We will use differential geometric methods for our considerations
and the notation is similar to the one in \cite{Giachetta}, where
the interested reader can find much more details about this geometric
machinery. To keep the formulas short and readable we will use tensor
notation and especially Einstein's convention on sums.

We use the standard symbol $\wedge$ for the exterior product (wedge
product), $\mathrm{d}$ is the exterior derivative, $\rfloor$ the
natural contraction between tensor fields. By $\partial_{\alpha}^{B}$
are meant the partial derivatives with respect to coordinates with
the indices $_{B}^{\alpha}$ and $[m^{\alpha\beta}]$ corresponds
to the matrix representation of the (second-order) tensor $m$ with
components $m^{\alpha\beta}$. E.g. taking a second-order tensor $m$
and a co-vector $\omega$, the components of the contraction $m\rfloor\omega$
read in local coordinates as $m^{\alpha\beta}\omega_{\alpha}$, where
the summation over $\alpha$ is performed (Einstein convention on
sums).

Furthermore $C^{\infty}(\cdot)$ denotes the set of the smooth functions
on the corresponding manifold. Moreover we will not indicate the range
of the used indices when they are clear from the context. Additionally,
pull backs and pull back bundles are only stated when necessary, when
they follow from the context they are not indicated to avoid exaggerated
notation.

Let us consider the bundle $\mathcal{Y}\rightarrow\mathcal{D},\,(X^{A},y^{\alpha})\rightarrow(X^{A})$.
The first jet manifold $\mathcal{J}^{1}(\mathcal{Y})$ possesses the
coordinates $(X^{A},y^{\alpha},y_{A}^{\alpha})$, where the capital
Latin indices $A,B$ are used for the base manifold $\mathcal{D}$
(independent coordinates) and $y_{A}^{\alpha}$ denote derivative
coordinates of first order (derivatives of the dependent coordinates
with respect to the independent ones) as well as
\[
\partial_{A}=\frac{\partial}{\partial X^{A}}\,,\,\partial_{\alpha}=\frac{\partial}{\partial y{}^{\alpha}}\,,\,\partial_{\alpha}^{A}=\frac{\partial}{\partial y{}_{A}^{\alpha}}.
\]
The jet structure also induces the so-called total derivative
\[
d_{A}=\partial_{A}+y_{A}^{\alpha}\partial_{\alpha}+y_{AB}^{\alpha}\partial_{\alpha}^{B}
\]
 acting on elements including first order derivatives and $y_{AB}^{\alpha}$
correspond to derivative coordinates of second order living in $\mathcal{J}^{2}(\mathcal{Y})$,
the second jet manifold. Based on the bundle structure $\mathcal{Y}\rightarrow\mathcal{D}$
let us introduce the vertical tangent bundle $\mathcal{V}(\mathcal{Y})$,
as well as
\[
\Lambda_{1}^{d}(\mathcal{Y})=\mathcal{T}^{*}(\mathcal{Y})\wedge(\overset{d}{\wedge}\mathcal{T}^{*}(\mathcal{D})),
\]
see also \cite{Giachetta}, with a typical element $\omega=\omega_{\alpha}\mathrm{d}y^{\alpha}\wedge\mathrm{d}V$
for $\Lambda_{1}^{d}(\mathcal{Y})$ where $\mathrm{d}V$ denotes the
volume element on the manifold $\mathcal{D}$, i.e. $\mathrm{d}V=\mathrm{d}X^{1}\wedge\ldots\mathrm{\wedge d}X^{d}$
with $\mathrm{dim}(\mathcal{D})=d$ and the functions $\omega_{\alpha}$
may depend on derivative coordinates. Furthermore, a typical element
for $\mathcal{V}(\mathcal{Y})$ reads as $v=v^{\alpha}\partial_{\alpha}$
and when $v^{\alpha}$ depends on derivative coordinates we call $v$
a generalized vertical vector field, see \cite{Olver}.

\section{Background material}

\subsection{Geometric preliminaries}

We will consider densities $\mathfrak{F}$ in the sequel (a quantity
that can be integrated), $\mathfrak{F}=\mathcal{F}\mathrm{d}V$ with
$\mathcal{F}\in C^{\infty}(\mathcal{J}^{1}(\mathcal{X}))$ (we restrict
ourselves to the first-order case). By $F=\int_{\mathcal{D}}\mathfrak{F}$
we denote the integrated quantity, where of course a section of the
bundle $\mathcal{Y}\rightarrow\mathcal{D}$, i.e. a map $y=\Phi(X)$
leading to $y_{A}=\partial_{A}\Phi(X)$ has to be plugged in to be
able to evaluate the integral properly.
\begin{prop}
\label{prop:Given-the-density}Given the density $\mathfrak{F}=\mathcal{F}\mathrm{d}V$
and a generalized vertical vector field $v:\mathcal{D}\rightarrow\mathcal{V}(\mathcal{Y})$,
together with its first jet-prolongation $j^{1}(v)=v^{\alpha}\partial_{\alpha}+d_{A}(v^{\alpha})\partial_{\alpha}^{A}$,
see \cite{Olver,Giachetta}, we obtain the decomposition
\begin{eqnarray}
\int_{\mathcal{D}}j^{1}(v)(\mathcal{F}\mathrm{d}V) & = & \int_{\mathcal{D}}v\rfloor\delta\mathfrak{F}+\int_{\partial\mathcal{D}}v\rfloor\delta^{\partial}\mathfrak{F}.\label{eq:LieDiffPde-1}
\end{eqnarray}
Here the map $\delta\mathfrak{F}=\delta{}_{\alpha}\mathcal{F}\,\mathrm{d}y^{\alpha}\wedge\mathrm{d}V$
(corresponding to the Euler Lagrange operator of $\mathfrak{F})$,
see \cite{Giachetta}, with the coefficients $\delta_{\alpha}\mathcal{\mathcal{F}}=\partial_{\alpha}\mathcal{F}-d_{A}\partial_{\alpha}^{A}\mathcal{F}$
(called the variational derivatives) is used, as well as the boundary
operator $\delta^{\partial}\mathfrak{F}=\partial_{\alpha}^{A}\mathcal{F}\,\mathrm{d}y^{\alpha}\wedge\mathrm{d}V_{A}$
with $\mathrm{d}V_{A}=\partial_{A}\rfloor\mathrm{d}V$ (the boundary
volume form).
\end{prop}

\begin{proof}
The proof follows by evaluating the Lie-derivative of the geometric
object $\mathfrak{F}$ with respect to the vector field $j^{1}(v)$
\begin{eqnarray}
j^{1}(v)(\mathcal{F}\mathrm{d}V) & = & \left(v^{\alpha}(\partial_{\alpha}\mathcal{F}-d_{A}\partial_{\alpha}^{A}\mathcal{F})+d_{A}(v^{\alpha}\partial_{\alpha}^{A}\mathcal{F})\right)\mathrm{d}V\nonumber \\
 & = & \left(v^{\alpha}\delta_{\alpha}\mathcal{F}+d_{A}(v^{\alpha}\partial_{\alpha}^{A}\mathcal{F})\right)\mathrm{d}V\label{eq:intbyparts}
\end{eqnarray}
and applying the Theorem of Stokes \cite{Olver} to (\ref{eq:intbyparts}).
\end{proof}
The relation (\ref{eq:LieDiffPde-1}) will be of key interest in the
forthcoming, since it provides a natural decomposition of the expression
$\int_{\mathcal{D}}j^{1}(v)(\mathcal{F}\mathrm{d}V)$ into a term
on the domain $\mathcal{D}$ and one on the boundary $\partial\mathcal{D}$.
Important is the case when the generalized vector-field $v$ is linked
to the solution of a pde system (via its semi-group, that $v$ may
generate), then the formal change of $F=\int_{\mathcal{D}}\mathcal{F}\mathrm{d}V$
along solutions of a pde system can be computed as $\int_{\mathcal{D}}j^{1}(v)(\mathcal{F}\mathrm{d}V)$
(provided all operations are admissible), which we denote by $\dot{F}=\int_{\mathcal{D}}j^{1}(v)(\mathcal{F}\mathrm{d}V)$
in this special case.

\subsection{Dirac structures}

Based on the space of power variables $\mathcal{F}\times\mathcal{E}$
(flows and efforts) and the symmetric bilinear pairing
\begin{equation}
\ll(f_{1},e_{1}),(f_{2},e_{2})\gg:=\left\langle e_{1},f_{2}\right\rangle +\left\langle e_{2},f_{1}\right\rangle \label{eq:DiracPairing}
\end{equation}
where $\left\langle \cdot,\cdot\right\rangle $ is the dual product
of the linear spaces $\mathcal{F}$ and $\mathcal{E}=\mathcal{F}^{*}$
a Dirac structure is a linear subspace $\mathbb{D}\subset\mathcal{F}\times\mathcal{E}$
such that $\mathbb{D}=\mathbb{D}^{\bot}$ with respect to the pairing
(\ref{eq:DiracPairing}). For $(f,e)\in\mathbb{D}$ one has $\left\langle e,f\right\rangle =0$
such that the Dirac structure preserves power. This concept can be
transfered to the case where $\mathcal{F}$ and $\mathcal{E}$ are
spaces of vector-valued functions over a spatial domain $\mathcal{D}$,
then infinite dimensional systems are the focus, and to allow for
non-zero energy flow through the boundary the so-called Stokes-Dirac
structure is introduced, see \cite{Maschke2005,SchaftMaschke,Machelli2004I,Machelli2005CDC}
and section \ref{sub:Approach-based-on}.

\section{Lagrangian Framework}

In this section we recapitulate the well-known Lagrangian framework
for first-order field theories, and we will derive the partial differential
equations as well as the boundary conditions in a geometric fashion.
Thus, we consider a bundle
\[
\mathcal{Q}\rightarrow\mathcal{D}_{L},\,(q^{\alpha},t^{0},X^{A})\rightarrow(t^{0},X^{A})
\]
where we use the shortcut $x^{i}=(t^{0},X^{A})$ such that the independent
variables are the time $t^{0}$ and the spatial ones $X^{A}\,,\, A=1,\ldots\dim(\mathcal{D}_{L})-1=n$
. A first order Lagrangian takes the form

\begin{equation}
\mathfrak{L}:\mathcal{J}^{1}(Q)\rightarrow\overset{n+1}{\wedge}\mathcal{T}^{*}(\mathcal{D}_{L})\label{eq:}
\end{equation}
$\mathfrak{L}=\mathcal{L}\omega$ with $\mathcal{L}\in C^{\infty}\left(\mathcal{J}^{1}\left(Q\right)\right)$
together with the volume element $\omega$ that meets
\[
\omega=\mathrm{d}t^{0}\wedge\mathrm{d}X^{1}\ldots\wedge\mathrm{d}X^{n}\,,\,\,\,\omega_{i}=\partial_{i}\rfloor\omega\,,\,\,\,\, i=0,\ldots,n
\]
The variational problem for a section $s:\mathcal{D}_{L}\mathcal{\rightarrow Q}$
is the following%
\footnote{At this point the pull back is essential and therefore indicated,
i.e. $j^{1}(\psi_{\epsilon}\circ s)^{*}\mathfrak{L}$ means, that
the first prolongation of $(\psi_{\epsilon}\circ s)$ has to be plugged
in into $\mathfrak{L}=\mathcal{L}\omega$ in order to evaluate the
integral properly.%
}
\begin{equation}
\left.\left(d_{\epsilon}\int_{\mathcal{D}_{L}}j^{1}(\psi_{\epsilon}\circ s)^{*}\mathfrak{L}\right)\right|_{\epsilon=0}=0,\label{eq:VarProb}
\end{equation}
where the flow $\psi_{\epsilon}$ is used to deform sections $s:\mathcal{D}_{L}\mathcal{\rightarrow Q}$
and whose generator is a vertical vector field $v_{L}:\mathcal{Q}\rightarrow\mathcal{V}(Q)$.

This is a well-known problem and treated for example in \cite{Giachetta,Olver}
and references therein. It is obvious that (\ref{eq:VarProb}) is
equivalent to
\begin{eqnarray}
\int_{\mathcal{D}_{L}}(j^{2}s)^{*}\left(j^{1}(v_{L})(\mathfrak{L})\right) & = & 0,\label{eq:VarProb1}
\end{eqnarray}
see for example \cite{Giachetta} and based on (\ref{eq:LieDiffPde-1})
where we replace $\mathcal{F}$ be $\mathcal{L}$ we obtain the decomposition
\begin{equation}
\int_{\mathcal{D}_{L}}v_{L}\rfloor\delta\mathfrak{L}+\int_{\partial\mathcal{D}_{L}}v_{L}\rfloor\delta^{\partial}\mathfrak{L}=0.\label{eq:VarSplitB}
\end{equation}
Consequently, the partial differential equations for a first order
Lagrangian follow as
\begin{equation}
\delta_{\alpha}(\mathcal{L})=0\,,\,\,\,\delta_{\alpha}=\partial_{\alpha}-d_{i}\partial_{\alpha}^{i}\label{eq:LAGEUL}
\end{equation}
and the boundary term is the second term in (\ref{eq:VarSplitB})
and reads in local coordinates as
\begin{equation}
\int_{\partial\mathcal{D}_{L}}v_{L}^{\alpha}\partial_{\alpha}^{i}\mathcal{L}\omega_{i}=0.\label{eq:BondGen}
\end{equation}
The boundary conditions can be fulfilled by either allowing for no
variations on (a part of) $\partial\mathcal{D}_{L}$, i.e. $v_{L}^{\alpha}=0$
or by $\partial_{\alpha}^{i}\mathcal{L}=0$ or by a combinations of
both approaches. Possible is also the inclusion of external boundary
variables $F_{e,\alpha}^{i}$ such that
\begin{equation}
\int_{\partial\mathcal{D}_{L}}v_{L}^{\alpha}(\partial_{\alpha}^{i}\mathcal{L}-F_{e,\alpha}^{i})\omega_{i}=0.\label{eq:BondGen-1}
\end{equation}
has to be met.

\section{Port-Hamiltonian Picture}

Now we turn to the Hamiltonian picture, where we discuss two different
port-Hamiltonian formulations. We will restrict ourselves to systems
without dissipation and without distributed control for simplicity,
but these properties can be included in both formalisms in a straightforward
manner, see \cite{SchaftMaschke,SchoeberlMCMDSPiezo2008}.

\subsection{Geometric approach based on underlying bundle structure\label{sub:Geometric-approach-based}}

We will introduce port-Hamiltonian systems described by pdes based
on a power balance relation, such that the power balance relation
together with the structure of the equations represent the physical
process.
\begin{defn}
\label{A-port-controlled}A port-Hamiltonian boundary control system
without dissipation on a bundle $\mathcal{X}\rightarrow\mathcal{D_{H}},\,(x^{\alpha},X^{A})\rightarrow(X^{A})$
takes the form of
\begin{equation}
\begin{array}{ccl}
\dot{x} & = & \mathcal{J}(\mathrm{\delta}\mathfrak{H)}\end{array}\label{eq:HamTimeInvStandardCF-1}
\end{equation}
with the Hamiltonian $\mathfrak{H}=\mathcal{H}\Omega\,,\,\,\Omega=\mathrm{d}X^{1}\wedge\ldots\wedge\mathrm{d}X^{d}$,
where $\mathcal{H}\in C^{\infty}(\mathcal{J}^{1}(\mathcal{X}))$ and
additional boundary conditions (possibly including boundary inputs,
optionally leading to so-called boundary ports). The map $\mathcal{J}$
is of the form $\mathbb{\mathcal{J}}:\Lambda_{1}^{d}(\mathcal{X})\rightarrow\mathcal{V\textrm{(}\mathcal{X}\textrm{)}}$
where $\mathcal{J}$ is a skew-symmetric map.
\end{defn}
In general the map $\mathcal{J}$ can be a differential operator,
see our paper \cite{SchoeberlLHMNL}, but within this contribution
we exclude this case (since in many examples, e.g. mechanics this
is not required).

Now we make use of proposition \ref{prop:Given-the-density} and replace
$\mathfrak{F}$ by $\mathfrak{H}$ in (\ref{eq:LieDiffPde-1}). Setting
$v=\dot{x}$ we obtain
\begin{eqnarray}
\dot{H} & = & \int_{\partial\mathcal{D_{H}}}\dot{x}\rfloor\delta^{\partial}\mathfrak{H}=\int_{\partial\mathcal{D_{H}}}\dot{x}^{\alpha}\partial_{\alpha}^{A}\mathcal{H}\Omega_{A}\label{eq:HdotPdeoO}
\end{eqnarray}
where $\Omega_{A}=\partial_{A}\rfloor\Omega$, which reflects the
power balance, since the total change of the functional $H$ along
solutions of (\ref{eq:HamTimeInvStandardCF-1}), is affected by a
boundary port (if it exists) depending on the boundary conditions.
See e.g. \cite{SchoeberlMCMDSPiezo2008} for a formal introduction
concerning the boundary ports.

\subsection{Approach based on underlying Stokes-Dirac Structure\label{sub:Approach-based-on}}

Following \cite{Machelli2005CDC} we shortly recapitulate the port-Hamiltonian
framework based on Stokes-Dirac structures. For more details we refer
to \cite{Maschke2005,SchaftMaschke}.

\begin{flushleft}
We consider the space of flows $\mathcal{F}$ and the space of efforts
$\mathcal{E}$, which are spaces of vector-valued functions over a
spatial domain $\mathcal{D}$. Given $\mathcal{J}_{SD}$ a skew-adjoint
matrix differential operator \cite{Machelli2005CDC}, the space
\[
\mathbb{D}=\left\{ \left.\left(f,e,w\right)\in\mathcal{F}\times\mathcal{E}\times\mathcal{W}\right|f=-\mathcal{J}_{SD}e,w=\mathcal{B}_{\mathcal{D}}(e)\right\}
\]
is a Stokes-Dirac structure, regarding the pairing
\begin{align*}
\begin{array}{c}
\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\ll\left(f_{1},e_{1},w_{1}\right),\left(f_{2},e_{2},w_{2}\right)\gg=\\
\qquad=\int_{\mathcal{D}}[e_{1}^{T}f_{2}+e_{2}^{T}f_{1}]\mathrm{d}V+\int_{\partial\mathcal{D}}\mathcal{B}_{\mathcal{J}_{SD}}(w_{1},w_{2})\mathrm{d}A
\end{array}
\end{align*}
where $\mathcal{B}_{\mathcal{J}_{SD}}$ is a boundary differential
operator induced by $\mathcal{J}_{SD}$ any by slight abuse of notation
$\mathrm{d}A$ corresponds to the boundary volume element. The map
$\mathcal{B}_{\mathcal{D}}$ is a boundary operator and the boundary
variables are $w$.
\par\end{flushleft}

If $(f,e,w)\in\mathbb{D}$ then
\begin{eqnarray}
0 & = & \int_{\mathcal{D}}e^{T}f\mathrm{d}V+\frac{1}{2}\int_{\partial\mathcal{D}}\mathcal{B}_{\mathcal{J}_{SD}}(w,w)\mathrm{d}A\label{eq:DiracEnergy}
\end{eqnarray}
holds. Given an energy density $\mathcal{H}\mathrm{d}V$ where $\mathcal{H}$
depends on the energy variables, a port-Hamiltonian boundary control
system without dissipation can be stated as
\begin{equation}
f=-\mathcal{J}_{SD}e\,,\,\, w=\mathcal{B}_{\mathcal{D}}(e)\label{eq:pdeDirac}
\end{equation}
where the energy variables $f$, are linked to the state variables
$\chi$ via $f=-\dot{\chi}$ and the efforts variables follow from
$e=\partial_{\chi}\mathcal{H}$.
\begin{rem}
Originally, in \cite{Machelli2005CDC} instead of $e=\partial_{\chi}\mathcal{H}$
the authors use $e=\delta_{\chi}H$ where $H=\int_{\mathcal{D}}\mathcal{H}\mathrm{d}V$,
but since $\mathcal{H}$ depends on energy variables, the variational
derivative degenerates to a 'partial' one.
\end{rem}
Furthermore, from $H=\int_{\mathcal{D}}\mathcal{H}\mathrm{d}V$ and
the relations (\ref{eq:DiracEnergy}) and (\ref{eq:pdeDirac}) one
has
\begin{equation}
\dot{H}=\frac{1}{2}\int_{\partial\mathcal{D}}\mathcal{B}_{\mathcal{J}_{SD}}(w,w)\mathrm{d}A\label{eq:BalanceDirac}
\end{equation}
since $\dot{H}=\int_{\mathcal{D}}\partial_{\chi}\mathcal{H}\dot{\xi}\mathrm{d}V=-\int_{\mathcal{D}}e^{T}f\mathrm{d}V$.
\begin{rem}
To derive this energy balance also proposition \ref{prop:Given-the-density}
can be applied, but since no jet-variables are included it simplifies
to $\dot{H}=\int_{\mathcal{D}}j^{1}(v)(\mathcal{H}\mathrm{d}V)=\int_{\mathcal{D}}v^{\alpha}\partial_{\alpha}\mathcal{H}\mathrm{d}V$
and to derive (\ref{eq:BalanceDirac}) a further integration by parts
must be performed, since $v$ corresponds to $-f$ which involves
the differential operator $\mathcal{J}_{SD}$.
\end{rem}

\section{The Mindlin Plate}

Let us consider a rectangular plate with lengths $l_{x},l_{y}$, where
$h$ will denote the thickness, which will be modeled based on the
hypothesis stated by Mindlin. Therefore, we choose as independent
coordinates the vertical deflection $w$ of the mid-plane as well
as the rotations of a transverse normal to the $X$ and $Y$ direction
termed $\psi$ and $\phi,$ respectively. The kinetic energy density
$\mathcal{K}$ and the potential energy density $\mathcal{V}$ can
be stated as

\begin{eqnarray*}
\mathcal{K} & = & \frac{\rho}{2}(\frac{h^{3}}{12}(\psi_{t}^{2}+\phi_{t}^{2})+hw_{t}^{2})
\end{eqnarray*}
and

\begin{eqnarray*}
\mathcal{V} & = & \frac{1}{2}kGh\left[(w_{X}-\psi)^{2}+(w_{Y}-\phi)^{2}\right]\\
 &  & +\frac{1}{2}D\frac{1-\nu}{2}(\psi_{Y}+\phi_{X})^{2}\\
 &  & +\frac{1}{2}(D(\psi_{X}^{2}+\nu\phi_{Y}\psi_{X})+D(\phi_{Y}^{2}+\nu\phi_{Y}\psi_{X})),
\end{eqnarray*}
where $\nu$ is the Poisson ratio, $k=\frac{\pi^{2}}{12},$ and $G,D$
are the plate stiffness and the plate module, respectively, see \cite{Machelli2005CDC}
and references therein.
\begin{rem}
The subscripts $t,X,Y$ correspond to the derivatives with respect
to these independent variables, according to the jet-bundle structure
in the Lagrangian framework. The subscripts $x,y$ to be used later,
correspond to quantities which are connected to the spatial variables
$X$ and $Y$ but they must not be confused with derivative variables.
Since we are in a time-invariant setting, we will use later on also
the $\dot{}$ notation, for time derivatives, instead of the subscripts
$t$.
\end{rem}
To derive the equations of motion we will use the variational principle
in a Lagrangian setting. Then given the partial differential equations,
we will interpret them in a Hamiltonian setting, either using the
approach presented in section \ref{sub:Geometric-approach-based}
and using an approach based on the Stokes-Dirac structure as in section
\ref{sub:Approach-based-on}.

\subsection{The Lagrangian picture}

In the Lagrangian framework we consider the bundle
\begin{equation}
\mathcal{Q}\rightarrow\mathcal{D}_{\mathcal{L}},\,(w,\psi,\phi,t,X,Y)\rightarrow(t,X,Y)\label{eq:BundleLag}
\end{equation}
together with the Lagrangian density $\mathfrak{L}=\mathcal{L}\omega$
with $\mathcal{L}=\mathcal{K}-\mathcal{V},\,\omega=\mathrm{d}t\wedge\mathrm{d}X\wedge\mathrm{d}Y.$
The variational derivatives follow form the chosen bundle structure
(\ref{eq:BundleLag}) and follow to
\begin{eqnarray*}
\delta_{w} & = & \partial_{w}-d_{t}\partial_{w}^{t}-d_{X}\partial_{w}^{X}-d_{Y}\partial_{w}^{Y}\\
\delta_{\psi} & = & \partial_{\psi}-d_{t}\partial_{\psi}^{t}-d_{X}\partial_{\psi}^{X}-d_{Y}\partial_{\psi}^{Y}\\
\delta_{\phi} & = & \partial_{\phi}-d_{t}\partial_{\phi}^{t}-d_{X}\partial_{\phi}^{X}-d_{Y}\partial_{\phi}^{Y}.
\end{eqnarray*}
From $\delta_{w}\mathcal{L}=0,\,\delta_{\psi}\mathcal{L}=0,\,\delta_{\phi}\mathcal{L}=0$
corresponding to (\ref{eq:LAGEUL}) we derive the partial differential
equations
\begin{eqnarray*}
\rho hw_{tt} & = & kGh(w_{XX}-\psi_{X})+kGh(w_{YY}-\phi_{Y})\\
\rho\frac{h^{3}}{12}\psi_{tt} & = & kGh(w_{X}-\psi)+D(\psi_{XX}+\nu\phi_{XY})\\
 &  & +\frac{1}{2}D(1-\nu)(\psi_{YY}+\phi_{XY})\\
\rho\frac{h^{3}}{12}\phi_{tt} & = & kGh(w_{Y}-\phi)+D(\phi_{YY}+\nu\psi_{XY})\\
 &  & +\frac{1}{2}D(1-\nu)(\psi_{XY}+\phi_{XX}).
\end{eqnarray*}
If we introduce
\begin{eqnarray}
M_{x} & = & D(\psi_{X}+\nu\phi_{Y})=-\partial_{\psi}^{X}\mathcal{L}\nonumber \\
M_{y} & = & D(\phi_{Y}+\nu\psi_{X})=-\partial_{\phi}^{Y}\mathcal{L}\nonumber \\
M_{xy} & = & D\frac{1-\nu}{2}(\psi_{Y}+\phi_{X})=-\partial_{\psi}^{Y}\mathcal{L=}-\partial_{\phi}^{X}\mathcal{L}\label{eq:MQ}\\
Q_{x} & = & kGh(w_{X}-\psi)=-\partial_{w}^{X}\mathcal{L}\nonumber \\
Q_{y} & = & kGh(w_{Y}-\phi)=-\partial_{w}^{Y}\mathcal{L}\nonumber
\end{eqnarray}
then the equations of motion take the familiar form%
\footnote{In order to be comparable with the literature we sometimes use $\partial_{X}$and
$\partial_{Y}$ although in a strict mathematical sense it should
be $d_{X}$ and $d_{Y}.$%
}
\begin{eqnarray}
\rho hw_{tt} & = & \partial_{X}Q_{x}+\partial_{Y}Q_{y}\nonumber \\
\rho\frac{h^{3}}{12}\psi_{tt} & = & Q_{x}+\partial_{X}M_{x}+\partial_{Y}M_{xy}\label{eq:PDELagMind}\\
\rho\frac{h^{3}}{12}\phi_{tt} & = & Q_{y}+\partial_{Y}M_{y}+\partial_{X}M_{xy}.\nonumber
\end{eqnarray}
The boundary conditions follow from
\[
\int_{\partial\mathcal{D_{L}}}\left(w_{L}\partial_{w}^{i}\mathcal{L}+\psi_{L}\partial_{\psi}^{i}\mathcal{L}+\phi_{L}\partial_{\phi}^{i}\mathcal{L}\right)\partial_{i}\rfloor(\mathrm{d}t\wedge\mathrm{d}X\wedge\mathrm{d}Y)=0
\]
where the variational vector field $v_{L}$ takes the form $v_{L}=w_{L}\partial_{w}+\psi_{L}\partial_{\psi}+\phi_{L}\partial_{\phi}$
and on the time-boundary no variation takes place (i.e. when $i=0$
then $v_{L}=0$) and (\ref{eq:MQ}) has to be used. Therefore we have
\begin{eqnarray}
\int_{\partial\mathcal{D_{L}}}\left(w_{L}Q_{x}+\psi_{L}M_{x}+\phi_{L}M_{xy}\right)\mathrm{d}t\wedge\mathrm{d}Y\nonumber \\
-\int_{\partial\mathcal{D_{L}}}\left(w_{L}Q_{y}+\psi_{L}M_{xy}+\phi_{L}M_{y}\right)\mathrm{d}t\wedge\mathrm{d}X=0\label{eq:LAGBondMind}
\end{eqnarray}
such that, e.g. if at $X=0$ we have that $w_{L}$ is arbitrary, then
$Q_{x}$ has to vanish or has to be compensated by an external boundary
term as in (\ref{eq:BondGen-1}), such that the familiar boundary
conditions are recovered.

\subsection{The Hamiltonian picture}

Based on the partial differential equations (\ref{eq:PDELagMind})
and the boundary conditions (\ref{eq:LAGBondMind}) we discuss the
two presented port-Hamiltonian formulations as well as the power balance
relations corresponding to the particular representation.

\subsubsection{Geometric approach}

Now we consider the bundle (which is different form the one in the
Lagrangian setting)
\[
\mathcal{X}\rightarrow\mathcal{D_{H}},\,(w,\psi,\phi,p_{w},p_{\psi},p_{\phi},X,Y)\rightarrow(X,Y).
\]
From the Legendre transform we derive the temporal momenta $p_{w}=\partial_{w}^{t}\mathcal{L},\, p_{\psi}=\partial_{\psi}^{t}\mathcal{L},\, p_{\phi}=\partial_{\phi}^{t}\mathcal{L}$
which read as
\begin{equation}
p_{w}=\rho hw_{t}\,,\,\, p_{\psi}=\rho\frac{h^{3}}{12}\psi_{t}\,,\,\, p_{\phi}=\rho\frac{h^{3}}{12}\phi_{t}\label{eq:Legendre}
\end{equation}
and the Hamiltonian follows as
\[
\mathcal{H}=\dot{w}p_{w}+\dot{\psi}p_{\psi}+\dot{\phi}p_{\phi}-\mathcal{L}.
\]
In the coordinates $(w,\psi,\phi,p_{w},p_{\psi},p_{\phi})$ together
with (\ref{eq:Legendre}) one has $\mathcal{H}=\mathcal{K}+\mathcal{V}$.

To derive the Hamiltonian formulation as in (\ref{eq:HamTimeInvStandardCF-1})
we set $x=(w,\psi,\phi,p_{w},p_{\psi},p_{\phi})$ and obtain
\begin{eqnarray}
\dot{x} & = & \mathcal{J}(\delta\mathfrak{H})\label{eq:pdeGeoEx}
\end{eqnarray}
which reads as
\[
\left[\begin{array}{c}
\dot{w}\\
\dot{\psi}\\
\dot{\phi}\\
\dot{p}_{w}\\
\dot{p}_{\psi}\\
\dot{p}_{\phi}
\end{array}\right]=\left[\begin{array}{cccccc}
0 & 0 & 0 & 1 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 1\\
-1 & 0 & 0 & 0 & 0 & 0\\
0 & -1 & 0 & 0 & 0 & 0\\
0 & 0 & -1 & 0 & 0 & 0
\end{array}\right]\left[\begin{array}{c}
\delta_{w}\mathcal{H}\\
\delta_{\psi}\mathcal{H}\\
\delta_{\phi}\mathcal{H}\\
\delta_{p_{w}}\mathcal{H}\\
\delta_{p_{\psi}}\mathcal{H}\\
\delta_{p_{\phi}}\mathcal{H}
\end{array}\right].
\]
The variational derivatives in this setting take the form
\begin{eqnarray*}
\delta_{w} & = & \partial_{w}-d_{X}\partial_{w}^{X}-d_{Y}\partial_{w}^{Y}\\
\delta_{\psi} & = & \partial_{\psi}-d_{X}\partial_{\psi}^{X}-d_{Y}\partial_{\psi}^{Y}\\
\delta_{\phi} & = & \partial_{\phi}-d_{X}\partial_{\phi}^{X}-d_{Y}\partial_{\phi}^{Y}
\end{eqnarray*}
and $\delta_{p_{w}}=\partial_{p_{w}}\,,\,\delta_{p_{\psi}}=\partial_{p_{\psi}}\,,\,\delta_{p_{\phi}}=\partial_{p_{\phi}}$because
of the different bundle structure compared to the Lagrangian approach.

The boundary ports follow from (\ref{eq:HdotPdeoO}) and we obtain
\begin{eqnarray*}
\dot{H} & = & \int_{\partial\mathcal{D}_{H}}\dot{x}^{\alpha}\partial_{\alpha}^{A}\mathcal{H}\partial_{A}\rfloor(\mathrm{d}X\wedge\mathrm{d}Y)\\
 & = & \int_{\partial\mathcal{D_{H}}}(\dot{x}^{\alpha}\partial_{\alpha}^{X}\mathcal{H}\mathrm{d}Y-\dot{x}^{\alpha}\partial_{\alpha}^{Y}\mathcal{H}\mathrm{d}X).
\end{eqnarray*}
From the special choice of the Hamiltonian $\mathcal{H}$ we observe
that the expressions $\partial_{\alpha}^{A}\mathcal{L}$ and $\partial_{\alpha}^{A}\mathcal{H}$
correspond (apart form the sign), and therefore we also have
\begin{eqnarray*}
\partial_{w}^{X}\mathcal{H}=Q_{x} & \,,\, & \partial_{w}^{Y}\mathcal{H}=Q_{y}\\
\partial_{\psi}^{X}\mathcal{H}=M_{x} & \,,\, & \partial_{\psi}^{Y}\mathcal{H}=M_{xy}\\
\partial_{\phi}^{X}\mathcal{H}=M_{xy} & \,,\, & \partial_{\phi}^{Y}\mathcal{H}=M_{y}
\end{eqnarray*}
which consequently leads to the power balance relation
\begin{eqnarray}
\dot{H} & = & \int_{\partial\mathcal{D_{H}}}(\dot{w}Q_{x}+\dot{\psi}M_{x}+\dot{\phi}M_{xy})\mathrm{d}Y\nonumber \\
 &  & -\int_{\partial\mathcal{D_{H}}}(\dot{w}Q_{y}+\dot{\psi}M_{yx}+\dot{\phi}M_{y})\mathrm{d}X\label{eq:HdotMinGeo}
\end{eqnarray}
which is based on (\ref{eq:HdotPdeoO}).
\begin{rem}
It should be noted that depending on the boundary conditions, along
$\partial\mathcal{D}$ a boundary port appears only if in the pairings
$\dot{x}^{\alpha}\partial_{\alpha}^{A}\mathcal{H}$ both 'players'
are not equal to zero. Furthermore it should be noted that in $\dot{x}^{\alpha}$
only $(\dot{w},\dot{\psi},\dot{\phi})$ remain, since in $\mathcal{H}$
only jet variables with respect to $w,\psi,\phi$ appear, i.e. there
is a $w_{X}$ present but no $(p_{w})_{X}$ and so on.
\end{rem}

\subsubsection{The Stokes-Dirac approach}

This approach is not based on a bundle structure, which distinguishes
dependent and independent variables strictly, but uses so-called energy
variables instead. Therefore, let us introduce the strain variables
as \cite{Machelli2005CDC}
\begin{eqnarray}
\Gamma_{x} & = & -\psi_{X}\nonumber \\
\Gamma_{y} & = & -\phi_{Y}\nonumber \\
\Gamma_{xy} & = & -(\psi_{Y}+\phi_{X})\label{eq:Strains}\\
\Gamma_{xz} & = & w_{X}-\psi\nonumber \\
\Gamma_{yz} & = & w_{Y}-\phi.\nonumber
\end{eqnarray}
Then one can introduce as state $\chi$ which consists of the momentum
variables and the strains
\[
\chi=(\rho h\dot{w},\Gamma_{xz},\Gamma_{yz},\rho\frac{h^{3}}{12}\dot{\psi},\rho\frac{h^{3}}{12}\dot{\phi},\Gamma_{x},\Gamma_{y},\Gamma_{xy})
\]
and the Hamiltonian $\mathcal{H}=\mathcal{K}+\mathcal{V}$ can be
rewritten as
\begin{eqnarray*}
\mathcal{K} & = & \frac{\rho}{2}(\frac{h^{3}}{12}(\dot{\psi}^{2}+\dot{\phi}^{2})+h\dot{w}^{2})
\end{eqnarray*}
and
\begin{eqnarray*}
\mathcal{V} & = & \frac{1}{2}(Q_{x}\Gamma_{xz}+Q_{y}\Gamma_{yz}-M_{xy}\Gamma_{xy}-M_{x}\Gamma_{x}-M_{y}\Gamma_{y})
\end{eqnarray*}
such that $e=\partial_{\chi}\mathcal{H}$ follows to
\[
e=(\dot{w},Q_{x},Q_{y},\dot{\psi},\dot{\phi},-M_{x},-M_{y},-M_{xy}).
\]
The partial differential equations can be stated as
\begin{equation}
\dot{\chi}=\mathcal{J}_{SD}e=\mathcal{J}_{SD}\partial_{\chi}\mathcal{H}\label{eq:pdeDiracEx}
\end{equation}
where $\mathcal{J}_{SD}$ takes the form
\[
\left[\begin{array}{cccccccc}
0 & \partial_{X} & \partial_{Y} & 0 & 0 & 0 & 0 & 0\\
\partial_{X} & 0 & 0 & -1 & 0 & 0 & 0 & 0\\
\partial_{Y} & 0 & 0 & 0 & -1 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0 & -\partial_{X} & 0 & -\partial_{Y}\\
0 & 0 & 1 & 0 & 0 & 0 & -\partial_{Y} & -\partial_{X}\\
0 & 0 & 0 & -\partial_{X} & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & -\partial_{Y} & 0 & 0 & 0\\
0 & 0 & 0 & -\partial_{Y} & -\partial_{X} & 0 & 0 & 0
\end{array}\right].
\]
The energy balance follows from equation (\ref{eq:BalanceDirac})
where (\ref{eq:DiracEnergy}) has to be evaluated by a further integration
by parts, since $\mathcal{J}_{SD}$ is a differential operator which
contributes to the boundary expression, i.e. $\mathcal{B}_{\mathcal{J}_{SD}}$
has to be constructed, see \cite{Machelli2005CDC}.

\section{Discussion and Comparison}

In this section we will discuss the main differences of the two presented
port-Hamiltonian scenarios as in sections (\ref{sub:Geometric-approach-based})
and (\ref{sub:Approach-based-on}) where we will highlight these aspects
by focusing on the presented example, the Mindlin plate.

\subsection{State variables}

The state variables in the geometric approach consist of the displacements/deflections
and the temporal momenta, i.e. $x=(w,\psi,\phi,p_{w},p_{\psi},p_{\phi})$
in our example where in many cases (mechanical systems) the temporal
momenta can be derived from a given Lagrangian by means of the Legendre
transformation. These temporal momenta are introduced mainly to obtain
explicit partial differential equations where the state variables
are differentiated with respect to a curve parameter, which is the
time.

In contrast to this, in the approach based on the Stokes-Dirac structures,
energy variables are used, such that the strains are introduced in
mechanical applications and one has $\chi=(\rho h\dot{w},\Gamma_{xz},\Gamma_{yz},\rho\frac{h^{3}}{12}\dot{\psi},\rho\frac{h^{3}}{12}\dot{\phi},\Gamma_{x},\Gamma_{y},\Gamma_{xy})$.
From a conceptional point of view the use of energy/power variables
may be beneficial since they are linked to the power balance relation
in a simple manner, but as can be seen already using the Mindlin plate
example, the five strain variables have to be derived from the three
independent deflection/displacement variables $(w,\psi,\phi)$ by
differentiation, see (\ref{eq:Strains}), such that additionally to
the partial differential equations (\ref{eq:pdeDiracEx}) also the
compatibility conditions (\ref{eq:Strains}) must be listed, such
that a constrained Hamiltonian representation is apparent.

\subsection{Control issues}

From a control point of view control methods like damping injection
or control by interconnection can be performed equivalently using
the two presented Hamiltonian representations, see for example \cite{SiukaActa,MachelliSIAM}
where the Timoshenko beam is analyzed. However as stated also above,
the use of energy variables allows for controlling the system for
instance to zero strain configuration, but the global position in
space cannot be controlled in a straightforward manner since the deflection/displacement
coordinates do not enter the formalism, in contrast to the approach
as in (\ref{sub:Geometric-approach-based}).

\subsection{The skew-symmetric operators $\mathcal{J}$ and $\mathcal{J}_{SD}$ }

By inspection it becomes apparent that $\mathcal{J}$ and $\mathcal{J}_{SD}$
differ significantly, since $\mathcal{J}$ is no differential operator
in contrast to $\mathcal{J}_{SD}$. This also has severe consequences
for the expressions $\delta\mathfrak{H}$ as in (\ref{eq:HamTimeInvStandardCF-1})
and $e=\partial_{\chi}\mathcal{H}$ as in (\ref{eq:pdeDirac}), where
it is vice versa, i.e. $\delta$ is a variational derivative and in
the Stokes-Dirac approach a partial derivative appears. Let us consider
for instance the fourth equation of (\ref{eq:pdeDiracEx}) which reads
as
\[
\rho\frac{h^{3}}{12}\ddot{\psi}=1\cdot Q_{x}-\partial_{X}(-M_{x})-\partial_{Y}(-M_{xy})
\]
as well as the fifth equation of (\ref{eq:pdeGeoEx}) which is
\[
\dot{p}_{\psi}=-\delta_{\psi}\mathcal{H}=-(\partial_{\psi}-d_{X}\partial_{\psi}^{X}-d_{Y}\partial_{\psi}^{Y})\mathcal{H}.
\]
From $\partial_{\psi}\mathcal{H}=-Q_{x}$ as well as from $\partial_{\psi}^{X}\mathcal{H}=M_{x}\,,\,\partial_{\psi}^{Y}\mathcal{H}=M_{xy}$
we easily observe that a part of the variational derivate $\delta_{\psi}=\partial_{\psi}-d_{X}\partial_{\psi}^{X}-d_{Y}\partial_{\psi}^{Y}$
is incorporated in $\mathcal{J}_{SD}$ (namely $(1,-\partial_{X},-\partial_{Y})$)
and that the variables $Q_{x},\, M_{x}$ and $M_{xy}$ are used directly.

\subsection{The energy balances}

From the relation (\ref{eq:HdotPdeoO}) the power balance is derived
easily once the state $x$ as well as the Hamiltonian density $\mathcal{H}$
is chosen - this is very simple, since $\mathcal{J}$ is no differential
operator and does not contribute to the boundary term. This is different
in the approach as in section (\ref{sub:Approach-based-on}) since
the boundary operator $\mathcal{B}_{\mathcal{J}_{SD}}$ has to be
derived from the special choice of $\mathcal{J}_{SD}$ and by an additional
integration by parts one ends up again by the same relation as in
(\ref{eq:HdotMinGeo}).

\section{Conclusion}

We have presented two different port-Hamiltonian representations based
on a given set of partial differential equations together with their
boundary conditions derived by the Lagrangian formalism using jet-bundles.
By means of the running example, the Mindlin plate, we have extensively
discussed and compared these two different Hamiltonian scenarios.
Further investigations should also include the field theoretic Hamiltonian
concepts coming from mathematical physics, like the polysymplectic
and/or the multisymplectic approach as in \cite{Giachetta,Gotay}
or in the spirit as in \cite{Nishida}.


\section{Introduction}
Conformal invariance is anticipated to
emerge in asymptotically free non-Abelian gauge theories
when the fermion species (flavor $N_f$)
exceeds a critical value $N_f={N_f}^*$
\cite{Caswell:1974gg,Banks:1981nn}.
The approach to the conformality from below
is in principle associated
with a pre-conformal (walking) behavior of the running
coupling, which has been advocated as a basis for strongly interacting
mechanisms of electroweak symmetry breaking~\cite{Sannino:2009za}.

Recent lattice studies~\cite{DelDebbio:2010zz}
focused on the computation of ${N_f}^*$
and the analysis of the conformal window itself,
either with fundamental fermions
or other representations.
Among the many interesting results with fundamental fermions, we single
out the observation
that the color SU(3) gauge theory with $N_f = 8$
is still in the hadronic phase
\cite{Appelquist:2009ty,Deuzeman:2008sc},
while $N_f = 12$ seems to be close to the critical number of flavors,
with some groups favoring conformality
\cite{Appelquist:2011dp,Appelquist:2009ty,
Deuzeman:2009mh,Hasenfratz:2010fi},
and others chiral symmetry breaking~\cite{Fodor:2011tu}.

In order to attack the walking and the conformal dynamics,
it is more informative beyond a fixed $N_f$
to investigate the vanishing or reducing chiral dynamics
with increasing varying $N_f$.
To this end,
we investigate the $N_f$ dependences of
the chiral phase transition at finite temperature ($T$)
based on our recent works~\cite{Deuzeman:2008sc,Miura:2011mc}:
The vanishing (reducing) finite T step-scaling
which attributes to the uniqueness of $T_c$ at each $N_f$ signals
the emergence of the (pre-)conformality.
Further, motivated by the recent
functional renormalization group (FRG) studies~\cite{BraunGies},
we examine the $N_f$ dependence of $T_c$,
whose vanishing (decreasing) behavior indicates
the (pre-)conformality.
This {\em thermodynamic} lattice study for
the large $N_f$ non-Abelian gauge theory
has played a crucial role to extract a notion of
more strongly interacting non-Abelian plasma~\cite{Liao:2012tw},
and it is expected to provide a new connection
between the lattice and the Gauge/Gravity duality
\cite{Gursoy:2010fj}.

\section{Simulation setups}
Simulations have been performed by utilizing
the publicly available MILC code~\cite{MILC}.
We use an improved version of the staggered action, the
Asqtad action, with a one-loop Symanzik
\cite{Symanzik}
and tadpole~\cite{LM1985} improved gauge action.
The tadpole factor $u_0$ is determined
by performing zero temperature simulations
on the $12^4$ lattice, and used as an input
for finite temperature simulations.
To generate configurations with mass degenerate
dynamical flavors,
we have used the rational hybrid Monte Carlo algorithm
(RHMC)~\cite{Clark:2006wq}.

Our observables are the chiral condensate
and the Polykov loop
\begin{equation}
a^3\langle\bar{\psi}\psi\rangle =
\frac{N_f}{4N_s^3N_t}
\Big\langle\mathrm{Tr\bigl[M^{-1}\bigr]}\Big\rangle
\ ,\quad
L =
\frac{1}{N_cN_s^3}\sum_{\mathbf{x}}
\mathrm{Re}
\bigg\langle
\mathrm{tr}_c\prod_{t=1}^{N_t}U_{4,t\mathbf{x}}
\bigg\rangle
\ ,\label{eq:PBPLOOP}
\end{equation}
where $N_s~(N_t)$ represents the number of lattice sites
in the spatial (temporal) direction and
$U_{4,t\mathbf{x}}$ is the temporal link variable.
The other important observable is
the ratio of
a scalar and a pseudo-scalar susceptibility~\cite{Kocic:1992is},
\begin{equation}
R_\pi \equiv \frac{\chi_\sigma}{\chi_\pi}
=
\frac{\partial \langle \bar \psi \psi\rangle/\partial m}
{\langle \bar \psi \psi\rangle/m}
=
\frac{\chi_{\mathrm{conn}} + \chi_{\mathrm{disc}}}
{\langle \bar \psi \psi\rangle/m}
\ ,\label{eq:R_pi}
\end{equation}
where
$
a^2\chi_{\mathrm{conn}} =
-N_f\langle \mathrm{Tr}[( MM )^{-1}] \rangle
/(4 N_s^3 N_t)
$
and
$
a^2\chi_{\mathrm{disc}}=
N_f^2
[\langle \mathrm{Tr} [M^{-1}]^2\rangle
-\langle \mathrm{Tr} M^{-1} \rangle^2]
/(16 N_s^3 N_t)
$.
Here, $R_{\pi}\sim\mathcal{O}(1)$ indicates
the scalar and pseudo-scalar degeneracy
attributed to the approximate chiral restoration,
while $R_{\pi}\ll 1$ indicates the chiral symmetry breaking
\cite{Deuzeman:2008sc,Kocic:1992is}.

\section{Results}
We evaluate the thermalized ensemble averages of
the chiral condensate (PBP) and Polyakov loop (PLOOP)
for various lattice couplings $\betaL$,
lattice sizes with the finite $T$ set up $N_s\gg N_t$,
and the number of flavors $N_f$.
All results have been obtained by using
a single value for a lattice bare fermion mass $am=0.02$.
Then we locate the lattice bare coupling $\betaLC$
associated with the chiral crossover
which is signaled by the drastic decrease (increase) of PBP (PLOOP)
as a function of $\betaL$.
In practice, the ratio of the scalar and pseudo-scalar susceptibility
$R_{\pi}$ gives a stronger signal of the chiral crossover,
owing to its renormalization invariant property.
In table \ref{Tab:bc},
we summarize the obtained critical lattice couplings
as a function of ($N_f,N_t$).
We have confirmed the (approximate) asymptotic scaling for
the normalized critical temperature
$T_c/\Lambda_{\mathrm{L}}$ varying $N_t$ at each $N_f$,
where $\Lambda_{\mathrm{L}}$ is so-called lattice Lambda.
This indicates that our $\betaLC$ have been determined
near to the continuum limit~\cite{Gupta:2000hr}.

\begin{table*}[ht]
\caption{
Summary of the (pseudo) critical lattice couplings $\betaLC$
for the theories with $N_f=0,~4,~6,~8$, $am=0.02$
and varying $N_t=4,~6,~8,~12$.
The result has partially been extracted from our recent paper
\protect\cite{Deuzeman:2008sc,Miura:2011mc}.
All results are obtained by using the action with
the same level improvements.
}\label{Tab:bc}
\begin{center}
\begin{tabular}{c|c|c|c|c}
\hline
$N_f\backslash N_t$ &
$4$&
$6$&
$8$&
$12$\\
\hline
$0$ &
$7.35\pm 0.1$&
$7.88\pm 0.05$&
$8.20\pm 0.1$&
$-$\\
$4$ &
$5.60\pm 0.1$&
$5.89\pm 0.05$&
$6.10\pm 0.1$&
$-$\\
$6$ &
$4.65\pm 0.05$&
$5.05\pm 0.05$&
$5.2\pm 0.05$&
$5.55\pm 0.1$\\
$8$ &
$-$&
$4.1125\pm 0.0125$&
$4.20\pm 0.1$&
$4.34\pm 0.04$\\
\hline
\end{tabular}
\end{center}
\end{table*}

We shall now discuss the emergence of the (pre-)conformality
by using our critical lattice $\gLC = \sqrt{10/\betaLC}$ collection.
The uniqueness of a critical temperature at each $N_f$,
$
T_c^{-1} = N_t\ a(\betaLC) = N_t^{\prime}\ a({\betaLC}^{\prime})
$
with $N_t\neq N_t^{\prime}$
gives a {\em thermal}
step-scaling $\Delta\betaLC = \betaLC - {\betaLC}^{\prime}$.
A vanishing (decreasing) $\Delta\betaLC$
can be the signal of the (pre-)conformality.

Our thermal step-scaling is a function of
$N_f$ and two lattice temporal extensions,
$\Delta\betaLC (N_f;N_t,N_t^{\prime})$,
and tends to be smaller with increasing $N_f$.
We here estimate the number of flavor satisfying
$\Delta\betaLC(N_f^*) = 0$
by extrapolating our $\betaLC$ collection into
the larger flavor region.
To this end, we plot the (pseudo) critical lattice coupling
$\gLC = \sqrt{10/\betaLC}$
as a function of $N_f$ in Fig.~\ref{Fig:MY},
which gives an extension of
Miransky-Yamawaki diagram \cite{Miransky:1997}
to finite $T$ cases.

Let us first pick up
the lattice critical couplings for
$N_f = 6$ and $8$,
and consider a ``constant $N_t$'' line.
As shown in the left panel of Fig.~\ref{Fig:MY},
$N_t = 6$ and $12$ lines
get to a joint at
$(\gLC, N_f^*) = (1.825\pm 0.02,\ 11.57\pm 0.17)$,
and $N_t = 6$ and $12$ lines
at $(\gLC, N_f^*) = (1.753\pm 0.02,\ 10.715\pm 0.17)$,
indicating the infra-red fixed point
with vanishing thermal step scalings.
Next we shall investigate
the critical lattice couplings at $N_t = 6$ and $8$
for whole range of $N_f = 0 - 8$.
They can be well fitted by assuming the functional expression
$
N_f(\gLC) = A\cdot \log~
\Bigl[B\cdot
\bigl(\gLC
- \gLC|_{N_f=0}
\bigr) + 1
\Bigr]
$
giving $(\gLC, N_f^*) = (1.88\pm 0.09,\ 11.39\pm 0.78)$
(right panel of Fig.~\ref{Fig:MY}).
Thus, the thermal step-scaling
with the use of our lattice critical couplings
supports the emergence of conformal window
near to $12$ flavor system,
whose (pre-)conformality is now under debate
in the recent lattice studies~\cite{DelDebbio:2010zz}.

\begin{figure}[ht]
\includegraphics[width=7.5cm]{./MY_fit_select.eps}
\includegraphics[width=7.5cm]{./MY_fit_log.eps}
\caption{(Pseudo) critical values of the lattice coupling
$\gLC=\sqrt{10/\betaLC}$ for theories with $N_f=0,~4,~6,~8$
and for several values of $N_t$
in the Miransky-Yamawaki phase diagram.
Left: We have picked up $\gLC = \sqrt{10/\betaLC}$ for
$N_f = 6$ and $8$,
and consider the ``constant $N_t$'' line
with $N_t = 6,\ 8,\ 12$. Candidates for the IRFP
have been estimated by the crossing points of
$N_t = 12$ line with $N_f = 6$ or $N_f = 8$ line.
Right:
The dashed line is a fit with the ansatz
$
N_f(\gLC) = A\cdot \log~
\bigl[B\cdot(\gLC- \gLC|_{N_f=0}) + 1\bigr]
$
for $N_t=6$ and $N_t=8$ results.}
\label{Fig:MY}
\end{figure}

As indicated by the constant $N_t$ line in Fig.~\ref{Fig:MY},
the critical coupling is an increasing function of $N_f$
for a fixed lattice temporal extension.
This behavior is the direct consequence
of enhanced fermion screening effects
due to the larger number of fermion species.
We have observed the thermal step
(or equivalently the asymptotic) scalings
for our $\betaLC$, thereby,
the enhancement of the screening effects
gets to a physical significance.
The test of the asymptotic scaling
has been the historical homework since
the pioneering finite $T$ study at large $N_f$
by J. Kogut and his collaborators~\cite{Kogut:1985pp}
and now we have completed it.

We shall now investigate the $N_f$ dependence
of the critical temperature.
In order to compare a physical quantity
such as a critical temperature
among different theories with a different $N_f$,
it is necessary to introduce
a $N_f$ independent reference scale by hand.
To set the reference scale, ideally speaking,
we would like to measure various sizes of Wilson loops
at various numbers of flavors in Monte-Carlo simulations
to obtain the running couplings $\bar{g}$
in the wide range of scales at each number of flavor.
In stead of performing such a massive simulation,
we approximately construct the renormalization flow
via the integral of the two-loop beta-function
$\beta(g)=-g^3(b_0 + b_1g^2)$
\begin{align}
\frac{T_c}{M(\gLR)}
= \frac{1}{N_t}\frac{a^{-1}(\gLC)}{M(\gLR)}
= \frac{1}{N_t}\int_{\gLR}^{\gLC}
\frac{dg}{-g^3(b_0 + b_1g^2)}
\ .\label{eq:Tc_M}
\end{align}
To specify the reference scale $M(\gLR)$,
we utilize our plaquette (tad-pole factor $u_0$) data
shown in the left panel of Fig.~\ref{Fig:TcM}.
Note that the plaquettes can be regarded as a kind of
renormalized couplings.
Let us consider a constant $u_0$ without $N_f$ dependences,
for instance $u_0 = 0.9$ in figure,
and read-off the corresponding bare lattice couplings at each $N_f$.
The obtained $\gL(N_f) = \sqrt{10/\betaL(N_f)}$ is used
as a reference coupling $\gLR$ in Eq.~(\ref{eq:Tc_M}).
This procedure imitates the scale setting
in the potential scheme renormalization,
and the use of $N_f$ independent $u_0$ is motivated by
the FRG scale setting method \cite{BraunGies}.

To be analogous to the FRG study,
we should choose $u_0$ so as to
get a UV $M(\gL)$ free from the chiral dynamics.
The middle panel of Fig.~\ref{Fig:TcM}
displays the $N_f$ dependence of
$T_c/M(\gLR)$
defined by Eq.~(\ref{eq:Tc_M}) with $u_0=0.9$.
Fitting $T_c/M(\gLR)$
with the FRG motivated ansatz
$T_c= K|N_f^* - N_f|^{(-2b_0^2/b_1)(N_f^*)}$,
we now read-off the lower edge of conformal window
$N_f^*\sim 9.47\pm 0.02$, which is somewhat smaller value
comparing to those obtained by the vanishing thermal scale settings.
In the middle panel of Fig.~\ref{Fig:TcM},
we find $T_c/M(\gLR)\ll 1$, indicating
the UV nature of the reference scale $M(\gLR)$.

To get more transparent view for the UV reference scale,
we here consider the particular reference coupling $\gLR$ -
the thermal critical coupling $\gTC$
which makes the reference scale $M(\gLR = \gTC)$
be equivalent to $T_c$ in Eq.~(\ref{eq:Tc_M}):
$T_c/M(\gTC) = 1$, giving a typical
interaction strength at $T_c$.
As shown in the right panel of Fig.~\ref{Fig:TcM},
the increasing nature of $\gTC$ indicates
a realization of
more strongly interacting non-Abelian plasma at larger $N_f$
as discussed in Ref.~\cite{Liao:2012tw}.
The criterion to set the UV reference scale
$M(\gLR)$ at every $N_f$ would be given by the condition
$\gLR\ll\gTC(N_f)$ for all $N_f$.
We find that $u_0\geq 0.84$ meets a requirement,
while the use of too large $u_0$ suffers from
the strong discritization errors.
In practice, we find that
the number of flavor giving the vanishing
$T_c/M(\gLR)$
is relatively stable within the range $0.84\leq u_0\leq 0.94$
which results in $9.85\geq N_f^*\geq 9.17$.

\begin{figure}[ht]
\includegraphics[width=5.3cm]{./u0_beta_all.eps}
\includegraphics[width=5.3cm]{./Tc_M_Nf_2_1_09.eps}
\includegraphics[width=5.3cm]{./gT_2_1.eps}
\caption{
Left:
The $\betaL$ dependences of the tadpole factor $u_0$
at zero temperature with the use of $12^4$ lattice.
Specifying a constant $u_0$ ({\em e.g.} $u_0=0.9$ in figure),
we read off the corresponding lattice couplings $\betaL$
which are used to define the scale $M$ at each
theory with $N_f$.
Middle:
The $N_f$ dependence of $T_c/M$ where
$M$ is the UV scale with $u_0=0.9$ at each theory with $N_f$.
The dashed line represents the fit for data
by using the FRG motivated ansatz
$T_c= K|N_f^* - N_f|^{(-2b_0^2/b_1)(N_f^*)}$
Right: The thermal critical coupling
with increasing $N_f$.}
\label{Fig:TcM}
\end{figure}

\section{Summary}
We have investigated
the (pre-)conformal dynamics
in color SU(3) gauge theories with
multi-species of fundamental fermions
by using the lattice Monte-Carlo simulation.
In order to study the conformality
beyond the fixed number of flavor $N_f$,
we have focused on
a reducing chiral dynamics at finite $T$
as a function of increasing $N_f$.
We have observed stronger fermion screenings
resulting from a larger fermion multiplicity
at larger $N_f$.
We have investigated a finite $T$ step-scaling
which follows from the uniqueness of
critical temperature ($T_c$) at each $N_f$,
then the vanishing step-scaling signals
the conformal dynamics at $N_f^* \sim 10 - 12$.
Further, motivated by the recent
FRG based studies \cite{BraunGies},
we have examined the $N_f$ dependence of $T_c$,
by introducing a UV $N_f$ independent reference scale $M(\gLR)$.
We have used the thermal critical coupling $\gTC$ as a criterion
to insure the UV nature of $M(\gLR)$
by imposing the condition $\gLR\ll\gTC(N_f)$ for all $N_f$.
We have found that
the number of flavor giving a vanishing $T_c/M(\gLR)$
is relatively stable within the range $0.84\leq u_0\leq 0.94$,
which results in $9.85\geq N_f^*\geq 9.17$.

As a future perspective,
we should measure various sizes of Wilson loops
at various numbers of flavors,
and perform more rigorous scale settings
in the potential scheme.
It is also mandatory to investigate the chiral limit
and the thermodynamic limit at large $N_f$.
This, together with a more extended set of flavor
numbers, will allow a quantitative analysis of the critical
behavior in the vicinity of the conformal IR fixed point.

\section*{Acknowledgments}
The author thanks Maria Paola Lombardo
and Elisabetta Pallante for
continuous discussions.
He thanks Edward Shuryak
for fruitful discussion during the xQCD workshop.
This work was in part based on the MILC Collaboration's public
lattice gauge theory code~\cite{MILC}.
The numerical calculations were carried out on the
IBM-SP6 and BG/P at CINECA,
Italian-Grid-Infrastructures in Italy,
and the Hitachi SR-16000 at YITP, Kyoto University in Japan.
\section*{References}


\section{Introduction}
Modern semiconductor nanodevices have reached critical device
dimensions in the range of several nanometers. These devices consist
of complicated two and three dimensional geometries composed of
multiple materials. Typically, about 10000 to 10 million atoms are
in the active device region with contacts controlling the current
injection. This finite extent suggests an atomistic, local and
orbital-based electronic structure representation. Quantitative
device design requires the reliable prediction of bandgaps and band
offsets within a few meV and effective masses at principal symmetry
points within a few percent.

Ab-initio methods that have no adjustable parameters offer such
atomistic representations. However, accurate models such as hybrid functionals~\cite%
{HSE06}, GW~\cite{Hybertsen_GW} and BSE
approximations~\cite{BSE_PRL} are computationally far too expensive
to be applied on multi-million atom devices. More approximate
ab-initio methods such as the local density approximations (LDA) and generalized gradient approximations (GGA)~\cite%
{RMP_DFT} do not reproduce band gaps, relative band offsets, and
effective masses accurately enough. Empirical methods such as the
empirical tight binding (ETB) method are numerically much more
efficient. The accuracy of ETB is hereby limited by the parameters
fitting. Previous ETB simulations in semiconductor nanodevices such
as resonant tunneling diodes~\cite{Klimeck_RTD}, quantum
dots~\cite{Klimeck_QuantumDot} and strained Si/SiGe quantum
wells~\cite{Valleysplitting_Boykin} showed quantitative agreement
with experiments.

The accuracy of the ETB method depends critically on the careful
calibration of the empirical parameters. The common way to determine
the ETB parameters is to fit ETB results to experimental band structures~\cite%
{Jancu_Tightbinding}~\cite{Boykin_TB_strain}.
One shortcoming of this method is its requirement of experimental
data that are often not available for new and exotic materials. In
addition, the ETB basis functions remain unknown, which makes it
notoriously difficult to predict wave function dependent quantities
with high precission. To overcome these shortcomings, some
approaches were developed to construct localized
basis functions from ab-initio results such as localized wannier functions~%
\cite{Marzari_Wannier_functions} or quasi-atomic orbitals~\cite%
{Qian_quasiatomic_orbitals}~\cite{Lu_quasiatomic_orbitals}~\cite%
{Urban_TBfromDFT}. Unfortunately, these functions are either not
reliably centered at atoms, or resulting Hamiltonian requires long
distance coupling
with large number of neighbors, which is numerically expensive. %
Nanoelectronic devices are increasingly based on complex
heterstructures in 1 or 2 dimensions including atomistc disorder and
strain.  In this realm models with one or two neighbors are
conceptually preferrable and simple to implement.  In contrast the
validity of many neighbor bulk basis states are very questionable in
this domain.

In this work, a DFT mapping method that constructs ETB parameters
from ab-initio calculations is presented. This method allows to
determine ETB basis functions that are centered at atoms and it
limits the interatomic coupling to the first or second nearest
neighbors. Since the method does not require experimental results,
the ETB parameters are less empirical. Two materials are considered
in this work: 1) the well known GaAs to validate the method and 2)
MgO which is recently used as a magnetic tunneling barrier material
in Magnetoresistive Random-Access-Memory devices. MgO lacks the
elaborate experimental analysis, but it is known to have small
spin-orbit interaction and a large band
gap~\cite{Fe_MgO_Fe_junctions}.

\begin{figure}[h]
\center
\includegraphics[width = 9cm]{flowchart}
\caption{The process of TB parameters construction from DFT
calculations.} \label{fig:fig_flowchart}
\end{figure}

Figure~\ref{fig:fig_flowchart} shows the flow chart of the
mapping method. The first step is to perform ab-initio calculations
of the band structure of a material. In general, any method that is
capable to calculate electronic band structures and wave functions
is suitable here. However, in this work, DFT calculations including
hybrid functionals corrections for 60 electronic bands performed
with VASP version 5.2~\cite{Vasp_HSE} are used. I
the second step, the ETB basis functions for each type of atom are
defined as
\begin{equation}
\Psi _{n,l,m}\left( \mathbf{r}\right) \equiv \Psi _{n,l,m}\left(
r,\theta ,\phi \right) =R_{n,l}\left( r\right) Y_{l,m}\left( \theta
,\phi \right) , \label{eq:definition_atomic_orbitals}
\end{equation}%
where the functions $Y_{l,m}$ are the complex spherical harmonics
and the
functions $R_{n,l}$ are exponentially damped plane waves%
\begin{equation}
R_{n,l}\left( r\right) =\sum_{i=1}^{N}\left[ a_{i}\sin \left( {\lambda _{i}r}%
\right) {+b_{i}\cos }\left( {{\lambda _{i}r}}\right) \right]
r^{n-1}\exp \left( -\alpha _{i}r\right) .
\label{eq:definition_basis_Radialpart}
\end{equation}%
Here, $\mathbf{r=}\left( r,\theta ,\phi \right) $ is the position
vector with the respective atom centered at the origin, $l$, $m$ are
the angular and magnetic quantum numbers of the orbital basis
function, and $n$ is the principal number of the atomic orbital as
used in H\"{u}ckel type basis functions~\cite{Cerd_ETH}. The
remaining parameters $a_{i},b_{i},\alpha _{i},\lambda _{i}$ are the
fitting parameters. The ETB basis functions are spin independent.
With a given set of ETB basis functions $\Psi _{TB}^{\mathbf{k}}$,

in the third step, a transformation matrix between these chosen DFT
basis functions
$\psi_{\textrm{DFT}}^{\mathbf{k}}$ and the
$\Psi _{ \textrm{TB}}^{\mathbf{k}}$ is calculated. Since the number
of the ETB basis functions is smaller than the DFT basis funcitons,
this transformation matrix is rectangular and represents a low rank
approximation.~\cite{LRA} Then, the DFT Hamiltonian is transformed
to the tight binding representation. The ETB Hamilton matrix
elements are approximated by two center integrals according to the
Slater-Koster
table~\cite{Slater_Tightbinding}~\cite{Podolskiy_TBElements}. Any
non-zero off-diagonal element of the overlap matrix is neglected.
ETB Hamilton matrix elements beyond either 1st or 2nd nearest
neighbor coupling are neglected. In Step 4, the band edges,
effective masses and eigen functions of the Hamiltonian at high
symmetry points are calculated and compared to the corresponding DFT
results. The overlaps of the ETB basis functions are also
determined. In the fifth step, all fitting parameters are adjusted
to improve 1) the agreement of the ETB\ results with the DFT ones
and 2) reduce the overlap matrix of the ETB basis functions to the
unity matrix. Steps 2 - 5 are repeated until the convergence
criterion is met, i.e. when the maximum difference of DFT and ETB band edges are within $%
10~\mathrm{meV}$, when the effective masses agree within 5\% and
when the eigenfunctions of DFT and ETB calculations agree by at
least 90\%. Step 6 requires to extract the converged ETB basis
functions and the ETB\ two center integrals.

\begin{figure}[h]
\center
\includegraphics[width = 11cm]{GaAs_TB_pub_apl_cmp}%
\caption{Band structure and density of states of GaAs by DFT(green
dashed lines), TB using parameters in
ref\protect\cite{Boykin_TB_strain} and TB using parameters by this
work.} \label{fig:fig_GaAsBand}
\end{figure}
For GaAs and MgO, eigenfunctions and eigen energies of the lowest 16
bands of the DFT calculations from $L$ to $\Gamma $ and $\Gamma $ to
$X$ are taken into account for the ETB fitting method. Here, the
wave functions for the topmost valence bands and lowest conduction
band valley are considered in the fitting of ETB eigenfunctions to
the DFT ones. The overlap of the ETB basis functions is partly
minimized in the fitting process. Most of the overlap matrix
elements vanish, but the maximum overlap i.e. in this case
the overlap of the $p$ orbitals of cations and $d$ orbitals of anions, e.g. $%
\langle p_{x}^{Ga}|d_{yz}^{As}\rangle =0.86$ at the $\Gamma $ point
remains comparably high.
\begin{figure}[t]
\center
\includegraphics[width = 10cm]{Ek_TE_GaAs_cmp}%
\caption{Comparison of E-k diagram of conduction bands (left figure)
and Transmission (right figure) of GaAs nanowire . The red lines are
results using TB parameters generated by this work; the blue line
are results using TB parameters by
ref\protect\cite{Boykin_TB_strain}.} \label{fig:fig_Ek_TE_GaAs}
\end{figure}
\begin{table}[t]
\begin{tabular}{c|cccc|ccc}
\hline
&  &  & GaAs &  &  & MgO &  \\
Quantity & DFT & TB & Error & ref\cite{Boykin_TB_strain} & DFT & TB
& Error
\\ \hline
$E_g(\Gamma)$ & 1.420 & 1.449 & $2.96\%$ & 1.413 & 7.831 & 7.499 &
$4.2\%$
\\
$E_g(X)$ & 1.973 & 1.947 & $0.9\%$ & 1.898 & 12.161 & 11.819 & $2.8\%$ \\
$E_g(L)$ & 1.728 & 1.718 & $0.6\%$ & 1.714 & 10.871 & 10.469 & $3.7\%$ \\
$m^*_{\Gamma}$ & 0.0692 & 0.0737 & $6.5\%$ & 0.0657 & 0.396 & 0.458 & $%
15.6\% $ \\
$m^*_{(X,l)}$ & 1.140 & 1.117 & $2.0\%$ & 1.881 & $-$ & $-$ & $-$ \\
$m^*_{(X,t)}$ & 0.219 & 0.231 & $5.5\%$ & 0.175 & $-$ & $-$ & $-$ \\
$m^*_{(L,l)}$ & 1.700 & 1.756 & $3.3\%$ & 1.728 & $-$ & $-$ & $-$ \\
$m^*_{(L,t)}$ & 0.133 & 0.138 & $3.8\%$ & 0.098 & $-$ & $-$ & $-$ \\
$\gamma_1 $ & 6.964 & 6.985 & $1.1\%$ & 7.388 & 0.952 & 0.889 & $6.6\%$ \\
$\gamma_2 $ & 2.084 & 2.151 & $3.6\%$ & 2.367 & 0.277 & 0.219 & $20.9\%$ \\
$\gamma_3 $ & 2.972 & 2.980 & $1.1\%$ & 3.098 & 0.376 & 0.234 & $37.7\%$ \\
\hline
\end{tabular}%
\caption{Comparison of important bandedges and effective masses of
GaAs and MgO by DFT and TB.} \label{tab:Target_comparison}
\end{table}GaAs is parameterized for the 1st nearest neighbor $sp^3d^5s^*$ ETB
model. The resulting parameters can be found in table
\ref{tab:Parameter_GaAs}.
The band structure and density of states of GaAs are shown in Fig.~\ref%
{fig:fig_GaAsBand}. Calculated bandstructures of the DFT method, of
the ETB method with parameters of the present mapping method and of
the ETB method with previously published
parameters~\cite{Boykin_TB_strain} are compared in
Fig.~\ref{fig:fig_GaAsBand} as well. The ETB calculations with new
parameters agree very well with the DFT results for energies below
$5$~$\mathrm{eV}$ whereas ETB calculations with previously published
parameters deviate from the DFT results already at about
$2$~$\mathrm{eV}$. Relevant band edges and effective masses are
compared in
table~\ref{tab:Target_comparison}, demonstrating a much better fit.

Figure~\ref{fig:fig_Ek_TE_GaAs} compares results of ETB calculations
of the band structure and the transmission coefficient of a GaAs $3.5~\mathrm{nm}%
\times 3.5~\mathrm{nm}$ squared nanowire when the new and the
previously published parameter set are used. The results agree for
the $\Gamma $ point at lower energies, but they deviate
significantly for the indirect conduction valley at about
$1.9~\mathrm{eV}$. The difference in the confinement energy of this
conduction valley originates from different transverse effective
masses at the L point of the two ETB parameter sets (see Table~\ref%
{tab:Target_comparison}). The modeling go the details of the high
symmetry points is particularly important for the emerging concept
of gamma-L transistors~\cite{LValley_Device}.
\begin{figure}[t]
\center
\includegraphics[width = 11cm]{MgO_TB_pub_apl_cmp}%
\caption{Band structure and density of states of MgO by DFT and TB.}
\label{fig:fig_Ek_TE_MgO}
\end{figure}

\begin{figure}[t]
\center
\includegraphics[width = 9cm]{TB_Basis_MgO
\caption{Radial part of Basis functions of Mg and O atoms in MgO.}
\label{fig:MgOBasisfunctions}
\end{figure}
MgO cyrstalizes in Rocksalt structure. Each Oxygen atom has 6
Magnesium atoms as 1st nearest neighbors and 12 Oxygen atoms as 2nd
nearest neighbors. The valence bands of MgO are formed by hybridized
orbitals of Oxygen atoms:
The $s$-orbitals and $p$-orbitals of Oxygen contribute to valence bands around $%
-17~\mathrm{eV}$ and $-1~\mathrm{eV}$ respectively, while the
$s^{\ast }$, $d$-orbitals of Oxygen and orbitals of Magnesium
contribute to the band structure for energies beyond $7~\mathrm{eV}$ (see Fig.~\ref%
{fig:fig_Ek_TE_MgO}). According to this, MgO is parameterized for a
2nd nearest neighbor $sp^3d^5s^*$ ETB\ model (for parameters see table \ref{tab:Parameter_MgO}). Within
this model, the interaction between two Oxygen atoms is required to
produce the correct valence bands. The interaction between two
Magnesium atoms is omitted, since the onsite energies of Magnesium orbitals are higher than $%
20~\mathrm{eV}$ so that the omission will only affect bands close to
that high and technically irrelevant energy. The ETB band structure
matches the DFT result well within the energies $5$ to
$20~\mathrm{eV}$. Important band properties calculated in the DFT
and the ETB method are listed in Table~\ref{tab:Target_comparison}.
The basis functions of Oxygen and Magnesium are shown in Fig.~\ref%
{fig:MgOBasisfunctions}. The basis functions of Oxygen are more
localized while the basis functions of Magnesium are more plane wave
like functions. This difference in the localization originates from
the fact that Magnesium orbital energies are higher than the Oxygen
ones.

In conclusion, a method to determine ETB parameters from density
function theory calculations is developed. The method is applied to
GaAs and MgO. First nearest neighbor ETB parameters and basis
functions for an $sp^3d^5s^*$ model of GaAs are presented. Results
of this parameterization agree well with the DFT calculations. The
new ETB parameters lead to lower indirect conduction bands when
applied to GaAs nanowires. Second nearest neighbor ETB parameters
and basis functions for an $sp^3d^5s^*$ model of MgO are also
obtained. The ETB\ results with this parameterization also agree
well with the DFT calculations.

nanoHUB.org computational resources operated by the Network for
Computational Nanotechnology funded by NSF are utilized in this
work. The research was funded by the Lockheed Martin Corporation and
NSF (Award No. 1125017)

\begin{table}[h]
\caption{\label{tab:Parameter_GaAs} Tight Binding Parameters for
bulk GaAs}
\begin{tabular}{|cc|cc|}
\hline
  Parameter & value & Parameter & value  \\
\hline
$E_{sa}$ & $-4.5863$ & $s^{*}_ap_c\sigma$ & $2.6877 $\\
$E_{pa}$ & $1.4694$ & $s^{*}_cp_a\sigma$ & $1.8335 $ \\
$E_{s*a}$ & $10.0480$ & $s_ad_c\sigma$ & $-2.1172$ \\
$E_{da}$ & $11.2878$ & $s_cd_a\sigma$ & $-2.9128$ \\
$E_{sc}$ & $-1.3323$ & $s^{*}_ad_c\sigma$ & $-0.4974 $\\
$E_{pc}$ & $9.5885$ &  $s^{*}_cd_a\sigma$ & $-2.9971$\\
$E_{s*c}$ & $25.6752$ & $pp\sigma$ &  $3.8065 $\\
$E_{dc}$ & $35.2863$ & $pp\pi$ & $-1.5010$ \\
$\Delta_{a}$ & $0.1259 $  & $p_ad_c\sigma$ & $-1.2077$ \\
$\Delta_{c}$ & $0.1235$  & $p_cd_a\sigma$ & $-1.9855$ \\
$ss\sigma$ &  $-1.7615$ & $p_ad_c\pi$ & $3.1547$ \\
$s^{*}s^{*}\sigma$ &  $-0.8374$ & $p_cd_a\pi$ & $2.3234$ \\
$s_a^{*}s_c\sigma$ & $-1.1173$ & $dd\sigma$ & $-1.9986$ \\
$s_as_c^{*}\sigma$ & $-2.9313$ & $dd\pi$ & $3.1681$  \\
$s_ap_c\sigma$ & $2.1768 $ & $dd\delta$ & $-2.3137 $ \\
$s_cp_a\sigma$ & $3.6705 $ & & \\
\hline
\end{tabular}
\end{table}

 $$
 $$

\begin{table}[h]
\caption{\label{tab:Parameter_MgO} Tight Binding Parameters for bulk
MgO}
\begin{tabular}{|cc|cc|}
\hline
Parameter &  Value &  Parameter & Value \\
\hline
$E_{sa}$ & $-7.1496$ & $p_ad_c\pi$ & $1.5284$ \\
$E_{pa}$ &  $5.8926$ & $p_cd_a\pi$ & $-4.0453$ \\
$E_{s*a}$ & $23.8138$ & $d_ad_c\sigma$ & $-1.0038$ \\
$E_{da}$ & $40.0285$ & $d_ad_c\pi$ & $5.0830$ \\
$E_{sc}$ & $38.4754$ & $d_ad_c\delta$ & $-0.6323$ \\
$E_{pc}$ & $32.1465$ & $s_as_a\sigma$  &   $-0.2718$ \\
$E_{s*c}$ & $45.5084$ & $s^{*}_as_a^{*}\sigma$  &   $-0.4690$ \\
$E_{dc}$ & $57.9865$ & $s_as_a^{*}\sigma$  &   $-0.0001$ \\
$\Delta_{a}$ & $0.0031$  & $s_ap_a\sigma$ &    $0.3388$ \\
$\Delta_{c}$ & $0.0298$  & $s^{*}_ap_a\sigma$   & $0.1965$ \\
$s_as_c\sigma$ & $-0.1192$ & $s_ad_a\sigma$ & $-0.3380$ \\
$s^{*}_as^{*}_c\sigma$ &  $1.6477$ & $s^{*}_ad_a\sigma$   & $-0.4407$  \\
$s_a^{*}s_c\sigma$ & $-0.6008$ & $p_ap_a\sigma$   & $0.4371$ \\
$s_as_c^{*}\sigma$ & $-0.6347$ & $p_ap_a\pi$ &   $ -0.0641$ \\
$s_ap_c\sigma$ & $-2.0013$ & $p_ad_a\sigma$ &   $ -0.4039$ \\
$s_cp_a\sigma$ &  $0.3283$ & $p_ad_a\pi$ &   $0.6986$ \\
$s^{*}_ap_c\sigma$ & $2.1584$ & $d_ad_a\sigma$ &   $-2.3768$ \\
$s^{*}_cp_a\sigma$ & $2.1282$ & $d_ad_a\pi$    &   $0.4556$ \\
$s_ad_c\sigma$ &  $0.6641$ & $d_ad_a\delta$ &   $0.0967$ \\
$s_cd_a\sigma$ & $-2.9483$ & & \\
$s^{*}_ad_c\sigma$ & $1.6890$ &  &\\
$s^{*}_cd_a\sigma$ & $3.1534$ &  &\\
$p_ap_c\sigma$ &  $0.1743$ & & \\
$p_ap_c\pi$ & $-0.4703$ &  &\\
$p_ad_c\sigma$ & $-1.9960$ &  &\\
$p_cd_a\sigma$ & $0.0519$ &  & \\
\hline
\end{tabular}
\end{table}

\section*{References}