Pump-power-driven mode switching in a microcavity device and its relation to Bose-Einstein condensation

We investigate the switching of the coherent emission mode of a bimodal microcavity device, occurring when the pump power is varied. We compare experimental data to theoretical results and identify the underlying mechanism to be based on the competition between the effective gain on the one hand and the intermode kinetics on the other. When the pumping is ramped up, above a threshold the mode with the largest effective gain starts to emit coherent light, corresponding to lasing. In contrast, in the limit of strong pumping it is the intermode kinetics that determines which mode acquires a large occupation and shows coherent emission. We point out that this latter mechanism is akin to the equilibrium Bose-Einstein condensation of massive bosons. Thus, the mode switching in our microcavity device can be viewed as a minimal instance of Bose-Einstein condensation of photons. We, moreover, show that the switching from one cavity mode to the other occurs always via an intermediate phase where both modes are emitting coherent light and that it is associated with both superthermal intensity fluctuations and strong anticorrelations between both modes.

We investigate the switching of the coherent emission mode of a bimodal microcavity device, occurring when the pump power is varied. We compare experimental data to theoretical results and identify the underlying mechanism to be based on the competition between the effective gain on the one hand and the intermode kinetics on the other. When the pumping is ramped up, above a threshold the mode with the largest effective gain starts to emit coherent light, corresponding to lasing. In contrast, in the limit of strong pumping it is the intermode kinetics that determines which mode acquires a large occupation and shows coherent emission. We point out that this latter mechanism is akin to the equilibrium Bose-Einstein condensation of massive bosons. Thus, the mode switching in our microcavity device can be viewed as a minimal instance of Bose-Einstein condensation of photons. We, moreover, show that the switching from one cavity mode to the other occurs always via an intermediate phase where both modes are emitting coherent light and that it is associated with both superthermal intensity fluctuations and strong anticorrelations between both modes.
In this article, we study the switching of the coherent emission mode in bimodal micropillar lasers occurring when the pump power is ramped up. By comparing experimental data to theoretical results based on a phenomenological model, we identify the basic mechanism underlying the mode switching to be the competition between effective gain on the one hand and the intermode kinetics on the other. Namely, the mechanism that selects which of the modes shows coherent emission is found to be fundamentally different for weak pumping (just above the threshold) and in the limit of strong pumping. For weak pumping, the selected mode (i.e. the * ham.leymann@gmail.com † dv@pks.mpg.de ‡ eckardt@pks.mpg.de mode selected for coherent emission) is characterized by the largest effective gain and coherent emission corresponds to lasing. In contrast, for strong pumping, the selected mode depends neither on the coupling to the gain medium nor on the loss rates of both modes. Instead it is determined completely by the intermode kinetics, i.e. by processes that transfer photons from one mode to the other. We show that this mechanism is formally identical to the one leading to Bose-Einstein condensation of an ideal gas of massive bosons (i.e. with a conserved particle number) in contact with a thermal bath. Therefore, the mode switching in our system can be viewed as a minimal instance of Bose-Einstein condensation of photons.
The question, whether a system of photons (or bosonic quasiparticles) with non-conserved particle number can undergo Bose-Einstein condensation in a similar way as a thermal gas of massive bosons has raised considerable interest in the last decade. Here the problem to be overcome is to achieve a quasi-equilibrium situation, where a single mode acquires a macroscopic occupation via a thermalizing kinetics, despite the non-equilibrium nature of the system resulting from particle loss to be balanced by pumping. Beautiful experiments, showing that such a situation can indeed be achieved, have been conducted in systems of exciton-polaritons [23][24][25][26][27][28][29] magnons [30][31][32][33], and photons in dye-filled cavities [34]. While the microcavity device investigated in this article is simpler than these systems, in the sense that it is described in terms of two relevant modes only, it captures one of the most important aspects of photon condensation in a minimal fashion, namely that the condensate mode is selected not l l l l l l l l l l l l l l l l l Illustration of the relevant processes in our bimodal system. The intermode kinetics between the modes h and l is determined by the rates Ri→j. The modes lose their photons with rates i. This loss is compensated by new excitation from the quantum dots (QDs) with gain rates gi. The quantum dots are pumped with rate P and lose their excitations spontaneously with rate τ −1 . by pumping but rather by the kinetics of the photons.
The starting point of our theoretical decription of a bimodal micropillar is a generic phenomenological master equation describing the relevant processes of the system [10,35]: the coupling between the pumped medium and the cavity modes, loss, and the intermode kinetics ( Fig. 1). Such birth-death models have also been used to study (analogs of Bose-Einstein condensation in) population dynamics [36,37], transport [38,39], and networks [40] as well as quantum gases of massive bosons [41,42]. This approach, which is different from the microscopic modeling pursued in former studies of bimodal (micro) lasers [8,10,19,[43][44][45][46][47] and interacting excitonpolariton systems, e.g., Refs. [48][49][50], provides excellent agreement with the experimental data (Fig. 2). In order to treat these equations analytically, we work out an approximation scheme that combines the theory of Bose selection, which was recently developed to describe (nonequilibrium) Bose condensation of ideal gases of massive bosons with conserved particle number [41], with particle loss and the coupling to a pumped reservoir (gain medium). We justify this approximation by exact numerical simulations. Somewhat counter-intuitively, our theory shows that it is the limit of strong pumping, where the selected mode is determined by the intermode kinetics, which is described by terms that are formally identical to those appearing in the description of massive bosons and their equilibrium condensation. We also find that the mode switching occurs via an intermediate phase where both modes are emitting coherently (see phase diagram in Fig. 3 below).
We investigate also the statistical properties of our system. In bimodal lasers, where one mode dominates the emission for all pump rates, the non-lasing mode exhibits usually super-thermal intensity fluctuations [7,10] and the emission of both modes is strongly anticorrelated [8,46,[51][52][53]. In the situation where mode switch- Measured and calculated microcavity input-output characteristics (high-effective-gain mode in red, low-effectivegain mode in gray). The injection current I is proportional to the pump rate P . Here, circles (connected to guide the eye), solid curves, and dashed curves refer to experimentally measured data, exact Monte-Carlo results, and asymptotic meanfield theory, respectively. Panel (a) depicts the mean occupations of the modes and the number of excited carriers (light green) in arbitrary units. The colored bars at the bottom of this panel mark the phases determined by the asymptotic theory (see Sec. IV C). Panel (b) shows the autocorrelation g (2) ii for both modes in comparison with the value 2 expected for thermal emission (dotted black line). Panel (c) depicts the crosscorrelation g (2) hl in comparison to the value 1 for correlated emission (dotted black line). For the theoretical curves, we use s = 0, G = 0.77 as well as (in units of the spontaneous loss rate τ −1 ) g h = 1.6 · 10 −3 , g l = 2.1 · 10 −3 , h = 2.2 · 10 −2 , l = 3.8 · 10 −2 , R l→h = 1.7 · 10 −4 , A h→l = 8.5 · 10 −6 (see App. A).
ing occurs, we find that the super-thermal intensity fluctuations of the non-selected mode and strong anticorrelations occur whenever a mode starts or ceases to be selected. We show that these experimentally observed statistical properties can be described theoretically by an effective reduction of the spontaneous inter-mode transitions caused by mode interactions.
This paper is organized as follows. In Sec. II the experimental setup is presented. The theoretical description in terms of a master equation is introduced in Sec. III. An analytical theory of the mode switching and its relation to Bose-Einstein condensation is then worked out in Sec. IV. Finally, in Sec. V, we investigate the statistical properties of the system, before coming to the conclusions in Sec. VI.

II. EXPERIMENT
Electrically-pumped quantum-dot micropillars are fabricated by etching of a planar AlAs/GaAs distributed Bragg reflector λ-cavity in which a single active layer of self-assembled In 0.3 Ga 0.7 As quantum dots is embedded centrally [54]. A detailed description can be found in Ref. [55]. The micropillar used in this study has a diameter of 3.0 µm. Due to the strong confinement of light, the micropillars exhibit a spectrum of discrete modes. The fundamental modes HE 1,1 are composed of two orthogonally linearly polarized components, which are ideally energetically degenerate. In reality, however, asymmetries in the manufacturing process, which results in slightly elliptical structures, lift the energetical degeneracy of the fundamental modes [1,56]. The resulting mode splitting of the micropillar used in this study is (42 ± 2) µeV . Besides a finite mode splitting the two fundamental modes also exhibit slightly different quality factors, 14000±1500 and 12500±1500, respectively. The different spectral and local overlaps of the modes with the gain medium and the modes polarization alters their coupling to the quantum dot emitters. As a consequence, the former mode (mode h) is characterized by a higher effective gain [i.e. gainloss ratio, see Eq. (4) below] than the latter one, which has a lower effective gain (mode l).
A high-resolution (25 µeV ) micro-electroluminescence setup is used to characterize the micropillars spectrally at cryogenic temperatures of 15 K. A linear polarizer as well as a λ/4-plate in front of the monochromator enables polarization-resolved spectroscopy. For statistical analysis via the autocorrelation function g (2) ii (τ = 0) with zero delay time of the emission a fiber-coupled Hanbury Brown and Twiss setup is used. The temporal resolution of the Hanbury Brown and Twiss configuration based on fast Si avalanche photon diodes is τ res = 40 ps. For measuring the equal-time crosscorrelation function g (2) hl (τ = 0) between the orthogonally polarized micropillar modes the emission is selected by a polarization maintaining 50/50 beamsplitter and the split beam is directed to a second identical spectrometer -with the polarizer in front of the second monochromator oriented orthogonally to the first. These equal-time correlations are defined by whereb i is the bosonic annihilator operator of mode i. The emission characteristics of the CW-pumped micropillar are shown in Fig. 2. In (a) the output characteristics of modes h (red) and l (gray) are shown as a function of the injection current I (quantifying the pumping strength). Up to about 13 µA (phase A) both modes are below threshold and show a small increase in output emission only. Between 13 µA and 36 µA (phase B) the mode h is selected and its emission increases strongly while that of mode l remains small. Between 36 µA and 50 µA (phase C) both modes are selected, but while the emission intensity of mode l increases that of mode h decreases. Beyond 50 µA (phase D) only mode l is selected. In Fig. 2 (b) and (c) the zero-delay autocorrelation g (2) ii and crosscorrelation g (2) hl functions are plotted. Note that, due to finite temporal resolution of the avalanche photodiode detectors and the low coherence time of the modes, we could not properly resolve g (2) ii experimentally for values above 1 in phase A and B [57].

III. MASTER EQUATION
Our starting point for the theoretical description of the system is a phenomenological master equation for the probabilities ρ n N to find the system in a state with N excited emitters and photon numbers n = (n h , n l ) in the high-and the low-effective-gain mode. It takes the form and shall be solved for the steady-state state obeying The first term on the right hand side of the master equation (2) describes how photons leave and enter the cavity modes via loss and coupling to the emitters. It is given by where P and τ −1 denote the pump and the loss rate of the emitters, respectively, g i quantifies the gain of cavity mode i from the emitters, and i is the loss rate of cavity mode i. An additional or removed photon in mode i is denoted by ±e i , [i.e., e h = (1, 0), e l = (0, 1)]. The modes h and l are defined by the higher and lower effective gain g i /l i , respectively, so that we obtain an effective-gain ratio The terms contained in C laser (ρ) are sufficient for a theoretical description of single-mode lasing in mode h [35]. The second term of the master equation (2) captures the intermode kinetics and reads It is characterized by the rates R i→j for a transition from mode i to mode j. The rate asymmetry of the direct intermode transitions is generally nonzero. The origin of this rate asymmetry was attributed to stimulated scattering due to carrier population oscillations, e.g., in coupled photonic crystal nanolasers [18]. Furthermore asymmetric backscattering of electromagnetic waves was observed in optical microcavities [58][59][60]. In our system the rate asymmetry A h→l is positive, so that the low-effective-gain mode l is favored by the intermode kinetics. The parameter s in Eq. (5) quantifies the ratio between spontaneous and induced intermode transitions. Its natural value is s = 1. However, a reduction of s to values s < 1 is a simple way to effectively capture inter-mode interactions that lead to a relative enhancement of transitions into strongly occupied modes. We will show below that, while s has (practically) no impact on the phase transition and the mean occupation(s) of the selected mode(s), it does affect the occupation of the non-selected mode. Only for s < 1, the master equation can describe the experimentally observed super-thermal fluctuations g (2) il > 2 of the nonselected mode. The numerical data shown in Figs. 2, 4, and 5 are obtained for s = 0.
The form of the term C kin (ρ) capturing the intermode kinetics is identical to that of the master equation for an ideal gas of massive bosons in contact with an environment [41]. If such a system of massive bosons is coupled to a thermal environment characterized by the temperature T , the intermode rates obey R j→i /R i→j = exp[−∆ ij /(k B T )] with Boltzmann constant k B and energy splitting ∆ ij = E i − E j between modes (singleparticle states) with energy E i . In the quantum degenerate regime of low temperature or high boson density, the system will form a Bose-Einstein condensate in the single-particle state of lowest energy. When increasing the total number of bosons N B in this regime, the occupation of an excited mode i approaches the finite value n i = (e Ei/k B T − 1) −1 , while the ground-state occupation increases linearly with N B . Even for a finite number of discrete energy levels i, in the limit N B → ∞ this behavior clearly describes Bose-Einstein condensation as the macroscopic occupation of one single-particle state (see e.g. Ref. [42]).
The most intriguing result of this paper, shown below, is that the behavior of the bimodal system in the limit of strong pumping strength P closely resembles that of a Bose-Einstein condensed gas of massive bosons in equilibrium. Remarkably here the selected mode does not depend on the effective gain, but is determined exclusively by the intermode transitions R i→j . This countner intuitive result is related to the fact that the intermode kinetics scales quadratically, but gain and loss only linearly with the mode occupations.
It is instructive to define the dimensionless parameter ε, which can be interpreted as the ratio ε = ∆ eff hl /(k B T eff ) of an effective energy splitting ∆ eff hl between both modes and an effective temperature T eff . In the limit of strong pumping, the intermode kinetics makes the photons condense into the mode corresponding to the lower effective energy. That is for ε > 0 or A h→l > 0 (ε < 0 or A h→l < 0) a Bose-Einstein condensate of photons is formed in mode l (mode h). In contrast, it is always the mode h, characterized by the higher effective gain, that starts lasing when the pump power P is ramped up. Thus, a rate asymmetry A h→l > 0 implies a switching from lasing in mode h to condensation in mode l, when the pump power is ramped up. In the following section, we derive an analytical theory, which describes this effect and establishes the analogy to equilibrium Bose-Einstein condensation for strong pumping.

A. Mean-field approximation
In order to obtain a closed set of kinetic equations for the mean mode occupations n i = n,N ρ n N n i and the mean number of excited emitters N = n,N ρ n N N , we perform the mean-field approximation n inj ≈ n i n j and Nn i ≈ N n i .
This approximation, which ignores non-trivial twoparticle correlations, is later justified by comparing it to exact solutions of the full master equation (2) obtained from Monte-Carlo simulations (see Fig. 2). Employing it, we derive kinetic equations of motion for the mean occupations:

B. Asymptotic theory
In a next step, for the sake of finding an analytical expression for the mean occupation(s) of the selected mode(s), in Eq. (10) we neglect spontaneous processes relative to corresponding stimulated ones, ( n i + a) n i with a = 1, s. This approximation is valid asymptotically in the limit of large occupations of the emitters and the selected mode(s). Note, for high-β cavities, where almost the entire spontaneous emission goes into the cavity modes, this assumption is plagued by the strong presence of spontaneous emission at the first threshold. Still, this assumption is valid in the coherent emission regime and, in particular, when describing transitions, where the selected modes change. Using the asymptotic approximation above, the stationary solution of Eq. (10) obeys This equation is solved by mean occupations that can be divided into two classes. For the non-selected modes i / ∈ S, one finds the trivial solution n i = 0, whereas the occupations of the selected modes i ∈ S obey a linear set of equations: The occupation of the non-selected states, which vanishes in leading order [Eq. (12)], can be computed in the next order of our approximation. For this purpose, we take into account those terms, which are linear in the number of excited emitters N or the occupations n i of the selected modes. We find: The dependence of the number of excited emitters on the pumping can be obtained by analogue reasoning, Following the strategy recently employed for massive bosons [41], the set of selected modes S can now be determined by the physical requirement that all modes (both selected and non-selected modes) must have positive occupations, For the non-selected modes this implies that the denominator of Eq. (14) has to be positive. Thus, the set of selected modes and their occupations are independent of the parameter s since s occurs neither in Eq. (13) nor in the denominator of Eq. (14).

C. Phase diagram of the bimodal microcavity system
Based on the asymptotic theory, we can now compute the non-equilibrium phase diagram. The concept of selected modes, which appeared naturally in the asymptotic theory, clearly separates the parameter space into Plotted versus the injection current I that is proportional to pump rate P and the ratio of the effective gain rates G = g l h /g h l with fixed g h l / h . Except for g l all parameters are the same as in Fig. 2. The dashed line indicates the value of g l used in Fig. 2. While for G < 1 all phases are present when increasing the pump rate P , for G > 1 the mode l, that is favored by the direct mode coupling becomes actually the high-effective-gain mode so that no mode switching occurs. different phases, where no, one, or both modes are selected. A transition, where one of the non-selected modes becomes selected and starts emitting coherent light, is indicated by the divergence of its occupation described by Eq. (14). In turn, a selected mode ceases to be selected when its occupation, obtained by solving Eq. (13), drops to zero. In this section we will use this argument to compute the phase boundaries analytically. The resulting phase diagram is depicted in Fig. 3.
In phase A, for small pumping power P , neither mode is selected, S = {}. According to Eq. (14) the mode occupations read and the number of excited emitters increases linearly with the pump rate, When the pumping is increased and reaches the occupation of the high-effective-gain mode diverges indicating the transition to a regime where this mode is selected and starts emitting coherently (see Eq. (14)).
Since before the transition no mode is selected yet, this asymptotic estimate for the critical pumping strength cannot be expected to mark precisely the threshold in the high-β limit. Note, however, that the estimates for further thresholds at P BC and P CD are not plagued by this problem, since they occur in the regime where at least one mode is selected already.
In phase B the high-effective-gain mode is selected, S = {h}, and the number of excited emitters is clamped at the threshold value [see Eq. (13) and Fig. 2 (a)] The excitation provided by increasing the pumping is directly transferred to the selected mode [61]. Consequently, the occupation of the selected high-effective-gain mode depends linearly on the pump rate [use Eq. (15)] The occupation in the non-selected low-effective-gain mode is given by [see Eq. (14)] In a case where the mode-coupling rates would favor the high-effective-gain mode (A h→l < 0), Eq. (22) would be valid for all pumping powers P > P AB . However, in our case, where the mode-coupling rates favor the loweffective-gain mode (A h→l > 0), increasing the pump rate (and with it also n h B ), eventually leads to the divergence of the right-hand side of Eq. (22). This occurs at the pump rate and indicates the transition to phase C. In phase C both modes are selected, S = {h, l}. The number of excited emitters increases again linearly with the pump rate P , The occupations of the high-and low-effective-gain mode de-and increase linearly with N C , respectively, When the number of the emitters N C reaches l /g l , the occupation n h C becomes zero, indicating the transition to the phase where this mode is no longer selected. The threshold pump rate P CD can be obtained from Eq. (24) analogously to the expression for P BC and reads Thus the extent of phase C is determined by the inverse effective gain ratio G [Eq. (4)].
In phase D, only mode l is selected, S = {l}. The number of excited emitters remains at the higher threshold value and the occupation of the selected mode l reads The occupation of the non-selected mode h is given by The crucial difference to Eq. (22) is that for A h→l > 0 an increase of n l D cannot produce a further root in the denominator of Eq. (30). Thus no further transition will occur [unless other parameters change as well when the pump power is ramped up]. Figure 3 shows the phase diagram resulting from the asymptotic theory with respect to the effective-gain ratio G (varied by varying the gain rate g l ) and the pumping strength P (proportional to the injection current I). While the precise shape of the phase boundaries depends on the parameters (and which of them are varied in order to modify the effective-gain ratio G), the topology of the phase diagram is generic. For G < 1 (and A h→l > 0), the system always undergoes a sequence of three transitions between the phases A, B, C, D when the pump rate is increased. For G > 1 (and A h→l > 0), where the mode labeled l actually becomes the high effective gain mode, only a single transition from phase A to phase D occurs. In summary, when the pump power P is ramped up, the system starts lasing in the mode characterized by the higher effective gain, whereas in the limit of strong pumping the selected mode is the one favored by the intermode kinetics. Moreover, the switching from selection of mode h to selection of mode l has to occur via an intermediate phase, where both modes are selected (unless the system is fine-tuned to G = 1).
The data plotted in Fig. 2 corresponds to a cut through this phase diagram following the dashed horizontal line in Fig. 3. The different phases obtained from the asymptotic theory are indicated by the colors at the bottom of panel (a). In Fig. 2(a), we can clearly see that the mean occupations obtained from the asymptotic theory (dashed lines) nicely reproduce the exact solution (solid line) of the master equation (2), which was obtained by Monte-Carlo simulations (for a detailed description of the method see Ref. [42]). This agreement justifies both the mean-field approximation and the asymptotic theory. More importantly, the theoretical curves also describe the experimental data (circles in Fig. 2) very well. In appendix A we explain how the parameters of our model are determined.

D. Relation to Bose-Einstein condensation
We have seen that in the limit of strong pumping strength P the selected mode is determined exclusively by the intermode kinetics, which is described by the rates R i→j . It is remarkable that in this limit neither the loss rates of the modes i nor their coupling g i to the emitters influences the selection of the mode. This counterintuitive result is related to the fact that gain and loss scale only linearly with the mode occupations, whereas the rates for the intermode kinetics have a quadratic dependence on the mode occupations. It implies that the mechanism leading to a macroscopic (or large) occupation of one of the modes is the same as the one that leads to the Bose-Einstein condensation of massive bosons in contact with a thermal environment. This is based on Bose-enhanced inter-mode kinetics (scattering) described, e.g., in Ref. [41] on the basis of a rate equation comprising the same terms as C kin [Eq. (5)]. Even though the system consists of two levels only, the notion of Bose condensation becomes sharp in the limit P → ∞, where the occupation of mode l, given by Eq. (29), approaches infinity, while that of mode h, given by Eq. (30), remains below a finite value. Note that also the number of excited emitters, given by Eq. (28), remains at a finite value in this limit. This has the important consequence that the occupation of the noncondensed mode h is determined completely by the intermode kinetics described by the rates R i→j . In the limit P → ∞ it approaches where ε is the parameter defined in Eq. (8). Thus, for s = 1 the occupation of mode l precisely corresponds to that of an excited state of energy εk B T in a Bose condensed system of massive bosons at temperature T . The fact that the "excited-state" occupation n h (the depletion) approaches a constant value for strong pumping, so that increasing the number of photons (bosons) in the system will only increase the condensate occupation, is another clear analogy to equilibrium Bose condensation. Irrespective of the value of s, for strong pumping the selection of the coherent emission mode is a result of the intermode kinetics.

V. PHOTON STATISTICS
In microlasers, where almost the entire spontaneous emission goes into a single mode, no sharp intensity jump is visible at the lasing threshold [5,35]. Instead the transition to coherent emission is indicated in the autocorrelation [62] of the emitted photons [57,[63][64][65][66][67][68]. To confirm the coherence properties of the emitted photons, we examine the photon statistics. The equal-time photon correlation functions (1), which we can write like measure the occupation number fluctuations of each mode and the crosscorrelation between the modes, respectively. In Fig. 2 (b) and (c) experimental and theoretical results for g (2) ii and g (2) hl are depicted. The theoretical values for g (2) ij are determined by a numerically exact Monte-Carlo simulation of the full master equation (2). As discussed above a change from g (2) hh = 2 to g (2) hh = 1 indicates the first threshold where the coherent emission in mode h sets in, as can be seen clearly in Fig. 2. For even stronger pumping, when entering or leaving phase C, in which both modes are selected, we observe pronounced anticorrelations between both modes (g hl < 1) as well as superthermal intensity flucutations (g ii > 2) of that mode i that changes its state from non-selected to selected or vice versa. Note that for our system parameters, phase C appears in a narrow interval of pump powers only so that its properties are overshadowed by those of the transitions BC and CD. As a result, the two minima of the crosscorrelations occurring at the transitions have merged to a single one.
In order to reproduce the measured superthermal intensity fluctuations at the transitions BC and CD theoretically, we have to choose s < 1, corresponding to the presence of intermode interactions. The best results are obtained for the value s = 0, which we used also in the simulations (the role of s will be discussed in more detail below and in Appendix C).
We will now investigate the signatures of the mode switching in the reduced two-mode photon distribution ρ n h ,n l = N ρ n h ,n l N , which gives the probability to find the system in a state with n h and n l photons in mode h and mode l, respectively. We compute this quantity by solving either the full master equation (phases A and B) or the reduced master equation for ρ n h ,n l (phases C and D, see Appendix B for details). Results for four different pump powers P , corresponding to the four phases A to D, are depicted in Fig. 4. The corresponding single-mode distributions ρ n h = n l ρ n h ,n l and ρ n l are shown as well.
In the non-selected phase A the distribution possesses a single maximum at n = (n h , n l ) = (0, 0). The selection of mode h (phases B and C) is associated with a maximum of the distribution at n = (n max h , 0) lying on the vertical axis, whereas the selection of mode l (phases C and D) with a maximum at n = (0, n max l ) lying on the horizontal axis. Thus, in phase C where both modes are selected, the distribution possesses two local maxima that are separated by a saddle point [69]. The emergence of the second maximum when entering phase C, is visible also in the occupation distribution of the mode that starts emitting coherently, as can be seen in Fig. 5 showing ρ n l directly before and after the transition from B to C. The build-up and later the presence of the sec- . These pump rates are also indicated by arrows in Fig. 2(b). The labels n h and n l mark the position of the expectation value for each mode. For a better comparison the single-mode statistics are depicted up to 0.025% in solid lines, while the remaining part is dashed. Note that plotted range of photon numbers increases from left to right. ond maximum in ρ n l is accompanied by the strong (superthermal) number fluctuations in this mode.
The impact of the effective parameter s is illustrated in Fig. 6, where the two-mode distribution is shown for s = 1 in phases B and C for the same parameters used in the corresponding panels in Fig. 4. For s = 1 states with zero occupation in one of the modes (situated along the axes of the plot) are much less attractive than for the value s = 0 considered before. As a striking consequence, the selection of both modes in phase C is not associated with two maxima in the distribution anymore, but rather with a single central maximum. The transition from phase B to phase C now corresponds to the shifting of the single maximum away from the horizontal axis. As a result, for s = 1 the system does not show the experimentally observed superthermal photon number fluctuations, as we discuss in more detail in Appendix C. It is an interesting observation that neglecting or taking into account the spontaneous emission between the modes has such a strong impact on the statistical proper- ties of the modeled bimodal system. Note, however, that despite this strong impact on the statistics, the mean occupations of the modes and the critical parameters for the phase transition do not show a strong dependence on s. This is also the result predicted by the asymptotic theory presented in the previous section.

VI. CONCLUSION
In this paper we investigated the pump-power driven mode switching in a bimodal microcavity. We presented experimental results and their explanation in terms of a transparent analytical theory. In particular we found that the transition has to occur via an intermediate phase where both modes are selected. Our theoretical description reveals, moreover, a close connection to the physics of equilibrium Bose-Einstein condensation in quantum gases of massive bosons. The mode switching can, there-fore, be viewed as a minimal instance of Bose-Einstein condensation of photons and its demarcation to lasing.
We also investigated the statistical properties of the system and pointed out that the mode switching is accompanied by superthermal intensity fluctuations as well as anticorrelations between both modes. This observation can be technically relevant, since a device producing a drastically increased occurrence rate of photon pairs (with a very narrow linewidth [45]) could be used to enhance two-photon excitation process used in fluorescence microscopy [70]. Furthermore a device that changes its predominantly emitting mode in dependence of the injection current could be applied for optical memories and other types of mode management [12,71].

ACKNOWLEDGMENTS
The authors thank A. Musia l for very helpful comments on the manuscript. TL, DV, and HAML contributed equally to this work. DV is grateful for support from the Studienstiftung des Deutschen Volkes. We acknowlege funding from the European Research  The asymptotic theory describes the generic form of the mode switching and its analytic expressions can thus be used to obtain the parameters of the master equation model. However, the theoretical parameters cannot be related directly to experimental parameters due to the unknown proportionality factor a i between the intensity of the emitted light J i and the occupation of modes, n i = a i J i , and the unknown excitation efficiency b of the pumping with respect to the injection current, P = bI.
The main properties of the switching are captured by the effective-gain ratio G which determines whether a switching occurs and the extent of phase C, P CD P BC = G −1 [see Eq. (27)]. This ratio can be obtained in the following way: Apply first a linear fit for the intensities of the selected modes i in each of the phases B, C, D, The ratio P CD P BC = I CD I BC is then determined either by the in- The parameters A h→l , g h , g l , h , l , are extracted for comparison between theory and experiment via the least squares method for all experimental data with I < 80 µA and are listed in the caption of Fig. 2. Since the time scale in Eq. (2) does not affect steady state properties, all parameters are measured in units of τ . The individual rates R l→h and R h→l do not affect the asymptotic theory, only the rate asymmetry A l→h does [as discussed above]. But the correlation function g (2) ll does depend on the individual rates, so that R l→h is chosen to reproduce this correlation function.
In the experiment the orientation of a polarization filter is chosen parallel to the passive cavity modes at the inversion point. Due to the interactions induced by the common gain medium the polarization-resolved spectrum exhibits a double peak structure, indicating that each polarization direction contains small portions of the other mode [45]. To make numerical and experimental fluctuations comparable we take into account that in each polarization a small fraction of the other mode is detected by introducing the mixing where n meas h,l denotes the quantity that is measured. The mode mixing prevents the experimental observed fluctuations of the non-selected mode from increasing monotonically in phase D [see Fig. 2(b)]. Comparing the theoretical results with the experimental data, we find the optimal value of the mixing parameter to be very small, c = 8 · 10 −3 . Its impact of the mixing is negligible in phases A, B, and C. It is taken into account in the numerical data presented in Figs. 2 and 8, where the index 'meas' is dropped.

Appendix B: Reduced density matrix
Under the (idealizing) assumption that all spontaneous emission goes into the selected modes i.e. τ −1 = 0 and β = 1 a reduced density matrix of the form ρ n h ,n l can be derived. To derive the equation of motion for the reduced density matrix we need to consider only the parts of Eq. (2) describing the pump and the photon emission into the cavity Here and in the following equation [. . . ] stands for terms that describe describe the loss of the individual modes as well as intermode transitions. When no spontaneous emission is lost into non-selected modes the carrier recombination and emission into the cavity modes is faster than any other process so that the term ij obtained from the numerical exact Monte-Carlo simulation of Eq. (2) (solid lines) with the ones obtained from the direct solution of the reduced density matrix Eq. (B1) (dashed lines). Note that the pump in the reduced density matrix is scaled to compensate that β = 1.
i g i N (n i + 1)ρ n N can be substituted by P ρ n N −1 . This means that whenever an emitter is excited by the pump its excitation is immediately emitted into the cavity, thus the emission into the modes can be described directly by P ρ n N −1 [35]. By the same reasoning or by simply shifting the indices (ρ n−e i N +1 ) one can find a substitute for all terms in Eq. (2) that correspond to the photon emission. Now we can trace over the emitter subspace ( N ), resulting in an equation of motion for the reduced density matrix, d dt ρ n = −P ρ n + i P g i n i ρ n−e i j g j (n j + 1) − g i + [. . . ].
As argued above, the reduced density matrix approach works under the assumption of β = 1. In our case β 0.2, so only a fraction of the pump effectively creates photons in the cavity. To be able to compare the reduced density matrix to the full model the pump power is scaled accordingly. Figure 7 shows the results of the Monte-Carlo simulations of the full equation compared to the results obtained from numerical solution of the reduced equation. The deviation of the reduced density matrix approach for small pump rates is not a problem since for low pump rates the full equation can still be solved numerically exactly (as it is done for Fig. 4 A and  B). Importantly, the reduced equation reproduces the results in the regime of high pump rates, where the exact solution for the full equation can no longer be obtained. ii ≤ 2. All parameters are chosen as in Fig. 2.

Appendix C: The role of spontaneous intermode transitions
In this appendix we investigate the impact of the effective parameter s, which quantifies spontaneous intermode transitions. Figure 8 shows the mode characteristics for the case with full spontaneous transitions between the modes (s = 1). In contrast to Fig. 2, the sharp kinks in the occupations [panel (a)] and the cross correlation [panel (c)] in phase C are less pronounced. However, the most significant deviation appears in the photon autocorrelations. For s = 1 the computed autocorrelation does not reproduce the experimentally observed superthermal fluctuations, g (2) ≤ 2.
This numerical observation can be backed analytically using the following argument. In phase B and D the joint distribution ρ n N factorizes approximately into two parts, one describing the non-selected mode i and the other one describing the selected mode j and the emitters, ρ n N ≈ ρ ni ρ nj N . This statement is confirmed by the numerical solution of the full master equation. Such a factorization allows one, to trace out both the selected mode and the emitters, to obtain an equation of motion for the photon number distribution of the non-selected mode: d dt ρ n = − N g[(n + 1)ρ n − nρ n−1 ] − [nρ n − (n + 1)ρ n+1 ] −R j→i n j (n + s)ρ n − (n − 1 + s)ρ n−1 −R i→j ( n j + s) nρ n − (n + 1)ρ n+1 .
Here, n denotes the photon number of the non-selected mode, n j , N the mean occupations of the selected mode and the emitters respectively, R j→i (R i→j ) the transition rate from (to) the selected mode j to (from) the non-selected mode, and g, the gain and loss rate of the non selected mode. If spontaneous transitions are fully included, s = 1, Eq. (C1) is solved by a distribution of the form ρ n = (1 − α)α n , which always yields g (2) = 2. This explains why we do not find superthermal fluctuations for s = 1.