\documentclass[11pt,dvips]{article} \textheight = 24truecm
\textwidth = 17truecm \hoffset = -2.3truecm \voffset = -2truecm
\renewcommand{\baselinestretch}{1.2}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{latexsym}
\usepackage{color}
\begin{document}
\thispagestyle{empty}
\begin{center}

\vspace{1.8cm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 {\bf xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx}\\
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\vspace{1.5cm}

{\bf xxxx}$^{a}${\footnote { email: {\sf xxxx}}} and {\bf
yyyy}$^{b,c}$ {\footnote { email: {\sf yyyy}}}

\vspace{0.5cm}
$^{a}${\it Department of Physics , Faculty of Sciences, University Ibnou Zohr,\\
 Agadir ,
Morocco}\\[1em]

$^{b}${\it LPHE-Modeling and Simulation, Faculty  of Sciences,
University
Mohammed V,\\ Rabat, Morocco}\\[1em]

$^{c}${\it Centre of Physics and Mathematics,
CPM, CNESTEN,\\ Rabat, Morocco}\\[1em]


\vspace{3cm} {\bf Abstract}
\end{center}
\baselineskip=18pt
\medskip

We study the evolution of geometric quantum discord (GMQD) of a two
qubits system coupled with two independent bosonic reservoirs. We
consider sub-ohmic, ohmic and super-ohmic. A special attention is
devoled to Dicke states and their superpositions.


\newpage
\section{Introduction}

  Entanglement has been attracting much attention from physicists in
both theory and experiment
\cite{11}\cite{22}\cite{33}\cite{44}\cite{55}, as it plays an
important role in quantum information processing such as quantum
communications \cite{22} and quantum cryptography \cite{66}.
However, there are some exceptions. For instance, quantum
computation \cite{77} based on one pure qubit does not employ
entanglement, but it needs other non-classical resources. This
suggests that entanglement is not the unique resource but just one
kind of quantum correlation. In fact, there is another kind of
quantum correlation named quantum discord.

Quantum discord, as a measure of bipartite non-classical
correlation, is a promising candidate and has generated a lot of
interest. For a quantum state $\rho $ in a composite Hilbert space
$H=H_{A}\otimes H_{B},$ the total amount of correlation can be
quantified by quantum mutual information \cite{88}:

\begin{equation}
I(\rho _{AB})=H(\rho _{A})+H(\rho _{B})-H(\rho _{AB})
\end{equation}

with $H(\rho )=-Tr[\rho \log _{2}\rho ]$ the von Neumann entropy and
$\rho _{A(B)}=Tr_{B(A)}\rho $ the reduced density matrix by tracing
system $B(A)$.

Quantum discord is a measure of nonclassical correlations that may
include entanglement but is an independent measure. We will document
with simple examples that the amounts of classical correlation,
quantum discord and entanglement bear no simple relationship to each
other. Taking system A as the apparatus, the quantum discord is
defined as follows:

\begin{equation}
D(\rho _{AB})=I(\rho _{AB})-C(\rho _{AB})
\end{equation}%
where $C(\rho _{AB})$ denotes the classical correlations of the state.%

The optimization procedure involved in the calculation of quantum
discord prevents one to write an analytical expression for quantum
discord even for simple two-qubit systems. Quantum discord is
analytically computed only for a few families of states including
the Bell-diagonal states \cite{99}\cite{111}, two-qubit X states
\cite{222}\cite{333}, two-qubit rank-2 states \cite{444}, a class of
rank-2 states of $4\otimes 2$ systems \cite{555}, and Gaussian
states of the continuous variable systems \cite{666}. Moreover,
based on the optimization of the conditional entropy, an algorithm
to calculate the quantum discord of the two-qubit states is
presented in . It is also important to have some computable bounds
on the quantum discord and some authors have obtained such bounds
\cite{777},\cite{888}.

This paper is organized as follows. In section 2, we give general
expression of GMQD for two qubits in Bosonic Reservoirs and discuss
the GMQD for a kind of X-state. In section 3, the Geometric quantum
discord. We study in section 4, the GMQD for Reservoirs with
Ohmic-Like spectral densities specially the GMQD for Dicke states
and their superpositions (generalized W states and GHZ states, the
superpositions of two Dicke states). In section 4, the conclusion.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Two spin  in independent bosonic reservoirs}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Bipartite states }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here, we
consider an ensemble of $N$ spin-$1/2$ particles described by the collective
angular momentum operators
\begin{equation}
J_{\alpha }=\frac{1}{2}\sum_{k=1}^{N}\sigma _{k\alpha },\text{ \ \ }\alpha
=x,y,z.  \label{collective}
\end{equation}%
Since the QD is adopted to
describe the nonlocal quantum correlation, it is necessary to derive
the pairwise density matrix. In this paper, we only study the
pairwise correlations.
%For the Dicke model and LMG model, the
%reduced pairwise matrix in the standard basis
For all collective spin models, including the Dicke model and the LMG model, the pairwise reduced density matrix
in the standard basis,
$\{|{\downarrow}{\downarrow}{\rangle},
|{\downarrow}{\uparrow}{\rangle},|{\uparrow}{\downarrow}{\rangle},
|{\uparrow}{\uparrow}{\rangle}\}$
(with ${\sigma_z}|{\uparrow}{\rangle}=|{%
\uparrow}{\rangle}$ and ${\sigma_z}|{\downarrow}{\rangle}=-|{\downarrow}{%
\rangle}$)~\cite{wang1},
can be derived as
\begin{equation}
~ \rho=\left(
\begin{array}{llll}
v_{+} & x_{+}^{*} & x_{+}^{*} & u^{*} \\
x_{+} & w & y & x_{-}^{*} \\
x_{+} & y & w & x_{-}^{*} \\
u & x_{-} & x_{-} & v_{-}
\end{array}
\right).  \label{rho:1}
\end{equation}
The detailed expressions for these elements are
\begin{eqnarray}
~
v_{\pm}&=&\frac{N^2-2N+4{\langle}J^{2}_{z}{\rangle}{\pm}4(N-1){\langle}
J_{z}{\rangle}}{4N(N-1)},  \label{rho:2} \\
x_{\pm}&=&\frac{(N-1){\langle}J_{+}{\rangle}{\pm}{\langle}[J_{+},J_{z}]_{+}{
\rangle}}{2N(N-1)},  \nonumber \\
w&=&\frac{N^2-4{\langle}J^2_{z}{\rangle}}{4N(N-1)},
y=\frac{{\langle}
J^2_x+J^2_y{\rangle}-N/2}{N(N-1)},  \nonumber \\
u&=&\frac{{\langle}J^2_+{\rangle}}{N(N-1)},  \nonumber
\end{eqnarray}
where $[A,B]_{+}=AB+BA$. $w=y$, for
$\sum_{\alpha=x,y,z}J^2_{\alpha}=J^2= \frac{N}{2}(\frac{N}{2}+1)$.
%For the symmetric states with parity conservation, we find
%$x_{\pm}=0$
In particular, since the Dicke and LMG models have symmetric ground states
with parity conservation, we find $x_\pm = 0$~\cite{wang1,julien1}.
%from Refs.~\cite{wang1,julien1}.
Hence the pairwise reduced density matrix in $X$ form is shown as
\begin{equation}
~ \rho=\left(
\begin{array}{llll}
v_{+} &0 & 0 & u^{*} \\
0 & w & y & 0 \\
0 & y & w & 0 \\
u & 0 & 0 & v_{-}
\end{array}
\right).  \label{rho:2}
\end{equation}
For the two-qubit states in $X$ form, the QD may be derived analytically, according
to Ref.~\cite{luo1}.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Coupling to bosonic reservoirs }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



The model under consideration is a quantum system of two qubits coupled to two independant bosonic reservoirs, the Hamiltonian of wich may be expressed as

