BEC–BCS-laser crossover in Coulomb-correlated electron–hole–photon systems

Many-body features caused by Coulomb correlations are of great importance for understanding phenomena pertaining to polariton systems in semiconductor microcavities, i.e. electron–hole–photon systems. Remarkable many-body effects are shown to exist in both thermal-equilibrium phases and non-equilibrium lasing states. We then show a unified framework for connecting the thermal-equilibrium and the non-equilibrium steady states based on a non-equilibrium Green's function approach. Bose–Einstein condensate (BEC)–Bardeen–Cooper–Schrieffer (BCS)-laser crossovers are investigated by using this approach.


Introduction
Semiconductor microcavity systems containing quantum wells (QWs) have attracted considerable attention in recent years as test beds for fundamental physics, as well as from the viewpoint of their potential applications [1,2]. One of the most exciting phenomena observed in this electron-hole-photon (e-h-p) system is a Bose-Einstein condensate (BEC) of an exciton-polariton gas, i.e. a macroscopic occupation of a single exciton-polariton state by a thermodynamic phase transition. Here, exciton-polaritons are composite quasi-bosonic particles consisting of excitons (electron-hole (e-h) bound pairs formed by the Coulomb interactions in a low-excitation regime) and microcavity photons coupled by their light-matter interactions [3,4]. Therefore, many fundamental features of exciton-polaritons arise from their original constituents, i.e. electrons, holes and photons.
In particular, features arising from photons play two important roles in the realization of the BEC. Firstly, the cavity photon mass is roughly four orders of magnitude smaller than that of the bare exciton. Therefore, the effective mass of the exciton-polariton also decreases to the same degree; this lends notable experimental advantages, such as a high critical temperature and low critical density. Secondly, the fast leakage of photons from the cavity results in the relatively short lifetime of the exciton-polariton state (∼10 −12 s). This is not an advantage but rather a disadvantage for the BEC because the system can be in thermal equilibrium 6 only when the 4 lifetime [5][6][7]. One typical example is the normal photon lasing 8 induced by stimulated emission of the cavity photons from the e-h plasma [29,30]: at the laser frequency, the stimulated emission becomes very fast and many carriers are expended. As a result, a dip is formed in the carrier distribution (kinetic hole burning) [29,30] because the thermalization processes can no longer supply the lost carriers at sufficient speed. This means that non-equilibrium features cannot be negligible at all.
In such a situation, the lasing operation can be well described by the Maxwellsemiconductor-Bloch equations (MSBE)-these are coupled kinetic equations with dampinginduced and thermalization-induced terms [30][31][32]. This approach explicitly treats the electrons and holes where the Coulomb interactions are included within the Hartree-Fock (HF) approximation. However, the applicable range of carrier densities is limited to high-density regimes because it is assumed that the actual distribution is driven to approach the Fermi distribution under the relaxation time approximation [29,30]. More importantly, MSBE cannot reproduce the results under thermal equilibrium even though the damping term is set to zero, and a sufficiently small but finite thermalization rate is assumed. This means that we cannot discuss the changes from the above-described thermal equilibrium results into non-equilibrium regimes because the respective theories assume incompatible situations. In this sense, the development of a more comprehensive theory now has increasing importance for understanding the e-h-p systems [33].
In this study, for the systematic investigation of the e-h-p system, we first explain the treatment of the condensed phase under thermal equilibrium. Then, conventional photon lasers are analyzed in the framework of the MSBE with the assumption that the system is in a non-equilibrium steady state. Although some aspects will overlap with our previous studies [19][20][21]30], we discuss the similarities and differences between the thermal equilibrium and the non-equilibrium steady state in particular, and several new results are shown. Thus, the condensed phases in thermal equilibrium and the photon lasing in a non-equilibrium steady state are separately analyzed, and their features are discussed in detail. Then, we present a new framework that can explain their intermediate regimes. This framework is an extension of a nonequilibrium Green's function approach developed in [34][35][36] in which the e-h system is simply modeled by two-level systems (disorder-localized excitons rather than electrons and holes in a clean system); we include an appropriate e-h picture as well as the Coulomb interactions within the HF approximation. Then, it is shown that the obtained coupled equations can successfully reproduce the thermal-equilibrium results of the condensed phases in all density regimes. Furthermore, the equations can also recover the results obtained from the MSBE when nonequilibrium features become important. This means that we can simultaneously discuss the thermal-equilibrium and non-equilibrium steady-state phenomena by using only one theory. The physical meanings of the conditions recovering the equilibrium and non-equilibrium results are clearly addressed, where a flux of particles has the effect of breaking e-h pairs. Here, it 8 In general, the laser physics cannot be described without non-equilibrium parameters because the steady state is determined by the balance of the pumping and loss. Furthermore, thermodynamic variables (the effective chemical potentials and temperatures) cannot be defined for the system and/or subsystems. This is in contrast to the quasithermal equilibrium 6 , and therefore, the laser physics cannot be described by the thermal equilibrium theory. In the polariton community, the terms 'lasing' and 'laser' are occasionally used even for a condensation dominated by the thermodynamics of the e-h-p system [6] if the interest is in fabricating a device, because the system is not in the true thermal equilibrium due to the pumping and loss. However, in this paper, these terms are used only when the condensation is inherently governed by non-equilibrium kinetics, according to [7,33]. 5 is interesting to note that the ability to simultaneously describe thermal-equilibrium and nonequilibrium steady-state (lasing) features is similar to that in the works [34][35][36] even though the two models, i.e. disorder-localized excitons and electrons and holes with Coulomb interactions, are different from each other. In this sense, there is a model-independent aspect in this type of physics. Finally, numerical results clarify the crossovers between the lasing and the thermal condensed phases (exciton BEC, e-h BCS, and exciton-polariton BEC and photonic BEC). However, unlike the normal photon lasing, it is shown that the obtained lasing is not described solely by the MSBE but is still influenced by the BCS gap equation, the features of which cannot be treated by the two-level models [34][35][36].

Thermal polariton condensates and kinetic semiconductor lasers
Thermal polariton condensates-equilibrium ground states and kinetic semiconductor lasers-non-equilibrium stationary states are similar phenomena because a macroscopic phase coherence of the polarization and the photon field is built up when the carrier density is increased above a threshold. However, the physical processes by which they reach the stationary states are significantly different, i.e. one realizes the minimization of the total free energy, whereas the other realizes the valance between the gain and the loss. As shown in a later section, these two different processes can be described in a unified theoretical framework. In this section, we first show two different theories for describing the stationary states of the thermal polariton condensates and kinetic semiconductor lasers. First, thermal polariton condensates are studied by the BCS-type theory for T = 0 in section 2.1. Then, kinetic semiconductor lasing is described by the MSBE in section 2.2. The similarities and differences between the two are discussed in detail.

