\section{Introduction}\label{sect_intro}

Strong magnetic fields emerge in various areas of physics, e.g.~in cosmology \cite{1,2}, in non-central heavy-ion collisions \cite{3,4,5,6}, in neutron stars physics \cite{7}; see \cite{8} for review on strongly interacting matter in magnetic fields.

The influence of a magnetic field (m.f.) on QCD thermodynamics was studied in many model approaches \cite{9}.
In particular, for the results within the Nambu--Jona--Lasinio (NJL) model and the holographic approach see \cite{10,11} and \cite{12}, respectively.

The principal problem of a strongly interacting system description is its non-perturbative dynamics.
In the present paper, we address the problem with the Field Correlator Method \cite{13} (see \cite{16} for a recent review).
The QCD thermodynamics within FCM was initially developed in \cite{17} with incorporation of Polyakov loops (for a comparison with lattice data, see \cite{18}).

Recently, the np QCD thermodynamics description at temperatures above deconfinement within FCM was improved by incorporating the color-magnetic confinement (CMC) \cite{20,19,21,22,23}.
CMC produces the effective gluon mass, the Debye screening mass \(m_D\), that linearly grows with temperature above deconfinement.

With this improvement, a solution to the famous Linde problem was suggested in \cite{19}.
A compatible with lattice data SU(3) thermodynamics description was constructed in \cite{20,21,22}.
An improved EoS for \(n_f=2+1\) QCD at nonzero baryon density   was derived in \cite{23}.
In the zero density limit, the results           \cite{20,19,21,22,23} are in good agreement with lattice data \cite{24},  using the Debye mass  found in  \cite{25}.
At nonzero density, the absence of the critical point in \((\mu,~T)\)-plane was demonstrated, which is in agreement with lattice data analysis.
In addition, the sound speed in QGP at zero and nonzero density was determined in \cite{26} and \cite{27}, respectively. For a summary of FCM results at  nonzero  density see \cite{Khaisim}.

The topic of the  np dynamics behind the Polyakov loop was suggested in \cite{28}.

The important subject of  thermodynamics in magnetic field was studied in detail in lattice and analytic calculations in \cite{an1,an2,an3,an4,an5,an6, ant1,ant2,thLatt2} at zero baryonic density.

In the present paper, we continue studying QCD thermodynamics in a uniform  magnetic field within FCM \cite{31,32,33,34} applying our methods to nonzero  density.
We start with the quark contribution to the dense deconfined QCD free energy with np interaction incorporated by means of CMC and Polyakov loops.
Then we replace the quark energy with the relativistic charge energy in a m.f to obtain an expression for the free energy in a m.f.
Finally, we derive a EoS for the finite baryon density QGP in a m.f., or expression for the QGP pressure, from the free energy expression;
we define the QGP magnetic susceptibility with the quark contribution to the QGP pressure.
We check our results in a strong m.f.~at the zero density limit by numerical comparison with lattice data.

The paper is organized as follows.
In section \ref{sect_press} we outline a dense QGP EoS derivation in the absence of a m.f.
In section \ref{sect_magpress} we generalize the EoS to the case of a nonzero m.f.
In section \ref{sect_susc} we calculate the QGP magnetic susceptibility.
In section \ref{sect_num} we present a numerical evaluation of our analytical results.
Section \ref{sect_discuss} contains conclusion and discussion.

\section{Quark pressure with Polyakov loops and CMC}\label{sect_press}

We start with a \(f\)-flavored quark (the index \(f=u,d,s...\) is suppressed unless otherwise stated) free energy in background color and electromagnetic fields $A_\mu (x)$ and $A^{(e)}_\mu (x)$ (for the gluon contribution to the QGP pressure, see \cite{20})

 \begin{eqnarray} \frac{1}{T} F_q (A,A^{(e)} ) = - \frac12 \operatorname{tr } \int^\infty_0 \xi (s)
 \frac{ds}{s} d^4 x\overline{(Dz)}^w_{xx} e^{-K-sm_q^2} \langle W_\sigma
 (C[z])\rangle, \label{1}\end{eqnarray}
where $K=\frac14 \int^s_0 \left( \frac{dz_\mu}{d\tau} \right)^2 d \tau$ is the kinetic kernel (we utilized the Feynman-Fock-Schwinger proper time formalism), $m_q$ is the current quark mass, and

\begin{eqnarray} W_\sigma (C[z]) = P_F P_A \exp \left(ig \int_{C[z]} A_\mu dz_\mu + ie \int_{C[z]}
A_\mu^{(e)} dz_\mu\right) \exp \int^s_0 \left(g\sigma_{\mu\nu} F_{\mu\nu} +
e\sigma_{\mu\nu} F_{\mu\nu}^{(e)}\right)  d\tau\label{2}\end{eqnarray}
is a path- and surface-ordered Wilson loop;  $\langle ... \rangle$ denotes the averaging over the stochastic color background field.
The integration measure implies periodic boundary conditions (\(\beta=T^{-1}\) is the inverse temperature) on a contour $C[z]$ (for details, see \cite{20,23})
\begin{eqnarray}
\overline{(Dz)}^w_{xy} = \prod^n_{m=1} \frac{d^4\Delta
z_k(m)}{(4\pi\varepsilon)^2}\sum_{n=0,\pm 1,\pm 2} (-)^n \frac{d^4p}{(2\pi)^4}
e^{ip_\mu(\sum\Delta z_\mu(m) - (x-y) - n \beta \delta_{\mu 4})}.\label{3}\end{eqnarray}

We write the series over the Matsubara frequencies in the form
  \begin{eqnarray}
P_q = 2 N_c\int^\infty_0  \frac{ds}{s} e^{-m^2_q s} \sum^\infty_{n=1}
(-)^{n+1}[S^{(n)} (s) + S^{(-n)} (s)],\label{4}\end{eqnarray}
where
\begin{eqnarray} S^{(n)} (s) = \int
(\overline{Dz} )^w_{on} e^{-K} \frac{1}{N_c} \operatorname{tr }  W_\sigma (C[z]).\label{5}\end{eqnarray}
The T-independent \(n=0\) term was subtracted.

At this point, we introduce the main np thermodynamics ingredients --- CMC and the Polyakov loop \(L\).
Within the FCM assumptions, color-electric and color-magnetic stochastic background fields are statistically independent.
Hence, the Polyakov loop and the spatial 3d projection of the quark Green's function \(S_3(s)\),
which is subjected to CMC, factorize
  \begin{eqnarray} S^{(n)}(s) =  \int (Dz_4)^w_{on}  e^{-K} L^nS_3(s) = \frac{e^{\frac{-n^2\beta^2}{4s}}L^n}{\sqrt{4\pi s}} S_3(s) .\label{6}\end{eqnarray}

We introduced the quark density by replacing the free energy in \eqref{1} with the thermodynamic potential \(\Omega=F-\mu N\)
(\(\mu\) is the f-flavored quark chemical potential).
The pressure definition is $P = -\left( \frac{\partial \Omega}{\partial V} \right)_{T,\mu}$, so
\begin{eqnarray}  P_q = \frac{ 4N_c}{\sqrt{4\pi}} \int^\infty_0 \frac{ds}{s^{3/2}} e^{- m_q^2s} S_3 (s) \sum_{n=1,2,...} (-)^{n+1} e^{-\frac{n^2}{4T^2s}}\cosh\left(\frac{\mu n}{T}\right) L^n,
 \label{7}\end{eqnarray}

We may extract the quark Polyakov loop $L=\exp \left(-\frac{V_1 (\infty, T)}{2T} \right)$ from lattice simulations (in which case, we may have to change its normalization) or calculate it within FCM \cite{28} (in which case, we still indirectly use lattice data on deconfinement dynamics).

Following \cite{20,23}, we utilize the approximation
\begin{eqnarray} S_3(s) \simeq \frac{1}{(4\pi s)^{3/2}} e^{- \frac{m^2_Ds}{4}}, ~~ m^2_D = c_D^2\sigma_s (T),\label{9}\end{eqnarray}
where \(c_D\sim 2\) \cite{25} is a free parameter in our framework, and the color-magnetic string tension $\sigma_s(T) \sim g^4 T^2$ calculated in \cite{25} defines magnetic Debye screening mass $m_D$.

Therefore, the EoS is
\begin{eqnarray} \frac{P_q}{T^4} = \frac{N_c}{4\pi^2} \sum^\infty_{n=1} \frac{(-)^{n+1}}{n^4} L^n\cosh \left(\frac{\mu n}{T}\right) \Phi_n (T), \label{10}\end{eqnarray}
where $K_2$ is the Macdonald function
\begin{eqnarray} \Phi_n (T) = \frac{8n^2\bar M^2}{T^2} K_2 \left( \frac{\bar M n }{T} \right), ~~ \bar M = \sqrt{m^2_q + \frac{m^2_D}{4}}.\label{11}\end{eqnarray}

To reveal the EoS analytical structure in the complex \(\mu\)-plane, we cast the series to the integral form following \cite{23}
\begin{eqnarray} \frac{P_q}{T^4} = \frac{N_c}{\pi^2} (\xi_+ + \xi_-),\label{13}\end{eqnarray}
\begin{align}
    \xi_\pm &= \sum_{n=1}^\infty \frac{(-)^{n+1}}{n^2} L^n e^{\pm \mu n\beta} \left(\frac{\bar M}{T}\right)^2 K_2 \left( \frac{\bar M n}{T}\right)\label{12}\\
            &= \frac{1}{12} \left( \frac{\bar M}{T} \right)^4 \int^\infty_0 \frac{u^4du}{\sqrt{1+u^2}}\frac{1}{1+\exp \left( \frac{\bar M}{T} \sqrt{1+u^2} + \frac{V_1}{2T} \mp \frac{\mu}{T} \right)}.\label{14}
\end{align}
Cuts in the complex plane immediately follow from \eqref{14}: \(\frac{|\Re\mu|}{T}\ge\frac{\bar M}{T}+\frac{V_1}{2T}\), \(\frac{\Im\mu}{T}=\pi(2n+1)\).

In what follows, we also use another form \cite{33} of EoS \eqref{10}
\begin{eqnarray} P_q = \frac{2 N_c}{\sqrt{\pi}} \int \frac{ d^3p}{(2\pi)^3} \sum^\infty_{n=1} (-)^{n+1} \sqrt{\frac{2}{\beta n}}
    L^n \cosh\left(\frac{\mu n}{T}\right)  \int^\infty_0 \frac{d\omega}{\sqrt{\omega}}
    e^{-\left( \frac{\bar M^2+\mbox{\boldmath${\rm p}$}^2}{2\omega} + \frac{\omega}{2} \right) n\beta}.\label{15}
\end{eqnarray}

\section{Quark pressure in a magnetic field}\label{sect_magpress}

Let us adjust \eqref{15} to the case of a uniform m.f.~directed along the z-axis.

The m.f.~alters quarks motion in the transverse plane.
Transverse motion appears in \eqref{15} through the energy spectra and the phase space.
To make up for the m.f., we introduce the Landau energy levels and modify the phase space following the standard prescription \cite{35}

\begin{eqnarray} \varepsilon\to\varepsilon^\sigma_{n_\bot}(p_z)= \sqrt{ p^2_z+ \bar M^2 + |e_q B | (2n_\bot +1 - \bar \sigma)},~~ \bar \sigma= \frac{e_q}{|e_q|} \sigma .\label{16}\end{eqnarray}
\begin{eqnarray} 2\int\frac{V_3d^3\mbox{\boldmath${\rm p}$}}{(2\pi)^3}\to \sum_{\sigma=\pm1}~\sum_{n_\bot=0}^{\infty}\int\frac{dp_z}{2\pi} \frac{|e_q B|}{2\pi} V_3.\label{18}\end{eqnarray}

The standard pressure definition may require an elaboration in case of a non-trivial phase space.
However, the quark pressure along the m.f.~is well-defined in the present case
\begin{eqnarray} P_z = -\left(\frac{\partial \Omega}{\partial V_z} \right)_{T,\mu},\quad dV_z=S_{\perp}dz.\label{pres} \end{eqnarray}
See Section 5 for the discussion of pressure anisotropy in magnetized QGP.
In what follows, we suppress the subscript, and write \(P\equiv P_z\).

With this reservation, we obtain
\begin{gather}
    P_q(B) = \sum_{n_\bot, \sigma} N_c   T \frac{|e_qB|}{2\pi}
    \frac{\chi(\mu)+\chi(-\mu)}{8},\label{19}\\
    \chi(\mu)=\int \frac{dp_z}{2\pi} \ln
    \left(1+\exp \left( \frac{\mu -V_1/2-\varepsilon^\sigma_{n_\bot}(p_z)}{T}\right)
    \right).  \label{20}
