\section{Introduction}

In the first decades of last century, physics encounters two great
revolutions; one in large scale structures' domain namely special and
general relativities and another in microscopic world i.e. quantum
mechanics. Since quantum mechanics deals with low speed microscopic
particles, physicists quickly thought about a new theory combining (at the
first step) special relativity and quantum mechanics called relativistic
quantum mechanics in order to describe high speed microscopic particles. In
this regime, we continuously encounter production and annihilation of
particles. Therefore, relativistic quantum mechanics was progressed and gave
birth to a more consistent theory i.e. quantum field theory (QFT) where the
Hilbert space is vast enough to include arbitrary number of particles. QFT
indicates impressive theoretical results confirmed up to many significant
digits by experiments in the case of weak interactions where it is allowed
to perform perturbation method \cite{Fey}. However, when coupling is strong,
one should follow other ways in order to compute physical quantities in the
context of QFT.

Fundamental particles physics is not the only place where QFT is employed.
In fact, the situations where one faces non-relativistic particles produced
and annihilated continuously are the other places that QFT methods can be
utilized. The most familiar situations with this behavior are condensed
matter systems. For instance, phonons in a lattice are produced and
annihilated. Also, electrons below or above Fermi level may pass it and
cause a hole production or annihilation. These events resemble the pair
production or annihilation. Another place where QFT methods are useful is
statistical field theory. When we have a system with continuous
configurations, the statistical mechanics of the system is formulated in the
context of statistical field theory. Accordingly, we deal with an
integration in the formula of central quantity of statistical mechanics
namely partition function instead of the usual summation. In the other hand,
the amplitude of the transition between two states in QFT interpreted as the
generating function is calculated through the path integral method and has
very similar formulation to the partition function of statistical field
theory. Hence, the methods of QFT can be effectively applied in these cases.
A very important example of such situation happens near a critical point
where we have a continuous phase transition.

As we mentioned above the common perturbation method does not work for
strong interactions in QFT and hence one should look for another methods to
calculate quantities in this case. One of these methods based on physics of
black holes (mostly thermodynamics of black holes) is gauge/gravity duality
(GGD). The first version of GGD is AdS/CFT correspondence presented by
Maldacena in 1998 \cite{Mald}. GGD connects the QFT (usually gauge field
theories) lives on $n$-dimensional boundary of the spacetime with black hole
solutions in $(n+1)$-dimensional gravity theories of the bulk. This
connection is performed by equality of generating function of QFT ($Z_{QFT}
) and the saddle-point approximation of bulk partition function in the case
of classical gravity. Through this, one can obtain the key quantity of QFT
namely generating function. Computation of almost all physical quantities in
QFT relies on generating function $Z_{QFT}$.

