*mps*

#  Probing Quark-Lepton Unification with Leptoquark and Higgs Decays

Pavel Fileviez Pérez<EMAIL_ADDRESS>Physics Department and Center for
Education and Research in Cosmology and Astrophysics (CERCA), Case Western
Reserve University, Cleveland, OH 44106, USA  Elliot Golias
<EMAIL_ADDRESS>Physics Department and Center for Education and
Research in Cosmology and Astrophysics (CERCA), Case Western Reserve
University, Cleveland, OH 44106, USA  Alexis D. Plascencia
<EMAIL_ADDRESS>Physics Department and Center for Education and
Research in Cosmology and Astrophysics (CERCA), Case Western Reserve
University, Cleveland, OH 44106, USA

###### Abstract

We point out unique relations between the decay widths for leptoquarks and
Higgs bosons that can be used to test the unification of quarks and leptons at
the TeV scale. We discuss the main predictions of the minimal theory for
quark-lepton unification and show how the different decays for leptoquarks and
Higgses are related by the symmetry of the theory.

##  I Introduction

After the discovery of the Standard Model (SM) Higgs boson at the Large Hadron
Collider (LHC) we know that the SM of Particle Physics can describe with high
precision the physics below the TeV scale. Nonetheless, there are reasons to
believe that the LHC could discover new forces and a new sector that could
help address some of the open issues in particle physics and cosmology.

It is well-known that the SM needs to be modified in order to explain the
origin of neutrino masses, the nature of dark matter and the baryon asymmetry
in the Universe. Unfortunately, we do not know what is the cut-off scale of
the SM and there is no certainty that the LHC will discover new physics. The
unification of forces in nature is one of the best ideas we have for physics
beyond the Standard Model. The simplest unified theories based on $\,{\rm
SU}(5)$ and $\,{\rm SO}(10)$ can describe physics at the high scale, $M_{\rm
GUT}\sim 10^{15-16}$ GeV, for reviews see Refs. Senjanović (2008); Nath and
Fileviez Perez (2007), and we cannot hope to directly test their predictions
at colliders because the fields are superheavy.

J. Pati and A. Salam Pati and Salam (1974) proposed the idea of matter
unification, where the SM quarks and leptons can be unified in the same
multiplet. This idea was crucial to understand the idea of grand unification.
The minimal Pati-Salam model predicts the same mass for neutrinos and up-
quarks, and the same for down-quarks and charged leptons. Therefore, if we use
the canonical seesaw mechanism Minkowski (1977); Yanagida (1979); Gell-Mann
_et al._ (1979); Mohapatra and Senjanovic (1980) the relevant scale is
$M_{R}\sim 10^{14}$ GeV.

In our view, the generic idea of quark-lepton unification is very appealing.
Some years ago, a simple theory for the unification of quarks and leptons
below the multi-TeV scale was proposed by one of the authors, P. Fileviez
Perez, and M. B. Wise in Ref. Fileviez Perez and Wise (2013). This theory is
based on $\,{\rm SU}(4)_{C}\otimes\,{\rm SU}(2)_{L}\otimes\,{\rm U}(1)_{R}$,
and neutrino masses are generated through the inverse seesaw mechanism
Mohapatra (1986); Mohapatra and Valle (1986) in order to have a consistent
theory where $\,{\rm SU}(4)_{C}$ is broken at the low scale. This theory tells
us that one can hope to test the idea of quark-lepton unification at current
or future colliders.

In this article, motivated by the possibility to test the idea of matter
unification at colliders, we investigate in detail the Higgs and leptoquark
decays in the theory proposed in Ref. Fileviez Perez and Wise (2013). This
theory predicts three leptoquarks and two Higgs doublets. The predicted
leptoquarks are $X_{\mu}\sim(\mathbf{3},\mathbf{1},2/3)_{\text{SM}}$,
$\Phi_{3}\sim(\mathbf{\bar{3}},\mathbf{2},-1/6)_{\text{SM}}$ and
$\Phi_{4}\sim(\mathbf{3},\mathbf{2},7/6)_{\rm SM}$. We point out relations
between the decay widths for Leptoquarks and Higgs bosons that can be used to
test the idea of quark-lepton unification at the TeV scale.

In the case of the scalar leptoquarks we find several unique relations for the
decay widths, for example:

$\displaystyle\sum_{i,j=1}^{3}\Gamma(\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j})$
$\displaystyle=$
$\displaystyle\left(\frac{M_{\phi_{3}^{1/3}}}{M_{\phi_{4}^{2/3}}}\right)\sum_{i,j=1}^{3}\Gamma(\phi_{4}^{2/3}\to\bar{e}_{i}d_{j}).$

Therefore, if these decays are measured it can be checked whether this
relation predicted by quark-lepton unification holds. In the above equation
the leptoquarks $\phi_{3}^{1/3}$ and $\phi_{4}^{2/3}$ are components of the
$\Phi_{3}$ and $\Phi_{4}$ fields, respectively.

The Higgs sector of this theory is special because there are two Higgs
doublets with only four different Yukawa couplings. One finds, for example,
that only two Yukawa couplings determine the masses for the down-quarks and
charged leptons. In the limit where there is a hierarchy between the two Higgs
vacuum expectation values, $\tan\beta=v_{2}/v_{1}\gg 1$, we find that the
quark-lepton unification predicts the following relation for the heavy CP-even
Higgs decay widths

$\sum_{i,j=1}^{3}\Gamma(H\to\bar{d}_{i}d_{j})=3\sum_{i,j=1}^{3}\Gamma(H\to\bar{e}_{i}e_{j}).$

This relation is quite unique because naively we expect the Higgs decays into
leptons and quarks to be different. Finally, we study the relations between
the leptoquark and Higgs decays when $\tan\beta\ll 1$. For example, we find:

$\displaystyle\sum_{i,j=1}^{3}\Gamma(\phi_{3}^{-2/3}\to\bar{d}_{i}e_{j})=4\left(\frac{M_{\phi_{3}^{-2/3}}}{M_{H}}\right)\sum_{i,j=1}^{3}\Gamma(H\to\bar{e}_{i}e_{j}).$

We discuss how the above relations and other relations can be used to test the
idea of quark-lepton unification at colliders.

This article is organized as follows: in Section II we discuss the minimal
theory for quark-lepton unification at the low scale. In Section III we
discuss the properties of the different leptoquarks and their decays. In
Section IV the new Higgses are discussed and we point out the main properties
of their decays into fermions. Our main findings are summarized in Section V.

##  II Minimal Quark-Lepton Unification

A simple renomalizable theory for quark-lepton unification was proposed in
Ref. Fileviez Perez and Wise (2013) which can be seen as a low energy limit of
the Pati-Salam theory. This theory is based on the gauge symmetry, $\,{\rm
SU}(4)_{C}\otimes\,{\rm SU}(2)_{L}\otimes\,{\rm U}(1)_{R}$ and the SM matter
fields are unified in three representations:

$\displaystyle F_{QL}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}u&\nu\\\
d&e\end{array}\right)\sim(\mathbf{4},\mathbf{2},0),$ (3) $\displaystyle F_{u}$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}u^{c}&\nu^{c}\end{array}\right)\sim(\mathbf{\bar{4}},\mathbf{1},-1/2),$
(5) $\displaystyle F_{d}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccc}d^{c}&e^{c}\end{array}\right)\sim(\mathbf{\bar{4}},\mathbf{1},1/2),$
(7)

while the gauge fields live in $A_{\mu}\sim(\mathbf{15},\mathbf{1},0)$. The
minimal Higgs sector has three scalar representations:
$\Phi\sim(\mathbf{15},\mathbf{2},1/2)$, $\chi\sim(\mathbf{4},\mathbf{1},1/2)$
and $H_{1}\sim(\mathbf{1},\mathbf{2},1/2)$. This minimal sector allows us to
write a full realistic theory for matter unification.

The gauge symmetry, $\,{\rm SU}(4)_{C}\otimes\,{\rm SU}(2)_{L}\otimes\,{\rm
U}(1)_{R}$, is spontaneously broken to the SM gauge group by the vacuum
expectation value (VEV) of the scalar field $\chi$, i.e.
$\langle\chi\rangle=\textrm{diag}(0,\,0,\,0,\,v_{\chi}/\sqrt{2})$, which gives
mass to the vector leptoquark $X_{\mu}$, defining the scale of matter
unification. See Appendix A for more details.

The Yukawa interactions for the charged fermions can be written as

$\displaystyle-{\cal L}$ $\displaystyle\supset$ $\displaystyle
Y_{1}\,{F}_{QL}F_{u}H_{1}\ +\ Y_{2}\,{F}_{QL}F_{u}\Phi$ (8) $\displaystyle+$
$\displaystyle Y_{3}\,H_{1}^{\dagger}{F}_{QL}F_{d}\ +\
Y_{4}\,\Phi^{\dagger}{F}_{QL}F_{d}+\mbox{h.c.},$

while for neutrinos one can implement the inverse seesaw mechanism using the
terms

$\displaystyle-{\cal L}$ $\displaystyle\supset$ $\displaystyle Y_{5}F_{u}\chi
S\ +\ \frac{1}{2}\mu SS+\mbox{h.c.}.$ (9)

Here the fields $S\sim(\mathbf{1},\mathbf{1},0)$ are SM fermionic singlets. In
order to achieve very small neutrino masses one needs a seesaw mechanism. The
minimal scenario to have the $\,{\rm SU}(4)_{C}$ symmetry broken at the low
scale and generate small neutrino masses without fine-tuning is using the
inverse seesaw mechanism Fileviez Perez and Wise (2013). This is a key idea
that allows us to realize matter unification below the multi-TeV scale.

The mass matrix for neutrinos in the basis ($\nu$, $\nu^{c}$, $S$) reads as

$\left(\nu\ \nu^{c}\ S\right)\left(\begin{array}[]{ccc}0&M_{\nu}^{D}&0\\\
(M_{\nu}^{D})^{T}&0&M_{\chi}^{D}\\\
0&(M_{\chi}^{D})^{T}&\mu\end{array}\right)\left(\begin{array}[]{c}\nu\\\
\nu^{c}\\\ S\end{array}\right),$ (10)

where $M_{\chi}^{D}=Y_{5}\,v_{\chi}/\sqrt{2}$ and $M_{\nu}^{D}$ is given in
Eq. (25).

The light neutrino masses are given by

$m_{\nu}\approx\mu\,(M_{\nu}^{D})^{2}/(M_{\chi}^{D})^{2},$ (11)

when $M_{\chi}^{D}\gg M_{\nu}^{D}\gg\mu$ and the heavy neutrinos form a
pseudo-Dirac pair.

In order to test the generic idea of quark-lepton unification we need to
understand the predictions of the different theories where this idea is
realized. In this article, we focus on the minimal theory for matter
unification that can describe physics below the multi-TeV scale, because we
can hope to test this idea at current or future colliders. For
phenomenological studies of this theory see Refs. Smirnov (1995); Faber  _et
al._ (2018, 2020); Fileviez Perez  _et al._ (2021).

##  III Leptoquark Decays

The theory discussed in the previous section predicts a vector leptoquark,
$X_{\mu}\sim(\mathbf{3},\mathbf{1},2/3)_{\text{SM}}$, associated to the
$\,{\rm SU}(4)_{C}$ symmetry, and four physical scalar leptoquarks. The scalar
leptoquarks $\Phi_{3}\sim(\mathbf{\bar{3}},\mathbf{2},-1/6)_{\text{SM}}$ and
$\Phi_{4}\sim(\mathbf{3},\mathbf{2},7/6)_{\rm SM}$ can be written in $\,{\rm
SU}(2)_{L}$ components as,

$\Phi_{3}=\begin{pmatrix}\phi_{3}^{1/3}\\\\[4.30554pt]
\phi_{3}^{-2/3}\end{pmatrix},\hskip 14.22636pt\text{and}\hskip
14.22636pt\Phi_{4}=\begin{pmatrix}\phi_{4}^{5/3}\\\\[4.30554pt]
\phi_{4}^{2/3}\end{pmatrix},$ (12)

where the numbers in the superscript denote the electric charge. The Yukawa
interactions for $\Phi_{3}$ and $\Phi_{4}$ are given by

$\displaystyle-\mathcal{L}$ $\displaystyle\supset$ $\displaystyle
Y_{2}\,\varepsilon_{ab}\,\ell_{L}^{a}\,\Phi_{4}^{b}\,(u^{c})_{L}+Y_{2}\,\varepsilon_{ab}\,Q_{L}^{a}\,\Phi_{3}^{b}\,(\nu^{c})_{L}$
(13) $\displaystyle+$ $\displaystyle
Y_{4}\,\Phi_{3}^{\dagger}\,\ell_{L}\,(d^{c})_{L}+Y_{4}\,\Phi_{4}^{\dagger}\,Q_{L}\,(e^{c})_{L}+\mbox{h.c.}\,,$

where $a$ and $b$ correspond to the $\,{\rm SU}(2)_{L}$ indices. Notice that
in this sector we only have two different Yukawa couplings because the $\,{\rm
SU}(4)_{C}$ symmetry relates the different Yukawa interactions in a unique
way.

The physical scalar leptoquarks in this theory are:

$\phi_{3}^{1/3},\,\,\,\phi_{4}^{5/3},\,\,\,\phi_{3}^{-2/3}\,\,\,\textrm{and}\,\,\,\phi_{4}^{2/3}.$

The leptoquarks $\phi_{3}^{-2/3}$ and $\phi_{4}^{2/3}$ can mix but the mixing
angle is determined by the electroweak scale, and hence, it is generically
very small. Consequently, in this work we ignore this mixing. For the
interactions of these fields see Appendix E.

Now, let us discuss the different decays of the leptoquarks and their decays
to understand how the quark-lepton symmetry predicts unique relations between
the different decay widths.

$X_{\mu}$ decays: The vector leptoquark $X_{\mu}$ can have the following
decays

$X_{\mu}\to\bar{e}_{i}d_{j},\,\,\bar{\nu}_{i}u_{j},$

where $i,j=1,2,3$ correspond to the family indices. Unfortunately, naively one
expects that the vector leptoquark mass must be above $10^{3}$ TeV to satisfy
the experimental bounds on rare decays such as $K_{L}\to e^{\pm}\mu^{\mp}$
Valencia and Willenbrock (1994), unless one uses the freedom on the mixings
between quarks and leptons. For the vector leptoquark there is a very simple
relation for the decay widths:

$\displaystyle\sum_{i.j=1}^{3}\Gamma(X_{\mu}\to\bar{e}_{i}d_{j})$
$\displaystyle=$ $\displaystyle
2\sum_{i,j=1}^{3}\Gamma(X_{\mu}\to\bar{\nu}_{i}u_{j}).$ (14)