\end{gather}
The integration yields (\(\varepsilon^\sigma_{n_\bot}\equiv\varepsilon^\sigma_{n_\bot}(0)\))
\begin{eqnarray} P_q (B,T) = \frac{N_c  |e_qB|T}{\pi^2} \sum_{n_\bot, \sigma} \sum^\infty_{n=1} \frac{(-)^{n+1}}{n}
L^n \cosh \left( \frac{\mu n}{T} \right)  \varepsilon^\sigma_{n_\bot} K_1 \left( \frac{n
\varepsilon^\sigma_{n_\bot}}{T}\right),\label{21}\end{eqnarray}
and the summation of Landau levels  over $n_\bot$
yields the expression for f-flavored quark pressure at finite density in a m.f.~(along the m.f.)
\begin{align}
    P_q (B) = &\frac{N_c e_q BT}{\pi^2} \sum^\infty_{n=1} \frac{(-)^{n+1}}{n}L^n \cosh\left(\frac{\mu n}{T}\right)
    \left(\bar M K_1 \left( \frac{n\bar M}{T}\right)+\right.\nonumber\\
    &+\left.\frac{2T}{n} \frac{e_qB+\bar M^2}{e_qB} K_2 \left( \frac{n}{T} \sqrt{e_q B + \bar M^2}\right)
    - \frac{ne_qB}{12T} K_0\left( \frac{n}{T} \sqrt{e_q B + \bar M^2}\right)\right).\label{31}
\end{align}

The accuracy of the $n_\bot$  summation in (\ref{21}) exemplified in (\ref{31}) is checked below in section V and Fig.1.

The total quark pressure is the sum over the flavors
\begin{eqnarray} P^{(tot)}_q (B,T) = \sum_f P_q^{(f)}(B,T).\label{27}\end{eqnarray}

As in the previous section, we cast the series to the integral form (see the Appendix of \cite{33} for details)
\begin{eqnarray} \sum_{n_\bot, \sigma}  \chi(\mu) = (I_1+ I_2+I_3),\label{23}\end{eqnarray}
where
\begin{gather}
    I_1 = \frac{1}{\pi T}\int^\infty_0 \frac{ p_z d p_z}{1+\exp \left( \frac{\sqrt{p^2_i+\bar M^2}-\mu+V_1/2}{T} \right)},\label{24}\\
    I_2 = \frac{1}{\pi T}\int^\infty_0 \frac{ p_z^2 d p_z} {|e_q B|} \int^\infty_0 \frac{ d\lambda}{ \sqrt{p^2_z+\bar M^2 +|e_q B|+\lambda}}
    ~\frac{1}{ 1+\exp \left( \frac{\sqrt{p^2_z+\bar M^2 +|e_q B|+\lambda}-\mu+V_1/2}{T}\right)}, \label{25}\\
    I_3 = -\frac{|e_q B|}{24\pi T }  \int^\infty_{-\infty}   \frac{ dp_z}{ \sqrt{p^2_z+\bar M^2 +|e_q B|}}
    ~\frac{1}{ 1+\exp \left( \frac{\sqrt{p^2_z+\bar M^2 +|e_q B|} - \mu+V_1/2}{T}\right)}. \label{26}
\end{gather}
The integral form of \eqref{31} is
\begin{eqnarray} P_q (B,T) =N_c T \frac{|e_qB|}{2\pi^2}  \sum^3_{i=1} (I_i(\mu) + I_i(-\mu)).\label{28}\end{eqnarray}

Finally, we probe the result \eqref{28} with the \(B\to 0\) limit.
$I_1(\mu)$ does not depend on $B$, $I_3(\mu, B)$ is $O(|B|)$, and $I_2(\mu, B)$ is $O(|B|^{-1})$, so \(B\) drops out
\begin{eqnarray} P_q = \frac{N_c T^4}{12\pi^2} \left[ \varphi\left( \frac{\mu - \frac{V_1}{2}}{T}, \frac{\bar M}{T}\right)
    + \varphi\left( -\frac{\mu +\frac{V_1}{2}}{T}, \frac{\bar M}{T} \right)\right]\label{29}\end{eqnarray}
\begin{eqnarray} \varphi(a, \nu) = \int^\infty_0 \frac{z^4 dz}{\sqrt{z^2+\nu^2}} \frac{1}{\exp (\sqrt{z^2+\nu^2} -a)+1}.\label{30}\end{eqnarray}
As expected, the expression coincides with the corresponding result \eqref{13} of section \ref{sect_press}.

\section{QGP magnetic susceptibility at nonzero baryon density}\label{sect_susc}

We use the following definition of magnetic susceptibility \cite{an3} of a thermal medium
\begin{eqnarray}  P(B,T) - P(0,T) = \frac{\hat\chi}{2}(eB)^2 + O((eB)^4). \label{34}\end{eqnarray}

The gluon part of the QGP free energy is sensitive to a m.f.~via the quark loops.
The corresponding contribution to \(\hat\chi\) is suppressed at least as \(\alpha_s^2\).
In the present paper, we neglect this contribution
\begin{eqnarray} \hat\chi\simeq\hat\chi_q^{(tot)}.\end{eqnarray}
In a sense, this approximation is similar to the staggered fermions approach.

We also introduce the f-flavored quark susceptibilities \(\hat\chi_q^{(f)}\)
\begin{eqnarray}  P_q^{(f)}(B,T) - P_q^{(f)}(0,T) = \frac{\hat \chi_q^{(f)}}{2}(e_q^{(f)}B)^2 + O((eB)^4), \label{34*}\end{eqnarray}
which are related to the total quark susceptibility as
\begin{equation}
    \hat\chi_q^{(tot)}= \sum_f\hat\chi_q^{(f)}(e_q^{(f)}/e)^2.\label{34**}

\end{equation}
We extract the susceptibilities from \eqref{31} by expanding it in powers of \(e_qB\)
\begin{eqnarray}  \hat \chi_q = \frac{   N_c}{3\pi^2} \sum^\infty_{n=1}  (-)^{n+1}
    L^n\cosh \left(\frac{\mu n}{T}\right) K_0 \left( \frac{n \bar M}{T} \right). \label{35}\end{eqnarray}

With the integral representation for the Macdonald function
\begin{eqnarray} K_0 (z) = \frac12 \int^\infty_0 \frac{dx}{x} e^{-\left( \frac{1}{x} + \frac{z^2 x}{4} \right)},\label{36}\end{eqnarray}
the integral form of the susceptibilities was obtained
\begin{eqnarray}  \hat \chi_q =\frac{   N_c}{3\pi^2}  (I_q(\mu) + I_q(-\mu)), \label{37}\end{eqnarray}
\begin{eqnarray} I_q (\mu) = \frac12 \int^\infty_0 \frac{dx }{x} \frac{L e^{(\mu/T)} e^{-\left( \frac{1}{x} + \frac{\bar M^2  x}{4T^2} \right)}}{1+ L e^{({\mu}/{T})} e^{-\left( \frac{1}{x} + \frac{\bar M^2  x}{4T^2} \right)}}\label{38}\end{eqnarray}

\section{Numerical analysis of the results}\label{sect_num}

We computed thermodynamic quantities of QGP with u-, d-, s-quarks by means of corresponding series partial summation.
The np inputs $V_1(\infty,~T)$ and \(m_D(T)\) are identical to these of \cite{Khaisim}.
They were adjusted in such a way that QGP pressure \eqref{10} at zero m.f.~and baryon density matches the lattice data \cite{qlatt2,latt3}.
The Debye mass is defined by the spatial string tension value \(m_D= c_D \sqrt{\sigma_s(T)}\); for \(\sigma_s(T)\) see (45) in \cite{25}.
The parameter \(c_D=1.6\) is close to the corresponding result of \cite{25}.

The major source of systematic error is $V_1(\infty,~T)$ uncertainty.
Within the parameters range, we estimate the total error of this section numerical results to be at most \(15\%\) (see section \ref{sect_discuss} for details).

Let us start with an accuracy check of the approximate Landau levels summation procedure.
The Fig.~\ref{fig_summ} demonstrates the difference between the result of \eqref{21} based on the partial sum  and   the integral form in \eqref{31}.
We find the difference to be negligible in the considered parameter range.

