text:
\section{\label{intro}Introduction}
Over the last years, there has been an increasing interest in the thermodynamic properties of degenerate electrons in the quantum mechanical regime.
Such information is vital for the description of highly compressed matter \cite{fletcher,kraus,regan},
including plasmas in laser fusion experiments \cite{lindl,hu,hurricane,nora,gomez,schmit}
and in compact stars and planet cores \cite{knudson,militzer,nettelmann}.
In addition, the widespread density functional theory (DFT) approach crucially depends on the availability
of accurate quantum Monte Carlo (QMC) data for the exchange correlation energy of the UEG, hitherto at zero temperature \cite{alder,fci}.
However, in recent years more and more applications with highly excited electrons have emerged, which require to go beyond ground state DFT.
Hence, there exists a high current need for an ab-initio thermodynamic description of the UEG at finite $T$.

The widely used path integral Monte Carlo (PIMC) method, e.g.\ \cite{cep}, is a powerful tool for the ab-initio simulation of both distinguishable
particles (often referred to as ``boltzmannons'', e.g.\ \cite{mil,clark}) and bosons and allows for quasi exact
results for up to $N\sim10^3$ particles at finite temperature \cite{bon,bon2}.
However, the application of PIMC to fermions is hampered by the notorious fermion sign problem (FSP), e.g.\ \cite{loh}, which might render
even small systems unfeasible for state of the art QMC methods and is known to be NP-hard for a given representation \cite{troyer}.
With increasing degeneracy effects, permutation cycles with opposite signs nearly cancel each other and the statistical uncertainty
grows exponentially. Hence, standard PIMC cannot provide the desired results without further improvement.
Brown \textit{et al.} \cite{brown} have presented the first finite temperature results for the UEG down to $r_s=1$ using
restricted PIMC (RPIMC) \cite{node},
a popular approach to extend PIMC to higher degeneracy, that is, lower temperature and higher density.
To avoid the FSP, this method requires explicit knowledge of the nodal surface of the density matrix, which is, in general, unknown
and one has to rely on approximations. The use of the ideal nodes for a nonideal system appears to be problematic,
as has been shown for the case of hydrogen \cite{hyd1,hyd2}. In addition, it has been shown analytically that RPIMC does not reproduce
the exact limit of the ideal Fermi gas ($r_s\to0$) \cite{vfil1,vfil2}. Therefore, the quality of the RPIMC data remains unclear.
Indeed, recent configuration PIMC (CPIMC) \cite{tim1,tim2} results for the highly degenerate UEG
by Schoof \textit{et al.} \cite{prl} have revealed a significant disagreement between the two methods at small $r_s$ and low temperature.
While the first application of a novel density matrix QMC (DMQMC) approach \cite{blunt} to the UEG for four particles reports
 excellent agreement with CPIMC \cite{malone},
additional simulations of larger systems are needed to resolve the discrepancy towards RPIMC.
For completeness, we mention that QMC results by Filinov \textit{et al.} \cite{vfil4} cannot be used as a benchmark
due to the different treatment of the homogeneous positive background and a different account of the long-range Coulomb interaction \cite{yakub1,yakub2} than the usual Ewald summation. In this situation an independent third first-principle method, capable to treat WDM parameters, would be highly desirable.

In this work we, therefore, investigate the applicability of the recently developed permutation blocking PIMC (PB-PIMC) approach \cite{dornheim} to the uniform electron gas.
The basic idea behind PB-PIMC is to combine antisymmetric imaginary time propagators  \cite{det1,det2,det3}, i.e., determinants, between all ``time slices'' with a higher order factorization of the density matrix \cite{ho1,ho2,ho3,ho4}.
This means that each particle is represented by a ``path'' consisting of $3 \times P$ coordinates (``beads''), where $P$ is the number of high-temperature factors (or propagators).
The application of determinants leads to a relieve of the FSP by an effective cancellation of positive and negative terms in the partition function,
which belong to permutation cycles of different parity in standard PIMC. However, since the blocking is most effective
if the thermal wavelength of a single propagator is of the same order as the mean interparticle distance,
it is crucial to employ a higher order factorization scheme which allows for sufficient accuracy with only a few time slices.

The details of our PB-PIMC scheme are described in section \ref{method}, after a brief introduction of the employed model \ref{model}.
In section \ref{sconvergence}, we present our simulation results starting with a detailed investigation of the convergence behavior with respect
to the factorization of the density matrix.
We proceed  by simulating $N=33$ spin-polarized electrons, which is a commonly used model system of the UEG, see section \ref{sdensity}.
Interestingly, our PB-PIMC approach allows us to obtain accurate results over the entire density range and, therefore, to make a comparison
with the pre-existing RPIMC and CPIMC results for the UEG.
Finally, in section \ref{td} we investigate the applicability of our method with respect to the temperature.
We find that PB-PIMC, in combination with CPIMC, allows for the simulation of the UEG over a broad parameter range, which includes
the physically most interesting regime of warm dense matter, cf.\ Fig.\ \ref{sketch}.

     \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{sketch.png}
      \caption{Density-temperature plain around the warm dense matter (WDM) regime. PB-PIMC significantly extends the range of applicability of standard PIMC (qualitatively shown by the red dashe line)
      towards lower temperature and higher density while CPIMC is applicable to the highly degenerate and weakly nonideal UEG \cite{prl}.
      RPIMC data \cite{brown} are available for $r_s \ge 1$.
      The orange area marks the conditions of WDM and inertial confinement fusion (ICF) \cite{hu}.
      }
      \label{sketch}
  \end{figure}

\section{\label{theory}Theory}
\subsection{\label{model}Model Hamiltonian}
The uniform electron gas, often referred to as ``Jellium'', is a model description of Coulomb interacting electrons with a neutralizing background
of positive charges which are uncorrelated and homogeneously distributed.
To describe an infinite system based on a finite number of particles, one implements periodic boundary conditions and
includes the interaction of the $N$ electrons in the main cell with all their images via Ewald summation.
Following the notation from \cite{fraser}, we express the Hamiltonian of the $N$ electron UEG (in atomic units) as
 \begin{eqnarray}
  \nonumber \hat{H} = -\frac{1}{2} \sum_{i=1}^N \nabla^2_i + \frac{1}{2} \sum_{i=1}^N\sum_{j\ne i}^N e^2 \Psi( \mathbf{r}_i, \mathbf{r}_j) + \frac{N e^2}{2}\xi \quad ,
 \end{eqnarray}
with $\xi$ being the Madelung constant and the periodic Ewald pair potential
\begin{eqnarray}
\nonumber \Psi(\mathbf{r}, \mathbf{s} ) &=& \frac{1}{V} \sum_{ \mathbf{G} \ne 0 } \frac{ e^{-\pi^2\mathbf{G}^2/\kappa^2} e^{2\pi i \mathbf{G}(\mathbf{r}-\mathbf{s})} }{ \pi\mathbf{G}^2}
 \\ \label{pair} &-& \frac{\pi}{\kappa^2 V} + \sum_\mathbf{R} \frac{ \textnormal{erfc}( \kappa | \mathbf{r}-\mathbf{s} + \mathbf{R} | ) }{ |\mathbf{r}-\mathbf{s}+\mathbf{R} | } \quad .
\end{eqnarray}
Here $\mathbf{R}=\mathbf{n}_1L$ and $\mathbf{G}=\mathbf{n}_2/L$ denote the real and reciprocal space lattice vectors, respectively, with the box length $L$ and volume $V=L^3$.
The specific choice of the Ewald parameter $\kappa$ does not influence the outcome of Eq.\ (\ref{pair}) and, therefore, can be used to optimize the
convergence. PB-PIMC requires explicit knowledge of all forces in the system, and the force between the electrons $i$ and $j$ can be obtained from
 \begin{eqnarray}
  \mathbf{F}_{ij} = - \nabla_i \Psi( \mathbf{r}_i, \mathbf{r}_j ) \quad . \label{iforce}
 \end{eqnarray}
The evaluation of Eq.\ (\ref{iforce}) is relatively straightforward and we find
\begin{eqnarray}
 \mathbf{F}_{ij} &=&  \frac{2}{V} \sum_{\mathbf{G}\ne 0}\left( \frac{ \mathbf{G} }{ \mathbf{G}^2 }  \textnormal{sin}\left[ 2\pi\mathbf{G}(\mathbf{r}_i - \mathbf{r}_j)\right]e^{- {\pi^2\mathbf{G}^2}/{\kappa^2}} \right)  \nonumber \\ &+&
 \sum_{\mathbf{R}} \frac{ \mathbf{r}_i - \mathbf{r}_j + \mathbf{R} }{ \alpha^3} \left( \textnormal{erfc}(\kappa\alpha) + \frac{2\kappa\alpha}{\sqrt{\pi}}e^{-\kappa^2\alpha^2} \right)
\quad , \nonumber
 \end{eqnarray}
 with the definition $\alpha = | \mathbf{r}_i - \mathbf{r}_j + \mathbf{R} |$.

 \subsection{\label{method}Simulation method}
 To calculate canonical expectation values with the PB-PIMC approach \cite{dornheim}, we write the partition function in coordinate representation
 as
 \begin{eqnarray}
\label{z} Z = \frac{1}{N!} \sum_{\sigma\in S_N} \textnormal{sgn}(\sigma) \int \textnormal{d}\mathbf{R}\ \bra{\mathbf{R}} e^{-\beta\hat{H}} \ket{\hat{\pi}_\sigma\mathbf{R}} \quad ,
\end{eqnarray}
 with $\mathbf{R}={\mathbf{r}_1,\dots,\mathbf{r}_N}$ containing the coordinates of all electrons, $\hat{\pi}_\sigma$ denoting the
 exchange operator which corresponds to a specific element $\sigma$ from the permutation group $S_N$ and $\beta=1/k_\textnormal{B}T$.
 For the next step, we make use of the usual group property of the density matrix in Eq.\ (\ref{z}) and arrive at an expression for $Z$ which requires
 the evaluation of $P$ density matrices at $P$ times higher temperature. However, instead of the primitive approximation $e^{-\epsilon\hat{H}}\approx e^{-\epsilon\hat{K}}e^{-\epsilon\hat{V}}$, with
 $\epsilon = \beta / P$ being the imaginary time step of a single propagator and the kinetic and potential contributions to the Hamiltonian $\hat{K}$ and $\hat{V}$, respectively, we
 use the fourth order factorization \cite{ho2,ho3}
 \begin{eqnarray}
\nonumber e^{-\epsilon\hat{H}} \approx & & e^{-v_1\epsilon\hat{W}_{a_1}} e^{-t_1\epsilon\hat{K}}  e^{-v_2\epsilon\hat{W}_{1-2a_1}} \\  & & \times e^{-t_1\epsilon\hat{K}} e^{-v_1\epsilon\hat{W}_{a_1}} e^{-2t_0\epsilon\hat{K}} \quad . \label{fop}
\end{eqnarray}
The $\hat{W}$ operators in Eq.\ (\ref{fop}) denote a modified potential, which combines $\hat{V}$ with double commutator terms of the form
 \begin{eqnarray}
 [[\hat{V},\hat{K}],\hat{V}] = \frac{\hbar^2}{m} \sum_{i=1}^N |\mathbf{F}_i|^2 \quad ,
\end{eqnarray}
and, therefore, requires the evaluation of all forces on each particle, $\mathbf{F}_i = -\nabla_i V(\mathbf{R})$.
Our final result for the partition function is given by
\begin{eqnarray}
\label{finalz} Z &=& \frac{1}{(N!)^{3P}} \int \textnormal{d}\mathbf{X} \prod_{\alpha=0}^{P-1} e^{-\epsilon\tilde V_\alpha}e^{-\epsilon^3u_0\frac{\hbar^2}{m}\tilde F_\alpha} \\
 & & \textnormal{det}(\rho_\alpha)\textnormal{det}(\rho_{\alpha A})\textnormal{det}(\rho_{\alpha B}) \quad , \nonumber
\end{eqnarray}
 with the definition of the potential and force terms
 \begin{eqnarray}
 \tilde V_\alpha &=& v_1 V(\mathbf{R}_\alpha) + v_2 V(\mathbf{R}_{\alpha A}) + v_1 V(\mathbf{R}_{\alpha B}) \quad ,\\
 \nonumber \tilde F_\alpha &=& \sum_{i=1}^N \left( a_1 |\mathbf{F}_{\alpha,i}|^2 + (1-2a_1) |\mathbf{F}_{\alpha A,i}|^2 + a_1 |\mathbf{F}_{\alpha B,i}|^2 \right)  \quad ,
\end{eqnarray}
 and the diffusion matrices
 \begin{eqnarray}
 \nonumber \rho_\alpha(i,j) = \lambda_{t_1\epsilon}^{-D} \sum_\mathbf{n} \textnormal{exp} \left( -\frac{\pi}{\lambda^2_{t_1\epsilon}} ( \mathbf{r}_{\alpha,j} - \mathbf{r}_{\alpha A,i} + \mathbf{n}L)^2 \right) \quad ,
