$z = z_h$: $$\begin{aligned}
A_t(z_h) &=& f(z_h)=0,\\
f(0) &=& 1,\\
A_t(0) &=& \mu + \rho z^2 +....\end{aligned}$$ Here, $\mu$ is the quark chemical potential and $\rho$ is quark number density.

Following the procedure in Ref.[@Li:2011hp], we can calculate the pressure of the system through the entropy density $$\begin{aligned}
s = \frac{e^{3A(z_h)}}{4z^3_{h}}.\end{aligned}$$ For fixed values of the chemical potential, the pressure density can be calculated by the integral $$\begin{aligned}
p = \int [s dT+\rho d\mu]=\int[s (\frac{\partial T}{\partial z_h}d z_h+\frac{\partial T}{\partial \mu}d \mu)+\rho d\mu]\end{aligned}$$ and the energy density of the system can be derived: $$\epsilon = -p + sT + \mu \rho.$$

At zero chemical potential $\mu=0$, by fitting the lattice results of equation of state for $N_f = 2$ QCD[@Burger2015The], one can fix parameters $\alpha = 1.8, \beta = 0.4\ {\rm GeV}, \gamma = 0.42\ {\rm GeV}$ . With these parameters, the pseudo critical temperature for the crossover is around $T_0=217~ {\rm MeV}$ [@Fang:2015ytf].

In order to describe the system at finite chemical potential, we also need to fix $h(\Phi)$ in Eq.(\[sb\]) and (\[S\_G\_E\]), which describes the chemical potential dependence of the gluon-dynamical potential of the system. From the experience in the PNJL model [@Li:2018ygx], we can use the higher order 