\begin{figure
    \center{\includegraphics[width=0.7\linewidth]{Comparison}}
    \caption{The  pressure $P/T^4$ from  \eqref{21}, the dash-dotted line,  compared with the integral  form, Eq. (\ref{31}) the dashed line.}
    \label{fig_summ}
\end{figure}

The set of trajectories for the pressure magnetic shift $\Delta P = P(eB,~\mu_B) - P(0,~\mu_B)$, computed with \eqref{31}, is presented in Fig.~\ref{fig1}
(baryon and quark chemical potentials are related as $\mu_f = \frac{1}{3}\mu_B$, $f = u,d,s$).
The results suggest strong paramagnetism of QGP, in agreement with the lattice study \cite{an2}.

\begin{figure}
    \center{\includegraphics[width=0.7\linewidth]{dataMFMU.pdf}}
    \caption{The magnetic shift $\Delta P = P(eB,\mu) - P(eB = 0,\mu)$ for dense QGP with u,~d,~s quark flavors.
            For a given m.f., each trajectory demonstrates triple splitting for different $\mu_B$.
            the Solid lines are for $\mu_B = 0$, the dash-dotted --- for $\mu_B = 0.2~GeV$, and the dashed --- for $\mu_B = 0.4~GeV$.}
    \label{fig1}
\end{figure}

With the expressions for gluon pressure from \cite{31} and for quark pressure \eqref{31}, we computed a special combination \(\Delta\) for dense QGP in an external m.f.
\begin{eqnarray}
    \Delta = \frac{\epsilon-3P_z+\mu n}{T^4}.
\end{eqnarray}
If pressure is considered as isotropic,   \(\Delta\) is the trace anomaly.

The resulting curves in Fig.~\ref{fig_anomaly} demonstrate all the qualitative features (the maximum position and width dependence on m.f., the curves intersections) of the same quantity computed with the lattice simulation \cite{an6} (see Fig.~12 of \cite{an6}) in the given temperature range.

\begin{figure}
    \includegraphics[width=0.7\linewidth]{Anomaly.pdf}
    \caption{The combination \(\Delta\) for QGP with u,~d,~s quark flavors as a function of \(T\) at various values of \(\mu_B[GeV]\) and \(eB[GeV^2]\).}
    \label{fig_anomaly}
\end{figure}

\section{Conclusion and Discussion}\label{sect_discuss}
In the present paper, pressure of QGP at nonzero baryon density in an external uniform m.f.~was calculated analytically for the first time

Within the parameters range (\(eB<0.5~GeV^2,~\mu_B<0.5~GeV\)), we neglect the vacuum properties alternation.
That is, we neglected the Polyakov loop and the Debye mass dependencies on the baryon density and the m.f.
According to the recent lattice data, the approximations \(L(T,\mu_B,B)\to L(T),~m_D(T,B)\to m_D(T),~m_D(T,\mu)\to m_D(T)\) yield up to 7\%
\cite{09547}, 5\% \cite{00842}, and 2\% \cite{09461} errors, respectively.
We consider the errors as independent, and estimate  the total EoS \eqref{31} error in ``the worst'' parameters region (low temperature, high m.f.~and density) as 15\%.
At larger densities, the Polyakov loop dependence on the chemical potential is to be accounted for.

The Polyakov loops interaction and the perturbative corrections were also neglected.
The non-perturbative numerical input data were set as to reproduce the correct QGP pressure values at the zero m.f.~and baryon density limit.

With the numerical analysis of our results, we observed QGP strong paramagnetism, which was previously found within lattice simulations.
The possible pressure anisotropy, on the other hand, is a matter of discussion.

QGP pressure anisotropy is observed in lattice simulations with a certain setup.
Basic methods for pressure calculation on lattice are \cite{an1}
\begin{enumerate}
\item $B=const$ method --- the variation of energy in a constant external m.f.~(the QGP is assumed to be in thermal equilibrium);
\item $\Phi=const$ method, in which the total magnetic flux $\Phi$ through the lattice is fixed.
\end{enumerate}
With the first method, pressure is isotropic, while with the second method
\begin{equation}
P_x = P_y = P_z - M \cdot B, \label{39}
\end{equation}
where $M = -\frac{1}{V}\frac{\partial F}{\partial B}$ is the QGP magnetization.

The ``anisotropic'' method is more subtle for real physical systems.
It requires a careful analysis of the energy exchange between longitudinal and transverse motion of quarks and the corresponding relaxation time.
We leave these details to a future publication.

The direct influence of a m.f.~on the field correlators was studied on the lattice in  \cite{07012}.
Since the effect is at most 10\% at \(eB<0.5~GeV^2\), we had assumed it to be negligible for our purposes.
The more exact accuracy analysis with the Polyakov loop and the screening mass variation justified this assumption.

\section{Acknowledgments}
The authors are grateful to B.O.~Kerbikov and M.S.~Lukashov for useful discussions, and to N.P.~Igumnova for technical help in the manuscript preparation.
This work was supported by the Russian Science Foundation Grant number 16-12-10414.


\section{\large Introduction}
\label{intro}

The discovery by the ATLAS and CMS collaborations at the CERN Large Hadron Collider (LHC)~\cite{ATLAS, CMS} of a bosonic particle, with properties which are compatible with those predicted for the Higgs boson of the Standard Model (SM), represents a significant progress in our understanding of the electroweak symmetry breaking mechanism. The SM is theoretically consistent with the inclusion of a $125~\rm{GeV}$ Higgs boson, in the sense that no Landau pole emerges, also if the model is extrapolated up to the Planck scale ($\Lambda_P\approx 10^{18}~\rm{GeV}$), where one has to accept a meta-stable vacuum and an unnatural high amount of fine-tuning ($10^{34}$) for the prediction of the Higgs boson mass at the electroweak scale ($\Lambda_{EW}\approx 10^2~\rm{GeV}$)~\cite{Isidori, Buttazzo, Bezrukov, Kniehl, kniehl2}. However, there are still several puzzles that remain unsolved by the SM dynamics. The hierarchy problem, the neutrino oscillation, the identification of the dark matter, the baryon asymmetry, among others, are all left unanswered and require new physics beyond the Standard Model. The minimal super-symmetric extension of the SM (MSSM) is the best motivated and the most intensively studied framework of new physics, providing a widely amount of precise predictions for experimental phenomena at the TeV scale~\cite{Nilles,HabernKane}. In most scenarios that are phenomenologically relevant~\cite{Stal, Carena, Bagnaschietal,HST} the LHC measured value, $M_h^{exp}=125.09\pm 0.24~\rm{GeV}$~\cite{MhMass,MhCMS,MhATLAS}, is associated with the lightest CP-even Higgs boson mass ($M_h$) which is theoretically predicted with great accuracy in the MSSM. Up to now, the dominant quantum corrections to $M_h$ have been computed at one-loop~\cite{Chankowski,Dabelstein,Pierce,Frank}, two-loop~\cite{alphaalphas,alphabalphas,alphatalphas,Carena2,CarenaHollik,MartinHiggs,Sophia,
DegrassiMSSM,SophiaPaser} and three-loop \cite{Harlander1,Harlander2,Kant,Harlander3,Edilson} level using the Feynman diagrammatic (FD) and the effective potential (EP) approaches. These MSSM predictions can accommodate the measured Higgs mass value of $125~\rm{GeV}$ and are consistent with the similarities of the measured Higgs couplings to those in the SM~\cite{SMCouplings}. Effective field theory (EFT) methods have been also considered to resum large-logarithms in case of a large mass hierarchy between $\Lambda_{EW}$ and the SUSY scale ($M_{SUSY}$)~\cite{Draper,Lee,Villadoro,Bagnaschi}. In particular, for values of $M_{SUSY}$ above a critical point where the fixed-order and EFT combined uncertainties are equal, the EFT computation is more accurate and therefore the usage of the SM~\cite{Allanach} or a two-Higgs-doublet-model (THDM)~\cite{Bahl2} as effective theories below the SUSY scale is preferred. Both the fixed-order and the EFT results are implemented in several publicly available codes. For the diagrammatic fixed-order calculations there are the programs SoftSUSY~\cite{Allanach2}, SUSPECT~\cite{SUSPECT}, CPSuperH~\cite{CPSuperH} and H3m~\cite{Harlander3}. Pure EFT calculations are implemented in SUSYHD~\cite{Villadoro} and MhEFT~\cite{MhEFT}. Moreover, different hybrid methods that combine both approaches have been recently developed in order to take profit of the features of each one. FlexibleSUSY~\cite{FlexibleSUSY}, based on SARAH~\cite{SARAH}, implements a hybrid method called Flexible-EFT-Higgs~\cite{Athron}. This approach was also included into the program SPheno~\cite{SPheno1,SPheno2}. \\ A hybrid method different from the one pursued in Flexible-EFT-Higgs has been implemented in FeynHiggs~\cite{FeynHiggs1,FeynHiggs2}. There are also in literature detailed numerical comparisons between the different diagrammatic, EFT and hybrid codes. In~\cite{Athron} it is discussed in details how the hybrid method Flexible-EFT-Higgs compares to the other EFT and diagrammatic codes. Finally, several numerical comparisons of the hybrid approach implemented in FeynHiggs to the pure EFT calculations have been studied in~\cite{Villadoro,Athron,SPheno2}. Those papers reported surprising non-negligible numerical differences between FeynHiggs and pure EFT codes for the prediction of $M_h$ at large SUSY scales. The observed differences come mainly from three sources. The scheme conversion of input parameters from OS to $\overline{DR}$, which can lead to large shifts due to uncontrolled higher-order terms. Unwanted effects from incomplete cancellations with subloop renormalization contributions in the determinations of the Higgs propagator pole and different
parametrizations of non-logarithmic terms. After performing the corresponding corrections, FeynHiggs results are in very good agreement with the results of SUSYHD~\cite{Bahl3}. \\ For the present study we  decided to use the fixed-order and EFT hybrid calculations currently included in FeynHiggs, which seems to be in a very good agreement with the other fixed-order and EFT codes and gives a reliable three-loop predictions of the Higgs boson mass for large SUSY scales, in order to provide a numerical comparison of our three-loop fixed-order predictions of the ligthest MSSM Higgs boson mass reported in~\cite{Edilson} with the fixed-order and EFT hybrid results found in literature. As the effects of the large logarithms are expected to become relevant when $M_{SUSY}$ grows, it is natural to ask how large $M_{SUSY}$ can be. We therefore provide in this article a phenomenological analysis about the compatibility of the experimental observations at the LHC for the Higgs boson mass and the region of parameters in the specific MSSM considered scenario to find an upper bound on the needed $M_{SUSY}$.

\section{ \large Three-Loop Fixed-Order Calculation of~${M_h}$}
\label{sec-1}

In contrast to the SM, the Higgs sector of the MSSM with real parameters (rMSSM) contains two complex doublets with opposite hyper-charges
\begin{eqnarray}
\tiny
H_{1}=\left(\begin{array}{c}
H_{1}^{0}+\frac{v_{1}}{\sqrt{2}}\\
H_{1}^{-}
\end{array}\right) & \tiny \rm{and} & \tiny H_{2}=\left(\begin{array}{c}
H_{2}^{+}\\
H_{2}^{0}+\frac{v_{2}}{\sqrt{2}}
\end{array}\right), \label{eq:Higgsdoublets}
\end{eqnarray}
where the neutral components, $H_{1,2}^0$ fluctuate around the vacuum expectation values (vevs) $v_{1,2}$. In the physical basis there are five Higgs bosons, three of them are neutral: the lightest ($h$) and heavy ($H$) CP-even Higgs bosons and the CP-odd Higgs boson ($A$). The other two, $H^{\pm}$, are charged and vev-less. Besides the SM electroweak boson masses, the rMSSM Higgs sector is parametrized in terms of two additional parameters: the mass of the CP-odd Higgs boson ($M_{A}$) and $tan\beta$, which is the ratio of the two vevs, $v_1/v_2$. The masses of the CP-even Higgs boson particles, $h$ and $H$, follow as predictions. \\ We focus in this section on the prediction of the lightest Higgs boson mass, $M_h$, at three-loop accuracy using a fixed-order FD computation which is based on the calculation of Higgs self-energy corrections at the given perturbative order. In this approach, the renormalized CP-even Higgs boson masses are obtained by finding the zeros of the determinant of the inverse propagator matrix (poles equation)
\begin{eqnarray}
\tiny
\left(\Delta_{H}\right)^{-1}=-i\left(\begin{array}{cc}
p^{2}-m_{H}^2 + {\displaystyle\sum_{l=1}^{3}} \widehat{\prod}_{HH}^{(l)} & {\displaystyle\sum_{l=1}^{3}}\widehat{\prod}_{hH}^{(l)}\\
{\displaystyle\sum_{l=1}^{3}}\widehat{\prod}_{hH}^{(l)} & p^{2}-m_{h}^2 + {\displaystyle\sum_{l=1}^{3}}\widehat{\prod}_{hh}^{(l)}
\end{array}\right), \label{eq:inversematrix}
\end{eqnarray}
where $m_{h}$ and $m_{H}$ denote the tree-level mass of $h$ and $H$ respectively and
\begin{eqnarray}
\tiny
\left.\widehat{\prod}\right._{ij}^{(l)}=\left.\prod\right._{ij}^{(l)}~-~\delta^{(l)}M_{ij}^{2} & \tiny ; & \tiny i,\: j = h,\: H, \label{eq:SEandRC}
\end{eqnarray}
are the corresponding $l$-loop renormalized self-energies. A particular feature of the rMSSM is the large size of the higher order quantum corrections to masses and couplings. They can lead to a considerably large shift on the value of the Higgs boson mass, where the bulk of the corrections comes from the SUSY-QCD sector of the Lagrangian. Thus, the dominant contributions to $\left.\widehat{\prod}\right._{ij}$ in eq. (\ref{eq:SEandRC}) involve the SM parameters $h_t$ (top Yukawa coupling), $M_t$ (top quark mass), $\alpha_s$ (strong coupling constant) and the MSSM parameters $M_{\tilde{g}}$ (gluino mass), $\theta_t$ (stop mixing angle), $\tilde{m}_{q_{1,2}}$ (squark masses) and $A_q$ (soft breaking parameters) where $q=u,d,t,b,c,s$. \\ Concerning the renormalization of the self-energy corrections, that is to say, the determination of the mass counter-terms $\delta^{(l)}M_{ij}^{2}$, we follow the mixed OS/$\overline{\rm{DR}}$ scheme defined in~\cite{Edilson}. Thus, the electroweak gaugeless limit at O($\alpha_t\alpha_s^2$) and the approximation of zero external momentum are assumed. As a consequence, we have avoided dealing with the Higgs wave function renormalization and also with the renormalization of $tan\beta$. Moreover, $v_{1,2}$ are defined as the minima of the full effective potential and therefore the tadpoles are renormalized on-shell according to the conditions:
\begin{eqnarray}
\tiny
T_{1,2}^{tree}=0, & \tiny \delta^{(l)}T_{1,2}=-T_{1,2}^{(l)},\label{eq:T(OS)}
\end{eqnarray}
where $T_{j}^{(l)}$ is the $l$-loop Higgs tadpole contribution. We have also imposed an on-shell renormalization to the $A$-boson mass,
\begin{eqnarray}
\tiny
\delta^{(l)}M_{AA}^2=Re\left[\left.\prod\right._{AA}^{(l)}\left(M_{A}^{2}\right)\right]. \label{eq:MA(OS)}
\end{eqnarray}
The three-loop corrections of O($\alpha_t\alpha_s^2$) also include the O($\alpha_s$) contributions to the one-loop counter-terms coming from the renormalization of the gluino mass, the top quark mass, the squark masses and the stop mixing angles in the $\overline{DR}$ scheme, as well as the two-loop $\overline{DR}$ renormalization of the top mass, the stop masses and stop mixing angles at O($\alpha_s^2$).\\ For the purposes of this article we have chosen a degenerate single-scale scenario  where all the super-symmetric masses are set equal to an effective scale $M_{SUSY}$,
\begin{eqnarray}
\tiny
M_{L,R} = M_{\tilde{g}} = M_A =\mu = M_{SUSY}. \label{eq:scales}
\end{eqnarray}
Here $\mu$ is the Higgsino mass and $M_{L,R}$ are the soft SUSY-breaking masses. We have also identified the lightest Higgs boson $h$ as the SM-like Higgs boson and therefore we have assumed the decoupling limit, $M_A = M_{SUSY} \gg M_t$. This degenerate scenario in the decoupling limit is known as the ``heavy SUSY" limit. As a consequence, the three-loop self-energy corrections to $m_{h,H}^{2}$ can approximately be obtained as a superposition of the $33$ vacuum integrals depicted in Fig. \ref{MasterIMS} with coefficients that are functions of the kinematic invariants and the space-time dimension.
\begin{figure}[h]
\centering
\includegraphics[width=18pc]{MasterIMS.pdf}
\caption{\small{Basis of three-loop Master Integrals. The dashed line represents a massless propagator. The thin solid line is the propagator with a mass at the electroweak scale $M_t$ and the thick solid line depicts the propagator involving the SUSY scale $M_{SUSY}$.}}
\label{MasterIMS}
\end{figure}
Each diagram of the basis in Fig. \ref{MasterIMS} represents a three-loop Master Integral of the form
\begin{eqnarray}
\tiny
I_{v_{1}\ldots v_{6}}=i\frac{e^{3\gamma_{E}\varepsilon}}{\pi^{3D/2}}\int{\displaystyle \prod_{l=1}^{3}}d^{D}q_{l}\left[{\displaystyle \prod_{j=1}^{6}}\frac{1}{P_{j}^{n_{j}}}\right], \label{eq:MI}
\end{eqnarray}
where
\begin{eqnarray*}
\tiny P_{1}=q_{1}^{2}-m_{1}^{2},~~~~~~~ & \tiny P_{2}=\left(q_{1}-q_{2}\right)^{2}-m_{2}^{2}, \\
\tiny P_{3}=\left(q_{2}-q_{3}\right)^{2}-m_{3}^{2}, & \tiny P_{4}=q_{3}^{2}-m_{4}^{2}, \\
\tiny P_{5}=q_{2}^{2}-m_{5}^{2},~~~~~~~ & \tiny P_{6}=\left(q_{1}-q_{3}\right)^{2}-m_{6}^{2}.
\end{eqnarray*}
There are two scales involved, the electroweak scale $M_t$, whose associated propagator is represented with a thin solid line and the super-symmetric scale $M_{SUSY}$ represented with a thick solid line. Massless propagators are represented with a dashed line. This basis was obtained using the integration by parts (IBP) method implemented in the code Reduze~\cite{Reduze}. Main part of the diagrams shown on Fig.~\ref{MasterIMS} have been analytically evaluated in \cite{MIs1,MIs2,MIs3,MIs4,MIs5,MIs6,MIs7,MIs8}. The numerical evaluation of the basis integrals was done with the programs TVID~\cite{Freitas1, Freitas2} and SecDec~\cite{SecDec}. In particular, the integral $I_{211100}$ requires a Laurent expansion up to first order in $\varepsilon$. The evanescent terms of O($\varepsilon^1$) was numerically evaluated with the help of SecDec.

\section{\large EFT Hybrid Calculation of~$\rm{M_h}$}
\label{sec-2}
When there is a large mass hierarchy between the electroweak scale and the scale of the SUSY particles, the fixed-order computations of the Higgs self-energy corrections contain large logarithms that can spoil the convergence of the perturbative expansion and yield unreliable predictions of the Higgs boson masses. A fixed-order computation is thus recommended for low values of $M_{SUSY}$ not separated too much from $M_t$. There is an alternative approach to calculate $M_h$ which yield accurate results for high SUSY scales. This approach is based on the EFT techniques~\cite{Villadoro,Heidi} and allows the resummation of the large logarithmic terms and the incorporation of higher-order contributions beyond the order of the fixed-order diagrammatic calculations. In the heavy SUSY limit the low-scale EFT below $M_{SUSY}$ is the SM. It requires just one EFT coupling, the effective Higgs self coupling $\lambda$, which correlates the high scale $M_{SUSY}$ and the low scale $M_t$ through the renormalization group equations (RGEs) and captures radiative corrections of the form
\begin{eqnarray}
\tiny
\alpha_j^{n+m-1}log^{n}\left(M_{SUSY}/M_t\right) & \tiny ; & \tiny j = \lambda,h_t,g_s,..., \label{eq:largelog}
\end{eqnarray}
for any $n$, by using the $m$-loop beta functions of $\alpha_j$, into the running coupling $\lambda(Q)$. In order to get a SM running Higgs mass in the $\overline{MS}$ scheme at the scale $M_t$, one has to multiply $\lambda(M_t)$ by $2v^2(M_t)$, where $v(M_t)\approx 246~\rm{GeV}$ is the $\overline{MS}$ vev evaluated at $M_t$. The physical Higgs mass requires to solve the pole equation
\begin{eqnarray}
\tiny
p^2 - 2\lambda(M_t)v^2(M_t) + \widetilde{\prod}_{hh}^{SM}\left(p^2\right) = 0, \label{eq:SMpoleEq}
\end{eqnarray}
with the SM Higgs boson self-energy,
\begin{eqnarray}
\tiny
\widetilde{\prod}_{hh}^{SM}\left(p^2\right) = \left[\left.\prod\right._{hh}^{SM}\left(p^2\right)-\frac{1}{\sqrt{2}v}T_h^{SM}\right]_{fin}, \label{eq:SMSE}
\end{eqnarray}
renormalized in the $\overline{MS}$ scheme but with the Higgs tadpoles renormalized to zero, i.e. $\delta T_h^{SM} = -T_h^{SM}$. \\ As higher dimensional operators are not included into the effective Lagrangian, the contributions suppressed by the heavy scale $M_{SUSY}$ are not considered. Consequently, the EFT calculation is less accurate than the fixed-order one for low SUSY scales. The fixed-order calculation is more accurate below a critical SUSY mass scale, estimated to be about $M_{SUSY}^{C}\approx 1.2~\rm{TeV}$ in~\cite{Allanach}, whereas above that scale the EFT calculation is more accurate. \\ In the latest released version of FeynHiggs~\cite{FeynHiggs3} both approaches are combined in order to supplement the full one-loop, leading and sub-leading two-loop diagrammatic results with a resummation of the leading + next to leading (LL+NLL)~\cite{THahn} and next to next to leading (NNLL)~\cite{Bahl1} logarithmic contributions coming from the top/stop sector.
For the resummation of large logarithms up to NLL two-loop RGEs and one-loop
matching conditions are needed, accordingly, the resummation up to NNLL requires three-loop RGEs and two-loop matching conditions. The hybrid results obtained from the combination of the two approaches are added into the pole equation of the full MSSM
\begin{eqnarray}
\tiny
p^2 - m_h^2 + \widetilde{\prod}_{hh}\left(p^2\right) + \Delta_{hh}^{log} = 0, \label{eq:MSSMpoleEq}
\end{eqnarray}
through the shift $\Delta_{hh}^{log}$ which contains the resummed large logarithms from the EFT as well as the logarithmic terms already present in the fixed-order Higgs self-energies,
\begin{eqnarray}
\tiny
\Delta_{hh}^{log} = - \left[2\lambda(M_t)v^2(M_t)\right]_{log} - \left[ \widetilde{\prod}_{hh}\left(m_h^2\right) \right]_{log}. \label{eq:hhShift}
\end{eqnarray}
The subscript $"log"$ means that only logarithmic terms are considered. The logarithms in the Higgs self-energy appear explicitly only after expanding in $v/M_{SUSY}$. This subtraction term ensures that the one- and two-loop logarithms, already contained in the fixed-order FD computation, are not counted twice. In general the higher-order logarithms obtained from the EFT and the hybrid approaches are not the same because the determination of the poles of the propagators (eq. \ref{eq:SMpoleEq} and eq. \ref{eq:MSSMpoleEq}) are performed in different models. However, this difference, which comes from the momentum dependence of the two-loop order non-SM contributions to the Higgs self-energy, cancels out with contributions coming from the subloop renormalization in the heavy SUSY limit, as was explicitly shown in~\cite{Bahl3}. Besides the unwanted effects from incomplete cancellations in the determination of the Higgs propagator pole, the effects due to non-logarithmic terms and its parametrization as well as the higher-order terms coming from the scheme conversion between $OS$ and $\overline{DR}$ parameters are all included into FeynHiggs~2.14~\cite{FeynHiggs3}.

\section{\large Numerical Results}
\label{sec-3}

In this section we present a numerical comparison of our three-loop fixed-order predictions of $M_h$ to the numerical predictions coming from the new version of FeynHiggs.
We have chosen a $\overline{DR}$ renormalization of the stop sector with the renormalization scale set to $\mu_r=M_{SUSY}$, which is equivalent to set $Q_t = -1$ in FeynHiggs. The one-/two-loop fixed-order and the EFT-hybrid FeynHiggs predictions are fixed such that the full MSSM is considered (mssmpart=4) in its real version (higgsmix=2, tlCplxApprox=0), no approximation is taken for the one-loop result (p2approx=4) and the $O(tan^n \beta)$ corrections are resummed (botResum=1). In particular, when the resummation of the large logarithms is included, we use the full LL, NLL and NNLL resummation (looplevel=2, loglevel=3). The top quark mass is renormalized in the SM $\overline{MS}$ scheme at NNLO (runningMT=1) since for loglevel different from zero a $\overline{DR}$ renormalization is not allowed. The input flags of FeynHiggs 2.14.3 are explicitly indicated, for more details the online manual of the code can be consulted at~\cite{manual}. To obtain the pole mass $M_h$ at three-loop level in the fixed-order approach, we have introduced the $O(\alpha_t\alpha_s^2)$ corrections as constant shifts in the FeynHiggs 1-loop + 2-loop Higgs renormalized self-energies~(looplevel=2 and loglevel=0) with the help of the function FHAddSelf but in this case we have used a $\overline{DR}$ renormalization of the top quark mass (runningMT=3).
\begin{figure}[H]
\centering
\includegraphics[width=18pc]{Plota.pdf}  \\
\includegraphics[width=17.9pc]{Plotb.pdf}
\caption{\small{Comparison of the $M_h$ predictions of FeynHiggs with the three-loop fixed-order computation of $M_h$ at $O(\alpha_t\alpha_s^2)$ in the heavy SUSY limit. The dot-dashed and the dashed lines are the fixed-order results of FeynHiggs at one and two -loop level respectively. The blue dotted line contains the NNLL resummation of the large logarithms in FeynHiggs. The blue band corresponds to the uncertainty in the NNLL prediction taken from FeynHiggs. The brown band is the CMS/ATLAS Higgs boson mass, $M_h^{exp}=125.09 \pm 0.24~\rm{GeV}$. The red solid line represents our three-loop fixed-order predictions. Up: Dependence of $M_h$ on the super-symmetric scale $M_{SUSY}$ for a vanishing stop mixing, $X_t/M_{SUSY}=0$. Down: Numerical differences between the FeynHiggs predictions and the three-loop fixed-order predictions of~$M_h$.}}
\label{Fig1}
\end{figure}
We start by considering the FeynHiggs fixed-order, FeynHiggs NNLL hybrid and three-loop $O(\alpha_t\alpha_s^2)$ predictions. The upper plot of Fig. \ref{Fig1} shows the dependence of $M_h$ on $M_{SUSY}$ for a vanishing stop mixing, $X_t/M_{SUSY}=0$, at the kinematic point $A_{e,\mu,\tau,u,d,c,s,b}=0$ and $tan\beta=10$, whereas the lower plot shows the numerical differences between all the considered FeynHiggs results and the $O(\alpha_t\alpha_s^2)$ prediction of $M_h$. In order to draw these plots we have adopted the heavy SUSY limit (eq. \ref{eq:scales}) and we have followed the next conventions. The one and two-loop fixed-order results of FeynHiggs are represented with the dot-dashed and the dashed lines respectively. The blue dotted line contains, in addition, the resummation of the large logarithms up to NNLL order. The blue band corresponds to the uncertainty in the NNLL prediction computed with the help of the FeynHiggs function FHUncertainties for the flag choise: mssmpart = 4, looplevel = 2, loglevel = 3, runningMT = 1.  In principle three effects are taken into account: i) the variation of the renormalization scale from $M_t/2$ to $2M_t$, ii) the use of $M_t^{pole}$ instead of $M_t^{run}$ in the two-loop corrections and iii) the exclusion of higher order resummation effects in $M_b$. The brown band is the experimental Higgs boson mass and its corresponding uncertainty, we have adopted the combined CMS/ATLAS result of the RUN 1 at the LHC, $M_h^{exp}=125.09 \pm 0.24~\rm{GeV}$~\cite{MhMass}, since there is not yet an official combined result for RUN~2~\cite{MhCMS,MhATLAS} observations. Finally, the red solid line contains our three-loop fixed-order corrections. \\
The first thing to note here (and also in Fig. \ref{Fig2}) is that the higher-order large logarithms coming from the EFT hybrid approach at NNLL level produce a growing positive shift on the two-loop predictions reaching a size of about 20 GeV for $M_{SUSY}=40~\rm{TeV}$. Additionally, the NNLL predictions are in a very good agreement with the three-loop $O(\alpha_t\alpha_s^2)$ results for $M_{SUSY}$ less than the value $M_{SUSY} \lesssim 10~{\rm TeV}$.
On the lower graph of Fig.~\ref{Fig1} one can see that in the region $2.2~{\rm TeV} \lesssim M_{SUSY} \lesssim 7.4~{\rm TeV}$ there is an approximately constant difference of about $0.2~\rm{GeV}$ between the red solid and the blue dotted line which is within the theoretical uncertainty (blue band) estimated to be about 0.6~GeV. Below this region the agreement is still good with a numerical difference of at most $1~\rm{GeV}$. However, for scales above 10~TeV the effects of the large logarithms in the red curve start to be relevant, the difference between the two results rapidly increases up to $\sim 21~\rm{GeV}$ when $M_{SUSY}$ grows to up to $20~\rm{TeV}$ and grows monotonically reaching 78 GeV at $M_{SUSY}=40~\rm{TeV}$. This pronounced behaviour depends crucially on our election of the input parameters $\mu_r,~M_{\tilde{g}}~{\rm and}~X_t$. The presence of $n$-loop logarithms of the form $log^{n}\left(M_{SUSY}/M_t \right)$ in the master integrals of Fig.~\ref{MasterIMS} can introduce additional large contributions in the three-loop predictions of $M_h$. \\
In Fig.~\ref{Fig2} the heavy SUSY limit has been smoothed to include an additional SUSY scale, the gluino mass $M_{\tilde{g}}$, in the fixed-order results. The NNLL resummation procedure is restricted to the case of $M_{\tilde{g}}$ equal to $M_{SUSY}$ since three-loop RGEs for an appropriate extension of the Standard Model with the gluino as additional fermion \cite { Bahl1}, for instance as a singlet of the gauge group, are not included in FeynHiggs. However, the main contributions sensitive to the gluino mass are captured by the two-loop result, the numerical effects due to a gluino threshold is numerically small and can be safely neglected, as was shown in~\cite{Bahl1}. We have considered a gluino mass of $M_{\tilde{g}}=1.5~\rm{TeV}$. The inclusion of this additional scale produces sizeable differences between the $O(\alpha_t\alpha_s^2)$ and the  NNLL results. For small SUSY scales below $\sim 3.5~\rm{TeV}$ the difference is always less than $1.3~\rm{GeV}$. For large SUSY scales ($M_{SUSY}>3.5~\rm{TeV}$) this difference grows to a maximum value of 4 GeV when $M_{SUSY}=20~\rm{TeV}$. Nevertheless, the numerical effect of the large logarithms in the red curve is reduced by a factor of around 5 regarding the results shown in Fig. \ref{Fig1}. Finally, we have studied the dependence of the NNLL and three-loop $M_h$ predictions on the stop mixing parameter $X_t$ in the heavy SUSY limit.
\begin{figure}
\centering
\includegraphics[width=18pc]{Plotc.pdf}  \\
\includegraphics[width=18pc]{Plotd.pdf}
\caption{\small{Numerical comparison of the $M_h$ predictions in a scenario where $M_{\tilde{g}}=1.5~\rm{TeV}$ and $X_t/M_{SUSY}=0$. These plots follow the same conventions as in the Figure \ref{Fig1}. Up: Evolution of $M_h$ as a function of $M_{SUSY}$. Down: Differences between the three-loop fixed-order and the FeynHiggs predictions.}}
\label{Fig2}
\end{figure}
In Fig. \ref{Fig3} we increased the value of $X_t/M_{SUSY}$ from 0.2 (thin curves) to 2.4 (thick curves). We observe that for high energy scales above $M_{SUSY} \gtrsim 10~\rm{TeV}$ the agreement between the two predictions and therefore the effect of the large logarithms on the red curves improves when $X_t/M_{SUSY}$ increases up to $1.5$, where the numerical size of the large logarithmic contributions in the $O(\alpha_t\alpha_s^2)$ results is reduced by a factor of about 7 regarding the non-mixing scenario ($X_t=0$). At the kinematic point $X_t/M_{SUSY}=1.5$, the curvature of $M_h$ as a function of $X_t$ changes its sign and therefore $\Delta M_h$ starts to increase again for even higher values ($X_t/M_{SUSY}>1.5$) reaching the maximum difference in the critical mixing $X_t/M_{SUSY}=2.4$ (thickest lines in Fig. \ref{Fig3}) which is another inflection point of $M_h(X_t)$ where the prediction of $M_h$ takes its higher value.
\begin{figure}
\centering
\includegraphics[width=18pc]{Plote.pdf} \\
\vspace*{0.5cm}
\includegraphics[width=18pc]{Plotf.pdf}
\caption{\small{Numerical comparison of the $M_h$ predictions for a non-vanishing stop mixing in the heavy SUSY limit. The blue dashed lines are the NNLL predictions of FeynHiggs and the red solid lines represent our three-loop fixed-order predictions. The brown band is the CMS/ATLAS Higgs boson mass, $M_h^{exp}=125.09 \pm 0.24~\rm{GeV}$. Up: $M_h$ as a function of $M_{SUSY}$ for different stop mixing values, $X_t/M_{SUSY}~=~0.2,~0.5,~0.7,~1.0,~1.5~\rm{and}~2.4$. Down: Absolute numerical differences between the three-loop fixed-order predictions and the NNLO results of FeynHiggs plotted in the upper figure. }}
\label{Fig3}
\end{figure}\\
We further explore the dependence of the Higgs boson mass on the SUSY input parameters $M_{SUSY}$, $X_t$ and $tan\beta$ in the heavy SUSY limit. The figures \ref{Fig1} - \ref{Fig3} show that the predicted value of $M_h$ grows when $M_{SUSY}$ increases and reach a maximum value at the critical point $X_t/M_{SUSY}=2.4$, whose location is independent of $M_{SUSY}$. It suggests that one can find boundaries for the region of rMSSM parameters which put further constraints on $M_h$. Fig. \ref{Fig5} shows the numerical values of $X_t/M_{SUSY}$ and $M_{SUSY}$ which produce the same Higgs mass prediction (gray curves). We have considered values of $M_h$ from $115~\rm{GeV}$ to $131~\rm{GeV}$ and set $tan\beta = 10$. We observed here that there is a minimum value of $M_{SUSY}$, located at the maximal point $X_t/M_{SUSY} = 2.4$, which is compatible with some election of the Higgs boson mass. Moreover, in the case of non stop mixing ($X_t=0$) one can find the higher value of $M_{SUSY}$ compatible with a given $M_h$. These extrema values grow when we consider higher values of $M_h$. This behaviour can also be seen at the intersection of the brown band with the blue dashed lines in Fig. \ref{Fig3} for a $125~\rm{GeV}$ Higgs mass.
\begin{figure}
\centering
\includegraphics[width=18pc]{Plot8.pdf}
\caption{\small{Dependence of $M_h$ on $M_{SUSY}$ and $X_t$ in the heavy SUSY limit. We have used $tan\beta = 10$. The gray lines represent the values of $M_{SUSY}$ and $X_t$ which produce the same Higgs boson mass. The predicted value of $M_h$ increases monotonically with $M_{SUSY}$. }}
\label{Fig5}
\end{figure}
If we use the combined CMS/ATLAS measured Higgs boson mass within the actual combined uncertainties, $M_h^{exp}=125.09 \pm 0.24~\rm{GeV}$, we will be able to fix upper and lower bounds on the SUSY scale $M_{SUSY}$ in the benchmark scenario considered in this work. \\ Fig. \ref{Fig6} shows the $125.09~\rm{GeV}$ contours (gray lines) as a function of $M_{SUSY}$, $tan\beta$ (Up: for values of $X_t/M_{SUSY}$ from 0 to 2.4) and $X_t/M_{SUSY}$ (Down: for values of $tan\beta$ from 4 to 30). The blue and the brown regions refer to the SUSY parameters compatible with $M_h^{exp}$. The purple lines represent the combined uncertainty for the cases enclosed inside.
Notice that if $tan\beta \leq 10$ then $M_{SUSY}$ strongly depends on $tan\beta$, moreover the parameter region of $tan\beta \lesssim 3$ is incompatible with the LHC observations of the Higgs boson mass if one considers SUSY scales below 20 TeV. For values above 10, the dependence is marginal and the curves flatten. As a consequence, at low $tan\beta$ values, independent of the election of $X_t$, it is not possible to find upper bounds on the required SUSY scale from the CMS/ATLAS Higgs mass value. There is still the possibility to use the vacuum stability of the Higgs potential to find upper bounds on $M_{SUSY}$ for small $tan\beta$ values \cite{Allanach} but this is beyond the scope of the present work. For higher values however ($tan\beta\gtrsim 10$), due to the curves are almost constant, one can identify a lower bound for $X_t/M_{SUSY}=2.4$ and an upper bound for a vanishing stop mixing parameter ($X_t=0$). When $tan\beta=10$, which is the point considered in all the above plots of this section, we find that $M_{SUSY}$ must be at most $12.5\pm1.2~\rm{GeV}$ (see purple line in upper plot) in order to be in agreement with the CMS/ATLAS Higgs mass value. $M_{SUSY}$ can be reduced up to $9.6~\rm{GeV}$ for $tan\beta=30$ and $X_t=0$.
\begin{figure}
\centering
\includegraphics[width=18pc]{Plot9.pdf} \\
\vspace*{0.5cm}
\includegraphics[width=18pc]{Plot10.pdf}
\caption{\small{Region of rMSSM parameters in the heavy SUSY limit which is compatible with the central value and the combined uncertainty of the CMS/ATLAS Higgs boson mass, $M_h^{exp}=125.09\pm0.24~\rm{GeV}$. Up: Gray lines represent the points ($M_{SUSY}$, $tan\beta$) compatible with a $125.09~\rm{GeV}$ Higgs mass for different values of the stop mixing parameter, $X_t/M_{SUSY}=0,~0.8,~1.2,~1.5,~2.4$. The purple line represents the combined uncertainty for the case of zero stop mixing. Down: Gray lines are the 125.09 GeV contours as a function of $M_{SUSY}$ and $X_t/M_{SUSY}$ for different values of the parameter $tan\beta$, $tan\beta=4,~5,~10,~30$. The purple lines are the points compatible with the combined uncertainty for the lowest value of $tan\beta$ considered.}}
\label{Fig6}
\end{figure}
One can significantly lower the required value of $M_{SUSY}$ to $1.2~\rm{TeV}$ when $|X_t/M_{SUSY}|$ increases up to $2.4$ and for $tan\beta=30$. The region $M_{SUSY}>12.5\pm1.2~\rm{TeV}$, where the three-loop fixed-order results blow up, is excluded by the combined CMS/ATLAS Higgs mass value in the simple scenario consider here. The coming combined result for RUN 2 by ATLAS and CMS will reduce the current uncertainty and therefore the upper bound on the SUSY scale (for higher values of $tan\beta$) could be reduced even more.