Bardeen-Cooper-Schrieffer (BCS)-type theory for the ground states of thermal polariton condensates
It is well known that a large number of particles created in the e-h-p system can condense into a single ground state thermodynamically if the particle lifetime is longer than the thermalization time [5][6][7]. Hence, the e-h-p system can be approximately regarded as a closed system when such a quasi-thermal equilibrium is achieved. In this subsection, we study the ground state properties obtained from the BCS-type variational approach at T = 0 [19][20][21] Here, ε e/h,k = k 2 /2m e/h + (E g − µ)/2, ω c , µ, and g are the electronic dispersion within an effective mass approximation, cavity frequency, chemical potential, and dipole coupling strength, respectively.ê k ,ĥ k andâ denote the annihilation operators of electrons, holes with momentum k, and cavity photons, respectively. In equation (2), the two colons indicate the normal ordering of the creation and annihilation operators. U q is the Coulomb interaction, and it is defined as Here, the elimination of the q = 0 term arises from the Coulomb interaction with a uniform positive charge background [32]. In this section,h = k B = 1 is used for simplicity. The ground state of the e-h-p system is realized at T = 0, and the e-h pairs and cavity photons are in a coherent state when we assume that the ground state is condensed. Such a coherent state can be described by a simple extension of a previous theory for excitonic insulators [17] to the e-h-p system. The BCS-type ground state is given by where u k , v k and λ are the variational parameters (which are real numbers [19]) carrying microscopic information about the polariton condensates [19-21, 37, 38]. For example, the expected value of the microscopic polarization, carrier distribution and amplitude of the cavity mode are given by , and â = λ, respectively. To consider the fermionic statistics of carriers, a constraint u 2 k + v 2 k = 1 should be imposed on the variational parameters.
The variational parameters in equation (6) can be determined by the given number of excitations N tot and variational equations for the mean-field free energy F MF = E MF + µN tot : where E MF is the mean-field energy defined as As a result, the following equations can be obtained: where ε k ≡ k 2 /(2m r ) + E g . In section 2.1.2, we first study these equations in the limit of low excitation density, i.e. v k 1, in order to understand their physical meanings. In section 2.1.3, numerical calculations are then shown to discuss the regimes including the limit of high excitation density. 2.1.2. Analytical forms in the low-density limit. If the excitation density is sufficiently low (v k 1), v 2 k ∼ 0 and u k ∼ 1 can be applied in equations (10) and (11), and a linearized set of equations in terms of small v k is obtained: This linearization corresponds to neglecting the fermionic nature of the carriers, i.e. the nonlinear saturation effects due to the Pauli exclusion principle (i.e. u 2 k + v 2 k = 1 within meanfield theory). Then, one obtains the equation where is the strength of the effective attraction between electrons and holes. Equation (14) is formally equivalent to the Schrödinger equation for the e-h relative motion. Then, we can substitute the microscopic polarization v k as v k = ηφ k , where φ k is the bound-state wavefunction of an e-h pair and η 2 = N eh is the number of e-h pairs (which is equal to the number of electrons at low-carrier density) [17,18]. Hence, equation (14) clearly shows that an e-h bound pair can be formed not only by Coulomb interactions but also by photon-mediated attraction.
Here, it would be instructive to consider one limiting case in which the wavefunction of the e-h bound pair is determined mainly by a Coulombic origin, i.e. U k−k g 2 /(ω c − µ). In such a situation, φ k reduces to the wavefunction of a (1s) exciton (≡ φ X,k ), which satisfies equation (14) with the substitution µ → ω X . Here, ω X is the energy level of the 1s exciton. Therefore, by substituting v k = ηφ X,k into equations (12) and (13) with k φ 2 X,k = 1, we obtain where the coupling parameter is scaled as g ≡ φ X (0)g with φ X (0) ≡ k φ X,k . The chemical potential µ is then given by one of the eigenvalues of these two coupled equations, which are the eigenenergies of the upper and lower polaritons: In our case, the solution µ = ω LP should be selected because the ground state is now considered and the situation corresponds to the BEC of e-h-p composite particles into a (zero-momentum) ground state of lower polariton branches. The selection of the lowest energy from several candidates of µ corresponds to energy minimization based on the thermodynamics at T = 0: the thermalized state at T = 0 should be in the ground state. Thus, a familiar result can be obtained in the low-density limit. Here, it should be noted that ω LP and ω UP in equation (17) are the well-known expressions obtained when excitons are treated as simple bosons [32]; however, such a treatment is correct only when the photon-mediated attraction is not a dominant factor for the determination of the wavefunction in equation (14), as described above.

