Thus, we propose the two-mode phase state $|\theta_{p},\theta_{q} \rangle$ which should be necessarily introduced as $$\begin{aligned}
\label{pq2mode}
  \hspace{-1cm} |\theta_{p},\theta_{q} \rangle = \frac{1}{s+1}\sum_{n=0}^{s}\sum_{m=0}^{s} \exp\left(i n \theta_{p}\right) \exp\left( i  m \theta_{q}\right) | n,m \rangle,
  \end{aligned}$$ with $$\begin{aligned}
\label{thetapq}
  \theta_{k}=\theta_{0}+\frac{2\pi k}{s+1},\;\;\;\;k=p,q,
  \end{aligned}$$ and $\theta_{0}$ is an arbitrary value. In the same manner which stated in (\[thetap\]), each of the indices, $p$ and $q$ can get the values $0,1,...,s$. By the way, we are able to define the two-mode Pegg-Barnett phase distribution function as follows $$\begin{aligned}
\label{pthetapq}
  \mathcal{P}_{\theta}(\theta_{p},\theta_{q})=\lim_{s\rightarrow +\infty} \left( \frac{s+1}{2\pi} \right)^{2} \langle \theta_{p},\theta_{q}|\hat{\rho}_{F}|\theta_{p},\theta_{q} \rangle,
  \end{aligned}$$ By putting the relation (\[pq2mode\]) into (\[pthetapq\]) and paying attention to the state vector of the whole system proposed in (\[say\]), which all of its time-dependent coefficients have been explicitly obtained, one may obtain the phase distribution function associated with the two-mode cavity field as $$\begin{aligned}
\label{pthetapqfinal}
  \hspace{-2cm}\mathcal{P}_{\theta}(\theta_{1},\theta_{2}) &=& \frac{1}{4 \pi^{2}} \left| \sum_{n_{1}=0}^{+\infty}\sum_{n_{2}=0}^{+\infty} q_{n_{1}} q_{n_{2}} A(n_{1},n_{2},t) \exp(- i n_{1} \theta_{1}) \exp(- i n_{2} \theta_{2}) \right|^{2} \nonumber \\
  \hspace{-2cm}&+& \frac{1}{4 \pi^{2}} \left| \sum_{n_{1}=0}^{+\infty}\sum_{n_{2}=0}^{+\infty} q_{n_{1}} q_{n_{2}} B(n_{1}+1,n_{2},t) \exp(- i n_{1} \theta_{1}) \exp(- i n_{2} \theta_{2})  \right|^{2} \nonumber \\
  \hspace{-2cm}&+& \frac{1}{4 \pi^{2}} \left| \sum_{n_{1}=0}^{+\infty}\sum_{n_{2}=0}^{+\infty} q_{n_{1}} q_{n_{2}} C(n_{1},n_{2}+1,t) \exp(- i n_{1} \theta_{1}) \exp(- i n_{2} \theta_{2})  \right|^{2}.
 \end{aligned}$$ Now, based on the 