\section{Introduction}
Spin-orbit coupling (SOC) is an essential ingredient in the study of exotic
phases in correlated electron systems,\cite{WitczakKrempa:2014hz}
 such as unconventional
superconductivity in $4d$ transition-metal oxides,\cite{Mackenzie:2003er,Wang:2011cm,Watanabe:2013ge,Kim:2014hn,Yang:2014bc,Meng:2014jo,Kim:2015kf,Borisenko:2015ex}
topological phases of matter in quantum spin-Hall insulators,\cite{Kane:2005hl,Hasan:2010ku,Qi:2011hb}
excitonic
insulators,~\cite{Khaliullin:2013du,Kunes:2014ea,Sato:2015hq,Kim:2017kt,Sato:2019by} and
 Kitaev-model-based insulators,\cite{Jackeli:2009hz,rau2016spin,savary2016quantum,winter2017models,hermanns2018physics,takagi2019concept}
to mention a few.
A prototypical minimal
 model that includes the interplay between spin-orbit coupling
and correlations is the relativistic  multiorbital Hubbard model.
Its nonrelativistic counterpart has been intensively investigated in the past and shows a rich phase
diagram.\cite{Werner2008,Medici2011,Kunes:2014ea,Hoshino2016}
The relativistic multiorbital Hubbard model,\cite{Watanabe:2013ge,Meng:2014jo,Sato:2015hq,Shinaoka:2015cg,Kim:2017kt,Sato:2019by} however, is much less understood since
the choice of algorithms is strongly limited
due to the extra computational complications associated with (multiple-)spin-orbit-coupled
 degrees of freedom.

One of the promising formalisms to investigate the
Hubbard model and its generalizations is the dynamical
mean-field theory (DMFT)~\cite{Georges:1996zz,Kotliar2004} that has
provided important insights into multiorbital physics
also in combination with
{\it ab initio} calculations for real materials.\cite{Kotliar:2006fl,yin2011kinetic,lechermann2016electron}
The continuous-time quantum Monte Carlo method (CTQMC),\cite{Gull:2011jd}
particularly the hybridization expansion algorithm (CTHYB),\cite{Werner:2006ko,Werner:2006iz} is the most widely used impurity solver in multiorbital DMFT calculations.
However, CTHYB suffers
 from the notorious sign problem when the SOC is included in the calculations.
The sign problem
grows exponentially with inverse
temperature~\cite{Loh:1990hj,Troyer:2005hv}
and
typically prevents the study of low-temperature
symmetry-broken phases.
Alleviating the sign problem in CTHYB would
help improve our
understanding of phenomena determined
by the interplay of spin-orbit coupling and correlations.\cite{Mackenzie:2003er,Kim:2015kf,Borisenko:2015ex,2015PhRvL.115x7001H,Kim:2018hj}

For quantum Monte Carlo (QMC) algorithms based on auxiliary fields,
there have been various successful advances which
unveiled the origin of the sign problem and suggested a solution in specific
cases,\cite{Chandrasekharan:ul,Wang:2015hm,Iglovikov:2015bj,Kung:2016he,Hann:2017ca,Yoo:2005ib,Alet:2016dd,Huffman:2016ih,Honecker:2016bs}
including the recently developed idea of Majorana
symmetry.\cite{Li:2015jf,Li:2016fi,Wei:2016cb}  The rotationally invariant
Hund coupling in the SO-coupled multiorbital Hubbard system,
 however, generates rather complex interaction terms.
Furthermore, the nonlocal-in-time expansion scheme of CTHYB makes it difficult to track the origin of the fermionic sign on the world line configuration.

In this paper, we systematically study the nature of the sign problem of the
CTHYB for the SO-coupled three-orbital Hubbard model and propose
a strategy to alleviate it.
We employ a numerical sign-optimization scheme, called spontaneous
perturbation stochastic approximation (SPSA),\cite{Spall:1992fh} to determine the optimal basis
in terms of average sign. Remarkably, this optimal basis can be well approximated by a simple one,
 which we denote as $V_{j3/2\mathrm{BA}}$.
The $V_{j3/2\mathrm{BA}}$ basis is obtained from the relativistic basis by a bonding-antibonding transformation.

\section{Model and Method}
Our model Hamiltonian is composed of the three terms $\mathcal{H} = \mathcal{H}_{\text{t}} + \mathcal{H}_{\text{soc}} + \mathcal{H}_{\text{int}}$, where $\mathcal{H}_{\text{t}}$, $\mathcal{H}_{\text{soc}}$, and $\mathcal{H}_{\text{int}}$ are, respectively,
the electron hopping, SOC, and Coulomb interaction terms.
For the noninteracting electron hopping part we assume a degenerate semicircular density of states $\rho(\omega) = (1/\pi D)\sqrt{1-(\omega/D)^2}$ and set the half-bandwidth $D$ as the unit of energy.
In the orbital-spin basis $V_{\mathrm{os}}$, $\mathcal{H}_{\mathrm{soc}}$ has the form
\begin{equation}
	\mathcal{H}_{\text{soc}} = -\lambda\sum^{}_{i}\sum^{}_{\substack{\alpha\alpha'\\ \sigma\sigma'}}c^{\dagger}_{i\alpha\sigma}\langle \alpha\sigma|\mathbf{L}_{\text{eff}}\cdot \mathbf{S}|\alpha'\sigma'\rangle c^{}_{i\alpha'\sigma'}~,
	\label{}