Numerical results.
In order to discuss the variational problem even for the higherdensity regimes, we need to numerically evaluate the variational equations (9)-(11) while maintaining the higher-order terms in v k because the fermionic nature of carriers arises. For later reference, we note that equations (10) and (11) can be combined into one equation by using u 2 k + v 2 k = 1: where Equation (18) is equivalent to the BCS gap equation for T = 0, where the attraction force between electrons and holes is again given by U eff k−k . This means that the photon-mediated attraction included in U eff k−k is always present at arbitrary density. Here, k in equation (18) is equivalent to the composite order parameter commonly used in a Feshbach resonance of a uniform Fermi atom gas [39,40]. The second term in k corresponds to the Coulomb enhancement, and BGR k comes from the renormalization of the single-particle energies, in the context of the band-gap renormalization [32].
In our previous papers [19,21], the variational problem was solved for the Coulomb interaction model with U q = 4 e 2 /(εq V ) (V is the quantization volume of the system) by an interpolation approximation method and it was shown that various types of condensed ground states can arise depending on the detuning parameter d ≡ ω c − E g and the total excitation density n tot ≡ N tot /V (the results are qualitatively the same as those without the approximation [20], which indicates that the interpolation method works well with allowing a significant reduction of numerical efforts). Here, we present the results for another model of the interacting carriers with a contact interaction U q = U (the contact potential model mimics the screening effect in the high-carrier-density regime [38]) in order to discuss qualitative behaviors. The results shown here are qualitatively the same as the previous results for the real Coulomb potential [19][20][21]. In the numerical analysis, we consider that the carriers are in a QW, i.e. in a two-dimensional (2D) system. For the contact potential in 2D systems, we have to introduce a high-energy cutoff (E c ) for the convergence of the mean-field energy; otherwise the bound-state energy diverges. We choose the coupling strength g, U and E c in order to achieve ω LP = E g − 2|E x | and ω X = E g − |E x | for d = −|E x |. This means that the excitonbinding energy (= |E x | ≡ 1 Ry) is equal to half of the vacuum Rabi splitting (= [ω UP − ω LP ]/2) under the resonance condition d = −|E x |.
The chemical potential µ is shown in figure 1 as a function of the total excitation density n tot for different values of the detuning parameter d. First, focusing on the cases of large blue detuning (purple and green lines in figure 1), the ground state of the e-h-p system should be the exciton BEC with no photonic excitations in the low-density limit. Hence, in this case, the chemical potential is found at the exciton level ω X . As the excitation density increases, the ground state changes to form a high-density coherent state, i.e. an e-h BCS state in which a large number of carriers form the Fermi sea while a macroscopic coherent polarization is still present (BEC-BCS crossover [17,18]). With a further increase in the excitation density, the Fermi energy reaches the cavity mode energy µ ∼ = ω c . Thus, the chemical potential increases monotonically with density and eventually saturates at µ = ω c in the high-density limit. This phenomenon is naturally understood from the principle of energy minimization by considering the fermionic nature of carriers and bosonic nature of photons: the carriers occupies the higherenergy (momentum) states to form the Fermi sea, while there is no upper limitation for the occupation number of a photon mode. As a result, it is favorable to increase the photon numbers rather than electrons and holes (photonic BEC [19][20][21]).
These interpretations can be directly supported by the plots of the carrier density n car (≡ k v 2 k /V ) and photon densities n ph (≡ λ 2 /V ) as a function of the total excitation density in figure 2(a) (purple and green lines for large blue detuning). The photon density starts increasing nonlinearly when n tot is 0.1-1 Bohr −2 of the value achieved when µ ∼ ω c in figure 1. Subsequently, only the photon density increases but the carrier density does not. Hence, the photon fraction (≡ n ph /n tot ) dominates the excitation density in the higher-density regime ( figure 2(b)). These features are fully consistent with the above explanations. Thus, the behaviors of µ, n car , and n ph are clearly understood in the cases of large blue detuning. Similar explanations can also be applied to the resonant case (ω c ∼ = ω X ). The difference due to the detuning is found in the chemical potential in the low-density limit: µ is equal to the lower polariton level ω LP for ω c ∼ = ω X (red and blue lines in figure 1). As a result, the coherent ground states of e-h-p systems change from exciton-polariton/exciton BEC to the photonic BEC as the excitation density increases (figure 2(b)). Depending on the detuning parameter, the change of the ground states to photonic condensation occurs smoothly or rapidly (figures 2(a) and (b)). The momentum distribution function of electrons n k is plotted in figure 3 for (a) the resonant case (d = −1Ry) and (b) large blue detuning (d = 10 Ry). Three curves in each figure are obtained for different densities n tot (the parameters correspond to the blank squares in figure 2(b)) where the ground states ranged from/to the exciton BEC, e-h BCS, exciton-polariton BEC, and photonic BEC. In the cases of exciton-polariton BEC (n tot = 0.1, 1 Bohr −2 in figure 3(a)) and exciton BEC (n tot = 0.1 Bohr −2 in figure 3(b)), the half-width of the distribution in k-space is larger in the resonant case (d = −1 Ry) than in the large blue detuning case (d = 10 Ry). This means that the e-h pair is tightly bound in real space by the Coulomb attraction and photon-mediated attraction in the resonant case because the shape of the distribution shows the wavefunction of e-h pairs through the relationship v 2 k = −η 2 φ 2 k at low densities (see also section 2.1.2). Then, at increased density for the large blue detuning case (n tot = 1 Bohr −2 in figure 3(b)), the shape of the distribution changes into the Fermi-Dirac function rounded at the Fermi surface. This change is typical in the context of the BEC-BCS crossover [17,18]. However, when the photonic fraction becomes large (n tot = 200 Bohr −2 in figures 3(a) and (b); see also figure 2(b)), the width of the distribution function becomes very large, implying the strong enhancement of the photon-mediated attraction, because in the presence of a large and coherent photon field, the dipole-coupling energy dominates the Coulomb interaction energy in E MF in equation (8). In the limit λ → +∞, the variational problem gives the solution u k v k → 1/2 and v 2 k → 1/2. Therefore, the width of the distribution is broadened in order to approach v 2 k = 1/2 in the high-density limit of the photonic BEC. Thus, we have discussed the behaviors of µ, n car , n ph , and n k and shown that the photonic nature plays an important role in the high-density regime because the photon-mediated attraction becomes actualized. Here, one can easily expect that the pair-breaking energy of an e-h pair would be strongly enhanced by this attraction force. Therefore, we finally show the minimal pair-breaking energy of an e-h pair (= min[2E k ], E k ≡ ξ 2 k + 2 k ) in figure 4 as a function of the carrier density n car , where the black line (d = +∞) is drawn as a reference for the absence of photon-mediated attraction. At low density, the ground states are the exciton-polariton BEC for the resonant case (d = −1 Ry) and the exciton BEC for large blue detuning. In the former case, the pair-breaking energy is determined by both the lower polariton energy and the excitonbinding energy (2 Ry dashed line in figure 4) because of the inherent photon-mediated attractive force. In the latter case, it is merely the exciton-binding energy |E x | (1 Ry dashed line in figure 4). At high carrier density, one can find a strong enhancement of the pair-breaking energy compared with the black line (absence of photon-mediated attraction), which implies that the photon-mediated attraction dominates U eff k−k . Here, it should be noted that the diverging behavior of the pair-breaking energy is quite natural because µ becomes close to ω c and U eff k−k diverges within the framework of the thermal equilibrium theory, which is consistent with the previous works [20,21,41]. As a result, it was shown that the e-h pairing becomes robust in the photonic BEC regime.
In summary, the coherent ground states of the e-h-p system are theoretically determined by the minimization of the free energy. The ground-state property varies with the excitation density and the detuning of the cavity mode. The variation is caused by the change in the nature of the particles, i.e. it depends on whether the particles are bosonic or fermionic and whether the particles are exciton-like or photon-like.

Maxwell-semiconductor-Bloch equations (MSBE) for kinetic semiconductor lasers
In the previous subsection, we have described the case in which the e-h-p system can be considered to be in thermal equilibrium. In contrast, the e-h-p system can also work as a semiconductor laser when electrons and holes are highly excited from an external source, and the coherent light is emitted into free space. Hence, in such a case, it cannot be regarded as a closed system but as an open system where the total energy is not conserved inside the e-h-psystem(seefootnote8).Thenon-equilibriumnaturecanbesimplydescribedbythe MSBE, which is a set of kinetic equations of motion for the cavity photon field λ, microscopic polarization p k , and population inversion N k = n e,k + n h,k − 1(n e,k ≡ e † k e k , n h,k ≡ h † −k h −k ), and not by the variational equations for the total free energy: Here, λ and p k are complex numbers unlike in the BCS theory. k and BGR k are the renormalized Rabi frequency and the band-gap renormalization, respectively, the definitions of which are the same as in equations (20) and (21) [32]. The terms κ, γ ⊥ , and γ represent the cavity loss, dephasing of the polarization and thermalization of the population inversion toward N 0 k , respectively. These terms are introduced phenomenologically within the relaxation time approximation [29,30]. In this approximation, it is assumed that the carriers are first excited following which intraband relaxations immediately thermalize the carriers (typically ∼100 fs) by emitting phonons and/or by scattering other electrons and holes. Hence, the time scale of the thermalization should be much shorter than the dynamics of the e-h-p system. As a result, the pumping processes can be effectively considered by the thermalization rate γ with a thermalized population inversion N 0 k in equation (24). It would be natural to use where f e,k and f h,k are the Fermi distribution functions when assuming 13 that the system thermalizes into the e-h plasma phase. In this paper, we focus on the steady state determined by the MSBE.
In the case of a single-mode laser, the cavity photon field λ and the microscopic polarization p k oscillate with a laser frequency L that is determined as given below. Therefore, in the nonequilibrium steady state, the following conditions should be satisfied: At the same time, the distribution function is time independent, In section 2.2.1, we first discuss the similarities and differences between the steady-state MSBE and the BCS-type theory at T = 0. Then, in section 2.2.2, we present numerical calculations.