\section{\large Conclusions}
\label{sec-4}

We have recently presented a fixed-order computation of the lightest rMSSM Higgs boson mass which extends the validity of the leading three-loop corrections to the whole parameter space of the rMSSM~\cite{Edilson}. This computation is in a very good agreement with the results of H3m~\cite{Harlander3} for low SUSY scales ($M_{SUSY} \lesssim 1.2~\rm{TeV}$). However for large $M_{SUSY}$ a numerical comparison with the available codes is missing. We have decided to filling this gap by checking our computation of $M_h$ with the three-loop results coming from the EFT hybrid approach implemented in FeynHiggs~2.14~\cite{FeynHiggs3} for the same observable. FeynHiggs includes the resummation of the large logarithms at high SUSY scales and is in a very good agreement with the other fixed-order and EFT codes. This allowed us to compare our results with a reliable three-loop $M_h$-prediction for $M_{SUSY}$ up to 20 TeV. We focused on a single SUSY scale scenario in the decoupling limit (heavy SUSY limit) where the SM is the low energy EFT. We specifically compared our $O(\alpha_t\alpha_s^2)$ and the FeynHiggs NNLL predictions of $M_h$ at the kinematical point $A_{e,\mu,\tau,u,d,c,s,b}=0$, $tan\beta=10$ and $\mu_r=M_{SUSY}$. We find a very good agreement between the two results for SUSY scales below 10~TeV in the case of vanishing stop mixing ($X_t=0$). This agreement can be improved for a different election of the parameters $M_{\tilde{g}}$ and $X_t$. The difference is estimated to be in the range $ 0.2~GeV \lesssim \Delta M_h \lesssim 1~\rm{GeV}$ for the region $M_{SUSY}\lesssim 10~\rm{TeV}$. Above $M_{SUSY}=10~\rm{TeV}$ we have observed significant differences that increase monotonically with $M_{SUSY}$. Such a behaviour is expected for high SUSY scales since the $O(\alpha_t\alpha_s^2)$ computation contain the effects of the large logarithms. Nevertheless, the region where the contributions of the large logarithms blow up is excluded by the combined CMS/ATLAS Higgs mass, $M_h^{exp}=125.09 \pm 0.24~\rm{GeV}$. We have derived an upper bound on the needed SUSY scale for the considered scenario. For values above $tan\beta=10$ the region $M_{SUSY}>12.5\pm1.2~\rm{TeV}$ is ruled out. \\