Notice that the total widths for a given decay channel are clean, meaning that
they are independent of the unknown mixing angles between quarks and leptons
determined by the matrix $V_{DE}$ defined in Appendix E. If the right-handed
neutrinos are much lighter than the $X_{\mu}$ boson, the decay widths of
$X_{\mu}$ into quarks and leptons are equal.

Scalar Leptoquarks: In the case of the scalar leptoquarks one can have the
decays

$\displaystyle\phi_{3}^{1/3}$ $\displaystyle\to$
$\displaystyle\bar{d}_{i}\nu_{j},\,\bar{d}_{i}N_{j};\hskip
34.14322pt\phi_{4}^{5/3}\to\bar{e}_{i}u_{j};$ $\displaystyle\phi_{3}^{-2/3}$
$\displaystyle\to$
$\displaystyle\bar{d}_{i}e_{j},\,\bar{u}_{i}\nu_{j},\,\bar{u}_{i}N_{j};\,\,\,\,\,\phi_{4}^{2/3}\to\bar{e}_{i}d_{j},\,\bar{\nu}_{i}u_{j},\,\bar{N}_{i}u_{j}.$

Even though the Feynman rules for the scalar leptoquarks listed in Appendix E
are involved, it is possible to define some total widths for different
channels that are independent of the unknown mixing angles entering in the
interactions. For example, by defining:

$\displaystyle\Gamma_{T}(\phi_{3}^{1/3}\to\bar{d}\nu)$
$\displaystyle\equiv\sum_{i,j=1}^{3}\Gamma(\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j})$
$\displaystyle=\frac{3M_{\phi_{3}^{1/3}}}{16\pi}{\rm
Tr}[Y_{4}^{\dagger}Y_{4}],$ (15)

which turns out to be independent of any mixing angle when the fermion masses
are neglected, see Appendix B for the details of the calculation. In order to
understand this simple result, notice that the individual decay width
$\Gamma(\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j})$ is a complicated function of
mixing matrices and Yukawa couplings. See the Feynman rules in Appendix E for
the details. However, since the leptoquarks are heavy we can neglect the
fermion masses and by summing over the family index we find a total decay
width that is a function only of the Yukawa couplings. This simple idea allows
us to find relations between the different decays.

We investigated the different leptoquark decays into fermions and found the
following set of relations predicted by quark-lepton unification:

$\displaystyle\frac{\Gamma_{T}(\phi_{3}^{1/3}\to\bar{d}\nu)}{M_{\phi_{3}^{1/3}}}$
$\displaystyle=$
$\displaystyle\frac{\Gamma_{T}(\phi_{4}^{2/3}\to\bar{e}d)}{M_{\phi_{4}^{2/3}}}$
(16) $\displaystyle=$
$\displaystyle\frac{\Gamma_{T}(\phi_{3}^{-2/3}\to\bar{d}e)}{M_{\phi_{3}^{-2/3}}}.$

Now, if the decay channels with the heavy pseudo-Dirac neutrinos are
available, neglecting all fermion masses one finds:

$\displaystyle\frac{\Gamma_{T}(\phi_{3}^{1/3}\to\bar{d}N)}{M_{\phi_{3}^{1/3}}}$
$\displaystyle=$
$\displaystyle\frac{\Gamma_{T}(\phi_{4}^{2/3}\to\bar{\nu}u)}{M_{\phi_{4}^{2/3}}}$
(17) $\displaystyle=$
$\displaystyle\frac{\Gamma_{T}(\phi_{3}^{-2/3}\to\bar{u}N)}{M_{\phi_{3}^{-2/3}}},$
$\displaystyle\Gamma_{T}(\phi_{4}^{5/3}\to\bar{e}u)$ $\displaystyle=$
$\displaystyle\frac{M_{\phi_{4}^{5/3}}}{M_{\phi_{3}^{-2/3}}}\left[\Gamma_{T}(\phi_{3}^{-2/3}\to\bar{u}N)\right.$
(18) $\displaystyle\left.+\,\Gamma_{T}(\phi_{3}^{-2/3}\to\bar{d}e)\right].$

Furthermore, if the pseudo-Dirac neutrinos are light and when $\phi_{3}^{1/3}$
and $\phi_{4}^{5/3}$ are the lightest elements of the $\Phi_{3}$ and
$\Phi_{4}$ respectively; then, their total widths will be dominated by the
$\ell+q$ decays, and hence, their lifetimes will satisfy

$M_{\phi_{3}^{1/3}}\,\,\tau_{\phi_{3}^{1/3}}=M_{\phi_{4}^{5/3}}\,\,\tau_{\phi_{4}^{5/3}}\,\,.$
(19)

Similarly, if $\phi_{3}^{1/3}$ and $\phi_{4}^{2/3}$ are the lightest elements
of the $\Phi_{3}$ and $\Phi_{4}$ we then have

$M_{\phi_{3}^{1/3}}\,\,\tau_{\phi_{3}^{1/3}}=M_{\phi_{4}^{2/3}}\,\,\tau_{\phi_{4}^{2/3}}\,\,.$
(20)

The same is true when $\phi_{3}^{-2/3}$ and either $\phi_{4}^{2/3}$ or
$\phi_{4}^{5/3}$ are the lightest elements

$\displaystyle
M_{\phi_{3}^{-2/3}}\,\,\tau_{\phi_{3}^{-2/3}}=M_{\phi_{4}^{2/3}}\,\,\tau_{\phi_{4}^{2/3}}\,\,,$
(21) $\displaystyle
M_{\phi_{3}^{-2/3}}\,\,\tau_{\phi_{3}^{-2/3}}=M_{\phi_{4}^{5/3}}\,\,\tau_{\phi_{4}^{5/3}}\,\,.$
(22)

Clearly, these relations are predictions from the unification of quark and
leptons and can be used to test this idea at particle colliders.

It is well-known that any theory, including the SM, does not predict the
values of the gauge and Yukawa couplings present in the interactions, but can
predict relations between the different physical quantities. In this case, the
minimal theory for quark-lepton unification is predicting a set of relations
for the decay widths that can be tested if these leptoquarks are discovered in
the near future.

It is important to mention that the branching ratios for the leptoquark decays
depend of the mass splittings of the different components of $\Phi_{3}$ and
$\Phi_{4}$. For example, $\phi_{3}^{1/3}$ can have the following decays with
the $\phi_{3}^{-2/3}$ in the final state: $\phi_{3}^{1/3}\to
W^{+}\phi_{3}^{-2/3},\,\,\pi^{+}\phi_{3}^{-2/3},\,\,\bar{e}_{i}\nu_{j}\phi_{3}^{-2/3}$.
Only when the mass splitting is small one can a large branching ratio for
$\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j}$. For a study about the relation between
the decays of a $\,{\rm SU}(2)_{L}$ doublet leptoquark see the studies in Ref.
Fileviez Perez  _et al._ (2009). See also Ref. Doršner  _et al._ (2016) for a
review about leptoquarks and Ref. Murgui and Wise (2021) for a recent
discussion about proton decay mediated by scalar Leptoquarks in this theory,
where the authors have shown that there are no dimension five contributions to
proton decay.

##  IV Higgs Decays

This theory predicts a unique Higgs sector with two Higgs doublets,
$H_{1}\sim(\mathbf{1},\mathbf{2},1/2)$ and
$H_{2}\sim(\mathbf{1},\mathbf{2},1/2)$, with only four independent Yukawa
couplings. The Yukawa couplings for the Higgses can be written as

$\displaystyle-\mathcal{L}_{Y}$ $\displaystyle=$
$\displaystyle\bar{u}_{R}\left(Y_{1}^{T}H_{1}+\frac{1}{2\sqrt{3}}Y_{2}^{T}H_{2}\right)Q_{L}$
(23) $\displaystyle+$
$\displaystyle\bar{N}_{R}\left(Y_{1}^{T}H_{1}-\frac{\sqrt{3}}{2}Y_{2}^{T}H_{2}\right)\ell_{L}$
$\displaystyle+$
$\displaystyle\bar{d}_{R}\left(Y_{3}^{T}H_{1}^{\dagger}+\frac{1}{2\sqrt{3}}Y_{4}^{T}H_{2}^{\dagger}\right)Q_{L}$
$\displaystyle+$
$\displaystyle\bar{e}_{R}\left(Y_{3}^{T}H_{1}^{\dagger}-\frac{\sqrt{3}}{2}Y_{4}^{T}H_{2}^{\dagger}\right)\ell_{L}+{\rm
h.c.}\,.$

Here we neglect the small mixing between the Higgs doublets and the $\chi$
field. Notice that in the general two Higgs doublet model (commonly referred
in the literature as the type-III 2HDM) there are eight different Yukawa
couplings, but in our case the $\,{\rm SU}(4)_{C}$ symmetry relates quarks and
leptons so there are only four independent Yukawa couplings. For reviews on
two-Higgs doublet models we refer the reader to Refs. Gunion  _et al._ (2000);
Branco  _et al._ (2012).

After symmetry breaking the charged fermions and the Dirac neutrino masses are
given by

$\displaystyle M_{U}$ $\displaystyle=$ $\displaystyle
Y_{1}\frac{v_{1}}{\sqrt{2}}+\frac{1}{2\sqrt{3}}Y_{2}\frac{v_{2}}{\sqrt{2}},$
(24) $\displaystyle M_{\nu}^{D}$ $\displaystyle=$ $\displaystyle
Y_{1}\frac{v_{1}}{\sqrt{2}}-\frac{\sqrt{3}}{2}Y_{2}\frac{v_{2}}{\sqrt{2}},$
(25) $\displaystyle M_{D}$ $\displaystyle=$ $\displaystyle
Y_{3}\frac{v_{1}}{\sqrt{2}}+\frac{1}{2\sqrt{3}}Y_{4}\frac{v_{2}}{\sqrt{2}},$
(26) $\displaystyle M_{E}$ $\displaystyle=$ $\displaystyle
Y_{3}\frac{v_{1}}{\sqrt{2}}-\frac{\sqrt{3}}{2}Y_{4}\frac{v_{2}}{\sqrt{2}}.$
(27)

Here the VEVs of the Higgs doublets are defined as $\langle
H^{0}_{1}\rangle=v_{1}/\sqrt{2}$, and $\langle
H^{0}_{2}\rangle=v_{2}/\sqrt{2}$. Notice that the above equations allow us to
have different masses for charged leptons and down-quarks. However, one needs
the inverse seesaw mechanism to generate small neutrino masses discussed in
Section II.

In our convention the mass matrices are diagonalized as

$\displaystyle U^{T}M_{U}U_{c}=M_{U}^{\rm diag},$ (28) $\displaystyle
D^{T}M_{D}D_{c}=M_{D}^{\rm diag},$ (29) $\displaystyle
E^{T}M_{E}E_{c}=M_{E}^{\rm diag}.$ (30)

Notice that the above relations tell us that the Yukawa coupling $Y_{4}$
defines the difference between $M_{E}$ and $M_{D}$, and one can use these
relations to write the decay widths for leptoquarks as a function of quark
masses. For example,

$\displaystyle\Gamma_{T}(\phi_{4}^{2/3}\to\bar{e}d)$ $\displaystyle=$
$\displaystyle\frac{3M_{\phi_{4}^{2/3}}}{16\pi}{\rm
Tr}[Y_{4}^{\dagger}Y_{4}]\simeq\frac{9m_{b}^{2}\,M_{\phi_{4}^{2/3}}}{32\pi\,v_{2}^{2}}.$
(31)

As in any two-Higgs doublet model, assuming CP-conservation, the physical
fields are: $h$ and $H$ the CP-even neutral fields, $A$ the CP-odd field, and
two charged Higgs bosons $H^{\pm}$. The $h$ field is identified as the SM-like
Higgs boson. We list all the Feynman rules in Appendix F. As in the case of
the leptoquarks, one can find a set of relations between the decay widths of
the Higgs bosons.

The physical Higgs fields are defined as

$\displaystyle\begin{pmatrix}H\\\
h\end{pmatrix}=\begin{pmatrix}\cos\alpha&\sin\alpha\\\
-\sin\alpha&\cos\alpha\end{pmatrix}\begin{pmatrix}H_{1}^{0}\\\
H_{2}^{0}\end{pmatrix},$ (32) $\displaystyle\begin{pmatrix}G^{0}\\\
A^{0}\end{pmatrix}=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}A_{1}^{0}\\\
A_{2}^{0}\end{pmatrix},$ (33) $\displaystyle\begin{pmatrix}G^{\pm}\\\
H^{\pm}\end{pmatrix}=\begin{pmatrix}\cos\beta&\sin\beta\\\
-\sin\beta&\cos\beta\end{pmatrix}\begin{pmatrix}H_{1}^{\pm}\\\
H_{2}^{\pm}\end{pmatrix},$ (34)

where $H^{0}_{i}$, $H_{i}^{\pm}$, $A_{i}^{0}$ are the neutral, charged, and
CP-odd components of the Higgs doublets, respectively, and $G^{0}$, $G^{\pm}$
are the Goldstone bosons. Furthermore, the mixing angle $\beta$ is related to
the VEVs of the Higgs doublets by $\tan{\beta}=v_{2}/v_{1}$. Since the $h$
field is the SM-like Higgs boson, in order to agree with the measured
properties of the SM Higgs field we will work in the limit
$\sin(\beta-\alpha)\to 1$ (or $\alpha\simeq\beta-\pi/2$). In this limit we
make sure the couplings of $h$ are SM-like.

The physical Higgs can have the usual decays to matter:

$\displaystyle H,\,A$
$\displaystyle\to\bar{u}_{i}u_{j},\,\,\bar{d}_{i}d_{j},\,\,\bar{e}_{i}e_{j},\,\,\bar{\nu}_{i}N_{j};$
$\displaystyle H^{+}$
$\displaystyle\to\bar{e}_{i}\nu_{j},\,\,\bar{e}_{i}N_{j},\,\,\bar{d}_{i}u_{j}.$

The decay widths of these Higgses are a function of the Yukawa couplings and
unknown mixing angles.

Following the same idea used for the leptoquark decays, we can define the
total widths for the different channels and then find simple relations between
the Higgs decays. In the limit of $\tan\beta\gg 1$, quark-lepton unification
predicts the following unique relation for the heavy Higgs

$\sum_{i,j=1}^{3}\Gamma(H\to\bar{d}_{i}d_{j})=3\sum_{i,j=1}^{3}\Gamma(H\to\bar{e}_{i}e_{j}).$
(35)

Notice that this is a striking relation for the decays into down-quarks and
charged leptons.