Similarities and differences between MSBE and BCS
-type theory at T = 0. We describe the similarities and differences between the MSBE and the BCS-type theory at T = 0. For this purpose, we first discuss how the steady state is determined, especially in terms of the optical susceptibility χ ( ), which gives a clear picture of the steady-state MSBE and BCS-type theory at T = 0. The susceptibility χ( ) is defined by the ratio of macroscopic polarization If one knows the stationary distribution N k , χ( ) is determined by the polarization equation in equation (23), which is the linear equation of p k and λ.
Using the definition of the susceptibility, we can rewrite equation (22) with equations (25) and (26) as These are the two balance equations that should be satisfied for the laser oscillation to be realized. In equation (28), the frequency of the cavity photons is modified by −Re[χ( L )] as a result of the refractive index change induced by the population inversion N k . Hence, equation (28) gives the frequency-matching condition between the induced polarization and the cavity photons under the frequency shift. In equation (29), G( L ) corresponds to the optical gain due to the population inversion N k > 0 at a certain k. In the stationary operation of lasers, the optical loss rate κ from the cavity is compensated for by the optical gain G( L ). Therefore, once a stationary distribution function N k is given, several candidates for the laser frequency L are found by the frequency-matching condition (equation (28)), and then, one of them is selected by the balance of the gain and loss conditions (equation (29)).
Such a picture can also be applied to the BCS-type theory at T = 0, described in section 2.1. From the real part of equation (10), one can obtain a frequency-matching condition identical to equation (29) except that L → µ. This is because the cavity photon field should be in resonance with the coherent polarization in the condensed phase. However, from the imaginary part of equation (10), no information can be obtained because the variational parameters are real numbers. This means that the state is not determined by the balance of the gain and loss in the case of a closed system. Instead, the lowest energy should be selected among the numbers 14 of candidates of µ satisfying the frequency-matching condition (see also section 2.1.2). This is because the fully thermalized state at T = 0 should be in the ground state, which is an important difference between the MSBE and the BCS-type theory at T = 0. Now, one should notice that equations (22) and (23) in the MSBE under the conditions of equations (25) and (26) have forms similar to those of equations (10) and (11) in the variational problems of the closed system. Actually, equations (10) and (11) in the variational problem can be obtained from equations (22), (23), (25) and (26) in the MSBE with κ = γ ⊥ = γ = 0 and L → µ. This is quite natural in the sense that the MSBE with κ = γ ⊥ = γ = 0 is equivalent to the Heisenberg equations of motion for a closed system, and equations (22) and (23) are equivalent to the eigenvalue problem under the conditions of equations (25) and (26). Hence, the BCS-type theory at T = 0 is formally recovered for a given total number of excitations N tot = k (N k + 1)/2 + |λ| 2 (corresponding to equation (9)) if (i) we additionally use a condition for the square length of a Bloch vector at momentum k, and if (ii) the ground state is subsequently selected from the candidates of the eigenstates by the energy minimization principle at T = 0. As shown below, operations (i) and (ii) are equivalent to introducing thermalization processes manually after neglecting the corresponding terms in the MSBE (the first term in equation (24)).
With regard to operation (i), we can see the effect of γ ⊥ and γ in s 2 k , which becomes less than that for the thermal-equilibrium cases. The equation of motion for s 2 k can be easily obtained from equations (23) and (24) as The right-hand side of equation (30) is negative for the Bloch sphere (s 2 k = 4| p k | 2 + N 2 k = 1) for −1 < N 0 k < 1 and 2γ ⊥ > γ ; these two inequalities are always fulfilled in physical systems (2γ ⊥ = 2γ > γ in our model presented in section 3.3). Hence, s 2 k becomes less than one by the effect of γ ⊥ and γ even if the system is initially prepared in a pure state (an eigenstate of the Hamiltonian; s 2 k = 1). Therefore, a non-equilibrium steady state for ds 2 k /dt = 0 is achieved inside the Bloch sphere. In contrast, in the BCS-type theory at T = 0, s 2 k = 1 is preserved by the normalization condition u 2 k + v 2 k = 1 (ds 2 k /dt = 0 with γ ⊥ = γ || = 0 in equation (30)) because the coherent state in equation (6) is an eigenstate of the Hamiltonian. Here, it should be noted that ds 2 k /dt = 0 can be obtained from the MSBE with κ = γ ⊥ = γ = 0, but the condition s 2 k = 1 cannot be obtained.
On the other hand, operation (ii) indicates that an energy minimization principle based on the thermodynamics at T = 0 is separately required in order to reproduce the BCS-type theory at T = 0 by using the MSBE with κ = γ ⊥ = γ = 0 and L → µ. This means that the thermalization processes described by γ in the MSBE are neglected once (γ = 0), and then, the thermalization is separately introduced manually. In this sense, one can notice that the additional condition s 2 k = 1 arises for the same reason. Hence, the MSBE should include the BCStype theory if the thermalization processes are correctly included. However, without the two artificial operations, the MSBE described in equations (22)-(24) cannot reproduce the results of BCS theory even though sufficiently small but finite γ and γ ⊥ are assumed with κ = 0. The problem would originate from the limitation of a sufficiently small γ (and γ ⊥ ), which conflicts with the initial assumption: the thermalization time (∼1/γ ) should be much shorter than the lifetime of the e-h-p system dynamics. This is similar to the difficulty of recovering thermalequilibrium results when describing interacting composite systems by the Markovian quantum master equation [42]. In section 3, we present a (non-Markov) unified theory for connecting the MSBE and the BCS-type theory in a true sense.