\vspace*{0.2cm}

\textbf{Acknowledgements.} This work is partially financially supported by the research grant N. 39844 of the call CONVOCATORIA 727 DE COLCIENCIAS PARA DOCTORADOS NACIONALES 2015.


\section{Introduction}
\label{sec:Int}

Exciton-polaritons are elementary excitations of semiconductor
microcavities and constitute hybrid quasi-particles of strongly
coupled light (cavity photons) and matter (quantum well excitons) excitations
\cite{DengRevModPhys2010,CarusottoRevModPhys2010}, retaining the properties
of both constituents. The excitons instill effective interactions inducing
polariton nonlinearities, whereas the small effective mass of photons
enables Bose-Einstein condensation (BEC) of polaritons
at high temperatures~\cite{YamaguchiPRL2013,SchneiderNature2013}.
First experimental observations of exciton-polariton BECs have been
reported more than a decade ago
\cite{KasprzakNature2006,BaliliScience2007}.
Since then, the remarkable properties of the exciton-polariton systems,
combined with their high condensation temperatures \cite{LerarioNatPhys2017},
have motivated much of the research activity in this field.
BECs of exciton polaritons have been studied in various geometries,
e.g., parabolic traps \cite{WoutersPRB2008,TralleroPRB2014},
double-well potentials \cite{WoutersPRL2007,LagoudakisPRL2010,RahmaniSciRep2016,AbbarchiNatPhys2013},
triple wells \cite{ZhangOQE2017}, and 1D and 2D polariton lattices
\cite{AmoPhysique2016,KlembtAPL2017,WhittakerPRL2018,OhadiPRB2018,PanPRB2019}.
In contrast to relatively stable atomic BECs, the polariton BECs
are open quantum systems; they are inherently strongly dissipative
and require continuous pumping, e.g., by external laser fields
\cite{BallariniNatCom2013}.