\begin{equation}\label{z}
H=\sum^2_{j=1}[\frac{\nu_j}{2}\sigma_{j,3}+\sum_k\omega_{j,k}b_{j,k}^\dagger
b_{j,k}+\sum_k\sigma_{j,3}(g_{j,k}b_{j,k}^\dagger+g_{j,k}^*b_{j,k}],
\end{equation}

where $\nu_j$ is the energy difference between the excited state $\left\vert 1\right\rangle _{j}$
and the ground state $\left\vert 0\right\rangle _{j}$, and $\sigma_{j,3}$ is the Pauli matrix of qubit j with
$\sigma_{j,3}\left\vert 1\right\rangle _{j}=\left\vert 1\right\rangle _{j}$ and
$\sigma_{j,3}\left\vert 0\right\rangle _{j}=-\left\vert 0\right\rangle _{j}$. $b_{j,k}^\dagger(b_{j,k})$
and $\omega_{j,k}$ denote the bosonic creation (annihilation) operator and the frequency of the kth mode of the reservoir of the qubit j,respectively.$g_{j,k}$
denotes the coupling strength between the qubit j and the kth mode. The Hamiltonian describes the spin-boson model without tunneling. A possible experiment
setup is a double-dot charge qubit placed in a freestanding semiconductor slab.

Suppose that the two qubits are initially writes in the X states,

\begin{equation}\label{v}
\rho_s(0)=\frac{1}{4}(I\otimes I+x_3 \sigma_3\otimes I+y_3 I\otimes \sigma_3+\sum_{i,j=1}^3 R_{ij}\sigma_i\otimes \sigma_j),
\end{equation}

the reservoirs are in thermal equilibrium states at temperature T,

\begin{equation}
\rho_{E{_j}}=exp(-\beta \sum_k\omega_{j,k}b_{j,k}^\dagger b_{j,k})/Z_{E_j},
\end{equation}

and the whole system is in the product state,


\begin{equation}
\rho(0)=\rho_s(0)\otimes \rho_{E_1}\otimes\rho_{E_2}
\end{equation}

where $c_i$ is a real number with $0\leq \left\vert c_{i}\right\vert \leq 1$ for each i,$Z_{E_j}$ is the partition function of the reservoir j with
$Z_{E_j}=Tr(exp(-\beta \sum_k\omega_{j,k}b_{j,k}^\dagger b_{j,k}))=\prod_k(1-e^{-\beta \omega_{j,k}})^{-1}$, and $\beta=1/(k_{B}T)$.
The whole system?s state at time t is governed by

\begin{equation}
\rho(t)=exp(-iHt)\rho(0)exp(iHt),
\end{equation}

and the state of the two qubits at time t can be obtained by
the partial trace

\begin{equation}\label{2}
\rho_s(t)=Tr_E[exp(-iHt)\rho(0)exp(iHt)].
\end{equation}

In order to calculate GMQD, we need to write $\rho_s(t)$ in the form of Eq \cite{1}. From Eq \cite{2}, the elements
of $\rho_s(t)$ can be expressed as

\begin{equation}\label{2}
\left\langle mm^{\prime }\right\vert \rho_s(t)\left\vert nn^{\prime
}\right\rangle =Tr[\xi_1^{mn}(t)\xi_1^{m'n'}(t)\rho_(0)].
\end{equation}

where $\xi_1^{mn}(t)=e^{iHt}\xi_1^{mn}e^{-iHt}$ is the Heisenberg operator of
qubit j, and $\xi_1^{mn}=\left\vert n\right\rangle _{j}\left\langle m\right\vert$. Noting that $[\xi_1^{mn},H]=0$, we have $\xi_1^{mn}(t)=\xi_1^{mn}$.
The operators $\xi_1^{mn}(t)$ for $m\neq\ n$ can be obtained by solving the Heisenberg equations of motion,




\begin{equation}\label{}
i\frac{d}{dt}b_{j,k}=\omega_{j,k}b_{j,k}(t)+g_{j,k}\sigma_{j,3},
\end{equation}
\begin{equation}\label{}
i\frac{d}{dt}\xi_{j}^{01}(t)=-\nu_j\xi_{j}^{01}(t)-2\sum_k[g_{j,k}b_{j,k}^\dagger(t)+g_{j,k}^*b_{j,k}(t)]\xi_{j}^{01}(t)
\end{equation}

The solutions to the above differential equations are

\begin{equation}\label{}
b_{j,k}(t)=e^{-i\omega_{j,k}t}[b_{j,k}+\frac{1}{2}\alpha_{j,k}(t)\sigma_{j,3}]
\end{equation}

\begin{equation}\label{}
\xi_{j}^{01}(t)=\xi_{j}^{01}exp\left\{ {{i\nu_jt-\sum_k[\alpha_{j,k}(t)b_{j,k}^\dagger-\alpha_{j,k}(t)b_{j,k}]}}\right\},
\end{equation}

with
\begin{equation}\label{}
\alpha_{j,k}=2g_{j,k}(1-e^{i\omega_{j,k}(t)})/\omega_{j,k}.
\end{equation}

with the help of Eqs.(8) and (12),we finally obtain

\begin{equation}
\rho_s(t)=\frac{1}{4}(I\otimes I+T_{30} \sigma_3\otimes I+T_{03} I\otimes \sigma_3+T_{11}\sigma_1\otimes \sigma_1+T_{12} \sigma_1\otimes \sigma_2+T_{21} \sigma_2\otimes \sigma_1+T_{22} \sigma_2\otimes \sigma_2
+T_{33} \sigma_3\otimes \sigma_3),
\end{equation}

Where
\begin{equation}
T_{11}=[R_{11}cos(\nu_{1}t)cos(\nu_{2}t)+R_{22}sin(\nu_{1}t)sin(\nu_{2}t)
-R_{12}cos(\nu_{1}t)sin(\nu_{2}t)-R_{21}sin(\nu_{1}t)cos(\nu_{2}t)]e^{-\gamma (t)}$$
$$T_{12}=[R_{11}cos(\nu_{1}t)sin(\nu_{2}t)-R_{22}sin(\nu_{1}t)cos(\nu_{2}t)
+R_{12}cos(\nu_{1}t)cos(\nu_{2}t)-R_{21}sin(\nu_{1}t)sin(\nu_{2}t)]e^{-\gamma (t)}$$
$$T_{21}=[R_{11}sin(\nu_{1}t)cos(\nu_{2}t)-R_{22}cos(\nu_{1}t)sin(\nu_{2}t)
-R_{12}sin(\nu_{1}t)sin(\nu_{2}t)+R_{21}cos(\nu_{1}t)cos(\nu_{2}t)]e^{-\gamma (t)}$$
$$T_{22}=[R_{11}sin(\nu_{1}t)sin(\nu_{2}t)+R_{22}cos(\nu_{1}t)cos(\nu_{2}t)
+R_{12}sin(\nu_{1}t)cos(\nu_{2}t)+R_{21}cos(\nu_{1}t)sin(\nu_{2}t)]e^{-\gamma (t)}$$
$$T_{30}=x_3$$
$$T_{03}=y_3$$
$$T_{33}=R_{33}
\end{equation}
and

\begin{equation}\label{a}
\gamma(t)=\sum_{j=1}^2\sum_k4\left\vert {g_{j,k}}\right\vert^2\omega_{j,k}^{-2}coth(\frac{\beta\omega_{j,k}}{2})[1-cos(\omega_{j,k}t)].
\end{equation}


\begin{equation}
\rho_s(t)=\left(%
\begin{array}{cccc}
  \rho_{00}(t) & 0 & 0 &  \rho_{03}(t) \\
  0 &  \rho_{11}(t) &  \rho_{12}(t) & 0 \\
  0 &  \rho_{21}(t) &  \rho_{22}(t) & 0 \\
   \rho_{30}(t) & 0 & 0 &  \rho_{33}(t) \\
\end{array}%
\right)
\end{equation}

where
\begin{equation}\label{t}
\rho_{00}(t)=1+T_{30}+T_{03}+T_{33}=\rho_{00}$$
$$\rho_{03}(t)=T_{11}-T_{22}-iT_{12}-iT_{21}=e^{-\gamma (t)}{\rho_{03}[cos(\nu_1 t +\nu_2 t))-isin(\nu_1 t +\nu_2 t )]}$$
$$\rho_{11}(t)=1+T_{30}-T_{03}-T_{33}=\rho_{11}$$
$$\rho_{12}(t)=T_{11}+T_{22}+iT_{12}-iT_{21}=e^{-\gamma (t)}{\rho_{12}[cos(\nu_1 t -\nu_2 t )-isin(\nu_1 t -\nu_2 t )]}$$
$$\rho_{21}(t)=T_{11}+T_{22}-iT_{12}+iT_{21}=e^{-\gamma (t)}{\rho_{21}[cos(\nu_1 t -\nu_2 t )+isin(\nu_1 t -\nu_2 t )]}$$
$$\rho_{22}(t)=1-T_{30}+T_{03}-T_{33}=\rho_{22}$$
$$\rho_{30}(t)=T_{11}-T_{22}+iT_{12}+iT_{21}=e^{-\gamma (t)}{\rho_{30}[cos(\nu_1 t +\nu_2 t )+isin(\nu_1 t +\nu_2 t )]}$$
$$\rho_{33}(t)=1-T_{30}-T_{03}+T_{33}=\rho_{33}
\end{equation}


These matrix elements can be written also as follows

\begin{equation}
\rho_{00}(t)= \rho_{00}$$
$$\rho_{03}(t)= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}\rho_{03}$$
$$\rho_{11}(t)= \rho_{11}$$
$$\rho_{12}(t)= e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} \rho_{12}$$
$$\rho_{21}(t)= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t} \rho_{21}$$
$$\rho_{22}(t)=\rho_{22}$$
$$\rho_{30}(t)= e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}  \rho_{30}$$
$$\rho_{33}(t)=\rho_{33}
\end{equation}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Closest classical states to two qubit $X$ states }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
To begin, we shall present the procedure leading to  the closest classically correlated state to the two-qubit $X$ state ({\ref{eqn1}}).
 The Fano-Bloch representation (\ref{fano-X})
