---
author:
- Makoto Yamaguchi and Tetsuo Ogawa
title: 'Equilibrium to nonequilibrium condensation in driven-dissipative semiconductor systems'
---

Introduction {#sec:intro}
============

In a semiconductor system, it is known that electron-hole (e-h) bound pairs can be formed by their Coulomb attraction when the conduction and valence band effectively reach an equilibrium state after the carriers are generated e.g. by laser excitation (Fig. \[RelaxationRedistribution\]). An exciton polariton is a quasi-bosonic particle composed of such a Coulomb-bound e-h pair (exciton) and a photon [@Weisbuch92; @Bloch98], the behaviors of which have attracted much attention due to their potential apprications through the Bose-Einstein condensation (BEC) [@Imamoglu96; @Deng02; @Kasprzak06], i.e. a macroscopic occupation of a single exciton-polariton state by a thermodynamic phase transition.

![ Excitation and thermalization process in a semiconductor. Electrons and holes generated by laser excitation subsequently undergo immediate intraband relaxations and redistributions in the conduction band (C.B.) and valence band (V.B.) to effectively reach an equilibrium state. Coulomb-bound e-h pairs (excitons) are formd when the equilibrium state is at sufficiently low temperature and low carrier density. []{data-label="RelaxationRedistribution"}](RelaxationRedistribution.eps){width="0.50\linewidth"}

A typical exciton-polariton system is shown in Fig. \[System\]. The system basically consists of semiconductor quantum wells (QWs) and a microcavity, the same structure as a vertical cavity surface emitting laser (VCSEL). In this context, a conventional lasing phase[^1] is involved in this system as well as the exciton-polariton BEC [@Snoke12]. At high densities, moreover, the Bardeen-Cooper-Schrieffer (BCS) -like ordered phase can potentially be caused where electrons and holes form the “Cooper pairs” [@Keldysh65; @Littlewood04], as is discussed in the BCS-BEC crossover in cold atom systems with Feshbach resonances [@Ohashi02; @Ohashi03]. These ordered phases are schematically shown in Fig. \[BEC-BCS-LASER\].

![ Schematic illustration of a typical exciton-polariton system. Exciton-polaritons are formed by the electrons and holes in the QWs and the photons confined between the two mirrors (microcavity). []{data-label="System"}](System.eps){width="0.50\linewidth"}

However, the BEC and BCS phases are in equilibrium, the situation of which is quite different from the semiconducotor laser in nonequilibrium. As a result, approaches for describing the BEC and BCS phases based on equilibrium statistical mechanics, e.g. the BCS theory [@Kamide10; @Byrnes10], are not applicable to the semiconductor laser because any nonequilibrium effects cannot be taken into account, such as pumping and loss. Conversely, past theories for describing the lasing operation, e.g. the Maxwell-Semiconductor-Bloch equations (MSBEs) [@Chow02; @Kamide11], cannot recover such equilibrium statistical approaches.[^2] The difficulty shown here has been one of problems to understand the underlying physics in exciton-polariton systems.

In such a situation, we have recently proposed a framework which can treat the phases of the BEC, BCS and laser in a unified way [@Yamaguchi12; @Yamaguchi13]. This framework is an extention of a nonequilibrium Green’s function approach developed in Ref. [@Szymanska06; @Szymanska07; @Keeling10] in which excitons are simply modeled by localized noninteracting two-level systems without internal e-h structures. Our formalism results in the BCS theory when the system can be regarded as in equilibrium, while it recovers the MSBE when nonequilibrium features become important. The internal e-h structures as well as the Coulomb interactions can also be taken into account within the mean-field approximation. In this contribution, we would like to give an introduction to such a “BEC-BCS-LASER crossover theory”.

![ Schematic illustration of several ordered phases involved in the exciton-polariton system. []{data-label="BEC-BCS-LASER"}](BEC-BCS-LASER.eps){width="0.90\linewidth"}

BCS theory and MSBE for exciton-polariton systems {#sec:BCS and MSBE}
=================================================

In exciton-polariton systems, the equilibrium phases (the BEC and BCS phases) can be described by the BCS theory while the nonequilibrium phase (the lasing phase) can be described by the MSBE. In this Section, we give an overview of the BCS theory and the MSBE to highlight their similarities and differences. For simplicity, we set $\hbar=k_{\mathrm{B}}=1$ in the followings.

Model {#subsec:Model}
-----

We first describe the Hamiltonian for the exciton-polariton system where electrons and holes in the QWs and photons in the microcavity are taken into account. The system Hamiltonian $\hat{H}_{\mathrm{S}}$ is then given by $\hat{H}_{\mathrm{S}} = \hat{H}_{0} + \hat{H}_{\mathrm{Coul}} + \hat{H}_{\mathrm{dip}}$. Here, $\hat{H}_{\mathrm{0}}$ is the Hamiltonian for free particles without interactions and written as $$\hat{H}_{\mathrm{0}} = \sum_{{\textbf{\itshape{k}}}} \left( {\epsilon_{\mathrm{e},{\textbf{\itshape{k}}}}}\oed_{{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}} + {\epsilon_{\mathrm{h},{\textbf{\itshape{k}}}}}\ohd_{{\textbf{\itshape{k}}}}\oh_{{\textbf{\itshape{k}}}} + {\epsilon_{\mathrm{ph},{\textbf{\itshape{k}}}}}\oad_{{\textbf{\itshape{k}}}}\oa_{{\textbf{\itshape{k}}}} \right),
\label{eq:HS}$$ where $\ope_{{\textbf{\itshape{k}}}}$, $\oh_{{\textbf{\itshape{k}}}}$, and $\oa_{{\textbf{\itshape{k}}}}$ are annihilation operators for electrons, holes, and photons with in-plane wave number ${\textbf{\itshape{k}}}$, respectively. ${\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}} = k^2/2m_{\mathrm{e(h)}} + E_{\mathrm{g}}/2$ is the energy dispersion of electrons (holes) with an effective mass $m_{\mathrm{e(h)}}$, while ${\epsilon_{\mathrm{ph},{\textbf{\itshape{k}}}}} = k^2/2m_{\mathrm{cav}} + E_{\mathrm{cav}}$ is that of photons with an effective mass $m_{\mathrm{cav}}$. $E_{\mathrm{g}}$ is the bandgap and $E_{\mathrm{cav}}$ is the energy of the cavity mode for ${\textbf{\itshape{k}}}=0$ [@Deng10Review].

In contrast, $\hat{H}_{\mathrm{Coul}}$ and $\hat{H}_{\mathrm{dip}}$ denote the Coulomb interaction and the light-matter interaction within the dipole approximation, respectively written as $$\begin{aligned}
\hat{H}_{\mathrm{Coul}} &=& \frac{1}{2}\sum_{{\textbf{\itshape{k}}},{\textbf{\itshape{k}}}',{\textbf{\itshape{q}}}}U'_{{\textbf{\itshape{q}}}}
    \left( \oed_{{\textbf{\itshape{k}}}+{\textbf{\itshape{q}}}}\oed_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{q}}}}\ope_{{\textbf{\itshape{k}}}'}\ope_{{\textbf{\itshape{k}}}} + ( \ope \leftrightarrow \oh ) 
    -2\oed_{{\textbf{\itshape{k}}}+{\textbf{\itshape{q}}}}\ohd_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{q}}}}\oh_{{\textbf{\itshape{k}}}'}\ope_{{\textbf{\itshape{k}}}} \right),
\label{eq:HCoul}\\
\hat{H}_{\mathrm{dip}} &=& -\sum_{{\textbf{\itshape{k}}},{\textbf{\itshape{q}}}} \left( g\oad_{{\textbf{\itshape{q}}}}\oh_{-{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}+{\textbf{\itshape{q}}}} + g^{*}\oa_{{\textbf{\itshape{q}}}}\oed_{{\textbf{\itshape{k}}}+{\textbf{\itshape{q}}}}\ohd_{-{\textbf{\itshape{k}}}} \right),
\label{eq:Hdip}\end{aligned}$$ where $U'_{{\textbf{\itshape{q}}}}=U'_{-{\textbf{\itshape{q}}}}$ and $U'_{{\textbf{\itshape{q}}}=0} \equiv 0$. Note that $[\hat{H}_{\mathrm{S}},\hat{N}_{\mathrm{S}}] = 0$ is satisfied when an excitation number of the system $\hat{N}_{\mathrm{S}}$ is defined as $\hat{N}_{\mathrm{S}} \equiv \sum_{{\textbf{\itshape{k}}}} [\oed_{{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}}/2 + \ohd_{{\textbf{\itshape{k}}}}\oh_{{\textbf{\itshape{k}}}}/2 + \oad_{{\textbf{\itshape{k}}}}\oa_{{\textbf{\itshape{k}}}}]$. For later convenience, therefore, we redefine $\hat{H}_{\mathrm{S}}-\mu\hat{N}_{\mathrm{S}}$ as $\hat{H}_{\mathrm{S}}$. This means that a grand canonical ensemble is assumed with a chemical potential $\mu$ if we are interested in equilibrium phases. In contrast, for time-dependent problems, this means that dynamics of physical quantities is captured on a rotating frame with the frequency $\mu$. Thus, $\mu$ is a given parameter identical to the chemical potential for the BEC and BCS phases (Subsection \[subsec:BCS\]), whereas it becomes a unknown variable equivalent to the lasing frequency for the semiconductor laser in a steady state (Subsection \[subsec:MSBE\]).

Mean-field approximation {#subsec:MF}
------------------------

The Hamiltonians shown in Subsection \[subsec:Model\] give a starting point for theories of the exciton-polariton system. However, in practice, it is difficult to exactly treat $\hat{H}_{\mathrm{Coul}}$ and $\hat{H}_{\mathrm{dip}}$ because these Hamiltonians cause many-body problems. In this Subsection, therefore, we discuss the mean-field (MF) approximation in order to reduce the problems to single-particle problems.

In general, the MF approximation is performed by writing a specific operator $\hat{O}$ as $\hat{O} = {\langle \hat{O} \rangle} + \mathnormal{\delta}\hat{O}$ and by neglecting quadratic terms with respect to $\mathnormal{\delta}\hat{O}$ in the Hamiltonians.[^3] Here, ${\langle \hat{O} \rangle} \equiv \mathrm{Tr}[\hat{O}\hat\rho]$ denotes the expectation value for the density operator $\hat\rho$ and the operator $\mathnormal{\delta}\hat{O}$ corresponds to a fluctuation around the expectation value. In our case, the interaction Hamiltonians of $\hat{H}_{\mathrm{Coul}}$ and $\hat{H}_{\mathrm{dip}}$ can easily be reduced to a single-particle problem by employing $\hat{O} \in \{ \oa_{{\textbf{\itshape{k}}}}, \oh_{-{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}'}, \oed_{{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}'}, \ohd_{{\textbf{\itshape{k}}}}\oh_{{\textbf{\itshape{k}}}'} \}$. As a result, with definitions of the photon field ${\langle \oa_{{\textbf{\itshape{k}}}} \rangle} \equiv \delta_{{\textbf{\itshape{k}}},0}\a0$, the polarization function ${\langle \oh_{-{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}'} \rangle} \equiv \delta_{{\textbf{\itshape{k}}},{\textbf{\itshape{k}}}'}\pk$, and the distribution functions of electrons ${\langle \oed_{{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}'} \rangle} \equiv \delta_{{\textbf{\itshape{k}}},{\textbf{\itshape{k}}}'}\nek$ and holes ${\langle \ohd_{{\textbf{\itshape{k}}}}\oh_{{\textbf{\itshape{k}}}'} \rangle} \equiv \delta_{{\textbf{\itshape{k}}},{\textbf{\itshape{k}}}'}\nhk$, the mean-field Hamiltonian $\hat{H}_{\mathrm{S}}^{\mathrm{MF}}$ is obtained as $$\begin{aligned}
\hat{H}_{\mathrm{S}}^{\mathrm{MF}} &=& \sum_{{\textbf{\itshape{k}}}} \left( {\tilde{\xi}_{\mathrm{e},{\textbf{\itshape{k}}}}}\oed_{{\textbf{\itshape{k}}}}\ope_{{\textbf{\itshape{k}}}} + {\tilde{\xi}_{\mathrm{h},{\textbf{\itshape{k}}}}}\ohd_{{\textbf{\itshape{k}}}}\oh_{{\textbf{\itshape{k}}}} - [{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}\oed_{{\textbf{\itshape{k}}}}\ohd_{-{\textbf{\itshape{k}}}} + \mathrm{H.c.}] \right) \nonumber\\
&&+\sum_{{\textbf{\itshape{k}}}} \left( {\xi_{\mathrm{ph},{\textbf{\itshape{k}}}}}\oad_{{\textbf{\itshape{k}}}}\oa_{{\textbf{\itshape{k}}}} - [g\pk\oad_0 + g^{*}\pk^{*}\oa_0] \right).
\label{eq:HMF}\end{aligned}$$ Here, constants are ignored because the following discussion is not affected. ${\tilde{\xi}_{\mathrm{e(h)},{\textbf{\itshape{k}}}}}$ and ${\xi_{\mathrm{ph},{\textbf{\itshape{k}}}}}$ are respectively defined as ${\tilde{\xi}_{\mathrm{e(h)},{\textbf{\itshape{k}}}}} \equiv {\tilde\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}} - \mu/2$ and ${\xi_{\mathrm{ph},{\textbf{\itshape{k}}}}} \equiv {\epsilon_{\mathrm{ph},{\textbf{\itshape{k}}}}} - \mu$, where ${\tilde\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}} \equiv {\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}} - \sum_{{\textbf{\itshape{k}}}'} U'_{{\textbf{\itshape{k}}}-{\textbf{\itshape{k}}}'} n_{\mathrm{e(h)},{\textbf{\itshape{k}}}'}$ denotes the energy dispersion of electrons (holes) renormalized by the repulsive electron-electron (hole-hole) Coulomb interaction, the first (second) term in Eq. (\[eq:HCoul\]). The well-known bandgap renormalization (BGR) in semiconductor physics is included in ${\tilde\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}}$. In contrast, ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \equiv g^{*}a_0 + \sum_{{\textbf{\itshape{k}}}'}U'_{{\textbf{\itshape{k}}}-{\textbf{\itshape{k}}}'}p_{{\textbf{\itshape{k}}}'}$ results from the attractive electron-hole Coulomb interaction, the third term in Eq. (\[eq:HCoul\]), and is called the generalized Rabi frequency [@HaugKoch]. ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$ has a role in forming e-h pairs as can be seen in Eq. (\[eq:HMF\]).

The mean-field Hamiltonian is thus obtained. However, note that the expectation values of ${\langle \hat{O} \rangle}$ (i.e. $\a0$, $\pk$, $\nek$, and $\nhk$) are included in $\hat{H}_{\mathrm{S}}^{\mathrm{MF}}$. For self-consistensy, therefore, the following relation should be satisfied : $${\langle \hat{O} \rangle} = {\mathrm{Tr}[\hat{O}\hat\rho^{\mathrm{MF}}({\langle \hat{O} \rangle})]}.
\label{eq:SCE}$$ Here, $\hat\rho^{\mathrm{MF}}$ is the density operator determiend by using $\hat{H}_{\mathrm{S}}^{\mathrm{MF}}$. The BCS theory and the MSBE shown below are obtained from this self-consistent equation.

BCS theory for exciton-polariton condensation {#subsec:BCS}
---------------------------------------------

First, we assume that the exciton-polariton system is in equilibrium. According to the equilibrium statistical mechanics, the density operator $\hat\rho^{\mathrm{MF}}$ at temperature $T$ can be described as $$\begin{aligned}
\hat\rho^{\mathrm{MF}}=\hat\rho_{\mathrm{eq}}^{\mathrm{MF}}\equiv\frac{1}{Z}\exp(-\beta\hat{H}_{\mathrm{S}}^{\mathrm{MF}}),
\label{eq:equilibrium}\end{aligned}$$ where $Z \equiv \mathrm{Tr}[\exp(-\beta\hat{H}_{\mathrm{S}}^{\mathrm{MF}})]$ and $\beta \equiv 1/T$. In this case, $\mu$ is a given parameter equivalent to the chemical potential, as mentioned above. With ${\epsilon_{\mathrm{e},{\textbf{\itshape{k}}}}} = {\epsilon_{\mathrm{h},{\textbf{\itshape{k}}}}}$ for simplicity, the self-consistent equations obtained from Eqs. (\[eq:HMF\])-(\[eq:equilibrium\]) are $$\begin{aligned}
\a0=\sum_{{\textbf{\itshape{k}}}'}\frac{g}{\xi_{\mathrm{ph},0}}p_{{\textbf{\itshape{k}}}'},\ \pk=\frac{{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}}{2{E_{{\textbf{\itshape{k}}}}}}\tanh \left( \frac{\beta {E_{{\textbf{\itshape{k}}}}}}{2} \right),
\label{eq:ap}\\
\nek=\nhk=\frac{1}{2} \left\{ 1-\frac{{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+}}{{E_{{\textbf{\itshape{k}}}}}}\tanh \left( \frac{\beta {E_{{\textbf{\itshape{k}}}}}}{2} \right)  \right\},
\label{eq:number}\end{aligned}$$ where ${\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{\pm} \equiv ({\tilde{\xi}_{\mathrm{e},{\textbf{\itshape{k}}}}} \pm {\tilde{\xi}_{\mathrm{h},{\textbf{\itshape{k}}}}})/2 $ and ${E_{{\textbf{\itshape{k}}}}} \equiv [({\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+})^2 + |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2]^{1/2}$. In the derivation, Bogoliubov transformations of $\ope_{{\textbf{\itshape{k}}}}$ and $\oh_{{\textbf{\itshape{k}}}}$ can be applied to the first line in Eq. (\[eq:HMF\]) for diagonalization, while a displacement of $\oa_{0}$ to the second line, because the Hilbert space of the first (second) line of Eq. (\[eq:HMF\]) is spanned only by the electron and hole (photon) degrees of freedom.

The gap equation, which is formally equivalent to the BCS theory for superconductors, can then be obtained by substituting Eq. (\[eq:ap\]) into the definition of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$: $$\displaystyle {\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} = \sum_{{\textbf{\itshape{k}}}'}U_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}}^{\mathrm{eff}}\frac{{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}{'}}}}{2{E_{{\textbf{\itshape{k}}}{'}}}} \tanh \left( \frac{\beta {E_{{\textbf{\itshape{k}}}{'}}}}{2} \right).
\label{eq:gap}$$ In this context, ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$ is an order parameter in the exciton-polariton system as well as in the superconducting system. $U_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}}^{\mathrm{eff}} \equiv |g|^2/\xi_{\mathrm{ph},0} + U'_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{k}}}}$ represents an effective attractive e-h interaction, from which one can find that photon-mediated process also contributes the attractive interaction. Notice that Eqs. (\[eq:number\]) and (\[eq:gap\]) are simultaneous equations with the unknown variables $\nek (=\nhk)$ and ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$. Especially for $T = 0$, this treatment is known to cover the equilibrium phases from the BEC to the BCS states [@Kamide10; @Byrnes10; @Comte82].