Numerical results.
In the remainder of this section, we present the numerical results of the stationary solution to the MSBE. The problem was presented previously for the Coulomb interaction model U q = 4 e 2 /( q V ) [30]. Here, the calculations are performed by a contact interaction model U q = U and equal electron and hole masses (therefore n e,k = n h,k = n k ) for a simple consideration of the many-body effects. In such a situation, the renormalized Rabi frequency is independent of the momentum: k = gλ + U k p k = . As a result, equations (22)-(24) can be a closed set of equations for macroscopic parameters: (we can choose it as real number without loss of generality), λ, L and n e ≡ k n e,k = k (N k + 1)/2. In the calculation, we assumed a positively detuned cavity mode ω c − E g = 2|E X | and κ = γ ⊥ = γ = 0.2|E X |.
In figure 5, the calculated photon density n ph = |λ| 2 and the electron density consumption δn ≡ n 0 e − n e (n 0 e ≡ k f e,k = k (N 0 k + 1)/2) are indicated by red and blue lines as a function of the injected electron density n 0 e at low temperature (panel (a)) T = 0.5|E x | and at high temperature (panel (b)) T = 3|E x |. In both figures, there exist threshold values of n 0 e above which the photon density increases. Similar curves are also obtained for the electron density consumption δn. This is easily understood by the following relationship obtained from equations (22) and (24) for the contact potential model: The threshold value at low temperature, ∼ = 0.08 Bohr −2 ( figure 5(a)), is much less than that at high temperature, ∼ = 0.8 Bohr −2 ( figure 5(b)). Moreover, the slope of the photon density curve is steeper during low-temperature operation (figure 5(a)) than during high-temperature operation ( figure 5(b)). These results indicate that the laser performance is highly efficient during low temperature operation, which is well known in the field of lasers [43].
Once the macroscopic parameters , λ, L and n e are obtained, the microscopic population n e,k and polarization | p k | are determined by equations (23) and (24). The results are shown in figure 6, where panels (A)-(E) show the low-temperature case with varying electron densities (from (A) the threshold density, ∼ = 0.08 Bohr −2 , up to (E) 0.38 Bohr −2 ), whereas panels (F)-(J) show the high-temperature case with varying densities (from (F) the threshold density, ∼ = 0.8 Bohr −2 , up to (J) 1.3 Bohr −2 ). In both plots, at the threshold density (panels (A) and (F) in figure 6), the distribution function coincides with the Fermi distribution function, and the polarization is negligible. As the density increases, the dip in the distribution function grows ((A) → (E) and (F) → (J)). The appearance of the dip in the population distribution is due to kinetic hole burning, which is observed above the threshold density and is caused by the stimulated recombination emission of the inverted carriers into the cavity mode. Because the carriers are strongly coupled with the lasing mode at the momentum k hole near the dip (kinetic hole), a large polarization is also found near k hole 9 . This feature is found irrespective of the temperature. The difference between low-and high-temperature operations lies in whether the burned hole and the polarization function are broad in k-space. At low temperature, we found that the polarization function is spread widely over k-space, indicating the strong (attractive) interaction between different e-h pairs with different momentums. This in turn means that BCS-like tightly bound e-h pairs (in real space) contribute to the low-temperature laser gain and correspond to a signature of non-equilibrium phase transition [44]. On the other hand, at high temperature, hole burning is observed only over a tiny area in k-space. This means that the e-h pairs related to the laser are almost unbound and are like e-h plasma.
Finally, we comment on the theoretical treatment of the optical gains caused by the strongly correlated e-h systems. As shown above, the methods of determining the ground state in closed systems and the non-equilibrium stationary state in open systems were quite different based on the two different approaches: the variational principle and MSBE. However, in order to understand new phenomena related to the non-equilibrium phase transition that can occur between the two limits, we need other frameworks that can be applied simultaneously to the non-equilibrium system and thermodynamic systems, where thermalization processes should be correctly included by non-Markov treatments [42]. In the next section, we show a unified theory for the e-h-p system based on the Keldysh-Green functions.

Non-equilibrium Green's function approach including both thermal polariton condensates and kinetic semiconductor lasers with Coulomb correlation effects
We have described the properties of the condensed phases in thermal equilibrium and the semiconductor lasers in non-equilibrium steady states in section 2, both of which are potentially obtainable under different conditions of the e-h-p system. However, the relationships between the two are ambiguous in the sense that we cannot discuss the change from the thermal equilibrium condensed phase into the non-equilibrium lasing state because the two approaches are not connected directly. Therefore, in this section, we present a unified theory based on a nonequilibrium Green's function approach, which is an extension of [34][35][36] for two-level systems. The many-body effects due to the Coulomb interactions are considered at the level of the HF approximation, which is sufficient for our discussion. The theory presented here can reproduce both the thermal equilibrium and the non-equilibrium results in all density regimes. We describe this section independently of section 2 as long as possible because several parameters change their physical meanings in the two different limits, and therefore, we should present definitions of parameters clearly. Hence, we expect that section 3 can be read independently of the information described in section 2.

Model
In our model, we assume a situation in which the e-h system is first excited and subsequent intraband relaxations are caused by carrier-carrier scatterings and interactions with phonons. These relaxation processes immediately redistribute the carriers in order to follow the Fermi distribution function. Then, these carriers interact with the cavity photons that are connected to free space. Hence, the e-h-p system is intrinsically an open quantum system. In such a system, it is common to describe the damping and pumping of a system (denoted by S) by its interaction with a reservoir (denoted by R) having a large number of degrees of freedom. Then, we focus on the evolutions of the system's variables by integrating out the reservoir's variables [45,46].
In this sense, the e-h-p system can be modeled as schematically shown in figure 7(a), where the reservoir for the photonic leakage is the free-space vacuum field and those for pumping are modeled by pumping baths. Here, the pumping baths are introduced somewhat artificially in order to consider the essential points regarding intraband relaxations without describing their specific mechanisms: electrons and holes are injected into (extracted from) the e-h system by their respective pumping baths in a time scale of 1/γ ( figure 7(b)). This modeling seems to be reasonable because a carrier relaxation from one state to the other can be regarded as a carrier extraction from the former and injection into the latter.
Using such an approach, the whole Hamiltonian is described asĤ =Ĥ S +Ĥ R +Ĥ SR , wherê are the system, reservoir and their interaction Hamiltonians, respectively. Here,ĉ c,k =ê k and c v,k (=ĥ † −k ) denote annihilation operators for electrons with wavenumber k in the conduction (c) and valence (v) bands, respectively, andâ q is an annihilation operator for photons with wavenumber q in the cavity. Similarly,b c,k andb v,k denote fermion annihilation operators of pumping baths andˆ p is a boson annihilation operator of free-space vacuum fields. α k and ζ k in equation (34) are the coupling constants between the system and each reservoir, satisfying with the following definitions of the density of states: Here, the dependence on the wavenumber is neglected in equation (35) for simplicity.Ĥ kin ,Ĥ e−e andĤ e−p in equation (32) are described aŝ respectively. ω c,k , ω v,k and ω ph,q denote the conduction-band, valence-band and photonic dispersion relationships, respectively. U q denotes the Coulomb interaction and is the same as the definition in equation (5).

Formulation
Based on the above Hamiltonians, we now present our formulation by assuming that a coherent photon field can be formed only for the q = 0 state and that the system is in a steady state described by With these assumptions, we derive a self-consistent closed set of coupled equations for the unknown variables of a 0 , p k , n c,k , n v,k and µ. Toward this end, it is convenient to perform a gauge transformation to eliminate the time dependence of â 0 (t) and ĉ † v,kĉ c,k (t) [35,36]. This is equivalent to a rotating frame and corresponds to simple energy shifts of (33), (36) and (37). Unless otherwise stated, we use this transformed Hamiltonian where â 0 = a 0 and ĉ † v,kĉ c,k = p k are achieved. In this situation, a kinetic equation of motion for the coherent photon field a 0 can be obtained by a straightforward Heisenberg-Langevin approach [36,45]: where equations (35) and (36) are used in the derivation. Equation (42) is a member of our self-consistent coupled equations. The other equations are systematically derived by the nonequilibrium Green's function approach, as described below.