reads
\begin{equation}\label{12X}
\rho_{12} = \frac{1}{4} \bigg[  \sigma_{0}\otimes \sigma_{0} +   T_{03}
\sigma_{0}\otimes \sigma_{3} +  T_{30}
\sigma_{3}\otimes \sigma_{0} + \sum_{kl} T_{kl}
\sigma_{k}\otimes \sigma_{l} \bigg]
\end{equation}
where the correlation matrix elements are obtainable from (\ref{Rij})
modulo some obvious substitutions.
The geometric measure of quantum discord is defined as the distance
the state $\rho_{12}$ and its closest classical-quantum state
presenting zero discord \cite{Dakic2010}
\begin{equation}\label{definition-discord}
D_{\rm g}(\rho_{12}) = \min_{\chi_{12}} \vert\vert \rho_{12} -
\chi_{12} \vert\vert^2
\end{equation}
where the Hilbert-Schmidt norm is defined by $ \vert\vert X
\vert\vert^2 = {\rm Tr} (X^{\dagger}X)$ and the minimization is
taken over the set of all classical states. When the measurement is
performed on the  qubit 1, the classical states
write
\begin{equation}
\chi_{12} = p_1 \vert \psi_1 \rangle \langle \psi_1 \vert \otimes \rho_1^{2} + p_2 \vert \psi_2 \rangle \langle \psi_2 \vert \otimes \rho_2^{2}
\end{equation}
where $\{ \vert \psi_1 \rangle , \vert \psi_2 \rangle \}$ is an
orthonormal basis related to the qubit $1$,  $p_i$ $(i=1,2)$ stands
for probability distribution and $\rho_i^{2}$ $(i=1,2)$ is the
marginal density of the qubit 2. The classically correlated states
$\chi_{12} $ can also be written as
\begin{equation}\label{xi12}
\chi_{12} = \frac{1}{4}\Bigg[ \sigma_{0}\otimes \sigma_{0} +
\sum_{i=1}^3 te_i ~~\sigma_i \otimes \sigma_{0} + \sum_{i=1}^3
(s_+)_{i} ~\sigma_0 \otimes \sigma_{i} + \sum_{i,j=1}^3  e_i
(s_-)_{j} ~\sigma_i \otimes \sigma_{j}\Bigg]
\end{equation}
where
$$ t = p_1-p_2, \qquad e_i = \langle \psi_1 \vert \sigma_i \vert \psi_1 \rangle,  \qquad
(s_{\pm})_{j} = {\rm Tr}\big((p_1\rho_1^{2} \pm p_2\rho_2^{2})
\sigma_{j} \big).$$
It follows that the distance between the density matrix $\rho_{12}$
and the classical state $\chi_{12}$, as measured by Hilbert-Schmidt
norm, is then given by
\begin{equation}\label{HS12-1}
|| \rho_{12} - \chi_{12}||^2 = \frac{1}{4} \bigg[  (t^2 -2te_3 T_{30} + T^2_{30}) +  \sum_{i=1}^{3}
(T_{0i} -  (s_+)_{i})^2 + \sum_{i,j=1}^3   ( T_{ij} -  e_i(s_-)_{j})^2  \bigg]
\end{equation}
The  minimization of the distance (\ref{HS12-1}),  with respect to the parameters $t$, $(s_+)_{i}$ and
$(s_-)_{i}$, gives
\begin{equation}\label{minHS12}
t= e_3T_{30} $$
$$ (s_+)_{1}= 0 \quad (s_+)_{2} = 0 \quad  (s_+)_{3} = T_{03} $$
$$ (s_-)_{i} = \sum_{j=1}^3 e_j T_{ji}.
%$$  (s_-)_{1} = e_1 T_{11} \quad (s_-)_{2} = e_2 T_{22} \quad (s_-)_{3} = e_3 T_{33}
\end{equation}
Inserting  these solutions in (\ref{HS12-1}), one has
\begin{equation}\label{HS12}
|| \rho_{12} - \chi_{12}||^2 = \frac{1}{4} \bigg[  {\rm Tr}K - \vec{e}^t K \vec{e} \bigg]
\end{equation}
where the matrix $K$ is defined by
\begin{equation}\label{matrixK12}
K = xx^{\dagger} + TT^{\dagger}
\end{equation}
with
$$ x^{\dagger}  = (0,0,T_{30})\qquad T = \left(%
\begin{array}{ccc}
T_{11} &  T_{12} & 0 \\
T_{21}   & T_{22} & 0 \\
0 & 0 &T_{33}  \\
\end{array}%
\right).$$
From equation (\ref{HS12}), one see that the minimal value of
Hilbert-Schmidt distance (\ref{HS12}) is reached for the largest
eigenvalue of the matrix $K$. We denote by $\lambda_1$, $\lambda_2$ and $\lambda_3$ the eigenvalues
of the matrix $K$ (\ref{matrixK12}) corresponding to the $X$ state
(\ref{eqn1})  or equivalently (\ref{12X}).  They
 are given by
\begin{equation}\label{vp:lambda12}
\lambda_1 = 4(\vert \rho_{14}\vert +\vert \rho_{23}\vert )^2,\qquad
\lambda_2 = 4(\vert \rho_{14} \vert - \vert \rho_{23}\vert )^2,
\qquad \lambda_3= 2[(\rho_{11}- \rho_{33})^2 + (\rho_{22}-
\rho_{44})^2].
\end{equation}
To get the minimal value of the Hilbert-Schmidt distance (\ref{HS12}) and subsequently the amount
of geometric quantum discord, one compares
$\lambda_1$, $\lambda_2$ and $\lambda_3$. As $\lambda_1$ is always greater than  $\lambda_2$,
the largest eigenvalue  $\lambda_{\rm max}$ is $\lambda_1$ or
$\lambda_3$. It follows that the geometric discord is given by
\begin{equation}
D_g(\rho_{12}) = \frac{1}{4}~ {\rm min}\{ \lambda_1 + \lambda_2 ,
\lambda_2 + \lambda_3\}.\label{eq:GMQD_new}
\end{equation}
To write down the  explicit  expressions of the closest classical
state $\chi_{12}$ to $\rho_{12}$, one has to determine the
eigenvector $\vec{e}_{\rm max}$ associated with the largest
eigenvalue $\lambda_{\rm max}$. In this respect, two cases ( $\lambda_{\rm max} =
\lambda_1$ and $\lambda_{\rm max} = \lambda_3$) are separately
discussed. We begin by density matrices $\rho_{12}$ (\ref{eqn1}) whose entries satisfy
the condition $\lambda_{\rm max} = \lambda_3$. The associated
eigenvector  is given by $\vec{e}_3 = (0,0,1)$. Replacing in the set
of  constraints (\ref{minHS12}), one has
\begin{equation}
\chi^3_{12} = \frac{1}{4}\Bigg[ \sigma_{0}\otimes \sigma_{0} +
T_{30} ~~\sigma_3 \otimes \sigma_{0} + T_{03}~\sigma_0 \otimes
\sigma_{3} + T_{33}~\sigma_3 \otimes \sigma_{3}\Bigg]
\end{equation}
In the second situation, the eigenvector corresponding to
$\lambda_1$ is given by $\vec{e}_1 = (\cos \frac{\phi}{2}, -
\sin\frac{\phi}{2} , 0)$ where $e^{i\phi} =
\frac{\rho_{14}\rho_{23}}{|\rho_{14}||\rho_{23}|} $. Reporting the
components of $\vec{e}_1$ in (\ref{minHS12}), one gets the closest
classical state
\begin{equation}
\chi^1_{12} = \frac{1}{4}\Bigg[ \sigma_{0}\otimes \sigma_{0} +
T_{30} ~~\sigma_3 \otimes \sigma_{0} + \sum_{i=1}^{2} \sum_{j=1}^{2}
\tilde{T}_{ij}~\sigma_i \otimes \sigma_{j} \Bigg]
\end{equation}
where
$$ \tilde{T}_{11} = \cos \frac{\phi}{2} (\cos \frac{\phi}{2} T_{11} - \sin \frac{\phi}{2} T_{21} )
\qquad  \tilde{T}_{12} =  \cos \frac{\phi}{2} (\cos \frac{\phi}{2}
T_{12} - \sin \frac{\phi}{2} T_{22} )$$
$$  \tilde{T}_{21} =  -\sin \frac{\phi}{2} (\cos \frac{\phi}{2} T_{11} - \sin \frac{\phi}{2} T_{21} )
\qquad \tilde{T}_{22} =  -\sin\frac{\phi}{2} (\cos \frac{\phi}{2}
T_{12} - \sin \frac{\phi}{2} T_{22} ).$$



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Evolution of geometric quantum discord}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




\section{Geometric measure of quantum discord under quantum
  decoherence channels\label{sec:GMQD_channel}}

A quantum channel can be described in the Kraus
representation
\begin{equation}
  \mathcal{E}(\rho)=\sum_{\mu}K_{\mu}\rho
  K_{\mu}^{\dagger},