\end{eqnarray}
with $D$ being the dimensionality, see e.g.\ \cite{det1}.
Eq.\ (\ref{finalz}) contains two free coefficients, $t_0$ and $a_1$, which can be used for optimization, cf.\ Fig.\ \ref{parameter},
and the integration is carried out over $3P$ sets of coordinates, $\textnormal{d}\mathbf{X} =\textnormal{d}\mathbf{R}_{0}\dots\textnormal{d}\mathbf{R}_{P-1}\textnormal{d}\mathbf{R}_{0A}\dots\textnormal{d}\mathbf{R}_{P-1A}\textnormal{d}\mathbf{R}_{0B}\dots\textnormal{d}\mathbf{R}_{P-1B}$.
Instead of explicitly sampling each permutation individually, as in standard PIMC, we combine configuration weights of both positive and negative sign
in the determinants, which leads to a cancellation of terms and, therefore, an effective blocking of permutations.
When the thermal wavelength of a single time slice, $\lambda_{t_1\epsilon} = \sqrt{2\pi\epsilon t_1\hbar^2/m}$, is comparable
to the mean interparticle distance, the effect of the blocking is most pronounced and the average sign in our simulations is significantly increased.
However, with an increasing number of propagators $P$, $\lambda_{t_1\epsilon}$ decreases and, eventually, the blocking will have no effect and the sign
converges towards the sign from standard PIMC. Hence, it is crucial to employ the high order factorization from Eq.\ (\ref{fop}), which allows
for reasonable accuracy even for only two or three propagators.
We simulate the canonical probability distribution defined by Eq.\ (\ref{finalz}) using the Metropolis algorithm \cite{metropolis}
and refer to \cite{dornheim} for a more detailed description of the PB-PIMC method.

\subsection{Energy estimator}
The consideration of periodicity in the diffusion matrices requires minor modifications in the energy estimator presented in \cite{dornheim},
which can be derived from the partition function via the familiar relation
\begin{eqnarray}
 E = - \frac{1}{Z} \frac{\partial Z}{\partial \beta} \quad . \label{tde}
\end{eqnarray}
Inserting the expression from Eq.\ (\ref{finalz}) into (\ref{tde}) and performing a lengthy but straightforward calculation leads to
\begin{eqnarray}
 \nonumber E &=& \frac{1}{P}\sum_{k=0}^{P-1}\left( \tilde V_k + 3\epsilon^2u_0\frac{\hbar^2}{m}\tilde F_k \right) + \frac{3DN}{2\epsilon}
 \\ &-& \nonumber \sum_{k=0}^{P-1}\sum_{\kappa=1}^N\sum_{\xi=1}^N\left( \frac{ \pi\eta^{k}_{\kappa\xi} }{ \epsilon P \lambda_{t_1\epsilon}^2 }
 + \frac{ \pi\eta^{kA}_{\kappa\xi} }{ \epsilon P \lambda_{t_1\epsilon}^2 } + \frac{ \pi \eta_{\kappa\xi}^{kB} }{ \epsilon P \lambda_{2t_0\epsilon}^2} \right)
 \quad ,
\end{eqnarray}
 with the definition
 \begin{eqnarray}
  \eta^k_{\kappa\xi} &=&  \frac{\left( \rho_k^{-1}\right)_{\kappa\xi} }{\lambda_{t_1\epsilon}^D } \nonumber
 \sum_\mathbf{n} \textnormal{exp}\left[ -\frac{\pi}{\lambda_{t_1\epsilon}^2}(\mathbf{r}_{k,\kappa}-\mathbf{r}_{kA,\xi}+L\mathbf{n})^2\right] \\
 & &  (\mathbf{r}_{k,\kappa}-\mathbf{r}_{kA,\xi}+L\mathbf{n})^2 \quad .
 \end{eqnarray}
 For completeness, we note that the total energy $E$ splits into the kinetic and potential contribution, $K$ and $V$, in precisely
 the same way as before \cite{dornheim}.

\section{Results\label{results}}

 \subsection{Convergence\label{sconvergence}}
\begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{param_potential.pdf}
      \caption{Influence of the relative interslice spacing $t_0$ for $N=4$, $r_s=4$ and $\theta=0.5$ on the convergence of the propagator.
      The exact result known from CPIMC (green line) is compared to the PB-PIMC results for $P=2$, $P=3$ and $P=4$ for
      the fixed free parameter $a_1=0.33$ over the entire $t_0$ range. The optimal value is located around $t_0=0.14$.
      }
      \label{parameter}\vspace*{-0.5cm}
  \end{figure}

    \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{convergence_potential.pdf}\\
         \vspace*{-0.8cm}
      \includegraphics[width=0.98\textwidth]{convergence_kinetic.pdf}
      \caption{Convergence of the potential (top) and kinetic (bottom) energy for $N=4$, $r_s=4$ and $\theta=0.5$ with $t_0=0.14$ and $a_1=0.33$.
      In the top panel, the potential energy $V$ is plotted versus the inverse number of propagators $P^{-1}\propto\epsilon$
      and the PB-PIMC results are compared to the exact value known from CPIMC. The bottom panel shows the same information for the kinetic energy $K$.
      }
      \label{convergence}\vspace*{-0.5cm}
  \end{figure}

  We begin the discussion of our simulation results by investigating the convergence of the energy with the number of imaginary time propagators $P$.
To enhance the performance, the free parameters from the propagator, $a_1$ and $t_0$, can be optimized.
  In Fig.\ \ref{parameter}, we choose $a_1=0.33$, which corresponds to equally weighted forces on all time slices, and plot the potential energy $V$, calculated with $P=2$, $P=3$ and $P=4$, versus $t_0$ over the entire possible range
  for a benchmark system of $N=4$ spin-polarized electrons with $\theta=0.5$ and $r_s=4$. To asses the accuracy, we compare these results with the exact energy known
  from CPIMC (green line). Evidently, the optimal choice for this free parameter is located around $t_0=0.14$, which is consistent with previous findings by Sakkos \textit{et al.} \cite{ho3}
  and the application of PB-PIMC to electrons in a quantum dot \cite{dornheim}.
  For completeness, we mention that the kinetic energy $K$
  exhibits the same behavior. Hence, we use the combination $a_1=0.33$ and $t_0=0.14$ for all
  presented simulations in this work. However, it should be noted that our method converges for all possible choices of the free parameters.
 In Fig.\ \ref{convergence}, we demonstrate the convergence of the energy with respect to the number of propagators for the same system as in Fig.\ \ref{parameter}.
 However, since $V$ and $K$ nearly cancel for this particular combination of $r_s$, $\theta$ and $N$, we investigate the convergence of
 both contributions separately.
 The top panel shows the potential energy versus the inverse number of propagators $P^{-1}\propto\epsilon$ and we compare the PB-PIMC results
 to the exact value (with the corresponding confidence interval) from CPIMC. We find that as few as two propagators allow for a
 relative accuracy $\Delta V/|V|\sim 10^{-4}$ and with $P=4$ the potential energy is converged within error bars.
In the bottom panel, we show the same information for the kinetic energy $K$. The variance of $K$ is one order of magnitude larger than that of $V$ and, for two propagators, we find the relative time step error $\Delta K/K\sim 10^{-3}$. With increasing $P$, the PB-PIMC results are fluctuating around the exact value, within error bars.

     \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{potential_difference.pdf}
      \caption{Accuracy of two and three propagators over a broad $r_s$ range for $N=4$ and $\theta=0.5$ with $t_0=0.14$ and $a_1=0.33$.
      We show the relative difference between the potential energy from PB-PIMC and CPIMC, $\Delta V/|V|$, for the optimal parameters from the fourth order
      propagator.
      }
      \label{pdiff}\vspace*{-0.2cm}
  \end{figure}

 Finally, we adress the $r_s$--dependence of the time step error by comparing PB-PIMC results for $V$ with $P=2$ (red crosses) and $P=3$ (blue squares)
 to the exact values from CPIMC.
 In Fig.\ \ref{pdiff}, the relative error of the potential energy $\Delta V/|V|$ is plotted versus $r_s$ for $N=4$ spin-polarized electrons at $\theta=0.5$.
 The increased errorbars for larger $r_s$ are a manifestation of the sign problem from CPIMC \cite{tim2}, while for the rest the statistical
 uncertainty from PB-PIMC predominates.
The time step error is smaller for three propagators over the entire $r_s$--range, as it is expected, and adopts a maximum
around $r_s=1$. This can be understood by recalling the source of the systematic error in PB-PIMC.
For $r_s\to0$, the UEG approaches an ideal system and the commutator error from $\hat{K}$ and $\hat{V}$ vanishes.
For $r_s\to\infty$, on the other hand, the particles are more separated and the system becomes more classical.
Therefore, the neglected commutator terms are most important at indermediate $r_s$, which is the case for the results in Fig.\ \ref{pdiff}.

We conclude that as few as two or three propagators provide sufficient accuracy to assess the discrepancy between CPIMC and RPIMC observed in previous studies \cite{prl}.
In particular, the selected benchmark temperature, $\theta = 0.5$, is even lower than for all other simulations to be presented in this work.
Hence, the observed time step error constitutes an upper bound for the accuracy of our results in the remainder of the paper.

 \subsection{Density parameter dependence\label{sdensity}}
  \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{N33_rs_sign.pdf}
      \caption{The average sign is plotted versus the density parameter $r_s$ for three different temperatures and $N=33$ spin-polarized electrons with $P=2$, $a_1=0.33$ and $t_0=0.14$.
      }
      \label{s1}\vspace*{-0.5cm}
  \end{figure}

    \begin{figure*}
      \centering
      \includegraphics[width=0.34\textwidth]{theta_delta_potential_theta1.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_potential_theta2.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_potential_theta4.pdf}
      \includegraphics[width=0.34\textwidth]{theta_delta_kinetic_theta1.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_kinetic_theta2.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_kinetic_theta4.pdf}
      \includegraphics[width=0.34\textwidth]{theta_delta_energy_theta1.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_energy_theta2.pdf}\hspace{-0.4cm}
      \includegraphics[width=0.34\textwidth]{theta_delta_energy_theta4.pdf}
      \caption{Comparison of PB-PIMC with CPIMC and RPIMC for $N=33$ spin-polarized electrons and three temperatures.
      In the top row, the relative deviation of the potential energy towards PB-PIMC with $P=2$, $t_0=0.14$ and $a_1=0.33$ is plotted
      versus $r_s$. The center and bottom rows display the same information for the kinetic and total energy, respectively.
      }
      \label{d1}
  \end{figure*}
 Among the most interesting questions regarding the implementation of PB-PIMC for the UEG is the range of applicability with respect to the density
 parameter $r_s$. To adress this issue, we simulate $N=33$ spin-polarized electrons, which corresponds to a closed momentum shell and is often used
 as a starting point for finite size corrections.
 In Fig.\ \ref{s1}, we show the average sign $S$ versus $r_s$ for three different temperatures over a broad density range.
 All curves exhibit a qualitatively similar behavior, that is, a smooth decrease of $S$ towards smaller $r_s$ until it saturates.
 At large $r_s$, the coupling induced particle separation mostly exceeds the extension of the single particle wavefunctions and quantum exchange effects do not play a dominant role.
 With decreasing $r_s$, the UEG approaches an ideal system and the particles begin to overlap, which leads to sign changes in the determinants.
 However, due the blocking, the average sign, instead of dropping exponentially, remains finite which implies that, for the three depicted temperatures, PB-PIMC is applicable over the entire density range. This is in stark
 contrast to standard PIMC, see e.g.\ supplement of \cite{brown}.
 Nervertheless, with decreasing temperature the sign drops and the FSP makes the simulations more involved, cf.\ section \ref{td}.\par
 In Fig.\ \ref{d1}, we compare the corresponding energies with RPIMC \cite{brown2} and CPIMC \cite{prl}, where they are available.