Non-equilibrium Green's function approach.
Here, we first outline the derivation of the remaining set of coupled equations in a somewhat general manner. With a definition of a closedtime non-equilibrium Green's function, and its behavior, Dyson's equation can be written as where τ i is a time variable on the closed-time path (denoted by C) and T C is the time ordering operator on that path [47]. Here, G 0 C denotes the Green's function with no interactions and C is the self-energy. Then, it is straightforward to obtain the corresponding equation in a real-time formalism [47], where the matrix form of the Green's functions and the self-energy is Here, G 0 , G and are described by X for notational simplicity and R, A and K denote the retarded, advanced and Keldysh parts, respectively. Because only the time difference t 1 − t 2 determines the time dependence of X R/A/K (t 1 t 2 ; α 1 α 2 ; k) in a steady state, a Fourier transform on t 1 − t 2 can be performed and Dyson's equation can be rewritten as by the matrix form of equation (46). This equation suggests that G(ν; k) can be obtained once a specific form of the self-energy (ν; k) is determined at some approximation level. Then, p k and n α,k (α = c, v) can be described as by equation (47) and In this manner, equations (42) and (48) with equation (47) can form a set of the coupled equations. However, it should be noted that the coupled equations can be closed only when G < (tt; αβ; k) in equation (48) can be written solely using our unknown variables a 0 , p k , n c,k , n v,k and µ. Therefore, our remaining tasks are to insert a specific self-energy into equation (45) and to check whether the obtained equations can be closed. Figures 8(a)-(c) show self-energy diagrams when only the lowestorder effects caused by the interaction Hamiltonians are considered (equations (34), (38), and (39)). Here, C,el−ph , C,el−el and C,SR are self-energies due to the light-matter coupling (equation (39)), Coulomb interaction (equation (38)) and system-reservoir coupling (equation (34)), respectively. Note that the diagrams in figures 8(a) and (b) are equivalent to the mean-field and HF approximations for light-matter coupling and Coulomb interaction, respectively. In this case, the self-energy C can be written as C = C,el-ph + C,el-el + C,SR ,

Self-energies.
where equations (5) and (40) are considered: the self-energy arising from the Hartree diagram ( figure 8(b)) vanishes, whereas the self-energy proportional to a 0 remains when Wick's theorem 22 is applied. δ C in equations (50) and (51) is defined as for for for where C 1 (C 2 ) in equation (53) denotes the forward (backward) branch of the closed time path C. Here, we note that the effects of photon mass and the center-of-mass motions of the e-h pairs are neglected by using these lowest-order self-energies (equations (50) and (51)). These effects would be included when higher-order self-energies are systematically taken into account. However, such a treatment makes the formalism complicated and therefore will be presented elsewhere. Fortunately, at least in thermal equilibrium regimes [17-21, 39, 40], these effects are not so large when the temperature is far below the critical temperature because the center-of-mass motions for e-h pairs are frozen and such e-h pairs are coupled only with the cavity q = 0 mode due to the law of momentum conservation. B C in equation (52) is defined as Within a Born approximation, the reservoir's Green function B C (in equations (52) and (54)) is not affected by the system's Green functions and can be determined only from the reservoirs' Hamiltonian (equation (33)). In such a situation, the reservoirs can be considered to be in thermal equilibrium and the relationship can be used in the description of equation (54). Here, β = 1/k B T denotes the pumping bath temperature and µ B c (µ B v ) is the chemical potential of the pumping bath for electrons in the conduction (valence) bands. Then, a straightforward conversion of equations (49)-(52) into the real-time formalism [47] and the subsequent Fourier transform yields a final form of (ν; k): where with the Fermi distribution function (with gauge-transformed energies) Here, k and BGR α,k in equation (56) are defined as respectively. As described in section 2, k is equivalent to the composite order parameter commonly used in a Feshbach resonance of a uniform Fermi atom gas [39,40] and is also identical to the renormalized Rabi frequency in the MSBE [32]. The second term inh k corresponds to the Coulomb enhancement and BGR α,k is equal to the renormalization of the single-particle energies [32].

Self-consistent closed set of coupled equations.
We have now obtained the explicit form of the self-energy (equation (56)). In addition, the inverse matrix of the free Green's function can be easily given by Therefore, equations (47) and (48) with an inverse Fourier transform yield the following equations: where L k (ν) is defined as and the notation is transformed into an e-h picture [32] by and In addition, the following definitions are also introduced: Thus, a set of coupled equations (equations (42), (61)-(63)) are obtained, which number one less than the number of unknown variables (a 0 , p k , n c,k , n v,k and µ). This is because the origin of the phase of k and a 0 can be taken arbitrarily. Therefore, without loss of generality, we can set k to be real for one arbitrary wavenumber k 0 : Thus, we successfully obtained a self-consistent set of the coupled equations (equations (42), (61)-(63), and (70); denoted by set (A) hereafter). Furthermore, one can easily check that these equations are evidently closed. This set of coupled equations is one of our main results.
In particular, the closeness of the equations is nontrivial, as noted before, because there is a possibility that it is mixed with variables other than a 0 , p k , n c,k , n v,k and µ by the Coulomb interaction. However, fortunately, the self-energy in equation (56) does not include any Green's functions that cannot be written by a 0 , p k , n c,k , n v,k and µ within our approximation level. Therefore, the equations are closed. Finally, we comment that equations (42) and (61) can be combined into one equation by the composite order parameter k (equation (59)) with the definition Hence, k can be regarded as unknown variables instead of a 0 and p k . As a result, a set consisting of equations (62), (63), (70) and (71) (denoted by set (A) * ) is also closed with new unknown variables k , n c,k , n v,k and µ. In this case, a 0 and p k can be calculated from equations (42), (59) and (61) after the solutions of k , n c,k , n v,k and µ are found. U eff,κ q,k is found to play the role of an effective potential in the thermal (quasi-thermal) equilibrium case. Thus, we have derived a self-consistent closed set of coupled equations. However, the obtained equations are complicated for understanding the physical meaning. As a first step, we consider two typical situations in the next subsection, which reproduce thermal (quasi-thermal) equilibrium and non-equilibrium steady-state results.