\end{equation}
where $\mathbf{L}_{\text{eff}}$ is the $l=1$ orbital angular momentum operator and $\mathbf{S}$ is the spin operator.
$c^{}_{i\alpha\sigma}$ ($c^{\dagger}_{i\alpha\sigma}$) is the electron annihilation (creation) operator of orbital $\alpha$ ($yz,zx,xy$) and spin $\sigma$ ($\uparrow,\downarrow$) at lattice site $i$.
We introduce the Slater-Kanamori form of the Coulomb interaction including the spin-flip and pair-hopping terms:
\begin{eqnarray}
	\mathcal{H}_{\text{int}} &=& U\sum^{}_{i,\alpha}n_{i\alpha\uparrow}n_{i\alpha\downarrow} + \sum^{}_{\substack{i,\alpha<\alpha' \\ \sigma\sigma'}}(U'-J_{\rm H}\delta_{\sigma\sigma'})n_{i\alpha\sigma}n_{i\alpha'\sigma'} \nonumber\\
	&& -J_{\rm H}\sum^{}_{i,\alpha<\alpha'}(c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\alpha'\downarrow}c^{}_{i\alpha'\uparrow}c^{}_{i\alpha\downarrow} + \text{H.c.})\nonumber\\
	&& +J_{\rm H}\sum^{}_{i,\alpha<\alpha'}(c^{\dagger}_{i\alpha\uparrow}c^{\dagger}_{i\alpha\downarrow}c^{}_{i\alpha'\downarrow}c^{}_{i\alpha'\uparrow} + \text{H.c.})~.
  \label{eqn:kanamori}
\end{eqnarray}
Here, $U$ ($U'$) is the
intraorbital (interorbital) Coulomb repulsion, and $J_{\text{H}}$ is the strength of the Hund coupling.
To make $\mathcal{H}_{\mathrm{int}}$ rotationally invariant, we choose $U'=U-2J_{\mathrm{H}}$~.
This model Hamiltonian covers a wide range of $4d$ and $5d$ materials listed in Appendix~\ref{app:material} by tuning the electron density, $\lambda$, $U$, and $J_{\mathrm{H}}$~.

We solve the model Hamiltonian within the framework of DMFT.\cite{Georges1996}
The lattice model is mapped onto a quantum impurity model
by means of a self-consistency relation.
The effective action of the quantum impurity model is written in terms of the local action $\mathcal{S}_{\mathrm{loc}}$ and the hybridization function $\Delta$,
\begin{eqnarray}
	\mathcal{S}_{\text{imp}} &=& \mathcal{S}_{\text{loc}} \nonumber\\
	&+& \sum^{}_{\substack{\alpha\alpha'\\\sigma\sigma'}}\int_{0}^{\beta}d\tau d\tau'~c^{\dagger}_{\alpha\sigma}(\tau)\Delta_{\alpha\sigma,\alpha'\sigma'}(\tau-\tau')c^{}_{\alpha'\sigma'}(\tau')~.\nonumber\\
	\label{}
\end{eqnarray}
Here, $c^{}_{\alpha\sigma}$ and $c^{\dagger}_{\alpha\sigma}$ are the fermionic
annihilation and creation operators
for orbital $\alpha$ and spin $\sigma$ for the impurity site.
CTHYB solves this impurity problem by performing an expansion of the partition function in powers of the hybridization function,
\begin{eqnarray}
	\mathcal{Z}_{\text{imp}} &=& \mathcal{Z}_{\mathrm{bath}}\sum^{\infty}_{k=0}\frac{1}{(k!)^2}\sum^{}_{\substack{\alpha_1\alpha'_1\\\sigma_1\sigma'_1}}\cdots\sum^{}_{\substack{\alpha_k\alpha'_k\\\sigma_k\sigma'_k}} \nonumber\\
	&\times&\int_{0}^{\beta}d\tau_1\int_{0}^{\beta}d\tau'_1\cdots \int_{0}^{\beta}d\tau_k\int_{0}^{\beta}d\tau'_k~\det\left[ \underline{\Delta}^{(k)} \right]\nonumber\\
	&\times& \left\langle\mathcal{T}_\tau c^{}_{\alpha'_1\sigma'_1}(\tau'_1)c^{\dagger}_{\alpha_1\sigma_1}(\tau_1) \cdots c^{}_{\alpha'_k\sigma'_k}(\tau'_k)c^{\dagger}_{\alpha_k\sigma_k}(\tau_k)\right\rangle_{\mathcal{S}_{\mathrm{loc}}}.\nonumber\\
	\label{eqn:Zimp}
\end{eqnarray}
Here, $\underline{\Delta}^{(k)}$ represents the $k\times k$ matrix whose elements are
\begin{equation}
	\underline{\Delta}^{(k)}_{ij} = \Delta_{\alpha_i\sigma_i,\alpha'_j\sigma'_j}(\tau_i-\tau'_j)~.
	\label{}
\end{equation}
The hybridization function contains the relevant information about the continuous bath degrees of freedom. Since $\Delta$ is directly updated in the DMFT self-consistency loop, CTHYB does not rely on a discretization or truncation of the bath, in contrast to other impurity solvers like the exact diagonalization or the numerical renormalization group.\cite{Georges:1996zz,Gull:2011jd}
In that sense, CTHYB is numerically exact in treating the bath degrees of freedom.

The resulting partition function, Eq.~(\ref{eqn:Zimp}), can be expressed as a sum of configuration weights,
\begin{equation}
	\mathcal{Z}_{\text{imp}} = \sum^{}_{\mathbf{x}} \omega(\mathbf{x})
	\label{}
\end{equation}
where
\begin{widetext}
\begin{equation}
	\omega(\mathbf{x}) = \mathcal{Z}_{\text{bath}}\frac{(d\tau)^{2k}}{(k!)^2}\left\langle\mathcal{T}_\tau c^{}_{\alpha'_1\sigma'_1}(\tau'_1)c^{\dagger}_{\alpha_1\sigma_1}(\tau_1) \cdots c^{}_{\alpha'_k\sigma'_k}(\tau'_k)c^{\dagger}_{\alpha_k\sigma_k}(\tau_k)\right\rangle_{\mathcal{S}_{\mathrm{loc}}}\det\left[ \underline{\Delta}^{(k)} \right]~.
	\label{}
\end{equation}
\end{widetext}
Here, $\mathbf{x}$ represents the configuration composed of $\left\{k; (\alpha_1,\sigma_1,\tau_1),(\alpha'_1,\sigma'_1,\tau'_1),\dots (\alpha_k,\sigma_k,\tau_k),(\alpha'_k,\sigma'_k,\tau'_k)\right\}$.
One can measure a general time-local observable $\mathcal{O}$ using MC sampling as
\begin{equation}
	\langle \mathcal{O}\rangle_{\text{MC}} = \frac{\sum^{}_{\mathbf{x}}\mathcal{O}(\mathbf{x})\omega(\mathbf{x})}{\sum^{}_{\mathbf{x}}\omega(\mathbf{x})}~.
	\label{}
\end{equation}
A potential sign problem appears when the weight $\omega(\mathbf{x})$ becomes negative for certain configurations.
In this case the observable can be sampled through the modified expression
\begin{equation}
	\langle \mathcal{O}\rangle_{\text{MC}} = \frac{\sum^{}_{\mathbf{x}}\text{sign}[\omega(\mathbf{x})]~\mathcal{O}(\mathbf{x})~|\omega(\mathbf{x})|}{\sum^{}_{\mathbf{x}}\text{sign}[\omega(\mathbf{x})]~|\omega(\mathbf{x})|}~.
	\label{}
\end{equation}
A severe sign problem occurs if the small average sign
\begin{equation}
	\langle \mathrm{sign}\rangle=\frac{\sum^{}_{\mathbf{x}}\mathrm{sign}[\omega(\mathbf{x})]|\omega(\mathbf{x})|}{\sum^{}_{\mathbf{x}}|\omega(\mathbf{x})|}
	\label{}
\end{equation}
amplifies the error propagation in the estimation of the observables.
Especially, for low temperatures, the relative error of the average sign becomes exponentially large as a function of the inverse temperature $\beta$.
Thus, reliable estimates of observables cannot be obtained in polynomial time.
An analogous phase problem appears if the configuration weights are complex.
However, since the average phase is always real, we will for simplicity use the terms sign problem and average sign in the following.

An important aspect to keep in mind is that the sign problem depends on the single-particle basis.
For example, it has been reported that if the electron operators $\tilde{c}^{}_{j,m}$ are expressed in the relativistic $j_{\text{eff}}$ basis $V_{j\mathrm{eff}}$,
\begin{eqnarray}
  \left(
  \begin{array}{c}
    \tilde{c}^{}_{\frac{1}{2},+\frac{1}{2}}\\
    \tilde{c}^{}_{\frac{3}{2},+\frac{1}{2}}\\
    \tilde{c}^{}_{\frac{3}{2},-\frac{3}{2}}
  \end{array}
  \right)
  &=& \frac{1}{\sqrt{6}}
  \left(
  \begin{array}{ccc}
	  -\sqrt{2} & +i\sqrt{2} & -\sqrt{2}\\
	  -1 & i & 2 \\
	  -\sqrt{3} & i\sqrt{3} & 0 \\
  \end{array}
  \right)
  \left(
  \begin{array}{c}
	  c^{}_{yz,\downarrow}\\
	  c^{}_{zx,\downarrow}\\
	  c^{}_{xy,\uparrow}\\
  \end{array}
  \right)~,\nonumber\\
  \left(
  \begin{array}{c}
    \tilde{c}^{}_{\frac{1}{2},-\frac{1}{2}}\\
    \tilde{c}^{}_{\frac{3}{2},-\frac{1}{2}}\\
    \tilde{c}^{}_{\frac{3}{2},+\frac{3}{2}}
  \end{array}
  \right)
  &=& \frac{1}{\sqrt{6}}
  \left(
  \begin{array}{ccc}
	  -\sqrt{2} & -i\sqrt{2} & +\sqrt{2}\\
	  +1 & i & 2 \\
	  +\sqrt{3} & i\sqrt{3} & 0 \\
  \end{array}
  \right)
  \left(
  \begin{array}{c}
	  c^{}_{yz,\uparrow}\\
	  c^{}_{zx,\uparrow}\\
	  c^{}_{xy,\downarrow}\\
  \end{array}
  \right),\hspace{8mm}
  \label{eqn:Ujeff}
\end{eqnarray}
the average sign is improved~\cite{Sato:2015hq} due to the diagonalized hybridization function matrix.\cite{Meng:2014jo}  (An average sign of 1 means no sign problem, while an average sign approaching 0 means a severe sign problem.)

To further optimize the basis beyond $V_{j\mathrm{eff}}$, we introduce the
numerical optimization scheme SPSA. In this scheme the average sign becomes the
objective function on the parameter space where each point corresponds to a
single-particle basis used in CTHYB.
The SPSA approximates the gradient for a given parameter point.
At each iteration, the objective function is measured at both the positively and negatively perturbed parameter points along a stochastically chosen direction.
The gradient is evaluated from this two-point measurement, and the parameter point is updated.
Compared to the finite-difference stochastic approximation which involves a number of measurements proportional to the dimension of the parameter space, the SPSA potentially converges faster in time if the parameter space is high-dimensional and the measurement of the objective function is computationally expensive.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.45\textwidth]{nldep.eps}
	\caption{
		Average sign for four different bases as a function of (a) the electron density $n$ and (b) and (c) the strength of SOC $\lambda$ for $T=0.03$, $U=2.0$, and $J_{\mathrm{H}}/U=0.15$~.
		$V_{\mathrm{opt}}$, $V_{j3/2\mathrm{BA}}$, $V_{j \mathrm{eff}}$, and $V_{\mathrm{os}}$ represent, respectively, the optimal basis found by the SPSA, the $j=3/2$ bonding-antibonding basis, the $j_{\mathrm{eff}}$ basis, and the orbital-spin basis.
		The average sign is calculated using the self-consistent paramagnetic solution for a given parameter set.
		The open squares and circles at $n=3$ correspond to the average signs of the Mott insulator solution.
}
	\label{fig:nldep}
\end{figure}
\section{Results}
In Fig.~\ref{fig:nldep} we plot the average sign for the four different bases considered in this work as a function of electron density
[Fig.~\ref{fig:nldep}(a)]  and spin-orbit coupling strength [Fig.~\ref{fig:nldep}(b) and (c)].
The orbital-spin basis $V_{\mathrm{os}}$ has a severe sign problem
since a drastic drop in the average sign of $V_{\mathrm{os}}$ appears as a function of the SOC strength both below and above half-filling.
Figure~\ref{fig:Udep} shows the evolution of the average sign as a function of $U$ ($J_{\mathrm{H}}$) at various electron fillings.
The suppressed average sign of $V_{\mathrm{os}}$ in the noninteracting limit
  implies that the source of the sign problem in this basis is the off-diagonal hybridization function generated by the SOC.
The basis $V_{j\mathrm{eff}}$, defined in Eq.~(\ref{eqn:Ujeff}), diagonalizes
the hybridization-function matrix of the Hamiltonian and thus recovers the
sign-free behavior in the noninteracting limit, as demonstrated in Fig.~\ref{fig:Udep}.  In the presence of a nonzero
Hund coupling $J_\mathrm{H}$, however, $V_{j\mathrm{eff}}$ shows a sign problem
as well, especially at intermediate interaction strengths and away from half-filling.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.45\textwidth]{udep.eps}
	\caption{
		Average sign of four different bases as a function of the intraorbital interaction strength $U$
		for $T=0.03$, $\lambda=0.25$, and $J_{\mathrm{H}}/U=0.15$.
		(a)-(d) correspond to the system with electron filling $n=2.5$, $2.0$, $3.5$, and $4.0$, respectively.
		Open squares and circles at $n=2$ and $4$ indicate that the self-consistent paramagnetic solution is Mott insulating.
}
	\label{fig:Udep}
\end{figure}

In order to further improve the average sign we introduce the basis $V_{j3/2\mathrm{BA}}$ in which $\tilde{c}^{}_{\frac{3}{2},+\frac{3}{2}}$ and $\tilde{c}^{}_{\frac{3}{2},-\frac{3}{2}}$ are mixed in the bonding-antibonding manner:
\begin{equation}
  \left(
  \begin{array}{c}
	  \tilde{a}^{}_{\frac{3}{2},+}\\
	  \tilde{a}^{}_{\frac{3}{2},-}
  \end{array}
  \right)
  = \frac{1}{\sqrt{2}}
  \left(
  \begin{array}{cccccc}
    1 & 1 \\
    1 & -1 \\
  \end{array}
  \right)
  \left(
  \begin{array}{c}
	  \tilde{c}^{}_{\frac{3}{2},+\frac{3}{2}}\\
	  \tilde{c}^{}_{\frac{3}{2},-\frac{3}{2}}
  \end{array}
  \right)~.\\
  \label{eqn:Uj3/2BA}
\end{equation}
A similar transformation was previously explored for an isolated trimer.\cite{2015PhRvB..92s5126S} It turns out that $V_{j3/2\mathrm{BA}}$ is superior to $V_{j\mathrm{eff}}$ for most parameter sets investigated.
In the noninteracting limit, $V_{j3/2\mathrm{BA}}$ (like $V_{j\mathrm{eff}}$) yields an average sign of unity (see Fig.~\ref{fig:Udep}).
As $U$ is switched on, $V_{j3/2\mathrm{BA}}$ results in a larger average sign than $V_{j\mathrm{eff}}$~.
In Fig.~\ref{fig:Udep}, the improvement of the average sign in $V_{j3/2\mathrm{BA}}$ is
 prominent in the itinerant regime with intermediate Coulomb interaction strength.
This improvement persists over the whole range of electron densities for $U=2$ as shown in Fig~\ref{fig:nldep}(a)~.
In the Mott insulating regime (marked by open squares and circles in Figs.~\ref{fig:nldep} and \ref{fig:Udep}),
on the other hand, the difference between the average signs of
$V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$ is smaller.  Since the CTHYB
is based on the expansion around the localized limit, those observations imply
that $V_{j3/2\mathrm{BA}}$ effectively prevents sign-problematic high-order
processes.

Note that neither $V_{j3/2\mathrm{BA}}$ nor $V_{j\mathrm{eff}}$ is an optimal basis for the non-SO-coupled model.
As we show in Figs.~\ref{fig:nldep}(b) and (c), at the
non-SO-coupled point ($\lambda=0$), $V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$
exhibit a sign problem, in contrast to $V_{\mathrm{os}}$~.
This demonstrates that Eq.~(\ref{eqn:kanamori}) includes dangerous interacting terms in the $V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$ basis.
One of these terms is the correlated-hopping (CH) term, which has the form
\begin{eqnarray}
	&&\frac{\sqrt{2}J_{\text{H}}}{3}\sum^{}_{s=\pm}\left[s(2\tilde{n}_{\frac{3}{2},s\frac{3}{2}} - \tilde{n}_{\frac{3}{2},\bar{s}\frac{3}{2}}- \tilde{n}_{\frac{3}{2},\bar{s}\frac{1}{2}})\tilde{c}^{\dagger}_{\frac{1}{2},s\frac{1}{2}}\tilde{c}^{}_{\frac{3}{2},s\frac{1}{2}} + \text{H.c.}\right]\nonumber\\
	\label{eqn:ch_Ujeff}
\end{eqnarray}
in $V_{j\mathrm{eff}}$ and
\begin{eqnarray}
	&&\frac{\sqrt{2}J_{\text{H}}}{6}\sum^{}_{s=\pm}\left[ s(\tilde{n}_{\frac{3}{2},+} + \tilde{n}_{\frac{3}{2},-}   - 2 \tilde{n}_{\frac{3}{2},\bar{s}\frac{1}{2}}) \tilde{c}^{\dagger}_{\frac{1}{2},s\frac{1}{2}}\tilde{c}^{}_{\frac{3}{2},s\frac{1}{2}}+\text{H.c.}\right]\nonumber\\
	&&+\frac{J_{\mathrm{H}}}{2}\left[(\tilde{n}_{\frac{1}{2},+\frac{1}{2}}-\tilde{n}_{\frac{1}{2},-\frac{1}{2}})\tilde{a}^{\dagger}_{\frac{3}{2},+}a^{}_{\frac{3}{2},-} + \text{H.c.}\right]
	\label{eqn:ch_Uj32BA}
\end{eqnarray}
in $V_{j3/2\mathrm{BA}}$~.
This CH term is the major source of the sign problem in the SO-coupled Hamiltonian.
Table~\ref{tab:mask} analyzes the effect of different terms in the local
 Hamiltonian on the average sign for both $V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$~.
There is a substantial drop in the average sign when the CH terms are introduced.

\begin{table}
	\centering
	\begin{tabular}{p{0.55\linewidth}p{0.20\linewidth}p{0.20\linewidth}}
		\hline\hline
		Terms& $\langle\mathrm{sign}\rangle_{V_{j\mathrm{eff}}}$ & $\langle\mathrm{sign}\rangle_{V_{j3/2\mathrm{BA}}}$\vspace{0.5pt}\\
		\hline
		Density-density (DD) & 1.00000 & 1.00000\\
		DD + spin flip (SF) & 1.00000 & 1.00000\\
		DD + four scattering (FS) & 0.99447(9) & 0.9516(3)\\
		DD + SF + pair hopping (PH) & 0.9760(2) & 0.9739(2)\\
		DD + PH & 0.9730(7) & 0.9705(2)\\
		DD + SF + PH + FS & 0.9308(3) & 0.9035(6)\\
		DD + \textit{correlated hopping} (CH) & 0.8212(9) & 0.8730(5)\\
		DD + CH + FS & 0.8061(6) & 0.7718(5)\\
		DD + PH + CH & 0.7728(6) & 0.8623(6)\\
		DD + SF + PH + CH & 0.7638(7) & 0.8814(5)\\
		DD + SF + PH + CH + FS & 0.6751(8) & 0.7406(4)\\
		\hline\hline
	\end{tabular}
	\caption{
		Average sign for the masked local Hamiltonians represented in $V_{j\mathrm{eff}}$ and $V_{j3/2\mathrm{BA}}$~.
		The results are for $n=3.5$, $T=0.03$, $\lambda=0.25$, $U=2.0$, and $J_\mathrm{H}/U$=0.15, and the
		self-consistent hybridization for the full local Hamiltonian.
	}
	\label{tab:mask}
\end{table}

In $V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$, there emerge other new terms involving four different flavors, which are not of the spin-flip type or pair-hopping type appearing in the Slater-Kanamori Hamiltonian in $V_{\mathrm{os}}$~.
Those terms are denoted as four-scattering terms in Table~\ref{tab:mask}~.
The average sign in both $V_{j3/2\mathrm{BA}}$ and $V_{j\mathrm{eff}}$ becomes even lower when the CH term is combined with the pair-hopping and the four-scattering terms.
The full interacting Hamiltonian in $V_{j\text{eff}}$ is explicitly written in Appendix~\ref{app:HintVjeff}.
Furthermore, we analyze the nature of the self-consistent solutions for the Hamiltonians without the CH or FS terms and we discuss the potential use of such {\it masked} Hamiltonians in Appendix~\ref{app:maskH}.

\begin{figure}[]
	\centering
	\includegraphics[width=0.45\textwidth]{spsa.eps}
	\caption{
		(a) Average sign of the positively and negatively perturbed parameter points during the SPSA search for the $SO(4)$ group ($n=3.5$, $T=0.03$, $\lambda=0.25$, $U=2.0$, and $J_{\mathrm{H}}/U=0.15$) as a function of iterations.
		The inset shows the convergence of the parameters defined in Eqs.~(\ref{eqn:ML}).
		(c) Landscape of the average sign in the $(\theta_1,\theta_2)$ subspace described by Eq.~(\ref{eqn:theta12space}).
		Green triangles represent the SPSA trajectory in the parameter space.
		(b) Enlarged plot of the SPSA trajectory of (c).
		The numeric labels in (b) denote iteration numbers during the SPSA search, and the arrow in (c) shows the direction in which the search proceeds.
	}
	\label{fig:SPSA}
\end{figure}

In what follows we will determine by the SPSA method the optimal basis in terms of the average sign for the
parameters with the most severe sign problem and show
that this optimal basis is nearly identical to the $V_{j3/2\mathrm{BA}}$ basis.
For that,
we search the basis space generated by the $SO(4)$ rotation group for $j=3/2$, whose $4\times 4$ matrix representation is denoted by $\mathbbm{M}$ and transforms the electron operators in $V_{j\mathrm{eff}}$ as follows:
\begin{gather}
	\tilde{a}^{}_{\frac{3}{2}} = \mathbbm{M}\cdot\tilde{c}^{}_{\frac{3}{2}}~,\\
	\tilde{a}^{\intercal}_{\frac{3}{2}} = \left(\tilde{a}^{}_{\frac{3}{2},1},\tilde{a}^{}_{\frac{3}{2},2},\tilde{a}^{}_{\frac{3}{2},3},\tilde{a}^{}_{\frac{3}{2},4}\right)~,\\
	\tilde{c}^{\intercal}_{\frac{3}{2}} = \left(\tilde{c}^{}_{\frac{3}{2},+\frac{1}{2}},\tilde{c}^{}_{\frac{3}{2},-\frac{1}{2}},\tilde{c}^{}_{\frac{3}{2},+\frac{3}{2}},\tilde{c}^{}_{\frac{3}{2},-\frac{3}{2}}\right)~.
	\label{}
\end{gather}
Since the off-diagonal hybridization function is a clear source of the severe sign problem as shown in the cases of $V_{\mathrm{os}}$ (Figs.~\ref{fig:nldep} and \ref{fig:Udep}), we exclude the mixing between the $j=1/2$ and $3/2$ subspaces to preserve the diagonal structure of the hybridization.
Without mixing between the $j=1/2$ and $3/2$ subspaces, one can fix the basis for the $j=1/2$ subspace using the rotational symmetry generated by the total angular momentum operator, $\mathbf{J}_{\mathrm{eff}}=\mathbf{L}_{\mathrm{eff}}+\mathbf{S}$ without loss of generality.

To parametrize the basis space, we introduce the isoclinic decomposition of $\mathbbm{M}$ as
	$\mathbbm{M} = \mathbbm{M}_{L}\mathbbm{M}_{R}$~,
where
\begin{align}
	&\mathbbm{M}_{L} {=}
	\left(
	\begin{array}{cccc}
		a & -b & -c & -d\\
		b & a & -d & c\\
		c & d & a & -b\\
		d & -c & b & a
	\end{array}
	\right),\hspace{2mm}
	\mathbbm{M}_{R}
	{=}
	\left(
	\begin{array}{cccc}
		p & -q & -r & -s\\
		q & p & s & -r\\
		r & -s & p & q\\
		s & r & -q & p
	\end{array}
	\right).
		\label{eqn:ML}
\end{align}
Here, $a^2+b^2+c^2+d^2 = 1$ and $p^2+q^2+r^2+s^2 = 1$~.
Under these constraints among $\{a,b,c,d\}$ and $\{p,q,r,s\}$ the dimension of the parameter space becomes $6$.

Figure~\ref{fig:SPSA} shows how the SPSA works while searching for the optimal basis in this parameter space.
The evolution of the average sign as a function of number of iterations
at both positively ($\langle\mathrm{sign}^+\rangle$) and negatively ($\langle\mathrm{sign}^-\rangle$) perturbed points in the parameter space is shown in Fig.~\ref{fig:SPSA} (a).
The average sign value converges to $\sim 0.74$ for $n=3.5$, $T=0.03$, $\lambda=0.25$, $U=2.0$, and $J_\mathrm{H}/U=0.15$~.
Within numerical accuracy, it is very close to the value of
$V_{j3/2\mathrm{BA}}$ defined in Eq.~(\ref{eqn:Uj3/2BA})~.  This shows that the
$V_{j3/2\mathrm{BA}}$ basis is at least near the local optimum in parameter
space.
The inset of Fig.~\ref{fig:SPSA}(a) illustrates the convergence of the parameters in Eqs.~(\ref{eqn:ML}).

Figure~\ref{fig:SPSA}(b) and (c) show the SPSA sequence in a small parameter subspace.
We introduce two parameters, $\theta_1$ and $\theta_2$, representing the restricted basis transformation
\begin{equation}
	\mathbbm{M}
	= \left(
	\begin{array}{cccc}
		\cos\theta_1 & \sin\theta_1 & 0 & 0\\
		-\sin\theta_1 & \cos\theta_1 & 0 & 0\\
		0 & 0 & \cos\theta_2 & \sin\theta_2\\
		0 & 0 & -\sin\theta_2 & \cos\theta_2\\
	\end{array}
	\right)~.
	\label{eqn:theta12space}
\end{equation}
 Figure~\ref{fig:SPSA}(c) shows that the landscape of the average sign is smooth, so that the
SPSA search based on the gradient approximation can successfully find the local optimum in that subspace.
Furthermore, $V_{j3/2\mathrm{BA}}$, corresponding to $\theta_1=0$ and $\theta_2=\pi/4$, is shown to be very close to the optimum found by the SPSA.
Figure~\ref{fig:SPSA}(b) plots the trajectory determined by the SPSA.

\begin{figure}[]
	\centering
	\includegraphics[width=0.45\textwidth]{Tdep.eps}
	\caption{
		Temperature dependence of the average sign
		for $T=0.03$, $\lambda=0.15$, $U=2.0$, and $J_{\mathrm{H}}/U=0.15$~.
		The electron filling for (a) and (b) is $3.5$ and $2.5$, respectively.
		The dashed (dotted) line corresponds to the fitting function of the form $\exp\left[-\beta\delta F + F_0\right]$ for the $V_{j3/2\mathrm{BA}}$ ($V_{j{\mathrm{eff}}}$) basis.
	}
	\label{fig:Tdep}
\end{figure}

We finally investigate the temperature scaling of the average sign for the different bases.
Figure~\ref{fig:Tdep} shows the exponential scaling $\langle \mathrm{sign}\rangle \sim \exp\left(-\beta\delta F + F_0\right)$
as a function of the inverse temperature $\beta$.
Here, $\delta F$ can be regarded as the free-energy difference from the auxiliary bosonic system without a sign problem, which determines the temperature scaling.
For both fillings, $n=3.5$ and $2.5$, the $V_{j3/2\mathrm{BA}}$ basis shows an improved temperature scaling exponent $\delta F$ and a comparable offset $F_0$~.

\section{Conclusion}
We have investigated the nature of the sign problem in the CTHYB for the SO-coupled multiorbital Hubbard model.
We found that the correlated hopping term that
appears in the $j_{\mathrm{eff}}$ basis is the major source of the sign problem,
 and we introduced a basis--the $j=3/2$ bonding-antibonding basis--which alleviates the effects of those terms.
By applying the stochastic optimization scheme, we found that
the $j=3/2$ bonding-antibonding basis is near a local optimum in the parameter space considered.
These results (i) provide useful guidelines for the choice
 of the single-particle basis in CTHYB simulations of SO-coupled systems,
(ii) introduce an algorithm to numerically determine the optimal basis, and
(iii) suggest using the intrinsic fermionic sign bound in a
systematic search of the full $SO(4)$ basis space
as a strategy for the further improvement of CTHYB for this class of systems.

In addition, our findings have relevance for descriptions beyond DMFT that include
nonlocal correlations. For example,
since the average sign is alleviated at the level of the partition function in $V_{j3/2\text{BA}}$,
we expect our results to be helpful for a reliable estimation of two-particle correlation functions,
which are the essential building blocks of diagrammatic extensions of DMFT.\cite{Rohringer2018}
For cluster extensions of DMFT, on the other hand, our numerical basis optimization scheme can be applied to an extended basis space including site (cluster) degrees of freedom.
A promising strategy might be to combine the $V_{j3/2\text{BA}}$ basis with a bonding-antibonding transformation that breaks loops on the cluster.\cite{2015PhRvB..92s5126S}

{\it Note added.} During the final stages of this work, a series of conceptually related studies appeared on  arXiv.\cite{Hangleiter2019,Torlai2019,Levy2019}

\begin{acknowledgments}
A.J.K. was supported by EPSRC through Grant No. EP/P003052/1.
P.W. acknowledges support from the Swiss National Science Foundation via NCCR Marvel.
R.V. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG)
through Grant No. VA117/15-1.
The computations were performed at the Centre for Scientific Computing (CSC)
in Frankfurt.
\end{acknowledgments}


\section{Supplemental Material}

\begin{table*}[h]
\centering
\begin{tabular}{@{}lllrllll@{}}
    \toprule
    Model                & $m_\chi$ [\si{GeV/c^2}]                       & $\sigma$ [\si{{cm}^2}]                    & \multicolumn{2}{c}{S2 ROI [PE]} & Observed & Expected signal  & Expected bg. \\ \colrule
    SI NR &  4 & \SI{6.4e-43}{} & {[}165,       & 271{]}       & 16    & $24.6 \pm 4.0$ & $3.7 \pm 0.6$ \\
    SI NR &  6 & \SI{4.8e-44}{} & {[}165,       & 276{]}       & 16    & $23.3 \pm 3.2$ & $4.0 \pm 0.6$ \\
    SI NR & 17 & \SI{1.8e-45}{} & {[}201,       & 578{]}       & 10    & $12.4 \pm 1.4$ & $5.9 \pm 0.8$ \\
    SI NR & 18 & \SI{5.2e-46}{} & {[}306,       & 1190{]}       & 3    & $4.4 \pm 0.5$ & $3.7 \pm 0.6$ \\
    SI-LM, $m_\phi = \SI{1}{MeV/c^2}$ & 4 & \SI{2.9e-38}{} & {[}165,       & 271{]}       & 16    & $24.8 \pm 4.2$ & $3.7 \pm 0.6$ \\
    SD (neutron-only) & 4 & \SI{9.8e-38}{} & {[}165,       & 271{]}       & 16    & $24.6 \pm 4.0$ & $3.7 \pm 0.6$ \\
    DM-electron scattering & 1 & $\SI{2.5e-40}{}$ & {[}165,       & 271{]}       & 16    & $24.7 \pm 4.1$ & $3.7 \pm 0.6$ \\ \colrule
    ALP & \SI{0.5e-6}{} & $g_\mathrm{ae} = \SI{7.7e-14}{}$ & {[}653,       & 897{]}       & 0    & $2.7 \pm 0.3$ & $0.5 \pm 0.1$ \\
    Dark photon & \SI{0.5e-6}{} & $\kappa = \SI{2.2e-16}{}$ & {[}653,       & 897{]}       & 0    & $2.7 \pm 0.3$ & $0.5 \pm 0.1$ \\
    \botrule
    \\
\end{tabular}
\caption{Results for specific DM models excluded at exactly 90\% confidence level. The columns show successively: DM model (abbreviated as in the text, LM denotes light mediator and ALP axion-like particle); DM mass $m_\chi$; interaction strength, specified by the cross section for elastic scattering of DM with a free nucleon (or an electron, for DM-electron scattering) at zero momentum transfer for most models, by the axio-electric coupling constant $g_\mathrm{ae}$ for ALPs, and the kinetic mixing parameter $\kappa$ for dark photons; the left and right bounds of the S2 ROI in photoelectrons (PE); the number of observed events in the ROI; the expected number of events from the excluded signal model; the expected number of background events from the modeled backgrounds. The final two columns also show aggregated systematic uncertainties.}
\label{tbl:intervals}
\end{table*}

\twocolumngrid

\clearpage
$\!$
\newpage
$\!$
\pagebreak
$\!$

\subsection{Details on the energy scale}

Figure \ref{fig:s2_resolution} shows the assumed response in the S2 region of this analysis to ER and NR signals at different energies and including all efficiencies. As stated in the text, for each dark matter model and mass, a specific S2 ROI is used, whose lower bound is never below 150 PE.

The average energy of events observed at a particular S2 depends on the underlying energy spectrum. For a flat energy spectrum, the majority of low-S2 events are underfluctuations, i.e.~events with an S2 below the median S2 response for their energy. Thus, the average energy of events from a flat spectrum observed at a given S2 is higher than a naive conversion based on the median S2 response of different energies would suggest. This is an instance of `Eddington bias', a systematic effect caused by statistical fluctuations. Similar effects are discussed in \cite{lux_dd}.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.49\textwidth]{s2_energy_resolution.pdf}
    \caption{S2 response after all efficiencies for ER (blue) and NR (red). Solid lines show the median S2 response to signals of a specific energy, considering S2s in [90, 3000] PE only. The shaded regions indicate the $\pm 1 \sigma$ spread of the response. No NR (ER) response is assumed below 0.7 (0.186) \si{keV}; the hypothetical S2 response to those regions is shown lightly shaded for illustration. Dotted lines show the mean energy of events observed at a particular S2, in a flat (in linear energy) spectrum with no $Q_y$ cutoff.}
    \label{fig:s2_resolution}
\end{figure}

\vfill\eject

\subsection{Details on the S2 width model}

The S2 waveform width is defined as
\begin{equation}
w = \big( \Delta t_{50\%} + \Delta t_r /2 \big) /2 ,
\end{equation}
where $\Delta t_{50\%}$ is the time between the first and third quartile of the S2 waveform, and $\Delta t_r$ the time between the start and median of the S2 waveform. This definition showed a slightly better correlation with depth than using $\Delta t_{50\%}$ alone.
We estimate the distribution of S2 widths using XENON1T's waveform simulator \cite{ap1}. Figure \ref{fig:width_validation} compares the S2 widths from simulation with those observed in neutron generator calibration data, showing that the median width agrees to within $\sim\!\SI{50}{ns}$. S2 widths of \ce{^{220}Rn}\xspace calibration data are in similarly good agreement. S1-tagged cathode events are broader than the model (by $\sim\! \SI{120}{ns}$), as expected from the higher diffusion constant and lower drift velocity in the low-field region at and immediately above the cathode.

The decrease of the mean S2 width with S2 at a fixed depth, clearly visible in Figures~\ref{fig:monster} and \ref{fig:cathode} for cathode events, is not discernible in Figure~\ref{fig:width_validation} due to a compensatory selection effect. At low S2, events with S1s are more common deep in the detector, where the charge loss to impurities and the photon detection efficiency are largest, but the average S2 width is also larger.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.49\textwidth]{wvsim_validation.pdf}
    \caption{First quartile (blue), median (green) and third quartile (red) of S2 widths in simulation (bands) and neutron generator calibration data (dots), as a function of depth $z$ (left) and S2 size (right), integrated over the other dimension (S2 on the left, $z$ on the right). The lower panels show the same, for the median only, with the nominal model subtracted. This figure uses events with S1s (so their $z$ is known) that pass event selections besides those on $z$ and S2 width; the multiple-S2 and SE-before-S2 cuts were relaxed to obtain additional statistics. The width of the bands denotes the uncertainty due to limited calibration data statistics (100 events/bin).}
    \label{fig:width_validation}
\end{figure}

\clearpage

\subsection{Details on the cathode background estimate}
Here we provide more details on our estimate of the background from decays on the cathode wires. As explained in the main text, most cathode events have detectable S1s even at very low S2. These cathode events with S1s are trivial to select, as they have a reconstructed depth $z$ within $\sim\! \SI{1}{cm}$ of the cathode, which sits at a drift length of $\SI{-96.9}{cm}$, with a spread that results from imperfectly corrected effects of the drift field inhomogeneity. In contrast, these events are homogeneously distributed in $(x,y)$, and their rate shows no significant time-dependence during SR1. We also observe \ce{^{214}Po} and \ce{^{210}Po} alpha decays at the cathode, and thus hypothesize that $\beta$ decays from \ce{^{214}Pb} and especially \ce{^{210}Pb} (due to its low $Q = \SI{64}{keV}$) contribute to the observed low-energy cathode events.

Figure~\ref{fig:cathode} illustrates the procedure for estimating the background from cathode events without S1s, after the width cut and inside the \SI{4}{GeV/c^2} NR analysis' S2 ROI, shown in \textbf{black} in figure \ref{fig:cathode}. First, we count cathode events with $S1 < \SI{20}{PE}$ that would appear in the S2 ROI if not for the width and z cuts. This region is shown in \textcolor{red}{red} in Figure~\ref{fig:cathode}. Next, we multiply this by the efficiency $\epsilon_w = (10.5 \pm 1.9) \%$ of the width cut on cathode events in the ROI, determined on cathode events with larger S1s ($20-\SI{200}{PE}$, not shown). Finally, we multiply by the ratio of cathode events without S1s and those with S1s $< \SI{20}{PE}$ in a control region at high S2 and high S2 width -- the count in the \textcolor{green}{green} divided by that in the \textcolor{orange}{orange} boxes in Figure~\ref{fig:cathode}. The count in the \textcolor{green}{green} box is corrected for (an insignificant) contamination by ordinary ER events. This is estimated by first counting events with non-cathode S1s in the control region -- the \textcolor{blue}{blue} box in Figure~\ref{fig:cathode} -- then scaling by one minus the S1 detection efficiency ($\epsilon_\mathrm{S1} > 90\%$) on ordinary ER events in the control region (determined as in \cite{ap1}). The cathode background estimate in the S2 ROI is thus:
\begin{equation}
\textcolor{red}{r} \cdot \epsilon_w \cdot [\textcolor{green}{g} - \textcolor{blue}{b} \cdot (1 - \epsilon_\mathrm{S1}) ] / \textcolor{orange}{o} ,
\end{equation}
where $\textcolor{red}{r}$, $\textcolor{green}{g}$, $\textcolor{blue}{b}$ and $\textcolor{orange}{o}$ are the event counts in the red, green, blue and orange boxes, respectively.

This procedure assumes that the ratio of cathode events with and without S1s is the same in the ROI as in the control region. However, events in the ROI have a lower S2, thus on average a lower energy, and are hence more likely to lack an S1 than events in the control region. Our estimate is therefore only a lower bound, and using it in the inference produces conservative limits.

We find $\textcolor{red}{r} = 46$, $\textcolor{orange}{o} = 336$, $\textcolor{green}{g} = 123$, $\textcolor{blue}{b} = 10$, leading to a cathode background estimate of $(1.62 \pm 0.46)$ events. The uncertainty includes statistical uncertainties on all the counts and the uncertainties on $\epsilon_w$ and $\epsilon_\mathrm{S1}$. The cathode background uncertainty is included in the aggregate background uncertainty described in the main text. We verified the estimate is stable well within its uncertainty under variations of the control region's S2 and S2 width bounds or the \SI{20}{PE} S1 threshold.

\onecolumngrid

\begin{figure}[b]
    \centering
    \includegraphics[width=\textwidth]{cathode_bg_estimate_summary.pdf}
    \caption{Illustration of the cathode background estimate. The scatterplots show (A) cathode events with $S1 < \SI{20}{PE}$, (B) events without S1s, and (C) events with $S1 < \SI{200}{PE}$ not from the cathode. Dot colors denote the S1 size in PE according to the colorbar in panel C. The colored boxes are used as described in the text to estimate the background in the DM ROIs after the width cut; the black box shows this region for \SI{4}{GeV/c^2} spin-independent NR DM. All events shown are from the search data, without the S2 width cut, and with the radial cut restricted to $R^2 < \SI{700}{cm}$ to avoid including misreconstructed events from the TPC walls in the blue and green boxes.}
    \label{fig:cathode}
\end{figure}
\twocolumngrid


\section{Introduction}
Entanglement in a few-atom system attracted enormous attention over last few decades \cite{leh1,leh2,ficek,entd1,entd2,lukin,eberly,tq_ent}.
Small qubit samples may form buildings blocks for even larger networks with huge potential applications for quantum technologies
\cite{zb_sc,QTBk,agarwal_2014,cirac,chk,wang,kimble,mek_ent,qt_les,qt_pl,qt_m}. Generally, a thorough description of the entanglement
creation in a two-atom system was given in \cite{fit}. From this point of view artificial atomic systems have been widely investigated as
well. Particularly, experimental realization of entanglement in two coupled charge qubits was performed in \cite{qcirc}. Entanglement of
two quantum dots inside a cavity injected with squeezed vacuum was predicted as well, in \cite{gao}. Ultra-strong dipole-dipole
interacting two-level superconducting flux qubits are naturally entangled through their corresponding environmental reservoir and maximum
coherence can be induced too \cite{kmj}. Furthermore, a pair of moderately dipole-dipole coupled and laser pumped two-level quantum
dots get maximal entangled via their environmental phonon thermostat which facilitates also the creation of a subradiant two-qubit state
\cite{ccim}. Recently, based on quantum dots systems, a relevant experimental realization of an interconnection among two qubits located
five meters apart from each other, via single photons, was reported in Ref.~\cite{qtsexp}.

Often to realize quantum states of matter or light one requires ground-state cooled individual or coupled quantum systems, respectively.
In this context, cooling of a quantum circuit via coupling to an independent or Dicke-like interacting multiqubit ensemble was demonstrated
in \cite{mmac}. The cooling of a nanomechanical resonator coupled to two interacting flux qubits via the corresponding subradiant Dicke
states was demonstrated as well, in Ref.~\cite{kxje}. A scheme for ground-state cooling of a mechanical resonator coupled to two
coupled quantum dots forming an effective $\Lambda$-type three-level structure was  presented in \cite{gxL}. Ground-state cooling of a
nanomechanical oscillator with $N$ spins was recently proposed in \cite{erem}.

Evidently, there is a considerable effort done to elucidate the relationship among
the entanglement and cooling phenomena. For instance, it was found that entanglement enhances cooling in microscopic quantum refrigerators
\cite{popq}. The atom-membrane cooling and entanglement using cavity electromagnetically induced transparency was investigated in
Ref.~\cite{ritsch}. Ground-state cooling enabled by critical coupling and dark entangled states was found in \cite{gray}. Further interesting
works on cooling or entanglement processes are presented in Refs.~\cite{c_ex,pt_mar,vitali,beige,popescu,jwang,aurelian}.

Here we shall take the opportunity of these advances and investigate the interconnection of the entanglement process in a coherently and
resonantly pumped dipole-dipole interacting two-qubit system, and the cooling effects of a boson mode with whom the qubits are longitudinally
coupled. Both the quantum subsystems are damped via their corresponding environmental reservoirs. We have found that the cooling of the
boson mode is accompanied by entanglement creation among the two-level qubit pair as long as the dipole-dipole frequency shift lies around
the boson mode frequency. Generally, the entanglement effect enhances during the cooling process while the amplitude of the applied coherent
field increases from zero. However, the maximum of the concurrence which was taken as entanglement measure, has slightly lower magnitudes
than those that would be obtained but in the absence of the boson mode coupling to the two qubits. Thus, the maximum boson mode cooling
efficiency realizes on the expense of the entanglement creation. Anyway, for lower bath temperatures, the concurrence will increase until
values which also can be reached without the boson mode coupling to the qubit pair. An intuitive explanation for the described effects is as
follows: Because we considered that the external coherent source frequency equals the qubit's transition frequency while its wave-vector
is perpendicular to the line connecting the two dipole-dipole interacting qubits, then a privileged way to excite the collective sample is via a
simultaneous absorption of a laser photon followed by a boson mode phonon absorption, respectively (see Fig.~\ref{fig-0ab}a). This is
because the dipole-dipole frequency splitting among the cooperative two-qubit states is close to the boson mode frequency. As a consequence,
the qubit pair excitation as a whole results in the boson mode cooling, whereas the two closely spaced emitters get entangled. Finally, these
effects occur for moderately weak external coherent driving fields which can avoid sample's deteriorations.
\begin{figure}[t]
\includegraphics[width=8cm]{fig0.eps}
\caption{\label{fig-0ab} (a) The two-qubit cooperative states $|\Phi_{i}\rangle$, where $i \in \{e,s,a,g\}$, stand for {\it excited}, {\it symmetrical},
{\it antisymmetrical} and {\it ground} Dicke states, respectively. A coherent laser source, with its wave-vector perpendicular to the line connecting
the two qubits, can excite the whole sample either directly, i.e. $|\Phi_{g}\rangle \leftrightarrow |\Phi_{s}\rangle \leftrightarrow |\Phi_{e}\rangle$,
or indirectly, that is $|\Phi_{g}\rangle \leftrightarrow |\Phi_{e}\rangle$. The latter path involves absorption of two photons simultaneously and is
less probable for weaker applied external fields - a situation considered here. Therefore, the first channel prevalent meaning that qubit's subsystem
excitation is taken place via absorption of a single laser photon of frequency $\omega_{L}$ followed, respectively, by a phonon absorption of
frequency $\omega$. A  process which leads to cooling of the single-mode boson field longitudinally coupled to the two qubits which become also
entangled. (b) The two-qubit dressed-states $|\Psi_{j}\rangle$ and the corresponding eigenvalues $\lambda_{j}$, $j \in \{1,2,3,4\}$, obtained
after diagonalization of the dipole-dipole coupled qubit-laser interaction Hamiltonian.}
\end{figure}

The article is organized as follows. In Sec. II we describe the system of interest together with the analytical approach used as well as
the boson mode features. In Sec. III we analyze the entanglement creation among the qubit pair and its relationship with the boson
mode cooling effects, respectively. We finalize the article with a summary given in Sec. IV.

\section{The master equation and equations of motion}
The system of interest consists of environmental vacuum mediated dipole-dipole coupled pair of identical two-level $\{|2\rangle_{j},|1\rangle_{j}\}$
quantum emitters $\{j \in q1,q2\}$, resonantly pumped by a coherent laser source. The laser wavelength is sufficiently bigger than the interqubit
spatial interval $|\vec r_{q_1q_2}|$, while the interparticle separation is larger than the linear size of the quantum emitter itself. Also, the laser
wave-vector is perpendicular to the line connecting the qubits. Hence, the qubits are in an equivalent position with respect to the driving field.
We have assumed that the transition frequencies of both qubits are equal and identical to the external pumping source frequency, respectively.
Furthermore, the two-level qubits interact with a single boson mode of frequency $\omega$ via a longitudinal coupling. The whole system
dampens via the interaction with the electromagnetic vacuum modes of the surrounding reservoir as well as through the corresponding boson
mode environmental thermostat, respectively. As an appropriate realistic model can be taken a coherently pumped two-qubit pair, formed of ions,
molecules, dimers or impurities as well as quantum dots or quantum wells, superconducting qubits etc., embedded in cavities, nanomechanical resonators
or quantum circuits ones etc. \cite{blat,les,keel,dimer,zol,sch0,schon,para,twqb_exp,monro,harm,qtsexp,qudits,discord,norri}.

In the following, we shall present the corresponding master equation describing the investigated model where all the involved parameters are included
properly.
\subsection{The master equation}
The master equation describing this global system in the Born-Markov approximations \cite{leh1,leh2,ficek,zb_sc,agarwal_2014} is given as follows
\begin{eqnarray}
\dot \rho &+&\frac{i}{\hbar}[\bar H,\rho]=-\sum_{\{j,l\} \in\{q1,q2\}}\bigl(\gamma_{jl}[S^{+}_{j},S^{-}_{l}\rho] + H.c.\bigr) \nonumber \\
&-& \frac{\kappa}{2}(1+\bar n)[b^{\dagger},b\rho] - \frac{\kappa}{2}\bar n[b,b^{\dagger}\rho] + H.c.,
\label{ros}
\end{eqnarray}
where an overdot denotes differentiation with respect to time. Here, the qubit operators $S^{+}_{j} = |2\rangle_{j}{}_{j}\langle 1|$,
$S^{-}_{j}=[S^{+}_{j}]^{\dagger}$ and $S^{(j)}_{z}=(|2\rangle_{j}{}_{j}\langle 2| - |1\rangle_{j}{}_{j}\langle 1|)/2$ obey the commutation relations:
$[S^{+}_{j},S^{-}_{l}]=2S^{(j)}_{z}\delta_{jl}$ whereas $[S^{(j)}_{z},S^{\pm}_{l}]=\pm S^{\pm}_{j}\delta_{jl}$. The respective boson mode creation,
$b^{\dagger}$, and annihilation, $b$, operators satisfy the following commutation relations: $[b,b^{\dagger}]=1$ and
$[b,b]=[b^{\dagger},b^{\dagger}]=0$. In Eq.~(\ref{ros}), the Hamiltonian characterizing the corresponding coherent quantum dynamics is
$\bar H=H + H_{i}$, where
\begin{eqnarray}
H = \hbar\omega b^{\dagger}b + \hbar g\sum_{j\in\{q1,q2\}}S^{(j)}_{z}(b+b^{\dagger}),
\label{hh}
\end{eqnarray}
and
\begin{eqnarray}
H_{i} = \hbar\Omega_{dd}\sum_{j \not = l\in \{q1,q2\}}S^{+}_{j}S^{-}_{l} + \hbar\Omega \sum_{j \in \{q1,q2\}}\bigl (S^{+}_{j} + S^{-}_{j}\bigr).
\label{hi}
\end{eqnarray}
Here, $g$ is the qubit-boson-mode coupling strength whereas $\Omega$ denotes the standard Rabi frequency, both assumed to be identical for each
qubit, respectively. $\gamma_{q1q1}=\gamma_{q2q2}=\gamma/2$ is the single-qubit spontaneous decay rate, while $\gamma_{q1q2}=\gamma_{q2q1}
=\gamma \chi_{r}/2$ describes the radiative coupling among the two-level qubits and $\Omega_{dd}$ corresponds to the dipole-dipole interaction potential,
respectively. The radiative coupling $\chi_{r}$ goes to zero (unity) for larger (smaller) interparticle separations $|\vec r_{q_1q_2}|$ in comparison to the
photon emission wavelength. Correspondingly, $\Omega_{dd}$ tends to zero or to the static dipole-dipole interaction potential. Finally, $\kappa$ is the
damping rate of the boson mode, while $\bar n$ gives its mean thermal phonon number corresponding to the frequency $\omega$ and environmental
temperature $T$.

For our further purposes, we diagonalize the Hamiltonian $(\ref{hi})$ describing the dipole-dipole coupled qubit pair interacting as well with an externally
applied coherent laser field, using the following two-qubit bare states: $|2_{q1}2_{q2}\rangle$, $|2_{q1}1_{q2}\rangle$, $|1_{q1}2_{q2}\rangle$ and
$|1_{q1}1_{q2}\rangle$. Hence, we arrive at the corresponding cooperative two-qubit eigenfunctions:
\begin{eqnarray}
|\Psi_{4}\rangle &=&-\bar a \{|2_{q1}2_{q2}\rangle + |1_{q1}1_{q2}\rangle\} + \bar b\{|2_{q1}1_{q2}\rangle + |1_{q1}2_{q2}\rangle\}, \nonumber \\
|\Psi_{3}\rangle &=&-\bar c \{|2_{q1}2_{q2}\rangle + |1_{q1}1_{q2}\rangle\} + \bar d\{|2_{q1}1_{q2}\rangle + |1_{q1}2_{q2}\rangle\}, \nonumber \\
|\Psi_{2}\rangle &=&\frac{1}{\sqrt{2}}\{|2_{q1}1_{q2}\rangle - |1_{q1}2_{q2}\rangle\}, \nonumber \\
|\Psi_{1}\rangle &=&\frac{1}{\sqrt{2}}\{|2_{q1}2_{q2}\rangle - |1_{q1}1_{q2}\rangle\}. \label{drst}
\end{eqnarray}
Here
\begin{eqnarray}
\bar a = \frac{(\Omega_{dd}-\lambda_{4})/\sqrt{2}}{\sqrt{(\Omega_{dd}-\lambda_{4})^{2}+4\Omega^{2}}}, ~~~
\bar b = \sqrt{\frac{2\Omega^{2}}{(\Omega_{dd}-\lambda_{4})^{2}+4\Omega^{2}}}, \nonumber \\
\bar c = \frac{(\Omega_{dd}-\lambda_{3})/\sqrt{2}}{\sqrt{(\Omega_{dd}-\lambda_{3})^{2}+4\Omega^{2}}}, ~~~
\bar d = \sqrt{\frac{2\Omega^{2}}{(\Omega_{dd}-\lambda_{3})^{2}+4\Omega^{2}}}, \nonumber
\end{eqnarray}
with
\begin{eqnarray*}
\lambda_{4} &=& \bigl(\Omega_{dd} - \sqrt{\Omega^{2}_{dd}+16\Omega^{2}}\bigr)/2,  \nonumber \\
\lambda_{3} &=& \bigl(\Omega_{dd} + \sqrt{\Omega^{2}_{dd}+16\Omega^{2}}\bigr)/2,
\end{eqnarray*}
whereas other eigenvalues are
\begin{eqnarray*}
\lambda_{2} = -\Omega_{dd}, ~~{\rm and}~~ \lambda_{1} = 0,
\end{eqnarray*}
respectively, see also Fig.~{\ref{fig-0ab}(b)}. Substituting the two-qubit dressed-state transformation (\ref{drst}) in the master equation (\ref{ros}),
while keeping the slowly varying terms only by assuming that $\omega > g$ with $\omega \approx \lambda_{3}$ as well as $\Omega_{dd} \gg \gamma$
and $|\lambda_{4}| \ll |\lambda_{2}|$, one arrives at the following main equation governing the quantum dynamics of the examined system
\begin{eqnarray}
\dot \rho(t) &+& \frac{i}{\hbar}[\bar H_{0},\rho] = -\frac{\gamma}{2}(1+\chi_{r})\biggl(2[\bar c\bar d R_{44} + \bar a\bar b R_{33} +
\frac{\bar c}{2\sqrt{2}}(R_{41}-R_{14}),\{4(\bar c\bar d R_{44} + \bar a\bar b R_{33}) \nonumber \\
&+& \sqrt{2}\bar c(R_{14} - R_{41})\}\rho] + 2(\bar a \bar d +\bar b \bar c)^{2}\{[R_{34},R_{43}\rho] +[R_{43},R_{34}\rho]\}
+ \bar a^{2}\{[R_{13},R_{31}\rho] \nonumber \\
&+& [R_{31},R_{13}\rho]\} - \sqrt{2}\bar a(\bar a\bar d + \bar b\bar c)\{[R_{43},R_{31}\rho] + [R_{13},R_{34}\rho] - [R_{34},R_{13}\rho]
- [R_{31},R_{43}\rho]\} \biggr ) \nonumber \\
&-& \frac{\gamma}{2}(1-\chi_{r})\biggl(\bar b^{2}\{[R_{32},R_{23}\rho] + [R_{23},R_{32}\rho]\}+\bar d^{2}\{[R_{24},R_{42}\rho]
+ [R_{42},R_{24}\rho]\} \nonumber \\
&+& \frac{1}{2}\{[R_{12},R_{21}\rho]+[R_{21},R_{12}\rho]\} - \frac{\bar d}{\sqrt{2}}\{ [R_{42},R_{21}\rho] + [R_{12},R_{24}\rho]
- [R_{24},R_{12}\rho] \nonumber \\
&-& [R_{21},R_{42}\rho]\} \biggr) -
\frac{\kappa}{2}(1+\bar n)[b^{\dagger},b\rho] - \frac{\kappa}{2}\bar n[b,b^{\dagger}\rho] + H.c.. \label{dmeq}
\end{eqnarray}
Here
\begin{eqnarray}
\bar H_{0} = \hbar \lambda_{4}R_{44} - \hbar \delta b^{\dagger}b - \hbar \bar g(R_{31}b + b^{\dagger}R_{13}), \label{h0}
\end{eqnarray}
where $\delta=\lambda_{3}-\omega$, whereas $\bar g =\sqrt{2}g\bar c$, see also Fig.~{\ref{fig-0ab}(b)}. The resulting two-qubit dressed-state operators
which enter in Eq.~(\ref{dmeq}) are defined as follows: $R_{\alpha\beta}=|\Psi_{\alpha}\rangle \langle \Psi_{\beta}|$, $\{\alpha,\beta \in 1, \cdots, 4\}$,
and satisfy the standard commutation relations $[R_{\alpha\beta},R_{\beta'\alpha'}]=R_{\alpha\alpha'}\delta_{\beta\beta'}$ -
$R_{\beta'\beta}\delta_{\alpha'\alpha}$. Note that $|\lambda_{4}| \ll |\lambda_{2}|$ means also that we deal with rather weaker applied laser fields, i.e.
$\Omega/\Omega_{dd} \ll 1$ or the Rabi frequency $\Omega$ is of the order of few $\gamma$'s or even less.

\subsection{The equations of motion}
Using the Master equation (\ref{dmeq}), one can obtain the following equations of motion describing the combined {\it laser pumped qubit pair plus boson
mode} sample where the corresponding pumping and damping effects are properly taken into account:
\begin{eqnarray}
\dot P^{(0)}_{n} &=& i\bar g(P^{(4)}_{n} - P^{(6)}_{n}) - \kappa\bar n\bigl((n+1)P^{(0)}_{n} - nP^{(0)}_{n-1}\bigr) \nonumber \\
&-& \kappa(1+\bar n)\bigl(nP^{(0)}_{n} - (n+1)P^{(0)}_{n+1}\bigr), \nonumber \\
\dot P^{(1)}_{n} &=& i\bar g P^{(4)}_{n} - \kappa\bar n\bigl((n+1)P^{(1)}_{n} - nP^{(1)}_{n-1}\bigr) \nonumber \\
&-&\kappa(1+\bar n)\bigl(nP^{(1)}_{n} - (n+1)P^{(1)}_{n+1}\bigr) + \gamma^{(1)}_{0}P^{(0)}_{n} \nonumber \\
&-& \gamma^{(1)}_{1}P^{(1)}_{n} - \gamma^{(1)}_{2}P^{(2)}_{n} - \gamma^{(1)}_{3}P^{(3)}_{n} +
\gamma^{(1)}_{11}P^{(11)}_{n}, \nonumber
\end{eqnarray}
\begin{eqnarray}
\dot P^{(2)}_{n} &=& - \kappa(1+\bar n)\bigl(nP^{(2)}_{n} - (n+1)P^{(2)}_{n+1}\bigr) \nonumber \\
&-&  \kappa\bar n\bigl((n+1)P^{(2)}_{n}- nP^{(2)}_{n-1}\bigr) + \gamma^{(2)}_{0}P^{(0)}_{n} \nonumber \\
&+& \gamma^{(2)}_{1}P^{(1)}_{n} - \gamma^{(2)}_{2}P^{(2)}_{n} + \gamma^{(2)}_{3}P^{(3)}_{n} - \gamma^{(2)}_{11}P^{(11)}_{n}, \nonumber \\
\dot P^{(3)}_{n} &=& -i\bar g P^{(6)}_{n} + \gamma^{(3)}_{0}P^{(0)}_{n} + \gamma^{(3)}_{1}P^{(1)}_{n} -  \gamma^{(3)}_{2}P^{(2)}_{n}  \nonumber  \\
&-& \gamma^{(3)}_{3}P^{(3)}_{n} - \gamma^{(3)}_{11}P^{(11)}_{n} -\kappa\bar n\bigl((n+1)P^{(3)}_{n}  \nonumber \\
&-& nP^{(3)}_{n-1}\bigr) - \kappa(1+\bar n)\bigl(nP^{(3)}_{n}-(n+1)P^{(3)}_{n+1}\bigr), \nonumber \\
\dot P^{(4)}_{n} &=& -i\delta P^{(5)}_{n} +2 i\bar g n (P^{(1)}_{n} - P^{(3)}_{n-1} ) - \kappa(1+\bar n)\bigl(2P^{(6)}_{n} \nonumber \\
&+&(2n-1)P^{(4)}_{n} -2 (n+1)P^{(4)}_{n+1}\bigr)/2 + \kappa\bar n\bigl(2nP^{(4)}_{n-1} \nonumber \\
&-& (2n+1)P^{(4)}_{n} \bigr)/2 - \gamma^{(4)}_{4}P^{(4)}_{n} + \gamma^{(4)}_{8}P^{(8)}_{n}, \nonumber \\
\dot P^{(5)}_{n} &=& -i\delta P^{(4)}_{n}  - \kappa(1+\bar n)\bigl(2P^{(7)}_{n} + (2n-1)P^{(5)}_{n} \nonumber \\
&-&2 (n+1)P^{(5)}_{n+1}\bigr)/2 - \kappa\bar n\bigl( (2n+1)P^{(5)}_{n} \nonumber \\
&-& 2nP^{(5)}_{n-1} \bigr)/2 - \gamma^{(5)}_{5}P^{(5)}_{n} + \gamma^{(5)}_{9}P^{(9)}_{n}, \nonumber \\
\dot P^{(6)}_{n} &=& -i\delta P^{(7)}_{n} +2 i\bar g(n+1)(P^{(1)}_{n+1} - P^{(3)}_{n}) - \kappa(1+\bar n) \nonumber \\
&\times& \bigl((2n+1)P^{(6)}_{n} -2 (n+1)P^{(6)}_{n+1}\bigr)/2 + \kappa\bar n\bigl(2nP^{(6)}_{n-1} \nonumber \\
&-& (2n+3)P^{(6)}_{n} + 2P^{(4)}_{n}\bigr)/2 - \gamma^{(6)}_{6}P^{(6)}_{n} + \gamma^{(6)}_{12}P^{(12)}_{n}, \nonumber \\
\dot P^{(7)}_{n} &=& -i\delta P^{(6)}_{n} - \kappa(1+\bar n)\bigl((2n+1)P^{(7)}_{n} -2 (n+1) \nonumber \\
&\times&P^{(7)}_{n+1}\bigr)/2 -\kappa\bar n\bigl( (2n+3)P^{(7)}_{n}- 2nP^{(7)}_{n-1}  \nonumber \\
&-&2P^{(5)}_{n} \bigr)/2  -\gamma^{(7)}_{7}P^{(7)}_{n} + \gamma^{(7)}_{13}P^{(13)}_{n}, \nonumber \\
\dot P^{(8)}_{n} &=&i(\lambda_{4}-\delta)P^{(9)}_{n} + i\bar g n P^{(11)}_{n} - \kappa\bar n\bigl( (2n+1)P^{(8)}_{n} \nonumber \\
&-&2nP^{(8)}_{n-1} \bigr)/2 - \kappa(1+\bar n)\bigl( (2n-1)P^{(8)}_{n} + 2P^{(12)}_{n}\nonumber \\
&-&2 (n+1)P^{(8)}_{n+1}\bigr)/2 + \gamma^{(8)}_{4}P^{(4)}_{n} - \gamma^{(8)}_{8}P^{(8)}_{n}, \nonumber \\
\dot P^{(9)}_{n} &=&i(\lambda_{4}-\delta)P^{(8)}_{n} + i\bar g n P^{(10)}_{n} - \kappa\bar n\bigl( (2n+1)P^{(9)}_{n} \nonumber \\
&-&2nP^{(9)}_{n-1} \bigr)/2 - \kappa(1+\bar n)\bigl((2n-1)P^{(9)}_{n}+2P^{(13)}_{n}\nonumber \\
&-&2(n+1)P^{(9)}_{n+1}\bigr)/2 + \gamma^{(9)}_{5}P^{(5)}_{n} - \gamma^{(9)}_{9}P^{(9)}_{n}, \nonumber \\
\dot P^{(10)}_{n} &=&i\lambda_{4}P^{(11)}_{n}+ i\bar g P^{(9)}_{n} - \kappa\bar n\bigl((n+1)P^{(10)}_{n} - nP^{(10)}_{n-1}\bigr) \nonumber \\
&-&\kappa(1+\bar n)\bigl(nP^{(10)}_{n} - (n+1)P^{(10)}_{n+1}\bigr) - \gamma^{(10)}_{10}P^{(10)}_{n}, \nonumber \\
\dot P^{(11)}_{n} &=&i\lambda_{4}P^{(10)}_{n}+ i\bar g P^{(8)}_{n} - \kappa\bar n\bigl((n+1)P^{(11)}_{n} - nP^{(11)}_{n-1}\bigr) \nonumber \\
&-&\kappa(1+\bar n)\bigl(nP^{(11)}_{n} - (n+1)P^{(11)}_{n+1}\bigr) + \gamma^{(11)}_{0}P^{(0)}_{n} \nonumber
\end{eqnarray}
\begin{eqnarray}
&-& \gamma^{(11)}_{1}P^{(1)}_{n} - \gamma^{(11)}_{2}P^{(2)}_{n} - \gamma^{(11)}_{3}P^{(3)}_{n} - \gamma^{(11)}_{11}P^{(11)}_{n}, \nonumber \\
\dot P^{(12)}_{n}&=& i(\lambda_{4}-\delta)P^{(13)}_{n} + i\bar g(n+1) P^{(11)}_{n+1} + \gamma^{(12)}_{6}P^{(6)}_{n} \nonumber \\
&-& \kappa(1+\bar n)\bigl((2n+1)P^{(12)}_{n}- 2(n+1)P^{(12)}_{n+1} \bigr)/2  \nonumber \\
&-& \kappa\bar n\bigl((2n+3) P^{(12)}_{n} -2nP^{(12)}_{n-1} - 2P^{(8)}_{n}\bigr)/2  \nonumber \\
&-& \gamma^{(12)}_{12}P^{(12)}_{n}, \nonumber \\
\dot P^{(13)}_{n}&=& i(\lambda_{4}-\delta)P^{(12)}_{n} + i\bar g(n+1) P^{(10)}_{n+1} + \gamma^{(13)}_{7}P^{(7)}_{n} \nonumber \\
&-& \kappa(1+\bar n)\bigl((2n+1)P^{(13)}_{n}- 2(n+1)P^{(13)}_{n+1} \bigr)/2  \nonumber \\
&-& \kappa\bar n\bigl((2n+3) P^{(13)}_{n} -2nP^{(13)}_{n-1} - 2P^{(9)}_{n}\bigr)/2  \nonumber \\
&-& \gamma^{(13)}_{13}P^{(13)}_{n}. \label{eqsm}
\end{eqnarray}
The system of equations (\ref{eqsm}) can be easily obtained if one first get the equations of motion for the variables
$\rho_{\alpha\beta}=\langle \alpha|\rho|\beta\rangle$, $\{\alpha,\beta \in 1, \cdots, 4\}$, see also \cite{quang}, using the Master Equation (\ref{dmeq}),
namely, $\rho^{(0)}=\rho_{11}+\rho_{22}+\rho_{33}+\rho_{44}$, $\rho^{(1)}=\rho_{11}$, $\rho^{(2)}=\rho_{22}$,
$\rho^{(3)}=\rho_{33}$, $\rho^{(4)}=b^{\dagger}\rho_{31}-\rho_{13}b$, $\rho^{(5)}=b^{\dagger}\rho_{31} +\rho_{13}b$,
$\rho^{(6)}=\rho_{31}b^{\dagger} - b\rho_{13}$, $\rho^{(7)}=\rho_{31}b^{\dagger} + b\rho_{13}$,
$\rho^{(8)}=b^{\dagger}\rho_{34}-\rho_{43}b$, $\rho^{(9)}=b^{\dagger}\rho_{34} + \rho_{43}b$,
$\rho^{(10)}=\rho_{14}-\rho_{41}$, $\rho^{(11)}=\rho_{14}+\rho_{41}$, $\rho^{(12)}=\rho_{34}b^{\dagger} - b\rho_{43}$,
and $\rho^{(13)}=\rho_{34}b^{\dagger} + b\rho_{43}$. The projection on the Fock states $|n\rangle$, i.e.,
$P^{(i)}_{n}=\langle n|\rho^{(i)}|n\rangle$, $\{i \in 0, \cdots, 13\}$, with $n \in \{0, \infty\}$, will lead us to Eqs.~(\ref{eqsm}).
The corresponding decay rates are given in the Appendix A.
\begin{figure}[t]
\includegraphics[width=4.56cm]{fig1.eps}
\includegraphics[width=4.58cm]{fig2.eps}
\begin{picture}(0,0)
\put(-160,70){(a)}
\put(-25,70){(b)}
\end{picture}
\caption{\label{fig-1ab} (a) The steady-state behavior of the scaled mean phonon number $\langle b^{\dagger}b\rangle/\bar n$
as well as (b) the steady-state behavior of the second-order phonon-phonon correlation function $g^{(2)}(0)$ versus $\Omega/\Omega_{dd}$,
respectively. The involved parameters are: $g/\gamma=2$, $\Omega_{dd}/\gamma=28$, $\omega/\gamma=30$, $\chi_{r}=0.98$,
$\bar n=20$ and $\kappa/\gamma=10^{-3}$.}
\end{figure}

Generally, in order to solve the infinite system of equations (\ref{eqsm}), one truncates it at a certain maximum value $n=n_{max}$ so that a further
increase of its value, i.e. $n_{max}$, does not modify the obtained results. As a consequence, the steady-state mean phonon number is expressed as:
\begin{eqnarray}
\langle b^{\dagger} b \rangle = \sum^{n_{max}}_{n=0}nP_{n}^{(0)}, \label{nfm}
\end{eqnarray}
with
\begin{eqnarray}
\sum^{n_{max}}_{n=0}P_{n}^{(0)}=1, \label{nrm}
\end{eqnarray}
while its steady-state second-order phonon-phonon correlation function is defined in the usual way \cite{glauber}, namely,
\begin{eqnarray}
g^{(2)}_{b}(0) &=& \frac{ \langle b^{\dagger}b^{\dagger}bb\rangle }{\langle b^{\dagger} b\rangle^{2} } \nonumber \\
&=&\frac{1}{\langle b^{\dagger}b\rangle ^{2}}\sum^{n_{max}}_{n=0}n(n-1)P_{n}^{(0)}. \label{gg2}
\end{eqnarray}

Based on Eqs.~(\ref{eqsm},\ref{nfm},\ref{nrm},\ref{gg2}), Figure \ref{fig-1ab}(a) shows the cooling of the boson mode while the pumping parameter
is being varying demonstrating an efficient cooling scheme. Respectively, Figure \ref{fig-1ab}(b) depicts the second-order phonon-phonon
correlation function, during the cooling process, demonstrating super-Poissonian phonon statistics, i.e. $g^{(2)}(0) > 2$, with only few phonons.
The cooling mechanism occurring in this system, when $\Omega \ll \Omega_{dd}$, can be intuitively understood if one refers to the two-qubit Dicke
states \cite{leh1,leh2,ficek}, namely, $|\Phi_{e}\rangle=|2_{q1}2_{q2}\rangle$, $|\Phi_{s}\rangle$=$\{|2_{q1}1_{q2}\rangle +
|1_{q1}2_{q2}\rangle\}/\sqrt{2}$, $|\Phi_{a}\rangle$=$\{|2_{q1}1_{q2}\rangle - |1_{q1}2_{q2}\rangle\}/\sqrt{2}$, and $|\Phi_{g}\rangle=
|1_{q1}1_{q2}\rangle$, see Fig.~\ref{fig-0ab}(a). When the external field frequency is in resonance with that of the qubit one, while its
wave-vector is perpendicular to the line connecting the two qubits, then the only way to laser excite the two-qubit sample is either via
$|\Phi_{g}\rangle \to |\Phi_{s}\rangle \to |\Phi_{e}\rangle$ or $|\Phi_{g}\rangle \to |\Phi_{e}\rangle$, respectively. The latter path involve
two-photon processes which are less probable for weaker external driving fields. Therefore, the first channel, i.e.
$|\Phi_{g}\rangle \to |\Phi_{s}\rangle \to |\Phi_{e}\rangle$, becomes active and involves an available phonon at a particular frequency from
simple reasons since the driving external coherent field is in resonance with the qubit's transition frequency, that is $\omega_{sg} \approx \omega_{L}
+ \omega$, see Fig.~\ref{fig-0ab}(a). Thus, one can conclude that the two-qubit system absorbs an external laser photon followed by a boson mode
phonon absorption in order to reach the symmetrical Dicke state $|\Phi_{s}\rangle$, and these processes lead to phonon cooling effects, respectively.
Figure \ref{fig-2ab}(a) demonstrates this statement in the sense that the population of the symmetrical two-qubit state, i.e. $|\Phi_{s}\rangle$,
increases in the presence of the boson mode, coupled to the two qubits, while compared to the case of its absence. In the dressed-state picture,
depicted in Fig.~\ref{fig-0ab}(b), cooling occurs evidently when the population in the dressed-state $|\Psi_{1}\rangle$ is larger than that residing
in the state $|\Psi_{3}\rangle$, respectively. If one inspects the explicit forms of these states, see Exps.~(\ref{drst}), then one can observe that
laser excitation involves directly all the two-qubit Dicke states, excepting the antisymmetrical one, i.e. $|\Phi_{a}\rangle$.
\begin{figure}[t]
\includegraphics[width=4.56cm]{fig3.eps}
\includegraphics[width=4.58cm]{fig4.eps}
\begin{picture}(0,0)
\put(-240,70){(a)}
\put(-25,70){(b)}
\end{picture}
\caption{\label{fig-2ab} (a) The steady-state behaviors of the population in the symmetrical two-qubit Dicke state
$|\Phi_{s}\rangle$=$\{|2_{q1}1_{q2}\rangle + |1_{q1}2_{q2}\rangle\}/\sqrt{2}$, i.e.
$\Pi_{s} =\langle|\Phi_{s}\rangle \langle \Phi_{s}|\rangle$, as a function of scaled pumping parameter $\Omega/\Omega_{dd}$.
(b) The corresponding behaviors but for the two-qubit concurrence $C$. In these plots the solid lines are plotted for $g\not =0$
while the dashed curves describe the same results but with $g=0$, respectively. All other parameters are the same as for
Fig.~(\ref{fig-1ab}).}
\end{figure}

Until now we have focused mainly on the boson mode properties. In the following Section using the analytical approach developed here, we shall
investigate the entanglement creation within the qubit subsystem and emphasize its connection to the boson mode cooling phenomena, respectively.
\section{Entanglement of the two-qubit system coupled with a single-mode boson field}
The entanglement and its definition is certainly a main topic within quantum computation theory. In this context, one of the widely accepted measures
of entanglement for a two qubit system is the concurrence, $C$, \cite{entd1,entd2}. Particularly, for a mixed state of qubits $\{q_{1},q_{2}\}$ with
density matrix $\tilde \rho_{q_{1}q_{2}}$, it is defined as
\begin{equation}
C=\rm{max}\{0,s_{1}-\sum^{4}_{\xi=2}s_{\xi}\}. \label{ccr}
\end{equation}
The quantities $s_{\xi}$, $\{ \xi \in 1,\cdots,4\}$, are the square roots of the eigenvalues of the following matrix product
\begin{equation}
Q=\tilde \rho_{q_{1}q_{2}}(\sigma_{q_{1}y}\otimes\sigma_{q_{2}y})\tilde \rho^{\ast}_{q_{1}q_{2}}(\sigma_{q_{1}y}\otimes\sigma_{q_{2}y}),
\label{mp}
\end{equation}
and, importantly, in descending order. Here, $\tilde \rho^{\ast}_{q_{1}q_{2}}$ denotes complex conjugation of $\tilde \rho_{q_{1}q_{2}}$, and
$\sigma_{jy}$ are Pauli matrices for the two-level systems ($j \in \{q_{1},q_{2}\}$). The values of the concurrence range from zero for an unentangled
state to unity for a maximally entangled two-particle state \cite{entd1,entd2}. The density matrix $\tilde \rho_{q_1q_2}$ can be represented in the basis
$|2_{q1}2_{q2}\rangle$, $|2_{q1}1_{q2}\rangle$, $|1_{q1}2_{q2}\rangle$ and $|1_{q1}1_{q2}\rangle$, which is symmetric under the exchange of the
sub-systems \cite{entd1,entd2,wang,mek_ent}. Hence, its elements are given as follows:
\begin{eqnarray}
\tilde \rho_{q_1q_2} = \left( \begin{array}{cccc}
\tilde \rho_{11} & \tilde \rho_{12} & \tilde \rho_{13} & \tilde \rho_{14}\\
\tilde \rho_{21} & \tilde \rho_{22} & \tilde \rho_{23} & \tilde \rho_{24}\\
\tilde \rho_{31} & \tilde \rho_{32} & \tilde \rho_{33} & \tilde \rho_{34} \\
\tilde \rho_{41} & \tilde \rho_{42} & \tilde \rho_{43} & \tilde \rho_{44}
\end{array} \right),
\label{rab}
\end{eqnarray}
where
\begin{eqnarray}
\tilde \rho_{11} &=& \frac{1}{4}(1+\rho_{11}-\rho_{22}) - \frac{\Omega_{dd}}{4\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}}(\rho_{33}-\rho_{44})
\nonumber \\
&-&\frac{\bar a}{\sqrt{2}}(\rho_{41}+\rho_{14}), \nonumber \\
\tilde \rho_{12} &=&\frac{\Omega}{\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}}(\rho_{33}-\rho_{44}) - \frac{\bar c}{\sqrt{2}}\rho_{41}, \nonumber \\
\tilde \rho_{13} &=&\tilde \rho_{12}, \nonumber \\
\tilde \rho_{14} &=&\frac{1}{4}(1+\frac{\Omega_{dd}}{\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}})\rho_{44} - \frac{\bar a}{\sqrt{2}}(\rho_{41}-\rho_{14})
\nonumber \\
&+& \frac{1}{4}(1-\frac{\Omega_{dd}}{\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}})\rho_{33} - \frac{1}{2}\rho_{11}, \nonumber \\
\tilde \rho_{21} &=&(\tilde \rho_{12})^{\dagger}, \nonumber \\
\tilde \rho_{22} &=& \frac{1}{4}(1+\rho_{22}-\rho_{11}) + \frac{\Omega_{dd}}{4\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}}(\rho_{33}-\rho_{44}),
\nonumber \\
\tilde \rho_{23} &=&\frac{1}{4}(1-\frac{\Omega_{dd}}{\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}})\rho_{44} - \frac{1}{2}\rho_{22}
\nonumber \\
&+& \frac{1}{4}(1 + \frac{\Omega_{dd}}{\sqrt{\Omega^{2}_{dd}+(4\Omega)^{2}}})\rho_{33}, \nonumber \\
\tilde \rho_{24} &=&\frac{\Omega}{\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}}(\rho_{33}-\rho_{44}) + \frac{\bar c}{\sqrt{2}}\rho_{14}, \nonumber \\
\tilde \rho_{31} &=&(\tilde \rho_{13})^{\dagger}, ~~~ \tilde \rho_{32} = \tilde \rho_{23}, ~~~
\tilde \rho_{33} = \tilde \rho_{22}, ~~~ \tilde \rho_{34} =\tilde \rho_{24}, \nonumber \\
\tilde \rho_{41} &=&(\tilde \rho_{14})^{\dagger}, ~~~ \tilde \rho_{42} =(\tilde \rho_{24})^{\dagger}, ~~~ \tilde \rho_{43}=(\tilde \rho_{34})^{\dagger},
\nonumber
\end{eqnarray}
\begin{eqnarray}
\tilde \rho_{44} &=& \frac{1}{4}(1+\rho_{11}-\rho_{22}) - \frac{\Omega_{dd}}{4\sqrt{\Omega^{2}_{dd} +(4\Omega)^{2}}}(\rho_{33}-\rho_{44})
\nonumber \\
&+&\frac{\bar a}{\sqrt{2}}(\rho_{41}+\rho_{14}). \label{roq12}
\end{eqnarray}
Inserting the matrix (\ref{rab}) in the expression (\ref{mp}) one can obtain the corresponding eigenvalues of $Q$ after some algebraic manipulations.
The corresponding steady-state behaviors for the concurrence $C$ are shown in the Figure~\ref{fig-2ab}(b) which were obtained with the help of
Eqs.~(\ref{eqsm}) as well as Exps.~(\ref{ccr}-\ref{roq12}). One can observe that the entanglement creation among the two-level qubit pair is
accompanied by cooling of the boson mode (compare the solid curves in Fig.~\ref{fig-1ab}(a) and Fig.~\ref{fig-2ab}(b), respectively). Moreover,
the maximal cooling effect is achieved when the entanglement is maximal as well. Note that even higher magnitudes for the concurrence $C$ can be
obtained but in the absence of the qubit's coupling to the single-mode boson field (compare the solid and dashed curves in Fig.~\ref{fig-2ab}b) meaning
that the maximal cooling efficiency realizes on the expense of the entanglement creation. Anyway, at the beginning of the steady-state evolution the
concurrence $C$ is larger than its value which would be obtained, however, in the absence of the boson mode, see Fig.~\ref{fig-2ab}(b). Furthermore,
for lower bath temperatures, the magnitude of the concurrence $C$ will reach the same values regardless of the boson mode presence. This is because
at those temperatures the corresponding two-qubit cooperative states, responsible for entanglement creation, are almost equal populated in both cases.

\section{Summary}
Summarizing, we have investigated the relationship among the entanglement creation in a laser-pumped dipole-dipole interacting two-level qubits
and the cooling effects of a boson mode which is longitudinally coupled with the both quantum emitters, respectively. We have found that cooling
occurs when the dipole-dipole frequency shift lies around the boson mode frequency. This happens because it was assumed that the driving coherent
field frequency is equal with the qubit's transition one, while the two qubits are in an equivalent position with respect to the pumping external field.
Hence, the only way to excite the two-qubit sample is via a concomitant absorption of a photon and a phonon, respectively, leading to cooling of the
boson mode. Furthermore, the quantum cooling process is accompanied by entanglement creation within the qubit sample which is demonstrated by
nonzero values for the concurrence, although its maximal magnitude, i.e. for $C$, is lower than that which would be obtained but in the absence of
the single-mode boson field. However, adjusting the external parameters one can optimize the entanglement as well. These effects are taken place
for rather weak external applied fields which may protect the sample from deteriorations.

This work was supported by grant No. 15.817.02.09F. Also, M.M.A. is grateful for the nice hospitality of the Theory Department of the
Horia Hulubei National Institute of Physics and Nuclear Engineering, Bucharest, Romania.