The top row displays the relative difference in the potential energy towards PB-PIMC with two propagators.
For $\theta=4$ and $\theta=2$, we find excellent agreement with CPIMC. For the lowest temperature, $\theta=1$,
the CPIMC values are systematically lower by $\Delta V/|V| \lesssim 10^{-3}$.
However, this discrepancy can be explained by the convergence behavior of the propagator, cf.\ Fig.\ \ref{pdiff}, since the potential (and kinetic)
energy is expected to converge from above towards the exact result.
To confirm this assumption, we also plot results for $P=3$ and $\theta=1$, visualized by the grey triangles.
Evidently, these points coincide with the CPIMC data everywhere within the errorbars and, thus, can be regarded as quasi-exact.
The RPIMC data for $V$, on the other hand, exhibit a systematic discrepancy with respect to PB-PIMC and CPIMC \cite{prl}.
At $r_s=1$, the energies approximately differ by $\Delta V/|V|\sim 0.02$, but the difference decreases with increasing $r_s$.
In the center row, we display the relative difference in the kinetic energy.
Again, all PB-PIMC results are in good agreement with CPIMC.
On the other hand, there is no clear systematic deviation between the PB-PIMC and RPIMC data, although most RPIMC-values for $\theta=1$ are lower while the opposite holds for most values for $\theta=4$.
Finally, the bottom row displays the relative difference in the total energy.
Interestingly, for $\theta=1$ the difference of RPIMC in $V$ and $K$ towards PB-PIMC nearly cancels, so that $E$ appears to be in good agreement.
In particular, even the value for $\theta=1$ and $r_s=4$, where the potential energy is an outlier, and both $V$ and $K$ exhibit a maximum deviation,
is almost within single error bars.
For completeness, we have also included the total energy for $\theta=1$ and $r_s=40$ from standard PIMC \cite{brown2}, cf.\ the black cirlce, which
is in excellent agreement with PB-PIMC as well.
For $\theta=2$ and $\theta=4$, most RPIMC values for $E$ are higher than PB-PIMC, although the deviation hardly exceeds twice the error bars.

 \subsection{Temperature dependece\label{td}}

   \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{N33_theta_sign.pdf}
      \caption{The average sign is plotted versus the temperature $\theta$ for $r_s=10$, $r_s=1$ and $r_s=0.1$ and $N=33$ spin-polarized electrons with $P=2$ and the free parameters $t_0=0.14$ and $a_1=0.33$.
      }
      \label{p1}\vspace*{-0.2cm}
  \end{figure}

 Finally, we investigate the performance of PB-PIMC with respect to the temperature.
 In Fig.\ \ref{p1}, the average sign is plotted versus $\theta$ for $N=33$ spin-polarized electrons at $r_s=10$, $r_s=1$ and $r_s=0.1$.
 All three curves exhibit a similar behavior, that is, a large sign $S$ at high temperature and a monotonous decay for
 $T\to0$. However, for $r_s=10$, the system is significantly less degenerate than for both other density parameters, and even
 at $\theta=0.5$, the average sign of $S\approx0.056$ indicates that the simulations are feasible.
 For $r_s=1$ and $r_s=0.1$, the decay of $S$ is more rapid and, at low temperature, the simulations are more involved.
 In particular, half the Fermi temperature seems to constitute the current limit down to which reasonable results can be achieved
 for such $r_s$--values (and this particle number) and, for $r_s=0.1$, the sign is zero within error bars, cf.\ the dashed line.
 Finally, we note that the average signs for the two smaller depicted $r_s$ parameters are more similar to each other than to $r_s=10$.
 We characterize the temperature in units of the ideal Fermi temperature, which is appropriate for weak coupling. However, for large $r_s$, the system becomes increasingly
 nonideal and, therefore, $\theta$ does not constitute an adequate measure for the degeneracy.

   \begin{figure}
      \centering
      \includegraphics[width=0.98\textwidth]{rs1_theta_potential_relative.pdf}\\
       \vspace*{-0.9cm}
      \includegraphics[width=0.98\textwidth]{rs1_theta_kinetic_relative.pdf}
      \caption{Comparison with CPIMC and RPIMC as a function of temperature. In the top panel, the relative deviation of the potential energy
      from the PB-PIMC result is plotted versus $\theta$ for $N=33$ spin-polarized electrons and $r_s=1$. The bottom panel displays the same information
      for the kinetic contribution.
      }
      \label{c1}\vspace*{-0.5cm}
  \end{figure}

 In Fig.\ \ref{c1}, we compare the energies of the $N=33$ electrons at $r_s=1$ from PB-PIMC both to RPIMC \cite{brown2} and CPIMC.
 The top panel displays the relative difference in the potential energy versus $\theta$.
 The CPIMC results for $V$ are in good agreement with PB-PIMC, while the RPIMC data are systematically higher, by about $2\%$.
 Interestingly, this behavior appears to be almost independent of the temperature.
 In the bottom panel, the same information is shown for the kinetic energy and, again, PB-PIMC agrees with CPIMC over the entire temperature range.
 The large statistical uncertainty at $\theta=0.5$ is a manifestation of the FSP in PB-PIMC, which prevents us from obtaining more precise kinetic
 energies with feasible computational effort. The RPIMC data for $K$ are
 slightly lower, at low temperature, which confirms the trend observed by Schoof \textit{et al.} \cite{prl}, and seems to converge towards the other methods for large $\theta$.

\section{Discussion}
 In summary, we have successfully extended the permutation blocking path integral Monte Carlo (PB-PIMC) method \cite{dornheim} to the uniform electron gas at finite temperature.
 We have started the discussion with a brief introduction of our simulation scheme, which combines a fourth-order factorization of the density matrix
 with the application of antisymmetric imaginary time propagators, i.e., determinants.
 This allows us to combine permutations, which appear as individual configurations with positive and negative sign in standard PIMC,
 into a single configuration weight. Therefore, the average sign in our simulations is significantly increased.

 To assert the quality of our numerical results, we have investigated the optimization of the free parameters of our propagator
 and demonstrated the convergence of both the potential and kinetic energy with respect to the number of imaginary time steps.
 We have found that even for the lowest considered temperature, $\theta=0.5$, as few as two propagators allow for a relative
 accuracy of $0.1\%$ and $0.01\%$ in the kinetic and potential energy, respectively. After this preparatory work,
 we have shown results for $N=33$ spin-polarized electrons, which is a commonly used model system as it is well suited to be a starting point for
 the extrapolation to the macroscopic limit (finite size corrections).
 Interestingly, PB-PIMC is feasible over the entire density range and, therefore, allows us to compare our results to both CPIMC and RPIMC data, where they are available.
 Our PB-PIMC data exhibit a very good agreement with CPIMC, for both the potential and kinetic energy,
 for all three investigated temperatures.
 %
 On the other hand, we observe deviations between PB-PIMC and RPIMC of up to $3\%$ in the potential energy, which
 decreases towards strong coupling. For the kinetic energy, we find no systematic trend although, for $\theta=1$, most of the RPIMC-values are smaller while, for $\theta=4$, most are larger than the PB-PIMC results. However, for both temperatures this deviation hardly exceeds twice the RPIMC errorbars.

 Finally, we have investigated the applicability of PB-PIMC to the $N=33$ spin-polarized electrons with respect to the temperature.
 With decreasing $\theta$, exchange effects lead to more negative determinants in the configuration weights and, therefore,
 a smaller average sign. For the physically most interesting density regime, $r_s\sim 1$, simulations are feasible above $\theta=0.5$ while
 for larger $r_s$ even lower temperatures are possible.
 A comparison of the energies for $r_s=1$ over the entire applicable temperature range has again revealed an excellent agreement with CPIMC. On the other hand, we observe a nearly $\theta$-independent relative deviation between PB-PIMC and RPIMC in the potential energy of approximately $2\%$, whereas differences in the kinetic energy are observed only towards low temperature.

 We conclude that our permutation blocking PIMC approach is capable to provide accurate results for the UEG over a broad parameter range. This approach is efficient above a minimum temperature of about $0.5 T_F$ and, thus, complements CPIMC.
 Even though PB-PIMC carries a small systematic error (which is controllable and depends only on the number of time slices),  we expect it to be useful for the development and test of other new techniques such as DMQMC \cite{blunt,malone} and other novel versions of fermionic PIMC, such as the approximate treatment of exchange cycles by DuBois \textit{et al.}\ \cite{dubois} or a variational approach to the RPIMC nodes, e.g.\ \cite{brownphd}.

 A natural follow-up of this work will be the extension of PB-PIMC to unpolarized systems which, in
 combination with CPIMC, should allow for a nearly complete description of the finite temperature UEG over the entire density range.
 In addition, we aim for the application or derivation of finite size corrections in order to extrapolate our results to the macroscopic limit \cite{fraser,drummond,lin}
 which could be followed by the construction of a new analytical fit formula for the UEG at finite temperature, e.g.\ \cite{brown3,karasiev}.

 Finally, since PB-PIMC allows for efficient simulations in the warm dense matter regime, applications to two-component
 plasmas, such as dense hydrogen \cite{bonitz,morales,proton}, are within reach.

\section*{Acknowledgements}
We acknowledge stimulating discussions with C. Hann (Kiel and Durham, North Carolina) and V.S. Filinov (Moscow). This work is supported by the Deutsche Forschungsgemeinschaft via SFB TR-24
project A9 and via project BO 1366/10 as well as by grant
SHP006 for CPU time at the Norddeutscher Verbund f\"ur
Hoch- und H\"ochstleistungsrechnen (HLRN).

\section*{References}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}

The purpose of this note, generalizing results of \cite{GM04,GLMZ05,GLM07}
for asymptotically constant-coefficient operators
and of Gesztesy--Makarov \cite{GM04} for
of periodic-coefficient Schr\"odinger operators,
is to give an explicit connection between
two objects related to the spectra of periodic-coefficient
operators on the line, namely, a Birman--Schwinger-type {\it characteristic
Fredholm determinant} $D(\lambda)$ introduced in \cite{JZ,BJZ} and the
Jost function-type {\it periodic Evans function} $E(\lambda)$
of Gardner \cite{G}.

Both of these objects are analytic functions whose zeros have
been shown in various contexts to agree in location and
multiplicity with the eigenvalues of the associated operator.
Thus, on the mutual domain for which these properties have been
established, it follows that they must agree up to a nonvanishing
analytic factor.
Here, we determine explicitly this nonvanishing factor, thus
illuminating the relation between the two functions while extending
the key property of agreement with eigenvalues to the union
of the domains on which it has been established for each
function separately.
Knowledge of this scaling factor is
useful also for numerical approximation of spectra as
in \cite{BJZ,BJNRZ1,BJNRZ2}, allowing convenient comparison
of different methods.

Our method of proof proceeds by comparison with
an interpolating first-order Birman--Schwinger determinant,
using a simple version of the Volterra integral computation of \cite{GM04}
together with the relation established in \cite{JZ} between
Birman--Schwinger determinants and Hill's method, a spectral Galerkin algorithm
for finite approximation of spectra.
We carry out the analysis here for the simplest interesting
case of a second-order elliptic operator; however, our arguments
extend readily
to the full class of operators discussed in \cite{JZ}.

Consider the eigenvalue problem for a general second-order
periodic-coefficient ordinary differential operator
$ L= B_0(x)^{-1}\big(\partial_x^2 + \partial_x A_1(x) + A_0(x)\big), $
written in form
\be\label{e:L}
\big(\partial_x^2 + \partial_x A_1(x) + A_0(x) -\lambda B_0(x)\big)U=0,
\ee
where $L$ is
defined on complex vector-valued functions $U\in L^2[0,X]$
with periodic boundary conditions,
$A_j, B_j\in L^2$ are matrix-valued and periodic on $x\in [0,X]$,
and $B_0$ is positive or negative definite in the sense that its
symmetric part $\Re B_0:=(1/2)(B_0+B_0^*)$ is positive or negative
definite.\footnote{The assumption of symmetry of $B_0$ made in \cite{JZ} is
not necessary for the arguments there.}
(Note that, by Lyapunov's Lemma,
definiteness can always be achieved by a variable--coefficient
change of coordinates,
provided that the eigenvalues of $B_0$ are of purely positive or purely
negative real part.)
Taking the Fourier transform, we may express \eqref{e:L}, equivalently,
as an infinite-dimensional matrix system
\be\label{infsys}
 (\D^2 +\D\A + \B)\U=0,
\ee
$ \B_{jk}= \hat A_{0,j-k} -\lambda \hat B_{0,j-k}$,
$\A_{jk}=\hat A_{1,j-k}$,
$\D_{jk}=(2\pi/X)\delta_j^k  ij$,
$\U_j= \hat U(j)$,
where $\hat f$ denotes the discrete Fourier transform of $f$.
(Here and elsewhere $i=\sqrt{-1}$.)
We shall alternate between these two representations as is convenient.

\section {The Birman--Schwinger Evans function}\label{s:review}
By a Birman--Schwinger-type procedure,
applying $(\partial_x^2-1)^{-1}$ to \eqref{e:L} on the left,
we obtain an equivalent problem
\be\label{equiv}
 (I+K( \lambda))U=0,
\ee
where $K=K_1+K_0$, with
$ K_1=\partial_x (\partial_x^2-1)^{-1} A_1,$
$K_0= (\partial_x^2-1)^{-1} (A_0+1 -\lambda B_0).$

\bl[\cite{JZ}]
For $A_j,B_j\in L^2$,
the operator $K$ is Hilbert-Schmidt.
\el