![ A schematic picture of the model including environments. The e-h system is excited and thermalized by the pumping baths by exchanging carriers. Photons in the system are lost into the vacuum. []{data-label="Model"}](Model.eps){width="0.50\linewidth"}

MSBE for semiconductor lasers {#subsec:MSBE}
-----------------------------

Next, a treatment based on the MSBE is explained for the discussion of the semiconductor laser, which is characterized by nonequilibrium. In contrast to the BCS theory, therefore, the effects of environments (Fig. \[Model\]) cannot be neglected for lasing; the excitation and thermalization of the e-h system and the loss of photons from the microcavity. For this reason, the dynamics of the total density operator $\hat\rho^{\mathrm{MF}}$ is discussed by writing the total mean-field Hamiltonian $\hat{H}^{\mathrm{MF}} \equiv \hat{H}_{\mathrm{S}}^{\mathrm{MF}}+\hat{E}$ with the couplings to the environments $\hat{E}$. Since $\ii\partial_t \hat\rho^{\mathrm{MF}} = [\hat{H}^{\mathrm{MF}},\hat\rho^{\mathrm{MF}}]$ in the Schrödinger picture, a time drivative of Eq. (\[eq:SCE\]) yields $$\ii\partial_t{\langle \hat{O} \rangle} = {\mathrm{Tr}[[\hat{O},\hat{H}_{\mathrm{S}}^{\mathrm{MF}}]\hat\rho^{\mathrm{MF}}]}+{\mathrm{Tr}[[\hat{O},\hat{E}]\hat\rho^{\mathrm{MF}}]},
\label{eq:SCE2}$$ where ${\mathrm{Tr}[\hat{A}\hat{B}]}={\mathrm{Tr}[\hat{B}\hat{A}]}$ is used. The MSBE is then obtained when the first term is derived from Eq. (\[eq:HMF\]) and the second term is replaced by phenomenological relaxation terms: $$\begin{aligned}
\partial_t\a0&=&-\ii \xi_{\mathrm{ph},0} \a0 + \ii g \textstyle{\sum_{{\textbf{\itshape{k}}}}} \pk - \kappa \a0,
\label{eq:MSBE_a0} \\
\partial_t\pk&=&-2\ii{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+}\pk-\ii{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} N_{{\textbf{\itshape{k}}}}-2\gamma (\pk - \pk^{0}),
\label{eq:MSBE_pk} \\
\partial_t n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}&=&-2\mathrm{Im}[{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} p_{{\textbf{\itshape{k}}}}^{*}]-2\gamma(n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}-n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}^{0}),
\label{eq:MSBE_nehk}\end{aligned}$$ where the last term in each equation is the relaxation term and $N_{{\textbf{\itshape{k}}}} \equiv \nek + \nhk -1$ denotes the degree of the population inversion.[^4] $\pk^{0}$ and $n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}^{0}$ are defined as $$\pk^{0} \equiv 0, \hspace{5pt} n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}^{0} \equiv f_{\mathrm{e(h)},{\textbf{\itshape{k}}}} 
\label{eq:RTA}$$ where $f_{\mathrm{e(h)},{\textbf{\itshape{k}}}} \equiv [1 + \exp \{ \beta ( {\tilde\epsilon_{\mathrm{e(h)},{\textbf{\itshape{k}}}}}-\mu^{\mathrm{B}}_{\mathrm{e(h)}}) \} ]^{-1}$ is the Fermi distribution with the chemical potential $\mu^{\mathrm{B}}_{\mathrm{e(h)}}$ of the electron (hole) pumping bath. The phenomenological approximation shown here is called the relaxation approximation [@Henneberger92]. Each relaxation term suggests that the photon field $a_{0}$ decays with a rate of $\kappa$, the distribution function $n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}$ is driven to approach the Fermi distribution $f_{\mathrm{e(h)},{\textbf{\itshape{k}}}}$ (Fig. \[Energy\]), i.e. thermalization (Fig. \[RelaxationRedistribution\]), and $\pk$ decays due to thermalization-induced dephasing.

