\section{Introduction}
  Nowadays, quantum chromodynamics (QCD)
  is established as the fundamental theory of the strong interaction.
  However, the non-perturbative properties of QCD
  are not fully understood, especially,
  on color confinement and chiral symmetry breaking.
  It is interesting issue to
  investigate the correspondence between these non-perturbative phenomena
  \cite{Suganuma:1995,Miyamura:1995,Gattringer:2006,
  Bruckmann:2007,Bilgici:2008,Synatschke:2008,Suganuma:2011,Gongyo:2012}.
  The lattice QCD calculations show
  the simultaneous chiral and deconfined phase transitions at finite temperature \cite{Karsch:2002}, which
  suggests a close relation between confinement and chiral
  symmetry breaking.

  As for chiral symmetry breaking,
  the chiral condensate $\langle \bar{q} q \rangle$
  is directly connected to the Dirac operator in QCD.
  The chiral condensate is proportional to the Dirac zero-mode density as
  \begin{equation}
    \langle \bar{q} q \rangle
    = - \lim_{m \rightarrow 0}
    \lim_{V\rightarrow \infty} \pi \langle \rho(0) \rangle,
    \label{}
  \end{equation}
  with the Dirac spectral density $\rho(\lambda)$,
  which is known as the Banks-Casher relation \cite{BanksCasher}.
  The Dirac zero-modes are also related to the topological charge
  via the Atiyah-Singer index theorem \cite{AtiyahSinger}.

  Therefore, it is interesting to investigate
  color confinement in terms of the Dirac-operator properties.
  Using Gattringer's formula \cite{Gattringer:2006},
  the Polyakov loop was analyzed by the Dirac spectrum sum
  with twisted boundary condition
  on lattice \cite{Bruckmann:2007,Bilgici:2008,Synatschke:2008}.
  In our previous studies \cite{Suganuma:2011,Gongyo:2012},
  we developed the Dirac-mode expansion method
  for the link-variable, and analyzed the role of the Dirac mode to
  the Wilson loop and the interquark potential.
  As for the hadron spectra, it is reported that
  the hadrons still exist as the bound state even without chiral symmetry breaking
  by removing low-lying Dirac-modes \cite{Lang:2011,Glozman:2012}.

  In this paper,
  based on the Dirac-mode expansion method \cite{Suganuma:2011,Gongyo:2012},
  we investigate the role of the Dirac mode to
  the Polyakov loop in both confined and deconfined phases
  at finite temperature in SU(3) lattice QCD.
  In Sec.2, we briefly review the formalism of the Dirac-mode
  expansion method in lattice QCD.
  In Sec.3, we show the lattice QCD results.
  Section 4 is devoted for the summary.

\section{Formalism}
  In this section, we briefly review the Dirac-mode expansion method
  in lattice QCD \cite{Suganuma:2011,Gongyo:2012},
  and formulation of the Dirac-mode projected Polyakov loop.

  \subsection{Dirac-mode expansion in lattice QCD}
  In lattice QCD, the Dirac operator $\Slash{D} = \gamma_\mu D_\mu$ is expressed as
  \begin{equation}
    \Slash{D}_{x,y} \equiv \frac{1}{2a}
    \sum_{\mu = 1}^4
    \gamma_\mu \left[ U_\mu(x) \delta_{x+\hat{\mu},y}
    - U_{-\mu}(x) \delta_{x-\hat{\mu},y}\right],
  \end{equation}
  using the link-variable $U_\mu(x) \in \mathrm{SU}(N_c)$
  and the lattice spacing $a$.
  Here, we use the convenient notation of $U_{-\mu}(x) \equiv U_\mu^\dagger(x-\hat{\mu})$,
  and $\hat{\mu}$ denotes for the unit vector on lattice in $\mu$-direction.
  In this paper, we adopt hermitian $\gamma$-matrices,
  i.e., $\gamma_\mu^\dagger = \gamma_\mu$, in the Euclidean space-time.
  Thus, the Dirac operator $\Slash{D}$ is anti-hermitian,
  and the Dirac eigenvalues are pure imaginary number.
  We introduce the normalized Dirac eigenstate $|n \rangle$ which satisfies
  \begin{equation}
    \Slash{D}|n\rangle = i\lambda_n | n \rangle,
  \end{equation}
  with the eigenvalue $\lambda_n \in \mathbf{R}$.
  The Dirac eigenfunction $\psi_n(x)$ defined by
  \begin{equation}
    \psi_n(x) \equiv \langle x | n \rangle
    \label{}
  \end{equation}
  satisfies $\Slash{D}\psi_n = i \lambda_n \psi_n$.

  We consider the operator formalism in lattice QCD \cite{Suganuma:2011,Gongyo:2012}.
  The link-variable operator $\hat{U}_\mu$ is defined
  by the matrix element
  \begin{equation}
    \langle x | \hat{U}_\mu | y \rangle = U_{\mu} (x) \delta_{x+\hat{\mu},y}.
    \label{}
  \end{equation}
  The Dirac-mode matrix element $\langle n | \hat{U}_\mu | m\rangle$
  is expressed as
  \begin{eqnarray}
    \langle n | \hat{U}_\mu | m \rangle
    &=& \sum_x \langle n | x \rangle \langle x |
    \hat{U}_\mu | x + \hat{\mu} \rangle
    \langle x + \hat{\mu} | m \rangle \nonumber \\
    &=& \sum_x \psi_n^\dagger(x) U_\mu(x) \psi_m(x+\hat{\mu}),
    \label{}
  \end{eqnarray}
  using the link-variable $U_\mu(x)$ and the Dirac eigenfunction $\psi_n(x)$.

  Using the completeness relation $\sum_n |n \rangle \langle n|=1$,
  any operator $\hat O$ can be expanded
  in terms of the Dirac-mode basis $|n \rangle$ as
  \begin{eqnarray}
    \hat O=\sum_n \sum_m |n \rangle \langle n|\hat O|m \rangle \langle m|,
    \label{eq:dirac-mode-expansion}
  \end{eqnarray}
  which is the mathematical basis of
  the Dirac-mode expansion \cite{Suganuma:2011,Gongyo:2012}.

  Based on the expansion in Eq.(\ref{eq:dirac-mode-expansion}),
  we introduce the Dirac-mode projection operator $\hat{P}$ as
  \begin{equation}
    \hat{P} \equiv \sum_{n \in \mathcal{A}} | n \rangle \langle n |,
    \label{}
  \end{equation}
  with the Dirac eigenstate $|n\rangle$,
  and arbitrary set of eigenmode subspace $\mathcal{A}$.
  For example, the IR and the UV Dirac-mode cut are given by
  \begin{equation}
    \hat{P}_{\rm IR} \equiv \sum_{|\lambda_n| \geq \Lambda_{\rm IR}} | n \rangle \langle n |,
    \qquad
    \hat{P}_{\rm UV} \equiv \sum_{|\lambda_n| \leq \Lambda_{\rm UV}} | n \rangle \langle n |,
  \end{equation}
  respectively, with the IR/UV cutoff parameter,
  $\Lambda_{\rm IR}$ and $\Lambda_{\rm UV}$.
  We define the Dirac-mode projected link-variable operator as
  \begin{equation}
    \hat{U}_\mu^P \equiv \hat{P} \hat{U}_\mu \hat{P}
    = \sum_{n \in \mathcal{A}} \sum_{m \in \mathcal{A}}
    | n \rangle \langle n | \hat{U}_\mu | m \rangle \langle m |,
    \label{}
  \end{equation}
  with the projection operator $\hat{P}$.
  Using the projected link-variable $\hat{U}_\mu^P$,
  we can analyze the individual contribution
  of each Dirac eigenmode to the various quantities,
  such as the Wilson loop \cite{Suganuma:2011,Gongyo:2012}.

  \subsection{Polyakov loop operator and Dirac-mode projection}
  Next, we formulate the Dirac-mode projection of the Polyakov loop.
  Hereafter, we consider the periodic SU(3) lattice of
  the space-time volume $V = L^3 \times N_t$ with lattice spacing $a$.
  In lattice QCD operator formalism,
  the Polyakov-loop operator is defined by
  \begin{equation}
    \hat{L}_P \equiv \frac{1}{3V} \prod_{i=1}^{N_t} \hat{U}_4
             =\frac{1}{3V} \hat{U}_4^{N_t}
    \label{}
  \end{equation}
  with the temporal link-variable operator $\hat{U}_4$.
  Taking the functional trace ``Tr'',
  we obtain the standard Polyakov loop $\langle L_P \rangle$ as
  \begin{eqnarray}
    \mathrm{Tr} \ \hat{L}_P
    &=& \frac{1}{3V} \mathrm{Tr} \
    \large\{ \prod_{i=1}^{N_t} \hat{U}_4 \large\} = \frac{1}{3V} \mathrm{tr}
    \sum_{\vec{x},t} \langle \vec{x},t | \prod_{i=1}^{N_t} \hat{U}_4 | \vec{x},t \rangle \nonumber  \\
    &=& \frac{1}{3V} \mathrm{tr} \sum_{\vec{x},t} \langle \vec{x},t | \hat{U}_4 | \vec{x},t + a \rangle
      \langle \vec{x},t+a | \hat{U}_4 | \vec{x},t+2a \rangle
      \cdots \langle \vec{x},t + (N_t-1)a | \hat{U}_4 | \vec{x},t \rangle \nonumber \\
    &=& \frac{1}{3V} \mathrm{tr}
      \sum_{\vec{x},t} U_4(\vec{x},t) U_4(\vec{x},t+a) \cdots U_4(\vec{x},t+(N_t-1)a)
      = \langle L_P \rangle,
  \end{eqnarray}
  where ``tr'' denotes the trace over SU(3) color index.

  Using the projection operator $\hat P$,
  we define the Dirac-mode projected Polyakov loop
  $\langle L_{P}^{\rm proj.} \rangle$ as
  \begin{eqnarray}
    \langle L_{P}^{\rm proj.} \rangle
    &\equiv& \frac{1}{3V} \mathrm{Tr} \ \large\{ \prod_{i=1}^{N_t} \hat{U}_4^P \large\} =
    \frac{1}{3V} \mathrm{Tr} \left\{
    \hat{P}\hat{U}_4 \hat{P} \hat{U}_4 \hat{P} \cdots \hat{P} \hat{U}_4 \hat{P} \right\} \nonumber \\
    &=& \frac{1}{3V}\mathrm{tr} \sum_{n_1,n_2,\dots, n_{N_t} \in \mathcal{A}}
      \langle n_1 | \hat{U}_4 | n_2 \rangle \langle n_2 | \hat{U}_4 | n_3 \rangle
      \cdots \langle n_{N_t} | \hat{U}_4 | n_1 \rangle.
    \label{eq:Lproj}
  \end{eqnarray}
  In particular, we consider the IR/UV Dirac-mode projected-Polyakov loop as
  \begin{eqnarray}
    \langle L_P \rangle_{\rm IR}
    &\equiv& \frac{1}{3V} \mathrm{tr}
    \sum_{|\lambda_{n_i}| \ge \Lambda_{\rm IR}}
    \langle n_1 | \hat{U}_4 | n_2 \rangle
    \cdots \langle n_{N_t} | \hat{U}_4 | n_1 \rangle, \\
    \langle L_P \rangle_{\rm UV}
    &\equiv& \frac{1}{3V} \mathrm{tr}
    \sum_{|\lambda_{n_i}| \le \Lambda_{\rm UV}}
    \langle n_1 | \hat{U}_4 | n_2 \rangle
    \cdots \langle n_{N_t} | \hat{U}_4 | n_1 \rangle,
    \label{}
  \end{eqnarray}
  with the IR/UV eigenvalue cutoff, $\Lambda_{\rm IR}$ and $\Lambda_{\rm UV}$.

\section{Lattice QCD calculation}
  In this section, we study
  the Polyakov loop in terms of the Dirac-mode
  in SU(3) lattice QCD at the quenched level.
  We use the LAPACK package for the full diagonalization of the Dirac operator \cite{LAPACK}.
  We use the Kogut-Susskind (KS) formalism for reduction of the computational costs
  \cite{Suganuma:2011,Gongyo:2012}.

  \subsection{The confined phase}
  First, we analyze the Polyakov loop properties in the confined phase.
  Here, we use $6^4$ lattice with $\beta = 5.6$, which
  corresponds to lattice spacing $a \simeq 0.25$fm \cite{Suganuma:2011,Gongyo:2012}.
  The total number of KS Dirac-modes is $L^3 \times N_t \times 3 = 3888$.
  Figure \ref{fig:DiracSpectrumCut} shows
  the lattice QCD result for the Dirac spectral density $\rho(\lambda)$
  and IR/UV-cut Dirac spectral density,
  \begin{equation}
    \rho_{\rm IR}(\lambda) \equiv \rho(\lambda)\theta(|\lambda|-
    \Lambda_{\rm IR}),
    \quad
    \rho_{\rm UV}(\lambda) \equiv \rho(\lambda)
    \theta( \Lambda_{\rm UV} - | \lambda | ),
  \end{equation}
  with $\Lambda_{\rm IR} = 0.5a^{-1}$ and $\Lambda_{\rm UV} = 2.0a^{-1}$.
  Both mode-cuts correspond to removing about 400 modes from full eigenmodes.

\begin{figure}
  \centering
  \includegraphics[width=0.325\textwidth,clip]{dirac_spectrum_confined.eps}
  \includegraphics[width=0.325\textwidth,clip]{dirac_spectrum_confined_ir_cut.eps}
  \includegraphics[width=0.325\textwidth,clip]{dirac_spectrum_confined_uv_cut.eps}
  \caption{\label{fig:DiracSpectrumCut}
    The Dirac spectral density in SU(3) lattice QCD
    on $6^4$ at $\beta$=5.6, i.e., $a\simeq$0.25fm.
  (a) The original spectral density $\rho(\lambda)$.
  (b) $\rho_{\rm IR}(\lambda)$ for IR-cut with $\Lambda_{\rm IR} = 0.5a^{-1}$.
  (c) $\rho_{\rm UV}(\lambda)$ for UV-cut with $\Lambda_{\rm UV} = 2.0a^{-1}$.
  }
\end{figure}

  We show in Figs.\ref{fig:PolyakovScatterConfined}(a)-(c)
  the scatter plot of the original Polyakov loop $\langle L_P \rangle$,
  $\langle L_P \rangle_{\rm IR}$ for low-lying Dirac-mode cut
  with $\Lambda_{\rm IR} = 0.5a^{-1}$,
  and $\langle L_P \rangle_{\rm UV}$ for UV Dirac-mode cut
  with $\Lambda_{\rm UV} = 2.0a^{-1}$, respectively.
  As shown in Fig.\ref{fig:PolyakovScatterConfined}(a),
  the Polyakov loop satisfies $\langle L_P \rangle \simeq 0$,
  which indicates the confined phase.

  By removing low-lying Dirac-modes,
  chiral symmetry breaking is effectively restored
  \cite{Gongyo:2012,BanksCasher,Lang:2011,Glozman:2012}.
  Actually, this IR Dirac-mode cut of
  $\Lambda_{\rm IR} = 0.5a^{-1} \simeq 0.4$GeV
  corresponds to about 98\% reduction of the quark condensate
  around the physical region $m_q \simeq 5$MeV \cite{Gongyo:2012}.
  However, as shown in Fig.\ref{fig:PolyakovScatterConfined}(b),
  the Polyakov loop $\langle L_P \rangle_{\rm IR}$ remains almost zero,
  which means unbroken $Z_3$ center symmetry.
  This result indicates that
  the single-quark energy is still extremely large, and the system remains
  in the confined phase even without chiral symmetry breaking.
  In the UV Dirac-mode cut, the chiral condensate is almost unchanged,
  and the Polyakov loop $\langle L_P \rangle_{\rm UV}$
  also remains almost zero,
  as shown in Fig.\ref{fig:PolyakovScatterConfined}(c).
  These results in Figs.\ref{fig:PolyakovScatterConfined}(b) and (c)
  show that the Polyakov loop is insensitive to the IR/UV Dirac-mode cut.