\begin{proof} (\cite{JZ})
Expressing $K$ in matrix form $\cK$ with respect to
the infinite-dimensional Fourier basis $\{e^{2\pi ijx/X}/\sqrt{X}\}$,
we find that
$
\cK_{1,jk}= \frac{(2\pi i/X)j}{1+(2\pi/X)^2j^2} \hat A_{1,j-k},
$
and similarly for $\cK_2$,
where $\hat A_j$ denotes Fourier transform
and $i:=\sqrt{-1}$.
Taking without loss of generality $X=2\pi$,
we find that
\be\label{calc}
\|K_1\|_{\cB_2}=\|\cK_1\|_{\cB_2}=
\sum_j
\frac{j^2}{(1+j^2)^2} \sum_k |A_1(j-k)|^2
=
\frac{j^2}{(1+j^2)^2} \|A_1\|_{L^2(x)}<+\infty,
\ee
and similarly for $\cK_2$.
\end{proof}

\begin{definition}\label{BSEvans}
\textup{
We define the {\it Birman--Schwinger Evans function} as
\be\label{e:evans2}
D(\lambda):={\det}_2(I-K( \lambda)),
\ee
where $\det_2$ denotes the $2$-modified Fredholm determinant,
defined for Hilbert--Schmidt perturbations of the identity;
see Appendix \ref{s:fredholm} for a review of the relevant theory.
}
\end{definition}

\bpr[\cite{JZ,BJZ}]\label{anlem}
For $A_j,B_j \in L^2$,
$D$ is analytic in $\lambda$; moreover, its zeros correspond in
location and multiplicity
with the eigenvalues of $L:=B_0^{-1}(\partial_x^2+\partial_x A_1+ A_0)$.
\epr

\begin{proof}
Analyticity follows by analyticity of the finite-dimensional Galerkin
approximations by which the Fredholm determinant is defined,
plus uniform convergence of the Galerkin approximations, a consequence
of \eqref{compare}.
Correspondence in location is immediate from the fact that
$\det_2(I-K)=0$ if and only if $(I-K)$ has a kernel, while
correspondence in multiplicity may be deduced by consideration
of a special sequence of Galerkin approximations on successive eigenspaces
of $L$, for which the Galerkin approximant can be seen to be a nonvanishing
multiple of the characteristic polynomial for the restriction of $L$
to these finite-dimensional invariant subspaces.
See \cite{JZ} for further details.
\end{proof}