\end{equation}
where $K_{\mu}$ are Kraus operators satisfying $\sum_{\mu}K_{\mu}^{\dagger}K_{\mu}=\sigma_0$.  As we
discussed in the previous section, to obtain the GMQD, we need to know
the expectation values of the Pauli matrices of the two qubits for the
state $\mathcal{E}(\rho)$. So we turn to the Heisenberg picture to
describe quantum channels via the map~\cite{Wang2010}
\begin{equation}
  \mathcal{E}^{\dagger}(A)=\sum_{\mu}K_{\mu}^{\dagger}AK_{\mu}
\end{equation}
with $A$ an arbitrary observable. Then the expectation value of $A$
can be obtained through $\langle
A\rangle=\mathrm{Tr}\left[A\mathcal{E}(\rho)\right]=\mathrm{Tr}\left[\mathcal{E}^{\dagger}(A)\rho\right]$.
Because an arbitrary Hermitian operator on $\mathbb{C}^{2}$ can be
expressed by $A=\sum_{i=0}^{3}r_{i}\sigma_{i}$ with
$r_{i}\in\mathbb{R}$, then a quantum channel for a qubit can be
characterized by the transmission matrix $M$ defined through
\begin{equation}
  \mathcal{E}^{\dagger}(\sigma_{i})=\sum_{j}M_{ij}\sigma_{j}\quad\mathrm{or\quad}M_{ij}=\frac{1}{2}\mathrm{Tr}\left[\mathcal{E}^{\dagger}(\sigma_{i})\sigma_{j}\right].
  \label{eq:Transmission_Matrix}
\end{equation}
Since $\mathrm{Tr}[\mathcal{E}^{\dagger}(\sigma_{i})\rho]=\sum_j
M_{ij}\mathrm{Tr}[\sigma_{j}\rho]$, $M_{ij}$ actually describes the
transformation of the polarized vector
$P_{i}\equiv\mathrm{Tr}[\sigma_{j}\rho]$.

Now we consider the case of two qubits under local decoherence
channels, i.e.,
$\rho=[\mathcal{E}_{A}\otimes\mathcal{E}_{B}](\rho_{0})$. To obtain
the GMQD of the output state $\rho$ through the channel, we need to
get the expectation matrix $\mathcal{R}$. With the Heisenberg picture,
we have
\begin{equation}
  \mathcal{R}_{ij}  = \mathrm{Tr}(\mathcal{E}_{A}^{\dagger}(\sigma_{i})\otimes\mathcal{E}_{B}^{\dagger}(\sigma_{j})\rho_{0})
  =(M_{A}\mathcal{R}_{0}M_{B}^{T})_{ij},
\end{equation}
where $\mathcal{R}_{0}$ is the expectation matrix under $\rho_{0}$,
i.e.,
$(\mathcal{R}_{0})_{ij}=\mathrm{Tr}(\sigma_{i}\otimes\sigma_{j}\rho_{0})$,
and $M_{A(B)}$ is the transformation matrix characterizing the quantum
channel $\mathcal{E}_{A(B)}$. So we obtain
$\mathcal{R}=M_{A}\mathcal{R}_{0}M_{B}^{T}$.

For simplicity, we assume $\mathcal{E}^{A}$ and $\mathcal{E}^{B}$ be
identical, hereafter. Next, we consider three typical kinds of
decoherence channels: the amplitude damping channel (ADC), the phase
damping channel (PDC), and the depolarizing channel (DPC). They are
described by the set of Kraus operators
respectively~\cite{PreskillLect,NielsenBook}:
\begin{eqnarray}
  K^{\mathrm{ADC}} & = & \left\{ \sqrt{s}|0\rangle\langle0|+|1\rangle\langle1|,\ \sqrt{p}|1\rangle\langle0|\right\} ,
  \label{eq:Kraus_ADC}\\
  K^{\mathrm{PDC}} & = & \left\{ \sqrt{s}\sigma_0 ,\ \sqrt{p}|0\rangle\langle0|,\ \sqrt{p}|1\rangle\langle1|\right\} ,
  \label{eq:Kraus_PDC}\\
  K^{\mathrm{DPC}} & = & \{\frac{1}{2}\sqrt{1+3s}\,\sigma_0 ,\ \frac{1}{2}\sqrt{p}\,\sigma_{x},\frac{1}{2}\sqrt{p}\,\sigma_{y},\ \frac{1}{2}\sqrt{p}\,\sigma_{z}\},
  \label{eq:Kraus_DPC}
\end{eqnarray}
with $s\equiv 1-p$. Here the real parameter $p\in[0,1]$ may be
time-dependent in some realistic
setup~\cite{NielsenBook,PreskillLect}. For instance, for the PDC, the
parameter $s$ may be like $\exp(-\gamma t)$ with $\gamma$ the rate of
damping.

From Eqs.~(\ref{eq:Transmission_Matrix}), (\ref{eq:Kraus_ADC}),
(\ref{eq:Kraus_PDC}), and (\ref{eq:Kraus_DPC}), the transmission
matrix $M$ of each channel can be got through the transformation of
the Pauli matrices in the Heisenberg picture~\cite{Wang2010} as
\begin{equation}
  M_{\mathrm{ADC}}  =  \left[\begin{array}{cccc}
      1 & 0 & 0 & 0\\
      0 & \sqrt{s} & 0 & 0\\
      0 & 0 & \sqrt{s} & 0\\
      -p & 0 & 0 & s\end{array}\right],\quad
  M_{\mathrm{PDC}}  =  \left[\begin{array}{cccc}
      1 & 0 & 0 & 0\\
      0 & s & 0 & 0\\
      0 & 0 & s & 0\\
      0 & 0 & 0 & 1\end{array}\right],\quad
  M_{\mathrm{DPC}}  =  \left[\begin{array}{cccc}
      1 & 0 & 0 & 0\\
      0 & s & 0 & 0\\
      0 & 0 & s & 0\\
      0 & 0 & 0 & s\end{array}\right].
\end{equation}







For the density matrix under consideration, the eigenvalues of the matrix $k$ are given by
$$ \lambda_1 \equiv \lambda_1 (t) = 4 (\vert \rho_{12}(t)\vert + \vert \rho_{03}(t)\vert)^2
$$
$$ \lambda_2  \equiv \lambda_2 (t) =  4 (\vert \rho_{12}(t) \vert - \vert \rho_{03}(t) \vert )^2 $$
$$ \lambda_3  \equiv \lambda_3 (t)  = 2\big[ (\rho_{00}(t) - \rho_{22}(t))^2 + (\rho_{11}(t) - \rho_{33}(t))^2\big]$$

finally, we find:

$$\lambda_1 (t)=4(\left\vert {e^{-\gamma
(t)}e^{-i(\upsilon_1-\upsilon_2)t}\rho_{12}}\right\vert+\left\vert
{e^{-\gamma
(t)}e^{-i(\upsilon_1+\upsilon_2)t}\rho_{03}}\right\vert)^2=e^{-2\gamma(t)}\lambda_1(0)$$
$$ \lambda_2 (t)=4(\left\vert {e^{-\gamma
(t)}e^{-i(\upsilon_1-\upsilon_2)t}\rho_{12}}\right\vert+\left\vert
{e^{-\gamma
(t)}e^{-i(\upsilon_1+\upsilon_2)t}\rho_{03}}\right\vert)^2=e^{-2\gamma(t)}\lambda_2(0)$$
$$\lambda_3 (t)=2[(\rho_{00}-\rho_{22})^2+(\rho_{11}-\rho_{33})^2]=\lambda_3(0)$$


Clearly, the eigenvalue $\lambda_1$ is larger than $\lambda_2$.\\







\subsection{GMQD for dicke and their superpositions} we consider a
state of two qubits which have been extracted from the whole
ensemble. Specifically, we concentrate on the states with exchange
symmetry and parity, whose two-qubit reduced density matrix can be
written as

\begin{equation}\label{}
    \rho^{AB}=\left(%
\begin{array}{cccc}
  \upsilon_+ & 0 & 0 & u^* \\
  0 & y & y & 0 \\
  0 & y & y & 0 \\
  u & 0 & 0 & \upsilon_- \\
\end{array}%
\right)
\end{equation}

in the basis of ${\left\vert 00\right\rangle,\left\vert 01
\right\rangle,\left\vert 10\right\rangle,\left\vert
11\right\rangle}$, with $\upsilon_+$, $\upsilon_-$. and $y$ real,
and $u^*$. the complex conjugate of u. The elements of the density
matrix can be represented by the expectation values of the spin
components
\begin{equation}\label{}
    \upsilon_{\pm}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle\pm 4 \left\langle J_{z}\right\rangle
    (N-1)}{4N(N-1)},
\end{equation}



\begin{equation}\label{}
    y=\frac{N^2-4\left\langle
    J^2_{z}\right\rangle}{4N(N-1)} \\\\\\\\\\\\\\, u=\frac{\left\langle
    J^2_{+}\right\rangle}{N(N-1)}
\end{equation}