We can find the following relation for the decay widths of the charged Higgs,
$H^{\pm}$, if the mass of the right-handed neutrinos, $M_{N_{i}}$, is much
smaller than $M_{H^{\pm}}$:

$\Gamma_{T}(H^{+}\to\bar{d}u)=3\left[\Gamma_{T}(H^{+}\to\bar{e}\nu)+\Gamma_{T}(H^{+}\to\bar{e}N)\right].$
(36)

A similar relation can be obtained for the CP-odd neutral Higgs:

$\sum_{i,j=1}^{3}\Gamma(A\to\bar{d}_{i}d_{j})=3\sum_{i,j=1}^{3}\Gamma(A\to\bar{e}_{i}e_{j}).$
(37)

In the $\tan\beta\ll 1$ limit the relations between the Higgs decays become

$\displaystyle\Gamma_{T}(H\to\bar{e}e)$ $\displaystyle=$ $\displaystyle 3\
\Gamma_{T}(H\to\bar{d}d),$ (38) $\displaystyle\Gamma_{T}(A\to\bar{e}e)$
$\displaystyle=$ $\displaystyle 3\ \Gamma_{T}(A\to\bar{d}d),$ (39)

and in this limit it is possible to relate the decay widths of leptoquarks
with the ones from the Higgs scalars as follows

$\displaystyle\Gamma_{T}(\phi_{4}^{2/3}\to\bar{\nu}u)$
$\displaystyle=\frac{12\,M_{\phi_{4}^{2/3}}}{M_{H}}\,\Gamma_{T}(H\to\bar{u}u),$
(40) $\displaystyle\Gamma_{T}(\phi_{3}^{-2/3}\to\bar{d}e)$
$\displaystyle=4\left(\frac{M_{\phi_{3}^{-2/3}}}{M_{H}}\right)\Gamma_{T}(H\to\bar{e}e).$
(41)

The following relation holds independent of the value of $\tan\beta$

$\Gamma_{T}(A\to\bar{e}e)=\frac{M_{A}}{M_{H}}\Gamma_{T}(H\to\bar{e}e).$ (42)

From our perspective, these relations are very unique. Notice that the
relations in Eqs. (35) - (42) tell us something interesting: the theory
predicts simple relations for the total decay widths into quarks and leptons.
Clearly, these are predictions from the unification of quarks and leptons. We
are not aware of any model for physics beyond the SM that can predict these
relations for the Higgs and leptoquark decays.

##  V Summary

The idea of quark-lepton unification is one of the best motivated ideas for
physics beyond the Standard Model. We have discussed the minimal gauge theory
for quark-lepton unification that can describe physics below the multi-TeV
scale. This theory predicts one vector leptoquark, four scalar leptoquarks and
a unique Higgs sector with only four independent Yukawa couplings.

In this article we pointed out unique relations between the decay widths for
leptoquarks that can be used to test the idea of quark-lepton unification. We
also discussed the Higgs sector of the theory. The theory has two Higgs
doublets with only four different Yukawa couplings determined by the symmetry
between quarks and leptons. We discussed the different Higgs decays and showed
how the total Higgs decay widths into quarks and leptons are related at large
(and small) values for the ratio between the vacuum expectation values. We
believe that these results should motivate new studies to test the idea of
matter unification at the LHC or future colliders.

Acknowledgments: The work of P.F.P. has been supported by the U.S. Department
of Energy, Office of Science, Office of High Energy Physics, under Award
Number DE-SC0020443.

##  Appendix A Gauge and Higgs Fields

In this theory the SM gluon fields, the vector leptoquark and the new neutral
gauge boson live in the adjoint representation of $\,{\rm SU}(4)_{C}$:

$\displaystyle A_{\mu}=\left(\begin{array}[]{cc}G_{\mu}&X_{\mu}/\sqrt{2}\\\
X_{\mu}^{*}/\sqrt{2}&0\\\ \end{array}\right)+T_{4}\
B_{\mu}^{{}^{\prime}}\sim(\mathbf{15},\mathbf{1},0),$ (45)

where $G_{\mu}\sim(\mathbf{8},\mathbf{1},0)_{\text{SM}}$ are the gluons,
$X_{\mu}\sim(\mathbf{3},\mathbf{1},2/3)_{\text{SM}}$ are vector leptoquarks,
and $B_{\mu}^{{}^{\prime}}\sim(\mathbf{1},\mathbf{1},0)_{\text{SM}}$. The
Higgs sector is composed of

$\displaystyle H^{T}_{1}$ $\displaystyle=$ $\displaystyle\left(H^{+}_{1}\
H^{0}_{1}\right)\sim(\mathbf{1},\mathbf{2},1/2),\,\,\,\,\,\,\chi=\left(\chi_{u}\
\chi_{R}^{0}\right)\sim(\mathbf{4},\mathbf{1},1/2),\,\,\,\,\,\,{\rm{and}}$
$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\Phi_{8}&\Phi_{3}\\\ \Phi_{4}&0\\\
\end{array}\right)+\sqrt{2}\,T_{4}\ H_{2}\sim(\mathbf{15},\mathbf{2},1/2).$
(48)

Here $H_{2}\sim(\mathbf{1},\mathbf{2},1/2)_{\text{SM}}$ is a second Higgs
doublet, $\Phi_{8}\sim(\mathbf{8},\mathbf{2},1/2)_{\text{SM}}$, and the scalar
leptoquarks $\Phi_{3}\sim(\mathbf{\bar{3}},\mathbf{2},-1/6)_{\text{SM}}$ and
$\Phi_{4}\sim(\mathbf{3},\mathbf{2},7/6)_{\rm SM}$. The $T_{4}$ generator of
$\,{\rm SU}(4)_{C}$ in the above equation is normalized as
$T_{4}=\frac{1}{2\sqrt{6}}\rm{diag}(1,1,1,-3).$

##  Appendix B Decay Widths

Here we present the details on how the expressions for the decay widths
simplify when we sum over family indices. Let’s consider the decay
$\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j}$ and sum over the family indices,

$\displaystyle\Gamma_{T}(\phi_{3}^{1/3}\to\bar{d}\nu)$
$\displaystyle\equiv\sum_{i,j=1}^{3}\Gamma(\phi_{3}^{1/3}\to\bar{d}_{i}\nu_{j})=\frac{3M_{\phi_{3}^{1/3}}}{16\pi}(N^{*}Y_{4}^{\dagger}D_{c}^{\dagger})^{ij}(NY_{4}^{T}D_{c}^{T})^{ij}$
$\displaystyle=\frac{3M_{\phi_{3}^{1/3}}}{16\pi}N^{\dagger\alpha
i}N^{i\sigma}D_{c}^{\dagger\beta
j}D_{c}^{j\rho}Y_{4}^{\dagger\alpha\beta}Y_{4}^{\rho\sigma}=\frac{3M_{\phi_{3}^{1/3}}}{16\pi}\delta^{\alpha\sigma}\delta^{\beta\rho}Y_{4}^{\dagger\alpha\beta}Y_{4}^{\rho\sigma}$
$\displaystyle=\frac{3M_{\phi_{3}^{1/3}}}{16\pi}{\rm
Tr}[Y_{4}^{\dagger}Y_{4}],$ (49)

due to the unitarity of mixing the matrices, the final result turns out to be
independent of any mixing parameter when the fermion masses are neglected. The
same simplification occurs for the decays of the vector leptoquark $X_{\mu}$
after summing over the family index.

##  Appendix C Heavy Gauge Boson Masses

The VEV of the scalar $\chi$ corresponding to
$\langle\chi\rangle=\textrm{diag}(0,\,0,\,0,\,v_{\chi}/\sqrt{2})$ is
responsible for the spontaneous breaking of $\,{\rm SU}(4)_{C}\otimes\,{\rm
SU}(2)_{L}\otimes\,{\rm U}(1)_{R}\to\,{\rm SU}(3)_{C}\otimes\,{\rm
SU}(2)_{L}\otimes\,{\rm U}(1)_{Y}$. After this spontaneous breaking, the gauge
boson corresponding to the $\,{\rm SU}(4)_{C}$ generator $T_{4}$
($B_{\mu}^{{}^{\prime}}$) mixes with the $\,{\rm U}(1)_{R}$ gauge boson
($Z^{\mu}_{R}$), which results in the massless $\,{\rm U}(1)_{Y}$ gauge boson
($B_{\mu}$) of the Standard Model and an orthogonal massive state
($Z^{{}^{\prime}}_{\mu}$) defined by

$\displaystyle\begin{pmatrix}Z_{R\mu}\\\
B^{{}^{\prime}}_{\mu}\end{pmatrix}=\begin{pmatrix}\cos\theta_{4}&\sin\theta_{4}\\\
-\sin\theta_{4}&\cos\theta_{4}\end{pmatrix}\begin{pmatrix}Z^{{}^{\prime}}_{\mu}\\\
B_{\mu}\end{pmatrix},$ (50)

where the mixing angle and the hypercharge gauge coupling are given by

$\displaystyle\sin{\theta_{4}}=\frac{g_{4}}{\sqrt{g_{4}^{2}+\frac{2}{3}g_{R}^{2}}}\hskip
14.22636pt{\rm and}\hskip
14.22636ptg_{Y}=\frac{g_{R}g_{4}}{\sqrt{g_{4}^{2}+\frac{2}{3}g_{R}^{2}}}.$
(51)

The mass of the heavy neutral gauge boson reads as

$\displaystyle
M_{Z^{{}^{\prime}}}^{2}\simeq\frac{1}{4}\left(g_{R}^{2}+\frac{3}{2}g_{4}^{2}\right)v_{\chi}^{2},$
(52)

while the mass of the vector leptoquark is given by

$M_{X}^{2}\simeq\frac{g_{4}^{2}}{4}v_{\chi}^{2}.$ (53)

##  Appendix D Scalar Potential and Leptoquark Masses

The scalar potential may be written as

$\displaystyle V$ $\displaystyle\supset$
$\displaystyle\mu^{2}_{H_{1}}H_{1}^{\dagger}{H_{1}}+\mu_{\chi}^{2}\chi^{\dagger}\chi+\mu_{\Phi}^{2}{\rm
Tr}\left[\Phi^{\dagger}\Phi\right]+\lambda_{1}H_{1}^{\dagger}{H_{1}}\chi^{\dagger}\chi+\lambda_{2}H_{1}^{\dagger}H_{1}{\rm
Tr}\left[\Phi^{\dagger}\Phi\right]+\lambda_{3}\chi^{\dagger}\chi{\rm
Tr}\left[\Phi^{\dagger}\Phi\right]$ $\displaystyle+$
$\displaystyle\left(\lambda_{4}H_{1}^{\dagger}\chi^{\dagger}\Phi\chi+\mbox{h.c.}\right)+\lambda_{5}H_{1}^{\dagger}{\rm
Tr}\left[\Phi\,\Phi^{\dagger}\right]H_{1}+\lambda_{6}\chi^{\dagger}\Phi\Phi^{\dagger}\chi+\lambda_{7}\left(H_{1}^{\dagger}H_{1}\right)^{2}+\lambda_{8}\left(\chi^{\dagger}\chi\right)^{2}$
$\displaystyle+$ $\displaystyle\lambda_{9}{\rm
Tr}\left[(\Phi^{\dagger}\Phi)^{2}\right]+\lambda_{10}\left({\rm
Tr}[\Phi^{\dagger}\Phi]\right)^{2}+\left(\lambda_{11}H_{1}^{a\dagger}{\rm
Tr}[\Phi^{a}\Phi^{b}]H_{1}^{b\dagger}+\lambda_{12}H_{1}^{a\dagger}{\rm
Tr}\left[\Phi^{a}\Phi^{b}\Phi^{b\dagger}\right]\right.$ $\displaystyle+$
$\displaystyle\left.\lambda_{13}H_{1}^{a\dagger}{\rm
Tr}\left[\Phi^{a}\Phi^{b\dagger}\Phi^{b}\right]+\mbox{h.c.}\right)+\lambda_{14}\chi^{\dagger}\Phi^{\dagger}\Phi\chi+\lambda_{15}{\rm
Tr}\left[\Phi^{a\dagger}\Phi^{b}\Phi^{b\dagger}\Phi^{a}\right]+\lambda_{16}{\rm
Tr}\left[\Phi^{a\dagger}\Phi^{b}\right]{\rm
Tr}\left[\Phi^{b\dagger}\Phi^{a}\right]$ $\displaystyle+$
$\displaystyle\lambda_{17}{\rm
Tr}\left[\Phi^{a\dagger}\Phi^{b\dagger}\right]{\rm
Tr}\left[\Phi^{a}\Phi^{b}\right]+\lambda_{18}{\rm
Tr}\left[\Phi^{a\dagger}\Phi^{b\dagger}\Phi^{a}\Phi^{b}\right]+\lambda_{19}{\rm
Tr}\left[\Phi^{a\dagger}\Phi^{b\dagger}\Phi^{b}\Phi^{a}\right]\,,$

where the trace is in $\,{\rm SU}(4)_{C}$ space, and $a,b$ are $\,{\rm
SU}(2)_{L}$ indices. Taking the limit $v_{\chi}\gg v_{1},v_{2}$ we find that

$\displaystyle M_{\Phi_{8}}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{3}\lambda_{4}}{4}\cot{\beta}-\frac{3}{8}(\lambda_{6}+\lambda_{14})\right)v_{\chi}^{2},$
(55) $\displaystyle M_{\Phi_{3}}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{3}\lambda_{4}}{4}\cot{\beta}+\frac{\lambda_{14}-3\lambda_{6}}{8}\right)v_{\chi}^{2},$
(56) $\displaystyle M_{\Phi_{4}}^{2}$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{3}\lambda_{4}}{4}\cot{\beta}+\frac{\lambda_{6}-3\lambda_{14}}{8}\right)v_{\chi}^{2},$
(57) $\displaystyle M_{H_{2}}^{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}\lambda_{4}}{4}\cot{\beta}\,v_{\chi}^{2}.$ (58)

These mass expressions imply the following tree-level sum rule Faber  _et al._
(2018):

$M_{\Phi_{8}}^{2}+2M_{H_{2}}^{2}=\frac{3}{2}\left(M_{\Phi_{3}}^{2}+M_{\Phi_{4}}^{2}\right),$
(59)

which implies that the masses of the scalars in the theory are related and of
the same scale.

##  Appendix E Leptoquark Feynman Rules

The vector leptoquark has the following interactions

$-\mathcal{L}\supset\frac{g_{4}}{\sqrt{2}}X_{\mu}(\bar{Q}_{L}\gamma^{\mu}\ell_{L}+\bar{u}_{R}\gamma^{\mu}\nu_{R}+\bar{d}_{R}\gamma^{\mu}e_{R})+{\rm
h.c.}\,.$ (60)