\begin{figure}
  \centering
  \includegraphics[width=0.27\textwidth,clip]{scatter-beta5.6-no-cut.eps}
  \includegraphics[width=0.27\textwidth,clip]{scatter-beta5.6-ir-cut.eps}
  \includegraphics[width=0.27\textwidth,clip]{scatter-beta5.6-uv-cut.eps}
\vspace{-0.2cm}
  \caption{\label{fig:PolyakovScatterConfined}
    The scatter plot of the Polyakov loop in the confined phase
    in SU(3) lattice QCD on $6^4$ at $\beta$=5.6, i.e.,
    $a\simeq$0.25fm and $T \equiv 1/(N_ta) \simeq$ 0.13GeV.
    (a) The original Polyakov loop $\langle L_P \rangle$.
    (b) $\langle L_P \rangle_{\rm IR}$ for low-lying Dirac-mode cut
        with $\Lambda_{\rm IR} = 0.5a^{-1}$.
    (c) $\langle L_P \rangle_{\rm UV}$ for UV Dirac-mode cut
        with $\Lambda_{\rm UV} = 2.0a^{-1}$.
  }
\end{figure}

\subsection{The deconfined phase at high temperature}
  Next, we study the role of the Dirac mode
  in the deconfined phase at high temperature.
  Here, we use $6^3 \times 4$ lattice at $\beta$=6.0.
  The total number of KS Dirac eigenmodes is $L^3 \times N_t \times 3$=2592.

  Figure~\ref{fig:PolyakovScatterDeconfined} shows
  the original Polyakov loop $\langle L_P \rangle$,
  $\langle L_P \rangle_{\rm IR}$ for low-lying Dirac-mode cut
        with $\Lambda_{\rm IR} = 0.5a^{-1}$, and
  $\langle L_P \rangle_{\rm UV}$ for UV Dirac-mode cut
        with $\Lambda_{\rm UV} = 2.0a^{-1}$, respectively.
  These mode-cuts correspond to removing about 200 modes from full eigenmodes.
  As shown in Fig.\ref{fig:PolyakovScatterDeconfined}(a),
  the Polyakov loop has a non-zero expectation value
  $\langle L_P \rangle \neq 0$,
  which shows the center group $Z_3$ structure on the complex plane.
  This property indicates the deconfined phase.

  After removing low-lying or UV Dirac-modes,
  as shown in Figs.\ref{fig:PolyakovScatterDeconfined}(b) and (c),
  the Dirac-mode projected Polyakov loop $\langle L_P \rangle_{\rm IR/UV}$
  still shows the non-zero value and the $Z_3$ center structure,
  which indicates the deconfined and $Z_3$ broken phase.
  In fact, in both cut cases of IR and UV Dirac modes,
  no drastic change occurs on the Polyakov loop,
  apart from a constant normalization factor.
  The Dirac-mode seems to be insensitive also for deconfinement properties of
  the Polyakov loop.

  \begin{figure}
    \centering
    \includegraphics[width=0.27\textwidth,clip]{scatter-beta6.0-no-cut.eps}
    \includegraphics[width=0.27\textwidth,clip]{scatter-beta6.0-ir-cut.eps}
    \includegraphics[width=0.27\textwidth,clip]{scatter-beta6.0-uv-cut.eps}
\vspace{-0.2cm}
    \caption{\label{fig:PolyakovScatterDeconfined}
    The scatter plot of the Polyakov loop in the deconfined phase
    in SU(3) lattice QCD on $6^3 \times 4$ at $\beta$=6.0, i.e.,
    $a\simeq$0.10fm and $T \equiv 1/(N_ta) \simeq$ 0.5GeV.
    (a) The original Polyakov loop $\langle L_P \rangle$.
    (b) $\langle L_P \rangle_{\rm IR}$ for low-lying Dirac-mode cut
        with $\Lambda_{\rm IR} = 0.5a^{-1}$.
    (c) $\langle L_P \rangle_{\rm UV}$ for UV Dirac-mode cut
        with $\Lambda_{\rm UV} = 2.0a^{-1}$.
    }
  \end{figure}

\subsection{$\beta$-dependence of Dirac-mode projected Polyakov loop}
  We also investigate the $\beta$-dependence of the absolute value
  of the Dirac-mode projected Polyakov loop,
  $|\langle L_{\rm P} \rangle_{\rm IR/UV}|$,
  at fixed $N_t$ and $L$.
  Here, we use $6^3 \times 4$ lattice with $\beta = 5.4 \sim 6.0$.

  Figure~\ref{fig:PolyakovDiracCut} (a) shows
  $|\langle L_{\rm P} \rangle_{\rm IR}|$
  with $\Lambda_{\rm IR} = 0.5a^{-1}$ and $1.0a^{-1}$, and
  Fig.\ref{fig:PolyakovDiracCut} (b) shows
  $|\langle L_{\rm P} \rangle_{\rm UV}|$
  with $\Lambda_{\rm UV} = 2.0a^{-1}$ and $1.7a^{-1}$.
  In terms of the removed number of Dirac modes,
  $\Lambda_{\rm IR} = 0.5a^{-1}$ and $1.0a^{-1}$
  approximately correspond to
  $\Lambda_{\rm UV} = 2.0a^{-1}$ and $1.7a^{-1}$, respectively.
  We have also added the original (no Dirac-mode cut) Polyakov loop
  $|\langle L_{\rm P} \rangle|$,
  which shows the phase transition around $\beta = 5.6 \sim 5.7$.
  Both IR and UV Dirac-mode projected Polyakov loop
  $|\langle L_{\rm P} \rangle_{\rm IR/UV}|$
  show the similar $\beta$-dependence as the original Polyakov loop
 $|\langle L_{\rm P} \rangle|$,
  apart from a constant normalization factor.

\begin{figure}
  \begin{center}
    \includegraphics[width=0.40\textwidth,clip]{polyakov-dirac-06-04-ir-cut.eps}
    \includegraphics[width=0.40\textwidth,clip]{polyakov-dirac-06-04-uv-cut.eps}
\vspace{-0.5cm}
  \end{center}
  \caption{ \label{fig:PolyakovDiracCut}
   The $\beta$-dependence of the Polyakov loop
   $|\langle L_{\rm P} \rangle_{\rm IR/UV}|$ in SU(3) lattice QCD on
   $6^3 \times 4$.
   (a) $|\langle L_{\rm P} \rangle_{\rm IR}|$ for IR Dirac-mode cut
       with $\Lambda_{\rm IR} = 0.5a^{-1}$ and $1.0a^{-1}$.
   (b) $|\langle L_{\rm P} \rangle_{\rm UV}|$ for UV Dirac-mode cut
       with $\Lambda_{\rm UV} = 2.0a^{-1}$ and $1.7a^{-1}$.
 The original Polyakov loop $|\langle L_{\rm P} \rangle|$ without cut is added.
  }
\end{figure}

\section{Summary and concluding remarks}
  In this paper, we have analyzed the direct relation between the Dirac eigenmodes
  and the Polyakov loop in SU(3) lattice QCD calculation at the quenched level.
  Using the Dirac-mode expansion method,
  we have carefully removed the relevant ingredient of chiral symmetry breaking
  from the Polyakov loop.

  In the confined phase, we have found that the Polyakov loop remains almost zero
  even without low-lying Dirac-mode.
  These low-lying modes are relevant for chiral symmetry breaking, as
  the Banks-Casher relation indicates.
  However, the Polyakov loop does not show any drastic changes,
  which indicates the system still remains in the confined phase.
  This result is consistent with the Wilson loop analysis,
  which shows the area law and the linear interquark potential
  even after removing low-lying Dirac-modes \cite{Suganuma:2011,Gongyo:2012}.
  We have also checked the UV Dirac-mode contribution to the Polyakov loop.
  By removing UV Dirac modes, the Polyakov loop remains almost zero.
  Thus, there seem to be no specific Dirac modes
  essential for the Polyakov loop in the confined phase.

  In addition to the confined phase,
  we have also analyzed the Polyakov loop properties in the deconfined phase
  at high temperature.
  In the deconfined phase, the Polyakov loop
  has a non-zero expectation value,
  which distributes in $Z_3$ direction in the complex plane.
  Even by removing low-lying or UV Dirac-modes,
  the behavior of the Polyakov loop $\langle L_P \rangle_{\rm IR/UV}$
  seems almost unchanged, apart from a constant normalization factor.

  These lattice QCD results suggest no direct connection
  between chiral symmetry breaking and color confinement
  through the Dirac eigenmodes, which indicates that
  one-to-one correspondence would not hold between them in QCD.
  If it is the case, the QCD phase diagram would exhibit more richer structure
  by mismatch of chiral and deconfinement phase transitions.

\section*{Acknowledgements}
  The lattice QCD calculations have been done on NEC-SX8 and NEC-SX9
  at Osaka University.
  This work is in part supported by a Grant-in-Aid for JSPS Fellows
  [No.23-752, 24-1458] and the Grant for Scientific Research
  [(C) No.23540306, Priority Areas ``New Hadrons'' (E01:21105006)]
  from the Ministry of Education, Culture, Science and Technology (MEXT)
  of Japan.