\subsubsection{Dicke states}

The so-called Dicke states under consideration are defined as
\begin{equation}
\left\vert n\right\rangle _{N}\equiv \left\vert N/2,-N/2+n\right\rangle ,%
\text{ \ \ \ \ \ \ \ \ \ \ }n=0,.....,N,
\end{equation}

As we have discussed in the previous section, to obtain the GMQD, we
need to calculate the expectation values of the spin components for
the stats $\left\vert n\right\rangle _{N}$, which are

\begin{eqnarray}
\left\langle J_{z}\right\rangle   = n-\frac{N}{2},  \quad
\left\langle J^2_{z}\right\rangle   = (n-\frac{N}{2})^{2},  \qquad
\left\langle J_{+}^{2}\right\rangle  = \left\langle
J_{-}^{2}\right\rangle =0.
\end{eqnarray}


%Two-qubit reduced density matrix can be written as
%\begin{equation}\label{}
 %   \rho^{AB}=\left(%
%\begin{array}{cccc}
%  \upsilon_+ & 0 & 0 & u^* \\
%  0 & y & y & 0 \\
%  0 & y & y & 0 \\
%  u & 0 & 0 & \upsilon_- \\
%\end{array}%
%\right)
%\end{equation}

%Where:
%\begin{equation}\label{}
%    \upsilon_{+}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle + 4 \left\langle J_{z}\right\rangle
%    (N-1)}{4N(N-1)}=\frac{n(n-1)}{N(N-1)},
%\end{equation}

%\begin{equation}\label{}
%    \upsilon_{-}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle - 4 \left\langle J_{z}\right\rangle
%    (N-1)}{4N(N-1)}=1+\frac{n^2-2nN+n}{N(N-1)},
%\end{equation}


%\begin{equation}\label{}
%    y=\frac{N^2-4\left\langle J^2_{z}\right\rangle}{4N(N-1)}=\frac{-n^2+nN}{N(N-1)},$$$$ u=\frac{\left\langle
%    J^2_{+}\right\rangle}{N(N-1)}=0
%\end{equation}

%Using the formula (\ref{t}) we can write the density matrix for each
%system at t after intercation.


%\begin{equation}
%\rho_{00}(t)= \rho_{00}=\upsilon_+=\frac{n(n-1)}{N(N-1)}$$
%$$\rho_{03}(t)= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}\rho_{03}= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}u^*=0$$
%$$\rho_{11}(t)= \rho_{11}=y=\frac{-n^2+nN}{N(N-1)}$$
%$$\rho_{12}(t)= e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} \rho_{12}=e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t}y
%=e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t}\frac{-n^2+nN}{N(N-1)} $$
%$$\rho_{21}(t)= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t} \rho_{21}= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t}y
%= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t}\frac{-n^2+nN}{N(N-1)}$$
%$$\rho_{22}(t)=\rho_{22}=y=\rho_{22}=\frac{-n^2+nN}{N(N-1)} $$
%$$\rho_{30}(t)= e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}\rho_{30}=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}u =0 $$
%$$\rho_{33}(t)=\rho_{33} =\upsilon_-=1+\frac{n^2-2nN+n}{N(N-1)}
%\end{equation}


%This  form will be more convenient for our purpose. For the density
%matrix under consideration, the eigenvalues of the matrix $k$ are
%given by
$$ \lambda_1 \equiv \lambda_1 (t) =e^{-2\gamma(t)}\lambda_1(0)=e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}
$$
$$ \lambda_2 \equiv \lambda_2 (t) = e^{-2\gamma(t)}\lambda_2(0)=e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}$$
$$ \lambda_3  \equiv \lambda_3 (t)= \lambda_3(0)=\frac{(N-2n)^2}{N^2}+\frac{[(N-2n)^2-N]^2}{N^2(N-1)^2}$$


$\lambda_1 (t)=\lambda_2 (t)$ and the maximum of the eigenvalues is
$\lambda_1 (t)$ or $\lambda_3 (t)$

%Therefore, for these X-states, equation (23) is reduced to a simpler
%form
%\begin{equation}\label{}
%    D_G(\rho)=\frac{1}{4}(min\{\lambda_1 (t),\lambda_3
%    (t)\}+\lambda_2
%    (t))
%\end{equation}

%Since they are functions of N and n, we will discuss the properties
%of the GMQD in one case with N and n fixed.

\subsubsection{Generalized GHZ states}

Now we discuss the GMQD of the generalized GHZ states which are
superpositions of two special Dicke states. The generalized GHZ
states which are important in the research of quantum mechanics are
states in the following form:

\begin{equation}\label{}
    \left\vert \psi\right\rangle_{GHZ}=cos\theta\left\vert
    0\right\rangle_N+e^{i\phi}sin\theta \left\vert
    N\right\rangle_N.
\end{equation}


For such states, the spin expectation values are


\begin{equation}\label{}
\left\langle J_{z}\right\rangle=-\frac{N}{2}cos 2\theta,$$
$$
\left\langle J^2_{z}\right\rangle=\frac{N^2}{4},$$
$$
\left\langle J_{+}^{2}\right\rangle=0.
\end{equation}


%Two-qubit reduced density matrix can be written as
%\begin{equation}\label{}
%    \rho^{AB}\left(%
%\begin{array}{cccc}
%  \upsilon_+ & 0 & 0 & u^* \\
%  0 & y & y & 0 \\
%  0 & y & y & 0 \\
%  u & 0 & 0 & \upsilon_- \\
%\end{array}%
%\right)
%\end{equation}

%Where:
%\begin{equation}\label{}
%    \upsilon_{+}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle + 4 \left\langle J_{z}\right\rangle
%    (N-1)}{4N(N-1)}=\frac{N-Ncos(2\theta)-1+cos(2\theta)}{2(N-1)},
%\end{equation}
%\begin{equation}\label{}
%    \upsilon_{-}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle - 4 \left\langle J_{z}\right\rangle
%    (N-1)}{4N(N-1)}=\frac{N+Ncos(2\theta)-1-cos(2\theta)}{2(N-1)},
%\end{equation}
%\begin{equation}\label{}
%y=\frac{N^2-4\left\langle J^2_{z}\right\rangle}{4N(N-1)}=0,$$$$ u=\frac{\left\langle
 %   J^2_{+}\right\rangle}{N(N-1)}=0.
%\end{equation}

%Using the formula (\ref{t}) we can write the density matrix for each
%system at t after intercation.

%\begin{equation}
%\rho_{00}(t)=
%\rho_{00}=\upsilon_+=\frac{N-Ncos(2\theta)-1+cos(2\theta)}{2(N-1)}$$
%$$\rho_{03}(t)= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}\rho_{03}= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}u^*=0$$
%$$\rho_{11}(t)= \rho_{11}=y=0$$
%$$\rho_{12}(t)= e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} \rho_{12}=e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} y=0 $$
%$$\rho_{21}(t)= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t} \rho_{21}=e^{-\gamma (t)} e^{i(\nu_1  -\nu_2 )t} y= 0 $$
%$$\rho_{22}(t)=\rho_{22}=y=0 $$
%$$\rho_{30}(t)=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}\rho_{30}=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}u=0 $$
%$$\rho_{33}(t)=\rho_{33} =\upsilon_-=\frac{N+Ncos(2\theta)-1-cos(2\theta)}{2(N-1)}
%\end{equation}

%This  form will be more convenient for our purpose. For the density
%matrix under consideration, the eigenvalues of the matrix $k$ are
%given by
$$ \lambda_1 \equiv \lambda_1 (t) = e^{-2\gamma(t)}\lambda_1(0)=0
$$
$$ \lambda_2  \equiv \lambda_2 (t) =e^{-2\gamma(t)}\lambda_2(0)=0 $$
$$ \lambda_3  \equiv \lambda_3 (t)  = 2\big[ (\rho_{00}(t) - \rho_{22}(t))^2 + (\rho_{11}(t) - \rho_{33}(t))^2\big]=\lambda_3(0)=1+cos^22\theta$$

which lead to
\begin{equation}\label{}
 D_g(\rho)=0
\end{equation}

\subsubsection{Superpositions of Dicke states}

In this subsection, weinvestigate a more general superposition of Dicke states, which
reads

\begin{equation}\label{}
    \left\vert \psi\right\rangle_{SD}=cos\theta\left\vert
    n\right\rangle_N+e^{i\phi}sin\theta \left\vert
    n+2\right\rangle_N, %
\text{ \ \ \ \ \ \ \ \ \ \ }n=0,.....,N-2
\end{equation}

with the angle $\theta\epsilon[0,\pi)$ and relative phase
$\varphi\epsilon[0,2\pi)$. Then the expressions of the relevant spin
expectation values are