For the scalar leptoquarks, after expanding the $\,{\rm SU}(2)_{L}$ indices of
the Yukawa interactions in Eq. (13), we can describe the interactions by

$\displaystyle-\mathcal{L}\supset$
$\displaystyle\,\,Y_{4}^{ji}\left[\bar{d}_{R}^{i}\left(\phi_{3}^{1/3}\right)^{*}\nu_{L}^{j}+\bar{d}_{R}^{i}\left(\phi_{3}^{-2/3}\right)^{*}e_{L}^{j}+\bar{e}_{R}^{i}\left(\phi_{4}^{5/3}\right)^{*}u_{L}^{j}+\bar{e}_{R}^{i}\left(\phi_{4}^{2/3}\right)^{*}d_{L}^{j}\right]$
$\displaystyle+Y_{2}^{ji}\left[\bar{N}^{i}\,\phi_{3}^{-2/3}\,u_{L}^{j}-\bar{N}^{i}\,\phi_{3}^{1/3}\,d_{L}^{j}+\bar{u}_{R}^{i}\,\phi_{4}^{2/3}\,\nu_{L}^{j}-\bar{u}_{R}^{i}\,\phi_{4}^{5/3}\,e_{L}^{j}\right],$
(61)

where $i,j$ correspond to family indices. The mass eigenstates
$\phi_{A}^{-2/3}$ and $\phi_{B}^{-2/3}$ are defined as

$\displaystyle\phi_{3}^{-2/3}$ $\displaystyle=$
$\displaystyle\cos\theta_{\text{LQ}}\,\phi_{A}^{-2/3}+\sin\theta_{\text{LQ}}\,\phi_{B}^{-2/3},$
(62) $\displaystyle(\phi_{4}^{2/3})^{*}$ $\displaystyle=$
$\displaystyle-\sin\theta_{\text{LQ}}\,\phi_{A}^{-2/3}+\cos\theta_{\text{LQ}}\,\phi_{B}^{-2/3}.$
(63)

The masses of the leptoquarks have to be above 1 TeV from LHC constraints and
the mixing is determined by the electroweak scale; therefore, the mixing angle
$\theta_{\text{LQ}}$ is very small in general.

The Feynman rules for the interactions of the leptoquarks with fermions are as
follows:

  * •

$X_{\mu}$:

$\displaystyle\bar{d}^{i}e^{j}X_{\mu}:\hskip
14.22636pti\frac{g_{4}}{\sqrt{2}}\left[(V_{DE})^{ij}P_{R}+(V_{2})^{ji}P_{L}\right]\gamma^{\mu},$
(64) $\displaystyle\bar{u}^{i}\nu^{j}X_{\mu}:\hskip
14.22636pti\frac{g_{4}}{\sqrt{2}}\left(K_{1}V_{\rm CKM}K_{2}V_{DE}K_{3}V_{\rm
PMNS}\right)^{ij}\gamma^{\mu}P_{L},$ (65)
$\displaystyle\bar{u}^{i}N^{j}X_{\mu}:\hskip
14.22636pti\frac{g_{4}}{\sqrt{2}}V_{1}^{ji}\gamma^{\mu}P_{R}.$ (66)

  * •

$\phi_{3}^{1/3}$:

$\displaystyle\bar{\nu}^{i}d^{j}\phi_{3}^{1/3}:\hskip
14.22636pti(V^{*}_{4})^{ij}\,P_{R},$ (67)
$\displaystyle\bar{N}^{i}d^{j}\phi_{3}^{1/3}:\hskip
14.22636pt-i(V_{3}^{T}\,K_{1}\,V_{\rm CKM}\,K_{2})^{ij}\,P_{L}.\hskip
142.26378pt$ (68)

  * •

$\phi_{4}^{5/3}$:

$\displaystyle\bar{u}^{i}e^{j}\phi_{4}^{5/3}:\hskip
14.22636pti\left[-(V_{5}^{T}\,V^{\dagger}_{\rm
PMNS}\,K_{3}^{*})^{ij}\,P_{L}+(V_{6}^{*})^{ij}\,P_{R}\right].\hskip
120.92421pt$ (69)

  * •

$\phi_{A}^{-2/3}$:

$\displaystyle\bar{N}^{i}u^{j}\phi_{A}^{-2/3}:\hskip
14.22636pti\cos{\theta_{\rm LQ}}(V_{3}^{T})^{ij}\,P_{L},$ (70)
$\displaystyle\bar{\nu}^{i}u^{j}\phi_{A}^{-2/3}:\hskip
14.22636pt-i\sin{\theta_{\rm LQ}}(V_{5}^{*})^{ij}\,P_{R},$ (71)
$\displaystyle\bar{e}^{i}d^{j}\phi_{A}^{-2/3}:\hskip
14.22636pti\left(\cos{\theta_{\rm LQ}}(K_{3}V_{\rm
PMNS}V_{4}^{*})^{ij}P_{R}-\sin{\theta_{\rm LQ}}(V_{6}^{T}K_{1}V_{\rm
CKM}K_{2})^{ij}P_{L}\right).$ (72)

  * •

$\phi_{B}^{-2/3}$:

$\displaystyle\bar{N}^{i}u^{j}\phi_{B}^{-2/3}:\hskip
14.22636pti\sin{\theta_{\rm LQ}}(V_{3})^{ji}\,P_{L},$ (73)
$\displaystyle\bar{\nu}^{i}u^{j}\phi_{B}^{-2/3}:\hskip
14.22636pti\cos{\theta_{\rm LQ}}(V_{5}^{*})^{ij}P_{R},$ (74)
$\displaystyle\bar{e}^{i}d^{j}\phi_{B}^{-2/3}:\hskip
14.22636pti\left(\sin{\theta_{\rm LQ}}(K_{3}V_{\rm
PMNS}V_{4}^{*})^{ij}P_{R}+\cos{\theta_{\rm LQ}}(V_{6}^{T}K_{1}V_{\rm
CKM}K_{2})^{ij}P_{L}\right),$ (75)

where $P_{L,R}$ are the chiral projection operators
$P_{L,R}=(1\mp\gamma^{5})/2$ and the interaction matrices are given by

$\displaystyle V_{1}$ $\displaystyle=N^{\dagger}_{c}U_{c},\hskip
59.75095ptV_{2}=E^{\dagger}_{c}D_{c},\hskip 48.36958ptV_{3}=U^{T}Y_{2}N_{c},$
(76) $\displaystyle V_{4}$ $\displaystyle=N^{T}Y_{4}D_{c},\hskip
45.52458ptV_{5}=N^{T}Y_{2}U_{c},\hskip 39.83368ptV_{6}=U^{T}Y_{4}E_{c},$
$\displaystyle U^{\dagger}D$ $\displaystyle=K_{1}V_{\rm CKM}K_{2},\hskip
19.91684ptE^{\dagger}N=K_{3}V_{\rm PMNS},\hskip
28.45274ptV_{DE}=D^{\dagger}E,$

$K_{1}$ and $K_{3}$ are diagonal matrices containing three phases, $K_{2}$ is
a diagonal matrix with two phases.

##  Appendix F Higgs Feynman Rules

The Feynman for the physical scalars in the two Higgs doublets rules
correspond to:

  * •

$H^{+}$:

$\displaystyle\bar{u}^{i}d^{j}H^{+}:\hskip
14.22636pti\left[(C_{Lud})^{ij}P_{L}+(C_{Rud})^{ij}P_{R}\right],$ (77)
$\displaystyle\bar{N}^{i}e^{j}H^{+}:\hskip 14.22636pti(C_{Ne})^{ij}P_{L},$
(78) $\displaystyle\bar{\nu}^{i}e^{j}H^{+}:\hskip 14.22636pt-i(C_{\nu
e})^{ij}P_{R}.$ (79)

  * •

$H^{-}$:

$\displaystyle\bar{d}^{i}u^{j}H^{-}:\hskip
14.22636pti\left[(C_{Rud}^{*})^{ji}P_{L}+(C_{Lud}^{*})^{ji}P_{R}\right],$ (80)
$\displaystyle\bar{e}^{i}N^{j}H^{-}:\hskip 14.22636pti(C_{Ne}^{*})^{ji}P_{R},$
(81) $\displaystyle\bar{e}^{i}\nu^{j}H^{-}:\hskip 14.22636pt-i(C_{\nu
e}^{*})^{ji}P_{L}.$ (82)

  * •

$h$:

$\displaystyle\bar{u}^{i}u^{j}h:\hskip
14.22636pti\left[(C_{uu}^{h})^{ij}P_{L}+(C_{uu}^{h*})^{ji}P_{R}\right],$ (83)
$\displaystyle\bar{N}^{i}\nu^{j}h:\hskip
14.22636pti\left[(C_{N\nu}^{h})^{ij}P_{L}+(C_{N\nu}^{h*})^{ji}P_{R}\right],$
(84) $\displaystyle\bar{d}^{i}d^{j}h:\hskip
14.22636pti\left[(C_{dd}^{h})^{ij}P_{L}+(C_{dd}^{h*})^{ji}P_{R}\right],$ (85)
$\displaystyle\bar{e}^{i}e^{j}h:\hskip
14.22636pti\left[(C_{ee}^{h})^{ij}P_{L}+(C_{ee}^{h*})^{ji}P_{R}\right].$ (86)

  * •

$H$:

$\displaystyle\bar{u}^{i}u^{j}H:\hskip
14.22636pti\left[(C_{uu}^{H})^{ij}P_{L}+(C_{uu}^{H*})^{ji}P_{R}\right],$ (87)
$\displaystyle\bar{N}^{i}\nu^{j}H:\hskip
14.22636pti\left[(C_{N\nu}^{H})^{ij}P_{L}+(C_{N\nu}^{H*})^{ji}P_{R}\right],$
(88) $\displaystyle\bar{d}^{i}d^{j}H:\hskip
14.22636pti\left[(C_{dd}^{H})^{ij}P_{L}+(C_{dd}^{H*})^{ji}P_{R}\right],$ (89)
$\displaystyle\bar{e}^{i}e^{j}H:\hskip
14.22636pti\left[(C_{ee}^{H})^{ij}P_{L}+(C_{ee}^{H*})^{ji}P_{R}\right].$ (90)

  * •

$A$:

$\displaystyle\bar{u}^{i}u^{j}A:\hskip
14.22636pt(C_{uu}^{A})^{ij}P_{L}-(C_{uu}^{A*})^{ji}P_{R}\,,$ (91)
$\displaystyle\bar{N}^{i}\nu^{j}A:\hskip
14.22636pt(C_{N\nu}^{A})^{ij}P_{L}-(C_{N\nu}^{A*})^{ji}P_{R}\,,$ (92)
$\displaystyle\bar{d}^{i}d^{j}A:\hskip
14.22636pt(C_{dd}^{A})^{ij}P_{L}-(C_{dd}^{A*})^{ji}P_{R}\,,$ (93)
$\displaystyle\bar{e}^{i}e^{j}A:\hskip
14.22636pt(C_{ee}^{A})^{ij}P_{L}-(C_{ee}^{A*})^{ji}P_{R}\,.$ (94)

where the interaction matrices are given by

$\displaystyle C_{Lud}$ $\displaystyle=U_{c}^{T}\left(Y_{1}^{T}\sin\beta-
Y_{2}^{T}\frac{\cos\beta}{2\sqrt{3}}\right)D,\hskip
31.2982ptC_{Rud}=-U^{\dagger}\left(Y_{3}^{*}\sin\beta-
Y_{4}^{*}\frac{\cos\beta}{2\sqrt{3}}\right)D_{c}^{*},$ $\displaystyle C_{Ne}$
$\displaystyle=N_{c}^{T}\left(Y_{1}^{T}\sin\beta+Y_{2}^{T}\frac{\sqrt{3}\cos\beta}{2}\right)E,\hskip
22.76228ptC_{\nu
e}=N^{\dagger}\left(Y_{3}^{*}\sin\beta+Y_{4}^{*}\frac{\sqrt{3}\cos\beta}{2}\right)E_{c}^{*},$
$\displaystyle C_{uu}^{H}$
$\displaystyle=U^{T}_{c}\left(Y_{1}^{T}\frac{\cos\alpha}{\sqrt{2}}+Y_{2}^{T}\frac{\sin\alpha}{2\sqrt{6}}\right)U,\hskip
36.98866ptC_{N\nu}^{H}=N_{c}^{T}\left(Y_{1}^{T}\frac{\cos\alpha}{\sqrt{2}}-Y_{2}^{T}\frac{3\sin\alpha}{2\sqrt{6}}\right)N,$
$\displaystyle C_{dd}^{H}$
$\displaystyle=D_{c}^{T}\left(Y_{3}^{T}\frac{\cos\alpha}{\sqrt{2}}+Y_{4}^{T}\frac{\sin\alpha}{2\sqrt{6}}\right)D,\hskip
36.98866ptC_{ee}^{H}=E_{c}^{T}\left(Y_{3}^{T}\frac{\cos\alpha}{\sqrt{2}}-Y_{4}^{T}\frac{3\sin\alpha}{2\sqrt{6}}\right)E,$
$\displaystyle C_{uu}^{h}$
$\displaystyle=U^{T}_{c}\left(-Y_{1}^{T}\frac{\sin\alpha}{\sqrt{2}}+Y_{2}^{T}\frac{\cos\alpha}{2\sqrt{6}}\right)U,\hskip
27.03003ptC_{N\nu}^{h}=N_{c}^{T}\left(-Y_{1}^{T}\frac{\sin\alpha}{\sqrt{2}}-Y_{2}^{T}\frac{3\cos\alpha}{2\sqrt{6}}\right)N,$
$\displaystyle C_{dd}^{h}$
$\displaystyle=D_{c}^{T}\left(-Y_{3}^{T}\frac{\sin\alpha}{\sqrt{2}}+Y_{4}^{T}\frac{\cos\alpha}{2\sqrt{6}}\right)D,\hskip
28.45274ptC_{ee}^{h}=E_{c}^{T}\left(-Y_{3}^{T}\frac{\sin\alpha}{\sqrt{2}}-Y_{4}^{T}\frac{3\cos\alpha}{2\sqrt{6}}\right)E,$
$\displaystyle C_{uu}^{A}$
$\displaystyle=U_{c}^{T}\left(Y_{1}^{T}\frac{\sin\beta}{\sqrt{2}}-Y_{2}^{T}\frac{\cos\beta}{2\sqrt{6}}\right)U,\hskip
28.45274ptC_{N\nu}^{A}=N_{c}^{T}\left(Y_{1}^{T}\frac{\sin\beta}{\sqrt{2}}+Y_{2}^{T}\frac{3\cos\beta}{2\sqrt{6}}\right)N,$
$\displaystyle C_{dd}^{A}$
$\displaystyle=D_{c}^{T}\left(-Y_{3}^{T}\frac{\sin\beta}{\sqrt{2}}+Y_{4}^{T}\frac{\cos\beta}{2\sqrt{6}}\right)D,\hskip
28.45274ptC_{ee}^{A}=E_{c}^{T}\left(-Y_{3}^{T}\frac{\sin\beta}{\sqrt{2}}-Y_{4}^{T}\frac{3\cos\beta}{2\sqrt{6}}\right)E.$