Solutions for the laser action can then be obtained by determining the unknown variables $\a0$, $\pk$, $\nek$, $\nhk$, and $\mu$ in Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) under a steady-state condition $\partial_t{\langle \hat{O} \rangle} = 0$. Again, we emphasize that $\mu$ is a unknown variable corresponding to the laser frequency in the steady-state MSBE, in contrast to the BCS theory. This is equivalent to find an appropriate frequency with which the lasing oscillation of $\a0$ and $\pk$ seems to remain stationary on the rotating frame.

![ Energy dispersions of electrons and holes. The distribution functions $n_{\mathrm{e},{\textbf{\itshape{k}}}}$ and $n_{\mathrm{h},{\textbf{\itshape{k}}}}$ are driven to approach the Fermi distributions by the respective pumping baths. []{data-label="Energy"}](EnergyDispersion.eps){width="0.50\linewidth"}

BEC-BCS-LASER crossover theory {#sec:BBL}
==============================

In the exciton-polariton system, as shown in Section \[sec:BCS and MSBE\], the BCS theory and the MSBE are theoretical frameworks starting from the common Hamiltonians with the same mean-field approximation. However, the difference is the way of deriving the self-consistent equations. In the case of the BCS theory, $\hat\rho^{\mathrm{MF}}$ is directly described by $\hat{H}_{\mathrm{S}}^{\mathrm{MF}}$ (Eq. (\[eq:equilibrium\])). In contrast, in the case of the MSBE, Eq. (\[eq:SCE2\]) is used to introduce the phenomenological relaxation terms. We note, however, that any assumption is not used for $\hat\rho^{\mathrm{MF}}$ in Eq. (\[eq:SCE2\]), which indicates that the MSBE may incorporate the BCS theory at least in principle.

In this context, an approach to derive the BCS theory from the MSBE should be discussed briefly. We first consider a situation where the effects of the environments are completely neglected, which is equivalent to set $\kappa = \gamma = 0$ in the MSBE. However, in this case, the BCS theory cannnot be derived because there is no term to drive the system into equilibrium in the MSBE.[^5] A natural condition to consider is physically a limit of $\gamma \rightarrow 0^{+}$ after $\kappa \rightarrow 0$ because the system should be thermalized even though the effects of environments are decreased. Unfortunately, however, the MSBE does not recover the BCS theory even by taking this limit. The relationship between the BCS theory and the MSBE is thus discontinuous in spite of the similarities of the two frameworks. Obviously, the phenomenological relaxation approximation causes such a problem.

In regard to this problem, we have recently constructed a unified framework [@Yamaguchi12] by using a nonequilibrium Green’s function approach [@Szymanska06; @Szymanska07; @Keeling10]. The framework, at first, takes an integral form of simultaneous equations and seems quite different from the MSBE (see also Appendix II). However, by rearranging the equations with particular attention to the problem mentioned above, all of important changes can successfully be incorporated in the relaxation terms in the MSBE [@Yamaguchi13]. The result simply replaces Eq. (\[eq:RTA\]) by $$\begin{aligned}
&\pk^{0} \equiv \ii\int\frac{\mathrm{d}\nu}{2\pi}\left[ {G^{\mathrm{R}}_{\mathrm{12},{{\textbf{\itshape{k}}}}}}(\nu)\{1-f^{\mathrm{B}}_{\mathrm{h}}(-\nu)\} - {G^{\mathrm{R*}}_{\mathrm{21},{{\textbf{\itshape{k}}}}}}(\nu) f^{\mathrm{B}}_{\mathrm{e}}(\nu) \right], \nonumber\\
&n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}^{0} \equiv \int\frac{\mathrm{d}\nu}{2\pi}f^{\mathrm{B}}_{\mathrm{e(h)}}(\nu)A_{\mathrm{11(22)}}(\pm\nu;{\textbf{\itshape{k}}}),
\tag{\ref{eq:RTA}$'$} \label{eq:RTA'}\end{aligned}$$ where $f^{\mathrm{B}}_{\mathrm{e(h)}}(\nu) \equiv [\exp \{ \beta(\nu - \mu^{\mathrm{B}}_{\mathrm{e(h)}} + \mu/2) \} + 1]^{-1}$ is the Fermi distribution of the electron (hole) pumping bath. ${G^{\mathrm{R}}_{\mathrm{\alpha\alpha'},{{\textbf{\itshape{k}}}}}}(\nu)$ is called the retarded Green’s function and described by elements of a matrix $$\begin{aligned}
G^{\mathrm{R}}_{{\textbf{\itshape{k}}}}(\nu)=\left(
\begin{array}{cc}
\nu - \tilde{\xi}_{\mathrm{e},{\textbf{\itshape{k}}}} + \ii\gamma & {\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \\
{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}^{*} & \nu + \tilde{\xi}_{\mathrm{h},{\textbf{\itshape{k}}}} + \ii\gamma \\
\end{array}\right)^{-1}.
\label{eq:GR}\end{aligned}$$ On the other hand, $A_{11(22)}(\nu;{\textbf{\itshape{k}}})$ is called the single-particle spectral function and defined as $$A_{\alpha\alpha'}(\nu;{\textbf{\itshape{k}}}) \equiv \ii (G^{\mathrm{R}}_{\alpha\alpha',{\textbf{\itshape{k}}}}(\nu) - G^{\mathrm{R}*}_{\alpha'\alpha,{\textbf{\itshape{k}}}}(\nu)).
\label{eq:A}$$ Here, $A_{11(22)}(\nu;{\textbf{\itshape{k}}})$ means the density of states for electron-like (hole-like) quasi-particles with the energy $\nu$ and wave number ${\textbf{\itshape{k}}}$.

Some readers might feel difficult to understand the formalism because the above definitions are unique to the Green’s function approach. However, all we have to do is the replacement of Eq. (\[eq:RTA\]) by Eq. (\[eq:RTA’\]). The unknown variables are still $\a0$, $\pk$, $\nek$, $\nhk$, and $\mu$, that is, the same as the MSBE. In this sense, the obtained equations are quite simple, which is one of strong points of this formalism. From the viewpoint of the Green’s funciton, it is relatively easy to understand the physical meaning of Eq. (\[eq:RTA’\]) due to the clear form; the energy integral of (distribution)$\times$(density of states).[^6] We refer to such a formalism as the BEC-BCS-LASER crossover theory.

Now, this formalism enables us to cleary understand the standpoint of the BCS theory. For this purpose, let us discuss the limit of equilibrium, based on the idea described above. In the followings, however, ${\epsilon_{\mathrm{e},{\textbf{\itshape{k}}}}} = {\epsilon_{\mathrm{h},{\textbf{\itshape{k}}}}}$ and a charge neutrality $\mu^{\mathrm{B}}_{\mathrm{e}} = \mu^{\mathrm{B}}_{\mathrm{h}}$ are assumed for simplicity. First, in the limit of $\kappa \rightarrow 0$, one can prove $\mu = \mu_{\mathrm{B}} (\equiv \mu^{\mathrm{B}}_{\mathrm{e}} + \mu^{\mathrm{B}}_{\mathrm{h}})$. This is the same as treating $\mu$ as a given parameter, and physically, means that the system reaches in chemical equilibrium with the pumping baths because there is no photon loss. The BCS theory is then derived after taking the limit of $\gamma \rightarrow 0^{+}$, where the integrals in Eq. (\[eq:RTA’\]) can be performed analytically. In this derivation, $\gamma \neq 0$ is required to be canceled down even though $\gamma$ does not appear in the final expression. This means that thermalization is essential to recover the equilibrium theory.

Thus, the BCS theory can be derived from the presented theory in the equilibrium limit. However, in some sense, this situation is physically trivial; the situation is not limited to such a trivial one for the system to be in equilibrium. Even under a condition where photons are continuously lost, it may be still possible to identify the system as being in equilibrium (quasi-equilibrium) as long as the e-h system is excited and thermalized. A true advantage of the above-presented framework becomes obvious in such a situation rather than in the trivial one. In this case, $\mu$ is still equivalent to the chemical potential but $\mu_{\mathrm{B}} > \mu$ because the system is influenced by the photon loss. As a result, $\mu$ becomes a unknown variable again. Furthermore, such a quasi-equilibrium condition can easily be obtained from Eqs. (\[eq:RTA’\])-(\[eq:A\]) as[^7]

(I) $\min[2{E_{{\textbf{\itshape{k}}}}}] \gtrsim \mu_{\mathrm{B}} - \mu + 2\gamma + 2T$.

Here, $\min[2{E_{{\textbf{\itshape{k}}}}}]$ is the minimum energy required for breaking e-h bound pairs and $\mu_{\mathrm{B}}-\mu > 0$ suggests that there is continuous particle flow from the pumping baths into the system.[^8] We can then interpret the condition (I); this is a condition that the particle flux, thermalization-induced dephasing (= $2\gamma$), and temperature effect (= $2T$), do not contribute to the dissociations of the e-h pairs.

However, the system can no longer be in quasi-equilibrium when nonequilibrium effect becomes significant. Let us therefore consider a situation where the MSBE, i.e. the physics of the semiconductor laser, becomes important. Such a condition can be found from Eqs. (\[eq:RTA’\])-(\[eq:A\]) as

(II) $\mu_{\mathrm{B}} - \mu \gtrsim \min[2{E_{{\textbf{\itshape{k}}}}}] + 2\gamma + 2T$,

because $f^{\mathrm{B}}_{\mathrm{e(h)}}(\pm\nu) \simeq f^{\mathrm{B}}_{\mathrm{e(h)}}({\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+})$ turns out to be a good approximation in Eq. (\[eq:RTA’\]) for ${\textbf{\itshape{k}}}$-resions satisfying

(II$'$) $~\mu_{\mathrm{B}} - \mu \gtrsim 2{E_{{\textbf{\itshape{k}}}}} + 2\gamma + 2T$.

Note that there are such ${\textbf{\itshape{k}}}$-resions whenever the condition (II) is fulfilled. As a result, we can obtain $p_{{\textbf{\itshape{k}}}}^0 \cong 0$ and $n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}^0 \cong f_{\mathrm{e(h)},{\textbf{\itshape{k}}}}$ which recovers the MSBE. However, we stress that the condition (II$'$) depends on the wave number ${\textbf{\itshape{k}}}$; there remain ${\textbf{\itshape{k}}}$-regions still described by the BCS theory. The MSBE and the BCS theory are, thus, coupled with each other in a strict sense. In this context, the lasing can be referred to as the BCS-coupled lasing when this viewpoint is emphasized. At the same time, the physical meaning of $\mu$ changes into the oscillating frequency of the laser action.

Second thresholds, band renormalization, and gain spectra {#sec:Results}
=========================================================

Figure \[Results\] shows the number of coherent photons in the cavity $|a_{0}|^2$ and the frequancy $\mu$ as a function of $\mu_{\mathrm{B}}$ calculated by our formalism.[^9] Plots are colour coded by red (blue) when the quasi-equilibrium condition (I) (the lasing condition (II)) is satisfied, while by green when neither of the conditions is satisfied. In Fig. \[Results\](a), $|a_{0}|^2$ arises with increasing $\mu_{\mathrm{B}}$, the point of which is called the first threshold. In this situation, the system is in quasi-equilibrium regime (red) and $\mu$ is around the lower polariton level[^10] $E_{\mathrm{LP}}$ in Fig. \[Results\](b). The first threshold therefore means that the exciton-polariton BEC is caused because the chemical potential of the system reaches the lowest energy of the exciton polariton, $E_{\mathrm{LP}}$.

With further increase of $\mu_{\mathrm{B}}$, the system changes from the quasi-equilibrium regime (red) into the lasing regime (blue) through a crossover regime (green). Around the crossover regime in Fig. \[Results\](a), a second threshold can be seen where the number of coherent photons grows rapidly again. $\mu$ is then blue-shifted from $E_{\mathrm{LP}}$ into the bare cavity level $E_{\mathrm{cav}}$. Furthermore, the kinetic hole burning can be seen in the distribution function of electrons $\nek$ (the blue arrow in the inset to Fig. \[Results\](a)). These resuls demonstrate that the exciton-polariton BEC has smoothly changed into the semiconductor laser with the second threshold.

In experiments [@Balili09; @Nelsen09; @Dang98; @Tempel12-1; @Tempel12-2; @Tsotsis12; @Kammann12], the second threshold and the blue shift has been reported since more than 10 years ago, the mechanism of which has been attributed to a shift into the weak coupling regime due to dissociations of Coulomb-bound e-h pairs (excitons); the lasing phase is then achieved as a result. However, there is no convincing discussion why such dissociations lead to nonequilibration essential for lasing.

![ Numerical results of (a) the coherent photon number in the cavity $|a_{0}|^2$ and (b) the frequency $\mu$ as a function of $\mu_{\mathrm{B}}$. Plots are colour coded by red and blue when satisfying the quasi-equilibrium condition (I) and the lasing condition (II), respectively. Green colours are used when neither of them are satisfied. $\mu$ represents the chemical potential in the quasi-equilibrium regime (red) but the laser frequency in the lasing regime (blue). Inset: the distribution function of electrons $n_{\mathrm{e},{\textbf{\itshape{k}}}}$ (black) and the polarization $\pk$ (red). In the lasing regime (B), a characteristic dip can be seen in the distribution (the bule arrow), which is known as one of the signatures of lasing and called the kinetic hole burning. []{data-label="Results"}](Results.eps){width="0.70\linewidth"}

According to our formalism, this empirical picture can be investigated and shown to be incorrect. This is because, even in the lasing regime, there are gaps around $\pm\mu/2$ in the renormalized band structure as shown in the left of Fig. \[Gain\]. An analytical form of $A_{11(22)}(\nu;{\textbf{\itshape{k}}})$, obtained by Eqs. (\[eq:GR\]) and (\[eq:A\]), enables us to conveniently study the renormalized band: $$A_{11(22)}(\nu;{\textbf{\itshape{k}}}) = 2|u_{{\textbf{\itshape{k}}}}|^2\frac{\gamma}{(\nu-{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{-} \mp {E_{{\textbf{\itshape{k}}}}})^2+\gamma^2}
            +2|v_{{\textbf{\itshape{k}}}}|^2\frac{\gamma}{(\nu-{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{-} \pm {E_{{\textbf{\itshape{k}}}}})^2+\gamma^2}.
\label{eq:A2}$$ Here, $u_{{\textbf{\itshape{k}}}}$ and $v_{{\textbf{\itshape{k}}}}$ are the Bogoliubov coefficients defined as $$u_{{\textbf{\itshape{k}}}} \equiv \sqrt{\frac{1}{2}+\frac{{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+}}{2{E_{{\textbf{\itshape{k}}}}}}}, \hspace{0.5cm}
v_{{\textbf{\itshape{k}}}} \equiv e^{\ii\theta_{{\textbf{\itshape{k}}}}}\sqrt{\frac{1}{2}-\frac{{\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{+}}{2{E_{{\textbf{\itshape{k}}}}}}},
\label{eq:Bogoliubov}$$ with $\theta_{{\textbf{\itshape{k}}}} \equiv \arg({\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}})$. These equations have remarkable similarities to the BCS theory in superconductors [@Abrikosov75; @Yamaguchi13]. Therefore, it is clear that the gaps are opened around $\pm\mu/2$ with the magnitude of $\min[2{E_{{\textbf{\itshape{k}}}}}]$ when ${\tilde{\xi}_{\mathrm{eh},{\textbf{\itshape{k}}}}}^{-}=0$ i.e. ${\epsilon_{\mathrm{e},{\textbf{\itshape{k}}}}} = {\epsilon_{\mathrm{h},{\textbf{\itshape{k}}}}}$ with $\mu^{\mathrm{B}}_{\mathrm{e}} = \mu^{\mathrm{B}}_{\mathrm{h}}$. Note, however, that the unknown variables contained in Eqs. (\[eq:A2\]) and (\[eq:Bogoliubov\]) are determined by the BEC-BCS-LASER crossover theory (Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) with Eqs. (\[eq:RTA’\])-(\[eq:A\])) rather than the BCS theory. In the BCS phase, the existence of the gap around $\pm\mu/2$ means the formation of Cooper pairs around the Fermi level because $\pm\mu/2$ is equivalent to the Fermi level. In contrast, in the lasing phase, $\pm\mu/2$ corresponds to the laser frequency.[^11] Thus, the gap indicates the formation of bound e-h pairs by mediating photons around the laser frequency. The semiconductor laser in Fig. \[BEC-BCS-LASER\] is drawn along this picture, where the e-h pairs are explicitly depicted.

![Left; A renormalized band structure in a lasing phase (the point B in Fig. \[Results\]). The gaps are opened around $\pm\mu/2$ with the magnitude of $\min[2{E_{{\textbf{\itshape{k}}}}}]$. Right; Optical gain spectra for the exciton-polariton BEC (A) and for the lasing phase (B). Panels (A) and (B) correspond to the point A and B in Fig. \[Results\], respectively. Aqua (pink) represents the gain (absorption). Panels (A) and (B) are reproduced from [@Yamaguchi13]. []{data-label="Gain"}](Gain.eps){width="0.80\linewidth"}

Such a “lasing gap” is, at least in principle, measureable in the optical gain spectrum $G(\omega)$ by irradiating probe light with frequency $\omega$ because $G(\omega)$ is strongly affected by the renormalized band structure in general. As a result, in the gain spectrum of the lasing phase (Fig. \[Gain\](B)), there appears a transparent region originating from the gap. The optical gain spectrum is thus one of important ways for the verification of the lasing gap.

In addition, we note that behaviors of the gain spectra vary drastically when the exciton-polariton BEC is changed to the laser phase. By comparing Figs. \[Gain\](A) and (B), for example, one can find only absorption but no gain in Fig. \[Gain\](A). This is mainly because there is no ${\textbf{\itshape{k}}}$-region with inverted population $N_{k}>0$ $(\Leftrightarrow n_{\mathrm{h},k} = n_{\mathrm{e},k} > 0.5)$ in Fig. \[Results\](A). In contrast, optical gain is caused in Fig. \[Gain\](B) because there are ${\textbf{\itshape{k}}}$-regions with $N_{{\textbf{\itshape{k}}}} > 0$ in Fig. \[Results\](B). Thus, the existence of the gain after the second threshold gives us important information to identify the phases in the system.[^12]

Conclusions and perspectives {#sec:Conclusions}
============================

In this contribution, we have presented a brief explanation of the BCS theory and the MSBE in the exciton-polariton system, to highlight their similarities and differences. We have then shown a framework of describing the BCS theory (the BEC and BCS phases) and the MSBE (the semiconductor laser) in a unified way. As a result, the existence of bound e-h pairs in the lasing phase as well as the lasing gap have been pointed out. The results presented here are the physics elucidated for the first time by considering the BEC, BCS, and Laser phases in a unified way. However, for example, effects of spontaneous emission [@Scully97] and pure dephasing [@Yamaguchi12-2] are still unclear. In this respect, further studies are needed for a full understanding of this system. Experimantal studies are also important, in particular, in a high density regime [@Balili09; @Nelsen09; @Dang98; @Tempel12-1; @Tempel12-2; @Tsotsis12; @Kammann12; @Horikiri13].

Although we have focused on the exciton-polariton system in this contribution, we finally would like to emphasize that this system has a close relationship with superconductors and the Feshbach resonance in cold atom systems because interacting Fermi and Bose particles play important roles in the formation of ordered phases. In this sense, it would be interesting to study the lasing gap by terahertz pulses in a manner similar to superconductors [@Matsunaga13; @Papenkort07]. Inclusions of the e-h center-of-mass fluctuations with mass imbalance are also important, as discussed in the cold atom systems [@Hanai13-1], because these effects cannot be taken into account within the mean-field approximation. We further note that fundamental problems of the nonequilibrium statistical physics are also included in this system in the sense of providing a bridge between the equilibrium and the nonequilibrium phases. We hope that our approach also stimulates new studies in a wide range of such fields.

The authors are grateful to K. Kamide, R. Nii, Y. Yamamoto, T. Horikiri, Y. Shikano, Y. Matsuo, T. Yuge, and M. Bamba for fruitful discussions. This work is supported by the JSPS through its FIRST Program, and DYCE, KAKENHI No. 20104008.

Appendix I: Excitonic effects in the low density limit {#appendix-i-excitonic-effects-in-the-low-density-limit .unnumbered}
======================================================

In semiconductor exciton-polariton systems, the excitonic effects play quite important roles in the formation of Coulomb-bound e-h pairs (excitons) and exciton-polaritons. In this Appendix, we, therefore, confirm the excitonic effects in our formalism [^13].

For this purpose, we now assume that the density of electrons and holes are sufficiently low ($n_{\mathrm{e},{\textbf{\itshape{k}}}}, n_{\mathrm{h},{\textbf{\itshape{k}}}} \ll 1$ or $N_{{\textbf{\itshape{k}}}} \cong -1$) with no pumping and loss ($\gamma = 0$ and $\kappa = 0$). Under this condition, $2\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}}$ in Eq. (\[eq:MSBE\_pk\]) can be written as $$\begin{aligned}
2\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}} &= \epsilon_{\mathrm{e},{\textbf{\itshape{k}}}} + \epsilon_{\mathrm{h},{\textbf{\itshape{k}}}} - \mu = \frac{k^2}{2m_{\mathrm{r}}}+E_g-\mu,
\label{eq:temp1}\end{aligned}$$ where $1/m_{\mathrm{r}}=1/m_{\mathrm{e}}+1/m_{\mathrm{h}}$. Then, Eqs. (\[eq:MSBE\_a0\]) and (\[eq:MSBE\_pk\]) can be described as $$\begin{aligned}
0&=-(\epsilon_{\mathrm{ph},0}-\mu)a_0+ g \sum_{{\textbf{\itshape{k}}}} p_{{\textbf{\itshape{k}}}},
\label{eq:a0_}\\
0&=-\left( \frac{k^2}{2m_{\mathrm{r}}}+E_g-\mu \right) p_{{\textbf{\itshape{k}}}} + g^{*}a_0 + \sum_{{\textbf{\itshape{k}}}'}U'_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{k}}}}p_{{\textbf{\itshape{k}}}'},
\label{eq:pk_}\end{aligned}$$ where the definition of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \equiv g^{*}a_0 + \sum_{{\textbf{\itshape{k}}}'}U'_{{\textbf{\itshape{k}}}-{\textbf{\itshape{k}}}'}p_{{\textbf{\itshape{k}}}'}$ is used. Especially, for $g = 0$ in Eq. , we obtain $$\begin{aligned}
\frac{k^2}{2m_{\mathrm{r}}} p_{{\textbf{\itshape{k}}}} - \sum_{{\textbf{\itshape{k}}}'}U'_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{k}}}}p_{{\textbf{\itshape{k}}}'}=-(E_g-\mu) p_{{\textbf{\itshape{k}}}},
\label{eq:Wannier}\end{aligned}$$ which is nothing but the Schr$\ddot{\mathrm{o}}$dinger equation in ${\textbf{\itshape{k}}}$-space for the single exciton bound state [@Yamaguchi12; @Comte82; @HaugKoch]. This means that the Coulomb-bound e-h pairs (excitons) can be formed in the low density limit in the presented formalism. In such a case, $p_{{\textbf{\itshape{k}}}}$ can be described by the bound state e-h pair wave-function $\phi_{{\textbf{\itshape{k}}}}$ ($p_{{\textbf{\itshape{k}}}} = \eta \phi_{{\textbf{\itshape{k}}}}$ with $\sum_{{\textbf{\itshape{k}}}}|\phi_{{\textbf{\itshape{k}}}}|^2 = 1$) with $\mu = E_{\mathrm{ex}}$, where $E_{\mathrm{ex}}$ is the energy level of the exciton and the binding energy corresponds to $E_{\mathrm{g}} - E_{\mathrm{ex}}$. The formation of the exciton is, thus, includeed in the theory.

Although $p_{{\textbf{\itshape{k}}}}$ is changed from the exciton wave-function $\phi_{{\textbf{\itshape{k}}}}$ by the photon-mediated attraction in the case of $g \neq 0$, it is instructive to consider the case where such an effect is not so large. In this limit, by substituting $p_{{\textbf{\itshape{k}}}} \cong \eta\phi_{{\textbf{\itshape{k}}}}$ into Eqs.  and , we obtain $$\begin{aligned}
0&=(\mu - E_{\mathrm{cav}})a_0 + g_{\mathrm{ex}}\eta, \\
0&=(\mu - E_{\mathrm{ex}})\eta + g_{\mathrm{ex}}^{*}a_0, 
\label{eq:eigen}\end{aligned}$$ where $g_{\mathrm{ex}} \equiv g\sum_{{\textbf{\itshape{k}}}}\phi_{{\textbf{\itshape{k}}}}=g\phi_{\mathrm{ex}}({\textbf{\itshape{r}}}=0)$ is the coupling constant renormalized by the exciton wave-function. Then, $\mu$ is given by one of the eigenvalues of these two coupled equations, which are the eigen-energies of the upper and lower polaritons: $$\begin{aligned}
E_{\mathrm{UP/LP}} = \frac{E_{\mathrm{cav}}+E_{\mathrm{ex}} \pm \sqrt{(E_{\mathrm{cav}}-E_{\mathrm{ex}})^2+4|g_{\mathrm{ex}}|^2}}{2}. 
\label{eq:UPLP}\end{aligned}$$ Here, $E_{\mathrm{UP}}$ and $E_{\mathrm{LP}}$ in Eq.  are the well-known expressions obtained when the excitons are treated as simple bosons [@HaugKoch]. This means that the formation of exciton-polaritons are also included in the theory. The excitonic effects are, thus, taken into account in our formalism within the mean-field approximation. We note that the procedure shown here is basically the same as Section 2.1.2 in Ref. [@Yamaguchi12].

Appendix II: Proof of equivalence {#appendix-ii-proof-of-equivalence .unnumbered}
=================================

In the main text, we have mentioned that the formalism in Ref. [@Yamaguchi12] seems quite different from the MSBE. This formalism can be described by the follwoing simultaneous equations with the unknown variables of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$, $\nek$, $\nhk$, and $\mu$: $$\begin{gathered}
{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} = \sum_{{\textbf{\itshape{k}}}'} U^{\mathrm{eff},\kappa}_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}} {\mathnormal{\Delta}_{{\textbf{\itshape{k}}}{'}}} \int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}'}(\nu) \\
          \times \left \{ (F_{\mathrm{e}}^{\mathrm{B}}(\nu) + F_{\mathrm{h}}^{\mathrm{B}}(\nu)) (\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}'}) 
          + (F_{\mathrm{e}}^{\mathrm{B}}(\nu) - F_{\mathrm{h}}^{\mathrm{B}}(\nu)) (\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}'} + \ii\gamma) \right\},
\label{eqa:Delta}\end{gathered}$$ $$n_{\mathrm{e(h)},{\textbf{\itshape{k}}}} = \frac{1}{2} \mp \int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}}(\nu)
        \left\{ F_{\mathrm{e(h)}}^{\mathrm{B}}(\nu)[(\nu\pm\tilde\xi_{\mathrm{h(e)},{\textbf{\itshape{k}}}})^2+\gamma^2] + F_{\mathrm{h(e)}}^{\mathrm{B}}(\nu)|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2 \right\},
\label{eqa:nehk}$$ where $U_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}}^{\mathrm{eff},\kappa} \equiv |g|^2/(\xi_{\mathrm{ph},0}-\ii\kappa)+U'_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{k}}}}$ and $$L_{{\textbf{\itshape{k}}}}(\nu) \equiv \frac{\gamma}{[(\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}}-{E_{{\textbf{\itshape{k}}}}})^{2}+\gamma^2][(\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}}+{E_{{\textbf{\itshape{k}}}}})^{2}+\gamma^2]}.
\label{eqa:L}$$ $F_{\mathrm{e}}^{\mathrm{B}}(\nu)$ and $F_{\mathrm{h}}^{\mathrm{B}}(\nu)$ are respectively defined as $$\begin{aligned}
F_{\mathrm{e}}^{\mathrm{B}}(\nu) & \equiv & \tanh \left( \frac{\beta[\nu - \mu^{\mathrm{B}}_{\mathrm{e}} + \mu/2]}{2} \right) = 1-2f^{\mathrm{B}}_{\mathrm{e}}(\nu),\\
\label{eqa:Fe}
F_{\mathrm{h}}^{\mathrm{B}}(\nu) & \equiv & \tanh \left( \frac{\beta[\nu + \mu^{\mathrm{B}}_{\mathrm{h}} - \mu/2]}{2} \right) = 2f^{\mathrm{B}}_{\mathrm{h}}(-\nu)-1.
\label{eqa:Fh}\end{aligned}$$ In this Appendix II, therefore, we prove that Eqs. (\[eqa:Delta\])-(\[eqa:Fh\]) are equivalent to Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) with Eqs. (\[eq:RTA’\])-(\[eq:A\]) under the steady-state condition $\partial_t{\langle \hat{O} \rangle} = 0$.