\begin{eqnarray}
\left\langle J_{z}\right\rangle  &=&(n-\frac{N}{2})cos^2\theta(n+2-\frac{N}{2})sin^2\theta,  \nonumber \\
\left\langle J^2_{z}\right\rangle  &=&(n-\frac{N}{2})^2cos^2\theta(n+2-\frac{N}{2})^2sin^2\theta,  \nonumber \\
\left\langle J_{+}^{2}\right\rangle &=&\frac{1}{2}e^{i\phi}sin
2\theta\sqrt{\mu_n}. \nonumber
\\
\end{eqnarray}

Where
\begin{equation}\label{}
    \mu_N=(n+1)(n+2)(N-n)(N-n-1).
\end{equation}


%Two-qubit reduced density matrix can be written as
%\begin{equation}\label{}
 %   \rho^{AB}\left(%
%\begin{array}{cccc}
%  \upsilon_+ & 0 & 0 & u^* \\
%  0 & y & y & 0 \\
%  0 & y & y & 0 \\
%  u & 0 & 0 & \upsilon_- \\
%\end{array}%
%\right)
%\end{equation}

%Where:
%\begin{equation}\label{}
 %   \upsilon_{+}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle + 4 \left\langle J_{z}\right\rangle
 %   (N-1)}{4N(N-1)}=\frac{n^2+4nsin^2\theta-n+2sin^2\theta}{N(N-1)}
%\end{equation}

%\begin{equation}\label{}
 %   \upsilon_{-}=\frac{N^2-2N+4\left\langle J^2_{z}\right\rangle - 4 \left\langle J_{z}\right\rangle
 %   (N-1)}{4N(N-1)}=\frac{N^2-N+n^2-2nN+n+4nsin^2\theta+4sin^2\theta-4Nsin^2\theta}{N(N-1)}
%\end{equation}


%\begin{equation}\label{}
%    y=\frac{N^2-4\left\langle J^2_{z}\right\rangle}{4N(N-1)}=\frac{-n^2+nN-4nsin^2\theta-4sin^2\theta+2Nsin^2\theta}{N(N-1)},$$$$ u=\frac{\left\langle
%    J^2_{+}\right\rangle}{N(N-1)}=\frac{1}{2}\frac{e^{i\phi}sin
%2\theta\sqrt{\mu_n}}{N(N-1)}.
%\end{equation}

%Using the formula (\ref{t}) we can write the density matrix for each
%system at t after intercation

%\begin{equation}
%\rho_{00}(t)=
%\rho_{00}=\upsilon_+=\frac{n^2-n+4nsin^2\theta+2sin^2\theta}{N(N-1)}$$
%$$\rho_{03}(t)= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}\rho_{03}= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}u^*= e^{-\gamma (t)} e^{-i(\nu_1 +\nu_2
%)t}e^{-i\phi}\frac{sin 2\theta}{\sqrt{2N(N-1)}}$$
%$$\rho_{11}(t)= \rho_{11}=y=\frac{-n^2+nN-4nsin^2\theta-4sin^2\theta+2Nsin^2\theta}{N(N-1)}$$
%$$\rho_{12}(t)= e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} \rho_{12}=e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} y=e^{-\gamma
%(t)} e^{-i(\nu_1  -\nu_2
%)t}\frac{-n^2+nN-4nsin^2\theta-4sin^2\theta+2Nsin^2\theta}{N(N-1)}
%$$
%$$\rho_{21}(t)= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t} \rho_{21}=e^{-\gamma (t)} e^{i(\nu_1  -\nu_2 )t}y=e^{-\gamma (t)} e^{i(\nu_1  -\nu_2 )t}
%\frac{-n^2+nN-4nsin^2\theta-4sin^2\theta+2Nsin^2\theta}{N(N-1)}
%$$
%$$\rho_{22}(t)=\rho_{22}=y=\frac{-n^2+nN-4nsin^2\theta-4sin^2\theta+2Nsin^2\theta}{N(N-1)}$$
%$$\rho_{30}(t)=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}\rho_{30}=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}u=e^{-\gamma (t)}
%e^{i(\nu_1  +\nu_2 )t} e^{i\phi}\frac{sin 2\theta}{\sqrt{2N(N-1)}}$$
%$$\rho_{33}(t)=\rho_{33} =\upsilon_-=\frac{N^2-N+n^2-2nN+n+4nsin^2\theta+6sin^2\theta-4Nsin^2\theta}{N(N-1)}
%\end{equation}

The GMQD for these states can be obtained by repeating the previous
processing step. AsDG is not a function of ?, we discuss the GMQD
versus ? in the following. Specifically, when n = 0


%\begin{equation}
%\rho_{00}(t)= \rho_{00}=\upsilon_+=\frac{2sin^2\theta}{N(N-1)}$$
%$$\rho_{03}(t)= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}\rho_{03}= e^{-\gamma (t)} e^{-i(\nu_1  +\nu_2 )t}u^*= e^{-\gamma (t)} e^{-i(\nu_1 +\nu_2
%)t}\frac{sin 2\theta}{\sqrt{2N(N-1)}}$$
%$$\rho_{11}(t)= \rho_{11}=y=\frac{2(N-2)}{N(N-1)}sin^2\theta$$
%$$\rho_{12}(t)= e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} \rho_{12}=e^{-\gamma (t)} e^{-i(\nu_1  -\nu_2 )t} y=e^{-\gamma
%(t)} e^{-i(\nu_1  -\nu_2 )t}\frac{2(N-2)}{N(N-1)}sin^2\theta
%$$
%$$\rho_{21}(t)= e^{-\gamma (t)}  e^{i(\nu_1  -\nu_2 )t} \rho_{21}=e^{-\gamma (t)} e^{i(\nu_1  -\nu_2 )t}y=e^{-\gamma (t)} e^{i(\nu_1  -\nu_2 )t}
%\frac{2(N-2)}{N(N-1)}sin^2\theta
%$$
%$$\rho_{22}(t)=\rho_{22}=y=\frac{2(N-2)}{N(N-1)}sin^2\theta$$
%$$\rho_{30}(t)=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}\rho_{30}=e^{-\gamma (t)} e^{i(\nu_1  +\nu_2 )t}u=e^{-\gamma (t)}
%e^{i(\nu_1  +\nu_2 )t} e^{i\phi}\frac{sin 2\theta}{\sqrt{2N(N-1)}}$$
%$$\rho_{33}(t)=\rho_{33} =\upsilon_-=\frac{N^2-N+6sin^2\theta-4Nsin^2\theta}{N(N-1)}
%\end{equation}


%This  form will be more convenient for our purpose. For the density
%matrix under consideration

$$ \lambda_1 \equiv \lambda_1 (t) =
e^{-2\gamma(t)}\lambda_1(0)=4e^{-2\gamma(t)}[\frac{2(N-2)}{N(N-1)}sin^2\theta+\frac{\vert
sin 2\theta \vert}{\sqrt{2N(N-1)}}]^2
$$
$$ \lambda_2  \equiv \lambda_2 (t) =
e^{-2\gamma(t)}\lambda_2(0)=4e^{-2\gamma(t)}[\frac{2(N-2)}{N(N-1)}sin^2\theta-\frac{\vert
sin 2\theta \vert}{\sqrt{2N(N-1)}}]^2$$
$$ \lambda_3  \equiv \lambda_3 (t)  = 2\big[ (\rho_{00}(t) - \rho_{22}(t))^2 + (\rho_{11}(t) - \rho_{33}(t))^2\big]=(1-\frac{4}{N}sin^2\theta)^2+[1-\frac{8(N-2)}{N(N-1)}sin^2\theta]^2$$


\section{numerical analysis}

\subsection{Dicke
states}

For the density matrix under consideration, the eigenvalues of the
matrix $k$ are given by
$$ \lambda_1 \equiv \lambda_1 (t) =e^{-2\gamma(t)}\lambda_1(0)=e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}
$$
$$ \lambda_2 \equiv \lambda_2 (t) = e^{-2\gamma(t)}\lambda_2(0)=e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}$$
$$ \lambda_3  \equiv \lambda_3 (t)= \lambda_3(0)=\frac{(N-2n)^2}{N^2}+\frac{[(N-2n)^2-N]^2}{N^2(N-1)^2}$$


$\lambda_1(t)=\lambda_2 (t)$ and the maximum of the eigenvalues is
$\lambda_1(t)$ or $\lambda_3 (t)$

Therefore, for these X-states, equation (23) is reduced to a simpler
form
\begin{equation}\label{}
    D_G(\rho)=\frac{1}{4}(min\{\lambda_1 (t),\lambda_3
    (t)\}+\lambda_2
    (t))
\end{equation}

Specifically, When $N=15$ \text{\\\\\\} $n=0,......,15$

\begin{itemize}
    \item For $n=0,...,4$ and $n=11,....,15$

    The geometric quantum discord $\rho_{12}$ is
\begin{equation}\label{}
    D_{\rho}=\frac{1}{2}(\lambda_1)=\frac{1}{2}e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}
\end{equation}
when $\lambda_3>\lambda_1$ (*)
    \item For $n=5$ and $n=10$ the condition (*) is satisfied for $t_0<t$. In this situation, the
geometric quantum discord is
\begin{equation}\label{}
D_G=\frac{1}{2}e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2}
\end{equation}