## References

  * Senjanović (2008) G. Senjanović, “Course on grand unification,” Springer Proc. Phys. 118, 137–179 (2008).
  * Nath and Fileviez Perez (2007) P. Nath and P. Fileviez Perez, “Proton stability in grand unified theories, in strings and in branes,” Phys. Rept. 441, 191–317 (2007), arXiv:hep-ph/0601023 .
  * Pati and Salam (1974) J. C. Pati and A. Salam, “Lepton Number as the Fourth Color,” Phys. Rev. D 10, 275–289 (1974), [Erratum: Phys.Rev.D 11, 703–703 (1975)].
  * Minkowski (1977) P. Minkowski, “$\mu\to e\gamma$ at a Rate of One Out of $10^{9}$ Muon Decays?” Phys. Lett. B 67, 421–428 (1977).
  * Yanagida (1979) T. Yanagida, “Horizontal gauge symmetry and masses of neutrinos,” Conf. Proc. C 7902131, 95–99 (1979).
  * Gell-Mann  _et al._ (1979) M. Gell-Mann, P. Ramond, and R. Slansky, “Complex Spinors and Unified Theories,” Conf. Proc. C 790927, 315–321 (1979), arXiv:1306.4669 [hep-th] .
  * Mohapatra and Senjanovic (1980) R. N. Mohapatra and G. Senjanovic, “Neutrino Mass and Spontaneous Parity Nonconservation,” Phys. Rev. Lett. 44, 912 (1980).
  * Fileviez Perez and Wise (2013) P. Fileviez Perez and M. B. Wise, “Low Scale Quark-Lepton Unification,” Phys. Rev. D 88, 057703 (2013), arXiv:1307.6213 [hep-ph] .
  * Mohapatra (1986) R. N. Mohapatra, “Mechanism for Understanding Small Neutrino Mass in Superstring Theories,” Phys. Rev. Lett. 56, 561–563 (1986).
  * Mohapatra and Valle (1986) R. N. Mohapatra and J. W. F. Valle, “Neutrino Mass and Baryon Number Nonconservation in Superstring Models,” Phys. Rev. D 34, 1642 (1986).
  * Smirnov (1995) A. D. Smirnov, “The Minimal quark - lepton symmetry model and the limit on Z-prime mass,” Phys. Lett. B 346, 297–302 (1995), arXiv:hep-ph/9503239 .
  * Faber  _et al._ (2018) T. Faber, M. Hudec, M. Malinsky, P. Meinzinger, W. Porod, and F. Staub, “A unified leptoquark model confronted with lepton non-universality in $B$-meson decays,” Phys. Lett. B 787, 159–166 (2018), arXiv:1808.05511 [hep-ph] .
  * Faber  _et al._ (2020) T. Faber, M. Hudec, H. Kolešová, Y. Liu, M. Malinsky, W. Porod, and F. Staub, “Collider phenomenology of a unified leptoquark model,” Phys. Rev. D 101, 095024 (2020), arXiv:1812.07592 [hep-ph] .
  * Fileviez Perez  _et al._ (2021) P. Fileviez Perez, C. Murgui, and A. D. Plascencia, “Leptoquarks and matter unification: Flavor anomalies and the muon g-2,” Phys. Rev. D 104, 035041 (2021), arXiv:2104.11229 [hep-ph] .
  * Valencia and Willenbrock (1994) G. Valencia and S. Willenbrock, “Quark - lepton unification and rare meson decays,” Phys. Rev. D 50, 6843–6848 (1994), arXiv:hep-ph/9409201 .
  * Fileviez Perez  _et al._ (2009) P. Fileviez Perez, T. Han, T. Li, and M. J. Ramsey-Musolf, “Leptoquarks and Neutrino Masses at the LHC,” Nucl. Phys. B 819, 139–176 (2009), arXiv:0810.4138 [hep-ph] .
  * Doršner  _et al._ (2016) I. Doršner, S. Fajfer, A. Greljo, J. F. Kamenik, and N. Košnik, “Physics of leptoquarks in precision experiments and at particle colliders,” Phys. Rept. 641, 1–68 (2016), arXiv:1603.04993 [hep-ph] .
  * Murgui and Wise (2021) C. Murgui and M. B. Wise, “Scalar leptoquarks, baryon number violation, and Pati-Salam symmetry,” Phys. Rev. D 104, 035017 (2021), arXiv:2105.14029 [hep-ph] .
  * Gunion  _et al._ (2000) J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson,  _The Higgs Hunter’s Guide_ , Vol. 80 (2000).
  * Branco  _et al._ (2012) G. C. Branco, P. M. Ferreira, L. Lavoura, M. N. Rebelo, M. Sher, and J. P. Silva, “Theory and phenomenology of two-Higgs-doublet models,” Phys. Rept. 516, 1–102 (2012), arXiv:1106.0034 [hep-ph] .


# Quantum thermodynamics of driven-dissipative condensates

Luísa Toledo Tude and Paul R. Eastham  School of Physics, Trinity College
Dublin, Dublin 2, Ireland  Trinity Quantum Alliance, Unit 16, Trinity
Technology and Enterprise Centre, Pearse Street, Dublin 2, Ireland

###### Abstract

Polariton condensates occur away from thermal equilibrium, in an open system
where heat and particles are continually exchanged with reservoirs. These
phenomena have been extensively analyzed in terms of kinetic equations. Based
on the collection of knowledge about polariton kinetics provided by these
simulations and by experimental works, we constructed a few-level model that
captures the main processes involved in the buildup of a ground-state
population of polaritons. This allows condensation to be understood as the
output of a heat engine and exposes the thermodynamic constraints on its
occurrence. The model consists of a three-level system interacting with a
field and connected to a hot and a cold thermal reservoir that represent a
non-resonant pump and the lattice phonons. This subsystem can drive a
condensate, through polariton-polariton scattering, which produces work in the
form of coherent light emission from the microcavity. We obtain a phase
diagram as a function of the temperatures of the two baths and investigate the
possible types of phase transition that lead to the condensate phase.

##  I Introduction

Bose-Einstein condensates are ordered states which show macroscopic quantum
effects like superfluidity. While originally an equilibrium phenomenon, there
are now many situations where condensation occurs in the non-equilibrium
steady-states of driven open systems. Examples include condensates of exciton-
polaritons in inorganic [1] and organic semiconductors [2], photons in dye-
filled cavities [3, 4], plasmon-polaritons [5], and magnons [6]. These
particles have a finite lifetime, and the condensate must be maintained
against losses by gain from an external pump.

Some condensates, such as those of photons and magnons, quickly reach a quasi-
equilibrium state, and so can to an extent be treated using equilibrium
thermodynamics. For others, however, this description does not apply. Progress
has been made by studying the dynamics and steady-states of different models,
using kinetic equations and field-theoretic approaches, among other
techniques. This diversity of approaches makes it difficult to identify the
general requirements for condensation. Here, we suggest a universal
description of non-equilibrium condensation as a heat engine. We show how
condensation can be described using a three-level laser model, whose
connection to the thermodynamics of heat engines, first pointed out by Scovil
and Schulz-DuBois, has been extensively studied [7, 8, 9, 10, 11, 12]. For
definiteness, we focus on the exciton-polariton condensate, and use the
collection of knowledge provided by theoretical and experimental works to
construct the three-level model. We investigate the thermodynamic properties
such as efficiency, and determine the conditions required for condensation. An
important result is that the occurrence of condensation is determined by two
temperatures whose difference controls the direction of energy flow in the
system, governing the formation of the condensate, and allowing the losses to
be overcome (Fig. 4). Our work further clarifies that non-equilibrium
condensation requires a population inversion, albeit of an unconventional
kind, and emphasizes the need for an effective coupling to a cold reservoir.
This last requirement – the need for a rapid depopulation of the lower state –
is well known in the context of lasers, but is rarely explicit in the
literature on driven-dissipative condensates.

The thermodynamics of three-level amplifiers has been studied before, and our
main contribution is to apply these results to driven-dissipative
condensation. However, this application does require two extensions of the
general framework. Firstly, we consider not just amplification, but the
balance of gain and loss in the full system comprising a condensate driven by
a gain medium. This allows us to compute non-equilibrium phase diagrams and
study how thermodynamic quantities such as efficiency vary across them.
Secondly, we allow for the possibility that the heat engine operates between
reservoirs at different chemical potentials, as well as different
temperatures, implying that it absorbs work as well as heat. This is necessary
in practice for driven-dissipative condensates, as discussed below, but also
relevant in cases such as electrically-driven lasers, where the voltage bias
is a source of work.

##  II Condensates as heat engines

Driven-dissipative condensates are typically understood as open quantum
systems, in which energy from an incoherent pump reservoir is converted into a
coherent condensate, which is in turn emitted into an environment. As
discussed further below, we argue that if the pump is an incoherent source it
corresponds to a heat bath, while the coherent emission from the condensate is
a source of work. With these identifications we identify a condensate as a
form of heat engine, converting heat, from the pump, into work. It follows
immediately, as a consequence of the second law of thermodynamics, that a
consistent description of condensation requires consideration of a cold
reservoir, in addition to the hot reservoir representing the pump.

To develop this idea further we construct a few-level model of condensation.
For definiteness, we consider inorganic microcavities, where the processes
involved in forming and maintaining the condensate are generally accepted. As
illustrated in Fig. 1(a), particles are created by a high-energy pump, which
produces electron-hole pairs that populate the exciton states at high
momentum. This population quickly begins to thermalize with the emission of
acoustic phonons [13], but such scattering becomes ineffective for the
polaritonic final states at low energy and momentum. As confirmed by kinetic
simulations [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26], this
bottleneck effect [27] can be overcome by polariton-polariton scattering. If
the density is large enough this scattering into low-energy states can exceed
their loss. It then becomes a stimulated process which leads to condensation
[28].

Figure 1: (a) Illustration of the main processes that lead to condensation of
inorganic microcavity polaritons. A non-resonant pump creates excitons that
relax to lower energy states through phonon emission. Pair-wise polariton
scattering leads to occupation in the lower energy modes if it can overcome
radiative losses. (b) The condensation process is modeled as the outcome of a
three-level heat engine with a hot bath connecting the ground state and the
most energetic state and a cold bath connecting the ground state to the middle
state. The condensate is a classical field interacting with the two excited
states of the working medium. (c) Model of a microcavity comprising $M$ three-
level systems that contribute to growth of the condensate, which competes with
loss due to the finite polariton lifetime.

To determine the thermodynamic constraints on condensation we consider a
minimal model, consisting of a heat engine whose working medium is a three-
level system with energies $\epsilon_{1}$, $\epsilon_{2}$, and $\epsilon_{3}$.
It can be argued that this is the simplest possible case [29]. As shown in
Fig. 1(b), we suppose that the highest energy state of the working medium,
$\ket{3}$, corresponds to a population of $n$ polaritons in the bottleneck
region. These can undergo pairwise scattering, adding one particle to the
condensate and promoting one polariton to a higher-energy exciton state. We
shall treat the condensate as a field external to the working medium, so this
takes the latter to an intermediate state, $\ket{2}$. The high-energy exciton
can then transition back to the bottleneck region by phonon emission, leading
to the state $\ket{1}$ with bottleneck population $n-1$. In the terminology of
the polariton parametric oscillator [28], the populated bottleneck states are
the pump states, and the high-energy exciton state generated by pairwise
scattering is the idler.

The three-level working medium drives the polariton condensate, which we treat
as a classical mean field. The relevant Hamiltonian is

$H_{s}=g_{c}(p_{c}^{\dagger}p_{i}^{\dagger}p_{p}p_{p}+p_{c}p_{i}p_{p}^{\dagger}p_{p}^{\dagger}),$
(1)

where $p_{c}$, $p_{p}$, and $p_{i}$ are the annihilation operators for
particles in the condensate state, pump state, and idler state, respectively,
and $g_{c}$ is the interaction strength. We consider a macroscopically
occupied condensate mode and approximate $g_{c}p_{c}\approx g_{c}\langle
p_{c}\rangle=e^{-i\omega t}\Omega/2$, where $\omega$ is the condensate
frequency, and for a condensate of $N$ particles

$\Omega=2g_{c}\sqrt{N}.$ (2)

In the basis of the three-level system, the final three terms in each product
in Eq. (1) become transition operators
$p_{i}^{\dagger}p_{p}p_{p}=\ket{2}\bra{3}$, and the Hamiltonian for the
scattering process reduces to that of the three-level system driven by a field
corresponding to the condensate,

$H_{s}=\frac{\Omega}{2}\big{(}e^{-i\omega t}\ket{3}\bra{2}+e^{i\omega
t}\ket{2}\bra{3}\big{)}.$ (3)

This implies that the energy transfers between the three-level system and the
condensate are work, in accordance with the standard partition of energy
currents [30]

$\frac{d\langle
E\rangle}{dt}=\Tr[\rho\frac{dH}{dt}]+\Tr[\frac{d\rho}{dt}H]=\langle\dot{W}\rangle+\langle\dot{Q}\rangle,$
(4)

in which work arises from the time-dependence of a Hamiltonian and heat from
that of the density-matrix. Using this along with Eq. (3), the power supplied
to the condensate is

$\displaystyle\langle\dot{W}\rangle=\left\langle\frac{\partial H_{s}}{\partial
t}\right\rangle$ $\displaystyle=\omega\Omega\mathfrak{Im}\left\\{\Tr[\rho
e^{i\omega t}\ket{2}\bra{3}]\right\\}$
$\displaystyle=\omega\Omega\mathfrak{Im}\rho_{32},$ (5)

where $\rho_{32}$ is an element of the density matrix in the time-dependent
basis introduced below.

Additionally, the three-level system exchanges energy and particles with a
reservoir of acoustic phonons and a reservoir of excitons that is generated by
the pump. The phonons act as the cold bath of the heat engine, extracting heat
from the working medium as the excitons relax to the bottleneck
($\ket{2}\rightarrow\ket{1}$). The excitons generated by the external pump act
as the hot bath. To stay within the framework of a three-level model we
suppose they feed the states at the bottleneck, causing transitions from
$\ket{1}$ to $\ket{3}$. The lifetime of a reservoir exciton is relatively
long, so that the hot reservoir forms a quasi-equilibrium state with
temperature $T_{h}$ and chemical potential $\mu$. This means that the heat
engine operates with a chemical potential drop between the hot and cold
reservoirs. Since particles flow through the engine across this chemical
potential drop, its energy input consists of both work and heat.