The dynamics of dissipative systems can be rendered pseudo-Hermitian
in the parity-time ($\mathcal{PT}$) symmetric setup
\cite{Bender1999,Bender2007,ElGanainy2018}, where the gain and loss
are exactly balanced in a complex potential with the reflection symmetric
real part (energy) and reflection antisymmetric imaginary part (gain/loss).
Systems with $\mathcal{PT}$-symmetry have recently attracted much interest
in various branches of physics,
extending from quantum mechanics \cite{Bender2002}
and field theory \cite{Bender2004} to optics \cite{Ganainy2007, Musslimani2008a,Peng2014}
and acoustics \cite{Zhu2014,Fleury2015}.
Since the experimental realization of the spontaneous $\mathcal{PT}$ symmetry
breaking \cite{Guo2009,Ruter2010}, there has been an enormous progress
in the field and a multitude of interesting phenomena have been observed,
such as, e.g., power oscillations \cite{Ruter2010},
double refraction \cite{Makris2010}, and non-reciprocal diffraction \cite{Makris2008}.
The presence of nonlinearity renders $\mathcal{PT}$ symmetric
systems even more remarkable, permitting, for example,
unidirectional \cite{Ramezani2010} and asymmetric \cite{Ambroise2012,Yang2015}
wave propagation in discrete and continuous structures
\cite{Kottos2010,Li2011,Li2013,Zezyulin2012}.

Despite the plethora of different settings, the study of $\mathcal{PT}$
symmetry in polariton structures is still in its infancy and most works
address this issue in the single species context
\cite{LienPRB2015,GaoNature2015}.
A natural next step is to consider polariton mixtures
\cite{SatijaPRA2009,ChestnovSciRep2016,Ohadi2017} in a $\mathcal{PT}$ symmetric
setting where the losses of the polaritons -- due to the semiconductor
exciton recombination and photon escape for the microcavity -- are
compensated by gain from the external laser pumping of the exciton population.
Here we study a two-species polariton mixture with self- and cross-interaction
nonlinearities in a double well structure. We identify the necessary condition
imposed on the pumping rate of the reservoirs and their losses to render
the polariton mixture $\mathcal{PT}$ symmetric. We construct the corresponding
model, analyze its static and dynamic properties, and show that nearly
perfect Rabi-like oscillations of the two polariton components between the
two wells can be observed in this system in the presence of moderate self-
and cross-interactions. This system can then formally be mapped onto
a system of two qubits, or spin-$1/2$ particles, coupled via
exchange ($XY$) interaction. We calculate the fidelity of an effective
 \textsc{swap} gate and quantum state transfer that the system can simulate,
despite being an essentially classical system of coupled BECs.
Our results can thus have interesting implications for
simulations of spin models with polariton lattices.

\section{The exciton-polariton system}
\label{II}

\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig1.pdf}
\end{center}
\caption{Schematic representation of a coupled double well potential
with a binary mixture of exciton-polaritons $\psi_{L,R}^{\pm}$.
Four exciton reservoirs $n_{L,R}^{\pm}$ replenish each polariton
species $(+)$, $(-)$ in each well ($L$), ($R$), respectively.
By balancing the polariton loss $\kappa_{L,R}^{\pm}$ and gain
$Q_{L,R}^{\pm}$ via appropriate reservoir pumping $P_{L,R}^{\pm}$,
we can render the system $\mathcal{PT}$ symmetric.}
\label{fig:fig1}
\end{figure}

We consider a mixture of two species of semiconductor cavity
exciton polaritons in a double-well potential,
as shown schematically in Fig.~\ref{fig:fig1}.
Each of the species $(+)$ or $(-)$ represents a polariton composed
of a circularly right- or left-polarized cavity photon coupled
to a semiconductor exciton transition ($\sigma_+$ or $\sigma_-$) with
the corresponding change of the magnetic quantum number $\Delta m = +1$ or $-1$.
Starting with the mean-field Gross-Pitaevskii equation for
the exciton-polariton BEC in a tight-binding double-well potential,
we follow the standard procedure \cite{RaghavanPRA1999} to derive
the discrete nonlinear Schr\"odinger equation for the four-component
system represented by the wavefunctions $\psi_{L}^{\pm}$ and $\psi_{R}^{\pm}$
of the $(\pm)$ species in the left ($L$) and right ($R$) wells,
respectively,
\begin{subequations}
\label{eqs:EPs}
\begin{eqnarray}
i \hbar \dot{\psi}_{L}^{\pm} &=&
\varepsilon_{L}^{\pm} \psi_{L}^{\pm} + g_s|\psi_{L}^{\pm}|^{2}\psi_{L}^{\pm}
+ g_c|\psi_{L}^{\mp}|^{2}\psi_{L}^{\pm} - J \psi_{R}^{\pm}
\nonumber \\ & &
 + i \frac{\hbar}{2} [Q_{L}^{\pm}(n_{L}^{\pm}) - \kappa_{L}^{\pm} ] \psi_{L}^{\pm} , \\
i \hbar \dot{\psi}_{R}^{\pm} &=&
\varepsilon_{R}^{\pm} \psi_{R}^{\pm} + g_s|\psi_{R}^{\pm}|^{2}\psi_{R}^{\pm}
+ g_c |\psi_{R}^{\mp}|^{2}\psi_{R}^{\pm} - J \psi_{L}^{\pm}
\nonumber \\ & &
 + i\frac{\hbar}{2} [Q_{R}^{\pm}(n_{R}^{\pm}) - \kappa_{R}^{\pm} ] \psi_{R}^{\pm} .
\end{eqnarray}
\end{subequations}
Here $\varepsilon_{L,R}^{\pm}$ are the single-particle energies in each well,
$g_{s}$ and $g_{c}$ are the self- and cross-interaction strengths of the
polaritons, and $J$ is the Josephson tunnel-coupling between the two wells,
assumed the same for both species. The complex wavefunctions can be written
as $\psi_{L,R}^{\pm} = \sqrt{N_{L,R}^{\pm}} \, e^{i \phi_{L,R}^{\pm}}$, where
$N_{L,R}^{\pm} \equiv |\psi_{L,R}^{\pm}|^2$ are the polariton populations
and $\phi_{L,R}^{\pm}$ are the phases.
In contrast to an atomic BEC \cite{RaghavanPRA1999}, the polariton BEC
is a non-conservative system, with the decay $\kappa^{\pm}_{R,L}$
due the exciton recombination losses and cavity photon losses.
The polariton population is continuously replenished via coupling to
exciton reservoirs with rates $Q_{R}^{\pm}(n_{R}^{\pm})$~\cite{WoutersPRL2007}.
In turn, the exciton populations in each reservoir $n_{R,L}^{\pm}$
obey the equations
\begin{subequations}
\label{eqs:reservoirs}
\begin{eqnarray}
\dot{n}_{L}^{\pm} &=& P_{L}^{\pm} - \Gamma_{L}^{\pm} n_{L}^{\pm} -
Q_{L}^{\pm}(n_{L}^{\pm}) N_{L}^{\pm} , \\
\dot{n}_{R}^{\pm} &=& P_{R}^{\pm}-\Gamma_{R}^{\pm} n_{R}^{\pm} -
Q_{R}^{\pm}(n_{R}^{\pm}) N_{R}^{\pm} ,
\end{eqnarray}
\end{subequations}
where $P^{\pm}_{R,L}$ are the rates of exciton creation, usually induced by
a suitable laser pumping \cite{BallariniNatCom2013},
$\Gamma^{\pm}_{L,R}$ are the decay rates of the excitons
in the corresponding reservoir, and $Q^{\pm}_{L,R}(n^{\pm}_{L,R})$ are the rates
of stimulated scattering of the reservoir excitons into the condensate
of $N_{L,R}^{\pm}$ polaritons, see Fig.~\ref{fig:fig1}.
For simplicity, the scattering rate can be approximated
by a linear function of the reservoir exciton population,
$Q^{\pm}_{L,R}(n^{\pm}_{L,R}) \simeq q^{\pm}_{L,R} n^{\pm}_{L,R}$ ~\cite{LienPRB2015}.

The model can be further simplified if we assume sufficiently
large reservoirs such that the exciton populations remain nearly
constant in time, $\dot{n}_{L}^{\pm} \simeq 0$, obtaining $n_{L,R}^{\pm} \simeq
\frac{P_{L,R}^{\pm}}{\Gamma_{L,R}^{\pm} + q_{L,R}^{\pm} N_{L,R}^{\pm}}$.
Assuming also that $q_{L,R}^{\pm} N_{L,R}^{\pm}  \ll \Gamma_{L,R}^{\pm}$,
i.e. the reservoirs are only weakly depleted by the coupling to
the polariton condensates, as compared to their strong pumping and decay,
we have $n_{L,R}^{\pm} = P_{L,R}^{\pm}/\Gamma_{L,R}^{\pm}$. Upon substitution
into Eqs.~(\ref{eqs:EPs}), the polariton gain/loss coefficients
become
\begin{equation}
\gamma_{L,R}^{\pm} = \frac{1}{2}[q_{L,R}^{\pm} P_{L,R}^{\pm}/\Gamma_{L,R}^{\pm} - \kappa_{L,R}^{\pm}] .
\end{equation}
Hence, for given polariton decay rates $\kappa_{L,R}^{\pm}$ -- determined
by photon losses from the cavity and the spontaneous decay (recombination) of the excitons --
and the stimulated scattering $q_{L,R}^{\pm}$ and decay $\Gamma_{L,R}^{\pm}$ rates of the
reservoir excitons, the polariton gain/loss coefficients $\gamma_{L,R}^{\pm}$ can be
precisely tuned by the laser pumping rates $P_{L,R}^{\pm}$.
In turn, the intensity and spatial distribution of the left- and right-circularly polarized radiation
leaking from the cavity is directly proportional to the population $N_{L,R}^{\pm}$ of the corresponding
polariton components and can thus serve for their continuous monitoring \cite{Ohadi2017}.

We note that interactions of the polaritons with the reservoir excitons
\cite{AmoPhysique2016,EstrechoPRB2019} can strongly affect the polariton dynamics, leading
to, e.g., instabilities \cite{DeveaudPhysique2016}, induce the polariton energy shifts
and even be used to optically engineer the polariton landscape \cite{FerrierPRL2011,Ohadi2017}.
But under our assumption of the nearly-constant populations of the exciton reservoirs,
their interaction with the polaritons merely leads to a modification of the tight-binding
trapping potentials and thereby to constant energy shifts which can be absorbed into $\varepsilon_{L,R}^{\pm}$.