when the transmission parameter $t$ satisfies $t_0>t$

we have $\lambda_3<\lambda_1$ and the geometric quantum discord is
given by
\begin{equation}\label{}
D_G=\frac{1}{4}(\lambda_3+\lambda_2)=\frac{1}{4}(\frac{(N-2n)^2}{N^2}+\frac{[(N-2n)^2-N]^2}{N^2(N-1)^2}+e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2})\end{equation}

    \item For $n=6,7$ and $n=8,9$

    the geometric quantum discord $\rho_{12}$ is
\begin{equation}\label{}
D_G=\frac{1}{4}(\lambda_3+\lambda_2)=\frac{1}{4}(\frac{(N-2n)^2}{N^2}+\frac{[(N-2n)^2-N]^2}{N^2(N-1)^2}+e^{-2\gamma(t)}\frac{4n^2(N-n)^2}{N^2(N-1)^2})\end{equation}
 when$\lambda_3<\lambda_1$
\end{itemize}

\begin{center}
  \includegraphics[width=3.5in]{n_0_sub_ohmic.eps}\\\label{1}
{\bf Figure 1:}  {(a) Dynamics of GMQD for the sub-Ohmic reservoirs
of Dicke states with $N$ fixed $(N=15)$ and $(n=1,..,4)$
$s=0.5,\lambda=0.1, \Omega \beta=1$ (numerical calculation with the
upper bound of m being $10^5$)}.
\end{center}

\begin{center}
  \includegraphics[width=3.5in]{n_5_sub_ohmic.eps}\\\label{1}
{\bf }  {}
\end{center}


\begin{center}
  \includegraphics[width=3.5in]{n_6_sub_ohmic.eps}\\\label{1}
{\bf }  {}
\end{center}


\begin{center}
  \includegraphics[width=3.5in]{n_7_sub_ohmic.eps}\\\label{1}
{\bf Figure 1:}  {(b) Dynamics of GMQD for the sub-Ohmic reservoirs
of Dicke states with $N$ fixed $(N=15)$ and $(n=5,6,7)$
$s=0.5,\lambda=0.1, \Omega \beta=1$ (numerical calculation with the
upper bound of m being $10^5$)}.
\end{center}

\begin{center}
  \includegraphics[width=3.5in]{22.eps}\\\label{2}
{\bf Figure 2:}  {(a) Dynamics of GMQD for the Ohmic reservoirs of
dicke states with $N$ fixed $(N=15)$ and $(n=1,..,4)$
$s=0.5,\lambda=0.1, \Omega \beta=1$ (numerical calculation with the
upper bound of m being $10^5$)}.
\end{center}



\begin{center}
  \includegraphics[width=3.5in]{n_5_ohmic.eps}\\\label{2}
{\bf}  {}.
\end{center}


\begin{center}
  \includegraphics[width=3.5in]{n_6_ohmic.eps}\\\label{2}
{\bf }  {}.
\end{center}

\begin{center}
  \includegraphics[width=3.5in]{n_7_ohmic.eps}\\\label{2}
{\bf Figure 2:}  {(b) Dynamics of GMQD for the Ohmic reservoirs of
dicke states with $N$ fixed $(N=15)$ and $(n=5,6,7)$
$s=0.5,\lambda=0.1, \Omega \beta=1$ (numerical calculation with the
upper bound of m being $10^5$)}.
\end{center}

\begin{center}
  \includegraphics[width=3.5in]{c.eps}\\\label{3}
{\bf Figure 3:}  {(a) Dynamics of GMQD for the super-Ohmic
reservoirs of Dicke states with $N$ fixed $(N=15)$ and $(n=1,..,4)$
$s=1.5,\lambda=0.2$ (numerical calculation with the upper bound of m
being $10^5$)}.
\end{center}


\begin{center}
  \includegraphics[width=3.5in]{d.eps}\\\label{3}
{\bf Figure 3:}  {(b) Dynamics of GMQD for the super-Ohmic
reservoirs of Dicke states with $N$ fixed $(N=15)$ and $(n=5)$
$s=1.5,\lambda=0.2$ (numerical calculation with the upper bound of m
being $10^5$)}.
\end{center}

\begin{center}
  \includegraphics[width=3.5in]{e.eps}\\\label{3}
{\bf Figure 3:}  {(c)Dynamics of GMQD for the super-Ohmic reservoirs
of Dicke states with $N$ fixed $(N=15)$ and $(n=0,..,15)$
$s=1.5,\lambda=0.2$ (numerical calculation with the upper bound of m
being $10^5$)}.
\end{center}


It can also be seen that there are two classes of the evolution of
GMQD. (a)$D_G$ is a monotonic decreasing function of time t with the
limit zero, if $\lambda_3>\lambda_1$. (b)$D_G$ is a piecewise
monotonic decreasing function with one turning point at the time
$t=\overline{t}$ defined by $\lambda_3=\lambda_1$ and then has the
limit zero as t tends to infinity, if $\lambda_3<\lambda_1$. These
two classes of the evolution of $D_G$ are plotted in fig 1 and 2.

In the figure 3, there are two possible classes of the evolution of
$D_G$.(a)and (b) $D_G$ is a monotonic decreasing function of time
with a constant limit $(1/2)\lambda_1$ or
$(1/4)(\lambda_3+\lambda_2)$, if $\lambda_3>\lambda_1$, or if
$\lambda_3<\lambda_1$. (b)$D_G$ is a piecewise monotonic decreasing
function with one turning point at the time  $t=\overline{t}$
defined by $\lambda_3=\lambda_1$, and then has a constant limit as t
tends to infinity, if $\lambda_3<\lambda_1$.


\subsection{Superpositions of Dicke states}



For the density matrix under consideration

\begin{equation}\label{x}
\lambda_1 \equiv \lambda_1 (t) =
e^{-2\gamma(t)}\lambda_1(0)=4e^{-2\gamma(t)}[\frac{2(N-2)}{N(N-1)}sin^2\theta+\frac{\vert
sin 2\theta \vert}{\sqrt{2N(N-1)}}]^2
$$
$$ \lambda_2  \equiv \lambda_2 (t) =
e^{-2\gamma(t)}\lambda_2(0)=4e^{-2\gamma(t)}[\frac{2(N-2)}{N(N-1)}sin^2\theta-\frac{\vert
sin 2\theta \vert}{\sqrt{2N(N-1)}}]^2$$
$$ \lambda_3  \equiv \lambda_3 (t)  = 2\big[ (\rho_{00}(t) - \rho_{22}(t))^2 + (\rho_{11}(t) - \rho_{33}(t))^2\big]=(1-\frac{4}{N}sin^2\theta)^2+[1-\frac{8(N-2)}{N(N-1)}sin^2\theta]^2
\end{equation}

It is noted that $D_G$ is a periodic function of $\theta$ with
period $\pi$. Accordingly, we plot the GMQD versus ? within one
period. To compare the eigenvalues $\lambda_1$ and $\lambda_3$, we
find that there are two cases.




%\begin{center}
%  \includegraphics[width=5in]{44.eps}\\\label{3}
%{\bf Figure 4:}  {Eigenvalues of K of the superpositions of Dicke
%states verus $\theta$ when $n=0$ with $N=3$ and $k=0.1,....0.9$}.
%\end{center}

\begin{center}
  \includegraphics[width=5in]{55.eps}\\\label{3}
{\bf Figure 5:}  {The GMQD of the superpositions of Dicke states
verus $\theta$ when $n=0$ with $N=3$ and $k=0.1,....0.9$}.
\end{center}

%\begin{center}
%  \includegraphics[width=5in]{66.eps}\\\label{3}
%{\bf Figure 6:}  {Eigenvalues of K of the superpositions of Dicke
%states verus $\theta$ when $n=0$ with $N=8$ and $k=0.1,....0.9$}.
%\end{center}

\begin{center}
  \includegraphics[width=3in]{77.eps}\\\label{3}
{\bf Figure 7:}  {The GMQD of the superpositions of Dicke states
verus $\theta$ when $n=0$ with $N=8$ and $k=0.1,....0.9$}.
\end{center}


When $N<8$, as can be seen in figure 4 (N=3), the maximum eigenvalue
is $\lambda_1$ in the middle region and $\lambda_3$ at the edge near
0 and $\pi$ for $k=0.6,0.7,0.8,0.9$ and for $k=0.1,0.2,0.3,0.4$ the
maximum eigenvalue is $\lambda_3$ in the middle region and at the
edge near 0 and $\pi$ and $\lambda_1$ between the middle region and
at the edge near 0.

In figure 5, there are three local maxima of $D_G$. The first and
third maxima display just at the point when $\lambda_2=\lambda_3$.
While the second one accurs at $\theta=\frac{\pi}{2}$ which is due
to the symmetry of equation (\ref{x}). We also observe that there
exist two special values of $\theta$, $\theta=\pi/3$ and
$\theta=2\pi/3$, where the GMQD is zero. These zero values of the
GMQD are induced by $\lambda_2=\lambda_3=0$ at these two points
(while $\lambda_1$ is the maximum).