In summary, the thermal machine operates in a mode in which the hot reservoir
inputs work and heat into the three-level system, which performs work by
interacting with the condensate and dumps the rest of the energy into a cold
bath of phonons. Combining the free Hamiltonian of the three-level system, the
coupling to the condensate, and the coupling to the baths, we have

$\displaystyle H_{0}+H_{sb}=$
$\displaystyle(\epsilon_{2}+\Delta)\ket{3_{R}}\bra{3_{R}}+\epsilon_{2}\ket{2_{R}}\bra{2_{R}}$
(6) $\displaystyle+$
$\displaystyle\epsilon_{1}\ket{1_{R}}\bra{1_{R}}+\frac{\Omega}{2}\big{(}\ket{3_{R}}\bra{2_{R}}+\ket{2_{R}}\bra{3_{R}}\big{)}$
$\displaystyle+$
$\displaystyle\sum_{k}g_{k}^{\text{ph}}\big{(}\ket{1_{R}}\bra{2_{R}}+\ket{2_{R}}\bra{1_{R}}\big{)}(b_{k}+b_{k}^{\dagger})$
$\displaystyle+$
$\displaystyle\sum_{q}g_{q}^{\text{p}}\big{(}\ket{3_{R}}\bra{1_{R}}x_{q}+\ket{1_{R}}\bra{3_{R}}x_{q}^{\dagger}\big{)},$

where the time dependence in Eq. (3) has been eliminated by using the rotating
frame $(\ket{1_{R}},\ket{2_{R}},\ket{3_{R}})=(\ket{1},\ket{2},e^{-i\omega
t}\ket{3})$. This would produce a time dependence in the coupling to the hot
bath, which, however, has been removed by using an interaction picture with
respect to a free Hamiltonian $\sum_{q}\omega x_{q}^{\dagger}x_{q}$ and making
the rotating wave approximation. The detuning
$\Delta=(\epsilon_{3}-\epsilon_{2})-\omega$ corresponds to the energy
difference between the low-energy final state of the scattering and the
condensate. We define the zero of energy to be $\epsilon_{1}=0$. For the other
energies we use values representative of a GaAs microcavity, $\epsilon_{2}=5$
meV and $\epsilon_{3}-\epsilon_{2}=1$ eV.

The third and fourth terms of Eq. (6) represent the interactions with the cold
and hot baths, with coupling strengths $g_{k}^{\text{ph}}$, and
$g_{q}^{\text{p}}$, and annihilation operators $b_{k}$ and $x_{q}$,
respectively. In the present case, they refer to phonons, and excitons in
states which repopulate the bottleneck. We note that one might expect the pump
to repopulate the high-energy states which would then relax to the bottleneck
by phonon emission. This could be incorporated by using a four-level model,
with an additional state $\ket{4}$ which is populated by the pump and relaxes
to $\ket{3}$ by phonon emission. Such a model could be further extended to
include higher phonon states, giving a form similar to that used for photon
condensation by Kirton and Keeling [4].

##  III Methods

To analyze Eq. (6) we use the standard method, previously applied to a three-
level heat engine in [10], in which one transforms to the eigenbasis of
$H_{0}$ and eliminates the heat baths within the Born-Markov approximation.
The transformation to the eigenbasis of $H_{0}$ is effected by the rotation

$\begin{pmatrix}\ket{3_{R}^{\prime}}\\\ \ket{2_{R}^{\prime}}\\\
\ket{1_{R}^{\prime}}\end{pmatrix}=\begin{pmatrix}\cos\theta/2&\sin\theta/2&0\\\
-\sin\theta/2&\cos\theta/2&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}\ket{3_{R}}\\\ \ket{2_{R}}\\\
\ket{1_{R}}\end{pmatrix}=U\begin{pmatrix}\ket{3_{R}}\\\ \ket{2_{R}}\\\
\ket{1_{R}}\end{pmatrix},$ (7)

which makes $H_{0}^{\prime}=U^{\dagger}H_{0}U$ diagonal when
$\tan\theta=\Omega/\Delta$. The resulting energies are $e_{1}=\epsilon_{1}=0$
and
$e_{2,3}=\big{(}\epsilon_{3}-\omega+\epsilon_{2}\mp\sqrt{\Delta^{2}+\Omega^{2}}\big{)}/2$.
We follow common practice and neglect the principal value terms, but do not
make the secular approximation. The resulting equations-of-motion are given in
Eqs. (14–16) and Eqs. (24–29).

This method can be extended, by introducing counting fields [31,
gasparinetti_heat-exchange_2014], to obtain the energy currents to the bath.
We find

$\displaystyle\langle\dot{Q}_{\text{c}}\rangle$
$\displaystyle=e_{3}R^{c}_{3}+e_{2}R^{c}_{2},$ (8)
$\displaystyle\langle\dot{E}_{\text{h}}\rangle$
$\displaystyle=(e_{3}+\omega)R^{h}_{3}+(e_{2}+\omega)R^{h}_{2}.$ (9)

The energy current to the cold bath is entirely heat and hence denoted
$Q_{c}$, whereas that to the hot bath comprises both heat and work, with the
latter arising from the flow of particles from the exciton reservoir, which in
general has a non-zero chemical potential. The rates which appear in the
energy currents are

$\displaystyle R^{c}_{2,3}$ $\displaystyle=\pi
J_{c}(e_{2,3})[\mp(n_{c}(e_{2,3})+1)\sin\theta\mathfrak{Re}[\rho_{23}]$
$\displaystyle+(1\pm\cos\theta)((n_{c}(e_{2,3})+1)\rho_{22}-n_{c}(e_{2,3})\rho_{11})]$
(10) $\displaystyle R^{h}_{2,3}=\pi J_{h}(e_{2,3}+$
$\displaystyle\omega)[\mp(n_{h}(e_{2,3}+\omega)+1)\sin\theta\mathfrak{Re}[\rho_{23}]$
$\displaystyle+(1\mp\cos\theta)((n_{h}(e_{2,3}+\omega)+1)\rho_{33}$
$\displaystyle\qquad\qquad\qquad-n_{h}(e_{2,3}+\omega)\rho_{11})].$ (11)

Here $\rho_{ij}$ are the elements of the density matrix in the basis
$(\ket{1_{R}},\ket{2_{R}},\ket{3_{R}})$. These rates can also be expressed in
terms of the elements of the density matrix in the diagonal basis,
$(\ket{1_{R}^{\prime}},\ket{2_{R}^{\prime}},\ket{3_{R}^{\prime}})$. Denoting
those elements by $\rho_{ij}^{\prime}$ the rates are

$\displaystyle R^{c}_{2,3}=\pi J_{c}(e_{2,3})\\{$
$\displaystyle(1\pm\cos\theta)[(n_{c}(e_{2,3})+1)\rho^{\prime}_{22,33}-n_{c}(e_{2,3})\rho^{\prime}_{11}]+(n_{c}(e_{2,3})+1)\sin\theta\mathfrak{Re}[\rho^{\prime}_{23}]\\}$
(12) $\displaystyle R^{h}_{2,3}=\pi J_{h}(e_{2,3}+\omega)\\{$
$\displaystyle(1\mp\cos\theta)[(n_{h}(e_{2,3}+\omega)+1)\rho^{\prime}_{22,{33}}-n_{h}(e_{2,3}+\omega)\rho^{\prime}_{11}]-(n_{h}(e_{2,3}+\omega)+1)\sin\theta\mathfrak{Re}[\rho^{\prime}_{23}]\\}.$
(13)

$J_{\text{c,h}}$ are the spectral densities for the cold and hot baths,
respectively. In the following we take a Lorentzian for the cold phonon bath,
$J_{\text{c}}(x)=(\alpha^{\text{c}}/2)\left[(x-\varepsilon_{2})^{2}+g_{\text{ph}}^{2}\right]^{-1}$,
and a step function for the hot pump bath
$J_{\text{h}}(x)=\alpha^{\text{h}}\Theta(x-E_{0})$. The parameters
$\alpha^{\text{c,h}}$ control the maximum values of the spectral density, and
$g_{\text{ph}}$ the width of the Lorentzian distribution, whose value needs to
be large enough to allow the Born-Markov approximation. $E_{0}$ is the energy
of the lowest exciton state in the pump reservoir. The bath occupation
function for the cold bath is
$n_{\text{c}}(E)=(e^{\beta_{\text{c}}E}-1)^{-1}$, and that for the hot bath is
the corresponding grand-canonical form, including the chemical potential,
$n_{\text{h}}(E)=(e^{\beta_{\text{h}}(E-\mu)}-1)^{-1}$.

##  IV Results

We begin by analyzing the dynamics of the cycle in terms of the occupations of
the states $\ket{n_{R}}$ ($n=1,2,3$), which we refer to as bare states. The
diagonal elements of the density matrix, in the bare-state basis, obey

$\displaystyle(\dot{\rho})_{11}=\left(R_{2}^{c}+R_{3}^{c}\right)+\left(R_{2}^{h}+R_{3}^{h}\right)\equiv
R^{c}-R^{h}$ (14)
$\displaystyle(\dot{\rho})_{22}=-\left(R_{2}^{c}+R_{3}^{c}\right)+R=-R^{c}+R$
(15)
$\displaystyle(\dot{\rho})_{33}=-\left(R_{2}^{h}+R_{3}^{h}\right)-R=R^{h}-R$
(16)

Here

$R=\langle\dot{W}\rangle/\omega=\Omega\mathfrak{Im}\rho_{32},$ (17)

and $R_{2,3}^{c,h}$ are given by Eqs. (10) and (11) or (12) and (13). From
these expressions we identify $R^{c}=(R_{2}^{c}+R_{3}^{c})$ as the rate of
population transfer from $\ket{2_{R}}$ to $\ket{1_{R}}$, due to the cold bath,
$R^{h}=-(R_{2}^{h}+R_{3}^{h})$ at that from $\ket{1_{R}}$ to $\ket{3_{R}}$,
due to the hot bath, and $R$ as that from $\ket{3_{R}}$ to $\ket{2_{R}}$, due
to the interaction with the condensate. The steady-state condition,

$R=R^{c}=R^{h},$ (18)

is that the rates around each part of the cycle are equal.

The energy fluxes to the baths, Eqs. (8,9), can be interpreted [11] in terms
of the eigenstates of $H_{0}$, $\ket{n_{R}^{\prime}}$, which differ from the
bare states, $\ket{n_{R}}$, around which the population circulates. Due to the
oscillating driving field each of the former is a Floquet state, or dressed
state, associated with a periodic quasi-energy. Fig. 2 illustrates the
formation of these states in the working medium. The bare spectrum, in the
absence of driving, becomes replicated, and the resulting level mixing
produces energy shifts.

Figure 2: Energy diagram of the three-level system with two Floquet zones
included in the diagonal ($e_{i}$) and bare basis ($\varepsilon_{i}$). Because
of the mixing generated by the condensate, the cold bath extracts from the
system energy contributions equivalent to the energies of the second and third
eigenvalues ($e_{2},e_{3}$). The hot bath acts in a similar way, but in the
second Floquet space, extracting $e_{2}+\omega$ and $e_{3}+\omega$ of energy
from the system. The power output is given by transitions from one zone to the
other, emitting $\omega$ of energy. The energy gap between the two excited
states of the diagonal Hamiltonian is defined as
$\Lambda=\sqrt{\Delta^{2}+\Omega^{2}}$.

The heat current to the cold bath, Eq. (8) is the sum of two contributions,
each the product of the quasi-energy of a dressed state $e_{2},e_{3}$, and one
of the two rates, $R^{c}_{2,3}$. Thus, we can interpret this heat flow as
arising from two transitions in the Floquet spectrum, as illustrated in Fig.
2. The hot bath acts similarly, but induces transitions from the ground state
to levels in the second quasi-energy band. The cycle is closed by the power
output, Eq. (5), that consists of quanta of energy $\omega$ emitted at the
rate $R$. This is illustrated in Fig. 2, as two independent transitions, each
between equivalent dressed states, suggesting the formation of two independent
cycles [11]. However, it should be noted that this interpretation implies a
secular approximation and only holds in a strong-driving limit. The energy
currents cannot in general be interpreted as coming from two separate cycles,
as $R^{c}_{2,3}\neq-R^{h}_{2,3}$. In fact, the rates in Eqs. (12,13) contain
contributions of the form of Boltzmann transition rates, and interference
terms weighted by the mixing angle $\theta$. The latter are responsible for
coupling the two fictitious cycles, that is, because of the sign difference
they constitute transitions to the ground state in inverted directions, which
is effectively an inter-cycle transition.

The first law of thermodynamics is
$\langle\dot{W}\rangle+\langle\dot{Q}_{c}\rangle+\langle\dot{E}_{h}\rangle=0$.
Using the forms for $\langle\dot{W}\rangle,\langle\dot{Q}_{c}\rangle$, and
$\langle\dot{E}_{h}\rangle$, and the steady-state condition, Eq. (18), it
becomes

$e_{2}(R^{c}_{2}+R^{h}_{2})+e_{3}(R^{c}_{3}+R^{h}_{3})=0.$ (19)

This condition is not exactly satisfied, but holds to a good approximation in
the parameter regimes used here, where the Born-Markov approximation is
appropriate.

We can describe a complete system, including the condensate, by noting that
the growth of the condensate comes from the output of a large number $M$ of
identical three-level heat engines. $M$ represents the number of states in the
bottleneck region that generate scattering to the ground state of the
polariton distribution, and is on the order of $M=10^{4}$ in typical cases.
The growth rate of a condensate, of population $N$, is then $MR(N)$, where the
dependence of the cycle rate $R$ on $\Omega$ and therefore $N$ has been made
explicit. This growth competes against the losses due to radiative polariton
decay, with rate $\gamma$. The steady-state condition for the condensate
number, $N_{c}$, is then that

$MR(N_{c})=N_{c}\gamma.$ (20)