A one-dimensional system, to be invariant under $\mathcal{PT}$ transformation
\cite{Bender1999,Bender2007,ElGanainy2018}, should be confined in a complex
potential which has reflection symmetric and reflection antisymmetric real
and imaginary parts, respectively.
For the coupled-mode equations of the polariton system considered
here, these requirements translate to the conditions
\begin{subequations}
\label{eqs:PTcond}
\begin{align}
\varepsilon^{\pm}_{L} &= \varepsilon^{\pm}_{R} , \label{eq:PTcond1} \\
\gamma_{L}^{\pm} &= - \gamma_{R}^{\pm} \label{eq:PTcond2}.
\end{align}
\end{subequations}
We assume that the system is initially prepared with equal number
of particles in each well, $N_{L}^{-} + N_{L}^{+} = N_{R}^{-} + N_{R}^{+} = N$.
We then switch on the tunnel coupling $J$ and increase
the reservoir pumping rate near the left well and reduce it
near the right, in order to balance gain and loss according
to the $\mathcal{PT}$ symmetry condition (\ref{eq:PTcond2}).
The dynamics of the system is governed by the equations
\begin{subequations}
\label{eqs:eps}
\begin{eqnarray}
i \dot{\psi}_{L }^{\pm} &=& g_s |\psi_{L}^{\pm}|^{2} \psi_{L}^{\pm}
+ g_{c} |\psi_{L}^{\mp}|^{2} \psi_{L}^{\pm}  - J \psi_{R}^{\pm}
+ i \gamma \psi_{L}^{\pm} , \quad \\
i \dot{\psi}_{R }^{\pm} &=& g_s |\psi_{R}^{\pm}|^{2} \psi_{R}^{\pm}
+ g_c |\psi_{R}^{\mp}|^{2} \psi_{R}^{\pm}  - J \psi_{L}^{\pm}
- i \gamma \psi_{R}^{\pm},
\end{eqnarray}
\end{subequations}
where we set $\hbar = 1$ and the zero-point energies
$\varepsilon^{\pm}_{L}=\varepsilon^{\pm}_{R}=0$,
and assumed the balanced gain/loss coefficients equal for both species,
$\gamma = \gamma_{L}^{\pm} = - \gamma_{R}^{\pm}$.

\section{$\mathcal{PT}$ symmetry breaking and fixed points}
\label{sec:III}

Let us first review the properties of a linear system with $g_{s}=g_{c}=0$
\cite{Guo2009,Ruter2010}.
The ($+$) and ($-$) polaritons decouple from each other and their equations
become equivalent. The $\mathcal{PT}$ symmetric and broken phases can be
determined from the eigenvalues of the corresponding Hamiltonian matrix,
$\Lambda_{\pm}= \pm \sqrt{J^{2}-\gamma^2}$.
For $\gamma \leq J$, we have the $\mathcal{PT}$ symmetric phase
and the system exhibits pseudo-Hermitian dynamics:
each polariton component coherently oscillates between the two wells,
\begin{gather}
\begin{bmatrix}
\psi_{L}(t) \\
\psi_{R}(t)
\end{bmatrix} =
U
\begin{bmatrix}
\psi_{L}(0) \\
\psi_{R}(0)
\end{bmatrix} , \nonumber \\
U =
\begin{bmatrix}
\cos (\Omega t) + \frac{\gamma}{\Omega} \sin (\Omega t) & i \frac{J}{\Omega} \sin (\Omega t)  \\
i \frac{J}{\Omega} \sin (\Omega t) & \cos (\Omega t) - \frac{\gamma}{\Omega} \sin (\Omega t)
\end{bmatrix} , \label{eq:PTlineval}
\end{gather}
with the effective Rabi frequency $\Omega = (\Lambda_+ -\Lambda_-)/2 = \sqrt{J^{2}-\gamma^2}$.
When $\gamma > J$, the eigenvalues $\Lambda_{\pm}$ become imaginary,
the $\mathcal{PT}$ symmetry breaks, and the population of the particles
in both wells will diverge exponentially.

In the presence of interactions, $g_{s},g_{c} \neq 0$, the $\mathcal{PT}$ phase
diagram cannot be determined from the eigenvalues of non-linear Hamiltonian,
but the relevant information can be extracted from the fixed points and their
stability properties \cite{Li2011,DuanmuPhilTrans2013}.
Assuming $N_{L}^{+} + N_{R}^{+} = N_{L}^{-} + N_{R}^{-} = N$,
we can recast Eqs.~(\ref{eqs:eps}) in terms of the polariton population
imbalances $z^{\pm}= (N_{L}^{\pm} - N_{R}^{\pm})/N$ and phase differences
$\Phi^{\pm}=\phi_{R}^{\pm}-\phi_{L}^{\pm}$ as
\begin{subequations}
\label{eqs:zphi}
\begin{align}
 \dot{z}^{\pm} &= - 2\sqrt{1-(z^{\pm})^2} \sin\Phi^{\pm}+2\frac{\gamma}{J}, \\
\dot{\Phi}^{\pm} &= \frac{g_{s}}{J} z^{\pm} + \frac{g_{c}}{J} z^{\mp}+2\frac{z^{\pm}}{\sqrt{1-(z^{\pm})^2}}\cos\Phi^{\pm},
\end{align}
\end{subequations}
where we use the dimensionless time $\tau = tJ$.
The fixed points correspond to the static solutions of these equations,
$\dot{z}^{\pm}=\dot{\Phi}^{\pm}=0$. By linearizing the above equations
in the vicinity of a fixed point, we find the stability eigenvalues
$\lambda_i$ which determine whether the fixed point is
(i) \textit{stable} if all eigenvalues have negative real parts,
(ii) \textit{unstable} if one or more eigenvalues have positive real part, and
(iii) \textit{elliptic} if all eigenvalues are imaginary.
The dependence of the stability eigenvalues on $\gamma$ thus
serves as the phase diagram of the nonlinear polariton system, since it
shows when the system diverges and when it follows Hermitian-like dynamics.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig2.pdf}
\end{center}
\caption{Real part of the stability eigenvalues $\lambda_{i}$ of
Eq.~(\ref{eq:lambdai}) around the fixed point (\ref{eq:fixed_point_1})
as a function of $\gamma$, for $g_{s} N =J$ and $g_c N =0.9J$.
Both axes are in units of $J$.}
\label{fig:fig2}
\end{figure}

A fixed point of the system corresponds to the trivial solution
\begin{equation}
z_{1}^{\pm} = 0 , \quad \Phi_{1}^{\pm}=\pm\arcsin\left(\frac{\gamma}{J}\right),
\label{eq:fixed_point_1}
\end{equation}
with equal population of the $(+)$ and $(-)$ polaritons in each well.
The analytic expression for the four stability eigenvalues is
\begin{equation}
\lambda_{i} = \pm \frac{\sqrt{-4(J^{2}-\gamma^{2})-2(g_{s} \pm g_{c})\sqrt{(J^{2}-\gamma^{2})}}}{J} . \label{eq:lambdai}
\end{equation}
In Fig.~\ref{fig:fig2} we show the real parts of the stability eigenvalues
$\lambda_{i}$ as a function of $\gamma$, for the system prepared in the vicinity
of the fixed point~(\ref{eq:fixed_point_1}). We observe only elliptic fixed
points, $\mathrm{Re}[\lambda_{i}]=0$, where any small deviation leads to
oscillations around the fixed point; and unstable fixed points,
$\mathrm{Re}[\lambda_{i}] \neq 0$, in the vicinity of which the system diverges.
Equations~(\ref{eq:lambdai}) indicate that as long as $g_{s}>g_{c}$ the
bifurcation to the unstable eigenvalues occurs when $\gamma$ reaches
a threshold value equal to $J$. For the values of $g_{c}$ sufficiently larger
than $g_{s}$ the bifurcation point is shifted to smaller values of $\gamma$.

We have examined the possibility of stable fixed points
with unequal populations of the two wells, corresponding to,
e.g., self-trapping of polariton population on one well
\cite{Albiez2005,AbbarchiNatPhys2013}.
Our analytical and numerical results do not support the existence
of such fixed points, as they break the $\mathcal{PT}$ symmetric
phase of the system \cite{Ramezani2010,GraefeJPA2012}.

\section{Dynamics of the system}
\label{sec:IV}

Before we examine the dynamics of our two-component $\mathcal{PT}$ symmetric
system, it is instructive to consider its Hermitian version with $\gamma=0$.
We assume an initial configuration with all ($+$) particles in the left well
and all ($-$) particles in the right. Figure~\ref{fig:fig3}(a) illustrates
the corresponding dynamics of population imbalances
$z^{\pm}= (N_{L}^{\pm} - N_{R}^{\pm})/N$. We observe that the population
of each species exhibits Rabi-like oscillations between the two wells
with frequency $J$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=\columnwidth]{fig3n.pdf}
\end{center}
\caption{Dynamics of the population imbalances $z^{\pm}$ of two polariton
species in a double well system, for
(a) Hermitian case $\gamma=0$,
(b) non-Hermitian case, with the net gain $\gamma_L=0.1J$ on the left well and
net loss $\gamma_R=-0.2J$ on the right well, and
(c,d,e) $\mathcal{PT}$ symmetric cases with balanced gain/loss
$\gamma_L = - \gamma_R= \gamma=0.1J$ on the two wells, and the interaction strengths
(c) $g_{s}N= J$, $g_{c}N=0.8 J$, (d) $g_{s}N = g_{c}N=gN=J$, and (e) $g_{s}N= J$, $g_{c}N=1.2 J$.
The dash-dotted (gray) line for the amplitude modulation in (b) is $\exp[ (\gamma_L + \gamma_R)t]$,
and in (d) is $\cos(\frac{\gamma g N}{2J} t)$.
Note the different time scales in (a,b) and (c,d,e).}
\label{fig:fig3}
\end{figure}

The exciton polariton system is inherently dissipative, which necessitates
continuous pumping of polaritons from the exciton reservoirs.
In Fig.~\ref{fig:fig3}(b) we illustrate the dynamics of the system
with the net loss from the right well, $\gamma_R=-0.2J$, being larger
than the net gain on the left well, $\gamma_L=0.1J$. As
expected, the polariton populations decay to zero with rate
$\gamma_L + \gamma_R$, as do the population imbalances $z^{\pm}$
normalized to the initial particle number $N$. In the opposite case
of gain being larger than loss, the polariton populations would diverge.

The Hermitian-like dynamics can be achieved with balanced loss and gain,
$\gamma_L = -\gamma_R = \gamma$, as per Eqs.~(\ref{eqs:PTcond}) corresponding
to the $\mathcal{PT}$ symmetric system. Note that the Hermitian system with
$\gamma_{L,R}=0$ can be considered as a special case of the $\mathcal{PT}$
symmetric system with zero imaginary part of the potential.
In Fig.~\ref{fig:fig3}(c,d,e) we illustrate the dynamics of the
interacting system.
In the case of $g_s \geq g_c$, the dynamics is characterized by slow
modulation of the amplitude of the Rabi-like oscillations between the two wells
(see Fig.~\ref{fig:fig3}(c)).
These modulations are not sinusoidal, exhibiting plateaus of maximal amplitude,
followed by rapid collapse and revival of the oscillations. In the case
of $g_s = g_c = g$, the modulation amplitude becomes harmonic, with the period
well approximated by $T = \frac{4\pi J}{\gamma g N}$ (see Fig.~\ref{fig:fig3}(d)).
But for moderate interaction strengths $g_s < g_c$,
the dynamics is nearly indistinguishable from that of the Hermitian system
(see Fig.~\ref{fig:fig3}(e)).

\section{Simulating two-qubit \textsc{swap} gate}
\label{sec:swap}

The two-species polariton system in a $\mathcal{PT}$ symmetric double well
can be formally mapped onto a system of two qubits, or spin-1/2 particles,
coupled via exchange interaction. This is an interesting analogy as it links
an essentially classical system described by coupled-mode equations
to a quantum system of coupled spins.
Here we consider a \textsc{swap} gate to simulate state transfer between two qubits.

The complete basis for a system of two qubits (spins) consists of the four
states $\{ \ket{\downarrow \downarrow}, \ket{\downarrow \uparrow},
\ket{\uparrow \downarrow}, \ket{\uparrow \uparrow} \}$,
and an arbitrary two-qubit state $\ket{\Psi}$ can be expanded as
\begin{equation}
\label{eq:2qubitstate}
\ket{\Psi} = C_{\downarrow \downarrow} \ket{ \downarrow \downarrow}
+ C_{\downarrow \uparrow} \ket{\downarrow \uparrow}
+ C_{\uparrow \downarrow} \ket{\uparrow \downarrow}
+ C_{\uparrow \uparrow}  \ket{\uparrow \uparrow}.
\end{equation}
We may associate the ($+$) component of the polariton in each well
$\psi_{L,R}^{+}$ with the amplitude of the $\ket{\uparrow }$ state
of the spin and the ($-$) component $\psi_{L,R}^{-}$ with the amplitude
of the $\ket{\downarrow}$ state of the spin.
Then, the superposition amplitudes in Eq.~(\ref{eq:2qubitstate}) are given by
$C_{\downarrow \downarrow} = c_{L}^{-} c_{R}^{-}$,
$C_{\downarrow \uparrow} = c_{L}^{-} c_{R}^{+}$,
$C_{\uparrow \downarrow} = c_{L}^{+} c_{R}^{-}$,
$C_{\uparrow \uparrow} = c_{L}^{+} c_{R}^{+}$, where
$c_{L}^{\pm} \equiv \frac{\psi_{L}^{\pm}}
{\sqrt{|\psi_{L}^{+}|^{2} + |\psi_{L}^{-}|^{2}}}$ and similarly for
$c_{R}^{\pm} \equiv \frac{\psi_{R}^{\pm}}
{\sqrt{|\psi_{R}^{+}|^{2} + |\psi_{R}^{-}|^{2}}}$