For this purpose, we here note that the following sum rule is satisfied for the single-partice spectral function: $$\int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi}A_{\alpha\alpha'}(\nu;{\textbf{\itshape{k}}}) = \delta_{\alpha,\alpha'}.
\label{eqa:SumRule}$$ This relation can be confirmed by the direct integration of $A_{\alpha\alpha'}(\nu;{\textbf{\itshape{k}}})$ described by elements of a matrix $$\begin{aligned}
A(\nu;{\textbf{\itshape{k}}}) &= \frac{-2}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}\left(
\begin{array}{cc}
\operatorname{Im}[D^{*}_{{\textbf{\itshape{k}}}}(\nu)(\nu+\tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}}+\ii\gamma)] & \operatorname{Im}[D_{{\textbf{\itshape{k}}}}(\nu)]{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \\
\operatorname{Im}[D_{{\textbf{\itshape{k}}}}(\nu)]{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}^{*} & \operatorname{Im}[D^{*}_{{\textbf{\itshape{k}}}}(\nu)(\nu-\tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}}+\ii\gamma)] \\
\end{array}\right) \nonumber \\
&= \frac{-2\gamma}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}\left(
\begin{array}{cc}
-[(\nu+\tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}})^2+\gamma^2+|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2] & 2{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}(\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}}) \\
2{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}^{*}(\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}}) & -[(\nu - \tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}})^2+\gamma^2+|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2] \\
\end{array}\right),
\label{eqa:A}\end{aligned}$$ which is obtained from the definition of Eq. (\[eq:A\]) with Eq. (\[eq:GR\]): $$G^{\mathrm{R}}_{{\textbf{\itshape{k}}}}(\nu) = \frac{1}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}\left(
\begin{array}{cc}
D^{*}_{{\textbf{\itshape{k}}}}(\nu)(\nu+\tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}}+\ii\gamma) & -D^{*}_{{\textbf{\itshape{k}}}}(\nu){\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \\
-D^{*}_{{\textbf{\itshape{k}}}}(\nu){\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}^{*} & D^{*}_{{\textbf{\itshape{k}}}}(\nu)(\nu-\tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}}+\ii\gamma) \\
\end{array}\right),
\label{eqa:GR}$$ where $$D_{{\textbf{\itshape{k}}}}(\nu) \equiv (\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}} + {E_{{\textbf{\itshape{k}}}}} + \ii\gamma)(\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}} - {E_{{\textbf{\itshape{k}}}}} + \ii\gamma).$$ The diagonal element $A_{11(22)}(\nu;{\textbf{\itshape{k}}})$ can then be described as Eq. (\[eq:A2\]). In the following, by using thse expressions, Eqs. (\[eqa:Delta\]) and (\[eqa:nehk\]) are derived from Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) with Eqs. (\[eq:RTA’\])-(\[eq:A\]).