In the following, we choose parameter values typical of GaAs microcavities.
The prefactor of the phonon spectral density is chosen in line with the
exciton relaxation times obtained numerically and experimentally [16]
$\alpha^{\text{c}}=0.1$ ps-1. We take the linewidth of the phonons to be
$g_{\text{ph}}=1.7$ ps-1. The temperature of the cold bath is the lattice
temperature of the semiconductor, which we take to be $10$ K. The spectral
density, as well as the chemical potential, and temperature of the hot bath,
depend non-trivially on the incoherent pump. We choose $\alpha^{\text{h}}=0.2$
ps-1, $E_{0}=1$ eV, $\mu-E_{0}=8$ meV and $T_{\text{h}}=200$ K initially, and
investigate the effects of varying these parameters. We assume a polariton
lifetime of $1/\gamma=1$ ps. For the polariton-polariton scattering strength
we use the expression derived in [21] to estimate $g_{c}\approx 0.048$ ps-1.
The detuning appears in our model with fixed condensate and transition
energies, and will determine the condensate in-scattering rate, i.e., the gain
from the working medium. In practice, however, the condensate, pump, and idler
states lie in a continuum, and a range of detunings are present. We expect
condensation to occur in the most favorable mode, and so focus on the resonant
case $\Delta=0$.

Figure 3: Phase diagram of polariton condensation as a function of the
temperatures of the hot and cold reservoir for different microcavity decay
rates. The points below each line correspond to the condensate region. The red
line corresponds to the case of infinite microcavity lifetime, where
condensation occurs at the onset of inversion. The parameter region of the
condensed phase is smaller for larger radiative decay, as it implies more
inversion is required to produce condensation.

Eq. (20) and the equations-of-motion can be solved numerically to determine
the condensate occupation $N_{c}$ in the steady-state. This solution can be
used to map the phase diagram, i.e., the region of parameter space in which
$N_{c}\neq 0$. Fig. 3 shows a numerically computed phase diagram for
condensation, in terms of the temperatures of the hot (pump) and cold
(lattice) heat baths, for several values of the decay rate $\gamma$. The
condensed state, with $N_{c}\neq 0$, lies below each curve, with the normal
state, $N_{c}=0$, above it. For these parameters the transition is continuous.
It can be seen that the phase boundaries approach the red diagonal line as
$\gamma\rightarrow 0$. This corresponds to the ideal thermodynamically-
reversible limit [7], in which condensation occurs when inversion is reached,
$\rho_{33}/\rho_{22}=1$. As noted by Scovil and Schulz-DuBois [7], the
requirement of inversion implies laser action is bounded by the Carnot
efficiency, and reaches it in the reversible limit where inversion first
appears. Generalizing their argument to include the chemical potential of the
hot reservoir, the condition for inversion

$\frac{\rho_{33}}{\rho_{22}}=\frac{\rho_{33}}{\rho_{11}}\frac{\rho_{11}}{\rho_{22}}=e^{-\beta_{h}(\varepsilon_{3}-\mu)}e^{\beta_{c}\varepsilon_{2}}>1,$
(21)

implies that the efficiency

$\displaystyle\eta=(W_{out}-W_{in})/Q_{in}=(\omega-\mu)/(\epsilon_{3}-\mu),$
(22)

is less than the Carnot efficiency $\eta_{C}=1-T_{c}/T_{h}$. Here we have used
$\varepsilon_{3}-\varepsilon_{2}=\omega$ and identified $\omega$ as the work
output per cycle, $\mu$ as the work input, and $\varepsilon_{3}-\mu$ as the
heat input (with $\varepsilon_{1}=0$).

The continuous transition can be understood using a perturbative expansion of
the steady-state in $\Omega$. In the normal state we have $\Omega=0$ so
$\theta\rightarrow 0$, $e_{2}=\varepsilon_{2},e_{3}=\varepsilon_{3}-\omega$.
The equations-of-motion reduce to the standard Lamb equations for a three-
level laser [11], given in the Appendix. The steady-state has the populations
$\rho_{33}/\rho_{11}$ and $\rho_{22}/\rho_{11}$ in equilibrium with the hot
and cold baths, respectively, and vanishing coherences. Expanding around this
solution gives the steady-state coherence to first order in $\Omega$, and
hence, using Eq. (17), the scattering rate

$R=\frac{\Omega^{2}\left(\rho_{33}-\rho_{22}\right)}{\gamma^{c}_{\downarrow}+\gamma^{h}_{\downarrow}}=\frac{4g_{c}^{2}N\left(\rho_{33}-\rho_{22}\right)}{\gamma^{c}_{\downarrow}+\gamma^{h}_{\downarrow}}.$
(23)

Here $\gamma^{c}_{\downarrow}$ ($\gamma^{h}_{\downarrow}$) are the emission
rates into the cold (hot) bath, respectively. Eq. (23) is the Fermi golden
rule expression for scattering into a final state, with a linewidth generated
by the emission into the baths. It defines, via Eq. (20), the critical
inversion at which gain exceeds loss, and hence the phase boundary. The
energies involved in the cycle are unaffected by the condensate to this order,
so the efficiency, Eq. (22), is constant.

The red line in Fig. 3 separates the diagram into two regions: above it, the
efficiency set by the energies is greater than the Carnot efficiency, and
condensation can never occur. It emerges immediately below this line for
$\gamma\rightarrow 0$, where condensation is supported by an infinitesimal
power flow from the working medium. For non-zero $\gamma$, condensation
requires finite, and therefore irreversible, power flows from the working
medium, which will produce a lower-than-Carnot efficiency for given bath
temperatures. The phase boundary therefore departs from the reversible line;
the non-zero $\gamma$ implies that the temperature difference $T_{h}-T_{c}$
required to drive condensation is increased, as it must overcome the loss. In
more practical terms, a larger $\gamma$ requires a larger inversion to
overcome the loss, and hence a higher (lower) temperature for the hot (cold)
bath.

Figure 4: Size of the condensate as a function of the temperature of the cold
(a) and hot (b) baths for various temperatures of the hot bath, and chemical
potentials, respectively. A second order phase transition is observed from the
normal state with $N_{c}=0$ to a condensed state with non-vanishing
occupation.

To determine the steady-state occupation in the condensed state we solve Eq.
(20) using the steady-state calculated numerically. This incorporates the
dependence of the dressed-state energies and wavefunctions on the condensate
occupation, which provides the nonlinearities, beyond Eq. (23), that are
needed to stabilize the condensate at a finite density. Results are shown in
Fig. 4, as functions of the two bath temperatures.

Considering first the dependence on $T_{h}$, we see that the condensate size
increases rapidly with $T_{h}$ once the threshold is crossed, and then
saturates or decreases. This is expected from the occupation of the hot bath
at the energy of the upper level ($\approx\varepsilon_{3}$). If we consider
the equilibrium of the upper level and the hot bath only, its population will
be $(1+e^{-\beta_{h}(\varepsilon_{3}-\mu)})^{-1}$, which is negligible until
$kT_{h}\sim(\varepsilon_{3}-\mu)$, at which point it rapidly grows before
saturating at $0.5$. The small decrease could be explained by noting that
$\gamma_{h}^{\downarrow}$, approximately proportional to
$(1+n_{h}(\varepsilon_{3}))$, increases with temperature due to the stimulated
emission into the hot reservoir. This will broaden the line and hence, as
expected from Eq. (23), reduce the gain. Such a mechanism has been shown to
produce an upper critical threshold, i.e., a maximum $T_{h}$, in the three-
level laser model [33].

The dependence of the condensate size on $T_{c}$ can be understood analogously
in terms of the thermal occupation of the cold bath, i.e., the phonons. Below
the critical temperature the condensate occupation increases smoothly, and
plateaus at very low temperatures. The critical $T_{c}$ arises here from the
requirements that the lower state of the $3-2$ transition should be
depopulated sufficiently to generate inversion – in other words, the phonons
must cool the high-momentum excitons, here represented by the pump and the
idler states, to a low temperature so that the idler population is small. The
saturation at low temperatures arises from the assumption of a single energy
gap for the phonon relaxation, which implies that the equilibrium idler
population is effectively zero for temperatures $kT_{c}\ll\varepsilon_{2}$.

Figure 5: Total in-scattering rate into the condensate as a function of its
size, with a Lorentzian (red) and constant (blue) phonon spectral density. The
rate reaches a maximum value for both spectral densities, however the
Lorentzian shape of the spectral density causes a reduction in $R$ for large
condensates.

Fig. 5 shows the total scattering rate into the condensate, $MR(N)$, in the
steady-state, as function of the condensate size $N$. This shows the existence
of a maximum power output of the three-level system. Such maxima are general
features of the three-level heat engine [11], caused by the energy-dependence
of the scattering rates, which are sampled at the dressed-state energies
$e_{2,3}$ with splitting $\sim\Omega$ when $\Omega\gg\Delta$. In Ref. 11 the
existence of a maximum power was attributed to the frequency-dependence of the
hot bath occupation function. Although this thermal factor will produce such
an effect in general, it cannot explain the particular peak in Fig. 5, since
here $kT_{h}$ is much larger than the energy shifts at the maximum power
point. It arises, instead, from the spectral function of the cold (phonon)
environment, which we have modeled as a Lorentzian of width $g_{ph}$. When
$N=0$ the dressed-states energies lie at the peak of this Lorentzian, but they
move away as $\Omega\sim\sqrt{N}$ increases, reducing the emission rate of
phonons into the environment and slowing the cycle. This interpretation is
consistent with the result for a flat phonon spectral density, shown for
comparison. Although the assumption of a Lorentzian spectral density for the
phonons is not necessarily realistic, the heat baths relevant for polariton
condensates do have spectral structure. This can be expected to lead to a
maximum power once the condensate nonlinearities become comparable to either
the scale of the structure or the temperatures.

Figure 6: In-scattering rates with $\Delta=0.3$ meV and $E_{0}-\epsilon_{3}$=0
(black solid), 0.33 meV (blue dotted), -0.2 meV (orange dotted). The red curve
is the loss rate. These parameters lead to in-scattering rates with a concave
region (black solid and orange dotted) and discontinuities (blue dotted and
orange dotted). Such forms give rise to first-order transitions and
bistability (see text).

In some cases there are additional phenomena which can arise from the non-
trivial spectral densities of the baths. Fig. 6 shows the gain as a function
of condensate size for three different sets of parameters, revealing two
unexpected effects. The first is the presence of discontinuities in the
steady-state, either connecting a region of zero gain to non-zero gain, or
causing a sudden drop in the gain. These arise from the interplay between the
energy shifts of the driving and the spectral density of the hot bath, which
we have assumed to have a step-function onset at energy $E_{0}$. The dashed-
blue curve arises in a situation where the upper level of the heat engine lies
just below the hot bath, $\varepsilon_{3}<E_{0}$, so that there is no
population of the upper level in the absence of the condensate, and hence no
gain. However, as the condensate size increases the upper dressed state,
$e_{3}$, moves to higher energies, and crosses $E_{0}$, at which point gain
appears. A related situation, shown by the orange dotted curve, appears if
$\epsilon_{3}$ starts just above $E_{0}$. In that case, while the two dressed
states start in the band, as the condensate occupation increases the lower-
energy one drops below $E_{0}$, and the gain suddenly decreases. The second,
more subtle, effect is also visible in this curve, as well as the black curve,
which corresponds to the case where $\epsilon_{3}$ starts exactly at $E_{0}$.
This is the presence of a concave part of the gain curve at small $N$. We
suggest this is because there is also a small positive detuning here, so that
the upper dressed state is composed mostly of $|3\rangle$ in the limit
$N\rightarrow 0$. However, as the splitting $\Omega$ reaches the detuning
$\Delta$, this component of the upper dressed state – which is the only one
pumped by the hot bath at this point in the curve – begins to reduce
significantly, suppressing the gain. We note this concave form of $R(N)$
implies the transition is first-order, with the crossing point solving
$MR(N)=\gamma N$ first appearing at a non-zero $N$. A first-order behavior can
also be expected from the physics of the step-like curve which, for the more
realistic case of a smoothly-increasing spectral density, will become a smooth
but concave $R(N)$. In addition, we expect these forms to lead to bistability
of the condensed states, since they have two intersections with the loss
curve.

Figure 7 shows some results for the power output and efficiency of the
condensate. The white regions denote the non-condensed phase, and the colors
depict the net power $P_{\text{out}}-P_{\text{in}}$, and net efficiency
$-P_{\text{out}}/Q_{\text{h}}$, in the steady-state. The two plots in the left
column show the effects of varying the hot bath temperature and chemical
potential, while the right column shows the corresponding effects of varying
the coupling strengths. The net power increases with increasing hot bath
temperature, in line with the condensate occupation shown in Fig. 4. However,
while the condensate occupation increases with increasing chemical potential,
the net power output decreases. This is because although increasing $\mu$
increases the cycle rate and hence the gross power, $\omega R=\omega\gamma
N/M$, it also increases the work input, giving an overall decrease in the net
power, $(\omega-\mu)R$. Considered as functions of the coupling strengths, the
net power output increases significantly with an increasing coupling to the
cold bath. However, the power output depends only weakly on the hot bath
coupling parameter, and can even decrease as it increases. We suggest that the
weak dependence arises because, for the temperatures chosen, the cycle rate is
limited by emission into the cold bath rather than absorption from the hot
one. This emphasizes the need for rapid phonon thermalization in the high-
momentum excitons, in order to drive condensation (see Fig. 1).

The efficiency, as a function of temperature and coupling constants, is shown
in the lower two panels of Fig. 7. It is shown relative to the Carnot
efficiency, which explains most of the dependence on the temperatures; the
energy shifts are small for these parameters, so that $\eta$ is well
approximated by Eq. (22). The efficiency at threshold is below the Carnot
efficiency due to the finite loss rate. Since the efficiency is determined
almost entirely by the energy levels, there is only a very small effect of the
coupling strengths.

Figure 7: Net output power (top row) and efficiency (bottom row) of a steady-
state condensate, as functions of the temperature and chemical potential of
the hot bath (left column), and of the coupling strengths to the two baths
(right column). The chemical potential is measured relative to the lowest
state of the pump reservoir, $E_{0}$, and the efficiencies are normalized by
the Carnot efficiency. The net output power grows if $T_{h}$ or $(E_{0}-\mu)$
increase and has stronger dependence on the coupling strength with the cold
than with the hot bath. The efficiency shows only a weak dependence on the
parameters considered.

##  V Conclusions

In summary, we have argued that a driven-dissipative condensate is a form of
heat engine, and as such requires consideration of both a hot and cold bath.
We considered in depth the case of a microcavity polariton condensate, and
constructed a minimal heat-engine model that captures the key processes
involved in condensation. Our results show how that condensation is determined
by the temperatures of both the hot (exciton) and cold (phonon) bath, and that
the temperature difference between these must exceed that required in the
reversible limit in order to overcome the polariton loss. They also emphasize
the importance of rapid cooling in the high-momentum exciton states, in order
to maintain effective scattering into the condensate. The maximum power output
of the condensate shows effects of the spectral densities of the environment,
and the formation of dressed states, which can also produce unusual phenomena
including a first-order phase boundary and bistability. Our results provide
guidance for extending the regimes and systems which support condensation, and
our methods could be extended to consider other examples such as the
condensation [3] and thermalization [34] of photons. It would also be
interesting to explore the reversed operation of the heat-engine and the
possibilities for phonon refrigeration [35].