Our aim is to realize the \textsc{swap} gate between the two spins representing
qubits. Consider the four input states
$\ket{\downarrow \downarrow}, \ket{\downarrow \uparrow},
\ket{\uparrow \downarrow}, \ket{\uparrow \uparrow}$,
which can cast as vectors
\begin{equation}
u_{\downarrow \downarrow}=\begin{bmatrix}
1  \\
0  \\
0  \\
0
\end{bmatrix} , \;
u_{\downarrow \uparrow}=\begin{bmatrix}
0  \\
1  \\
0  \\
0
\end{bmatrix} , \;
u_{\uparrow \downarrow}=\begin{bmatrix}
0  \\
0  \\
1  \\
0
\end{bmatrix} , \;
u_{\uparrow \uparrow}=\begin{bmatrix}
0 \\
0  \\
0  \\
1
\end{bmatrix} .
\end{equation}
In the simplest, non-interacting and Hermitian case,
$g_{s}=g_{c}=\gamma=0$ and $J \neq 0$, the dynamics of the system is
analytically solvable and for each input state the time-dependent
amplitudes of the output state (\ref{eq:2qubitstate}) are given in
Table~\ref{table1}. We observe periodic oscillations of the excitation
between the two spins, or qubits, with frequency $J$. Choosing for the
interaction time $t$ half the period of the oscillations, $tJ = \pi/2$,
yields the transformation
\begin{equation}
\label{swap} U_{\textsc{swap}}=-\begin{bmatrix}
1 & 0 & 0 & 0  \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix} ,
\end{equation}
which is precisely the \textsc{swap} gate (with an overall minus sign).
For comparison, we recall the dynamics of the coupled spin system
governed by the spin-exchange ($XY$) Hamiltonian
$H = - J \hat{\sigma}_L^+ \hat{\sigma}_R^- + \mathrm{H.c.}$,
with $\hat{\sigma}^+ = \ket{\uparrow}\bra{\downarrow}$ and
$\hat{\sigma}^- = \ket{\downarrow}\bra{\uparrow}$, for which the initial
states $\ket{\downarrow \downarrow}$ and $\ket{\uparrow \uparrow}$ remain
unchanged, while the other two basis states evolve as
\begin{eqnarray*}
\ket{\downarrow \uparrow} & \to &\cos(Jt) \ket{\downarrow \uparrow}
 + i \sin(Jt) \ket{\uparrow \downarrow} , \\
\ket{\uparrow \downarrow} & \to &\cos(Jt) \ket{\uparrow \downarrow}
 + i \sin(Jt) \ket{\downarrow \uparrow} .
\end{eqnarray*}
The precise dynamics of quantum spin system differs from that in Table I,
but at time $tJ=\frac{\pi}{2}$ we obtain the equivalent $i$\textsc{swap} gate
\begin{equation}
\label{iswap} U_{i\textsc{swap}}= \begin{bmatrix}
1 & 0 & 0 & 0  \\
0 & 0 & i & 0 \\
0 & i & 0 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix} .
\end{equation}

\begin{table}[t]
\centering
\setlength\tabcolsep{5pt}
\begin{tabular}{|c |c |c |c |c|}
\hlin
$~$ & $u_{\downarrow \downarrow} $ & $u_{\downarrow \uparrow}$ & $u_{\uparrow \downarrow}$ & $u_{\uparrow \uparrow}$    \\ [1.5ex]
\hline
$C_{\downarrow \downarrow} $  & $e^{2iJt}$ &$ i \sin(Jt) \cos(Jt)$ & $ i \sin(Jt) \cos(Jt)$   & 0   \\ [2ex]  \hline
$C_{\downarrow \uparrow}$  & $0$ & $\cos^2(Jt)$ & $-\sin^2(Jt)$ & $0$  \\  [2ex] \hline
$C_{\uparrow \downarrow}$& $0$ & $-\sin^2(Jt)$ & $\cos^2(Jt)$ & $0$ \\  [2ex] \hline
$C_{\uparrow \uparrow}$ & 0 & $ i \sin(Jt) \cos(Jt)$ & $ i \sin(Jt) \cos(Jt)$  & $e^{2iJt}$  \\   [2ex]
\hline
\end{tabular}
\caption{ \label{table1} The time dependent amplitudes
$C_{\downarrow \downarrow},C_{\downarrow \uparrow},C_{\uparrow \downarrow},C_{\uparrow \uparrow}$
of the two-qubit state in Eq.~(\ref{eq:2qubitstate}), for four different
inputs $u_{\downarrow \downarrow},u_{\downarrow \uparrow},u_{\uparrow \downarrow},u_{\uparrow \uparrow}$,
for a non-interacting, Hermitian system, $g_{s}=g_{c}=\gamma=0$.
The table represents the elements of the transformation matrix $U$.}
\end{table}

Consider now the dynamics of the interacting, non-Hermitian $\mathcal{PT}$
symmetric system. For each input state $u_i$ at $t=0$, the output state
$\ket{\Psi_o}$ at time $t >0$ is
\begin{equation}
\label{eq:2qubitstatei}
\ket{\Psi_o^{(i)}}= C^{(i)}_{\downarrow \downarrow} \ket{\downarrow \downarrow} +
C^{(i)}_{\downarrow \uparrow} \ket{\downarrow \uparrow} +
C^{(i)}_{\uparrow \downarrow} \ket{\uparrow \downarrow} +
C^{(i)}_{\uparrow \uparrow} \ket{\uparrow \uparrow} .
\end{equation}
We can construct the transformation matrix $U$ as follows:
For each input state $u_{i}$, the amplitudes of the output state
$C^{(i)}_{o}$ ($i,o \in \{ \downarrow \downarrow, \downarrow \uparrow,
\uparrow \downarrow,\uparrow \uparrow \}$) constitute the $i$th column of $U$.
The complete matrix is then obtained by evaluating the dynamics of the system
for the four input basis states.

We characterize gate transformation using the fidelity~\cite{Molmer2007},
\begin{equation}
\label{fidelity}
F=\frac{1}{n(n+1)} \left[ \textrm{Tr}(MM^{\dagger}) + |\textrm{Tr}(M)|^2 \right],
\end{equation}
where $n = 4$ is the size of the Hilbert space (matrix dimension), and
$M=U_{0}^{\dagger}U$ with $U_{0}$ the unitary transformation corresponding
to the  desired quantum gate, e.g., $U_{\textsc{swap}}$, and $U$ is the actual
transformation that we obtain from our simulations. In other words, for
any initial state $\ket{\Psi}$, the desired final state is $U_{0} \ket{\Psi}$
while the actual state is $U\ket{\Psi}$, and Eq.~(\ref{fidelity}) gives
the fidelity of the transformation averaged over all the input states.

The linear $\mathcal{PT}$ case, $g_s=g_c=0$, is analytically solvable
via Eq.~(\ref{eq:PTlineval}). At half-period $tJ=\frac{\pi}{2}$ we then obtain
the fidelity $F=0.992$ for $\gamma=0.1J$. The small reduction of fidelity is
due to incomplete population transfer,
$\min\limits_{t} |C_{\downarrow \uparrow (\uparrow \downarrow)}| \simeq \frac{\gamma^2}{J^2-\gamma^2}$,
and a slight shift (stretching) of the oscillation frequency by $\gamma$, i.e.,
the effective Rabi frequency is $\Omega \equiv \sqrt{J^2-\gamma^2} \simeq J( 1 - \gamma^2/2J^2)$.
For larger $\gamma=0.3J$, the fidelity at $tJ=\frac{\pi}{2}$ is further reduced
to $F=0.935$. The obtained fidelities can be slightly improved by choosing
appropriately delayed \textsc{swap} time $tJ=\frac{\pi}{2 -\gamma^2/J^2}$
to compensate for the reduction of the effective Rabi frequency of the oscillations.

\begin{figure}[t]
\centerline{\includegraphics[width=1\columnwidth]{fig4.pdf}}
\caption{Time dependence of the complex amplitudes $C_{o}$ of the output
state $\ket{\Psi_o}$ for each input state $u_{i}$, up to time $t=\pi/J$,
for the system with $g_sN=J$, $g_cN=0.5J$ and $\gamma=0.1J$.
In each panel, the blue (solid) line is the absolute value of the
corresponding amplitude, $|C_o^{i}|$, and the red (dashed) line is the phase
$\phi_C = \arg (C_o^{i})$.}
\label{fig:fig4}
\end{figure}

We next consider the interacting system with equal self and cross-interaction
strengths $g_sN=g_cN=J$ and $\gamma=0.1J$. This system is not analytically
solvable and we numerically calculate the time-dependent amplitudes $C^{(i)}_{o}$.
We obtain fidelity $F=0.991$ of the \textsc{swap} gate at time $tJ=\frac{\pi}{2}$.
Further increasing the loss/gain rate decreases the fidelity,
e.g., for $\gamma=0.3J$ we obtain $F =0.922$, and for $\gamma=0.5J$
we obtain $F=0.799$, which can again be improved by using the
optimal \textsc{swap} times $t = \frac{\pi}{2 \Omega}$.

Finally, we examine the general case of unequal self- and cross-interaction
strengths. Figure~\ref{fig:fig4} shows the evolution of the complex amplitudes
$C^{(i)}_{o}$ for each input state $u_{i}$ for the system with
$g_sN=J$, $g_cN=0.5J$ and $\gamma=0.1J$. Each panel of the figure is a direct
visualization of the corresponding time-dependent matrix element of
transformation $U$, up to time $tJ=\pi$ of one complete period of population
oscillation for the corresponding linear, Hermitian system.
At half the period $tJ=\frac{\pi}{2}$, we obtain the  \textsc{swap} fidelity
$F=0.982$, which is only slightly lower than in the previous cases. With
increasing the difference between $g_{s}$ and $g_{c}$ the fidelity decreases,
albeit rather slowly. For instance, for $g_sN=J$, $g_cN=0.1J$ the fidelity is
$F=0.963$ for $\gamma=0.1J$ and $F=0.867$ for $\gamma=0.3J$.

\section{Conclusions}

To summarize, we have considered a two-species polariton mixture with self- and
cross-interaction nonlinearities in a $\mathcal{PT}$ symmetric double well structure and examined
the static and dynamic properties of the system. We have shown that this essentially classical
system described by the mean-field coupled-mode equations can nevertheless mimic quantum
dynamics of exchange-coupled spins implementing quantum gates.

Our work can be relevant to analog simulations of few-body quantum systems with coupled
exciton-polaritons in lattice potentials \cite{AmoPhysique2016,KlembtAPL2017,WhittakerPRL2018,OhadiPRB2018,PanPRB2019}.
Thus, one-dimensional polariton lattices can be used to study the dynamics of a single spin excitation,
e.g. for the characterization of faithful state or excitation transfer in spin chains
with engineered couplings and tunable disorder \cite{Petrosyan2010,Nikolopoulos2014}.
Furthermore, the quantum problem of two interacting particles or spin excitations
on a one-dimensional lattice can be mapped onto an appropriately designed two-dimensional lattice of polaritons,
similarly to the simulations of the one-dimensional two-particle Hubbard model with two-dimensional arrays
of coupled optical waveguides \cite{Krimer2011,Corrielli2013,Rai2015,Mukherjee2016}.
Finally, tunable polariton lattices can serve as versatile simulators of classical spin-lattice models,
such as Ising or $XY$ \cite{LagoudakisNJP2017,BerloffNatMater2017,KalininNJP2018},
which in turn can tackle certain NP-hard optimization problems.

\acknowledgments
We thank E. Paspalakis for fruitful discussions.
This work was co-financed by Greece (General Secretariat for Research and Technology),
and the European Union (European Regional Development Fund), in the framework of
the bilateral Greek-Russian Science and Technology collaboration on Quantum Technologies
(POLISIMULATOR project).