Situations recovering thermal (quasi-thermal) equilibrium and non-equilibrium results
Here, the obtained coupled equations (set (A) or equivalently (A) * ) are examined in situations where thermal (quasi-thermal) equilibrium and non-equilibrium steady-state results would be reproduced. For this purpose, ω e,k = ω h,k and a charge neutrality µ B e = µ B h are assumed for simplicity and the following two situations are discussed: In both cases, the integral calculations of equations (61)-(63) and (71) are simplified because F B e (ν) and F B h (ν) can be treated as constants around ν =ω − eh,k ± E k , which is the position of two peaks in L k (ν) (see also figures 9(a) and (b)). 3.3.1. Thermal and quasi-thermal equilibrium results. First, the condition (I) min[2h E k ] µ B − µ + 2hγ + 2k B T is discussed. In this case, higher and lower plateaus of F B e (ν) and F B h (ν) lie at the two peak positions of L k (ν) (see also the inset of figure 9). Hence, F B e (ν) and F B h (ν) in equation (71) can be approximated by the values at ν =ω − eh,k ± E k : whereω − eh,k is replaced by −h −1 µ/2 because ω e,k = ω h,k in equation (68) and the last approximation is based on condition (I) and the actual shape of tanh(x). As a result, the following equation is obtained: is used in the derivation. Here, equation (75) is equivalent to the BCS gap equation for electrons and holes (equation (18) in section 2.1.3). In this sense, U eff,κ q,k and µ denote the effective potential between electrons and holes and the chemical potential of the e-h-p system, respectively. In the same manner, min[2h E k ] is equal to the minimum energy required for breaking an e-h pair. Now, we discuss the physical meaning of condition (I). First, it should be noted that µ = µ B is achieved when κ is set to zero because the imaginary part of equation (71) for k = k 0 should be zero [34][35][36]. The limit κ → 0 means that the e-h-p system is thermalized solely by the pumping baths. Hence, µ = µ B is a natural result because the e-h-p system should be in chemical equilibrium with the pumping baths. Therefore, the situation is undoubtedly a thermal equilibrium. In contrast, µ B = µ implies that the system is affected by the photon leakage (κ = 0). In this case, the system should have some aspects of non-equilibrium. However, if (I) is satisfied, the behavior of the e-h-p system is still described by the thermal-equilibrium BCS gap equation (equation (75)), the temperature of which is the same as that of the pumping baths but the chemical potential is not (µ = µ B ). Hence, this situation can still be classified as quasi-thermal equilibrium. In this context, (I) is interpreted as a condition for achieving quasithermal equilibrium. Then, the physical meaning of (I) also becomes clear: (I) is a condition for the e-h pairs not to be broken by the flux of particles due to the chemical non-equilibrium (= µ B − µ), thermalization-induced dephasing (= 2hγ ) and temperature effect (= 2k B T ).
Here, the significance of the particle flux was numerically pointed out by Szymanska et al [34,35] by regarding the e-h system as simple two-level systems. However, it is now clarified that the particle flux caused by µ B − µ has the effect of breaking e-h pairs by the explicit form of (I).

3.3.2.
Non-equilibrium steady-state results. Next, we discuss condition (II): µ B − µ 2h E k + 2hγ + 2k B T . In this case, the lower plateau of F B e (ν) and the upper plateau of F B h (ν) play a role at the two peak positions of ν =ω − eh,k ± E k (see also the inset of figure 9). Because the values of F B e (ν) and F B h (ν) do not change much between the two peaks, they can be approximated by the value at ν =ω − eh,k . Then, the integrals of equations (61)-(62) can be explicitly calculated by using equation (76) and As a result, the equations can be obtained by defining the population inversions of the system and pumping baths as N k = n e,k + n h,k − 1 and N 0 eh,k 2, respectively. After rewriting equation (78) as with the relationship − 1 2 the following equations are finally obtained: Here, in the derivation, we used the relationshipω + eh, (80). Equations (81) and (82) combined with equation (42) are identical to the MSBE described in section 2.2 (γ = γ ⊥ = 2γ ) under the steady-state assumptions of equation (41). Hence, we can see that the thermalization after the pumping is indeed described by γ (equation (82)), but it inevitably involves thermalization-induced dephasing (equation (81)) [34][35][36]. Thus, MSBE is successfully reproduced from our coupled equations (set (A)) under condition (II).
Here, condition (II) implies that the behavior of the e-h-p system is governed by (i) a large degree of non-equilibrium due to the photon leakage and (ii) a strong excitation, for the following two reasons: firstly, a large value of µ B − µ means that the system is significantly affected by the photon leakage because µ B − µ = 0 when the leakage is blocked (κ → 0), as described above. This suggests that a large degree of non-equilibrium is achieved. Secondly, it is easily confirmed that a complete population inversion of the pumping baths (N 0 k ∼ 1) is always achieved for wavenumber k that satisfies (II) because F B e (ω − eh,k ) ≈ −1 and F B h (ω − eh,k ) ≈ 1 (see also the inset of figure 9). This implies that the system is strongly excited and in a high-density regime. Therefore, it is quite natural that the MSBE is reproduced under condition (II) because the high-density regime achieved by the balance of the pumping and loss is the very situation assumed in the MSBE.
Here, we stress that condition (II) depends on the wavenumber k. Therefore, in a strict sense, the MSBE is not correct for all k. However are assumed to be the values estimated from the Fermi distribution function, MSBE is obtained for all k. In other words, MSBE is often used for all k by assuming that N 0 k is the value of the Fermi distribution function. In contrast, this assumption can be eliminated when the presented equations (set (A) or equivalently (A) * ) are directly solved.
Thus, the non-equilibrium steady-state result, i.e. MSBE, can be reproduced under condition (II). However, the condition recovering the MSBE is not limited to the presented one, e.g. a high-temperature limit (k B T h E k +hγ ) can also reproduce MSBE formally. Although this condition does not contain direct information about the non-equilibrium features, the MSBE is naturally recovered because the Markov approximation is justified [36,42] (see also section 2.2.1).