###### Acknowledgements.

P.R.E. acknowledges support from Science Foundation Ireland (21-FFP-10142).

## Authors declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Luisa Toledo Tude: Conceptualization (supporting), formal analysis (lead),
methodology (lead), software (lead), visualization (lead), writing – original
draft preparation (lead), writing – review and editing (equal). Paul Eastham:
Conceptualization (lead), supervision (lead), funding acquisition (lead),
methodology (supporting), supervision (lead), writing – original draft
preparation (supporting), writing – review and editing (equal).

## Data availability

The code used to generate the data that support the findings of this study is
openly available in Zenodo at http://doi.org/. In order to meet institutional
and research funder open access requirements, any accepted manuscript arising
shall be open access under a Creative Commons Attribution (CC BY) reuse
license with zero embargo.

##  Appendix A Equations-of-motion

The equations-of-motion for the diagonal elements of the density matrix, in
the bare basis, are given in Eqs. (14–16). We give here the remaining terms,
considering the contributions from the two baths separately. We use an
abbreviated notation in which subscripts denote the frequencies at which the
spectral densities and bath occupations are sampled, so that in expressions
for the cold (hot) bath we have $J_{i}=J_{c}(e_{i}),n_{i}=n_{c}(e_{i})$
($J_{i}=J_{h}(e_{i}+\omega),n_{i}=n_{h}(e_{i}+\omega)$). For the cold bath we
have

$\displaystyle\dot{\rho}_{23}|_{c}=-\frac{\pi}{2}$
$\displaystyle\\{[\rho_{11}(J_{2}n_{2}-J_{3}n_{3})+\rho_{33}(J_{3}(1+n_{3})-J_{2}(1+n_{2}))]\sin\theta$
$\displaystyle+\rho_{23}[J_{2}(1+n_{2})(1+\cos\theta)+J_{3}(1+n_{3})(1-\cos\theta)]\\},$
(24) $\displaystyle\dot{\rho}_{12}|_{c}=-\frac{\pi}{2}\\{($
$\displaystyle\rho_{12}-\rho_{21})[J_{2}(1+2n_{2})(1+\cos\theta)+J_{3}(1+2n_{3})(1-\cos\theta)]$
$\displaystyle+$
$\displaystyle(\rho_{13}-\rho_{31})[J_{3}(1+n_{3})-J_{2}(1+n_{2})]\sin\theta\\},$
(25) $\displaystyle\dot{\rho}_{13}|_{c}=-\frac{\pi}{2}$
$\displaystyle\\{\rho_{13}[J_{2}n_{2}(1+\cos\theta)+J_{3}n_{3}(1-\cos\theta)]$
$\displaystyle+\rho_{21}[J_{2}n_{2}-J_{3}n_{3}]\sin\theta\\}.$ (26)

The corresponding expressions for the hot bath are

$\displaystyle\dot{\rho}_{23}|_{h}=-\frac{\pi}{2}$
$\displaystyle\\{[\rho_{11}(J_{2}n_{2}-J_{3}n_{3})+\rho_{22}(J_{3}(1+n_{3})-J_{2}(1+n_{2}))]\sin\theta$
$\displaystyle+\rho_{23}[J_{2}(1+n_{2})(1-\cos\theta)+J_{3}(1+n_{3})(1+\cos\theta)]\\},$
(27)
$\displaystyle\dot{\rho}_{13}|_{h}=-\frac{\pi}{2}\\{\rho_{13}[J_{3}(1+2n_{3})(1+\cos\theta)+J_{2}(1+2n_{2})(1-\cos\theta)]$
$\displaystyle+\rho_{12}[J_{3}(1+n_{3})-J_{2}(1+n_{2})]\sin\theta\\},$ (28)
$\displaystyle\dot{\rho}_{12}|_{h}=-\frac{\pi}{2}$
$\displaystyle\\{\rho_{12}[J_{3}n_{3}(1+\cos\theta)+J_{2}n_{2}(1-\cos\theta)]\\}.$
(29)

In the limit $\theta\rightarrow 0$, the energies of the dressed states reduce
to the original energies,
$e_{2}=\varepsilon_{2},e_{3}=\varepsilon_{3}-\omega$, and the equations-of-
motion become the standard Lamb equations for a three-level laser [11]. The
population transfer rates become

$\displaystyle R^{c}$ $\displaystyle=2\pi J_{2}(1+n_{2})\rho_{22}-2\pi
J_{2}n_{2}\rho_{11}$
$\displaystyle=\gamma^{c}_{\downarrow}\rho_{22}-\gamma^{c}_{\uparrow}\rho_{11},$
(30) $\displaystyle-R^{h}$ $\displaystyle=2\pi J_{3}(1+n_{3})\rho_{33}-2\pi
J_{3}n_{3}\rho_{11}$
$\displaystyle=\gamma^{h}_{\downarrow}\rho_{33}-\gamma^{h}_{\uparrow}\rho_{11}.$
(31)

The dissipative contributions to the equations-of-motion for the off-diagonal
elements of the density matrix become

$\displaystyle\dot{\rho}_{23}=-\frac{1}{2}(\gamma^{c}_{\downarrow}+\gamma^{h}_{\downarrow})\rho_{23},$
(32)
$\displaystyle\dot{\rho}_{13}=-\frac{1}{2}(\gamma^{h}_{\uparrow}+\gamma^{h}_{\downarrow}+\gamma_{c}^{\uparrow})\rho_{13},$
(33)
$\displaystyle\dot{\rho}_{12}=-\frac{1}{2}(\gamma^{c}_{\uparrow}+\gamma^{c}_{\downarrow}+\gamma_{h}^{\uparrow})\rho_{12}.$
(34)

where a counterrotating term from the cold bath has been neglected in Eq.
(34). These expressions describe the decay of the coherences, with the rates
related in the expected way to those of the populations.

## References

## References

  * Keeling and Kéna-Cohen [2020] J. Keeling and S. Kéna-Cohen, “Bose–Einstein Condensation of Exciton-Polaritons in Organic Microcavities,” Annu. Rev. Phys. Chem. 71, 435–459 (2020).
  * Deng, Haug, and Yamamoto [2010] H. Deng, H. Haug, and Y. Yamamoto, “Exciton-polariton Bose-Einstein condensation,” Rev. Mod. Phys. 82, 1489–1537 (2010).
  * Klaers  _et al._ [2010] J. Klaers, J. Schmitt, F. Vewinger, and M. Weitz, “Bose–Einstein condensation of photons in an optical microcavity,” Nature 468, 545–548 (2010).
  * Kirton and Keeling [2013] P. Kirton and J. Keeling, “Nonequilibrium Model of Photon Condensation,” Phys. Rev. Lett. 111, 100404 (2013).
  * Hakala  _et al._ [2018] T. K. Hakala, A. J. Moilanen, A. I. Väkeväinen, R. Guo, J.-P. Martikainen, K. S. Daskalakis, H. T. Rekola, A. Julku, and P. Törmä, “Bose–Einstein condensation in a plasmonic lattice,” Nature Phys 14, 739–744 (2018).
  * Demokritov  _et al._ [2006] S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, “Bose–Einstein condensation of quasi-equilibrium magnons at room temperature under pumping,” Nature 443, 430–433 (2006).
  * Scovil and Schulz-DuBois [1959] H. E. D. Scovil and E. O. Schulz-DuBois, “Three-Level Masers as Heat Engines,” Phys. Rev. Lett. 2, 262–263 (1959).
  * Geusic, Schulz-DuBios, and Scovil [1967] J. E. Geusic, E. O. Schulz-DuBios, and H. E. D. Scovil, “Quantum Equivalent of the Carnot Cycle,” Phys. Rev. 156, 343–351 (1967).
  * Kosloff [1984] R. Kosloff, “A quantum mechanical open system as a model of a heat engine,” The Journal of Chemical Physics 80, 1625–1631 (1984).
  * Geva and Kosloff [1994] E. Geva and R. Kosloff, “Three-level quantum amplifier as a heat engine: A study in finite-time thermodynamics,” Phys. Rev. E 49, 3903–3918 (1994).
  * Geva and Kosloff [1996] E. Geva and R. Kosloff, “The quantum heat engine and heat pump: An irreversible thermodynamic analysis of the three-level amplifier,” The Journal of Chemical Physics 104, 7681–7699 (1996).
  * Mitchison [2019] M. T. Mitchison, “Quantum thermal absorption machines: Refrigerators, engines and clocks,” Contemporary Physics 60, 164–187 (2019), arxiv:1902.02672 .
  * Imamog¯lu  _et al._ [1996] A. Imamog¯lu, R. J. Ram, S. Pau, and Y. Yamamoto, “Nonequilibrium condensates and lasers without inversion: Exciton-polariton lasers,” Phys. Rev. A 53, 4250–4253 (1996).
  * Kasprzak  _et al._ [2006] J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, “Bose–Einstein condensation of exciton polaritons,” Nature 443, 409–414 (2006).
  * Balili  _et al._ [2007] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West, “Bose-Einstein Condensation of Microcavity Polaritons in a Trap,” Science 316, 1007–1010 (2007).
  * Kavokin and Malpuech [2003] A. Kavokin and G. Malpuech,  _Cavity Polaritons_ , 1st ed., Thin Films and Nanostructures No. 32 (Elsevier, Acad. Press, Amsterdam, 2003).
  * Kavokin [2017] A. Kavokin, ed.,  _Microcavities_ , second edition ed., Series on Semiconductor Science and Technology No. 16 (Oxford University Press, Oxford ; New York, NY, 2017).
  * Carusotto and Ciuti [2013] I. Carusotto and C. Ciuti, “Quantum fluids of light,” Rev. Mod. Phys. 85, 299–366 (2013).
  * Cao  _et al._ [2004] H. T. Cao, T. D. Doan, D. B. Tran Thoai, and H. Haug, “Condensation kinetics of cavity polaritons interacting with a thermal phonon bath,” Phys. Rev. B 69, 245325 (2004).
  * Doan  _et al._ [2005] T. D. Doan, H. T. Cao, D. B. Tran Thoai, and H. Haug, “Condensation kinetics of microcavity polaritons with scattering by phonons and polaritons,” Phys. Rev. B 72, 085301 (2005).
  * Tassone and Yamamoto [1999] F. Tassone and Y. Yamamoto, “Exciton-exciton scattering dynamics in a semiconductor microcavity and stimulated scattering into polaritons,” Phys. Rev. B 59, 10830–10842 (1999).
  * Kasprzak  _et al._ [2008] J. Kasprzak, D. D. Solnyshkov, R. André, L. S. Dang, and G. Malpuech, “Formation of an Exciton Polariton Condensate: Thermodynamic versus Kinetic Regimes,” Phys. Rev. Lett. 101, 146404 (2008).
  * Piermarocchi  _et al._ [1996] C. Piermarocchi, F. Tassone, V. Savona, A. Quattropani, and P. Schwendimann, “Nonequilibrium dynamics of free quantum-well excitons in time-resolved photoluminescence,” Phys. Rev. B 53, 15834–15841 (1996).
  * Hartwell and Snoke [2010] V. E. Hartwell and D. W. Snoke, “Numerical simulations of the polariton kinetic energy distribution in GaAs quantum-well microcavity structures,” Phys. Rev. B 82, 075307 (2010).
  * Malpuech  _et al._ [2002] G. Malpuech, A. Kavokin, A. Di Carlo, and J. J. Baumberg, “Polariton lasing by exciton-electron scattering in semiconductor microcavities,” Phys. Rev. B 65, 153310 (2002).
  * Porras  _et al._ [2002] D. Porras, C. Ciuti, J. J. Baumberg, and C. Tejedor, “Polariton dynamics and Bose-Einstein condensation in semiconductor microcavities,” Phys. Rev. B 66, 085304 (2002), arxiv:cond-mat/0206276 .
  * Tassone  _et al._ [1997] F. Tassone, C. Piermarocchi, V. Savona, A. Quattropani, and P. Schwendimann, “Bottleneck effects in the relaxation and photoluminescence of microcavity polaritons,” Phys. Rev. B 56, 7554–7563 (1997).
  * Baumberg  _et al._ [2002] J. J. Baumberg, P. G. Savvidis, P. Lagoudakis, M. Martin, D. Whittaker, R. Butte, M. Skolnick, and J. Roberts, “Polariton traps in semiconductor microcavities,” Physica E: Low-dimensional Systems and Nanostructures 13, 385–389 (2002).
  * Linden, Popescu, and Skrzypczyk [2010] N. Linden, S. Popescu, and P. Skrzypczyk, “How Small Can Thermal Machines Be? The Smallest Possible Refrigerator,” Phys. Rev. Lett. 105, 130401 (2010).
  * Adesso  _et al._ [2018] G. Adesso, J. Anders, F. Binder, L. A. Correa, and C. Gogolin, eds.,  _Thermodynamics in the Quantum Regime: Fundamental Aspects and New Directions_, 1st ed., Fundamental Theories of Physics No. 195 (Springer International Publishing : Imprint: Springer, Cham, 2018).
  * Esposito, Harbola, and Mukamel [2009] M. Esposito, U. Harbola, and S. Mukamel, “Nonequilibrium fluctuations, fluctuation theorems, and counting statistics in quantum systems,” Rev. Mod. Phys. 81, 1665–1702 (2009).
  * Gasparinetti  _et al._ [2014] S. Gasparinetti, P. Solinas, A. Braggio, and M. Sassetti, “Heat-exchange statistics in driven open quantum systems,” New J. Phys. 16, 115001 (2014), arxiv:1404.3507 .
  * Li  _et al._ [2017] S.-W. Li, M. B. Kim, G. S. Agarwal, and M. O. Scully, “Quantum statistics of a single-atom Scovil-Schulz-DuBois heat engine,” Phys. Rev. A 96, 063806 (2017), arxiv:1710.00902 [quant-ph] .
  * Busley  _et al._ [2023] E. Busley, L. Espert Miranda, C. Kurtscheid, F. Wolf, F. Vewinger, J. Schmitt, and M. Weitz, “Sunlight-pumped two-dimensional thermalized photon gas,” Phys. Rev. A 107, 052204 (2023).
  * Klembt  _et al._ [2015] S. Klembt, E. Durupt, S. Datta, T. Klein, A. Baas, Y. Léger, C. Kruse, D. Hommel, A. Minguzzi, and M. Richard, “Exciton-Polariton Gas as a Nonequilibrium Coolant,” Phys. Rev. Lett. 114, 186403 (2015).