The idea of GGD has been applied frequently in order to analyse condensed
matter systems. For instance quantum Hall effect \cite{Hall}, fractional
quantum Hall effect \cite{FQH}, Nernst effect \cite{Ner} and superconductors
\cite{SupCond,SupCondLif} have been studied by this method. Furthermore,
there are interesting strongly correlated electronic and atomic systems and
also non-relativistic ones which possess Schrodinger symmetry \cite{Schr}.
Nevertheless, the near a critical point dynamics of such systems can be
described by a relativistic conformal field theory or a more subtle scaling
theory respecting the Lifshitz symmetry \cite{Lif
\begin{equation}
t\rightarrow \lambda ^{z}t,\text{ \ \ \ \ }\vec{\mathbf{x}}\rightarrow
\lambda \vec{\mathbf{x}}.
\end{equation
The spacetime which supports above symmetry on its $r$-infinity boundary is
\cite{Lif}
\begin{equation}
ds^{2}=-\frac{r^{2z}}{l^{2z}}dt^{2}+\frac{l^{2}dr^{2}}{r^{2}}+r^{2}d\vec
\mathbf{x}}^{2},  \label{lifmet}
\end{equation
known as Lifshitz spacetime. In (\ref{lifmet}), $z$ is dynamical critical
exponent. As regards the fact that thermodynamics of black holes plays a
central role in GGD, the literature has encountered increasing interests on
this issue. Solutions and thermodynamics of asymptotic Lifshitz black
objects (usually called Lifshitz solutions) in the presence of massive gauge
fields have been studied in \cite{masssol} and \cite{Deh1} respectively.
Also, thermodynamics of Lifshitz solutions of gravity models containing
higher curvature corrections has been considered under study \cite{hcc}.
Lifshitz solutions with toroidal event horizons and their thermal properties
in the framework of gravity theory coupled to Abelian gauge fields with
negative cosmological constant have been studied in \cite{toplif}. Moreover,
in Einstein-dilaton gravity in the presence of massless gauge fields,
thermal behavior of uncharged \cite{peet} and linearly charged \cite{tario}
Lifshitz black solutions have been investigated. In this context, also
nonlinearly charged Lifshitz solutions with power-law Maxwell field have
been taken under consideration \cite{Deh2}. Furthermore, thermodynamics of
Gauss-Bonnet-dilaton Lifshitz black branes has been discussed in \cite{Deh3
. The idea of holography has been applied for systems with Lifshitz symmetry
in different contexts including dilaton gravity from different points of
view \cite{hololif} and behaviors of some quantities such as viscosity \cit
{vislif} and different kinds of conductivity \cit
{SupCondLif,CondLif,1106,1211,1311,1402,1508} have been analyzed. The
authors of Ref. \cite{SupCondLif}, constructed an Abelian Higgs model in a
gravity background which is dual to a strongly coupled system at a
Lifshitz-fixed point and finite temperature. They also explored the
conductiviy via holographic techniques. The holographic conductivity in the
presence of an additional term in gauge field action including the Weyl
tensor coupled to Maxwell fied strengths in a domain-wall background whose
near horizon IR geometry is Lifshitz black hole and asymptotic geometry is
AdS has been calculated by performing both memberane paradigm and Kubo's
formula \cite{1106}. Employing the gauge/gravity duality the holographic
superconductors were explored with $z=2$ Lifshitz scaling \cite{1211}. The
effects of Lifshitz dynamical exponent $z$ on the holographic superconductor
models including $s$-wave and $p$-wave models have been disclosed both
numerically and analytically in \cite{1311}. Via holography, superconducting
phase transitions of a system dual to Yang-Mills field coupled to an axion
as probes of black hole with arbitrary Lifshitz scaling have been studied
with $p_{x}+ip_{y}$ condensate \cite{1402}. While the latter phase is known
as unstable in the relativistic case ($z=1$) and absence of axion field, in
the non-relativistic case ($z\neq 1$\textbf{) }with axion field, stability
of it has been studied and behavior of Hall conductivity in
non-superconducting pahase has been studied numerically as a function of
Lifshitz scaling \cite{1402}. Furthermore, an analytic computation for
longtitudinal DC conductivity corresponding to Lifshitz-like fixed points
based on linear response theory has been performed in the presence of chiral
anomalies \cite{1508} and an appropriate holographic set up has been
constructed to calculate the Lifshitz sector of the DC conductivity at
strong coupling and low charge density limit.

It is well-known that the behavior of conductivity have a close relation to
electrodynamics model under consideration. There are several nonlinear
electrodynamics models in addition to linear Maxwell one. The first one of
these models is Born-Infeld electrodynamics presented in 1930's \cite{BI}.
This model which arises in open superstrings and D-branes \cite{BI1,BI2},
removes the divergency of charged particles' self energy. The nonlinearity
of this model is determined by a parameter $\beta $ called nonlinear
parameter. Born-Infeld electrodynamics recovers the linear Maxwell model for
large $\beta $'s. Recently, another nonlinear electrodynamics model has been
introduced which is called exponential electrodynamics \cite{hendiexp}. The
large $\beta $ behavior of this model is similar to Born-Infeld
electrodynamics. The advantage of exponential electrodynamics compared to
Born-Infeld theory is that while it does not fully remove the divergency of
electric field at the source, it does make this divergency much weaker in
comparison with Maxwell field \cite{hendi2,hendi3}. This is more reasonable
compared to the Born-Infeld case, since \textit{near the origin} where
r\rightarrow0$, the electric field of a point-like charged particle should
be an increasing function. Besides, it has been observed that the
exponential nonlinear electrodynamics has crucial effects on condensation
and critical temperature of a holographic superconductor \cite{ZPCJ.NN}.
Furthermore, in recent studies on holographic superconductors, it has been
shown that the exponential nonlinear electrodynamics can increase the
critical values of the external magnetic field as the temperature goes to
zero \cite{SheAb}. In this paper, we consider Einstein-dilaton gravity in
the presence of exponential electrodynamics coupled to dilaton field. We
construct exact asymptotically Lifshitz black hole/brane solutions and
discuss the thermal behavior of them. Finally, we apply GGD to obtain
viscosity and conductivity corresponding to our solutions for ($3+1$) and (
2+1$)-dimensional systems on boundary respectively.

The layout of this paper is as follows. In Sec. \ref{Thermo}, we first find
the exact asymptotic Lifshitz topological black hole solutions of
Einstein-dilaton gravity in the presence of exponential electrodynamics
coupled to dilaton field. Then, we calculate the conserved and
thermodynamics quantities in order to check the satisfaction of
thermodynamics first law. At the end of Sec. \ref{Thermo}, we disclose the
effects of parameters of model on thermal stability of our solutions. In
Sec. \ref{GGD}, we first apply the GGD through pole method to obtain the
viscosity of a system lives on three-dimensions. In continue of Sec. \re
{GGD}, we study the holographic conductivity of two-dimensional systems for
both linear Maxwell and nonlinear exponential elctrodynamics. Finally, we
figure out the behavior of real and imaginary parts of conductivity for
asymptotic AdS and Lifshitz cases. The last section is devoted to present
summary and concluding remarks.

\section{Thermodynamics of Asymptotic Lifshitz Solutions}

\label{Thermo} In this section, we first find asymptotic Lifsthitz
solutions. Then, we study thermodynamics of our solutions and show that the
first law of thermodynamics is satisfied on the horizon. Finally, we
investigate thermal stability of the solutions in both canonical and
grand-canonical ensembles.

\subsection{Action and Higher-Dimensional Solutions}

The action of ($n+1$)-dimensional ($n\geq 3$) Einstein-dilaton gravity
coupled to a nonlinear electrodynamics and two linear Maxwell fields can be
written as
\begin{gather}
S=\int_{\mathcal{M}}d^{n+1}x\mathcal{L},  \label{Act} \\
\mathcal{L}=\frac{\sqrt{-g}}{16\pi }\left( \mathcal{R}-\frac{4}{n-1}(\nabla
\Phi )^{2}-2\Lambda +L(F,\Phi )-\sum\limits_{i=1}^{2}e^{-4/(n-1)\lambda
_{i}\Phi }H_{i}\right) ,  \notag
\end{gather
where $\mathcal{R}$ and $\Phi $ are Ricci scalar and dilaton scalar field,
respectively. Here $F=F^{\mu \nu }F_{\mu \nu }$ and $H_{i}=(H_{i})_{\mu \nu
}(H_{i})^{\mu \nu }$, where $F_{\mu \nu }=\partial _{\lbrack \mu }A_{\nu ]}
, $(H_{i})_{\mu \nu }=\partial _{\lbrack \mu }(B_{i})_{\nu ]}$, and $A_{\nu
} $ and $(B_{i})_{\nu }$ are electromagnetic vector potentials. $\lambda
_{i} $ and $\Lambda $ are some constants and $L(F,\Phi )$ is the Lagrangian
of nonlinear electrodynamic matter source. Varying the action (\ref{Act})
with respect to the metric $g_{\mu \nu }$, scalar field $\Phi $, and
electromagnetic vector potentials $A_{\nu }$ and $(B_{i})_{\nu }$, leads to
the following equations of motion,
\begin{eqnarray}
\mathcal{R}_{\mu \nu } &=&\frac{g_{\mu \nu }}{n-1}\left\{ 2\Lambda
+2L_{F}F-L(F,\Phi )-\sum\limits_{i=1}^{2}H_{i}e^{-4\lambda _{i}\Phi
/(n-1)}\right\}  \notag \\
&+&\frac{4}{n-1}\partial _{\mu }\Phi \partial _{\nu }\Phi -2L_{F}F_{\mu
\lambda }F_{\nu }^{\text{ \ }\lambda }+2\sum\limits_{i=1}^{2}e^{-4\lambda
_{i}\Phi /(n-1)}\left( H_{i}\right) _{\mu \lambda }\left( H_{i}\right) _{\nu
}^{\text{ \ }\lambda },  \label{FE1}
\end{eqnarray
\begin{eqnarray}
&&\nabla ^{2}\Phi +\frac{n-1}{8}L_{\Phi }+\sum\limits_{i=1}^{2}\frac{
\lambda _{i}}}{2}e^{-{4\lambda _{i}\Phi }/({n-1})}H_{i}=0,  \label{FE2} \\
&&\triangledown _{\mu }\left( L_{F}F^{\mu \nu }\right) =0,  \label{FE3} \\
&&\triangledown _{\mu }\left( e^{-{4\lambda _{i}\Phi }/({n-1})}\left(
H_{i}\right) ^{\mu \nu }\right) =0,  \label{FE4}
\end{eqnarray
where we have used the convention $X_{Y}=\partial X/\partial Y$. In this
paper, we intend to consider exponential nonlinear electrodynamics \cit
{hendiexp}. The Lagrangian of such type of nonlinear electrodynamics coupled
to the dilaton field in higher dimensions can be written as \cit
{SheyHaj,SheyKaz}
\begin{equation}
L(F,\Phi )=4\beta ^{2}e^{4\lambda \Phi /(n-1)}\left[ \exp \left( \frac
-e^{-8\lambda \Phi /(n-1)}F}{4\beta ^{2}}\right) -1\right] .  \label{ELag}
\end{equation
The behavior of exponential nonlinear electrodynamics coupled to dilaton
field (END) for large $\beta $, namely
\begin{equation}
\lim_{\beta \rightarrow \infty }L(F,\Phi )=-e^{-4\lambda \Phi /(n-1)}F+\frac
e^{-12\lambda \Phi /(n-1)}F^{2}}{8\beta ^{2}}+O\left( \frac{1}{\beta ^{4}
\right) ,  \label{Larbeta}
\end{equation
is similar to large $\beta $ behavior of the Born-Infeld electrodynamics
coupled to dilaton field (BID
\begin{equation}
\lim_{\beta \rightarrow \infty }L_{BI}(F,\Phi )=-e^{-4\lambda \Phi /(n-1)}F
\frac{e^{-12\lambda \Phi /(n-1)}F^{2}}{8\beta ^{2}}+O\left( \frac{1}{\beta
^{4}}\right) ,
\end{equation
wher
\begin{equation*}
L_{BI}(F,\Phi )=4\beta ^{2}e^{4\lambda \Phi /(n-1)}\left[ 1-\sqrt{1+\frac
e^{-8\lambda \Phi /(n-1)}F}{2\beta ^{2}}}\right] .
\end{equation*
This similarity makes the form chosen for END justifiable and also implies
that one can consider both END and BID as the Lagrangian of nonlinear
electrodynamics coupled to dilaton field.

In order to construct asymptotic Lifshitz topological black hole solutions
in higher dimensions, we use the metric
\begin{equation}
ds^{2}=-\frac{r^{2z}f(r)}{l^{2z}}dt^{2}+{\frac{l^{2}dr^{2}}{r^{2}f(r)}
+r^{2}d\Omega _{n-1}^{2},  \label{met}
\end{equation
where $z\geq 1$ is dynamical critical exponent, $d\Omega _{n-1}^{2}$ is an (
n-1$)-dimensional hypersurface with constant curvature $(n-1)(n-2)k$ and
volume $\omega _{n-1}$ and $f(r)\rightarrow 1$ as $r\rightarrow \infty $.
Lifshitz dilaton topological black hole solutions in the presence of linear
Maxwell and nonlinear power-law Maxwell fields have been studied in \cit
{tario} and \cite{kzsd}, respectively. The first term of (\ref{Larbeta})
recovers the linear case of electrodynamic Lagrangian of \cite{kzsd}.
Therefore, we expect that our solutions recover the results of this case.
Using the metric (\ref{met}), one can immediately integrate (\ref{FE3}) and
\ref{FE4}) and find
\begin{equation}
F_{rt}=\frac{qe^{4\lambda \Phi /(n-1)}}{r^{n-z}}\exp \left[ -\frac{1}{2
L_{W}\left( \varrho \right) \right] ,  \label{Frt}
\end{equation
\begin{equation}
\left( H_{i}\right) _{rt}=\frac{q_{i}e^{4\lambda _{i}\Phi /(n-1)}}{r^{n-z}},
\label{Hrt}
\end{equation
where $\varrho \equiv q^{2}l^{2z-2}/(\beta ^{2}r^{2n-2})$ and $L_{W}(x)
\mathrm{Lambert}W(x)$ is the Lambert function which satisfies the identity
\cite{abram}
\begin{equation}
L_{W}(x)e^{L_{W}(x)}=x,
\end{equation
and has the series expansion
\begin{equation}
L_{W}(x)=x-x^{2}+\frac{3}{2}x^{3}-\frac{8}{3}x^{4}+\cdots ,  \label{serexp}
\end{equation
which converges provided $\left\vert x\right\vert <1$. For $\beta
\rightarrow \infty $, $F_{rt}$ given in Eq. (\ref{Frt}) reduces to the
linear Maxwell case \cite{kzsd}
\begin{equation}
\lim_{\beta \rightarrow \infty }F_{rt}=\frac{qe^{4\lambda \Phi /(n-1)}}
r^{n-z}}-\frac{q^{3}l^{2z-2}e^{4\lambda \Phi /(n-1)}}{2r^{3n-z-2}\beta ^{2}
+O\left( \frac{1}{\beta ^{4}}\right) .
\end{equation
Substituting Eqs. (\ref{Frt}) and (\ref{Hrt}) into Eqs. (\ref{FE1}) and (\re
{FE2}), one can find the following field equations
\begin{gather}
\frac{(n-1)^{2}rf^{\prime }+(n-1)^{2}nf+4r^{2}f\Phi ^{\prime 2}}{2(n-1)l^{2}
+\Lambda -\frac{(n-1)(n-2)k}{2r^{2}}  \notag \\
+\sum\limits_{i=1}^{2}\frac{q_{i}^{2}e^{4\lambda _{i}\Phi /(n-1)}}
l^{2(1-z)}r^{2(n-1)}}+4\beta ^{2}e^{4\lambda \Phi /(n-1)}\Theta =0,
\label{fe1}
\end{gather}

\begin{gather}
\frac{(n-1)^{2}rf^{\prime }+(n-1)^{2}(n-2)f+2(n-1)^{2}fz-4r^{2}f\Phi
^{\prime 2}}{2(n-1)l^{2}}+\Lambda -\frac{(n-1)(n-2)k}{2r^{2}}  \notag \\
+\sum\limits_{i=1}^{2}\frac{q_{i}^{2}e^{4\lambda _{i}\Phi /(n-1)}}
l^{2(1-z)}r^{2(n-1)}}+4\beta ^{2}e^{4\lambda \Phi /(n-1)}\Theta =0,
\end{gather}

\begin{gather}
\frac{r^{2}f^{\prime \prime }+(2n+3z-3)rf^{\prime }+\frac{4}{n-1}r^{2}f\Phi
^{\prime 2}+(2z^{2}+2(n-2)z+(n-1)(n-2))f}{2l^{2}}+\Lambda  \notag \\
-\frac{(n-3)(n-2)k}{2r^{2}}-\sum\limits_{i=1}^{2}\frac{q_{i}^{2}e^{4\lambda
_{i}\Phi /(n-1)}}{l^{2(1-z)}r^{2(n-1)}}+4\beta ^{2}e^{4\lambda \Phi
/(n-1)}\Psi =0,
\end{gather}

\begin{equation}
\frac{r^{2}f^{\prime }\Phi ^{\prime }+\Phi ^{\prime }frz+nrf\Phi ^{\prime
}+fr^{2}\Phi ^{\prime \prime }}{l^{2}}-\sum\limits_{i=1}^{2}\frac
q_{i}^{2}\lambda _{i}e^{4\lambda _{i}\Phi /(n-1)}}{l^{2(1-z)}r^{2(n-1)}
-4\beta ^{2}\lambda e^{4\lambda \Phi /(n-1)}\Theta =0,  \label{fe4}
\end{equation
where

\begin{equation*}
\Theta =\frac{1}{2}+\frac{ql^{z-1}}{2\beta r^{n-1}}\left( \sqrt
L_{W}(\varrho )}-\frac{1}{\sqrt{L_{W}(\varrho )}}\right) ,
\end{equation*
and

\begin{equation*}
\Psi =\frac{1}{2}-\frac{ql^{z-1}}{2\beta r^{n-1}\sqrt{L_{W}(\varrho )}}.
\end{equation*
The solutions of Eqs. (\ref{fe1})-(\ref{fe4}) can be obtained as
\begin{equation}
{\Phi (r)}{=}\frac{(n-1)\sqrt{z-1}}{2}{\ln \left( \frac{r}{b}\right) ,}
\label{Phi}
\end{equation
\begin{eqnarray}
f(r) &=&1-\frac{m}{r^{n+z-1}}+\frac{(n-2)^{2}kl^{2}}{(n+z-3)^{2}r^{2}}
\notag \\
&&-\frac{8\beta ^{2}l^{2}b^{2z-2}}{(n-1)(n-z+1)r^{2z-2}}\left\{ \frac{1}{2}
\frac{ql^{z-1}(n-z+1)}{2\beta }r^{z-n-1}\int r^{1-z}\left( \sqrt
L_{W}(\varrho )}-\frac{1}{\sqrt{L_{W}(\varrho )}}\right) dr\right\} ,  \notag
\\
&&  \label{f}
\end{eqnarray
where $m$ is a parameter which is related to the total mass of black hole.
The above solutions will fully satisfy the field equations (\ref{fe1})-(\re
{fe4}), provided
\begin{gather}
\lambda =-\sqrt{z-1},\text{ }\lambda _{1}=\frac{n-1}{\sqrt{z-1}},\text{
\lambda _{2}=\frac{n-2}{\sqrt{z-1}},  \notag \\
q_{1}^{2}=\frac{(n+z-1)(z-1)b^{2(n-1)}}{2l^{2z}},  \notag \\
q_{2}^{2}=\frac{k(n-1)(n-2)(z-1)b^{2(n-2)}}{2(z+n-3)l^{2(z-1)}},  \notag \\
\Lambda =-\frac{(n+z-1)(n+z-2)}{2l^{2}}.  \label{constants}
\end{gather
Since in the case of $k=-1$, $q_{2}$ is imaginary (except for $z=1$), we
exclude this case and focus on the black hole ($k=1$) and black brane ($k=0
) solutions in the remaining part of this paper. It is notable to mention
that althogh at first glance it seems that constants (\ref{constants})
diverge for $z=1$, one should note that $\lambda _{i}\Phi $\ is finite for
this case and $H_{i}=0$. Hence, for $z=1$, action (\ref{Act}) reduces to
(A)dS action in the presence of a nonlinear electrodynamics field. One can
perform the integration in Eq. (\ref{f}) by using MATHEMATICA software. One
obtains
\begin{eqnarray}
f(r) &=&1-\frac{m}{r^{n+z-1}}+\frac{(n-2)^{2}kl^{2}}{(n+z-3)^{2}r^{2}}
\notag \\
&&-\frac{4\beta ^{2}l^{2}b^{2z-2}}{(n-1)(n-z+1)r^{2z-2}}+\frac{8\beta
^{2}l^{2}b^{2z-2}(2n-2)^{\left( z-5n+3\right) /\left( 2n-2\right) }}
(z-2)^{(z+3n-5)/(2n-2)}r^{n+z-1}}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}
\right) ^{\left( n-z+1\right) /\left( 2n-2\right) }  \notag \\
&&\times \left\{ 4(n-1)^{2}\left[ \Gamma \left( \frac{3n+z-5}{2n-2},\frac{2-
}{2n-2}L_{W}(\varrho )\right) -\Gamma \left( \frac{3n+z-5}{2n-2}\right)
\right] \right.  \notag \\
&&\left. -(z-2)^{2}\left[ \Gamma \left( \frac{z-n-1}{2n-2},\frac{2-z}{2n-2
L_{W}(\varrho )\right) -\Gamma \left( \frac{z-n-1}{2n-2}\right) \right]
\right\} ,  \label{ff}
\end{eqnarray
where $\Gamma (x,y)$ and $\Gamma (x)$ are gamma functions and are related to
each other with the relation
\begin{equation}
\Gamma (x,y)=\Gamma (x)-\frac{y^{x}}{x}\mathbf{F}(x,1+x,-y),
\label{relation}
\end{equation
where $\mathbf{F}(a,b,c)$ is the hypergeometric function \cite{abram}. At
the first look, it seems that the function (\ref{ff}) diverges for $z=2$.
However, the factor $(z-2)$ in denominator is removed when one uses Eq. (\re
{relation}). Indeed, one can reexpress (\ref{ff}) by using relation (\re
{relation}) in terms of the hypergeometric functions as
\begin{eqnarray}
f(r) &=&1-\frac{m}{r^{n+z-1}}+\frac{(n-2)^{2}kl^{2}}{(n+z-3)^{2}r^{2}}
\notag \\
&&-\frac{4\beta ^{2}l^{2}b^{2z-2}}{(n-1)(n-z+1)r^{2z-2}}+\frac{4\beta
^{2}l^{2}b^{2z-2}}{(n-1)r^{n+z-1}}\left( \frac{q^{2}l^{2z-2}}{\beta
^{2}L_{W}(\varrho )}\right) ^{\left( n-z+1\right) /\left( 2n-2\right) }
\notag \\
&&\times \left\{ \frac{L_{W}^{2}(\varrho )}{3n+z-5}\mathbf{F}\left( \frac
3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}(\varrho )\right)
\right.  \notag \\
&&\left. +\frac{1}{n-z+1}\mathbf{F}\left( \frac{z-n-1}{2n-2},\frac{z+n-3}
2n-2},\frac{z-2}{2n-2}L_{W}(\varrho )\right) \right\} .  \label{ff1}
\end{eqnarray
Using the fact that $f(r_{+})=0$ where $r_{+}$ is the outermost event
horizon, one can obtai
\begin{equation}
{m(r_{+})}=r_{+}^{n+z-1}+\frac{(n-2)^{2}kl^{2}}{(n+z-3)^{2}r_{+}^{3-n-z}}
\frac{8\beta ^{2}l^{2}b^{2z-2}\Xi }{(n-1)(n-z+1)r_{+}^{z-n-1}}  \label{mr+}
\end{equation
where
\begin{eqnarray}
\Xi &=&\frac{1}{2}-\frac{(n-z+1)}{2r_{+}^{n-z+1}}\left( \frac{q^{2}l^{2z-2}}
\beta ^{2}L_{W}\left( \varrho _{+}\right) }\right) ^{\left( n-z+1\right)
/\left( 2n-2\right) }  \notag \\
&&\times \left\{ \frac{L_{W}^{2}\left( \varrho _{+}\right) }{3n+z-5}\mathbf{
}\left( \frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left(
\varrho _{+}\right) \right) \right.  \notag \\
&&\left. +\frac{1}{n-z+1}\mathbf{F}\left( \frac{z-n-1}{2n-2},\frac{z+n-3}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \varrho _{+}\right) \right) \right\} ,
\end{eqnarray
and $\varrho _{+}=q^{2}l^{2z-2}/(\beta ^{2}r_{+}^{2n-2})$. The large $\beta $
limit of $f(r)$ can be obtained by using the series expansion (\ref{serexp
). The result is
\begin{eqnarray}
f(r) &=&1-\frac{m}{r^{n+z-1}}+\frac{(n-2)^{2}kl^{2}}{(n+z-3)^{2}r^{2}}
\notag \\
&&+\frac{2q^{2}b^{2z-2}l^{2z}}{\left( n-1\right) \left( n+z-3\right)
r^{2n+2z-4}}-\frac{q^{4}b^{2z-2}l^{4z-2}}{2\left( n-1\right) \left(
3n+z-5\right) \beta ^{2}r^{4n+2z-6}}+O\left( \frac{1}{\beta ^{4}}\right) .
\notag \\
&&  \label{flarbeta}
\end{eqnarray
As we mentioned before, for $\beta \rightarrow \infty $ the NED Lagrangian
reduces to linear Maxwell Lagrangian. Therefore, we expect that the large
\beta $ limit of $f(r)$ recover the metric function of topological Lifshitz
black holes in Maxwell theory \cite{kzsd}. One should note that this is
indeed the case and Eq. (\ref{flarbeta}) recover the linear case of $f(r)$
presented in \cite{kzsd} when $\beta \rightarrow \infty $. The behavior of
f(r)$ is depicted in Figs. \ref{fig1a} and \ref{fig1b} for different values
of $\beta $. Figs. \ref{fig1a} and \ref{fig1b} correspond to black brane (
k=0$) and black hole ($k=1$) solutions respectively.
\begin{figure*}[t]
\centering
\subfigure[$k=0$, $n=4$]{
   \label{fig1a}\includegraphics[width=.46\textwidth]{fig1a.eps}\qquad}}
\subfigure[$k=1$, $n=5$]{
   \label{fig1b}\includegraphics[width=.46\textwidth]{fig1b.eps}}
\caption{The behavior of $f(r)$ versus $r$ for $l=1.5$, $b=0.3$, $q=1.7$,
z=1.5$ and $m=1.6$.}
\label{fig1}
\end{figure*}
Inserting Eqs. (\ref{Phi}) and (\ref{constants}) into (\ref{Frt}), one find
\begin{equation}
F_{rt}=\frac{qb^{2z-2}}{r^{n+z-2}}\exp \left[ -\frac{1}{2}L_{W}\left(
\varrho \right) \right] .
\end{equation
The fact that our solutions are static implies $A_{t}=A_{t}\left( r\right)
. Hence, we can obtain the gauge potential using relation $A_{t}\left(
r\right) =\int F_{rt}dr$ a
\begin{eqnarray}
A_{t} &=&\mu +\frac{qb^{2z-2}(2n-2)^{(z-n-1)/(2n-2)}}
(2-z)^{(-5+3n+z)/(2n-2)}}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}}\right)
^{\left( 3-n-z\right) /\left( 2n-2\right) }  \notag \\
&&\times \left\{ 2(n-1)\left[ \Gamma \left( \frac{3n+z-5}{2n-2},\frac{2-z}
2n-2}L_{W}\left( \varrho \right) \right) -\Gamma \left( \frac{3n+z-5}{2n-2
\right) \right] \right.  \notag \\
&&\left. -(z-2)\left[ \Gamma \left( \frac{z+n-3}{2n-2},\frac{2-z}{2n-2
L_{W}\left( \varrho \right) \right) -\Gamma \left( \frac{z+n-3}{2n-2}\right)
\right] \right\} .  \label{At}
\end{eqnarray
It is remarkable to note that since from one side $L_{W}\left( \varrho
\right) \rightarrow 0$ as $r\rightarrow \infty $ and from another side
\Gamma \left( x,0\right) =\Gamma \left( x\right) $, one can easily check
that $A_{t}$ is finite at infinity and reduces to the constant $\mu $.
Although (\ref{At}) seems divergent for $z=2$, we can remove the factor
(z-2)$ in denominator by using (\ref{relation}) and restate $A_{t}$ a
\begin{eqnarray}
A_{t} &=&\mu -qb^{2z-2}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}L_{W}\left(
\varrho \right) }\right) ^{\left( 3-n-z\right) /\left( 2n-2\right) }  \notag
\\
&&\times \left\{ \frac{L_{W}\left( \varrho \right) }{3n+z-5}\mathbf{F}\left(
\frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left( \varrho
\right) \right) \right.  \notag \\
&&\left. +\frac{1}{n+z-3}\mathbf{F}\left( \frac{n+z-3}{2n-2},\frac{3n+z-5}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \varrho \right) \right) \right\} .
\label{At1}
\end{eqnarray
We require that $A_{t}$ vanishes at horizon $r=r_{+}$. Thus, $\mu $ can be
calculated a
\begin{eqnarray}
\mu &=&qb^{2z-2}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}L_{W}\left( \varrho
_{+}\right) }\right) ^{\left( 3-n-z\right) /\left( 2n-2\right) }  \notag \\
&&\times \left\{ \frac{L_{W}\left( \varrho _{+}\right) }{3n+z-5}\mathbf{F
\left( \frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left(
\varrho _{+}\right) \right) \right.  \notag \\
&&\left. +\frac{1}{n+z-3}\mathbf{F}\left( \frac{n+z-3}{2n-2},\frac{3n+z-5}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \varrho _{+}\right) \right) \right\} .
\end{eqnarray
In the next subsection we study thermodynamics of the solutions by computing
thermodynamic and conserved quantities.

\subsection{Thermodynamics of Lifshitz black holes/branes}

In this section, we study thermodynamics of lifshitz black hole and brane
solutions. Since, the mass is a fundamental quantity in studying the
thermodynamics, we start with its calculation. Using the modified
subtraction method of Brown and York (BY) \cite{modBY}, one can find mass
per unit volume $\omega _{n-1}$ as (see Ref. \cite{kzsd})
\begin{equation}
M=\frac{(n-1)m}{16\pi l^{z+1}},  \label{Mass}
\end{equation
where $m$ is the geometrical mass given by (\ref{mr+}) in terms of outer
horizon radius $r_{+}$. Using the Gauss law, one can also compute the charge
as
\begin{equation}
Q=\frac{\,{1}}{4\pi }\int r^{n-1}L_{F}F_{\mu \nu }n^{\mu }u^{\nu }d{\Omega },
\label{chdef}
\end{equation
where $n^{\mu }$ and $u^{\nu }$ are the unit spacelike and timelike normals
to a sphere of radius $r$ given as
\begin{equation*}
n^{\mu }=\frac{1}{\sqrt{-g_{tt}}}dt=\frac{l^{z}}{r^{z}\sqrt{f(r)}}dt,\text{
\ \ \ \ }u^{\nu }=\frac{1}{\sqrt{g_{rr}}}dr=\frac{r\sqrt{f(r)}}{l}dr.
\end{equation*
Therefore, the charge per unit volume $\omega _{n-1}$ is
\begin{equation}
Q=\frac{ql^{z-1}}{4\pi }.  \label{charge}
\end{equation
The electric potential $U$ is defined a
\begin{equation}
U=A_{\mu }\chi ^{\mu }\left\vert _{r\rightarrow \infty }-A_{\mu }\chi ^{\mu
}\right\vert _{r=r_{+}},  \label{Pot}
\end{equation
where $\chi =\partial _{t}$ is the null generator of the horizon. Therefore,
it is a matter of calculation to obtain $U$ by using (\ref{At1}) as
\begin{eqnarray}
U &=&qb^{2z-2}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}L_{W}\left( \varrho
_{+}\right) }\right) ^{\left( 3-n-z\right) /\left( 2n-2\right) }  \notag \\
&&\times \left\{ \frac{L_{W}\left( \varrho _{+}\right) }{3n+z-5}\mathbf{F
\left( \frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left(
\varrho _{+}\right) \right) \right.  \notag \\
&&\left. +\frac{1}{n+z-3}\mathbf{F}\left( \frac{n+z-3}{2n-2},\frac{3n+z-5}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \varrho _{+}\right) \right) \right\} .
\label{elecpot}
\end{eqnarray
The area law of the black hole entropy states that the entropy of a black
hole is the quarter of event horizon area \cite{Beck}. The entropy of almost
all kinds of black holes in Einstein gravity including dilaton black holes
is calculated by using this near universal law \cite{hunt}. Thus, the
entropy per unit volume $\omega _{n-1}$ of the Lifshitz black holes can be
found as
\begin{equation}
s=\frac{r_{+}^{n-1}}{4}.  \label{entropy}
\end{equation
The Hawking temperature is
\begin{equation*}
T=\frac{r_{+}^{z+1}f^{\prime }\left( r_{+}\right) }{4\pi l^{z+1}},
\end{equation*
which can be calculated as
\begin{eqnarray}
T &=&\frac{(n+z-1)r_{+}^{z}}{4\pi l^{z+1}}+\allowbreak \frac
(n-2)^{2}kl^{1-z}}{4\pi (n+z-3)r_{+}^{2-z}}  \notag \\
&&-\frac{2\beta ^{2}l^{1-z}b^{2z-2}}{\pi (n-1)r_{+}^{z-2}}\left[ \frac{1}{2}
\frac{ql^{z-1}}{2\beta r_{+}^{n-1}}\left( \sqrt{L_{W}\left( \varrho
_{+}\right) }-\frac{1}{\sqrt{L_{W}\left( \varrho _{+}\right) }}\right)
\right] .  \label{Temp}
\end{eqnarray
The large $\beta $ limit of temperature reproduces the temperature of
Einstein-Maxwell dilaton Lifshitz black holes \cite{kzsd}:
\begin{eqnarray}
T &=&\frac{(n+z-1)r_{+}^{z}}{4\pi l^{z+1}}+\allowbreak \frac
(n-2)^{2}kl^{1-z}}{4\pi (n+z-3)r_{+}^{2-z}}-\frac{q^{2}l^{z-1}b^{2z-2}}{2\pi
(n-1)r_{+}^{2n+z-4}}+\frac{q^{4}l^{3z-3}b^{2z-2}}{8\pi
(n-1)r_{+}^{4n+z-6}\beta ^{2}}+O\left( \frac{1}{\beta ^{4}}\right) ,  \notag
\\
&&
\end{eqnarray
as one expects. We need the Smarr-type mass formula in order to check the
satisfaction of first law of thermodynamics. Using (\ref{mr+}), (\ref{Mass
), (\ref{charge}) and (\ref{entropy}), the Smarr-type mass can be written as
\begin{eqnarray}
{M}\left( s,Q\right) &=&\frac{(n-1)\left( 4s\right) ^{(n+z-1)/(n-1)}}{16\pi
l^{z+1}}+\frac{\left( n-1\right) (n-2)^{2}k(4s)^{(n+z-3)/(n-1)}}{16\pi
l^{z-1}\left( n+z-3\right) ^{2}}  \notag \\
&&-\frac{\beta ^{2}b^{2z-2}(4s)^{(n-z+1)/(n-1)}\Pi }{2\pi l^{z-1}(n-z+1)},
\end{eqnarray
where
\begin{eqnarray}
\Pi &=&\frac{1}{2}-\frac{(n-z+1)}{2(4s)^{(n-z+1)/(n-1)}}\left( \frac{16\pi
^{2}Q^{2}}{\beta ^{2}L_{W}\left( \zeta \right) }\right) ^{\left(
n-z+1\right) /\left( 2n-2\right) }  \notag \\
&&\times \left\{ \frac{L_{W}\left( \zeta \right) ^{2}}{3n+z-5}\mathbf{F
\left( \frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left(
\zeta \right) \right) \right.  \notag \\
&&\left. +\frac{1}{n-z+1}\mathbf{F}\left( \frac{z-n-1}{2n-2},\frac{n+z-3}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \zeta \right) \right) \right\} ,
\end{eqnarray
and $\zeta =\pi ^{2}Q^{2}/\left( \beta ^{2}s^{2}\right) $. As $\beta
\rightarrow \infty $, the behavior of Smarr-type mass is
\begin{eqnarray}
M\left( s,Q\right) &=&\frac{(n-1)(4s)^{(n+z-1)/(n-1)}}{16\pi l^{z+1}}+\frac
(n-1)(n-2)^{2}k(4s)^{(n+z-3)/(n-1)}}{16\pi (n+z-3)^{2}l^{z-1}}  \notag \\
&&+\frac{2\pi Q^{2}b^{2z-2}(4s)^{(3-n-z)/(n-1)}}{(n+z-3)l^{z-1}}-\frac{8\pi
^{3}Q^{4}b^{2z-2}(4s)^{(5-3n-z)/(n-1)}}{(3n+z-5)l^{z-1}\beta ^{2}}+O\left(
\frac{1}{\beta ^{2}}\right) .  \notag \\
&&
\end{eqnarray
This satisfy our expectation that the mass of the linear Maxwell case should
be recovered \cite{kzsd}. Now, in order to check the satisfaction of the
first law of thermodynamics, we take $S$ and $Q$\ as a complete set of
extensive quantities for mass $M(s,Q)$. Then, we define their conjugate
intensive quantities as temperature $T$\ and electric potential $U$ that
implies
\begin{equation}
T=\left( \frac{\partial M}{\partial s}\right) _{Q}\text{ \ \ \ \ and \ \ \ \
}U=\left( \frac{\partial M}{\partial Q}\right) _{s}.  \label{intqua}
\end{equation
Intensive quantities computed by (\ref{intqua}) are in coincidance with ones
obtained by (\ref{elecpot}) and (\ref{Temp}). Hence, the first law of
thermodynamics
\begin{equation}
dM=Tds+UdQ.
\end{equation
is satisfied. In the next subsection we investigate the stability of this
thermodynamic system under thermal perturbations.

\subsection{Thermal Stability in the Canonical and Grand-Canonical Ensembles}

\begin{figure*}[t]
\centering
\subfigure[~]{
   \label{fig2a}\includegraphics[width=.46\textwidth]{fig2a.eps}\qquad}}
\subfigure[~]{
   \label{fig2b}\includegraphics[width=.46\textwidth]{fig2b.eps}}
\caption{The behaviors of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$
and $T$ versus $r_{+}$ for $k=0$ with $l=b=1$, $q=1.1$, $z=3$ and $n=5$.}
\label{fig2}
\end{figure*}
\begin{figure*}[t]
\centering
\subfigure[~]{
   \label{fig3a}\includegraphics[width=.46\textwidth]{fig3a.eps}\qquad}}
\subfigure[~]{
   \label{fig3b}\includegraphics[width=.46\textwidth]{fig3b.eps}}
\caption{The behaviors of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$
and $T$ versus $z$ for $k=0$ with $l=b=1$, $q=0.5$, $r_{+}=0.7$ and $n=5$.}
\label{fig3}
\end{figure*}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig4.eps}}
\caption{The behavior of $10^{4}\mathbf{H}_{s,Q}^{M}$ and $T$ versus $r_{+}$
for $k=0$ with $l=2$, $b=1$, $q=0.2$, $n=5$, $z=6$ and $\protect\beta =0.1$.}
\label{fig4}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig5.eps}}
\caption{The behavior of $\mathbf{H}_{s,Q}^{M}$ and $10T$ versus $z$ for
k=0 $ with $l=b=1$, $q=0.2$, $n=5$, $r_{+}=0.7$ and $\protect\beta =0.1$.}
\label{fig5}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig6.eps}}
\caption{The behavior of $\mathbf{H}_{s,Q}^{M}$ and $T$ versus $\protec
\beta $ for $k=0$ with $l=b=1$, $q=0.2$, $n=5$, $z=6.5$ and $r_{+}=0.7$.}
\label{fig6}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig7.eps}}
\caption{The behavior of $10^{2}(\partial ^{2}M/\partial s^{2})_{Q}$, $1
\mathbf{H}_{s,Q}^{M}$ and $T$ versus $r_{+}$ for $k=1$ with $l=7$, $b=1$,
q=0.5$, $n=4$, $z=1.5$ and $\protect\beta =3$.}
\label{fig7}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig8.eps}}
\caption{The behavior of $10^{-1}(\partial ^{2}M/\partial s^{2})_{Q}$,
10^{-2}\mathbf{H}_{s,Q}^{M}$ and $T$ versus $z$ for $k=1$ with $l=b=1$,
q=0.1$, $n=5$, $r_{+}=0.6$ and $\protect\beta =3$.}
\label{fig8}
\end{figure}
\begin{figure*}[t]
\centering
\subfigure[~]{
   \label{fig81a}\includegraphics[width=.46\textwidth]{fig81a.eps}\qquad}}
\subfigure[~]{
   \label{fig81b}\includegraphics[width=.46\textwidth]{fig81b.eps}}
\caption{The behaviors of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$
and $T$ versus $r_{+}$ for $k=1$ with $l=b=1$, $q=1.1$, $z=3$ and $n=5$.}
\label{fig8+1}
\end{figure*}
\begin{figure*}[t]
\centering
\subfigure[~]{
   \label{fig82a}\includegraphics[width=.46\textwidth]{fig82a.eps}\qquad}}
\subfigure[~]{
   \label{fig82b}\includegraphics[width=.46\textwidth]{fig82b.eps}}
\caption{The behaviors of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$
and $T$ versus $z$ for $k=1$ with $l=b=1$, $q=1.1$, $r_{+}=0.8$ and $n=5$.}
\label{fig8+2}
\end{figure*}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig9.eps}}
\caption{The behavior of $(\partial ^{2}M/\partial s^{2})_{Q}$ and $T$
versus $\protect\beta $ for $k=1$ with $l=1$, $b=0.3$, $q=0.1$, $n=5$,
r_{+}=0.5$ and $z=1.3$.}
\label{fig9}
\end{figure}
\begin{figure}[h]
\centerline{\includegraphics[width=.46\textwidth]{fig10.eps}}
\caption{The behavior of $10^{-2}\mathbf{H}_{s,Q}^{M}$ and $T$ versus
\protect\beta $ for $k=1$ with $l=0.72$, $b=0.8$, $q=0.4$, $n=6$, $r_{+}=0.5$
and $z=3$.}
\label{fig10}
\end{figure}
In the first part of this section, we showed that one can regard a black
hole solution as a thermodynamic system. On the other hand, it is necessary
to investigate the stability of a thermodynamic system under thermal
perturbations. Therefore, this subsection is devoted to study the thermal
stability of the obtained solutions of the previous sections. The formal way
to analyze the stability of a thermodynamic system with respect to small
variations of the thermodynamic coordinates is by investigating the behavior
of the entropy $s(M,Q)$\ around equilibrium. This analysis can also be
performed in terms of the Legendre transformation of entropy namely $M(s,Q)
.\ In any ensemble, the local stability requires that the energy $M(s,Q)$\
be a convex function of its extensive variable \cite{Cal2,Gub}. The number
of thermodynamic variables is ensemble-dependent. In the canonical ensemble,
since the charge is fixed, the positivity of $(\partial ^{2}M/\partial
s^{2})_{Q}$\ in the ranges where temperature $T$ is positive suffices to
ensure the local stability. For the Lifshitz black solutions we have,
\begin{eqnarray}
\left( \frac{\partial ^{2}M}{\partial s^{2}}\right) _{Q} &=&\frac
z(n-1+z)r_{+}^{z-n+1}}{(n-1)\pi l^{z+1}}+{\frac{k\left( n-2\right)
^{2}(z-2)r_{+}^{z-n-1}}{\pi l^{z-1}(n-1)\left( z+n-3\right) }}  \notag \\
&&+\frac{8(z-2)\beta ^{2}r_{+}^{3-n-z}}{\pi {(n-1)}^{2}l^{z-1}b^{2-2z}
\left\{ \frac{1}{2}+\frac{ql^{z-1}r_{+}^{1-n}}{2\beta (z-2)}\left[ (n+z-3
\sqrt{L_{W}\left( \varrho _{+}\right) }-\frac{z-2}{\sqrt{L_{W}\left( \varrho
_{+}\right) }}\right] \right\} .  \notag \\
&&  \label{dMSS}
\end{eqnarray
The behavior of (\ref{dMSS}) for large $\beta $ is
\begin{eqnarray}
\left( \frac{\partial ^{2}M}{\partial s^{2}}\right) _{Q} &=&\frac
z(n-1+z)r_{+}^{z-n+1}}{(n-1)\pi l^{z+1}}+{\frac{k\left( n-2\right)
^{2}(z-2)r_{+}^{z-n-1}}{\pi l^{z-1}(n-1)\left( z+n-3\right) }}+\frac
2q^{2}(2n+z-4)b^{2z-2}r_{+}^{5-3n-z}}{\pi (n-1)^{2}l^{1-z}}  \notag \\
&&-\frac{q^{4}(4n+z-6)b^{2z-2}r_{+}^{7-5n-z}}{2\pi (n-1)^{2}l^{3-3z}\beta
^{2}}+O(\frac{1}{\beta ^{4}}),  \label{lbetadMdS}
\end{eqnarray
which recovers the result of \cite{kzsd} for the linear Maxwell
electrodynamics, as expected. Before we turn to grand canonical case it is
remarkable to find $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$for
highly nonlinear case i.e. $\beta \rightarrow 0$. In this case the behavior
of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$is
\begin{eqnarray}
\lim_{\beta \rightarrow 0}\left( \frac{\partial ^{2}M}{\partial s^{2}
\right) _{Q} &=&\frac{z(n-1+z)r_{+}^{z-n+1}}{(n-1)\pi l^{z+1}}+{\frac
k\left( n-2\right) ^{2}(z-2)r_{+}^{z-n-1}}{\pi l^{z-1}(n-1)\left(
z+n-3\right) }}  \notag \\
&&+\frac{4q\beta r_{+}^{4-2n-z}(n+z-3)}{\pi {(n-1)}^{2}b^{2-2z}}\sqrt
L_{W}\left( \frac{q^{2}l^{2z-2}}{\beta ^{2}r_{+}^{2n-2}}\right) }.
\label{sbetadMdS}
\end{eqnarray
Soon we use (\ref{sbetadMdS}) to discuss thermal stability of highly
nonlinear solutions. In the grand-canonical ensemble, since $Q$\ is not a
fixed parameter, the local stability requires the positivity of determinant
of Hessian matrix $\mathbf{H}_{sQ}^{M}=\left[ \partial ^{2}M/\partial
s\partial Q\right] >0$. The determinant of the Hessian matrix $\mathbf{H
_{s,Q}^{M}$ for the solutions under consideration can be calculated as
\begin{eqnarray}
\mathbf{H}_{s,Q}^{M} &=&\frac{4\beta b^{2z-2}l^{3-3z}}{q(n-1)^{2}r_{+}^{n-3}
\sqrt{L_{W}\left( \eta _{+}\right) }\left[ \Upsilon +\frac{4q\beta
(z-2)b^{2z-2}l^{z-1}}{(n-1)r_{+}^{n+2z-3}}\left( \sqrt{L_{W}\left( \eta
_{+}\right) }-\frac{1}{\sqrt{L_{W}\left( \eta _{+}\right) }}\right) \right]
\notag \\
&&-\frac{16b^{2z-2}l^{2-2z}(z-2)}{r_{+}^{n-z-1}(2n-2)^{2}}\left( \frac{\beta
^{2}L_{W}\left( \eta _{+}\right) }{q^{2}l^{2z-2}}\right) ^{(n+z-3)/(2n-2)}
\notag \\
&&\times \left[ \Upsilon +\frac{4q\beta b^{2z-2}l^{z-1}}{(n-1)r_{+}^{n+2z-3}
\left( (n+z-3)\sqrt{L_{W}\left( \eta _{+}\right) }-\frac{z-2}{\sqrt
L_{W}\left( \eta _{+}\right) }}\right) \right]  \notag \\
&&\times \left\{ \frac{L_{W}\left( \eta _{+}\right) }{3n+z-5}\mathbf{F
\left( \frac{3n+z-5}{2n-2},\frac{5n+z-7}{2n-2},\frac{z-2}{2n-2}L_{W}\left(
\zeta \right) \right) \right.  \notag \\
&&\left. +\frac{1}{n+z-3}\mathbf{F}\left( \frac{n+z-3}{2n-2},\frac{3n+z-5}
2n-2},\frac{z-2}{2n-2}L_{W}\left( \zeta \right) \right) \right\} ,
\end{eqnarray
wher
\begin{equation*}
\Upsilon =\frac{z(n+z-1)}{l^{2}}+\frac{(z-2)(n-2)^{2}k}{(n+z-3)r_{+}^{2}}
\frac{4(z-2)\beta ^{2}b^{2z-2}}{(n-1)r_{+}^{2z-2}}.
\end{equation*
The large $\beta $ behavior of $\mathbf{H}_{sQ}^{M}$ is
\begin{eqnarray}
\mathbf{H}_{s,Q}^{M} &=&\frac{4z(n+z-1)b^{2z-2}r_{+}^{4-2n}}
(n-1)(n+z-3)l^{2z}}+\frac{4(n-2)^{2}(z-2)kb^{2z-2}r_{+}^{2-2n}}
(n-1)(n+z-3)^{2}l^{2z-2}}  \notag \\
&&-\frac{8q^{2}(z-2)b^{4z-4}r_{+}^{8-4n-2z}}{(n-1)^{2}(n+z-3)}+\frac
2q^{2}b^{2z-2}}{(n-1)(3n+z-5)\beta ^{2}r_{+}^{4n-4}}  \notag \\
&&\times \left( \frac{q^{2}(z-2)b^{2z-2}l^{2z-2}}{(n-1)(n+z-3)r_{+}^{2n+2z-6
}-\frac{3k(n-2)^{3}}{n+z-3}-\frac{3z(n+z-1)r_{+}^{2}}{l^{2}}\right) +O\left(
\frac{1}{\beta ^{4}}\right) ,  \notag \\
&&
\end{eqnarray
which coincides with the one of linear Maxwell case \cite{kzsd}. In what
follows we discuss the stability for $k=0$ (black branes) and $k=1$
(spherical black holes), separately.

$k=0$: In this case, calculations show that the solutions are always stable
in the canonical ensemble. Although it is difficult to see this for
arbitrary $\beta $ from (\ref{dMSS}), we can show this fact for large and
small $\beta $'s. For large $\beta $ where the linear Maxwell regime is
dominant, this fact that the solutions are always stable in canonical
ensemble has been pointed out in \cite{kzsd}. For small $\beta $, the
behaviour of $\left( \partial ^{2}M/\partial s^{2}\right) _{Q}$ has been
given in (\ref{sbetadMdS}). It is obvious from (\ref{sbetadMdS}) that
\lim_{\beta \rightarrow 0}\left( \partial ^{2}M/\partial s^{2}\right) _{Q}>0$
for $k=0$. Figs. \ref{fig2a} and \ref{fig3a} depict the fact that black
branes are always stable in canonical ensemble. Of course, one should check
the positivity of temperature for these choices. This can be seen in Figs.
\ref{fig2b} and \ref{fig3b} which show that the temperature is positive for
these choices and therefore we have black branes. In grand canonical
ensemble, the system is thermally stable provided that the radius of the
black hole is larger than $r_{+\min }$ (Fig. \ref{fig4}). Figure \ref{fig5}
shows that there is a $z_{\max }$ that for values greater than it we
encounter instability. We also have a $\beta _{\max }$ that black branes are
stable under thermal perturbations for values lower than it as one can see
in Fig. \ref{fig6}.

$k=1$: As one can see from Fig. \ref{fig7}, there is a Hawking-Page phase
transition between small and large black holes in both ensembles. In terms
of $z$ we have an ensemble-dependent minimum $z_{\min }$ which depends on
the parameters of the system and black holes are unstable for values lower
than it as Fig. \ref{fig8} shows. The value of $z_{\min }$ is greater in
grand-canonical ensemble compare to that in canonical ensemble. For highly
nonlinear case (small $\beta $) and also linear case (large $\beta $) black
holes are stable for $z\geq 2$ in canonical ensemble as one can see from
\ref{lbetadMdS}) and (\ref{sbetadMdS}). Hence, one can conclude that there
are always stable black holes for $z\geq 2$ in canonical ensemble (Figs. \re
{fig8+1} and \ref{fig8+2}). Fig. \ref{fig9} shows that black holes are
unstable in canonical ensemble for $\beta <\beta _{\min }$ in the case of
z<2$. In grand-canonical ensemble we have stable black holes for $\beta
<\beta _{\max }$ (Fig. \ref{fig10}).

\section{Gauge/Gravity Duality\label{GGD}}

In this section, we would like to perform the gauge/gravity duality idea on
our obtained solutions. We first try to calculate zero-frequency shear
viscosity for a hydrodynamic system. Then, we employ this idea to study the
behavior of the holographic conductivity for a two dimensional system for
both asymptotic AdS ($z=1$) and asymptotic Lifshitz cases and present
experimental observations matched with obtained results.

\subsection{Holographic Viscosity}

Here we intend to calculate the ratio of shear viscosity to entropy $\eta /s$
in the zero frequency limit. For this aim, we use the pole method \cit
{vis1,vis2}. We consider the five-dimensional planar metri
\begin{equation}
ds^{2}=-\frac{r^{2z}f(r)}{l^{2z}}dt^{2}+{\frac{l^{2}dr^{2}}{r^{2}f(r)}
+r^{2}(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}),  \label{plmet}
\end{equation
where the QFT lives on the $4$-dimensional $r$-infinity boundary. Defining
y=1-r_{+}^{2}/r^{2}$, one can rewrite (\ref{plmet}) a
\begin{equation}
ds^{2}=-\frac{F(y)r_{+}^{2z}}{l^{2z}(1-y)^{z}}dt^{2}+{\frac{l^{2}dy^{2}}
4F(y)(1-y)^{2}}}+\frac{r_{+}^{2}}{1-y}(dx_{1}^{2}+dx_{2}^{2}+dx_{3}^{2}),
\label{plmet2}
\end{equation
where

\begin{eqnarray*}
F(y) &=&1-\frac{m(1-y)^{(3+z)/2}}{r_{+}^{3+z}}-\frac{4\beta
^{2}l^{2}b^{2z-2}(1-y)^{z-1}}{3(5-z)r_{+}^{2z-2}}+\frac{4\beta
^{2}l^{2}b^{2z-2}(1-y)^{(3+z)/2}}{3r_{+}^{3+z}}\left( \frac{q^{2}l^{2z-2}}
\beta ^{2}L_{W}(\varrho _{y})}\right) ^{\left( 5-z\right) /6} \\
&&\times \left\{ \frac{L_{W}^{2}(\varrho _{y})}{7+z}\mathbf{F}\left( \frac
7+z}{6},\frac{13+z}{6},\frac{z-2}{6}L_{W}(\varrho _{y})\right) +\frac{1}{5-z
\mathbf{F}\left( \frac{z-5}{6},\frac{z+1}{6},\frac{z-2}{6}L_{W}(\varrho
_{y})\right) \right\} ,
\end{eqnarray*
$\varrho _{y}\equiv q^{2}l^{2z-2}(1-y)^{3}/(\beta ^{2}r_{+}^{6})$ and
horizon located at $y=0$. In order to use the pole method, we first apply an
off-shell perturbatio
\begin{equation}
dx_{i}\rightarrow dx_{i}+\varepsilon e^{-i\omega t}dx_{j},
\end{equation
where $\varepsilon $ is an infinitesimal positive parameter. Then, by using
the perturbed metric and the expansion of $F(y)$ near the horizon, namely
\begin{eqnarray}
F(y) &=&\left[ \frac{z+3}{2}-\frac{2\beta ^{2}b^{2z-2}l^{2}}{3r_{+}^{2z-2}}
\frac{2q\beta b^{2z-2}l^{z+1}}{3r_{+}^{2z+1}}\left( \sqrt{L_{W}\left(
\varrho _{+}\right) }-\frac{1}{\sqrt{L_{W}\left( \varrho _{+}\right) }
\right) \right] y \\
&&-\left[ \frac{(z+3)(z+1)}{8}-\frac{(3z+2)q^{2}b^{2z-2}l^{2z}}{6r_{+}^{2z+4
}\exp \left[ -\frac{1}{2}L_{W}\left( \varrho _{+}\right) \right] \right.
\notag \\
&&\left. +\frac{(3z-1)\beta ^{2}b^{2z-2}l^{2}}{6r_{+}^{2z-2}}\left( \exp
\left[ -\frac{1}{2}L_{W}\left( \varrho _{+}\right) \right] -1\right) \right]
y^{2}+O\left( y^{3}\right) ,
\end{eqnarray
we calculate the residue of the pole at $y=0$ in the Lagrangian densit
\begin{equation}
\text{\textrm{Res}}_{y=0}\mathcal{L}=-\frac{\varepsilon ^{2}\omega
^{2}r_{+}^{3}e^{-2i\omega t}}{128\pi ^{2}T},  \label{ResL}
\end{equation
where $T$ is given by Eq. (\ref{Temp}). Now, via the formula \cite{vis1,vis2
\begin{equation}
\eta =-8\pi T\lim_{\varepsilon ,\omega \rightarrow 0}\frac{\text{\textrm{Res
}_{y=0}\mathcal{L}}{\varepsilon ^{2}\omega ^{2}},
\end{equation
the shear viscosity for zero frequency case can be obtained a
\begin{equation}
\eta =\frac{r_{+}^{3}}{16\pi }.
\end{equation
Therefore, the ratio $\eta /s$ is obtained as
\begin{equation}
\frac{\eta }{s}=\frac{1}{4\pi },  \label{shvis}
\end{equation
where we have used Eq. (\ref{entropy}). Therefore, we regain the well-known
value $1/4\pi $ for $\eta /s$ as many cases in Einstein gravity. This
implies that neither non-AdS symmetry of the system on boundary nor the
presence of the dilaton and additional gauge fields cannot affect the value
of the shear viscosity of the system.

\subsection{Holographic Conductivity}

Our aim in this subsection is to calculate the holographic conductivity by
performing gauge/gravity duality for both linear Maxwell (infinite $\beta $)
and nonlinear cases. Finally, we depict the behaviour of conductivity for
linear and nonlinear electrodynamics cases and show that our results are
supported by some experimental observations.

\subsubsection{Linear Maxwell case ($\protect\beta$-infinity)}

We consider a four-dimensional planar metri
\begin{equation}
ds^{2}=-\frac{\mathcal{F}_{\infty }(u)r_{+}^{4z}}{l^{2z}u^{2z}}dt^{2}+{\frac
l^{2}du^{2}}{\mathcal{F}_{\infty }(u)u^{2}}}+\frac{r_{+}^{4}}{u^{2}
(dx_{1}^{2}+dx_{2}^{2}),  \label{plmet3}
\end{equation
where

\begin{equation}
\mathcal{F}_{\infty }(u)=1-\frac{mu^{z+2}}{r_{+}^{2z-4}}+\frac
q^{2}b^{2z-2}l^{2z}u^{2z+2}}{zr_{+}^{4z+4}},  \notag
\end{equation
which can be obtained by defining $u=r_{+}^{2}/r$ in the metric (\ref{met}).
The horizon locates at $u=r_{+}$, while the QFT lives at $u=0$ boundary. We
turn on perturbations $g_{tx_{1}}\left( u\right) e^{-i\omega t}$ and
A_{x_{1}}\left( u\right) e^{-i\omega t}$ and therefore receive one Einstein
and one electrodynamic additional equations of motion as {\small
\begin{equation}
A_{x_{1}}^{\prime \prime }\left( u\right) +\left[ \frac{\mathcal{F}_{\infty
}^{\prime }(u)}{\mathcal{F}_{\infty }(u)}+\frac{3(1-z)}{u}\right]
A_{x_{1}}^{\prime }\left( u\right) +\frac{l^{2z+2}u^{2z-2}}{r_{+}^{4z
\mathcal{F}_{\infty }(u)^{2}}\left[ \omega ^{2}A_{x_{1}}\left( u\right)
\frac{qb^{2z-2}u^{z}\mathcal{F}_{\infty }(u)}{l^{2}r_{+}^{2z}}\left(
2g_{tx_{1}}\left( u\right) +ug_{tx_{1}}^{\prime }\left( u\right) \right)
\right] =0,  \label{E10}
\end{equation
} an
\begin{equation}
2g_{tx_{1}}\left( u\right) +ug_{tx_{1}}^{\prime }\left( u\right)
=4qr_{+}^{2z-4}u^{2-z}A_{x_{1}}\left( u\right) ,  \label{E2}
\end{equation
where the prime denotes derivative with respect to $u$. Combining Eqs. (\re
{E10}) and (\ref{E2}), one can easily find a decoupled equation for
A_{x_{1}}$

\begin{equation}
A_{x_{1}}^{\prime \prime }\left( u\right) +\left[ \frac{\mathcal{F}_{\infty
}^{\prime }(u)}{\mathcal{F}_{\infty }(u)}+\frac{3(1-z)}{u}\right]
A_{x_{1}}^{\prime }\left( u\right) +\frac{l^{2z+2}u^{2z-2}}{r_{+}^{4z
\mathcal{F}_{\infty }(u)^{2}}\left[ \omega ^{2}-\frac{4q^{2}b^{2z-2}u^{2
\mathcal{F}_{\infty }(u)}{l^{2}r_{+}^{4}}\right] A_{x_{1}}\left( u\right) =0.
\label{E3}
\end{equation
Near the boundary $u=0$, Eq. (\ref{E3}) reduces approximately to
\begin{equation}
A_{x_{1}}^{\prime \prime }\left( u\right) +\frac{3(1-z)}{u}A_{x_{1}}^{\prime
}\left( u\right) +\cdots =0,  \label{nearbound}
\end{equation
which has the solutio
\begin{equation}
A_{x_{1}}\left( u\right) =A^{0}+A^{1}u^{3z-2}+\cdots ,  \label{nearboundsol}
\end{equation
where $A^{0}$ and $A^{1}$ are constants of integration. Following the
procedure presented in appendix \ref{app1}, we can calculate the
conductivity as

\begin{equation}
\sigma =\frac{(3z-2)r_{+}^{6z-4}A^{1}}{4\pi i\omega b^{2z-2}l^{z+1}A^{0}}
\end{equation}

\subsubsection{Nonlinear electrodynamics}

Here we follow the above procedure to calculate the holographic conductivity
in the case of nonlinear electrodynamics. We consider the four-dimensional
planar metri
\begin{equation}
ds^{2}=-\frac{\mathcal{F}(u)r_{+}^{4z}}{l^{2z}u^{2z}}dt^{2}+{\frac
l^{2}du^{2}}{\mathcal{F}(u)u^{2}}}+\frac{r_{+}^{4}}{u^{2}
(dx_{1}^{2}+dx_{2}^{2}),
\end{equation
where

\begin{eqnarray}
\mathcal{F}(u) &=&1-\frac{mu^{z+2}}{r_{+}^{2z+4}}-\frac{2\beta
^{2}l^{2}b^{2z-2}u^{2z-2}}{(4-z)r_{+}^{4z-4}}+\frac{2\beta
^{2}l^{2}b^{2z-2}u^{z+2}}{r_{+}^{2z+4}}\left( \frac{q^{2}l^{2z-2}}{\beta
^{2}L_{W}(\varrho _{u})}\right) ^{\left( 4-z\right) /4}  \notag \\
&&\times \left\{ \frac{L_{W}^{2}(\varrho _{u})}{4+z}\mathbf{F}\left( \frac
4+z}{4},\frac{8+z}{4},\frac{z-2}{4}L_{W}(\varrho _{u})\right) +\frac{1}{4-z
\mathbf{F}\left( \frac{z-4}{4},\frac{z}{4},\frac{z-2}{4}L_{W}(\varrho
_{u})\right) \right\} ,  \notag
\end{eqnarray
and $\varrho _{u}\equiv q^{2}l^{2z-2}u^{4}/r_{+}^{8}\beta ^{2}$. Turning on
perturbations $g_{tx_{1}}\left( u\right) e^{-i\omega t}$ and
A_{x_{1}}\left( u\right) e^{-i\omega t}$, we find below equations of motion

\begin{gather}
A_{x_{1}}^{\prime \prime }\left( u\right) +\left[ \frac{\mathcal{F}^{\prime
}(u)}{\mathcal{F}(u)}+\frac{3(1-z)}{u}+\frac{2L_{W}\left( \varrho
_{u}\right) }{u(1+L_{W}\left( \varrho _{u}\right) )}\right]
A_{x_{1}}^{\prime }\left( u\right)  \notag \\
+\frac{l^{2z+2}u^{2z-2}}{r_{+}^{4z}\mathcal{F}(u)^{2}}\left[ \omega
^{2}A_{x_{1}}\left( u\right) -\frac{\beta b^{2z-2}u^{z-2}\mathcal{F}(u)}
l^{z+1}r_{+}^{2z-4}}\sqrt{L_{W}\left( \varrho _{u}\right) }\left(
2g_{tx_{1}}\left( u\right) +ug_{tx_{1}}^{\prime }\left( u\right) \right)
\right] =0,  \label{E1}
\end{gather
an
\begin{equation}
2g_{tx_{1}}\left( u\right) +ug_{tx_{1}}^{\prime }\left( u\right)
=4qr_{+}^{2z-4}u^{2-z}A_{x_{1}}\left( u\right) ,  \label{E21}
\end{equation
The result of combining Eqs. (\ref{E1}) and (\ref{E21}) is decoupled
equation for $A_{x_{1}}
\begin{gather}
A_{x_{1}}^{\prime \prime }\left( u\right) +\left[ \frac{\mathcal{F}^{\prime
}(u)}{\mathcal{F}(u)}+\frac{3(1-z)}{u}+\frac{2L_{W}\left( \varrho
_{u}\right) }{u(1+L_{W}\left( \varrho _{u}\right) )}\right]
A_{x_{1}}^{\prime }\left( u\right)  \notag \\
+\frac{l^{2z+2}u^{2z-2}}{r_{+}^{4z}\mathcal{F}(u)^{2}}\left( \omega ^{2}
\frac{4q\beta b^{2z-2}\mathcal{F}(u)}{l^{z+1}}\sqrt{L_{W}\left( \varrho
_{u}\right) }\right) A_{x_{1}}\left( u\right) =0.  \label{E31}
\end{gather
Near the $u=0$ boundary, Eq. (\ref{E31}) reduces approximately to Eq. (\re
{nearbound}) with the solution (\ref{nearboundsol}). Performing the
procedure of appendix \ref{app1}, the holographic conductivity can be
obtained a
\begin{equation}
\sigma =\frac{(3z-2)r_{+}^{6z-4}A^{1}}{4\pi i\omega b^{2z-2}l^{z+1}A^{0}
\left. \exp \left[ -\frac{1}{2}\left( \frac{l^{z}\omega
A^{0}r_{+}^{2z-6}u^{3-z}}{\beta b^{2z-2}}\right) ^{2}\right] \right\vert
_{u=0}.  \label{holcond}
\end{equation
As one can see from (\ref{holcond}), the behaviour of conductivity is
different for different ranges of $z$

\begin{equation}
\sigma =\left\{
\begin{array}{ll}
\frac{(3z-2)r_{+}^{6z-4}A^{1}}{4\pi i\omega b^{2z-2}l^{z+1}A^{0}}, & \text
for }z<3 \\
\frac{(3z-2)r_{+}^{6z-4}A^{1}}{4\pi i\omega b^{2z-2}l^{z+1}A^{0}}\exp \left[
-\frac{1}{2}\left( \frac{l^{z}\omega A^{0}r_{+}^{2z-6}}{\beta b^{2z-2}
\right) ^{2}\right] , & \text{for }z=3 \\
0, & \text{for }z>
\end{array
\right.
\end{equation}

\subsubsection{Behavior of the conductivity and experimental results \labe
{num}}

\begin{figure}[t]
\centering
\subfigure[]{
   \label{fig11a}\includegraphics[width=.46\textwidth]{fig11a.eps}\qquad}}
\subfigure[]{
   \label{fig11b}\includegraphics[width=.46\textwidth]{fig11b.eps}}
\caption{The behaviors of real and imaginary parts of conductivity $\protec
\sigma $ versus $\protect\omega /T$ for $z=1$ and different values of $q$
with $l=b=r_{+}=1$.}
\label{fig11}
\end{figure}

\begin{figure*}[t]
\centering
\subfigure[]{
   \label{fig12a}\includegraphics[width=.46\textwidth]{fig12a.eps}\qquad}}
\subfigure[]{
   \label{fig12b}\includegraphics[width=.46\textwidth]{fig12b.eps}}
\caption{The behaviors of real and imaginary parts of conductivity $\protec
\sigma $ versus $\protect\omega /T$ for $z=1$ and different values of
\protect\beta $ with $l=b=r_{+}=1$.}
\label{fig12}
\end{figure*}

\begin{figure*}[t]
\centering
\subfigure[]{
   \label{fig13a}\includegraphics[width=.46\textwidth]{fig13a.eps}\qquad}}
\subfigure[]{
   \label{fig13b}\includegraphics[width=.46\textwidth]{fig13b.eps}}
\caption{The behaviors of real and imaginary parts of conductivity $\protec
\sigma $ versus $\protect\omega /T$ for $z=1.1$ and different values of $q$
with $l=b=r_{+}=1$.}
\label{fig13}
\end{figure*}

\begin{figure*}[t]
\centering
\subfigure[]{
   \label{fig14a}\includegraphics[width=.46\textwidth]{fig14a.eps}\qquad}}
\subfigure[]{
   \label{fig14b}\includegraphics[width=.46\textwidth]{fig14b.eps}}
\caption{The behaviors of real and imaginary parts of conductivity $\protec
\sigma $ versus $\protect\omega /T$ for $z=1.1$ and different values of
\protect\beta $ with $l=b=r_{+}=1$.}
\label{fig14}
\end{figure*}

In order to illustrate the behavior of the conductivity, we should solve the
decoupled differential equation for $A_{x_{1}}$. By definin
\begin{equation}
A_{x_{1}}\left( u\right) =\mathcal{F}(u)^{-i4\pi \omega /T}S(u),  \label{Ax}
\end{equation
near the horizon, where $S(u)=1+a(u-r_{+})+b(u-r_{+})^{2}+\cdots $, we
remove the oscillations for numerical stability. It is remarkable to note
that at the horizon where $S(u)=1$, $A_{x_{1}}\left( u\right) \propto
\mathcal{F}(u)^{\pm i4\pi \omega /T}$. However, we choose $\mathcal{F
(u)^{-i4\pi \omega /T}$ in (\ref{Ax}) in order to perform ingoing boundary
condition. The coefficients $a$, $b$,$\cdots $ can easily be found by
looking for Taylor series expansion of differential equation (\ref{E3}) at
the horizon. These coefficients are initial values necessary for solving Eq.
(\ref{E3}) numerically to obtain the conductivity.

Figs. \ref{fig11} and \ref{fig12} depict the real and imaginary parts of the
conductivity for $z=1$ in terms of $\omega /T$. It is notable to mention
that the real part of the conductivity is the dissipative part while the
imaginary part is the reactive one. Fig. \ref{fig11} shows that $\sigma
_{DC}\equiv \mathrm{Re}\left[ \sigma \left( \omega =0\right) \right] $ is
greater for smaller values of $q$. The same behavior can be seen for the
nonlinear parameter $\beta $ i.e. $\sigma _{DC}$ increases when $\beta $
decreases. The behaviors of the real and imaginary parts of the conductivity
are in excellent agreement with the experimental results reported in \cit
{exp1}. These experimental results show the increase in $\mathrm{Re}\left[
\sigma \right] $\ near the zero frequency due to impurities and ionic
lattice. The holographic conductivity computed by employing a gravity dual
model with linear Maxwell electrodynamics cannot produce this behavior \cit
{hart}. However, as one can see in Figs. \ref{fig11} and \ref{fig12}, for
suitable choices of parameters the nonlinearity of electrodynamics model can
result this behavior for real part of conductivity.

The behaviors of the real and imaginary parts of the conductivity in terms
of $\omega /T$ for $z=1.1$ are illustrated in Figs. \ref{fig13} and \re
{fig14}. In this case, $\sigma _{DC}$ changes with respect to $q$ and $\beta
$ in the similar manner that it does for $z=1$. However, in contrast to the
asymptotic AdS spaces in which the real part of the conductivity tends to a
constant in large frequencies, in this case ($z=1.1$) \textrm{Re}$\left[
\sigma \right] $ grows. Similar behavior for \textrm{Re}$\left[ \sigma
\right] $ has been reported recently for the optical conductivity of
single-layer graphene induced by mild oxygen plasma exposure \cite{exp2}. It
is remarkable to note that our numerical calculations show that there is no
significant difference between the behaviour of conductivity for $\beta =10$
and larger.

\section{Summary and Concluding Remarks}

In this paper, we constructed a new class of Einstein-dilaton-Lifshitz black
solutions in the presence of exponential nonlinear electrodynamics. Our
solutions respect the Lifshitz symmetry $t\rightarrow \lambda ^{z}t$ and
\vec{\mathbf{x}}\rightarrow \lambda \vec{\mathbf{x}}$ at $r$-infinity
boundary where $z(\geq 1)$ is the dynamical critical exponent. The
exponential electrodynamics behaves as Born-Infeld electrodynamics for large
values of the nonlinear parameter $\beta $ where they reduce to linear
Maxwell regime, as expected. It is worth mentioning that, while the
Born-Infeld nonlinear electrodynamics removes divergences in the electric
field and has finite value \textit{near the origin} where $r\rightarrow0$,
the exponential form of the nonlinear electromagnetics does not cancel the
divergency of the electric field exactly at $r = 0$, however, its
singularity is much weaker than Maxwell theory. This is more reasonable
compared to the Born-Infeld case, since \textit{near the origin} where
r\rightarrow0$, the electric field of a point-like charged particle should
be an increasing function. The behavior of the electric field of the
exponential nonlinear electrodynamics \textit{near} the origin in the
absence and in the presence of the dilaton field was explicitly shown in
table A of Refs. \cite{hendi2} and \cite{SheyKaz}, respectively. Besides, it
was argued that in applications of the AdS/CFT correspondence to
superconductivity, exponential nonlinear electrodynamics, makes crucial
effects on the condensation as well as the critical temperature of the
superconductor \cite{ZPCJ.NN}. It was also recently observed that, in the
holographic superconductor, the exponential nonlinear electrodynamics can
increase the critical values of the external magnetic field as the
temperature goes to zero \cite{SheAb}.

We considered topological black holes with zero ($k=0$), positive ($k=1$)
and negative ($k=-1$) horizon curvatures, but the reality of charge of
asymptotic Lifshitz supporting Maxwell matter field imposes that there is no
allowed solution with negative horizon curvature except for asymptotic AdS
case ($z=1$). Next, we calculated the conserved and thermodynamical
quantities namely mass, charge, electric potential, entropy and temperature.
We obtained the mass in terms of the extensive quantities entropy $s$ and
charge $Q$ i.e. $M\left( s,Q\right) $ and checked that the intensive
quantities temperature $T$ and electric potential $U$ calculated from it
match with those obtained from the geometry of the black hole. Thus, the
first law of thermodynamics is satisfied.

Subsequently, we turned out to study thermal stability of the obtained
solutions. We showed that black brane solutions ($k=0$) are always thermally
stable in canonical ensemble. In grand-canonical ensemble black branes with
horizon radius $r_{+}$ larger than $r_{+\min }$ are stable while there is a
maximum $z$ ($\beta $) that black branes are unstable for values greater
than it. We observed that black holes ($k=1$) encounter the Hawking-Page
phase transition between small and large black holes in both canonical and
grand-canonical ensembles. For $z\geq 2$, our black holes are always stable
in canonical ensemble. For other values of $z$ in canonical ensemble as well
as grand-canonical ensemble, there is a minimum for dynamical critical
exponent that black holes are stable for values greater than it. Moreover,
black holes are stable for nonlinear parameters greater than $\beta _{\min }$
in canonical ensemble while they are unstable for values greater than $\beta
_{\max }$ in grand-canonical ensemble.

Afterward, we performed the gauge/gravity duality to compute the ratio of
shear viscosity to entropy for a three-dimensional hydrodynamic system by
using the pole method. We obtained the well-known $1/4\pi $ result which
shows that even the non-AdS symmetry of hydrodynamic system cannot affect
the value of viscosity. Finally, we turned to study the behavior of the
holographic conductivity for two-dimensional systems. We examined the issue
for both linear Maxwell and nonlinear exponential electrodynamics. Our
investigations revealed that the effect of nonlinearity is vanishing the
conductivity for $z>3$. We depicted the behaviors of real and imaginary
parts of conductivity for asymptotic AdS ($z=1$) and Lifshitz cases and
pointed out some two-dimensional graphene systems taken under consideration
experimentally which have behaviors resemble our numerical results.

\acknowledgments{We thank from the Research Council of Shiraz
University. This work has been supported financially by Research
Institute for Astronomy \& Astrophysics of Maragha (RIAAM), Iran.}


\section{Introduction}

At large temperatures and densities, strongly-interacting matter is expected to change from a confined hadronic phase to a deconfined phase of quarks and gluons where chiral symmetry is restored. As the partition function of quantum chromodynamics (QCD) cannot be solved perturbatively in the region of the phase transition, we have to rely on other techniques to study the QCD phase diagram. Lattice QCD has successfully discovered the crossover nature of the transition \cite{Aoki:2006we,Borsanyi:2010bp} and established continuum-extrapolated results for the QCD equation of state, both at small baryochemical potential $\mu_{\rm B}$ \cite{Borsanyi:2013bia,Bazavov:2014pvz}. Due to the fermionic sign problem, standard lattice QCD methods become computationally ineffective in the regime of large densities. A couple of methods have been developed to extend the region of current lattice QCD calculations to finite baryochemical potential \cite{Fodor:2001au,deForcrand:2002ci,Endrodi:2011gv}, but quantitative conclusions cannot yet be drawn. A QCD critical point has been excluded up to values of $\mu_{\rm B}/T\lesssim1$. Another approach to the QCD phase diagram, which reproduces lattice results at vanishing baryochemical potential and can be extended to larger densities, comes from solving a coupled set of Dyson-Schwinger equations for the quark and gluon propagators \cite{Fischer:2014ata, Eichmann:2015kfa}. Here, a potential critical point is found to be located at $(T^c, \mu^c_{\rm q})=(115,168) \textrm{ MeV}$.

Ratios of cumulants of conserved quantities like the net-baryon number or net-electric charge are sensitive to a critical point  \cite{Asakawa:2000wh,Stephanov:1998dy,Stephanov:1999zu} signaling the singularity of thermodynamic quantities via their relation to susceptibilities \cite{Cheng:2008zh,Asakawa:2009aj,Gupta:2011wh,Friman:2011pf}. In this context higher-order cumulants are of special interest because they are more sensitive to the correlation length of fluctuations \cite{Stephanov:2008qz,Athanasiou:2010kw}. From universality arguments it has been demonstrated that the critical contributions to the kurtosis, in particular, may become negative approaching the critical point from the crossover side in heavy-ion collision experiments \cite{Stephanov:2011pb}.
This leads to the expectation of measuring a decreasing kurtosis lowering the beam energies, followed by a more complicated non-monotonic structure depending on the interplay of the location of the QCD critical point and freeze-out conditions for fluctuation measures.

Measurements of the net-proton and net-charge kurtosis and skewness have been reported by the STAR collaboration \cite{Aggarwal:2010wy,Adamczyk:2013dal}, where significant deviations from the hadron resonance gas and UrQMD calculations were found at lower beam energies.

To understand the experimental data, it is important to develop dynamical models which are able to describe nonequilibrium effects of the QCD phase transition. Even if thermalization times are small during the evolution of the system created in a heavy-ion collision and local equilibrium is thus achieved, near the critical point the thermalization time diverges with a certain power of the correlation length given by the dynamical universality class \cite{HALPERIN:1969zza}. This phenomenon is called critical slowing down and limits the divergence of fluctuations due to finite-time effects. A phenomenological approach has been applied in \cite{Berdnikov:1999ph} to understand the growth of the correlation length in an evolving system. It was found that the correlation length does not grow beyond $1.5-2$~fm, but memory effects let the system remain correlated for a longer period than expected in equilibrium. Expanding Fokker-Planck dynamics in terms of powers of the correlation length over the system size in the scaling regime, the importance of memory effects was underlined in \cite{Mukherjee:2015swa}. In a real-time evolution of non-Gaussian moments the magnitude and the sign of the critical contributions could be significantly different from the equilibrium expectations. In this model a simplified and homogeneous expansion was assumed and the back-reaction of the order parameter fluctuations on the surrounding matter was not taken into account.

In this work, we focus on the real-time evolution of the fluctuations of the order parameter for chiral symmetry, the sigma field, as obtained in a coupled dynamics. While the sigma field is propagated explicitly via a stochastic-relaxational equation, it interacts with a fermionic heat bath which expands fluid dynamically \cite{Nahrgang:2011mg}. It has been demonstrated in \cite{Herold201414,Nahrgang:2011mv,Herold:2013bi,Herold:2014zoa}, that this model is able to describe critical slowing down as well as spinodal decomposition within a dynamical setup.

We follow the evolution of the system over various hypersurfaces of constant energy density of the coupled system and compare the fluctuations in the sigma field to the fluctuations in net-proton numbers which are obtained from a Cooper-Frye particlization prescription. We furthermore give a comparison to a mean-field evolution in order to pin down the fluctuations stemming from the initial state versus dynamical fluctuations near the phase transition.

We begin with a description of nonequilibrium chiral fluid dynamics (N$\chi$FD) in Sec.~\ref{sec:model}, including the equations of motion and a brief description of the implementation of initial state and particlization. In Sec.~\ref{sec:results}, we investigate the dynamics of the net-proton and sigma field kurtosis for an evolution in the crossover regime near the critical point. We conclude and give a brief outlook in Sec.~\ref{sec:summary}.

\section{Nonequilibrium chiral fluid dynamics}
\label{sec:model}

We study the fluid dynamical evolution in a heavy-ion collision using a quark-meson model with dilaton field,
\begin{subequations}
\label{eq:Lagrangian}
\begin{align}
{\cal L}&=\overline{q}\left(i \gamma^\mu \partial_\mu-g_{\rm q} \sigma\right)q + \frac{1}{2}\left(\partial_\mu\sigma\right)^2
+ \frac{1}{2}\left(\partial_\mu\chi\right)^2\\\nonumber
        &\phantom{=} + {\cal L}_A- U_{\sigma}-U_{\chi}~, \\
\label{eq:LagrangianA}
 {\cal L}_A&=-\frac{1}{4}A_{\mu\nu}A^{\mu\nu}+\frac{1}{2}g_A^2\left(\frac{\chi}{\chi_0}\right)^2 A_\mu A^\mu~, \\
U_{\sigma}&=\frac{\lambda^2}{4}\left[\sigma^2-f_{\pi}^2\left(\frac{\chi}{\chi_0}\right)^2\right]^2-h\left(\frac{\chi}{\chi_0}\right)^2\sigma~, \\
U_{\chi}&=\frac{1}{4}B\left(\frac{\chi}{\chi_0}\right)^4\left[\ln\left(\frac{\chi}{\chi_0}\right)^4-1\right]~,
\end{align}
\end{subequations}
as introduced in \cite{Sasaki:2011sd}. In addition to the breaking and restoration of chiral symmetry it accounts for scale symmetry via the dilaton or glueball field $\chi$ that is identified with a gluon condensate. In the present version of this model we consider light quarks only, so $q=(u,d)$. The two condensates $\sigma$ and $\chi$ dynamically generate masses for the constituent quarks and gluons, thus allowing us to fix the coupling parameters $g_{\rm q}=3.37$ and $g_A=850 \textrm{ MeV}$ from the ground state nucleon and glueball masses. The additional parameters of the chiral sigma model are standard values: the pion decay constant of $f_\pi=93$~MeV, the pion mass $m_\pi=138$~MeV, the explicit symmetry breaking term $h=f_\pi m_\pi^2$ and the self-coupling $\lambda^2=\frac{m_\pi^2-m_\sigma^2}{2f_\pi^2}$. For more details, the reader is referred to \cite{Sasaki:2011sd,Herold:2014zoa}.

In mean-field approximation, the effective thermodynamic potential reads
\begin{equation}
\label{eq:effpot}
 V_{\rm eff}=\Omega_{q\bar q}+\Omega_{A}+U_{\sigma}+U_{\chi}+\Omega_0~,
\end{equation}
with the quark and gluon contributions
\begin{align}
\Omega_{\rm q\bar q}&=-2 N_f N_c T\int\frac{\mathrm d^3 p}{(2\pi)^3} \left\{\ln\left[1+\mathrm e^{-\frac{E_{\rm q}-\mu}{T}}\right]\right.\\\nonumber
                    &\phantom{-2 N_f N_c T\int\frac{\mathrm d^3 p}{(2\pi)^3}}+\left.\ln\left[1+\mathrm e^{-\frac{E_{\rm q}+\mu}{T}}\right]\right\}~, \\
\Omega_{A}&=2 (N_c^2-1) T\int\frac{\mathrm d^3 p}{(2\pi)^3} \left\{\ln\left[1-\mathrm e^{-\frac{E_A}{T}}\right]\right\}~,
\end{align}
which depend on temperature $T$ and quark chemical potential $\mu=\mu_{\rm B}/3$,
$\Omega_0$ in Eq.~(\ref{eq:effpot}) is an unimportant constant to set the total energy to zero in the ground state. The quasiparticle energies of constituent quarks and gluons are generated via their effective masses as $E_{\rm q}=\sqrt{p^2+m_{\rm q}^2}$ and $E_A=\sqrt{p^2+m_A^2}$.

The mean-field values of the condensate fields, $\langle\sigma\rangle$ and $\langle\chi\rangle$, are obtained by minimizing the effective thermodynamic potential $V_{\rm eff}$ via
\begin{equation}
\label{eq:equilibrium}
\left. \frac{\partial V_{\rm eff}}{\partial\sigma}\right |_{\sigma=\langle\sigma\rangle}=0~,~~\left. \frac{\partial V_{\rm eff}}{\partial\chi}\right |_{\chi=\langle\chi\rangle}=0~.
\end{equation}
In what we will call the mean-field evolution, the order parameter fields are set to their mean-field values neglecting fluctuations and the pressure is given by $p=-V_{\rm eff}$. The energy and quark number density of the system are thus evaluated as $e=T\partial p/\partial T +\mu n_{\rm q}-p$ and $n_{\rm q}=\partial p/\partial \mu$. This is equivalent to conventional deterministic fluid dynamical calculations using a chiral equation of state (EoS).

From the curvature of the effective potential at the equilibrium value, the mass of the sigma field and thus the inverse correlation length are obtained as
\begin{equation}
 \label{eq:corrl}
 m_\sigma^2=\frac{1}{\xi_{\rm eq}^2}=\frac{\partial^2 V_{\rm eff}}{\partial\sigma^2}\bigg|_{\sigma=\langle\sigma\rangle}\, .
\end{equation}

In order to study nonequilibrium effects, we follow our previous works \cite{Nahrgang:2011mg,Nahrgang:2011vn,Herold:2013bi,Herold:2014zoa}, and propagate both order parameters explicitly. For the chiral condensate we derive a stochastic relaxation equation from the two-particle irreducible effective action as
\begin{equation}
\label{eq:eomsigma}
 \partial_\mu\partial^\mu\sigma+\eta_{\sigma}\partial_t \sigma+\frac{\delta V_{\rm eff}}{\delta\sigma}=\xi~,
\end{equation}
which takes into account interactions with the surrounding quark heat bath via a dissipative term and a stochastic noise field $\xi$. In the simplest approximation the noise is Gaussian

\begin{equation}
\label{eq:dissfluctsigma}
 \langle\xi(t,\vec x)\xi(t',\vec x')\rangle_\xi=\delta(\vec x-\vec x')\delta(t-t')m_\sigma\eta_{\sigma}\coth\left(\frac{m_\sigma}{2T}\right)~,
\end{equation}
 and has a vanishing expectation value $\langle\xi(t,\vec x)\rangle=0$. Due to the discretizing of the space-time delta function the noise term will be dependent on the lattice spacing \cite{CassolSeewald:2007ru}. In order to avoid this numerical cut-off dependence, we coarse-grain the noise term over the spatial extension of the equilibrium estimate for the correlation length as given in Eq. (\ref{eq:corrl}).

The damping coefficient $\eta_{\sigma}$ depends on temperature and chemical potential,
\begin{equation}
\label{eq:dampingcoeff}
  \eta_{\sigma}=\frac{12 g^2}{\pi}\left[1-2n_{\rm F}\left(\frac{m_\sigma}{2}\right)\right]\frac{1}{m_\sigma^2}\left(\frac{m_\sigma^2}{4}-m_{\rm q}^2\right)^{3/2}~,
\end{equation}
 and vanishes near the critical point, where the mass of the sigma field becomes zero and the constituent quarks massive.
Below the phase transition we use a damping coefficient of $\eta_{\sigma}=2.2/\textrm{fm}$ as has been estimated for the $\sigma-\pi$ interaction in \cite{Greiner:1996dx}.

In the temperature regime of interest the dilaton field only fluctuates minimally around its equilibrium value, as the restoration of scale symmetry occurs at much higher temperatures only. We, therefore, propagate small fluctuations according to the classical Euler-Lagrange equation of motion
\begin{equation}
\label{eq:eomchi}
 \partial_\mu\partial^\mu\chi+\frac{\delta V_{\rm eff}}{\delta\chi}=0~.
\end{equation}

The relaxation times of the constituent quarks and gluons are assumed to be much smaller than the long-wave length sigma mode and can thus be treated in local thermal equilibrium with the pressure

\begin{equation}
\label{eq:pressure}
 p(T,\mu; \sigma, \chi) = -\Omega_{q\bar q}-\Omega_{A}~.
\end{equation}

In order to conserve the total energy and momentum of the coupled system, the divergence of the energy-momentum tensor of the quark-gluon fluid $T^{\mu\nu}$ equals a source term from the sigma and dilaton fields

\begin{align}
\label{eq:fluidT}
\partial_\mu T^{\mu\nu}&=-\partial_\mu\left(T_\sigma^{\mu\nu}+T_\chi^{\mu\nu}\right)~,\\
\label{eq:fluidN}
\partial_\mu N_{\rm q}^{\mu}&=0~.
\end{align}
Thus, the fluid dynamical fields become stochastic as the evolution of the sigma field follows a stochastic differential equation.

Recently, we have used this model to study the dynamical evolution through a first-order phase transition, where spinodal decomposition plays an important role \cite{Mishustin:1998eq,Sasaki:2007db,Randrup:2009gp,Randrup:2010ax,Steinheimer:2012gc}. We have demonstrated the formation of non-uniform structures in the energy and baryon density \cite{,Herold201414} and the dynamical enhancement of fluctuations in the medium \cite{Herold:2014zoa}. Such effects are especially interesting for upcoming experiments at FAIR \cite{Friman:2011zz} and NICA \cite{nica:whitepaper} which will make the region of high baryon densities in the phase diagram accessible.

\subsection{Initial state}

In this paper, we use event-by-event initial conditions, as opposed to previous publications, where we used an averaged initial state, usually a smooth sphere or ellipsoid. The initial energy and baryon density profiles are obtained from the UrQMD transport model \cite{Bass:1998ca,Bleicher:1999xi} run at a  center-of-mass energy per nucleon pair of $\sqrt{s_{\rm NN}}=19.7\textrm{ GeV/c}$. However, as the underlying EoS is different from that of the effective chiral model used here, we have to scale the resulting quantities such that we can investigate the region around the critical point during the fluid dynamical evolution.
With this set of initial conditions from UrQMD the fields are initialized at their local equilibrium values according to the temperature and baryochemical potential profiles.
All events are generated with zero impact parameter.

The UrQMD initial state has been used in recent hybrid model calculations at lower beam energies \cite{Auvinen:2013sba} looking at observables like elliptic and triangular flow. It has been noted that it becomes less reliable at lower energies where different space-time regions might not thermalize along a contour of proper time, but gradually during the evolution. Also the impact of the core-corona separation is more important. We apply the UrQMD initial state on the crossover side to the right of the critical point, where the system traverses the critical region, but does not extend to extreme baryonic densities. Our study is of exploratory nature to investigate the effect of initial state fluctuations versus dynamical fluctuations near the critical point.

The second- and the fourth-order moments of the event-by-event volume averaged sigma-field fluctuations at initial proper time $\tau_0$ is $\langle\Delta\sigma^2\rangle_{0}=4\textrm{ MeV}^2$ and  $\langle\Delta\sigma^4\rangle_{0}=64\textrm{ MeV}^4$ respectively. The initial kurtosis is thus given by $\kappa\sigma^2=4\textrm{ MeV}^2$.

\subsection{Hypersurfaces of constant energy density and particlization}

In order to follow the evolution of the initial and the dynamical fluctuations we average the sigma field over hypersurfaces of constant energy density and look at the event-by-event fluctuations. The energy density is given by the sum of the local energy density of the fluid and of the order parameter fields
\begin{align}
\label{eq:edens}
e=&\,\,e_{\rm fluid}+\frac{1}{2}\left(\frac{\partial\sigma}{\partial t}\right)^2+\frac{1}{2}\left(\nabla\sigma\right)^2+U_{\sigma} \\\nonumber &+\frac{1}{2}\left(\frac{\partial\chi}{\partial t}\right)^2+\frac{1}{2}\left(\nabla\chi\right)^2+U_{\chi}~.
\end{align}
In order to make a first qualitative connection to experimental observables we apply a Cooper-Frye particlization prescriptions \cite{Cooper:1974qi,Cooper:1974mv} to produce protons and antiprotons from the fluid dynamical fields, in particular from the energy density in Eq. (\ref{eq:edens}) by the help of the Cornelius hypersurface finder developed in \cite{Huovinen:2012is}. Besides (anti-)protons we produce all non-strange particles implemented in the UrQMD model, such that the fully integrated energy, momentum, net-charge and net-baryon number are conserved exactly in each event. In general, the effect of sigma field fluctuations should couple to particle production, in particular (anti-)protons and pions, via an interaction term like $g\bar{p}\sigma p$ \cite{Stephanov:1999zu, Stephanov:2008qz, Athanasiou:2010kw}. These contributions will be considered in future work. In this study we do not apply a subsequent hadronic cascade, but work in this direction is underway.

\section{Results}
\label{sec:results}

The goal of this work is to study the real-time evolution of the kurtosis in a system which follows a trajectory on the crossover side near the critical point as seen in Fig.~\ref{fig:traj}. We have calculated the respective temperatures and quark chemical potentials as volume averages over different hypersurfaces of constant energy density and then averaged these values over a set of events. Along the phase boundary we notice a bending which typically occurs in the crossover fluid dynamical trajectories and isentropes of chiral effective models in mean-field approximation \cite{Kahara:2008yg,Nakano:2009ps,Herold201414}. We compare mean-field, where $\sigma=\langle\sigma\rangle$, $\chi=\langle\chi\rangle$, from Eq.~(\ref{eq:equilibrium}), and nonequilibrium evolutions, i.e. with the order parameter fields evolved according to Eqs.~(\ref{eq:eomsigma}, \ref{eq:eomchi}), of the system to disentangle initial fluctuations and dynamically evolved critical fluctuations. The hypersurface- and event-averaged quantities $T$ and $\mu$ for the trajectory do not differ significantly between the mean-field and the nonequilibrium evolution.

\begin{figure}[t]
{
\centering
    \includegraphics[width=0.5\textwidth]{traj.pdf}
}
\caption[trajectory]{(Color online) Event-averaged trajectory near the critical point (black dot). Lines of constant energy density are drawn to indicate the position of the particlization procedures.}
\label{fig:traj}
\end{figure}

\begin{figure}[t]
{
\centering
    \includegraphics[width=0.5\textwidth]{netproton_noneq_eq.pdf}
}
\caption[kurtosis]{(Color online) Kurtosis of the net-proton number as function of freeze-out energy for mean-field and nonequilibrium evolution.}
\label{fig:kurt}
\end{figure}

\begin{figure}[t]
{
\centering
    \includegraphics[width=0.5\textwidth]{netproton_noneq_c4c2.pdf}
}
\caption[kurtosis]{(Color online) Kurtosis of the net-proton number as function of particlization energy for nonequilibrium evolution compared with the thermodynamic equilibrium values calculated from generalized susceptibilities $c_4/c_2$.}
\label{fig:c4c2}
\end{figure}

First we extract the net-proton number kurtosis, a quantity that is also studied in experiment. It is calculated as
\begin{equation}
\kappa\sigma^2=\frac{\langle\Delta N^4\rangle}{\langle\Delta N^2\rangle}-3\langle\Delta N^2\rangle~,
\end{equation}
where $\Delta N= N-\langle N\rangle$ is the event-wise fluctuation in the net-proton number and $\langle\dots\rangle$ denotes an average over events.
For each event we extract the multiplicity with the above described particlization procedure over several hypersurfaces of constant energy density. Then, in order to observe fluctuations, we apply a cut in the rapidity of $|y|<0.5$ and in transverse momentum $0.4 \textrm{ GeV}<p_T<0.8 \textrm{ GeV}$. This is similar to the published data from the STAR Collaboration \cite{Aggarwal:2010wy,Adamczyk:2013dal}. The results are shown in Fig.~\ref{fig:kurt}, with a comparison between the mean-field and a nonequilibrium evolution. We can immediately see that for the mean-field scenario the values drop only slightly beyond unity, while in nonequilibrium we obtain a clear minimum at $e=2.5e_0$ where the kurtosis reaches a value of about $0.3$. Although there are event-by-event initial state fluctuations, they only result in  a flat behavior of the kurtosis as function of the particlization energy density. In the case of nonequilibrium, where we explicitly propagate the order parameter and allow for energy and momentum exchange between the field and the quark-gluon fluid, fluctuations build up when the systems passes the crossover region, resulting in a dip in the net-proton kurtosis.
To compare this to the thermodynamic net-quark number kurtosis, we calculated the ratio of the generalized susceptibilities $c_4/c_2=\kappa\sigma^2$ along the trajectory in Fig.~\ref{fig:traj}. These susceptibilities are defined as \cite{Skokov:2010uh}
\begin{equation}
c_n=\frac{\partial^n (p/T^4)}{\partial(\mu_{\rm q}/T)^n}~.
\end{equation}
The result is shown in Fig.~\ref{fig:c4c2}, together with the nonequilibrium net-proton kurtosis. Note that in order to make a graphical comparison possible, we have scaled $c_4/c_2$ by a factor of $10^{-5}$. We see a sharp minimum with a negative quark-number kurtosis at $e=2.6e_0$, around the same point where we have the minimum in the net-proton kurtosis. We note two things: First, the equilibrium signal for criticality survives even in the dynamical environment representing the situation in a heavy-ion collision, though clearly less pronounced. Second, the resulting suppression of the net-proton kurtosis is spread out over a larger range of energy densities as a result of the inhomogenous medium and critical slowing down. It is important to remember that the values of $T$ and $\mu$ for the calculation of $c_4/c_2$ are averaged over the whole volume of the fireball. Therefore, a stronly negative kurtosis does not only occur at $2.6e_0$, but also around this value, whereas in smaller regions of space-time.

\begin{figure}[t]
{
\centering
    \includegraphics[width=0.5\textwidth]{netproton_noneq_sigma.pdf}
}
\caption[proton]{(Color online) Net-proton kurtosis and kurtosis of the volume-averaged sigma field in comparison. The sigma field kurtosis is shown in unit $\textrm{MeV}^2$.}
\label{fig:sigma}
\end{figure}

In Fig.~\ref{fig:sigma} we compare the net-proton kurtosis of the nonequilibrium evolution to the kurtosis of the sigma field. We extract the latter one by volume-averaging the sigma field on hypersurfaces of constant energy density and subsequently calculating the event-by-event fluctuations of the obtained values $\langle\sigma\rangle_V$. Here we see a similar course in the two curves with a minimum at nearly the same energy density of about $2.5e_0$.
In addition we observe a maximum of the sigma-field kurtosis at lower energy densities, which occurs below the equilibrium expectation for the phase transition and can be attributed to the importance of memory effects in a nonequilibrium evolution.

It is important to note that while the sigma field is not, the net-proton fluctuations are generally subject to global net-baryon number conservation \cite{Bleicher:2000ek,urqmdkurtosis,Bzdak:2012an}. At the presently investigated range of baryochemical potentials, however, baryon stopping should only have a negligible effect.

\section{Conclusions}
\label{sec:summary}

We have studied the net-proton kurtosis within the model of nonequilibrium chiral fluid dynamics, including a particlization procedure. This model captures the essential nonequilibrium dynamics of the order parameters at the QCD phase transition and critical point. We evaluated the kurtosis of both the sigma field and the net-proton number along an evolution near the critical point as a function of the energy density on a hypersurface. Here we compared mean-field fluid dynamical calculations to those with an explicit propagation of the order parameter fields, the chiral and gluon condensates. This takes into account the nonequilibrium evolution of the fluctuations via a stochastic relaxation equation. In the nonequilibrium case we found a minimum in both sigma and net-proton kurtosis at the same energy density. This minimum occurs around the same energy density as the minimum in the thermodynamic net-quark number kurtosis, implying that the suppression of the net-proton kurtosis is a remnant of the negative thermodynamic kurtosis. In comparison to that, a mean-field evolution without propagation of the order parameters shows a flat kurtosis as function of the energy density on the hypersurface.

The aim of future work will be to use a more realistic EoS, possibly including both hadronic and quark degrees of freedom. Models including quarks and hadrons have been studied in \cite{Dexheimer:2009hi,Steinheimer:2010ib,Turko:2013taa}, its parameters constrained by both lattice QCD data at small baryochemical potentials as well as neutron star properties at small temperatures. We will then study the kurtosis as a function of beam energy to compare results with the beam energy scan program at STAR. Finally, it is necessary to consider the evolution of produced fluctuations in a hadronic cascade to account for effects such as for example isospin randomization and charge diffusion processes \cite{Kitazawa:2011wh,Kitazawa:2012at,Nahrgang:2014fza,Sakaida:2014pya}.

\section*{Acknowledgements}

This work is funded by Suranaree University of Technology (SUT) and CHE-NRU (NV.10/2558) project. The authors thank Xiaofeng Luo and Volker Koch for fruitful discussions and Dirk Rischke for providing the SHASTA code that was used for the fluid dynamical simulation. The computing resources have been provided by the National e-Science Infrastructure Consortium of Thailand, the Center for Computer Services at SUT and the Frankfurt Center for Scientific Computing. MN acknowledges support from the U.S. Department of Energy under grant DE-FG02-05ER41367 and a fellowship within the Postdoc-Program of the German Academic Exchange Service (DAAD).

\section*{References}
\bibliographystyle{apsrev4-1}