Numerical results for the BEC-BCS-laser crossover
Thus, we have described situations for recovering the thermal (quasi-thermal) equilibrium and non-equilibrium results. However, the presented conditions, (I) and (II), cannot be discussed without finding solutions for the coupled equations in section 3.2.3 because the conditions depend on the solved values of k and µ. Moreover, the relationships between several condensed phases (e.g. e-h BCS, exciton BEC and photonic BEC) and lasing remain ambiguous. Therefore, we show typical numerical results.
In our calculations, the cavity resonance is set to 20 meV greater than the band-gap energy (hω ph,0 = E g + 20 meV) and the k dependence of k is eliminated by using a contact potential of U q = U . In addition,hω e,k =hω h,k =h 2 k 2 /2m + E g /2 are assumed with a charge neutrality µ B e = µ B h and T = 0 K for simplicity. In this context, the discussion presented here is qualitative even though the parameters are determined to be realistic as long as possible 10  in the range of µ B <hω ph,0 (µ B − E g < 20 meV). This effect can be explained by the e-h pair breaking due to the thermalization-induced dephasing, as described in section 3.3.1, becausē hγ increases to as much ash . As a result, in figure 10(c), min[2h E k ] with large γ (green points) becomes comparable to and less than E Gap (green dashed line) for µ B <hω ph,0 , which is in contrast to the results with small γ (red data in figure 10(c)). Hence, the e-h-p system for µ B <hω ph,0 is described by thermal equilibrium theory when γ is small.
On the other hand, when µ B >hω ph,0 (µ B -E g > 20 meV), both the results (red and green points) are different from the thermal equilibrium results in figure 10(a). At the same time, µ is fixed at around E g + 20 meV =hω ph,0 in figure 10(b) even though µ B is increased (i.e. µ = µ B ). These two results suggest that the e-h-p system is out of thermal equilibrium with the pumping bath. However, the e-h-p system is still in quasi-thermal equilibrium because min[2h E k ] E Gap is achieved, as shown in figure 10(c) (see also section 3.3.1). Here, the quasi-thermal equilibrium is more easily achieved for larger γ because min[2h E k ] becomes large when γ is increased. This is in contrast to the case for µ B <hω ph,0 , where the system can be viewed as being in thermal equilibrium for small γ . In the present case (µ B >hω ph,0 ), it is important that µ B − µ takes a large value (figure 10(b)), because this means that the photon leakage has a considerable effect, unlike what happens in the case of µ B <hω ph,0 . Hence, the e-h-p system cannot be in thermal equilibrium when the thermalization speed is not sufficient, as mentioned in section 1. The thermalization speed of the e-h-p system is determined by that of the e-h system (= 1/γ ) and the speed of the light-matter coupling (= 1/g). As a result, the system approaches quasi-thermal equilibrium when γ is increased because the thermalization time is shortened by large γ .
Next, the effect of κ is discussed by focusing on the red and blue results in figures 10(a)-(c). Here,h is largely decreased for µ B >hω ph,0 but not for µ B <hω ph,0 when κ is increased from 1.0 eV (red) to 100 eV (blue) in figure 10(a). This result is natural because the effect of photon leakage becomes more important for µ B >hω ph,0 . Consequently, min[2h E k ] becomes smaller than E Gap at µ B − E g ∼ 70 meV and goes below E MSBE at µ B − E g ∼ 95 meV for large κ (blue data in figure 10(c)). This means that there are several k states that do not follow the gap equation for µ B − E g 70 meV. Then, the number of states that can be described by the MSBE begins to increase after the transition regime (70 meV µ B − E g 95 meV). As a result, nonequilibrium features become dominant for µ B − E g 95 meV, the signature of which appears as a kink in the plot ofh and min[2h E k ].
In order to discuss the connection between thermal (quasi-thermal) equilibrium and non-equilibrium features in more detail, the electron distribution function n e,k (= n h,k ) and n e,k once approaches the Fermi distribution (figures 11(A)-(C)) and becomes broad (figures 11(C)-(E)) as the excitation is increased. Simultaneously, p k shows a peak around the Fermi momentum (figures 11(A)-(C)) and then spreads out to form a plateau p k ∼ 0.5 (figures 11(C)-(F)). In this situation, it is also found that the photon fraction of the condensed phase increases sharply (inset of figure 11). These behaviors of n e,k , p k and the photonic fraction are consistent with the results shown in section 2 and previous works [19][20][21]. This is quite natural because the results shown in figures 11(A)-(E) are obviously in thermal (quasi-thermal) equilibrium (min[2h E k ] > E Gap in figure 10(c)). Hence, the relevant condensed phases can be easily identified as (A) exciton BEC, (B) e-h BCS, (C) e-h polariton BCS, and (E) photonic BEC.
In the case of large κ, major differences can be observed for µ B >hω ph,0 : although n e,k and p k in figure 11(F) are still close to those in figure 11(C), they gradually go out of thermal (quasi-thermal) equilibrium (figure 11(G)), and finally, characteristic dips appear in n e,k and p k ( figure 11(H)). Simultaneously, the MSBE-dominated (red shaded) area is broadened from the cavity resonance. Hence, the dip in n e,k is equivalent to kinetic hole burning [29,30], which means that a lasing operation is obviously achieved. Thus, crossovers between the exciton BEC, e-h BCS and lasing operation can be found in the calculation. However, we note that the obtained lasing is still different from the conventional one because the gain is not caused by the e-h plasma but rather by a Coulomb-correlated condensed phase described by the BCS gap equation. This effect is reflected in the shape of | p k | where the distribution is widely spread in k space ( figure 11(H)) due to the Coulomb correlations between the various k states. In this sense, BCS-coupled lasing is achieved in figure 11(H), the features of which cannot be treated by the simple non-interacting models [34][35][36]. One can notice that the behaviors of n e,k and | p k | are close to those in panel (E) of figure 6, which is already influenced by BCS-type e-h pairing, as described before. The other aspects, e.g. spectral behaviors, are beyond the scope of this paper, but we stress that this somewhat new mechanism can be justified because our approach can be applicable for both thermal equilibrium and non-equilibrium steady states with Coulomb correlations.

Conclusions
We have discussed various phenomena that can be achieved in the e-h-p system. Thermalequilibrium phases are first discussed by the variational approach at arbitrary density, and the crossovers between the exciton BEC, e-h BCS, exciton-polariton BEC and photonic BEC phases are presented. In particular, a large photonic fraction and photon-mediated interactions between electrons and holes are obtained in the photonic BEC phase. Then, non-equilibrium lasing operations are studied by the MSBE approach in steady states, where the lasing is evidenced by the kinetic hole burning. In such a situation, the e-h-p system is far from thermal equilibrium and its behavior is governed by the balance of the gain and the loss. However, as noted above, the situations covered by the variational approach and the MSBE are not connected with each other and their intermediate regimes cannot be discussed. Therefore, we finally presented a unified framework based on the non-equilibrium Green's function approach. Self-consistent coupled equations have been derived with clear explanations of their physical meanings. In particular, we have shown that the loss of the cavity photons leads to a difference between the chemical potentials of the e-h-p system and the pumping baths, leading to e-h pair breaking. The e-h-p system then moves into the non-equilibrium regime, and a lasing operation is achieved when the difference between the chemical potentials becomes sufficiently large. As a result, BEC-BCS-laser crossovers can be obtained. The lasing obtained as a result of the crossovers, however, is not identical to the conventional one in the sense that the lasing is achieved not by the e-h plasma phase but by a Coulomb-correlated condensed phase where the Cooper(-like) pairing plays an important role. Therefore, we referred to this type of lasing as the BCS-coupled lasing in order to distinguish it from the normal photon lasing. However, stabilities and spectral features, for example, are still unclear although these behaviors are important for discussing the characteristics of lasing. In this respect, further studies are needed for a full understanding of this phenomenon. An inclusion of the photon mass and the e-h center-of-mass fluctuation [18,39,40] is also critical when discussing the connections to the BEC phases, in particular in high-temperature regimes around the critical temperature. It is also important to use the real Coulomb potential when the theory is quantitatively compared with relevant experiments although the contact potential is employed to discuss the qualitative behaviors in this paper.