\section{Hill's method}\label{hills}

Hill's method consists of truncating \eqref{infsys} at wave number
$J$, i.e., considering the $(2j+1)$-dimensional minor $|j|\le J$,
and solving the resulting finite-dimensional system to obtain
approximate eigenvalues for $L$: that is, the eigenvalues of
the $(2J+1)\times (2J+1)$ matrix
\be\label{LJ}
 L_J=(\B_{0,J})^{-1}(\D_J^2 +\D_J\A_J + \B_J) ,
\ee
$\B_{0,J,jk}=\hat B_{0,j-k}$,
where subscripts $J$ indicate truncation at wave number $J$, or, equivalently,
restriction to the $J$th centered minor as described above.

Left-multiplying by $(\cD_J^2-1)^{-1}$, similarly as above,
we may rewrite the truncated eigenvalue equation
 $(L_J-\lambda)U=(\D_J^2 +\D_J\A_J + \B_J)\cU=0 $ equivalently as
$ (I+ \cK_J)\cU=0, $
where $\cK_J=\cK_{1,J}+\cK_{2,J}$ is the truncation of the
Fourier representation $\cK=\cK_1+\cK_2$ of operator $K$, that is,
 $\cK_1=\D_J(\D_J^2-I)^{-1}\A_{1,J}$,
 $\cK_2=(\D_J^2-I)^{-1}(\A_{0,J}+1-\lambda \cB_{0,J})$.
We define the
{\it truncated Birman--Schwinger Evans function}, accordingly, as
\be\label{truncevans}
D_J(\lambda):={\det}_2(I-\cK_J)=
\det (\cD_J^2-1) \det \cB_{0,J} e^{\tr \cK_J(\lambda)} \det (L_J-\lambda).
\ee

\bpr[\cite{JZ}]\label{trunccorr}
Each $D_J$ is analytic in $\lambda$; moreover,
the zeros of $D_{J}$ correspond
in location and multiplicity with those of $L_{J}$.
\epr

\begin{proof}
Immediate,
from \eqref{truncevans}, the assumed uniform positive definiteness
of $B_0$, hence also of $\cB_{0}$ and $\cB_{0,J}$,
and properties of the characteristic polynomial.
\end{proof}

\bpr[\cite{JZ}] \label{convthm}
For $A_j, B_j\in L^2$, $D_{J}(\lambda)\to D(\lambda)$
as $J\to \infty$, uniformly on $|\lambda|\le R$.
\epr

\begin{proof}
By diagonality of $\cD$, $\cD_J\cA_J=
(\cD\cA)_J$,
hence $D_{J}$ is a sequence of Galerkin
approximations of the Fredholm determinant $D=\det_2(I-\cK)$
on a complete set of successively larger subapaces.
By definition of the Fredholm determinant, therefore,
$D_{J}\to D$ as $J\to \infty$.
\end{proof}

\bc
[Convergence of Hill's method \cite{JZ}] \label{Lconvthm}
For $A_j, B_j\in L^2$,
 the eigenvalues of $L_{J}$
approach the eigenvalues of $L$ in location and multiplicity
as $J\to \infty$, uniformly on $|\lambda|\le R$.
\ec

\begin{proof}
Immediate, from Lemma \ref{anlem}, Propositions \ref{trunccorr}
and \ref{convthm}, and properties of uniformly convergent analytic functions.
\end{proof}

\section{The Evans function of Gardner}

Rewriting \eqref{e:L} as a first-order system
\be\label{e:1L}
\partial_x W=A(\lambda)W,
\quad
W:=\bp U\\ U_x\ep,
\quad
\mA(\lambda):=\bp 0 & I\\
(\lambda B_0-A_0 -\partial_x A_1) & -A_1
\ep,
\ee
we may alternatively define an Evans function by
shooting, as done by Gardner.

\begin{definition}\label{def:GEvans}
\textup{
We define the {\it Jost-type Evans function} of Gardner, following \cite{G}, as
\be\label{e:GEvans}
E(\lambda):=
\det (\Psi(X)-I),
\ee
where $\Psi$ denotes the fundamental solution of \eqref{e:1L},
satisfying
$\Psi(0)=I$, and $X$ is the period.
}
\end{definition}

\bpr[\cite{G}]\label{gardlem}
For $A_j,B_j \in C^1$,
$D$ is analytic in $\lambda$; moreover, its zeros correspond in
location and multiplicity
with the (periodic) eigenvalues of $L:=B_0^{-1}(\partial_x^2+\partial_x A_1+ A_0)$.
\epr

\begin{proof}
Analyticity and agreement in location follow immediately
by analytic dependence of solutions of ODE and the fact that $\lambda$
is an eigenvalue precisely if there exists a solution
$W(x)=\Psi(x)W_0$ of \eqref{e:1L} for which $W(X)=\Psi(X)W_0= W_0$.
Agreement in multiplicity may be shown by a more detailed argument
based on choice of a Jordan basis as in \cite{G}.
\end{proof}

\section{Connections}

Define now the intermediate
{\it first-order Birman--Schwinger Evans function}
\be\label{int}
F(\lambda):={\det}_2
\big( (\partial_x-1)^{-1}(\partial_x-\mA)\big)
={\det}_2(I-\mK),
\ee
where $\mK:=(\partial_x-1)^{-1}(\mA - 1)$.

\bt\label{conn1}
For $A_j,B_j \in L^2$,
\be\label{FE}
F(\lambda)=
\gamma (e^X-1)^{-2n} E(\lambda),
\ee
where
\be\label{e:gamma}
\gamma:= e^{ \frac{e^X}{1-e^X}
\int_0^X(\tr(\mA)-2n)(y)dy}=e^{\frac{e^X}{e^X-1}(\tr(A_{1,ave})+2n)X},
\ee
and $A_{1,ave}$ denotes the mean over one period of $A_1$.
\et

\br\label{Xinf}
\textup{
From \eqref{FE}-\eqref{e:gamma}, we obtain
$
F(\lambda)\sim
E(\lambda) e^{ -\int_0^X \tr(\mA)(y)dy}
=
\det(I- \Psi(X)^{-1})
$
as $X\to \infty$, so that $F$ is asymptotic to a ``backward''
version of Gardner's Evans function.
}
\er

\begin{proof}
By a straightforward computation taking into account periodicity
together with the jump condition $G(y^+,y)-G(y_-,y)=1$ at $x=y$,
we find that the Green kernel associated with $(\partial_x-1)^{-1}$ is
semiseparable, of form
\be\label{G}
G(x,y)=\begin{cases}
\frac{e^{x-y}}{1-e^X},&x>y,\\
\frac{ e^{X+x-y}}{1-e^X},&x<y.\\
\end{cases}
\ee

Following \cite{GM04}, we may thus decompose $G$ as the sum
$G(x,y)=f(x)g(y)+J(x,y)$ of a separable kernel
\be
f(x)g(y):\quad
f(x)= e^x,
\quad
g(y)= \frac{e^X e^{-y}}{1-e^X},
\ee
and a Volterra kernel
\be\label{H}
J(x,y)=
\begin{cases}
e^{x-y},&x>y,\\
0,&x<y.\\
\end{cases}
\ee
Likewise, the kernel $\kappa(x,y)$ of $\mK$ decomposes as
\be\label{decomp}
\kappa(x,y)=f(x)g(y)(\mA(y)-1) +H(x,y),
\ee
where $H(x,y=J(x,y)(\mA(y)-1)$ is again a Volterra kernel,
vanishing for $x<y$.

Using the Fredholm determinant facts \eqref{mult}, \eqref{commute}
compiled in Appendix \ref{s:fredholm},
together with the fact that ${\det}_2(I-H)=1$ for a Volterra operator
$H$,\footnote{
Heuristically a strictly lower triangular perturbation of the identity;
$H$ has no nonzero eigenvalues.}
we find (here, with slight abuse of notation, denoting operators
by their associated kernels) that
\ba\label{reduct1}
{\det}_2 \big(I- &(\partial_x-1)^{-1}(\mA-1)\big)\\
& = {\det}_2 \big(I- H -fg(\mA-1))\\
&= {\det}_2 \big(I- H\big) {\det}_2 \big(I- (I- H)^{-1}fg(\mA-1)\big)
e^{-\tr\big(H(I-H)^{-1}fg(\mA-1) \big)}\\
&= {\det}_2 \big(I- (I- H)^{-1}fg(\mA-1)\big)
e^{\tr\big(fg(\mA-1) \big)}
e^{-\tr\big((I-H)^{-1}fg(\mA-1) \big)}\\
&= \gamma {\det}_{\CC^{2n}} \big(I- \langle g(\mA-1),(I- H)^{-1}f \rangle \big),
\ea
where $\langle \cdot, \cdot\rangle$ denotes $L^2[0,X]$ inner
product and
$\gamma:= e^{\tr\big(\langle g(\mA-1), f\rangle \big)}=
e^{\frac{e^X}{e^X-1}(\tr(A_{1,ave})+2n)X}$.

Computing
$$
(I-H)\Psi(x)= \Psi(x) - \int_0^x e^{x-y}(\mA-1)\Psi(y)dy,
$$
and recalling that $(\partial_x -1)\Psi=(\mA-1)\Psi$, where
$\Psi$ denotes the fundamental solution of $(\partial_x -\mA)\Psi=0$,
$\Psi(0)=I$,
so that, by Duhamel's principle/variation of constants,
$\Psi(x)=e^x + \int_0^x e^{x-y}(\mA-1)\Psi(y)dy$,
we find that
$(I-H)\Psi= e^x =f$, or $(I-H)^{-1}f=\Psi$.
Substituting into \eqref{reduct1}, we obtain,
noting that $\Theta:=e^{-y}\Psi(y)$ satisfies
$\partial_y \Theta= (\mA-1)\Theta$,
\ba\label{finalred}
{\det}_2 \big(I- (\partial_x-1)^{-1}(\mA-1)\big)&=
 \gamma{\det}_{\CC^{2n}} \Big(I- \Big(\frac{e^X}{1-e^X}\Big)
\int_0^X(\mA-1) e^{-y}\Psi(y)dy \Big)\\
&= \gamma {\det}_{\CC^{2n}}
\Big(I- \Big(\frac{e^X}{1-e^X}\Big) e^{-y}\Psi(y)|_0^X \Big)\\
&=\gamma {\det}_{\CC^{2n}}
\Big(I- \Big(\frac{e^X}{1-e^X}\Big) (e^{-X}\Psi(X)-I) \Big)\\
&= \gamma (e^X-1)^{-2n} {\det}_{\CC^{2n}} (\Psi(X)-I),
\ea
yielding the result.
\end{proof}

Now, define the
{\it truncated first-order Birman--Schwinger Evans function} as
\be\label{1truncevans}
F_J(\lambda):={\det}_2(I-\hat \cK_J),
\ee
where $\hat \cK_J$ denotes the $J$th Fourier truncation of the
Fourier transform $\hat \cK$ of $\mK$.
\bl\label{oneconv}
For $A_j, B_j\in L^2$, $F_{J}(\lambda)\to F(\lambda)$
as $J\to \infty$, uniformly on $|\lambda|\le R$.
\el

\begin{proof}
As in the argument of Proposition \ref{convthm},
this follows essentially by the definition of the Fredholm determinant.
\end{proof}

\bl\label{sim}
Similarly as in \eqref{truncevans}, we have
\be\label{eval1truncevans}
F_J(\lambda)
= \det(\cD_J-1)^2
\det \cB_{0,J} e^{\tr \hat \cK_J(\lambda)} \det (L_J-\lambda).
\ee
\el

\begin{proof}
Using $(I-\mK)=(\partial_x-1)^{-1}(\partial_x - \mA )$,
$\mA(\lambda)=\bp 0 & I\\ (\lambda B_0-A_0 -\partial_x A_1) & -A_1 \ep $
to obtain
$$
(I-\hat \cK_J)=\bp \cD_J-I_J & 0\\
0& \cD_J-I_J \ep^{-1}
\bp
\cD_J & -I_J\\
\cB_J +\widehat {(\partial_x A_1)}_J &\cD_J +\cA_J
\ep,
$$
and, using the block determinant formula
$\det \bp a & -I\\b&c\ep=
\det(ca+b)$,
$$
\begin{aligned}
\det \bp
\cD_J & -I_J\\
\cB_J +\widehat {(\partial_x A_1)}_J &\cD_J +\cA_J
\ep
&=
\det\big((\cD_J+\cA)\cD_J +\widehat {(\partial_x A_1)}_J+ \cB_J \big)\\
&=
\det\big(\cD_J^2+ \cD_J \cA_J + \cB_J \big)\\
&=
\det \cB_{0,J} \det (L_J-\lambda).
\end{aligned}
$$
\end{proof}

\bc\label{conn2}
For $A_j,B_j \in L^2$, $ F(\lambda)= \epsilon
e^{\hat \delta- \delta} D(\lambda) , $
where, for $X=2\pi$,
\ba\label{deltagamma}
\hat \delta&:= \tr (\hat A_{1,0}+2I)\Big(1+ \sum_{1\le |j|}
\frac{1}{j^2+ij}\Big),\\
\delta&:= -\tr(\hat A_{1,0}+I -\lambda \hat B_{0,0})\sum_{j}\frac{1}{j^2+1},\\
\epsilon&:=\lim_{J\to \infty}\prod_{|j|\le J} \Big(1+\frac{2}{ij+1}\Big).
\ea
\ec

\begin{proof}
We may readily calculate that
\ba\label{trK}
\tr \cK&= \tr(\hat A_{1,0})\sum_{|j|\le J}\frac{j}{-j^2-1}
+
\tr(\hat A_{1,0}+I -\lambda \hat B_{0,0})\sum_{|j|\le J}\frac{1}{-j^2-1}
\\
&=
-\tr(\hat A_{1,0}+I -\lambda \hat B_{0,0})\sum_{|j|\le J}\frac{1}{j^2+1}
\ea
and
\ba\label{trhatK}
\tr \hat \cK&=
\bp -(\cD_J-I_J)^{-1} & I\\
*
&
(\cD_J-I_J)^{-1}(-A_1-I) \ep \\
& = -\tr (\hat A_{1,0}+2I)\sum_{|j|\le J} \frac{1}{ij-1}\\
& = -\tr (\hat A_{1,0}+2I)\Big(-1+ \sum_{1\le |j|\le J}
\Big( \frac{1}{ij-1}-\frac{1}{ij}\Big)\Big)\\
& = \tr (\hat A_{1,0}+2I)\Big(1+ \sum_{1\le |j|\le J}
\frac{1}{j^2+ij}\Big),
\\
\ea
both uniformly convergent by absolute convergence of $\sum 1/j^2$.
Similarly,
$$
\prod_{|j|\le J}\frac{(ij-1)^2}{-j^2-1}=
\prod_{|j|\le J} \Big(1+\frac{2}{ij+1}\Big)
\sim e^{\sum_{|j|\le J} \frac{2}{ij+1} + O(j^{-2})}
\sim e^{\sum_{|j|\le J} O(j^{-2})}
$$
may be seen to be uniformly convergent.
Comparing
\eqref{truncevans} and \eqref{eval1truncevans}
and taking the limit as $J\to \infty$, we obtain the result.
\end{proof}

Collecting information, we have our final result, relating
$D$, $E$, and $F$.

\bc\label{finalconn}
For $A_j,B_j \in L^2$, and $X=2\pi$,
\be\label{finalrel}
D(\lambda)=
\frac{e^{\delta-\hat \delta}}{\epsilon}
F(\lambda) =
\frac{e^{\delta-\hat \delta}}{\epsilon} \gamma (e^X-1)^{-2n} E(\lambda),
\ee
where $\gamma$ is as in \eqref{e:gamma}
and $\delta$, $\hat \delta$, and $\epsilon$
are as in \eqref{deltagamma}.
\ec

\br\label{gardrmk}
\textup{
Since $\frac{e^{\delta-\hat \delta}}{\epsilon} \gamma (e^X-1)^{-2n}$
is analytic in $\lambda$ and nonvanishing,
Corollary \ref{finalconn} together with
Proposition \ref{anlem} gives an alternative
proof of Proposition \ref{gardlem}.
}
\er


text:
\section{Introduction}
\label{intro}

The recent cosmological observations have revealed that the present universe is accelerating in expansion \cite{Perl99,Riess98} due to the presence of an exotic component called
dark energy which constitutes more than seventy percent of the total cosmic density. Earlier to this discovery, there were indications
from weak \cite{RefregierA03} and strong\cite{TysonJA98} lensing phenomena, large scale structure \cite{AllenSW03}, galaxy rotation
\cite{ZwickyF93,RubinVC70,RubinVC80}, on the presence of another exotic component, the dark matter, which constitutes nearly
23$\%$ of the total density. The other known components are radiation and baryonic matter, comprising only about 1$\%$ and 5$\%$
respectively\cite{RefregierA03,TysonJA98,AllenSW03,ZwickyF93,RubinVC70}. Till date there is no
clear understanding regarding the nature of dark energy and dark matter. Many models have been proposed to explain the nature and evolution of these
two hitherto unknown components. In the most successful model of the universe, the standard $\Lambda$CDM model, the dark matter is
considered
as some non-relativistic matter of zero pressure and dark energy is considered as cosmological constant.
Due to the severe discrepancy between the observed and predicted values of the cosmological constant, called the cosmological constant problem, which doesn't
have any natural
explanation in the $\Lambda$CDM model, many alternative models have been considered in the recent literature with varying dark energy densities instead of
pure cosmological constant\cite{SolaJ13}. Regarding dark matter also, there were no general consensus at all about its nature.

Even though the current data supports the exotic components like dark matter and dark energy, it is not really ruling out the existence of other
exotic components. For instance, there exist speculations regarding the existence of a component called "dark radiation"\cite{Dutta11}. Another
exotic component of interest is stiff fluid. In some early papers,
Zel'dovich\cite{Zel'dovich62,Zel'dovich72} proposed a model in which the early stage of the universe is composed of stiff gas of baryonic particles
with equation of state $p_z=\rho_z,$ where $p_z$ is the pressure and $\rho_z$ is the density. This equation of state implies that the velocity of
sound in such fluids will become equal to that of light\cite{Zel'dovich62}. Later many others studied the nature and evolution of such a stiff fluid
in an expanding universe.
Steili et al\cite{Steili10} have shown that the self
interaction field between dark matter components will behave like
stiff fluid  and can be taken as an indication that the so called dark matter can be of stiff nature at least in the early
stages of the universe.
On investigating the early stages of
the universe, Barrow\cite{Barrow86} was able to establish the possible
existence of Zel'dovich fluid (or stiff fluid).
In dissipative cosmological models also, the presence of stiff fluid in the early stage of the universe has been speculated\cite{Cruz1}.
The equation of state of the Zel'dovich fluid directly implies that the energy density of it decreases as,
$\rho_z \propto a^{-6}$ as the universe expands, where $a$ is the scale factor of the universe.  Because of this, the normal
Zel'dovich fluid would have decoupled comparatively earlier to other components like
 radiation, which evolves as $a^{-4}$ and  matter evolving as $a^{-3}.$
 Consequently the presence of normal Zel'dovich fluid would have a dominant effect only during the early stage of the universe.
 Inspired by these, the effect of stiff
fluid on the primordial abundances of light elements was numerically
computed by Dutta and Scherrer\cite{Dutta10}. By comparing the prediction with observational abundance of light elements, especially with the
abundance primordial helium-4, these
authors have obtained a constraint on the Zel'dovich fluid density as,  $\rho_z / \rho_c < 30,$ where $\rho_c$ is the critical density of the universe.
Recently a simple cosmological model consisting of stiff fluid, non-relativistic matter and cosmological constant as dark energy was studied
by Chavanis\cite{Chavanis1} and he has shown that, the universe began with a big-bang and will end in a de Sitter phase.

Even though the conventional Zel'dovich fluid would not have any effect on the late stage of the universe, the Zel'dovich fluid with bulk
viscosity could have strong influence on the later stage. Effect of bulk viscosity, especially on the late evolution of the universe was
studied by many\cite{Brevik1,Athira15,Jerin1}. It is worthwhile to note the analayses of the bulk viscosity in the context of early inflation in
reference\cite{Bamba1}.
A cosmological model with bulk viscous Zel'dovich fluid as the dominant cosmic component was analysed in
reference\cite{TitusAswatiManoj14}.
This model predicts the late acceleration of the universe when the weighted bulk viscous coefficient is in the range $4 < \displaystyle 3\zeta /H_0 <6.$
The authors have further studied the evolution of the model in the statefinder, $(r,s)$ plane, and have shown that the current state of the universe
is arguably different from the standard $\Lambda$CDM model.
 The statefinder parameters are the most suitable diagnostics to contrast the model with other relevant dark energy models, expecially $\Lambda$CDM.
The evolution of the parameters
$r$ and $s$ and their present values are used to contrast the model. The $r-s$ plane is a two dimensional plane with $r$ and $s$ as the coordinates,
used to show the comprehensive evolution of (these parameters)
the model. For more details on statefinder analysis see \cite{Sahni1}.
The authors also proved that the model satisfies the generalized second law of
thermodynamics at the Hubble horizon. This work was further extended in reference\cite{TitusRajagopal16}, where the authors have used the Union
2 supernovae data to constrain
the possible value of the viscous coefficient and obtained $\zeta/H_0 \sim 5.25$ for $H_0 \sim 70$km/s/Mpc, where $H_0$ is the Hubble parameter of the present epoch
of the universe. However the age of
the universe was found to be around $10$ GY only. The over all behavior of this model is such that, it asymptotically evolves to a de
Sitter one with deceleration
parameter, $q \to -1.$
The authors\cite{TitusRajagopal16} have substantiated their results by performing a dynamical system analysis and have shown that the asymptotic
de Sitter epoch is stable. Another important point, evident from their work
is that the effective equation of state of the Zel'dovich fluid would start evolving from stiff nature ($p/\rho=1$), pass through a state of
pressureless matter ($p=0$) and eventually tend towards, $p=-\rho$ corresponding to de Sitter epoch. So the inclusion of the bulk viscosity in
the stiff matter would naturally take the fluid to manifest as pressureless matter and subsequently as a pure cosmological constant during the successive
stages of the evolution of the universe. In spite of this reasonable evolution of the universe in the bulk viscous Zel'dovich fluid, the odd man out is
the prediction regarding the age of the universe, which is less compared to the current observational results. In the present work we try to alleviate this
drawback of the model by incorporating one more component, a varying vacuum energy as dark energy.
We found that with such an addition in the cosmic components, the age of the universe can be enhanced to around $12$ GY.

A scalar field model of the universe with both stiff fluid and an effective interacting vacuum were proposed
by Cataldo et al.\cite{Cataldo1}. The authors have considered a scalar field in FLRW universe and shown that, the field
can effectively manifest as a mixture two barotropic fluids, one with an equation of state $p=\rho,$ corresponding to the
stiff fluid and the other with an equation of state $p=-\rho$ mimicking the vacuum. The self interaction of the field, which
manifests as the interaction between the components, would allow the cosmological vacuum to become a dynamical quantity. We will
further describes regarding this work in a later section.

In later  works, many have found that there exists an alternate route to the decaying vacuum, which comes from
quantum field theory techniques based on the renormalization group approach in curved spacetimes. The effective
action in such theories inherits quantum effects from the matter sector. In general, the
renormalization group techniques in curved spacetimes leads to a dependence of vacuum
energy on the Hubble parameter H of the form, $\rho_{\Lambda}(H,\dot H) = M_{Pl}^2 \Lambda(H,\dot H)$
\cite{Sola3,Sola4,Shapiro1} and references therin).

In the present work we found that the addition of the decaying vacuum energy will not affect the asymptotic
properties of the model, such that the end de Siter phase is still a stable one.
The paper is organized as follows. In section 2, we discuss the scalar field representation of a model with a mixture of
stiff fluid and interacting scalar field, which motivates the discussions in the following sections. In section 3, we present the
analytical solutions of the Hubble parameter, scale factor and the other cosmological
parameters. In section 4, we obtain the model paramters by constraining the model with supernovae observational data and then compute the evolution of the
different cosmological parameters. We also evaluate the age of the universe in this section. Section 5 is devoted to the study of the asymptotic properties
of this model followed by the conclusions in section 6.

\section{Scalar field approach to stiff fluid and interacting vacuum}

Let us consider the evolution of a self interacting scalar field, $\phi.$ which is minimally coupled to gravity in flat
isotropic and homogeneous universe. In describing the evolution of the scalar field we mainly follow the reference
\cite{Cataldo1}, so for more details see that reference. The evolution of the scalar field is governed by
the equations
\begin{equation}\label{eqn:phi1}
 3H^2 = 8\pi G \rho_{\phi}
\end{equation}
and
\begin{equation}\label{eqn:phi2}
 \dot \phi \ddot \phi + 3 H {\dot \phi}^2 = \dot \phi \frac{dV(\phi)}{d\phi}
\end{equation}
where $H$ is the Hubble parameter, $V(\phi)$ is the potential, dot represents a derivative with respect to cosmic time and
equation(\ref{eqn:phi2}) represents the dynamical evolution of the field. The equation of state of the field is,
\begin{equation}
 \omega_{\phi}=\frac{p_{\phi}}{\rho_{\phi}},
\end{equation}
in which the pressure and density of the scalar field are given by,
\begin{equation}
 \rho_{\phi}=\frac{{\dot \phi}^2}{2} + V(\phi), \, \, \, \, \, \, p_{\phi}=\frac{{\dot \phi}^2}{2} - V(\phi)
\end{equation}

The self interacting scalar field can be effectively treated as the mixture of two interacting fluids, with densities
$\rho_1$ and $\rho_2$ which are having the equation of state $\omega_1$ and $\omega_2$ respectively. Then the effective
pressure is,
\begin{equation}
 P_{eff}=\omega_{eff} \rho_{eff} = \frac{\omega_1 \rho_1 +\omega_2 \rho_2}{\rho_1+\rho_2} \rho_{eff},
\end{equation}
where $\rho_{eff}=\rho_1+\rho_2.$
Hence it is now possible to have two fluids, one with equation of state, $\omega_1=1,$ corresponding to stiff fluid and
the other with equation of state, $\omega_2=-1$ corresponding to vacuum energy, if one identifies the corresponding
densities and pressures as,
\begin{equation}\label{eqn:rhop1}
 \rho_1=\frac{{\dot \phi}^2}{2}, \, \, \, \, \, p_1=\frac{{\dot \phi}^2}{2}
\end{equation}
and
\begin{equation}\label{eqn:rhop2}
 \rho_2=V(\phi), \, \, \, \, \, p_2=-V(\phi).
\end{equation}
One can either take these components as isolated from each other, so that each one of them satisfying separate
conservation laws or can be taken as interacting components following the conservation equations,
\begin{equation}\label{eqn:con123}
 {\dot \rho}_1 + 3 H \left(\rho_1+p_1\right)= Q, \, \, {\dot \rho}_2 + 3 H \left(\rho_2+p_2\right)= -Q,
\end{equation}
for $Q>0$ the energy flows for $\rho_2$ to $\rho_1$ and for $Q<0$ the energy flow in the reverse direction. From equations
(\ref{eqn:rhop1}),(\ref{eqn:rhop2}) and (\ref{eqn:con123}), it follows that,
\begin{equation}
 {\dot \phi} \ddot \phi + 3 H {\dot \phi}^2 = Q(t), \, \, \, \, \dot \phi \frac{dV(\phi)}{d\phi}=-Q.
\end{equation}
These equations are the equivalent evolution equation of the scalar field. Hence it is possible to consider a scalar field
as a mixture of a stiff fluid interacting with an effective vacuum energy.

In reference \cite{Cataldo1}, the authors have described the evolution of such a universe. Our aim is not in line with that.
Instead we consider a phenomenological form for the decaying vacuum along with the bulk viscous stiff fluid (Zel'dovich fluid) and
analyze both the background evolution of the universe, particularly in finding the age and also the asymptotic behavior.

\section{Background evolution with viscous Zel'dovich fluid and decaying vacuum $\Lambda(t)$}

We consider a flat FLRW universe with standard metric
\begin{equation}\label{eqn:metric}
ds^2=-c^2dt^2+a^{2}(t)\left(dr^2+r^2 d\theta^2+r^2sin^2\theta \,
d\phi^2 \, \right)
\end{equation}
where $t$ is the cosmic time, $(r,\theta,\phi)$ are the co-moving
coordinates and $a(t)$ is the scale factor. The cosmic components are the
bulk viscous Zel'dovich and a varying cosmological constant. A non-viscous ideal Zel'dovich fluid
obeys the equation of state $p_z=\rho_z$\cite{Zel'dovich62}, which will change once the fluid is assumed to be viscous.
Bulk viscosity can be arised due to the deviations of the system from local thermodynamic equilibrium\cite{Zimdahl1}. It arises as an
effective pressure to restore the system back to its equilibrium, whenever the cosmic fluid expands or contracts too fast.
In incorporating the viscosity into the analysis we adopt
the Eckart's formulation\cite{Eckart40}.  Landau and Lifshitz also discussed a formulation\cite{Landau58} equivalent to Eckart's.
There had been works showing that the thermodynamical equillibria in Eckart's theory as unstable\cite{Hiscock85}  and signals could propagate through the
fluid at superluminal velocities\cite{Israel76}. But these difficulties can be overcome by taking account of the higher order terms. Such a more general
formalism was developed by Israel et al \cite{Israel791,Israel792}, to which Eckart's theory would appear as a first order limit. However owing
to the simplicity many authors favor Eckart's formalism for a first step analysis.
Later Hiscock and Salmonson \cite{Hiscock91} had shown that Eckart's formalism can be reasonably applied to  FLRW universe with late acceleration.

 For a first order deviation from thermodynamical equilibrium
the energy momentum tensor can be written as,
\begin{equation}
T^{\mu\nu}=\rho u^{\mu} u^{\nu} + \left(p+\Pi\right)h^{\mu\nu}
\end{equation}
where
 $\rho$ is the density of the fluid component, $u^{\mu}$ is the four velocity of an observer in Hubble flow,
 $\left(p+\Pi\right)$ is the effective pressure and  $h^{\mu \nu}=u^{\mu} u^{\nu}+g^{\mu \nu},$ with $g^{\mu\nu}$ as the metric coefficients.
 The above mentioned approach results in the equation for
 effective pressure as $p^{\prime}=p+\Pi$ for the bulk viscous Zel'dovich fluid and it can be expressed as,
\begin{equation}\label{eqn:pzdash}
p^{\prime}_z=p_z-3\zeta H
\end{equation}
where we have $p=p_z$ the normal pressure, $\Pi=-3\zeta H,$ the viscous pressure and $H$ is the Hubble parameter.
 In cosmology, bulk viscosity arises as an effective
pressure to restore the system back to its thermal equilibrium, which is broken when the cosmological fluid expands
(or contracts) too fast. This bulk viscosity pressure generated
ceases as soon as the fluid reaches the thermal equilibrium\cite{Wilson1}.

The second cosmic component in the present model is the time varying cosmological parameter given by\cite{Bessada1},
\begin{equation}
\Lambda(t)=3\alpha H^2
\end{equation}
where $\alpha$ is a free parameter, value of which would be less than one.
 Earlier
introduction of this kind of decaying vacuum was consdiered by Carvalho and Lima\cite{Carvalho1}, where the authors restricted to $\alpha \leq 1/2.$
A higher
value of $\alpha$ resulted into incompatible age for the universe as claimed by many authors like \cite{Bessada1}, so the values of $\alpha$ is usually
restricted to below one.
Since this is effectively a form of time varying
vacuum energy,
its equation of state is taken as, $ \displaystyle \omega_{\Lambda}=(p_{\Lambda}/\rho_{\Lambda})=-1.$ Basically the $\Lambda(t)$
models have been originated from curved space quantum field theories\cite{SolaJ11}. Often there appears a constant additive term along
with the time varying part in the equation for $\Lambda(t),$ which, as argued by many
\cite{SolaJ13,Paxy1} facilitate the transition from the decelerating to an accelerating epoch of the expanding universe. But in the present model such
an additive constant is not needed due to the presence of viscosity in the Zel'dovich fluid component, which will otherwise guarantee such a transition
from an early deceleration to a later accelerating phase of expansion.

The Friedmann metric along with the standard Einstein's field equation will give the Friedmann equation for a flat universe as,
\begin{equation}\label{eqn:FrRaEq1a}
3H^2(t)=\rho_z + \rho_{\Lambda}
\end{equation}
where $\rho_{\Lambda}=\Lambda(t)$ (in standard units, $8\pi G=1, \, c=1$), is the time varying cosmological parameter,
equivalent to the standard dark energy density. These components together satisfy the conservation law (in the absence of
any source, i.e. $Q=0$),
\begin{equation}
\dot{ \rho}_{\Lambda} + \dot{\rho}_z+3 H \left(\rho_z + p'_z \right)=0.
\end{equation}
where a dot represents a derivative with respect to cosmic time.
Combining the Friedmann equation and the conservation equation,  leads to
\begin{equation}\label{eqn:dotRho2}
\dot{\rho}_z+\dot{\rho}_{\Lambda}= - 6 H\left(\rho_z-\frac{3}{2}\zeta H\right).
\end{equation}
But from Friedmann equation, $\dot{\rho}_z+\dot{\rho}_{\Lambda} = 6 H \dot H.$ On substituting this, the above equation becomes,
\begin{equation}\label{eqn:dotH}
\dot{H}=- \left(\rho_z-\frac{3}{2}\zeta H\right).
\end{equation}
Again from Friedmann equation one can substitute for $\rho_z$ as,
\begin{equation}\label{eqn:rhoz2}
\rho_z=3\left(1-\alpha \right)H^2.
\end{equation}
From equations (\ref{eqn:dotH}) and (\ref{eqn:rhoz2}),
\begin{equation}\label{eqn:diffEqH}
\dot{H}+3H\left((1-\alpha)H-\frac{\zeta}{2}\right)=0.
\end{equation}
Solving equation (\ref{eqn:diffEqH}) we obtain the Hubble parameter as,
\begin{equation}\label{eqn:Hubble0}
H(t)=\eta \left[1+\coth\left(3(1-\alpha)\eta(t-t_0)+\phi \right) \right]
\end{equation}
where $\displaystyle \eta=\frac{\zeta}{4(1-\alpha)}$, $\displaystyle \phi=\coth^{-1}\left(\frac{H_0}{\eta}-1\right)$ and $H_0$ is the current
value of Hubble parameter. Integrating the above equation, we obtain the equation for the scale factor as
\begin{equation}\label{eqn:aoft0}
a(t)= e^{\eta(t-t_0)}\left(\frac{\sinh[3\eta(1-\alpha)(t-t_0)+\phi]}{\sinh(\phi)}\right)^{\frac{1}{3(1-\alpha)}}.
\end{equation}
Using this equation, the Hubble parameter in equation (\ref{eqn:Hubble0}) can be recast as,
\begin{equation}\label{eqn:Hubble1}
H=\frac{\zeta}{2(1-\alpha)}+\left(H_0-\frac{\zeta}{2(1-\alpha)}\right)a^{-3(1-\alpha)}.
\end{equation}
In the asymptotic limit $a(t) \to \infty $ the Hubble parameter becomes a constant,
$\displaystyle H \to \frac{\zeta}{2(1-\alpha)}$ which corresponds to the de Sitter phase with exponential increase in the scale factor,
while in the limit $a(t) \to 0,$ the Hubble parameter
evolves as, $H \sim a^{-3(1-\alpha)},$ which points to an earlier decelerated epoch dominated with Zel'dovich fluid with density
$\rho_z \sim H^2 \sim a^{-6(1-\alpha)}.$ Thus the existence of the transition from an early decelerated to a late
accelerated epoch is guaranteed. In these limits the scale factor will evolve as follows.
As $t \to \infty$ the scale factor will evolve as,
$a(t) \to e^{2\eta(t-t_0)},$ this exponential increase corresponds to the de Sitter epoch, while in the early stage of the evolution, corresponding to
$3\eta(1-\alpha)(t-t_0)<1$ the above form of
$a(t)$ almost implies that, $a(t) \sim \left((1+3\eta(1-\alpha)(t-t_0)\right)^{1/3(1-\alpha)},$ representing a decelerating epoch.
 What is important here is that the transition occurs
  without the aid of an additive constant in the $\Lambda(t).$ In non-viscous models like
entropic dark energy\cite{Basilokos1} or Ricci dark energy\cite{Paxy1}, the presence of a bare constant cosmological term is essential for having a
transition from the early decelerated epoch to the late accelerated epoch.

The evolution of the cosmological parameters, like Hubble parameter, scale factor etc are depending upon the numerical values of the model
parameters $\alpha$ and
viscous coefficient $\zeta.$ However it is clear from the expression of scale factor in equation(\ref{eqn:aoft0}) that for a constant $\alpha$
the beginning of the universe corresponding to
 $a=0$ would have occurred earlier into the past of the universe as $\zeta$ assumes higher values.
For constant $\zeta$ and increasing $\alpha,$ the situation will be the same too. In both the cases the age of the universe increases
compared to a model with only Zel'dovich fluid as the cosmic component. However, only with an extraction
of these parameters, a final conclusion regarding the age of the universe can be made.

  We have considered a flat universe ($k=0$), since observations
 strongly indicate that our universe is flat \cite{Komatsu1}.
 The inflationary models theoretically propose a very small value for curvature around $\Omega_{k0} \sim 10^{-5}$ while observations
 favor a value of the order of $10^{-2}.$
 Basically
 for non-flat universe, the Friedmann equation becomes,
 \begin{equation}\label{eqn:frw11}
  3H^2 = \rho_z + \rho_{\Lambda}+\rho_k.
 \end{equation}
 where $\rho_k=-ka^{-2}.$
Since the
interaction is only between Zel'dovizh fluid and the vacuum, the conservation law is,
\begin{equation}\label{eqn:cons11}
 \dot\rho_z + \dot{\rho}_{\Lambda} + 3 H \left(\rho_z+p^{\prime}_z\right)=0.
\end{equation}
Using equations (\ref{eqn:pzdash}) and (\ref{eqn:rhoz2}) and through simple algebra we can rewrite the above equation as,
\begin{equation}
 \dot H +3H\left( (1-\alpha)H-\frac{\zeta}{2}\right)=\dot{\rho}_k.
\end{equation}
In our original work we took, $\rho_k=0$ consequently the rhs of the above equation is zero. But for non-flat universe the contribution due to the
rhs term, $\dot{\rho}_k=-2H\rho_k$ is extremly small especially in the late stage, first of all due to the decreasing nature of $H$ and secondly
due to the extremely low magnitude of $\rho_k.$  Hence the solution of the above equation (our model) would be almost close to the solution of the
corresponding homogeneous equation with zero curvature.

\section{Extraction of model parameters and evolution of cosmic parameters }
The best fit values for $\zeta,\alpha$ and $H_0$  are estimated using type Ia supernova observational data. Union data
containing 307 data points \cite{Kovalsky08}
in the red shift range $0.01<z<1.55$ has been used here. For $i^{th}$  supernova  at a red shift $z,$ having
an apparent magnitude $m$ and absolute magnitude
$M,$ the distance modulus is,
\begin{figure}[h]
\centering
\includegraphics[scale=0.5]{Chi_Square_Alpha.pdf}
\caption{Contour plot for the parameters $H$ and $\bar\zeta$ for fixed $\alpha=0.14.$}
\label{fig:cont1}
\end{figure}
\begin{equation}\label{eqn:distancemodulus}
\mu_i (z)=(m-M) = 5\log_{10} d_L(z)+25,
\end{equation}
where $d_L(z)$ is its luminosity distance, which is depending on the model parameters also and is given by
\begin{equation}\label{eqn:dLz}
d_L(z,\bar{\zeta},\alpha)=\frac{c(1+z)}{H_0} \int_0^z \frac{dz^{'}}{h(z^{'},\bar{\zeta},\alpha)},
\end{equation}
where $\displaystyle h(z^{'},\bar\zeta,\alpha)=\frac{H(z^{'},\bar\zeta,\alpha)}{H_0}.$
Here we have redefined the viscosity coefficient by weighting it with the present Hubble parameter as, $\displaystyle \bar{\zeta}=\frac{\zeta}{H_0}.$
Equation for $h$ in terms of $z,$ the cosmological red shift is obtained by substituting for scale factor using $a=\frac{1}{1+z}$ in equation(\ref{eqn:Hubble1})
as
\begin{equation}\label{eqn:h}
h(z^{'},\bar\zeta,\alpha)=\frac{\bar{\zeta}}{2(1-\alpha)}+\left(1-\frac{\bar{\zeta}}{2(1-\alpha)}\right)(1+z)^{3(1-\alpha)}.
\end{equation}

The theoretical distance moduli for various red shifts are obtained using equation (\ref{eqn:dLz}) and are
compared with the corresponding observational data. The statistical $\chi^2$ function for comparing the theoretical and observational values of
the distance moduli is defined as,
\begin{equation}
\chi^2=\sum_{i=1}^n {\left[\mu_{ith}-\mu_{iob} \right]^2 \over \sigma_i^2}
\end{equation}
where $\mu_{ith}$ is the theoretical value of the distance modulus of the $i^{th}$ supernova for a given redshift and $\mu_{iob}$ is its observed distance
modulus corresponding to the same redshift.  $\sigma_i$ is the variance of the measurement of $i^{th}$ supernova and $n=307$ is the number of data points.
$\mu_{ith}$ being dependent on $(\alpha,\bar{\zeta},H_0),$ the best estimates of the parameters $(\bar\zeta,\alpha,H_0)$ are obtained by
minimizing the $\chi^2$ function. The minimum of the $\chi^2$ indicates the
goodness-of-fit of the model apart from giving the
best estimates of the model parameters.

For obtaining the $\chi^2$ function we also used
Background (CMB) data from the WMAP 7-yr observation and the Baryon Acoustic Oscillation
(BAO) data from Sloan Digital Sky Survey(SDSS). The BAO signal has been directly
detected by SDSS survey at a scale $\sim$100MPc. The BAO peak parameter value was first proposed by
Eisenstein, D. J. et al\cite{eisenstein} and is defined as
\begin{equation}\label{equ:17}
\mathcal{A}=\frac{\sqrt{\Omega_{m}}}{h(z_{1})^{\frac{1}{3}}}\left(\frac{1}{z_{1}}\int_{0}^{z_{1}}\frac{dz}{h(z)}\right)^{\frac{2}{3}}
\end{equation}
Here  h(z) is the Hubble parameter, $z_{1} = 0.35$ is the red shift of
the SDSS sample. Using SDSS data from luminous red
galaxies survey the value of the parameter $\mathcal{A}$(for flat
universe) is given by $\mathcal{A} = 0.469 \pm 0.017$\cite{eisenstein}. The
$\chi^{2}$ function for the BAO measurement takes the form
\begin{equation}\label{equ:18}
\chi^{2}_{BAO}=\frac{(\mathcal{A}-0.469)^{2}}{(0.017)^{2}}
\end{equation}
The CMB shift parameter is the first peak of CMB power
spectrum\cite{bond} which can be written as
\begin{equation}\label{equ:19}
\mathcal{R}=\sqrt{\Omega_{m}}\int_{0}^{z_{2}}\frac{dz}{h(z)}
\end{equation}
Here $z_{2}$ is the red shift at the last scattering surface. From
the WMAP 7-year data, $z_{2}=1091.3$. At this red shift
$z_{2}$, the value of shift parameter would be $\mathcal{R}=1.725\pm
0.018$\cite{Komatsu1}. The $\chi^{2}$ function for the CMB measurement can be
written as
\begin{equation}\label{equ:20}
\chi^{2}_{CMB}=\frac{(\mathcal{R}-1.725)^{2}}{(0.018)^{2}}
\end{equation}
Considering three cosmological data sets together, i.e. (SNe+BAO+CMB),  the total $\chi^{2}$ function is
then given by
\begin{equation}\label{equ:21}
\chi^{2}_{total}=\chi^{2}_{SNe}+\chi^{2}_{BAO}+\chi^{2}_{CMB}
\end{equation}

 By minimising the $\chi^2,$ we found parameter values as, $\alpha=0.14, H_0=70.3\textrm{km/s/Mpc}, \, \textrm{and} \,
 \bar\zeta=1.446$ for $\chi^2$ per degrees
of freedom, $\chi^2_{dof}= \frac{\chi^2_{min}}{n-m}=1.016,$ where $n$ is the number of data points and $m=3$ the number of free parameters.
We have constructed the confidence interval plane for the parameters $H$ and $\bar\zeta$ keeping $\alpha=0.14,$ its best estimated value. The confidence
intervals corresponding to 68.4$\%$ and 95.4$\%$ show fairly good behavior and are given in figure.\ref{fig:cont1}. The best
fit values for the parameters $H_0$ and $\bar \zeta$ with corrections for a confidence of 64.8$\%$ are
$H_0=70.03_{-0.46}^{+0.54}$ and $\bar\zeta=1.446_{-0.023}^{+0.018}. $ For 95.4$\%$ probability the corrected parameter values are
$H_0=70.3^{+0.565}_{-0.47}$ and $\bar\zeta=1,446^{+0.095}_{0.032}$ for $\alpha=0.14.$
It may be noted that in reference
\cite{Bessada1}, the authors have extracted an upper limit for the parameter $\alpha$ by constraining a model with a decaying
vacuum, $\Lambda=\Lambda_0+3\alpha H^2,$
as $\alpha \leq 0.15.$

In discussing the evolution of different cosmological parameters, it is better to start with the equation of state parameter. As it was shown in some
of the earlier works\cite{TitusAswatiManoj14}, the equation of state of the viscous Zel'dovich fluid has natural evolution from its extreme
stiff nature (corresponds $\omega=1$)
to the de Sitter type behavior through radiation (corresponds to $\omega=1/3$) like and matter like (corresponds to $\omega=0$) natures.
First we will consider the net equation of state, comprising both the Zel'dovich fluid and the decaying vacuum, which can be obtained by
the standard procedure as,
\begin{equation}\label{eqn:Omega1}
\omega(z)=-1-\frac{1}{3}\frac{d}{dx}(\ln {h^2})
\end{equation}
From the equations (\ref{eqn:h})and (\ref{eqn:Omega1}), the equation of state can be expressed as,
\begin{equation}\label{eqn:Omega2}
\omega(z)=-1+\left(2(1-\alpha)-\bar\zeta\right)\frac{1}{h}(1+z)^{3(1-\alpha)}.
\end{equation}
For  $\alpha=0$, and $\bar\zeta=0$ equation of state tends to $\omega(z) \to 1$ for very large redshift, which corresponds to the early epoch dominated with
Zel'dovich fluid with negligible viscosity.
\begin{figure}[h]
	\centering
	\includegraphics[scale=1]{Omega_alpha_fx.pdf}
	\caption{Evolution of $\omega$ with time for constant $\alpha=0.14$ and varied values of the viscous coefficient, $\bar\zeta.$}
	\label{fig:OAfx}
\end{figure}
For the chosen value of the parameter, $\alpha=0.14$ a transition into the late accelerating phase would occur for a range of values of the
viscous coefficient, $1.72>\bar\zeta >0.$   For the best fit values of the parameters, $(\alpha,\bar\zeta,H_0)=(0.14,1.445,70.03)$
 the evolution of the equation of state parameter is as expected (see figure \ref{fig:OAfx}), that is,
 starting with a value corresponding to stiff fluid and gradually approaching the value corresponding to de Sitter epoch.
The current value of $\omega$ corresponding to the best fit values of the parameters is found to be  $\omega \sim -0.72,$
very much close to the range of values of the equation of state from WMAP data\cite{Komatsu1}.

The equation of state of the viscous Zel'dovich fluid can be obtained using a similar approach as,
\begin{equation}
 \omega_z = -1 - \frac{1}{3} \frac{d\Omega_z}{dx}=-1 - \frac{(1-\alpha)}{3} \frac{d\ln h^2}{dx}.
\end{equation}
The difference in $\omega_z$ compared to $\omega(z)$ is that the second term on the right hand side of the above equation
contains an extra term, $(1-\alpha),$  so that the final expression for $\omega_z$ becomes,
\begin{equation}
 \omega_z=-1+(1-\alpha) \left(2(1-\alpha)-\bar\zeta\right)\frac{1}{h}(1+z)^{3(1-\alpha)}.
\end{equation}
The general evolution of $\omega_z$ is hence similar to that of $\omega(z)$ except in the particular numerical values corresponding to
different epochs. But both
will approach the de Sitter value as $a \to \infty.$

The deceleration parameter $q(z)$ is a measure of the acceleration and can be obtained from the basic equation,
\begin{equation}\label{eqn:q0}
q=-1-\frac{\dot H}{H^2}
\end{equation}
\begin{figure}[h]
\centering
\includegraphics[scale=0.7]{q_alpha_fx.pdf}
\caption{Evolution of deceleration parameter with $\alpha=0.14$ and varying $\bar\zeta.$ }
\end{figure}
is also evaluated. Using equation(\ref{eqn:Hubble1}),  deceleration parameter takes the form
\begin{equation}
q=-1+\left(\frac{3(1-\alpha)}{1+(\frac{1}{\frac{2(1-\alpha)}{\bar\zeta}-1})(1+z)^{-3(1-\alpha)}}\right)
\end{equation}
When both the model parameters $\alpha$ and $\bar\zeta$ are equal to zero, the cosmic component becomes pure Zel'dovich fluid and $q=2.$
For the best estimated values of the model parameters, the transition is found to occur at a redshift, $z\approx 0.61$
which is again close to the observed value. Irrespective of the values of $\bar\zeta$ the deceleration parameter asymptotically
approaches the value, $q=-1.$

Finally we will discuss about the evolution of the scale factor. The age of the universe can be directly obtained for the evolution of the scale factor.
The evolution of the scale factor is given in equation(\ref{eqn:aoft0}). We have already shown its asymptotic behavior in a previous section, that
in the early epoch it evolves as in the decelerated phase and in the extreme future it evolves as in de Sitter epoch. In general,
the form of $a(t)$ indicates
the presence of big-bang as $(t-t_0) \to -\infty.$
The evolution of it as shown in Figure (\ref{fig:a}) for the best fit values  of the model parameters indicate it. But for higher values of $\bar\zeta$ it
is found that the big-bang occur at earlier times as evident from figure (\ref{fig:aAlphaFx}).
\begin{figure}[h]
	\centering
	\includegraphics[scale=0.5]{aVStheta_2.pdf}
	\caption{Evolution of scale factor with time. The profile corresponds to best fit values, $\alpha=0.14$,
	$\bar\zeta=1.445$ and $H_0=70.03$ $km^{-1}s^{-1}Mpc^{-1}.$}
	\label{fig:a}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{a_alpha_fx.pdf}
\caption{The spectrum of curves on evolution of $a(t)$ with $\alpha=0.14$ and varying $\bar\zeta.$
There is big-bang for $\bar\zeta<1.72.$}
\label{fig:aAlphaFx}
\end{figure}
It is also clear that, for extremely high values of the viscous parameter, the scale factor would have no-zero values in the beginning
indicating the absence of big-bang. It is found that, there is no big-bang in the model for $\bar\zeta>1.72.$
This
means that age of the universe defined only for $\bar\zeta<1.72.$

The age of the universe in the present model can be obtained by equating the scale factor to zero, which leads to,
\begin{equation}
 Sh \left[3\eta(1-\alpha)(t_B-t_0)+\phi\right]=0.
\end{equation}
where we took, $t=t_B,$ as the big bang time.
This leads to the equation of the age of the universe as,
\begin{equation}
 t_0-t_B=\frac{\phi}{3\eta(1-\alpha)}.
\end{equation}
On substituting the expressions for $\phi$ and $\eta,$ the above equation can be simplified into,
\begin{equation}
 t_0-t_B=\frac{2}{3\bar\zeta} \ln \left(\frac{1}{1-2\bar\eta}\right) H_0^{-1},
\end{equation}
where $\bar\eta=\eta/H_0.$ For the best estimated values of the model parameters, it is found that, the above equation gives an age of the universe
in the range,
 $11.39-12.18$ GY. First of all, this age is higher than the age predicted by the model in which
 Zel'dovich fluid is the dominant component\cite{TitusRajagopal16}. Secondly it is near to the age of the universe obtained from the
 data on oldest globular clusters\cite{Pont1,CarrettaAge00, SalarisM97,SukyongYi01}. In this sense the model seems to solve the problem of age which
 existed in the model with Zel'dovich fluid as the only cosmic component. So it can be concluded that, the inclusion of a varying dark energy component
 along the bulk viscous Zel'dovich fluid is essential for the consistency with the age determination of the present universe.

 The analysis so far reveals that the inclusion of the additional component, the decaying vacuum to the viscous Zel'dovich fluid, the model gives
a reasonable back ground evolution of the universe. Apart from this, the age of the resulting universe will be high compared to the
model with Zel'dovich fluid
alone as the cosmic component. The viscous Zel'dovich fluid component evolves in such a way that, during the very early period the matter
component is a stiff
fluid, compatible with many theoretical speculations\cite{Dutta10}. But as the universe evolves, the equation of state is smoothly evolving
towards that of pressureless fluid,
corresponding to the non-relativistic matter. Hence the model supporting the speculation that in the early period the matter
would have existed as a stiff fluid. In
the late stage, the evolution is compatible with the standard $\Lambda$CDM model, in predicting the observational parameters,
including the age of the universe. In next section
we do a dynamical system analysis of the model, which may throw more light on the viability of the present model.

\section{Dynamical system analysis}
 Dynamical system analysis is an effective method to extract the useful information about the stability of the asymptotic behavior of the model.
 For this, we have to express the cosmological equations governing the evolution of the model as a set of autonomous differential equations.
 Then the concerned information can be obtained by finding
 \begin{figure}[h]
\centering
\includegraphics[scale=0.4]{PSHR.pdf}
\centering
\caption{Plot of vector field of the phase space around the critical point(0.607,0.457). The arrowhead of
trajectories are tilted towards the critical point.}
\label {fig:PSHR}
\end{figure}
 critical points and analyzing the nature of the trajectories in the neighborhood of the critical points. Eventually it will
 become clear from this analysis  whether the model in question is consistent with the presently accepted cosmological paradigm.
 In the present case this means
 whether the model predicts a stable evolution from a pre-decelerated epoch to a later accelerated epoch.

The first step in the dynamical system analysis is to identify proper phase-space variables. In the present case,
we define,
\begin{equation}
u=\Omega_z=(1-\alpha)h^2, \, \,  \, \, \, \,  v=\frac{1}{1+\frac{1}{h}}
\end{equation}
as the phase-space variables where the quantities have their usual meaning. These variables will take the range,
$0\leq u \leq 1$ and $0 \leq v \leq 1.$  The resulting coupled autonomous differential equations can then be formed using the
Friemdmann equations and they are,
\begin{equation}\label{eqn:2DEquationsOFmotion1}
\dot{u}=6H_0\frac{u}{1-\alpha}\left(\frac{\bar\zeta}{2}-\frac{(1-\alpha)v}{1-v} \right)
\end{equation}
and
\begin{equation}\label{eqn:2DEquationsOFmotion2}
\frac{\dot{v}}{\sqrt{1-\alpha}}=3H_0 v^2\left(\frac{\bar\zeta}{2\sqrt{u}}-\sqrt{1-\alpha}\right).
\end{equation}
\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{VHR_Vector_Field.pdf}
\caption{Plot of vector field of the phase space close to the critical point.}
\label{fig:vhrVF}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{VHR_Stream_Plot.pdf}
\caption{Stream plot of the trajectories around the critical point(0.607,0.457)}
\label{fig:vhrSP}
\end{figure}
The points in the phase space exhibit isomorphism with the exact solutions of cosmological field equations.
The ODEs are easier to be solved when the derivatives in them are written in terms of $\tau=ln(a).$ The behavior of ODEs in
the linear, closeby regions of critical points can be expressed in terms of a matrix equation accommodating $\displaystyle u'=\frac{du}{d \tau}$ and
$\displaystyle v'=\frac{dv}{d \tau}$  enabling flux analysis in terms of $\tau$ parameter of the equivalent autonomous ODE system.
We then establish the correspondence of the dynamics of $\rho_z=\rho_z(t),\rho_{\Lambda}=\rho_{\Lambda}(t)$ and $H=H(t)$ with the mentioned
simplified flux analysis in phase space for small perturbations in the linear limit $u\to u+\delta u(\tau)$ and $v\to v+\delta v(\tau),$ where
$\delta u$ and $\delta v$ are the perturbations.  The critical points are the solutions of the algebraic equations $P(u,v)=0$ and $Q(u,v)=0$ where
$\displaystyle P(u,v)=6H_0\frac{u}{1-\alpha}\left(\frac{\bar\zeta}{2}-\frac{(1-\alpha)v}{1-v}\right) $ and
$\displaystyle Q(u,v)=3H_0 v^2\left(\frac{\bar\zeta}{2\sqrt{u}}-\sqrt{1-\alpha}\right).$
The perturbations around the critical points satisfy the equation,
\begin{eqnarray}
\left(
\begin{array}{c}
  \delta u^{'} \\ \delta v^{'}
\end{array}
\right)
=
\left(
\begin{array}{cc}
\left(\frac{\partial P}{\partial u} \right)_0 & \left(\frac{\partial P}{\partial v} \right)_0  \\
\left( \frac{\partial Q}{\partial u} \right)_0 & \left(\frac{\partial Q}{\partial v}\right)_0
\end{array}
\right)
\left(
\begin{array}{c}
\delta u\\ \delta v
\end{array}
\right)
\end{eqnarray}
where the suffix,$'0'$ denotes the value at the critical point and $2\times2$ matrix in the above equation is the Jacobian.
The nature of the eigen values of the Jacobian matrix determine the behavior of the system near the critical points.

The critical point of interest corresponding to the equations (\ref{eqn:2DEquationsOFmotion1}) and (\ref{eqn:2DEquationsOFmotion2}) is
$(u_c,v_c)=(0.607,0.457).$ It can be seen that the critical value $v_c$ corresponds to the end de Sitter epoch, where the Hubble parameter
becomes, $h \to \frac{\bar{\zeta}}{2(1-\alpha)}.$  Then using the relation  $v=\frac{1}{1+\frac{1}{h}}$ and the best fit values of the model parameter,
it can be
seen that $v_c=0.457$ for the end de Sitter phase. The value $u_c=0.607$, the mass parameter of the bulk viscous Zel'dovich fluid also corresponds
to the end de Sitter epoch. That is, in the end de Sitter epoch, $u \to \frac{\bar\zeta^2}{4(1-\alpha)}$, which for the
best estimated values of the model parameters will become equal to 0.607 and is $u_c.$
The eigen values corresponding to this critical point is
 $(0,-4.11).$

The first eigenvalue, $0$ apparently suggests the absence of any isolated
critical point. But it is seen that such a situation does not arise originally from the ODEs by setting $\dot{u}=0$ and $\dot{v}=0.$
So the apparent discrepancy suggestive of lack of an isolated critical point arises from the errors of linear approximation of trajectory flux
of the system in its immediate neighborhood.
A low resolution view of vector field in phase-space is depicted in the figure (\ref{fig:PSHR}). Higher resolution
vector field plot as in figure (\ref{fig:vhrVF}) makes the view of the critical point as an attractor.
The high resolution stream plot as in figure(\ref{fig:vhrSP}) also is a clear indicator that the critical point is an attractor.

\section{Conclusion}

Many have speculated that matter present in the early stage of the universe were of stiff nature, with equation of state, $p/\rho=1.$
But owing to the fast decrease in its
energy density, it would have effect only on processes like primordial nucleosynthesis. Later works which studied the bulk viscous
stiff fluid, found that they can even
cause the late acceleration of the universe. But the main drawback of such models were that, they predicted less age for the universe.
In the present work we have studied
a model with Bulk viscous stiff fluid and decaying vacuum energy as cosmic components. We found that the model possesses
reasonably good back ground evolution, so as
to produce a late acceleration at about a redshift compatible with the observational results. The model also predicts a
de Sitter epoch as the end phase. It is found that
the acceleration is mainly due to the effect of bulk viscosity, because for the decaying vacuum to produce a transition
into the late acceleration, there has to be a constant
in the vacuum energy density. But the decaying vacuum we have considered doesn't have such a constant.
The effect of varying  vacuum energy is reflected in the age of the
universe. It was found that the age of the universe was increased compared to the model with viscous Zel'dovich fluid alone
as the component. Age obtained is in agreement
with the age deduced from the observations of the oldest globular clusters.

It was found that the equation of state of the fluid starts form the stiff nature, but eventually reduces to that of the matter and finally
goes over to that of a
pure cosmological constant corresponding to the de Sitter epoch. The dynamical system analysis shows that, the end de Sitter phase
is a stable one. During this stable end
de Sitter phase, the density of the Zel'dovich fluid is found to be around $0.6,$ which itself confirms that the late
phase of the universe in this model is dominantly
controlled by the viscous nature of the Zel'dovich fluid rather than the decaying vacuum energy.

\noindent{Acknowledgement}
One of the authors (KRN) thanks CUSAT Kochi, India for financial support in the form of a research fellowship. We also wish to
thank IUCAA, Pune, India for the local hospitality during our visits, where part of the work had been carried out.