First, we discuss $\nu$-integral forms of the population inversion $N_{{\textbf{\itshape{k}}}} \equiv \nek + \nhk -1$ and the polarization function $\pk$ because Eqs. (\[eqa:Delta\]) and (\[eqa:nehk\]) are described by the integration with respect to $\nu$. From Eqs. (\[eq:MSBE\_nehk\]) and (\[eq:RTA’\]) with $\partial_t n_{\mathrm{e(h)},{\textbf{\itshape{k}}}}=0$, we obtain $$N_{{\textbf{\itshape{k}}}} = -\frac{2}{\gamma}\operatorname{Im}[{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}\pk^{*}]
            + \int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi} \left\{ f^{\mathrm{B}}_{\mathrm{e}}(\nu)A_{11}(\nu;{\textbf{\itshape{k}}})
            -(1-f^{\mathrm{B}}_{\mathrm{h}}(-\nu))A_{22}(\nu;{\textbf{\itshape{k}}})   \right\}.
\label{eqa:PopInv1}$$ where Eq. (\[eqa:SumRule\]) is used. In a similar manner, from Eqs. (\[eq:MSBE\_pk\]) and (\[eq:RTA’\]) with $\partial_t \pk=0$, $$\begin{aligned}
\pk = &- \frac{{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}}{2(\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}} - \ii\gamma)}N_{{\textbf{\itshape{k}}}} \nonumber \\
      &+ \frac{\gamma}{\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}} - \ii\gamma} 
        \int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi}\left\{ [1-f^{\mathrm{B}}_{\mathrm{h}}(-\nu)]G^{\mathrm{R}}_{12}(\nu;{\textbf{\itshape{k}}})
        -f^{\mathrm{B}}_{\mathrm{e}}(\nu)G^{\mathrm{R}*}_{21}(\nu;{\textbf{\itshape{k}}}) \right\}.