\section{#1}\setcounter{equation}{0}}
\newcommand{\newappendix}[1]{\section*{#1}\setcounter{equation}{0}}

\newcommand{\scs}{\scriptscriptstyle}
\newcommand{\f}{\frac}
\newcommand{\fm}[2]{{\textstyle \frac{#1}{#2}}}

\def\spins{\hbox{\tiny spins}}

\begin{document}
\title{Phase space analysis for  three and four massive particles in  final states}
\author{H.M.~Asatrian}
\author{A.~Hovhannisyan}
\author{A.~Yeghiazaryan}
\affiliation{Yerevan Physics Institute, 0036 Yerevan, Armenia}
\begin{abstract}
We propose  formulae for computing the phase space integrals of
$1\to 3$ and $1\to 4$ processes with  massive particles in final
states. As an application of these formulae we study the final state
mass effects in some interesting phenomenological cases, giving
fully integrated analytic results for the corresponding phase
spaces. We consider also the $B_s-\bar{B}_s$ process at NNLO and
calculate one of the most complicated master integrals, which
contributes to the $\Delta\Gamma_{B_s}$ at $O(\alpha_s^2)$.
\end{abstract}
\pacs{12.38.Bx, 13.20.He, 13.20.-v}
\maketitle

\newsection{Introduction \label{sec:intro}}
The $B_{s,d}$  -meson rare decays and oscillations are known to be a
unique source of indirect information about physics at scales of
several hundreds GeV. In the Standard Model all these processes
proceed through loop diagrams and thus are relatively suppressed. In
the extensions of the Standard Model the contributions stemming from
the diagrams with "new" particles in the loops can be comparable or
even larger than the contribution from the SM. Thus getting
experimental information about rare decays puts strong constraints
on the extensions of the SM or, if we are lucky, can lead to a
strong disagreement with the SM predictions providing an evidence of
some "new physics". To make a rigorous comparison between experiment
and theory, one has to get refined theoretical predictions for the
rare decay at hand. In particular,  perturbative QCD corrections
through next-to-next-to leading order (NNLO) in $\alpha_S$ are
needed. When calculating higher order QCD corrections, along with
the calculation of virtual corrections, it is necessary to take into
account real emission of gluons and quark-antiquark pairs. Then the
phase space integrals contain infrared and collinear singularities,
which can be regularized using dimensional regularization scheme.

Sometimes however it is more convenient to introduce a small mass as
an infrared regulator which can regularize both the soft and
collinear divergences. This can be for instance massive s quark for
$b\to s \gamma$ \cite{Kapustin:1995fk} and $b\to s \gamma \gamma$
\cite{Asatrian:2011ta}  or massive light quarks for $b\to
s\bar{q}q\gamma$ $(q=u, d, s)$ \cite{Kaminski:2012eb}. Another
example is the calculation of CP-asymmetry and width
difference-$\Delta\Gamma$ in $B_{s,d}-\bar{B}_{s,d}$ mixing at NLO
and NNLO in $\alpha_s$, when the calculation gets simplified by the
usage of a small gluon mass as an infrared regulator.

In this article we derive formulae for the $1\to 3$
\cite{Scharf:1993ds} and $1\to 4$ phase spaces, where the final
state particles are massive. Following the approach of
\cite{Asatrian:2002va} for the $1\to 3$ process  we derive formulae
for the phase space in the rest frame of two particles in the final
state and obtain a factorized parameterization of the  phase-space
for the case when all three particles in final state are massive. In
case of $1\to 4$ process we consider similar method, deriving the
phase space formula in the rest frame of three particles in the
final state. Here we come to the factorized parameterization of
phase-space  when only one particle in the final state is massive.

The paper is organized as follows. In Section II we discuss
three-particle phase space. Section III is devoted to the
investigation of  four-particle phase space. In Section IV we apply
our technique to the calculation of master integral connected with
the $B-\bar{B}$ mixing at $O(\alpha_s^2)$. The Appendix contains
formulae of $1\to 3$ and $1\to 4$ integrated phase spaces for some
particular cases.

\newsection{Three-Particle Phase Space}

Here we consider the three particle decay, when particle with
momentum $p^{'}$ decays into three particles with momenta $p^{'}_i$
and masses $m_i$, $i=1,2,3$; $p^{'}=p^{'}_1+p^{'}_2+p^{'}_3$. We
start from the well-known expression for the differential decay
width:
\begin{eqnarray}
\label{phN} \nonumber d\Gamma=\frac{1}{2m} \overline{|M|^2} D\Phi \,
,
\end{eqnarray}
where $\overline{|M|^2}$ is the squared matrix element, summed and
averaged over spins and colors of the particles in the final and
initial states respectively and $m$ is the mass of decaying
particle.

The phase space formula in the rest frame of particles 1 and 3 with
momenta $p_1^{'}$ and $p_3^{'}$ was derived in
\cite{Asatrian:2002va}:
\begin{eqnarray}
\label{phasefactorsleptrest}
&& \hspace{-0.3cm} D\Phi^{}=m^{2(d-3)} D\Phi^{}_1\, D\Phi_2 \, ds_{13}\, ,\\
&&\nonumber \hspace{-0.3cm} D\Phi^{}_1= (2\pi)^d \,
s_{13}^{(d-2)/2}\, \frac{d^{d-1} p}{2 p^0}
\frac{d^{d-1}p_2}{(2\pi)^{d-1}2p_2^0} \, \delta^d(p-p_2-q),\\
\nonumber && \hspace{-0.3cm} D\Phi_2=
\frac{d^{d-1}p_1}{(2\pi)^{d-1}2p_1^0} \,
\frac{d^{d-1}p_3}{(2\pi)^{d-1}2p_3^0} \, \delta^d(q-p_1-p_3) \, .
\end{eqnarray}
where $d=4-2\epsilon$ and   "dimensionless" momenta $p=p^{'}/m$,
$p_i=p_i^{'}/m$, i=1,2,3 are introduced. In
(\ref{phasefactorsleptrest}) we introduced an additional integration
over $s_{13}=q^{2}$, $q=p_1+p_3$ (in the considered frame
$q=(\sqrt{s_{13}},\vec{0})$ ). Now we can integrate over  $d-1$
components of $p_1$ and $p_2$ using the spatial parts of  two
$d$-dimensional $\delta$ functions. To carry out the remaining
integrations we choose the coordinate axes in a way that  particle
momenta have the following components:
\begin{eqnarray} \label{vectorform}
p &=& (E,|\vec{p}|,0,0,....) \, , \nonumber \\
p_3 &=& (E_3,|\vec{p}_3| \cos \vartheta, |\vec{p}_3| \sin
\vartheta,0,....)\, ,
\end{eqnarray}
where the dots correspond to the components of extra space
dimensions, which are all zero. Making use of the remaining two
one-dimensional $\delta$-functions we  express $E$ and $E_3$ in the
following way:
\begin{eqnarray}
E = \frac{1+s_{13}-x_2}{2\sqrt{s_{13}}} \quad ; \quad E_3 =
\frac{x_3+s_{13}-x_1}{2\sqrt{s_{13}}} \, ,
\end{eqnarray}
where $x_i=m^2_i/m^2$, $i=1,2,3$ are dimensionless variables. After
integration over those angles on which $\overline{|M|^2}$ does not
depend, we obtain a factorized formula for the three massive
particles phase space (unlike \cite{Scharf:1993ds}, where the
corresponding formula is not factorized):
\begin{widetext}
\begin{eqnarray}
\label{PS3} \nonumber D\Phi(1 \to 3) &=& m^{2(d-3)}
\frac{\Omega_{d-1} \, \Omega_{d-2} }{2^{4} \, (2\pi)^{2d-3} }
\,s_{13}^{(d-4)/2} \, \left(|\vec{p}| \, |\vec{p}_3| \,
\sin\vartheta \right)^{d-3} \, d\vartheta \, ds_{13}
\\
&=& m^{2(d-3)} \frac{\Omega_{d-1} \, \Omega_{d-2} }{2^{4} \,
(2\pi)^{2d-3} } \,s_{13}^{(d-4)/2} \, \left(|\vec{p}| \, |\vec{p}_3|
\, \right)^{d-3} (1-\cos^2\vartheta)^{\frac{d-4}{2}} \,
d\cos\vartheta \, ds_{13}
\end{eqnarray}
where $\Omega_d$ is the solid angle in $d$ dimensions,
\end{widetext}
\begin{eqnarray}
\Omega_d=\frac{2\pi^{d/2}}{\Gamma(d/2)}\,,
\end{eqnarray}
and
\begin{eqnarray}
\nonumber && |\vec{p}| =
\frac{1}{2\sqrt{s_{13}}}\sqrt{(1+x_2-s_{13})^2-4x_2}\, ,
\\
\quad && |\vec{p}_3| =
\frac{1}{2\sqrt{s_{13}}}\sqrt{(x_1+x_3-s_{13})^2-4x_1x_3} \,
\end{eqnarray}
are 3-momenta expressed by $K\ddot{a}ll\acute{e}n$ functions.

For the integration limits we have
\begin{eqnarray}
&& \nonumber (\sqrt{x_1}+\sqrt{x_3}\,)^2\leq s_{13}\leq
(1-\sqrt{x_2}\,)^2 ~,
\\
&& \hspace{1.5cm} -1\leq \cos\vartheta \leq1.
\end{eqnarray}

We represent all  scalar products via variables $s_{13}$ and
$\cos\vartheta$:
\begin{eqnarray}
&& \nonumber p.p_2=\frac{1+x_2-s_{13}}{2}
\\
&& \nonumber p.p_1=\frac{2(1+x_1+x_3)-x_2-X}{4}
\end{eqnarray}
\begin{eqnarray}
&& \nonumber p.p_3=\frac{-2(x_1+x_3)-x_2+2s_{13}+X}{4}
\\
&& \nonumber p_2.p_1=\frac{2+4x_3-x_2-2s_{13}-X}{4}
\\
&& \nonumber p_2.p_3=\frac{-x_2-4x_3+X}{4}
\\
&& p_1.p_3=\frac{s_{13}-x_1-x_3}{2}
\end{eqnarray}
where
\begin{eqnarray}
\nonumber X&=& -4|\vec{p}| \,
|\vec{p}_3|\cos\vartheta+1+2(x_1+x_3)-s_{13}
\\ &&-\frac{(1+s_{13}-x_2)(x_1-x_3)}{s_{13}}\,,
\end{eqnarray}
and $|\vec{p}|$ and $|\vec{p}_3|$ are defined in (2.6).

From the formula (\ref{PS3}) for the three-particle phase space one
can get a simpler formula assuming that  the masses of two particles
in the  final state are equal to zero: $m_{2}=m_{3}=0$. In that case
we make the replacements: $\cos\vartheta=2\lambda_1-1$ and
$s_{13}=\lambda_2 (1-x_1)+x_1$, where $0<\lambda_{1,2}<1$. Then we
get:
\begin{widetext}
\begin{eqnarray}
\label{PS3N} D\Phi(1 \to 3) &=& m^{2(d-3)} \frac{\Omega_{d-1} \,
\Omega_{d-2} }{2^{d+1} \, (2\pi)^{2d-3} }\, d \lambda_1 \,
d\lambda_2
\\
&& \nonumber \times \, (1-x_1)^{2d-5}\,\left [\lambda_2 (1-x_1)+x_1
\right ]^{(2-d)/2} \, \left[\lambda_2\, (1-\lambda_2)\right]^{d-3}\,
\left[\lambda_1(1-\lambda_1)\right]^{\frac{d-4}{2}}
\end{eqnarray}
\end{widetext}
If we put $x_1=0$ (i.e. $m_1=0$) in the expression (\ref{PS3N}), we
get the formula (2) of \cite{Anastasiou:2003gr}.

In the Appendix A we study the final state mass effects in some
interesting phenomenological cases, providing fully integrated
analytic formulae for the $1\to 3$ processes.

\newsection{Four-Particle Phase Space}

Now we proceed to the case of four particle decay. For this case one
can derive the expression for the phase space  in the similar way as
in \cite{Asatrian:2002va}. The only difference is that instead of
two particle's rest frame introduced in \cite{Asatrian:2002va} here
the rest frame of three particles is used (see \cite{GehrmannDe
Ridder:2003bm}). We consider a particle with momentum $p^{'}$, which
decays into four particles with momenta $p_i^{'}$ and masses $m_i$,
$i=1,...,4$ and $p^{'}=p_1^{'}+p_2^{'}+p_3^{'}+p_4^{'}$. Here also
we introduce "dimensionless" momenta $p=p^{'}/m$, $p_i=p_i^{'}/m$,
i=1,...,4. In the rest frame of particles 2, 3 and 4 with momenta
$p_2^{'}$, $p_3^{'}$, $p_4^{'}$ (where
$\vec{p}_2+\vec{p}_3+\vec{p}_4=0$) four-particle phase space can be
written in the following way
\begin{eqnarray}
&&d\Phi^{}=m^{3d-8}(2\pi)^dd\Phi^{}_1\, d\Phi_2 \, ds_{234}\, ,\\
&&\nonumber d\Phi^{}_1=(s_{234})^{(d-2)/2} \, \frac{d^{d-1} p}{2
p^0} \, \frac{d^{d-1}p_1}{(2\pi)^{d-1}2p_1^0} \, \delta^d(p-p_1-q)\,
,
\\
\nonumber &&D\Phi_2=\frac{d^{d-1}p_2}{(2\pi)^{d-1}2p_2^0} \,
\frac{d^{d-1}p_3}{(2\pi)^{d-1}2p_3^0} \,
\frac{d^{d-1}p_4}{(2\pi)^{d-1}2p_4^0} \,
\\ \nonumber && \hspace{1cm} \times \delta^d(q-p_2-p_3-p_4)
\, ,
\end{eqnarray}
where $s_{234}$=$(p_2+p_3+p_4)^2$ and $q$ is a d-dimensional vector
with components $(\sqrt{s_{234}},\vec{0})$ in the rest frame of
particles 2, 3 and 4. We choose the coordinate axes in such a way,
that the particle momenta have the following components:
\begin{eqnarray} \label{vectorform}
&& \hspace{-0.45cm} p_2 = (E_2,|\vec{p}_2|,0,0,....) \, , \nonumber
\\
&& \hspace{-0.45cm} p_1 = (E_1,|\vec{p}_1| \cos \vartheta_1,
|\vec{p}_1| \sin \vartheta_1,0,....)\, ,
\\
&& \hspace{-0.45cm} p_3 = (E_3,|\vec{p}_3| \cos \vartheta_3,
|\vec{p}_3| \sin \vartheta_3\cos \phi_3,|\vec{p}_3| \sin
\vartheta_3\sin \phi_3,....) , \nonumber
\end{eqnarray}
where the dots, as before, correspond to the components of extra
space dimensions, which are all zero. We can immediately integrate
over those angles on which none of the scalar products depends and
obtain
\begin{eqnarray}
&& d^{d-1}p_2= \Omega_{d-1}\,|\vec{p}_2|^{d-2}~d|\vec{p}_2|\,,
\\ \nonumber
&& d^{d-1} p_1=
\Omega_{d-2}\,(\sin\vartheta_1)^{d-4}~d\cos\vartheta_1~|\vec{p}_1|^{d-2}~d|\vec{p}_1|\,,~~~~~
\\
&& \nonumber d^{d-1}p_3=
\Omega_{d-3}\,(\sin\phi_3)^{d-5}~d\cos\phi_3~(\sin\vartheta_3)^{d-4}
\\ \nonumber && \hspace{1.3cm} \times d\cos\vartheta_3~|\vec{p}_3|^{d-2}~d|\vec{p}_3|\,.
\end{eqnarray}
$p_4$ and $p$ can be integrated out immediately using the spatial
parts of the $\delta$ functions in (3.1). Thus for the phase space
we get
\begin{widetext}
\begin{eqnarray}
&& \hspace{-0.6cm} D\Phi_1=  \frac{(s_{234})^{(d-2)/2}}{2E} \,
\frac{\Omega_{d-2}(\sin\vartheta_1)^{d-4}~d\cos\vartheta_1~|\vec{p}_1|^{d-2}~d|\vec{p}_1|}{(2\pi)^{d-1}2E_1}
\, \delta\left(E-E_1-\sqrt{s_{234}}\right),
\\
&& \hspace{-0.6cm} D\Phi^{}_2=
\frac{\Omega_{d-1}|\vec{p}_2|^{d-2}d|\vec{p}_2|}{(2\pi)^{d-1}2 E_2}
\, \frac{\Omega_{d-3}
(\sin\phi_3)^{d-5}d\cos\phi_3(\sin\vartheta_3)^{d-4}d\cos\vartheta_3
|\vec{p}_3|^{d-2}d|\vec{p}_3|}{(2\pi)^{d-1}2E_3}
\frac{\delta\left(\sqrt{s_{234}}-E_2-E_3-E_4\right)\,}{(2\pi)^{d-1}2E_4},
\,
\end{eqnarray}
where
\end{widetext}
\begin{eqnarray}
&&E_4=
\sqrt{|\vec{p}_2|^2+|\vec{p}_3|^2+2|\vec{p}_2||\vec{p}_3|\cos\vartheta_3+x_4}\,,~~~~~
\\
&&E=\sqrt{|\vec{p}_1|^2+1}\, ,
\\&& x_i=m_i^2/m^2, i=1,...,4   \nonumber
\end{eqnarray}
where $|\vec{p_i}|=\sqrt{E_i^2-x_i}$, $i=1,2,3$. Making use of the
remaining two one-dimensional $\delta$-functions, we can express
$E_1$ and $E_3$ in terms of the other variables. Using the $\delta$
function in $d\Phi_1$ and (3.7) we get:
\begin{eqnarray}\label{E1E3}
E_1 = \frac{1-x_1-s_{234}}{2\sqrt{s_{234}}} \,\, , \quad E =
\frac{1-x_1+s_{234}}{2\sqrt{s_{234}}} \, .
\end{eqnarray}
From the $\delta$ function in $d\Phi_2$ together with (3.6) we find
the expression for $E_3$
\begin{widetext}
\begin{eqnarray}\label{theta1}
\nonumber\ E_3&=&\frac{1}
   {-2\left(\sqrt{s_{234}}-E_2\right)^2+2 |\vec{p_2}|^2 \cos ^2 \vartheta_3}\,\left[\left(E_2-\sqrt{s_{234}}\right)
   \left(s_{234}-2 E_2 \sqrt{s_{234}}+x_2+x_3-x_4\right)\right.
   \\
  &+& \left.|\vec{p}_2| \cos \vartheta_3 \sqrt{\left(s_{234}-2 E_2
\sqrt{s_{234}}+x_2-x_3-x_4\right)^2-4 x_3
\left(x_4+|\vec{p}_2|^2\sin ^2 \vartheta_3)\right)}\right].
\end{eqnarray}
Finally, introducing new variables: $z_1=\cos\vartheta_1$,
$z_{31}=\cos\vartheta_3$, $z_{32}=\cos\phi_3$, we get for the
four-particle phase space:
\begin{eqnarray}\label{PS4}
\nonumber D\Phi(1\to 4)&=&\frac{m^{3d-8}}{32} \,
\frac{\Omega_{d-1}\Omega_{d-2}\Omega_{d-3}}{(2\pi)^{3d-4}} \,dz_1
~dz_{31}~dz_{32}~dE_2~ds_{234}
\\
&\times&
(1-z_1^2)^{\frac{(d-4)}{2}}~(1-z_{31}^2)^{\frac{(d-4)}{2}}~(1-z_{32}^2)^{\frac{(d-5)}{2}}
{\frac{
W^{d-3}}{\sqrt{s_{234}}-E_2+{\frac{E_3}{|\vec{p_3}|}|\vec{p_2}|~z_{31}}}}
\end{eqnarray}
\end{widetext}
where
\begin{eqnarray}
W=|\vec{p_1}|~|\vec{p_2}|~ |\vec{p_3}|~ \sqrt{s_{234}}\,.
\end{eqnarray}
The integration variables  $s_{234},E_2,z_1,z_{31},z_{32}$ in
(\ref{PS4})
 have the following limits
\begin{eqnarray}
(\sqrt{x_2}+\sqrt{x_3}+\sqrt{x_4})^2\leq s_{234}\leq
(1-\sqrt{x_1}\,)^2\,,
\end{eqnarray}
\begin{eqnarray}
\nonumber &&\sqrt{x_2}\leq E_2\leq \frac{s_{234}+x_2-(\sqrt{x_3}+
\sqrt{x_4})^2}{2\sqrt{s_{234}}}\,,
\\
&& \hspace{2cm} -1\leq z_i\leq1\,.
\end{eqnarray}

Different possible scalar products between 4-momenta of the
particles can be obtained automatically from the expressions for
their components (3.2) using  the conservation of momentum. All the
scalar products can be expressed through the 5 integration variables
in (\ref{PS4}):
\begin{eqnarray} && \nonumber
p_1.p_2=E_1E_2-|\vec{p_1}||\vec{p_2}|\cos{\vartheta_1}\,,
\\
&& p_2.p_3=E_2E_3-|\vec{p_2}||\vec{p_3}|\cos{\vartheta_3}\,,
\\
&& \nonumber p_1.p_3 =
E_1E_3-|\vec{p_1}||\vec{p_3}|(\cos{\vartheta_1}\cos{\vartheta_3}
\\ \nonumber && \hspace{1cm} +\sin{\vartheta_1}\sin{\vartheta_3}\cos{\phi_3})\,,
\end{eqnarray}
and
\begin{eqnarray}
&&\nonumber \hspace{-1cm} p_1.p_4 =
E_1\sqrt{s_{234}}-p_1.p_3-p_1.p_2\,,
\\
&&\nonumber \hspace{-1cm} p_2.p_4=E_2\sqrt{s_{234}}-x_2-p_2.p_3 \,,
\\
&&\hspace{-1cm} p_3.p_4 = E_3\sqrt{s_{234}}-x_3-p_2.p_3\,,
\end{eqnarray}
In right-hand side of (3.15) the expressions given in (3.14) should
be used. Note that in (3.14) and (3.15) $E_1$ and $E_3$ are
determined by formulae (3.8) and (3.9).

We can also obtain a factorized formula for the 4 particle phase
space with 1 massive and 3 massless particles in the final state,
taking $x_{2}=x_{3}=x_{4}=0$. Making the following substitutions in
(\ref{PS4})
\begin{eqnarray}
\nonumber && s_{234}=(1-\sqrt{x_1})^2 \lambda_1,\,
E_2=\frac{1}{2}\sqrt{\lambda_1}(1-\lambda_2) (1-\sqrt{x_1}),
\\ &&
z_{31}=\frac{\lambda_2(1-\lambda_4)-\lambda_4}
{\lambda_2(1-\lambda_4)+\lambda_4},\, z_1=2\lambda_3-1,
\\
&& \nonumber z_{32}=2\lambda_5-1,
\end{eqnarray}
the differential phase space can be expressed in the following form:
\begin{widetext}
\begin{eqnarray}\label{PS4N}
\nonumber D\Phi(1\to 4)&=&\frac{m^{3d-8}}{128} \,
\frac{\Omega_{d-1}\Omega_{d-2}\Omega_{d-3}}{(2\pi)^{3d-4}}
\,d\lambda_1 d\lambda_2d\lambda_3d\lambda_4d\lambda_5
\\
&\times& \nonumber
(1-\sqrt{x_1})^{3d-7}\left[(1-\lambda_1)((1+\sqrt{x_1})^2-\lambda_1(1-\sqrt{x_1})^2)\right]^{(d-3)/2}
\\&\times&
\left[\lambda_1(1-\lambda_2)\right]^{d-3}
\left[\lambda_2(1-\lambda_3)\lambda_3(1-\lambda_4)\lambda_4\right]^{(d-4)/2}
\left[\lambda_5(1-\lambda_5)\right]^{(d-5)/2}.
\end{eqnarray}
\end{widetext}
where $\lambda_i=[0,1]$, $i=1...5$. The corresponding
scalar products in terms of the $\lambda_i$ will have the following
form
\begin{eqnarray}
\nonumber s_{234}&=&(p_2+p_3+p_4)^2=(1-\sqrt{x_1})^2 \lambda_1,
\\
\nonumber s_{34}&=&(p_3+p_4)^2=(1-\sqrt{x_1})^2 \lambda_1 \lambda_2,
\\
\nonumber s_{23}&=&(p_2+p_3)^2=(1-\sqrt{x_1})^2 \lambda_1
(1-\lambda_2) \lambda_4,
\end{eqnarray}
\begin{eqnarray}
\nonumber
s_{134}&=&(p_1+p_3+p_4)^2=(s_{134}^+-s_{134}^-)\lambda_3+s_{134}^-,
\\
s_{13}&=&(p_1+p_3)^2=(s_{13}^+-s_{13}^-)\lambda_5+s_{13}^-,
\end{eqnarray}
where
\begin{widetext}
\vspace{-0.2cm}
\begin{eqnarray}
\nonumber && \hspace{-0.4cm}
s_{134}^{\pm}=\frac{1}{2}\left[1+\lambda_2-\lambda_1(1-\lambda_2)(1-\sqrt{x_1})^2+x_1(1-\lambda_2)
\pm
(1-\sqrt{x_1})(1-\lambda_2)\sqrt{(1-\lambda_1)((1+\sqrt{x_1})^2-\lambda_1(1-\sqrt{x_1})^2)}\right],
\\
\nonumber && \hspace{-0.4cm}
s_{13}^{\pm}=x_1+\frac{1}{2}\,(1-\sqrt{x_1})\left[\left(\lambda_2\left(1-\lambda_4\right)+\lambda_4\right)
\left(1+\sqrt{x_1}-\lambda_1(1-\sqrt{x_1})\right)\right.
\\
\nonumber && \hspace{0.2cm}
+\left.\left(\lambda_2(1-\lambda_4)-\lambda_4\right)\left(1-2\lambda_3\right)\sqrt{\left(1-\lambda_1\right)
\left((1+\sqrt{x_1})^2-\lambda_1(1-\sqrt{x_1})^2\right)}\right]
\\
&& \hspace{0.2cm} \pm
2(1-\sqrt{x_1})\sqrt{\lambda_2(1-\lambda_3)\lambda_3(1-\lambda_4)\lambda_4(1-\lambda_1)
\left((1+\sqrt{x_1})^2-\lambda_1(1-\sqrt{x_1})^2\right)}\,.
\end{eqnarray}
\end{widetext}
As it is expected if we put $x_1=0$ in (\ref{PS4N}) we get the
formula (21) of \cite{Anastasiou:2003gr} derived for massless
particles in the final state.

In the Appendix B we consider the final state mass effects in some
interesting phenomenological cases, providing fully integrated
analytic formulae for the $1\to 4$ processes.

\newsection{Contribution to $B-\bar{B}$ mixing
at $\mathcal{O}(\alpha_s^2)$}

In this Section we illustrate how one can use formulae derived in
previous sections for calculation of  the absorptive part of a
$B-\bar{B}$ mixing at $\mathcal{O}(\alpha_s^2)$, which contributes
to the width differences $\Delta \Gamma$ and the CP-asymmetry.

An example of a diagram that contributes to the $B-\bar{B}$ mixing
at $O(\alpha_s^2)$ is shown in Fig.~1. This diagram contains both
3-particle and 4-particle cuts. The light quark masses we put equal
to zero.

We introduce a small mass $m_g$ for a gluon as an infrared
regulator. In this way, in the matching one doesn't need $\epsilon$
and $\epsilon^2$ parts of NLO and LO Wilson coefficients
correspondingly. Moreover, in the Effective Theory side the
renormalization is much easier, because the $(1/\epsilon^n)$ terms
appear only as UV singularities and so their cancelation can be
tracked easier. Further in order to eliminate IR finite terms coming
from the interference of IR $(1/\epsilon)$ with UV structures of
$O(\epsilon)$ a special class of operators $(\propto\epsilon)$ has
to be introduced in  case of dimensional regularization. Whereas
with introduction of a gluon mass those terms vanish automatically.
This is the major difference in using dimensional regularization or
gluon mass for IR singularities.

The reduction to the master integrals is done by means of a program
FIRE \cite{Smirnov:2008iw}.

Here we consider one of the most complicated master integrals (MI),
shown on the right side of the Fig.~1. This MI has one
three-particle cut and one four-particle cut. While in the case of
3-particle cut $i$ on the right side of the Fig.~1 denotes either
heavy $m_c$, $m_b$ or  light quarks, in the case of a 4-particle cut
$i$  stands for only  light quarks. The corresponding cuts can be
calculated analytically as an expansion over a gluon mass. The
coefficient of a corresponding MI contains inverse masses of $m_g$
up to the forth power.
\begin{figure}[t]\label{fig:PS_pic}
\vspace*{-4.8cm}
\begin{center}
\hspace*{-2.2cm}
\includegraphics[width=15cm,angle=0]{PS_pic.ps}
\vspace*{-13.5cm} \caption{Left: an example of an infrared divergent
diagram contributing to the $B-\bar{B}$ process at NNLO.
 The virtual gluons are represented by curly lines. Right: one of the
 master integrals appearing after reduction of the left diagram with two
 cuts A and B ($c$- is a charm, $g$- is a gluon and $i$- is either charm or
  beauty or a light quark).}
\end{center}
\vspace*{-3.1mm}
\end{figure}

The 3-particle cut (cut A) can be factorized into a product of a
$q\bar{q}$ fermion loop and a 3-particle phase space, which is
presented in the Appendix, formula (\ref{PS3mmml}). In the case of
4-particle cut we get:
\begin{widetext}
\begin{eqnarray}
r_B &=&n_l\frac{m_b^2}{98304 \pi ^5} \left[16\pi^2(2-3x_c)x_c+192
x_c^2 \text{Li}_2(\sigma )-96 (2-x_c) x_c \text{Li}_2\left(\sigma
^2\right)\right.
   \\
   && \nonumber \hspace{1.8cm}
   +12 \sqrt{1-4 x_c} (2 x_c+1) (8 \log (1-\sigma )-4 \log (\sigma )-5)
      \\
&& \nonumber +48 \log (\sigma ) \left(x_c^2 \log (\sigma )-4 (2-x_c)
x_c \log (1-\sigma )+(1-x_c)^2+2 (2-x_c) x_c \log (1-4 x_c)\right)
\\
   && \nonumber +x_g\cdot \left(-8 \pi ^2 (5-6 x_c)
   +96 (1-2 x_c) \text{Li}_2(\sigma )+48 (3-2 x_c) \text{Li}_2\left(\sigma ^2\right)\right.
   \\
   && \nonumber \hspace{1cm} +24 \log (\sigma ) (4 (3-2 x_c) \log (1-\sigma )
   +(1-2 x_c) \log (\sigma )-2 (3-2 x_c) \log (1-4 x_c))
\\
&& \left. \nonumber  \hspace{1cm} +8 \sqrt{1-4 x_c} \left(7\pi
^2-6-48 \log ^2(1-\sigma )-3 \log ^2(\sigma )+48 \log (1-\sigma )
\log (\sigma )\right) \right)
    \\
   && \nonumber + x_g^2\cdot\frac{4}{(1-4 x_c)^{3/2}} \left(-28 \pi ^2 x_c^2-3 \left(6 x_c^2-10
   x_c+5\right)\right.
\\
&&  \nonumber \hspace{2cm} +12 \log (\sigma ) \left(1-2 x_c-2
x_c^2-\sqrt{1-4 x_c} (1-x_c)+x_c^2 \log (\sigma )\right)
\\
&& \left. \nonumber \hspace{2cm} -24 \log (1-\sigma ) \left(1-2
x_c-2 x_c^2-8 x_c^2 \log (1-\sigma )+8 x_c^2 \log (\sigma )\right)
\right)
   \\
   && \nonumber +\log (x_g) \cdot \left(-24 \sqrt{1-4 x_c} (2 x_c+1)-96 x_c (1-x_c) \log (\sigma
   )\right.
   \\
   && \nonumber \hspace{2.2cm} +x_g \left(96 (1-x_c)\log (\sigma )+96 \sqrt{1-4 x_c} \left(2 \log (1-\sigma )-\log (\sigma )\right)\right)
   \\
   && \left. \nonumber \hspace{2.2cm} +x_g^2 \frac{24}{(1-4 x_c)^{3/2}} \left(1-2 x_c-2 x_c^2-16 x_c^2 \log (1-\sigma )+8 x_c^2 \log (\sigma )\right) \right)
   \\
  && \nonumber \left. + \log ^2(x_g) \cdot
\left(-x_g 24 \sqrt{1-4 x_c}+x_g^2\frac{48 x_c^2}{(1-4
x_c)^{3/2}}\right) \right]+\mathcal{O}\left(x_g^3\right),
\end{eqnarray}
where $n_l$ is the number of light quarks ($n_l=3$),
$x_c=(m_c/m_b)^2$, $x_g=(m_g/m_b)^2$ and $\sigma=\frac{1-\sqrt{1-4
x_c}}{1+\sqrt{1-4 x_c}}$.\\
\end{widetext}

In the expression (4.2) the infrared singularities  appear as
$\log^n(x_g)$, n=1,2. The infrared singularities of the diagram
shown on the left side of Fig.~1, together with the $\log^n(x_g)$-s
of the other diagrams get canceled in the matching with the
corresponding Effective Theory diagrams with $Q$, $\tilde{Q}_S$
insertions \cite{Lenz:2006hd,Beneke:1996gn}. It must be mentioned
also, that the diagrams both in the full and Effective Theory have
to be renormalized. In addition, in the Effective Theory
corresponding Evanescent operators have to be taken into account as
well \cite{Herrlich:1994kh,Beneke:1998sy}.

\newsection{Conclusions \label{sec:concl}}

To conclude, we have presented new formulae for computing $1\to 3$
and $1\to 4$ processes, with massive particles in the final states.
We have also obtained a factorized formula for the 4 particle phase
space with 1 massive and 3 massless particles in the final state. On
the example of $B-\bar{B}$ mixing we  demonstrate  the capability of
our technique for  NNLO calculations, considering one of the most
complicated master integrals. Further we  study the final state mass
effects in some interesting phenomenological cases providing fully
integrated analytic formulae for 3  and 4 particle phase spaces.

\newappendix{ACKNOWLEDGMENTS}

This work was supported by the State Committee of Science of Armenia
program number 11-1c014 and VolkswagenStiftung program number 86426.
A. H. was partially supported by the EU contract MRTN-CT-2006-035482
(Flavianet) in early stage of this work.

\newappendix{APPENDIX A: FULLY INTEGRATED 3-PARTICLE PHASE-SPACE \label{app:interm}}
\def\theequation{A.\arabic{equation}}

Here, as an application of the formula (\ref{PS3}) we present fully
integrated 3-particle phase space. For the case of $m_1=m_2$,
$m_3=m_{light}$ ($m_{light}\ll m_{1,2}$), expanding (\ref{PS3}) over
$m_{light}/m$ and $\epsilon$ and performing integration,
 we get in the $\overline{MS}$  scheme
\begin{widetext}
\begin{eqnarray}\label{PS3mmml}
&& PS_3(m_1,m_1,m_{light})=
\left(\frac{m^2}{\mu^2}\right)^{-2\epsilon}
\frac{m^2}{128\pi^3}\left[ \frac{1}{2}\sqrt{1-4x_1}(1+2x_1)+2
x_1\left(1-x_1\right) \log (\sigma)\right.
   \\
   && \nonumber \hspace{1.4cm} +x_{\ell} \left(-2 \left(1-x_1\right) \log (\sigma)+\sqrt{1-4x_1}
   \left(\log (x_{\ell})-4 \log (1-\sigma)+2 \log (\sigma)\right)\right)
   \\
&& \nonumber \hspace{1.4cm} +x_{\ell}^2 \frac{1}{2(1-4x_1)^{3/2}}
\left(-1+2 x_1+2 x_1^2-4 x_1^2 \left(\log (x_{\ell})-4 \log
(1-\sigma)+2 \log (\sigma)\right)\right)
     \\
&& \nonumber +~ \epsilon\cdot \left[ -\frac{4 \pi ^2}{3} (1-x_1)
x_1+4 (1-x_1) x_1 \left(\text{Li}_2(\sigma )+\text{Li}_2\left(\sigma
^2\right)\right)\right.
\\
&& \nonumber \hspace{1cm}+(1+2 x_1) \sqrt{1-4 x_1}
\left(\frac{13}{4}-\frac{1}{2} \log (1-4 x_1)-2 \log (1-\sigma
)+\log (\sigma )\right)
\\
&& \nonumber \hspace{1cm} - (1-x_1) \log (\sigma ) \left(1-5 x_1+6
x_1 \log (1-4 x_1)-12 x_1 \log (1-\sigma )+x_1 \log (\sigma )\right)
   \\
   && \nonumber +x_{\ell} \left(\frac{4 \pi ^2}{3} (1-x_1) -4 (1-x_1)
   \left(\text{Li}_2(\sigma )+\text{Li}_2\left(\sigma ^2\right)\right)\right.
    \\
    && \hspace{1cm} \nonumber +(1-x_1) \log (\sigma )
    \left(-6+6 \log (1-4 x_1)-12 \log (1-\sigma )+\log (\sigma
    )\right)
    \\
    && \nonumber \hspace{1cm} +\sqrt{1-4 x_1} \left(
   4 \log (1-\sigma ) (2 \log (1-\sigma )-2 \log (\sigma )-3)+(6-\log (\sigma )) \log
   (\sigma )\right.
   \\
   && \left. \left. \nonumber \hspace{1cm} +\log (1-4 x_1) (4 \log (1-\sigma )-2 \log (\sigma )-\log (x_{\ell}))-\frac{1}{2} \log
   ^2(x_{\ell})+3 \log (x_{\ell})-2 \pi ^2+2\right)\right)
   \\
   && \nonumber +x_{\ell}^2\,
   \frac{1}{(1-4x_1)^{3/2}}\left(4 \pi ^2 x_1^2+\frac{1}{4} \left(22 x_1^2-14 x_1+7\right)
   \right.
   \\
   && \nonumber \hspace{1cm} +\log (x_{\ell}) \left(x_1^2 \log (x_{\ell})+\frac{1}{2} \left(-6 x_1^2+2 x_1-1\right)+2 x_1^2 \log (1-4 x_1)\right)
   \\
   && \nonumber \hspace{1cm} -2 \log (\sigma ) \left(-x_1^2 \log (\sigma )+2 x_1^2-2 x_1+1-\sqrt{1-4 x_1} (1-x_1)\right)
\\
&& \nonumber \hspace{1cm} + \frac{1}{2} \log (1-4 x_1) \left(-16
x_1^2 \log (1-\sigma )+8 x_1^2 \log (\sigma )-2 x_1^2-2 x_1+1\right)
   \\
   && \nonumber \left.\left.\left. \hspace{1cm} +\log (1-\sigma )
   \left(4 \left(2 x_1^2-2 x_1+1\right)-16 x_1^2 \log (1-\sigma )
   +16 x_1^2 \log (\sigma )\right)\right)\right] \right]
   +\mathcal{O}\left(\epsilon^2,x_{\ell}^3\right)\,.
\end{eqnarray}
For the case $m_1>m_3$ and $m_2=0$ we have
\begin{eqnarray}
&& PS_3(m_1,0,m_3)=\frac{m^2}{128\pi^3}\left[ \frac{1}{2}
(1+x_1+x_3) \sqrt{1-2 (x_1+x_3)+(x_1-x_3)^2}\right.
\\
&& \nonumber \hspace{1cm} -2 (x_1+x_3-2 x_1 x_3) \coth
^{-1}\left(\sqrt{\frac{1-\left(\sqrt{x_1}-\sqrt{x_3}\right)^2}
{1-\left(\sqrt{x_1}+\sqrt{x_3}\right)^2}}\right)
\\
&& \nonumber +(x_1-x_3) \left(2 \log
\left(\sqrt{1-\left(\sqrt{x_1}+\sqrt{x_3}\right)^2}
\left(\sqrt{x_1}-\sqrt{x_3}\right)+\sqrt{1-\left(\sqrt{x_1}-\sqrt{x_3}\right)^2}
   \left(\sqrt{x_1}+\sqrt{x_3}\right)\right)\right.
   \\
   && \left.\left. \nonumber \hspace{3cm} -\frac{1}{2} \log (16 x_1 x_3)\right)\right]
   +\mathcal{O}\left(\epsilon\right)\,.
\end{eqnarray}
For $m_1\gg m_3=m_{light}$, expanding over $m_{light}/m$, we get
also the terms proportional to $\epsilon^1$
\begin{eqnarray}
&& PS_3(m_1,0,m_{light})=
\left(\frac{m^2}{\mu^2}\right)^{-2\epsilon}
\frac{m^2}{128\pi^3}\left[ \frac{1}{2}
\left(1-x_1^2\right)+x_1\log(x_1) \right.
\\
&&\nonumber
+x_{\ell}\,((1-x_1)\log(x_{\ell})-2(1-x_1)\log(1-x_1)-x_1\log (x_1))
+x_{\ell}^2\left(\frac{1}{2}-\frac{1}{1-x_1}\right)
   \\
&& \nonumber +\epsilon\cdot\left[\frac{13}{4}
\left(1-x_1^2\right)-\frac{\pi^2x_1}{3}+\frac{1}{2}x_1(6-x_1-\log
(x_1))\log(x_1)\right.
\\
\nonumber && \hspace{1cm} -2\log(1-x_1)\left(1-x_1^2+x_1
\log(x_1)\right) -2x_1(\text{Li}_2(1-x_1)-\text{Li}_2(x_1))
\\
\nonumber && +x_{\ell} \left( (1-x_1) \left(6 \log ^2(1-x_1)-6 \log
(1-x_1)+3 \log (x_{\ell})+2\right)\right.
\\
&& \nonumber \hspace{1cm} -\frac{1}{2} (1-x_1) \left(4 \log (1-x_1)
\log (x_{\ell})+\log ^2(x_{\ell})\right) -\pi ^2 \left(1-\frac{4
x_1}{3}\right)
\\
&& \left. \nonumber \hspace{1cm} -x_1 \left(-2 \log
(1-x_1)-\frac{\log (x_1)}{2}+3\right) \log (x_1)+2
(\text{Li}_2(1-x_1)-x_1 \text{Li}_2(x_1))\right)
\\
&& \nonumber \left.\left. +x_{\ell}^2 \left(\frac{1+x_1}{1-x_1}
\left(3 \log (1-x_1)-\frac{\log
(x_{\ell})}{2}+\frac{7}{4}\right)-\frac{x_1 \log
(x_1)}{1-x_1}\right)\right]\right]+\mathcal{O}\left(\epsilon^2,x_{\ell}^3\right)\,,
\end{eqnarray}
where the following notations were introduced
\begin{eqnarray}\label{sigma}
x_{\ell}=(m_{light}/m)^2\,,\hspace{0.5cm}
x_{1,3}=(m_{1,3}/m)^2\,,\hspace{0.5cm} \sigma=\frac{1-\sqrt{1-4
x_1}}{1+\sqrt{1-4 x_1}}\,.
\end{eqnarray}
\end{widetext}

\newappendix{APPENDIX B: FULLY INTEGRATED 4-PARTICLE PHASE-SPACE}
\def\theequation{B.\arabic{equation}}

As an application of the formula (\ref{PS4}) we give fully
integrated 4-particle phase space. For the case when $m_1=m_2$,
$m_3=m_{light}$, $m_4=0$ ($m_{light}\ll m_{1,2}$), expanding
(\ref{PS4})  over $m_{light}/m$ and $\epsilon$ and performing
integration we get
\begin{widetext}
\begin{eqnarray}
&& \nonumber PS_4(m_1,m_1,m_{light},0)=\frac{m^4}{8192 \pi ^5}
\left[-\frac{4 \pi ^2 x_1^2}{3} +\frac{1}{3} \sqrt{1-4 x_1}
\left(1+20x_1+12 x_1^2\right)+4 x_1 \left(1+x_1-2 x_1^2\right) \log
(\sigma )\right.
\\
&& \hspace{1.3cm} +4 x_1^2 \log (\sigma ) (2 \log (1-4x_1)-4 \log
(1-\sigma )+\log (\sigma ))- 16 x_1^2 \text{Li}_2(-\sigma)
\\
&& \nonumber \hspace{0.5cm} + x_{\ell}\cdot \left(-16 x_1^2
\text{Li}_2(\sigma )+8x_1 (2-x_1) \text{Li}_2\left(\sigma
^2\right)+4 \left(x_1^2-1+\sqrt{1-4 x_1} (2x_1+1)\right) \log
(\sigma )\right.
\\
   && \nonumber\hspace{1.3cm} -4 x_1^2 \log ^2(\sigma )-\log (1-\sigma ) \left(8 \sqrt{1-4 x_1} (2 x_1+1)-16 (2-x_1) x_1
   \log (\sigma )\right)
\\
   && \nonumber \hspace{1.3cm}-8 (2-x_1) x_1 \log (\sigma ) \log (1-4
   x_1)+\log (x_{\ell}) \left(8 (1-x_1) x_1 \log (\sigma )+2 \sqrt{1-4
   x_1} (2 x_1+1)\right)
\\
   && \nonumber \left. \hspace{1.3cm}-\frac{4}{3} \pi ^2 (2-3 x_1) x_1+3 \sqrt{1-4 x_1}
   (2 x_1+1)\right)
\\
&& \nonumber \hspace{0.5cm} \left. + x_{\ell}^2\cdot \left(-4
(1-x_1) \log (\sigma )-\sqrt{1-4 x_1}\, (3+8 \log (1-\sigma )-4 \log
(\sigma)-2 \log
   (x_{\ell}))\right)\right]+\mathcal{O}\left(\epsilon,x_{\ell}^3\right)\,.
\end{eqnarray}
For $m_1>m_3$, $m_2=m_4=0$ we have
\begin{eqnarray} \label{ff2}
&& \nonumber PS_4(m_1,0,m_3,0)=\frac{m^4}{8192 \pi ^5} \left[4
\left(x_1^2-x_3^2\right) \tanh
   ^{-1}\left(\sqrt{\frac{1-\left(\sqrt{x_1}+\sqrt{x_3}\right)^2}
   {1-\left(\sqrt{x_1}-\sqrt{x_3}\right)^2}}
   \frac{\sqrt{x_1}-\sqrt{x_3}}{\sqrt{x_1}+\sqrt{x_3}}\right)\right.
\\
&& \hspace{1.0cm} +\frac{1}{3} \sqrt{1-2 (x_1+x_3)+(x_1-x_3)^2}
\left(1+x_1^2+x_3^2+10 (x_1+x_3+x_1 x_3)\right)
   \\
   && \nonumber \hspace{1.0cm}+2 (x_1 - x_3) \log (4 x_1)+\left(x_1^2 (1-2 x_3)+2 x_3 (1-x_1 x_3)+x_3^2+4 x_1 x_3 \log (x_3)\right)
      \\
   && \nonumber \hspace{1.2cm} \times \left(\log (4 x_1 x_3)-2 \log \left(1-x_1-x_3+\sqrt{1-2
(x_1+x_3)+(x_1-x_3)^2}\right)\right)
   \\
   && \nonumber \hspace{1.0cm}-4 (x_1-x_3) \log \left(1+x_1-x_3+\sqrt{1-2
(x_1+x_3)+(x_1-x_3)^2}\right)
   \\
   && \nonumber \hspace{1.0cm}-8 x_1 x_3 \text{Li}_2\left(\frac{2 x_1}{x_1+x_3-1-\sqrt{1-2
(x_1+x_3)+(x_1-x_3)^2}}\right)
   \\
   && \nonumber \left. \hspace{1.0cm}+8 x_1 x_3 \text{Li}_2\left(\frac{x_1+x_3-1-\sqrt{1-2
(x_1+x_3)+(x_1-x_3)^2}}{2
   x_3}\right)\right]+\mathcal{O}\left(\epsilon\right)\,.
\end{eqnarray}
For $m_1\gg m_3=m_{light}$, expanding the expression (\ref{ff2})
over $m_3$, we get
\begin{eqnarray}
&& PS_4(m_1,0,m_{light},0)=\frac{m^4}{8192 \pi ^5} \left[\frac{1}{3}
(1-x_1) \left(x_1^2+10x_1+1\right)+2 x_1 (x_1+1) \log (x_1)\right.
\\
&& \nonumber \hspace{0.5cm} +x_{\ell} \left(\left(1-x_1^2\right)
(3-4 \log (1-x_1)+2 \log (x_{\ell})) +2 x_1 \log (x_1) (-x_1-4 \log
(1-x_1)+2 \log (x_{\ell})) \right.
\\
&& \left. \nonumber  \hspace{1.3cm}-8 x_1 \text{Li}_2(1-x_1)\right)
\\
&& \nonumber  \hspace{0.5cm} \left.+x_{\ell}^2 (-(1-x_1) (3+4 \log
(1-x_1)-2 \log (x_{\ell}))-2 x_1 \log
(x_1))\right]+\mathcal{O}\left(\epsilon,x_{\ell}^3\right)\,.
\end{eqnarray}
\end{widetext}


\section{Introduction}\label{sec_intro}
Due to the high degree of controllability and isolatedness, trapped
ultra-cold
atoms serve as an ideal system for observing many-body quantum phenomena
\cite{many_body_physics_ultrac_atoms_zwerger_dalibard_bloch} and can even
be employed to simulate quantum systems of quite a broad physical context
\cite{quantum_simulators_with_ultracold_quantum_gases_Bloch2012}.
In particular, there is a growing
interest  in the meanwhile accessible regime where a mean-field description
\cite{Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008,
Pitaevskii_Stringari_Bose-Einstein_Condensation2003} as
given by the Gross-Pitaevskii equation fails. Such states
can be realized in
e.g. optical lattices
\cite{quantum_phase_transition_SF_to_MI_Greiner_Nature2002}. Feshbach
\cite{Feshbach_resonances_ultracold_gase_Grimm_2010} or
confinement induced resonances
\cite{Olshanii_PRL98_quasi_1d_scattering,
observation_of_1d_TG_gas_Kinoshita_Science2004,
confinement_induced_molecules_in_1d_fermigas_Esslinger_PRL2005} can be employed
to tune
the inter-atomic interaction strength. In particular, quasi-one-dimensional
trapping geometries can enhance correlation effects in the strong interaction
regime, leading to fascinating novel phases
\cite{girardeau_bose_fermi_map_1960,
observation_of_1d_TG_gas_Kinoshita_Science2004,
TG-gs_of_atoms_in_opt_lattice_Bloch_Nature2004,
realization_of_sTG_gas_Haller_Science2009} and quantum phase transitions
\cite{pinning_quantum_phase_transisition_Haller2010}. As the mean-field theory
becomes exact for weak interactions and large particle numbers
\cite{proof_BEC_for_dilute_trapped_gases_Lieb_Seiringer,
Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008,
Pitaevskii_Stringari_Bose-Einstein_Condensation2003}, beyond
mean-field physics can also be expected for small ensembles and finite
interaction strengths. The latter regime is experimentally explored in
e.g. arrays of decoupled one-dimensional
tubes typically containing two to 60 atoms
\cite{pinning_quantum_phase_transisition_Haller2010}. Therefore, the
transition from the few- to many-body behaviour is, in particular for the
strongly correlated quantum dynamics, a subject of immediate interest.

Moreover it is meanwhile also experimentally routinely achievable to trap
and manipulate different
components or
species\footnote{In the following, we will use the term {\it species}
irrespectively of whether it refers to different elements, isotopes or
internal states of an isotope.}, which allows for studying distinguishable
subsystems with indistinguishable constituents. Such mixtures can be
realized e.g. by preparing alkali atoms in different hyperfine states
\cite{production_of_2_overlapping_BECs_symopathic_cooling_WiemanPRL1997} or by
trapping different elements
\cite{bec_of_Potassium_by_sympathic_cooling_Inguscio_Science2001}. Due to
the interplay between the intra- and inter-species interaction strengths, these
systems show a number of intriguing features such as phase separation
\cite{dynamics_of_comp_separation_in_binary_BEC_Cornell_PRL1998} including
symbiotic excitations like interlacing vortex lattices with mutually filled
cores \cite{vortex_lattice_dynamics_in_rotating_spinor_BEC_Cornell_PRL2004} and
dark-bright solitons
\cite{oscillations_and_int_of_dark_and_darkbright_solitons_Becker_Nature2008},
spin-charge separation
\cite{spin-charge_sep_in_2comp_Bosegase_Schollwoeck_PRA2008}, various tunneling
effects
\cite{dyn_of_2_spec_BEC_in_double_well_Pindzola_PRA2009,
josephson_osc_in_spinor_BEC_Sanpera_PRA2009,
symmetry_breaking_and_restoring_in_mixture_BEC_double_well_Clark_PRA2009,
quantum_bose_josephson_junct_with_binary_BEC_Citro_JPB2010,
material_barrier_tunneling_1d_fewboson_mixtures_Pflanzer_Zoellner,
interspecies_tunneling_1d_bose_mixtures_Pflanzer_Zoellner,
correl_mixture_ tunneling_double_well_Chatterjee_PRA2012,
impact_of_spatial_correl_on_tunneling_mixture_Capo2012}, collective
excitations
\cite{collect_osc_of_2_colliding_BECs_Inguscio_PRL2000,
dyn_of_2_colliding_BECs_Minardi_PRA2000} and counterflow and paired
superfluidity
\cite{counterflow_SF_of_2spec_commens_lattice_Svistunov_PRL2003,
counterflow_and_paird_S F_in_1d_Bose_mixtures_Clark_PRA2009,
detecting_paired_counterflow_SF_via_dipole_osc_Clark_PRA2011}. The purpose of
the present work is to develop a broadly applicable and efficient {\it ab
initio} method for the quantum dynamics of such mixtures in order to explore
the fundamental dynamical processes in trapped ultra-cold multi-species setups
and to study the few- to many-body transition.

Simulating the quantum dynamics of an interacting many-body system, however,
is a tough task in general due to the exponential\footnote{For distinguishable
particles. Indistinguishable particles lead to a binomial scaling.} scaling of
the state space with the number of particles.
Besides e.g. the time-dependent density matrix renormalization group
approach \cite{real_time_evo_using_DMRG_White_Feiguin_PRL2004,
td_DMRG_using_adaptive_eff_Hilbert_spaces_Vidal2004,
td-DMRG_methods_Schollwoeck2004, DMRG_at_age_of_MPS_Schollwoeck2010}, a
promising concept to soften
this scaling is based on a many-body wave function expansion with respect to a
time-dependent, with the system comoving basis. This idea has been
incorporated in the multi-configuration time-dependent Hartree method (MCTDH)
\cite{mctdh_approach_Meyer_Manthe_Cederbaum,MCTDH_BJMW2000}. Being based
on time-dependent Hartree products as the many-body basis, MCTDH is designed
for distinguishable particles, but has also been applied to bosonic few-body
systems (e.g.
\cite{few_boson_dyn_doubl_wells_Zoellner,
impact_of_spatial_correl_on_tunneling_mixture_Capo2012}). Later the MCTDH
theory has been generalized and extended in several ways: There is the
multi-layer MCTDH (ML-MCTDH) method
\cite{Multilayer_formulation_Wang_Thoss_2003,
Multilayer_multiconfig_Manthe_2008,
 ML_implementation_Vendrell_Meyer_JChemPhys_2011}, which takes correlations
between various subsystems into account and is thus
particularly suitable for system-bath problems with distinguishable degrees
of freedom (e.g. \cite{wang_shao12}). Taking
the fermionic or bosonic
particle exchange \mbox{(anti-)} symmetry in the time-dependent many-body basis
into account, MCTDH has been specialized to
treat larger fermionic (MCTDHF)
\cite{MCTDHF_approach_to_multi_electr_dynamics_in_laser_fields_Scrinzi2003} or
bosonic systems (MCTDHB)
\cite{
role_of_excited_states_in_splitting_of_BEC_by_time_dep_barrier_Steltsov_PRL2007,
MCTDHB_many_body_dynacmics_of_bosonic_systems_phys_rev_a_Alon_Streltsov_Cederbau
m}. Furthermore, a direct
extension of MCTDHB and MCTDHF to treat bose-bose, bose-fermi and fermi-fermi
mixtures has been
developed \cite{MCTDHX_recursively_Cederbaum2012}, including the possibility
of particle-conversions
\cite{MCTDHB_with_particle_conversion_Alon_Streltsov_Cederbaum_PRA2009}. An
alternative approach to systems of indistinguishable particles is the so-called
ML-MCTDH method in second-quantization representation, which employs
the factorization of the many-body Hilbert space into a direct product of
Fock spaces
\cite{MLMCTDH_in_2nd_quantization_representation_Wang_Thoss_Chem_Phys_2009}.

In this work, we derive and apply a novel {\it ab initio} approach to the
non-equilibrium dynamics of ultra-cold correlated bosonic mixtures, which takes
all correlations of the
many-body system into account. We call this method the {\it multi-layer
multi-configuration time-dependent
Hartree method for bosons} (\mbox{ML-MCTDHB}). The
multi-layer
structure of our many-body wave function ansatz allows us to adapt our many-body
basis to system specific inter- and intra-species correlations, which leads to a
beneficial scaling. Moreover, the bosonic exchange symmetry is directly employed
for an efficient treatment of the indistinguishable bosonic subsystems. We apply
ML-MCTDHB to simulate the correlated tunneling dynamics of a
mixture of three bosonic species in a double well trap. It is shown that the
dynamics of the population imbalances of the species significantly differ
for ultra-weak and vanishing inter-species interaction strengths. In
particular, bosons of different kind show a bunching tendency
and the inter-species interaction strengths allow for tuning these correlations
up to a certain degree.

This paper is organized as follows: In section \ref{sec_method}, the
derivation and properties of the \mbox{ML-MCTDHB} method for bosonic
mixtures are presented. \mbox{ML-MCTDHB} is then applied to simulate the
complex tunneling behaviour of a mixture of three bosonic species in
section \ref{sec_appl}. Finally, we summarize our results and embed the
presented \mbox{ML-MCTDHB} theory for mixtures into a more general framework in
section
\ref{sec_concl}.
\section{The \mbox{ML-MCTDHB} method}\label{sec_method}
Let us consider an ensemble of $S$ bosonic species. In the
ultra-cold
regime, the interaction between neutral atoms can be modelled
by a contact interaction \cite{Huang_Yang_paper_PhysRev1957,
Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008,
Pitaevskii_Stringari_Bose-Einstein_Condensation2003}. For simplicity, we
restrict ourselves to
one-dimensional settings, which can be prepared by energetically freezing out
the transversal degrees of freedom
\cite{Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008}. The
Hamiltonian of such a mixture
with $N_\sigma$ bosons of species $\sigma=1,...,S$
reads:
\begin{equation}\label{hamil_tot}
 \hat H =\sum_{\sigma=1}^S \left(\hat \mathcal{H}_\sigma+\hat V_\sigma\right)
	 +\sum_{1\leq \sigma<\sigma'\leq S} \hat W_{\sigma\sigma'}.
\end{equation}
Here, $\hat \mathcal{H}_\sigma$ denotes the one-body Hamiltonian of the species
$\sigma$
containing a in general species dependent trapping potential $U_\sigma$:
\begin{equation}\label{hamil_1b}
 \hat \mathcal{H}_\sigma
\sum_{i=1}^{N_\sigma}
\left(
    \frac{\hat p^\sigma_i\phantom{)\!\!}^2}{2m^\sigma}+
	  U_\sigma(\hat x_i^\sigma)\right),
\end{equation}
and $\hat V_\sigma$, $\hat W_{\sigma\sigma'}$ refer to the intra-species
interaction of species $\sigma$ and to the inter-species interaction between
$\sigma$ and $\sigma'$ bosons, respectively:
\begin{eqnarray}\label{hamil_intra}
 \hat V_\sigma =g_\sigma\sum_{1\leq i<j \leq N_\sigma} \delta(\hat x_i^\sigma-
		\hat x_j^\sigma),\\\label{hamil_inter}
 \hat W_{\sigma\sigma'} = g_{\sigma\sigma'}
		\sum_{i=1}^{N_\sigma}\sum_{j=1}^{N_{\sigma'}}
		\delta(\hat x_i^\sigma-\hat x_j^{\sigma'}).
\end{eqnarray}
Please note that the intra- and inter-species interaction strengths $g_\sigma$,
$g_{\sigma\sigma'}$ have to be properly renormalized with respect to their 3d
values as a consequence of dimensional reduction
\cite{Olshanii_PRL98_quasi_1d_scattering}. We remark that the Hamiltonian may be
explicitly time-dependent for studying driven systems.

\subsection{Wavefunction ansatz}\label{ssec_ansatz}
The \mbox{ML-MCTDHB} method is an {\it ab initio} approach to the time-dependent
Schr\"odinger equation for systems like (\ref{hamil_tot}). To reduce the
number of basis states necessary for a fair representation of the total
wave function $|\Psi(t)\rangle$, we employ a time-dependent, with the system
comoving basis and restrict ourselves to the following class of ansatzes: For
each species $\sigma$, we take $M_\sigma$ time-dependent orthonormal species
states
$|\psi^{(\sigma)}_i(t)\rangle$ ($i=1,...,M_\sigma$), i.e. states of all the
$N_\sigma$ bosons of species $\sigma$, into account. Due to the
distinguishability of bosons of different species, the total wave function is
expanded in terms of Hartree products of these many-body states:
\begin{equation}\label{exp_top}
 |\Psi(t)\rangle = \sum_{i_1=1}^{M_1}...\sum_{i_S=1}^{M_S} A_{i_1,...,i_S}(t)\;
  |\psi^{(1)}_{i_1}(t)\rangle ... |\psi^{(S)}_{i_S}(t)\rangle.
\end{equation}
Each species state $|\psi^{(\sigma)}_i(t)\rangle$ refers to a system of
$N_\sigma$ indistinguishable bosons and should therefore be expanded in terms
of bosonic number states $|\vec n\rangle^\sigma_t$:
\begin{equation}\label{exp_spec}
  |\psi^{(\sigma)}_{i}(t)\rangle = \sum_{\vec n|N_\sigma}
      C^\sigma_{i;\vec n}(t)\;|\vec n\rangle^\sigma_t,
\end{equation}
where we allow each $\sigma$ boson to occupy $m_\sigma$ time-dependent single
particle
functions (SPFs) $|\phi^{(\sigma)}_j(t)\rangle$, indicated by the
time-dependence of bosonic number states $|\vec n\rangle^\sigma_t$. The integer
vector $\vec
n=(n_1,...,n_{m_\sigma})$ contains the occupation number $n_j$ of the $j$-th
SPFs such that all $n_j$'s sum up to $N_\sigma$, indicated by the symbol
``$\vec n|N_\sigma$'' in the summation.

Summarizing, our wave function ansatz consists of three layers: The expansion
coefficients $A_{i_1,...,i_S}(t)$ form the top layer. Then we have the
$C^\sigma_{i;\vec n}(t)$'s on the species layer which allow the species states
to move with the
system and, finally on the particle layer, the SPFs
$|\phi^{(\sigma)}_j(t)\rangle$ allow for rotations of the single particle
basis. It is crucial to notice that, in contrast to the standard method for
solving the time-dependent Schr\"odinger equation by propagating expansion
coefficients while keeping the basis time-independent, ML-MCTDHB is based on an
expansion with respect to a comoving basis with a two-fold time-dependence in
terms of the species states $|\psi^{(\sigma)}_{i}(t)\rangle$ and the SPFs
$|\phi^{(\sigma)}_j(t)\rangle$. This two-fold time-dependence allows for
significantly reducing the number of basis states leading to a very efficient
algorithm. Please also note that our \mbox{ML-MCTDHB} approach to mixtures
conceptually differs from ML-MCTDH in second-quantization representation
\cite{MLMCTDH_in_2nd_quantization_representation_Wang_Thoss_Chem_Phys_2009} by
the facts that we only employ two layers, one for the whole species and one for
the single bosons, but allow for a time-dependent single
particle basis.

Having the number of grid points for representing the SPFs fixed, the
numbers of species states $M_\sigma$ and SPFs $m_\sigma$ serve
as numerical control parameters: Taking $m_\sigma$ to be equal to the number of
grid points and $M_\sigma$ equal to the number of number state configurations,
i.e. $(N_\sigma+m_\sigma-1)!/[N_\sigma!(m_\sigma - 1)!]$, the ansatz
(\ref{exp_top},\ref{exp_spec}) proves to be numerically exact. Opposite to this
full CI limit, the choice $m_\sigma=M_\sigma=1$ leads to the mean-field or
Gross-Pitaevskii approximation
\cite{Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008,
Pitaevskii_Stringari_Bose-Einstein_Condensation2003}. In between
these two limiting case, any choice with $m_\sigma$ equal to or smaller than the
number of
grid points and $M_\sigma\leq(N_\sigma+m_\sigma-1)!/[N_\sigma!(m_\sigma - 1)!]$
is possible, which allows us to adapt our ansatz to system specific intra- and
inter-species correlations. If, for instance, the inter-species interactions are
relatively weak compared to the intra-species interactions, a ``species
mean-field'' ansatz with $M_\sigma=1$ but $m_\sigma>1$ might be sufficient.

\subsection{Equations of motion}\label{ssec_eqm}
Our final task is to find appropriate equations of motion for the ansatz
constituents $A_{i_1,...,i_S}(t)$, $C^\sigma_{i;\vec n}(t)$ and
$|\phi^{(\sigma)}_j(t)\rangle$, whose time dependence we will omit in the
notation from now on. In order
to find the variationally optimal wave function
$|\Psi(t)\rangle$ within  our class of ansatzes for given $M_\sigma$,
$m_\sigma$, we can employ the McLachlan variational principle, which enforces
the minimization of the error of our equations of motion with respect to the
exact Schr\"odinger equation \cite{mc_lachlan_1963}. In practice, however, it
is easier to work with the Dirac-Frenkel variational principle $\langle
\delta\Psi|(i\partial_t-\hat H)|\Psi\rangle=0$ with $|\delta\Psi\rangle$ being
a variation within our ansatz class ($\hbar\equiv1$)
\cite{dirac_variational,frenkel_variational}, which turns out to be equivalent
to McLachlan's variational principle on our
manifold of wave function ansatzes\cite{equivalence_of_td_var_princ}.

The variation of the top layer coefficients $A_{i_1,...,i_S}$ gives us the
usual linear equation of motion known from matrix mechanics:
\begin{equation}\label{eqm_top}
 i\partial_t A_{i_1,...,i_S} = \sum_{j_1=1}^{M_1}...\sum_{j_S=1}^{M_S}
      \langle \psi^{(1)}_{i_1}...\,\psi^{(S)}_{i_S}|\hat H
      |\psi^{(1)}_{j_1}...\,\psi^{(S)}_{j_S}\rangle\;A_{j_1,...,j_S},
\end{equation}
where the Hamiltonian matrix with respect to Hartree products of species
states becomes time-dependent due to the coupling to the $C^\sigma_{i;\vec n}$
coefficients and to the SPFs. Its explicit form is given in
\cite{ml-mctdhb_unpub}.

Varying the species state expansion coefficients $C^\sigma_{i;\vec n}$, we
obtain the following equations of motion on the species layer:
\begin{eqnarray}\label{eqm_spec}
\eqalign\nonumber
 i\partial_t C^\sigma_{i;\vec n} &=&  {^{\sigma\!}\langle}\vec
n|(\mathds{1}-\hat P^{1;\sigma})\;\sum_{\vec m|N_\sigma}
    \Big( \sum_{j,k=1}^{m_\sigma} {[h_\sigma]}_{jk}\,\hat a_{\sigma
j}^\dagger\hat a_{\sigma k} \;|\vec m\rangle^\sigma\;C^\sigma_{i;\vec m} + \\
    &+& \frac{1}{2}\sum_{j,k,q,p=1}^{m_\sigma}
{[v_{\sigma}]}_{jkqp}\,
\hat a_{\sigma j}^\dagger \hat a_{\sigma k}^\dagger \hat a_{\sigma q} \hat
a_{\sigma p} \;|\vec m\rangle^\sigma\;C^\sigma_{i;\vec m}+\\\nonumber
&+& \sum_{\sigma'\neq\sigma}\sum_{s,t=1}^{M_{\sigma}}
\sum_{u,v=1}^{M_{\sigma'}}\sum_{j,k=1}^{m_\sigma}
{[\eta_{1,\sigma}^{-1}]}_{is}
  \;{[ \eta_{ 2, \sigma\sigma'}]}_{sutv}
\;{[w_{\sigma\sigma'}]}^{jk}_{uv}\;
  \,\hat a_{\sigma j}^\dagger\hat
a_{\sigma k} \;|\vec m\rangle^\sigma\;C^\sigma_{t;\vec
m}\Big),
\end{eqnarray}
where $\hat a^{(\dagger)}_{\sigma i}$ denotes the bosonic annihilation
(creation) operator corresponding to the SPF $|\phi^{(\sigma)}_i\rangle$,
obeying the canonical commutation relations $[\hat a_{\sigma i},\hat
a_{\sigma' j}]=0$ and $[\hat a_{\sigma i},\hat
a^\dagger_{\sigma' j}]=\delta_{\sigma\sigma'}\delta_{ij}$. ${[h_\sigma]}_{jk}$
and ${[v_{\sigma}]}_{jkqp}$ represent the matrix elements of the one-body
Hamiltonian and the intra-species interaction potential
with respect to the SPFs, respectively. The inter-species interaction leads to
the
mean-field matrix $[w_{\sigma\sigma'}]$ coupling both SPFs and species states.
The reduced density matrix of the species
$\sigma$ and the reduced density matrix of the subsystem constituted by the
species $\sigma$ and $\sigma'$ ($\sigma\neq\sigma'$) enter (\ref{eqm_spec}) as
$[\eta_{1,\sigma}]$ and
$[ \eta_{ 2, \sigma\sigma'}]$, respectively (cf. (\ref{rho_1spec}),
(\ref{rho_2spec})). The orthonormality of the species
states is ensured by the projector $\hat P^{1;\sigma}=\sum_{s=1}^{M_\sigma}
|\psi^{(\sigma)}_s\rangle\!\langle\psi^{(\sigma)}_s|$. Formulas for the above
ingredients are given in \ref{app_spec_layer} and an efficient scheme for
applying the annihilation and creation operators to the number states can be
found in \cite{ml-mctdhb_unpub} (see also
\cite{
general_mapping_for_bos_and_ferm_ops_in_fock_space_Streltsov_Alon_Cederbaum} in
this context).

Finally, the variation of the SPFs leads to the following non-linear
integro-differential equations:
\begin{eqnarray}\label{eqm_part}
\eqalign\nonumber
 i\partial_t|\phi^{(\sigma)}_i\rangle = (\mathds{1}-\hat P^{2;\sigma})\;
    &\Big(&\hat h_\sigma|\phi^{(\sigma)}_i\rangle + \\
    &+& \sum_{j,k,q,p=1}^{m_\sigma}
  {[\rho_{1,\sigma}^{-1}]}_{ij}
  \;{[ \rho_{ 2, \sigma\sigma}
]}_{jkqp}
  \;{[ \hat v_{\sigma}]}_{kq}\;
  |\phi^{(\sigma)}_p\rangle +\\\nonumber
&+& \sum_{\sigma'\neq\sigma}\sum_{j,q=1}^{m_\sigma}\sum_{k,p=1}^{m_{\sigma'}}
{[\rho_{1,\sigma}^{-1}]}_{ij}
  \;{[ \rho_{ 2, \sigma\sigma'}
]}_{jkqp}
  \;{[\hat
w_{\sigma\sigma'}]}_{kp}\;
  |\phi^{(\sigma)}_q\rangle	\Big).
\end{eqnarray}
Here, $[\rho_{1,\sigma}]$ denotes the
reduced density matrix of a $\sigma$ boson and $[ \rho_{ 2, \sigma\sigma}]$,
$[\rho_{2, \sigma\sigma'}]$ ($\sigma\neq\sigma'$) refer to the reduced two-body
density matrix of
two $\sigma$ bosons, a $\sigma$ and a $\sigma'$ boson, respectively (cf.
(\ref{rho_1part})-(\ref{rho_1A1Bpart})). $\hat h_\sigma$
corresponds to the one-body Hamiltonian $\frac{\hat
p_\sigma^2}{2m^\sigma}+
	  U_\sigma(\hat x_\sigma)$ and the intra- and
inter-species interactions enter these equations of motion in the form of the
mean-field operator matrices $[\hat v_{\sigma}]$ and $[\hat
w_{\sigma\sigma'}]$, respectively. All these ingredients are explicated in
\ref{app_part_layer}.
The projector $\hat P^{2;\sigma}=\sum_{s=1}^{m_\sigma}
|\phi^{(\sigma)}_s\rangle\!\langle\phi^{(\sigma)}_s|$ again ensures the
orthonormality of the SPFs. So we have arrived at a set of highly coupled
evolution equations
(\ref{eqm_top}, \ref{eqm_spec},
\ref{eqm_part}), whose general properties we analyse in the following section.

\subsection{Properties of the \mbox{ML-MCTDHB} theory}\label{ssec_prop}
Derived from the Dirac-Frenkel variational principle, the \mbox{ML-MCTDHB}
evolution
equations preserve both norm and energy \cite{MCTDH_BJMW2000}. Moreover, one
can show that for a Hamiltonian with a (single particle) symmetry,
\mbox{ML-MCTDHB}
respects both the symmetry of the SPFs and the symmetry of the many-body
state, given that initially the SPFs and the many-body
state have a well-defined symmetry \cite{ml-mctdhb_unpub}.

In the full CI limit, i.e. $m_\sigma$ equal to the number of grid points and
$M_\sigma=(N_\sigma+m_\sigma-1)!/[N_\sigma!(m_\sigma - 1)!]$, the projectors in
(\ref{eqm_spec},\ref{eqm_part}) turn into unit operators such that both the
species states and the SPFs become time-independent. In this numerically
exact limit, \mbox{ML-MCTDHB} becomes equivalent to the standard method of
solving the
time-dependent Schr\"odinger equation by propagating only the $A$-coefficients.
The full CI limit, however, is numerically only manageable for extremely
small
particle numbers, whereas ML-MCTDHB being based on a smaller but with the
system optimally comoving basis can treat much larger ensembles.
In the
opposite mean-field limit $m_\sigma=M_\sigma=1$, the time-dependence of the $A$-
and the $C$-coefficients is given by trivial phase factors. With all the various
reduced
density matrices being equal to the c-number one, the equations
(\ref{eqm_part}) just differ from the coupled Gross-Pitaevskii equations of the
mean-field theory for mixtures
\cite{Bose-Einstein_Condensation_Dilute_Gases_Pethick_Smith_2008,
emerg_phenom_BEC_08} by a physically irrelevant phase factor as a consequence of
the projector $\hat P^{2;\sigma}$.

A converged \mbox{ML-MCTDHB} calculation takes all correlations into account.
These
correlations can be studied by means of reduced density matrices of various
subsystems, which the \mbox{ML-MCTDHB} method provides for free.
Single-particle coherence as well as correlations between two bosons of the same
or of different species can be unravelled with the help of
$[\rho_{1,\sigma}]$, $[ \rho_{ 2, \sigma\sigma}]$ and $[\rho_{2,
\sigma\sigma'}]$.
The entropy of a species as well as correlations between two species can be
deduced from $[\eta_{1,\sigma}]$ and $[\eta_{ 2, \sigma\sigma'}]$, for example.
Moreover, an analysis of the natural populations and orbitals of various
subsystems both serves as an internal convergence check (see
below)\cite{MCTDH_BJMW2000} and can give physical insights
\cite{red_density_matrices_and_coherence_of_trapped_bosons_Sakmann_PRA2008}.

In the case of just one species and in the full CI limit on the species
level $M_\sigma=(N_\sigma+m_\sigma-1)!/[N_\sigma!(m_\sigma - 1)!]$, the
\mbox{ML-MCTDHB} theory becomes equivalent to MCTDHB
\cite{
role_of_excited_states_in_splitting_of_BEC_by_time_dep_barrier_Steltsov_PRL2007,
MCTDHB_many_body_dynacmics_of_bosonic_systems_phys_rev_a_Alon_Streltsov_Cederbau
m} and its generalization to mixtures
\cite{MCTDHX_recursively_Cederbaum2012}, respectively. If, however, less
species states are sufficient for a converged
simulation, \mbox{ML-MCTDHB} proves to have a better scaling. With $n$ being the
number of grid points, one has to pay:
\begin{equation}
 \prod_{\sigma=1}^S M_\sigma + \sum_{\sigma=1}^S  \left ( M_\sigma\,
  {N_\sigma+m_\sigma-1 \choose m_\sigma-1} + m_\sigma\,n
\right)
\end{equation}
complex coefficients for storing a \mbox{ML-MCTDHB} wave function, which
should be compared with the costs for a corresponding MCTDHB expansion:
\begin{equation}
 \prod_{\sigma=1}^S {N_\sigma+m_\sigma-1 \choose m_\sigma-1} +
 \sum_{\sigma=1}^S m_\sigma\,n.
\end{equation}
For a detailed scaling comparison of the MCTDH type methods, we refer to
\cite{ml-mctdhb_unpub}.
\section{Application to correlated tunneling dynamics}\label{sec_appl}

Let us now explore the tunneling dynamics of three bosonic
species, refered to as the A, B and C species in the following, in a double
well trap. This setup both unravels interesting correlation
effects and illustrates the beneficial scaling of \mbox{ML-MCTDHB} by
introducing the
extra species layer.

In the following, we assume that the three species are realized as different
hyperfine states of an alkali element resulting in equal masses $m^\sigma$ for
all the bosons. Furthermore, each species shall consist of $N_\sigma=6$ bosons
and shall experience the very same trapping potential made of a harmonic trap
superimposed with a Gaussian at the trap centre, i.e.
$U_{\sigma}(x)\equiv U(x)=x^2/2+
h/\sqrt{2\pi s^2}\;\exp(-x^2/2s^2)$ in harmonic oscillator units
$\hbar=m^\sigma=\omega=1$. We choose $h=3$ and $s=0.2$ for the
height and width of the barrier, respectively, which leads to three
bands below the barrier, each consisting of two single particle eigenstates. The
lowest band is
separated by an energy difference of $1.63$ from the first excited band and its
level spacing amounts to $\Delta E\approx 0.23$ leading to a tunneling period
of $T=2\pi/\Delta E\approx 27$ for non-interacting particles. For the contact
interaction strengths, we take $g_A (N_A-1)=0.2$, $g_B=0.75\,g_A$ and $g_{AB} =
0.05\,\sqrt{g_Ag_B}$. Furthermore, the C bosons are assumed to have no
intra-species interaction, i.e. $g_C=0$, but an attractive, vanishing or
repulsive coupling to the bosons of species A and B: $g_{XC}\equiv
g_{AC}=g_{BC}\in\{-0.5\,g_{AB},0.0,0.5\,g_{AB}\}$. Anticipating the results, we
will show that
this very weak interaction of strength $g_{XC}$ has
a significant impact on the correlation between the A and B bosons.

As the particle numbers are the same
for all species and because of the not too
different interaction strengths, we provide for each species the
same number of species states, $M_\sigma\equiv M$, and SPFs\footnote{The SPFs
are represented by
means of a harmonic discrete variable representation (DVR)
\cite{MCTDH_BJMW2000}.},
$m_\sigma\equiv m$. For preparing the initial state of the mixture, we
block the right well by means of a high step function potential. All
bosons are then put into
the ground state of the resulting single particle Hamiltonian and, afterwards,
we let the interacting many-body system relax to its ground state by propagating
the
\mbox{ML-MCTDHB} equations of motion in imaginary time. Ramping down the
step function potential instantaneously, the resulting many-body
state
is finally propagated in real time in the original double well trap.
Afterwards, we infer the probability of a particle to be in some well and the
probability of finding two particles of the same or different species in the
same well from the corresponding reduced one-body and two-body
density matrices.

Here we would like to point out that we do not aim at an
exhaustive study of this setup. Rather than showing a systematic parameter
scan, we would like to present one striking
example of multi-species non-equilibrium dynamics hardly being accessible in
this precision by other methods but \mbox{ML-MCTDHB}, thereby illustrating
the beneficial scaling and efficiency of the method. As we shall see, this setup
shows very interesting correlation effects.

\subsection{Short time tunneling dynamics}
Let us firstly focus on the tunneling dynamics for an attractive coupling of
the C bosons to the bosons of the A and B species, i.e. $g_{XC}<0$, up to
time $t=100$. From
figure
\ref{fig_6x6x6_tunnel_prob}, we
see that the A, B and C bosons exhibit Rabi tunneling with respect to the
tunneling period on this time interval. The amplitude of the probability
oscillations, however, decreases in the course
of time for the A and the
B bosons. This decrease can be interpreted as a temporal equilibration of the
occupation probability of the left well as one can infer from the inset of
figure \ref{fig_6x6x6_tunnel_prob} showing a somewhat lower accuracy long time
propagation (see below). We also clearly see that
the decrease of the probability amplitude
is a genuine many-body property, not present in the mean-field description via
coupled Gross-Pitaevskii equations. In
contrast to this, the tunneling
amplitude of the C bosons is not damped and its dynamics in the many-body
calculation coincides with the mean-field description. This is a consequence
of the vanishing intra-species and the very weak inter-species interaction
strength.

A further phenomenon not
captured in the mean-field picture is unravelled in figure
\ref{fig_6x6x6_tunnel_jointprob}: The probabilities for finding two bosons
of the same species in the same well oscillate between 0.5 and
1.0 with the frequency $2/T$ in the mean-field calculation. In the many-body
calculation, however, the probability for finding two A or B bosons in
the same well features damped oscillations leading to a
saturation of 0.73, which indicates a bunching tendency, while the probability
of finding two C bosons stays oscillating between 0.5 and 1.0.

\begin{figure}
 \centering
 \includegraphics[width=0.8\textwidth]
{1b_prob_A_B_C_fig1.eps}
 \caption{Tunneling dynamics of 6 A, 6 B and 6 C bosons, initially loaded into
the left well, for the interaction strengths
$g_A (N_A-1)=0.2$, $g_B=0.75\,g_A$, $g_C=0.0$, $g_{AB} =
0.05\,\sqrt{g_Ag_B}$ and $g_{XC}\equiv g_{AC}=g_{BC}=-0.5\,g_{AB}$: The
probabilities for
finding
an A boson (blue line), a B boson (red) and a C boson (green) in the left well.
Solid lines: \mbox{ML-MCTDHB} data for $m=M=3$. Crosses: $m=M=4$. Dotted
lines: mean-field results ($m=M=1$). The first three
Rabi-tunneling periods are represented by
the dashed vertical lines. Inset: Longtime propagation of the probability to
find an A boson left ($m=3$, $M=5$) with $g_A$, $g_B$ as above. Four cases: (i)
$g_{AB}=g_{XC}=0$ (black line), $g_{AB}=0.05\,\sqrt{g_Ag_B}$ with (ii)
$g_{XC}=0$ (green), (iii) $g_{XC}=0.5\,g_{AB}$ (red), (iv)
$g_{XC}=-0.5\,g_{AB}$ (blue).}
\label{fig_6x6x6_tunnel_prob}
\end{figure}

\begin{figure}
 \centering
 \includegraphics[width=0.8\textwidth]
{2b_prob_AA_BB_CC_fig2.eps}
 \caption{Time
evolution of the probabilities for finding two A bosons (blue line), two B
bosons (red), two
C bosons (green) in the same well for the setting of figure
\ref{fig_6x6x6_tunnel_prob}. Solid lines: \mbox{ML-MCTDHB} data with $m=M=3$.
Crosses: $m=M=4$. Dotted lines: mean-field results.}
\label{fig_6x6x6_tunnel_jointprob}
\end{figure}

For discussing the convergence of the simulation, the \mbox{ML-MCTDHB}
calculation for $m=M=3$ is compared with the results for $m=M=4$ in
figures \ref{fig_6x6x6_tunnel_prob} and \ref{fig_6x6x6_tunnel_jointprob}.
The single particle probabilities show an excellent
agreement. Only for the joint probability of finding two particles of the same
or different (not shown)
kind in the same well, there are marginal deviations.
Hence, we can definitely regard the simulation as being converged. This
judgement is also
supported by the time evolution of the natural populations: From figure
\ref{sfig_6x6x6_natpop_spec}, we infer that most of the time only two natural
orbitals
significantly contribute to the reduced density matrix of the whole species A
and, hence, to the total wave function. For times larger
than $\sim70$, a third species state gains a weight more than 1\%. Thus,
much less species states than $M=(N+m-1)!/[N!(m - 1)!]$, i.e. the full CI
limit on the species layer, are enough for a fair representation of the total
wave function.
Figure
\ref{sfig_6x6x6_natpop_partA}
shows that the initially fully condensed state of the A bosons evolves into
a two-fold fragmented state. Increasing the
number of particle
SPFs from $m=3$ to $4$ just leads to a reshuffling of the third-highest natural
population without affecting the results. The natural populations corresponding
to a B boson show a similar behaviour due to the similar intra-species
interaction strengths (not shown). In contrast to this, the C bosons stay in a
condensed state and become depleted only by $1.2\%$ in the long time
propagation up to $t=300$ (not shown). Please note that the extra species layer
is crucial for this
convergence check: Our $m=M=4$ simulation lasted roughly a week\footnote{For
$n=250$ gridpoints on an Intel\textregistered$\,$ Xeon\textregistered$\,$  CPU
E5530  with 2.40GHz.}, while a
corresponding MCTDHB calculation would require to propagate 146 times
more coefficients.
\begin{figure}[h]
\centering
\subfigure[]{
\centering
\includegraphics[width=0.5\textwidth]
{natpop_spec_fig3a.eps}
 \label{sfig_6x6x6_natpop_spec}
}
\subfigure[]{
 \centering
 \includegraphics[width=0.5
\textwidth]{natpop_partA_fig3b.eps}
 \label{sfig_6x6x6_natpop_partA}
}
 \caption{Time evolution of the eigenvalues of the reduced
density matrices of (a)
the species A, $[\eta_{1,A}]$, and (b) a single A boson, $[\rho_{1,A}]$,
(cf. (\ref{rho_1spec}), (\ref{rho_1part})). Solid lines: $m=M=3$
\mbox{ML-MCTDHB}
results. Crosses: $m=M=4$. The horizontal dashed line in (b)
indicates an ordinate value of 0.5, i.e. the perfect twofold fragmented state.
Parameters as in figure
\ref{fig_6x6x6_tunnel_prob}.}
\label{fig_6x6_natpop_part}
\end{figure}

\subsection{Long time propagation and build up of correlations}
Now let us explore the tunneling on longer time scales with a somewhat
lower
accuracy calculation choosing $m=3$ and $M=5$. A comparison with a $m=3$, $M=4$
simulation shows only very small, quantitative deviations in the observables
under consideration (plots not shown). In the inset of figure
\ref{fig_6x6x6_tunnel_prob}, the time-evolution
of the probability to find an A boson in the left well is shown comparing four
different situations, namely $g_{AB}> 0$ with $g_{XC}=0$,
$g_{XC}>0$ or $g_{XC}<0$ and
$g_{AB}=g_{XC}=0$. Although all the inter-species interaction strengths
are much smaller than $g_A$ and $g_B$, their concrete values have a
strong influence on the tunneling dynamics: For no inter-species interactions,
there happens to be a partial revival of the tunneling oscillation after a
temporal equilibration. In the case of $g_{AB}>0$, only for a
vanishing or attractive coupling of the C species to the other species one can
observe such a temporal equilibration state with subsequent partial
tunneling revival of, however, smaller amplitude in comparison to the
former case. A repulsive coupling between the C and the other species does
not
result in a complete temporal equilibration to a probability of $0.5$ but rather
leads to a reduction of the amplitude of the oscillations around the
equilibration value to $0.12$. While the B bosons show a
dynamics similar to the A bosons, the C bosons tunnel almost unaffected by the
inter-species interactions and exhibit mean-field tunneling
oscillations (plots not shown).

In order to measure the correlations between different
species, we compare the conditional probability of finding a $\sigma$ boson in
e.g. the left well given that a $\sigma'$ boson has already been found there
with the marginal probability of finding a $\sigma$ boson in the left well:
Let $P(\sigma, L;\sigma', L)$ denote the probability for finding a $\sigma$ and
a $\sigma'$ boson in
the left well and let $P(\sigma, L)$ ($P(\sigma', L)$) be the probability for
finding a $\sigma$ ($\sigma'$) boson in the left well. Then the above mentioned
correlation measure reads $f_{LL}(\sigma,\sigma'):=P(\sigma, L;\sigma', L)
/[P(\sigma, L)P(\sigma', L)]$ and we similarly define $f_{RR}(\sigma,\sigma')$
for the right well. Please note that $(f_{LL},f_{RR})$ is a
straightforward extension of the
diagonal elements of the $g_2$ coherence / correlation measure
\cite{quantum_theory_of_opt_coherence_Glauber1963} to spatially discrete systems
with
distinguishable components. The dynamics of the centre of mass positions of the
$\sigma$ and the $\sigma'$
species has an impact on $f_{LL}$ and $f_{RR}$, of course.
In order to diminish this impact, we
finally construct our correlation measure for finding a $\sigma$ and a $\sigma'$
boson in the same well as
$f(\sigma,\sigma'):=([f_{LL}(\sigma,\sigma')^2+f_{RR}(\sigma,\sigma')^2]
/2)^{1/2}$. If the $\sigma$ and the $\sigma'$ bosons tunnel independently,
$f(\sigma,\sigma')$ will be unity. A value of $f(\sigma,\sigma')$ greater
(smaller) than one indicates an overall bunching (anti-bunching) tendency.

In figure \ref{fig_corr}, we find that a bunching tendency between an A and a B
boson clearly builds up with a maximal correlation measure $f(A,B)$ of up to
$40\%$ above unity. For  $75\lesssim t \lesssim
225$, this bunching tendency turns out to be most intense for
the repulsive coupling of the C bosons to the other species and becomes least
intense for an attractive coupling. In the absence of inter-species interactions
with the C species, i.e. if $g_{AB}>0$ is the only non-vanishing inter-species
interaction strength, the correlation measure $f(A,B)$ lies mostly inbetween
these two curves. For a large
fraction of the propagation time, the coupling of the C bosons to the other
two species can thus control the inter-species correlations between species A
and B up to a certain degree. Due to the fact that the C bosons
approximately
perform Rabi-oscillations with respect to the occupation probability of the
left well, one might come to the conclusion that the C bosons provide a
time-dependent potential for the other two species hardly experiencing
a backaction on the considered time-scale. That this descriptive
picture can only be approximately valid up to a certain time, can be inferred
from the second largest natural population of $[\eta_{1,C}]$, which monotonously
increases up to $6.2\%$ ($g_{XC}<0$) and $14\%$ ($g_{XC}>0$), respectively.
\begin{figure}
 \centering
 \includegraphics[width=0.8\textwidth]
{interspec_corr_fig4.eps}
 \caption{Time evolution of the inter-species correlation measure $f(A,B)$
for $g_{AB}=0.05\,\sqrt{g_Ag_B}$ and  $g_{XC}=0$ (black line),
$g_{XC}=0.5\,g_{AB}$ (red) and $g_{XC}=-0.5\,g_{AB}$ (blue). $m=3$, $M=5$.
Corresponding to
absent inter-species interactions, the dash-dotted
serves as a reference. Other parameters as in figure
\ref{fig_6x6x6_tunnel_prob}.}
\label{fig_corr}
\end{figure}
\section{Conclusion and outlook}\label{sec_concl}
We have presented a novel {\it ab initio} method for
simulating the non-equilibrium dynamics of mixtures of ultra-cold bosons. In
particular, ML-MCTDHB is suitable for dealing with explicitly time-dependent
systems, which will be explored in future works. Being based on an expansion in
terms of permanents and on a multi-layer
ansatz, our \mbox{ML-MCTDHB} method takes optimally and efficiently the bosonic
exchange symmetry within each species into account and allows for adapting the
ansatz to
system specific intra- and inter-species correlations. Hereby, the numbers of
provided single particle functions and species states serve as control
parameters for ensuring convergence. For any choice of these numbers of basis
functions, \mbox{ML-MCTDHB} rotates the species states and single particle
functions such that one obtains a variationally optimal representation of the
many-body wave function at any instant in time. This allows to achieve
convergence with a much smaller basis than methods being based on a
time-independent basis. Moreover, if the inter-species interactions are not too
strong, i.e. do not require to consider as many species states as there are
number state configurations for a given number of single particle functions,
ML-MCTDHB proves to have a much better scaling than the best state-of-the-art
method MCTDHB
\cite{
role_of_excited_states_in_splitting_of_BEC_by_time_dep_barrier_Steltsov_PRL2007,
MCTDHB_many_body_dynacmics_of_bosonic_systems_phys_rev_a_Alon_Streltsov_Cederbau
m, MCTDHX_recursively_Cederbaum2012}.
In the case of only a single species state and one
single particle function, \mbox{ML-MCTDHB}
reduces to coupled Gross-Pitaevskii equations.

Employing ML-MCTDHB for a tunneling scenario of three species, we have entered
a parameter regime which is hardly accessible by other methods in such a
controlled precision. Our simulations show that the imbalances of the
populations can feature a temporal equilibration with subsequent revival of the
population oscillations, where the duration of and the fluctuations around the
equilibration state as well as the degree of completeness of the revival
crucially depend on the inter-species interaction strengths. In our setup, we
have furthermore found two-body bunching correlations between
the first two species. The strength of this correlation can be tuned by a weak
attractive or repulsive coupling of the third species - with no intra-species
interaction - to the first two species without significantly altering the
tunneling dynamics of that third species.

In this paper, \mbox{ML-MCTDHB} has been formulated for systems confined by
quasi-one-dimensional traps and interacting via contact interactions. A
direct generalization to arbitrary dimensions and interaction potentials is
possible, of course. Moreover, it is also feasible to generalize
\mbox{ML-MCTDHB} further by applying the multi-layering concept on the level of
the single particle functions, which allows for optimally describing bosons in
quasi-one- or -two-dimensional traps embedded in three-dimensional space with or
without internal degrees of freedom \cite{ml-mctdhb_unpub}.
Incorporating internal degrees of freedom on the level of SPFs then allows
for taking particle converting interactions into account.
\ack
The authors would like to thank Hans-Dieter Meyer and Jan Stockhofe for
fruitful discussions on MCTDH methods and symmetry conservation.
Particularly, the authors would like to thank Jan Stockhofe for the DVR
implementation of the \mbox{ML-MCTDHB}
code. S.K. gratefully acknowledges financial support by the Studienstiftung des
deutschen Volkes. L.C. and P.S. gratefully acknowledge funding by the
Deutsche Forschungsgemeinschaft in the framework of the SFB 925 ``Light induced
dynamics and control of correlated quantum systems''.