However,When $N \geq 8$ as is shown in fig 6, $\lambda_3>\lambda_1$
no matter what $\theta$ is. Consequently, $D_G$ has only one maximum
at $\theta=\pi/2$, as is shown in figure 7, and it can analytically
be written as

\begin{equation}\label{}
    D_G=e^{-2\gamma(t)}[\frac{8(N-2)^2sin^4\theta}{N^2(N-1)^2}+\frac{sin^22\theta}{N(N-1)}]
\end{equation}


\section{Conclusion}
In this paper, in order to investigate the pairwise quantum
correlations, we have studied the dynamics of GMQD The model under
consideration consists of two qubits coupled with two independent
bosonic reservoirs described by Ohmic-like spectral densities. The
Hamiltonian of the model is described as Eq (\ref{z}), and the
two-qubit system is initially in the X states decoupled from the
environments, described by Eq (\ref{v}). We examine the evolution of
GMQD for the three types of reservoirs, sub-ohmic, ohmic, and
super-ohmic particularly the Dicke states and their superpositions.
$D_G$ is a monotonic decreasing function of time t or a poiecewise
monotonic decreasing function with one turning point before becoming
frozen phenomenon.

 \begin{thebibliography}{99}
\bibitem{11}C. H. Bennett, Phys. Rev. Lett. 68, 3121 (1992).

\bibitem{22} C. H. Bennett et al., Phys. Rev. Lett. 70, 1895 (1993).

\bibitem{33} A. Peres, Phys. Rev. Lett. 77, 1413 (1996).

\bibitem{44} W. K. Wootters, Phys. Rev. Lett. 80, 2245 (1998).

\bibitem{55} M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum
Information, (Cambridge University Press, Cambridge, (2000), Chap.
3.

\bibitem{66} N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Rev. Mod. Phys.
74, 145 (2002).

\bibitem{77} E. Knill and R. Laflamme, Phys. Rev. Lett. 81, 5672 (1998).

\bibitem{88} Groisman B, Popescu S and Winter A 2005 Phys. Rev. A 72 032317.

\bibitem{99} Dakic B, Vedral V and Brukner C 2010 Phys. Rev. Lett. 105190502.

\bibitem{111} Daoud M and Ahl Laamara R 2012 Phys. Lett. A 376 2361.

\bibitem{222} H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901
(2001).

\bibitem{333} L. Henderson and V. Vedral, J. Phys. A 34, 6899
(2001).

\bibitem{444} V. Vedral, Phys. Rev. Lett. 90, 050401 (2003).

\bibitem{555} J. Maziero, L. C. Cel\`eri, R. M. Serra, and V. Vedral, Phys,
Rev A 80, 044102 (2009).

\bibitem{666} S. Luo, Phys. Rev. A 77, 042303 (2008).

\bibitem{777} T. Werlang, S. Souza, F. F. Fanchini, and C. J. Villas Boas,
Phys. Rev. A 80, 024103 (2009).

\bibitem{888} R. F. Werner, Phys. Rev. A 40, 4277 (1989).

\bibitem{999} D. Girolami, and G. Adesso, Phys. Rev. A 83, 052108 (2011).

\end{thebibliography}
\end{document}





\section{Appendix: Ohmic-Like spectral densities}

Equation (\ref{eq:GMQD_new}) gives the GMQD of two qubits in the
bosonic reservoirs. $D[\rho_s(t)]$ is a time-dependent function
through the parameter $\gamma(t)$ defined by Eq (\ref{a}). In the
continuum limit, $\sum_k 4 g_{j,k}^2$ can be replaced by $\int
d\omega J_j(\omega)\delta(\omega_{j,k}-\omega)$. We further let
$J_1(\omega)=J_2(\omega)=J(\omega)$, which means that the two
resevoirs have the same spectral density. In this case, there is
\begin{equation}\label{b}
    \gamma(t)=2\int_0^\infty d\omega J(\omega)\omega^{-2}\coth(\beta
    \omega/2)[1-\cos(\omega t)].
\end{equation}

For Ohmic-like spectral densities, $J(\omega)$ reads
\begin{equation}\label{}
    J(\omega)=\lambda \Omega^{1-s} \omega^s e^{-\omega/\Omega}
\end{equation}

where $\lambda$ is a dimensionless coupling constant and $\Omega$ is
the cutoff frequency. $0<s<1$, $s=1$, and $s>1$ are respectively
corresponding to sub-Ohmic, Ohmic, and super-Ohmic reservoirs. In
the following paragraphs, we discuss the evolution of $D[\rho_s(t)]$
in three kinds of Ohmic-like reservoirs, respectively.

\subsection{Sub-Ohmic reservoirs}

For the case of sub-ohmic reservoirs, $0<s<1$, the integral
(\ref{b}) can be calculated
\begin{equation}\label{}
    \gamma(t)=\frac{2\lambda
    \Gamma(s)}{s-1}\left\{ {1-(1+\Omega^2t^2)^{(1-s)/2}cos[(s-1)arctan(\Omega
    t)]}\right\}$$$$+\frac{4\lambda \Gamma(s)}{s-1}\sum_{m=1}^\infty(1+m\Omega
    \beta)^{1-s}\times \left\{ {1-[1+(\frac{\Omega t}{1+m\Omega
    \beta})^2]^{(1-s)/2}\times cos[(s-1)arctan(\frac{\Omega t}{1+m\Omega
    \beta})]}\right\}
\end{equation}

where $ \Gamma(s)$ is the function. To investigate the evolution of
$D[\rho_s(t)]$, we calculate the derivative of $\gamma(t)$ for time
t,

\begin{equation}\label{}
    \frac{\partial \gamma(t)}{\partial t}=2 \lambda \Gamma(s)\Omega \left\{ {(1+\Omega^2t^2)^{-s/2}\sin[s arctan(\Omega t)
    ]+2\sum_{m=1}^\infty [(1+m\Omega \beta)^2+\Omega^2t^2]^{-s/2}\times \sin[s arctan(\frac{\Omega t}{1+m\Omega \beta})]}\right\}.
\end{equation}
 Since $0<\arctan[\Omega t/(1+m\Omega \beta)]<\arctan(\Omega
 t)<\pi/2$, there $\partial \gamma(t)/\partial t>0$ for all time t.
Hence, $\gamma(t)$ is a monotonic increasing function of t, and it
tends to infinity as t goes to infinity.



\subsection{Ohmic reservoirs}

For the case of Ohmic reservoirs, s=1, the integral (\ref{b})can be
worked out as

\begin{equation}\label{}
\gamma(t)=\lambda[ln(1+\Omega^2t^2)+4ln\Gamma(1+\frac{1}{\Omega
t})-2ln\left\vert {\Gamma(1+\frac{1}{\Omega
\beta}+i\frac{t}{\beta})}\right\vert^2]
\end{equation}

The derivative of $\gamma(t)$ is
\begin{equation}\label{}
    \frac{\partial \gamma(t)}{\partial t}=2\lambda \Omega^2t[\frac{1}{1+\Omega^2t^2}+\sum_{m=1}^\infty\frac{2}{(1+m\Omega
    \beta)^2+\Omega^2t^2}].
\end{equation}

Since $\partial \gamma(t)/\partial t>0$, $\gamma(t)$ is monotonic
increasing for all time t,and it tends to infinity as T goes to
infinity.



\subsection{Super-Ohmic reservoirs}
For the case of super-Ohmic reservoirs, $s > 1$, the integral
(\ref{b}) can be worked out as

\begin{equation}\label{}
\gamma(t)=2\lambda
\Gamma(s-1)\left\{{1-(1+\Omega^2t^2)^{(1-s)/2}}\right\}\times
cos[(s-1)arctan(\Omega t)$$$$+4\lambda \Gamma(s-1)\sum_{m=1}^\infty
(1+m \Omega \beta)^{(1-s)}\times \left\{{1-[1+(\frac{\Omega
t}{1+m\Omega \beta})^2]^{(1-s)/2}\times cos[(s-1)arctan(\frac{\Omega
t}{1+m\Omega \beta})]}\right\}
\end{equation}

and the time derivative of it reads


\begin{equation}\label{}
    \frac{\partial \gamma(t)}{\partial t}=2 \lambda \Gamma(s)\Omega \left\{ {(1+\Omega^2t^2)^{-s/2}\sin[s arctan(\Omega t)
    ]+2\sum_{m=1}^\infty [(1+m\Omega \beta)^2+\Omega^2t^2]^{-s/2}\times \sin[s arctan(\frac{\Omega t}{1+m\Omega \beta})]}\right\}.
\end{equation}
 Since $0<\arctan[\Omega t/(1+m\Omega \beta)]<\arctan(\Omega
 t)<\pi/2$, there $\partial \gamma(t)/\partial t>0$ for all time t.
Hence, $\gamma(t)$ is a monotonic increasing function of t, and it
tends to infinity as t goes to infinity.