\label{eqa:pk1}\end{aligned}$$ Therefore, with Eqs. (\[eqa:A\]) and (\[eqa:GR\]), substitution of Eq. (\[eqa:pk1\]) into Eq. (\[eqa:PopInv1\]) yields $$\begin{gathered}
N_{{\textbf{\itshape{k}}}}=\int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi}\frac{2\gamma}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}
            \left \{ f^{\mathrm{B}}_{\mathrm{e}}(\nu) [(\nu+\tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}})^2-|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2+\gamma^2] \right. \\
            \left. +[f^{\mathrm{B}}_{\mathrm{h}}(-\nu)-1] [(\nu-\tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}})^2-|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2+\gamma^2] \right \}.
\label{eqa:PopInv2}\end{gathered}$$ By substituting Eq. (\[eqa:PopInv2\]) into Eq. (\[eqa:pk1\]), we also find $$\pk={\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi} \frac{2\gamma}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}
            \left \{ [f^{\mathrm{B}}_{\mathrm{h}}(-\nu)-1](\nu - \tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}} - \ii\gamma)
            - f^{\mathrm{B}}_{\mathrm{e}}(\nu)(\nu + \tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}} + \ii\gamma)   \right \}.
\label{eqa:pk2}$$ Although the derivation of Eqs. (\[eqa:PopInv2\]) and (\[eqa:pk2\]) is straightforward, the following equations would be useful in the derivation: $$\begin{aligned}
\pm|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2\operatorname{Im}[D_{{\textbf{\itshape{k}}}}(\nu)(\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}}-\ii\gamma)]
 + &([\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}}]^2+\gamma^2)\operatorname{Im}[D^{*}_{{\textbf{\itshape{k}}}}(\nu)(\nu \mp \tilde\xi_{\mathrm{e(h)},{\textbf{\itshape{k}}}}+\ii\gamma)] \nonumber \\
 &= -\gamma({E_{{\textbf{\itshape{k}}}}}^2 + \gamma^2)((\nu \mp \tilde\xi_{\mathrm{e(h)},{\textbf{\itshape{k}}}})^2-|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2+\gamma^2), \\
(\nu+\tilde\xi_{\mathrm{h}})^2 - |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2 + \gamma^2 &-D_{{\textbf{\itshape{k}}}}(\nu) = 2 (\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}} - \ii\gamma)(\nu + \tilde\xi_{\mathrm{h}} + \ii\gamma), \\
(\nu-\tilde\xi_{\mathrm{e}})^2 - |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2 + \gamma^2 &-D^{*}_{{\textbf{\itshape{k}}}}(\nu) = -2 (\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}} - \ii\gamma)(\nu - \tilde\xi_{\mathrm{e}} - \ii\gamma). \end{aligned}$$ The $\nu$-integral forms of $N_{{\textbf{\itshape{k}}}}$ and $\pk$ are thus obtained as Eqs. (\[eqa:PopInv2\]) and (\[eqa:pk2\]), respectively. These expressions are helpful to find the $\nu$-integral forms of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$, $\nek$, and $\nhk$, which turn out to be the same as Eqs. (\[eqa:Delta\])-(\[eqa:Fh\]), as shown below.

From the definition of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} \equiv g^{*}a_0 + \sum_{{\textbf{\itshape{k}}}'}U'_{{\textbf{\itshape{k}}}-{\textbf{\itshape{k}}}'}p_{{\textbf{\itshape{k}}}'}$ and Eq. (\[eq:MSBE\_a0\]) with $\partial_t\a0 = 0$, ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$ can be described as $${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} = \sum_{{\textbf{\itshape{k}}}'} \left\{ \frac{|g|^2}{\xi_{\mathrm{ph},0}-\ii\kappa}+U'_{{\textbf{\itshape{k}}}'-{\textbf{\itshape{k}}}} \right\}p_{{\textbf{\itshape{k}}}'}
       = \sum_{{\textbf{\itshape{k}}}'} U_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}}^{\mathrm{eff},\kappa} p_{{\textbf{\itshape{k}}}'}.
\label{eqa:temp2}$$ Therefore, after the substitution of Eq. (\[eqa:pk2\]) into Eq. (\[eqa:temp2\]), we obtain $$\begin{gathered}
{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} = \sum_{{\textbf{\itshape{k}}}'} U_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}}^{\mathrm{eff},\kappa} \Delta_{{\textbf{\itshape{k}}}'} \int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}'}(\nu)
    \left \{ [F^{\mathrm{B}}_{\mathrm{h}}(\nu)-1][\nu - \tilde\xi_{\mathrm{e},{\textbf{\itshape{k}}}'} - \ii\gamma] \right. \\
    \left. + [F^{\mathrm{B}}_{\mathrm{e}}(\nu)-1][\nu + \tilde\xi_{\mathrm{h},{\textbf{\itshape{k}}}'} + \ii\gamma] \right\},\end{gathered}$$ where the definitions of Eqs. (\[eqa:L\])-(\[eqa:Fh\]) are used with $L_{{\textbf{\itshape{k}}}}(\nu) = \gamma/|D_{{\textbf{\itshape{k}}}}(\nu)|^2$. This equation can be rewritten as $$\begin{gathered}
{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}} = \sum_{{\textbf{\itshape{k}}}'} U^{\mathrm{eff},\kappa}_{{\textbf{\itshape{k}}}',{\textbf{\itshape{k}}}} {\mathnormal{\Delta}_{{\textbf{\itshape{k}}}{'}}} \int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}'}(\nu) 
            \left \{ (F_{\mathrm{e}}^{\mathrm{B}}(\nu) + F_{\mathrm{h}}^{\mathrm{B}}(\nu) - 2) (\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}'}) \right. \\
          + \left. (F_{\mathrm{e}}^{\mathrm{B}}(\nu) - F_{\mathrm{h}}^{\mathrm{B}}(\nu)) (\tilde\xi^{+}_{\mathrm{eh},{\textbf{\itshape{k}}}'} + \ii\gamma) \right\}.
\label{eqa:Delta2}\end{gathered}$$ By noting $\int \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}'}(\nu) (\nu-\tilde\xi^{-}_{\mathrm{eh},{\textbf{\itshape{k}}}'}) = 0$, we thus find that Eq. (\[eqa:Delta2\]) is equivalent to Eq. (\[eqa:Delta\]).

Our remaining task is now to derive the $\nu$-integral forms of $\nek$ and $\nhk$. By multiplying ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$ by the complex conjugate of Eq. (\[eqa:pk2\]), $$\frac{1}{\gamma}\operatorname{Im}[{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}\pk^{*}] = |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2 \int^{\infty}_{-\infty}\frac{\dd\nu}{2\pi} \frac{2\gamma}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2} \{ f^{\mathrm{B}}_{\mathrm{e}}(\nu) + f^{\mathrm{B}}_{\mathrm{h}}(-\nu) -1 \},
\label{eqa:temp1}$$ can be obtained. The substitution of Eq. (\[eqa:temp1\]) into Eq. (\[eq:MSBE\_nehk\]) with $\partial_t n_{\mathrm{e(h)},{\textbf{\itshape{k}}}} = 0$, then, yields $$\begin{gathered}
n_{\mathrm{e(h)},{\textbf{\itshape{k}}}} = \int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} \frac{2\gamma}{|D_{{\textbf{\itshape{k}}}}(\nu)|^2}
    \left \{ [(\nu \pm \tilde\xi_{\mathrm{h(e)},{\textbf{\itshape{k}}}})^2 + \gamma^2]f^{\mathrm{B}}_{\mathrm{e(h)}}(\pm\nu) \right. \\
    \left. -|{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2[ f^{\mathrm{B}}_{\mathrm{h(e)}}(\mp\nu) -1] \right\},\end{gathered}$$ which can be rewritten as $$\begin{gathered}
n_{\mathrm{e(h)},{\textbf{\itshape{k}}}} = \int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}}(\nu)
    \left \{ [(\nu \pm \tilde\xi_{\mathrm{h(e)},{\textbf{\itshape{k}}}})^2 + \gamma^2 + |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2] \right. \\
    \left. \mp [(\nu \pm \tilde\xi_{\mathrm{h(e)},{\textbf{\itshape{k}}}})^2 + \gamma^2]F^{\mathrm{B}}_{\mathrm{e(h)}}(\nu) \mp |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2F^{\mathrm{B}}_{\mathrm{h(e)}}(\nu) \right\}.
\label{eqa:nehk2}\end{gathered}$$ We then find that Eq. (\[eqa:nehk2\]) is identical to Eq. (\[eqa:nehk\]) because, from Eqs. (\[eqa:SumRule\]) and (\[eqa:A\]), $$\int^{\infty}_{-\infty} \frac{\dd\nu}{2\pi} L_{{\textbf{\itshape{k}}}}(\nu)[(\nu \pm \tilde\xi_{\mathrm{h(e)},{\textbf{\itshape{k}}}})^2 + \gamma^2 + |{\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}|^2] = \frac{1}{2}.$$ Thus, we have shown that Eqs. (\[eqa:Delta\]) and (\[eqa:nehk\]) are derived from Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) with Eqs. (\[eq:RTA’\])-(\[eq:A\]). This means that the formalism in Ref. [@Yamaguchi12] is equivalent to Eqs. (\[eq:MSBE\_a0\])-(\[eq:MSBE\_nehk\]) with Eqs. (\[eq:RTA’\])-(\[eq:A\]).

[^1]: In this contribution, the terms ‘lasing’ and ‘laser’ are used only when the condensation is inherently governed by non-equilibrium kinetics, according to [@Bajoni08; @Deng03]. In other words, thermodynamic variables of the system, such as temperatures, cannot be defined for lasing phases. However, we note that these terms are occasionally used even for a condensation dominated by the thermodynamics of the electron-hole-photon system [@Kasprzak08] if the interest is in fabricating a device.

[^2]: We note that theories for describing dynamics of equilibrium phases, e.g. the Gross-Pitaevskii equation, can asymptotically be derived from Maxwell-Bloch equations [@Berloff13] even though these theories still do not recover equilibrium statistical approaches.

[^3]: We note, however, that physical guesses are required for the determination of what operator(s) should be chosen as $\hat{O}$, e.g. from experiments.

[^4]: Here, $-1 \le N_{{\textbf{\itshape{k}}}} \le +1$ because $0 \le n_{\mathrm{e(h)},{\textbf{\itshape{k}}}} \le 1$. Population inversion is formed in ${\textbf{\itshape{k}}}$-resions with $N_{{\textbf{\itshape{k}}}}>0$.

[^5]: For $\kappa = \gamma = 0$, the steady state of the MSBE becomes identical to the BCS theory if the solution of the BCS theory is chosen as an initial condition in the MSBE because $[\hat{\rho}_{\mathrm{eq}}^{\mathrm{MF}},\hat{H}_{\mathrm{S}}^{\mathrm{MF}}]=0$. However, this is a special case.

[^6]: The retarded Green’s funcition is also seen as a kind of density of states.

[^7]: Under the condition (I), $f^{\mathrm{B}}_{\mathrm{e(h)}}(\nu)$ in Eq. (\[eq:RTA’\]) can be approximated by the values at $\nu = \pm {E_{{\textbf{\itshape{k}}}}}$ because $A_{\alpha\alpha'}(\nu;{\textbf{\itshape{k}}})$ and $G^{\mathrm{R}}_{\alpha\alpha',{\textbf{\itshape{k}}}}(\nu)$ have peaks around $\nu = \pm {E_{{\textbf{\itshape{k}}}}}$, as seen in Eqs. (\[eq:A2\]) and (\[eqa:GR\]).

[^8]: This means that the system is chemically non-equilibrium with the pumping baths even if the system is in quasi-equilibrium.

[^9]: In the numerical calculations, the ${\textbf{\itshape{k}}}$-dependence of ${\mathnormal{\Delta}_{{\textbf{\itshape{k}}}}}$ is eliminated by using a contact potential $U'_{{\textbf{\itshape{q}}} \neq 0} = U = 2.66 \times 10^{-10}$ eV with cut-off wave number $k_{\mathrm{c}} = 1.36 \times 10^{9}$ $m^{-1}$. The other parameters are $m_{\mathrm{e}} = m_{\mathrm{h}} = 0.068m_0$ ($m_0$ is the free electron mass), $\mu^{\mathrm{B}}_{\mathrm{e}} = \mu^{\mathrm{B}}_{\mathrm{h}}$, $T = 10 $K, $g = 6.29 \times 10^{-7}$eV, $\gamma = 4$meV, and $\kappa = 100 \mathrm{\mu eV}$. In this context, our calculations are not quantitative but qualitative even though the parameters are taken as realistic as possible. In this situation, the exciton level ($\equiv E_{\mathrm{ex}}$) is formed at 10 meV below $E_{\mathrm{g}}$ ($E_{\mathrm{ex}} = E_{\mathrm{g}} - 10$ meV) and the lower polariton level $E_{\mathrm{LP}}$ is created at 20 meV below $E_{\mathrm{g}}$ ($E_{\mathrm{LP}} = E_{\mathrm{g}} - 20$ meV) under the resonant condition $E_{\mathrm{cav}} = E_{\mathrm{ex}}$ [@Yamaguchi13].

[^10]: For $E_{\mathrm{LP}}$ and $E_{\mathrm{ex}}$, see also Appendix I. Excitonic effects are discussed in the low density limit.

[^11]: The origin of the gap is analogous to the Rabi splitting in resonance fluorescence [@Scully97; @Schmitt-Rink88; @Henneberger92; @Yamaguchi13].

[^12]: In fact, the second threshold and the blue-shift can also be caused by a different mechanism even if the system remains in quasi-equilibrium. In this situation, however, the gain spectrum shows only absorption [@Yamaguchi12].

[^13]: Discussion in Appendix I is reproduced from the supplemental material in Ref. [@Yamaguchi13].
