Bimodal behavior of microlasers investigated with a two-channel photon-number-resolving transition-edge sensor system

We explore the photon-number distribution of bimodal quantum-dot micropillar lasers with a two-channel transition-edge sensor (TES) detection system. The two channels of the photon-number-resolving TES system simultaneously detect light emission of two orthogonal components of the micropillar’s fundamental emission mode. The applied experimental scheme provides unprecedented access to the joint photon-number distribution and enables a profound insight into the dynamics and photon statistics of the gain-coupled mode components. In particular, the two-channel TES measurements reveal an optical bistability of the correlated laser modes leading to temporal hopping between emission associated with Poissonian and thermal-like emission statistics. The experimental data and theoretical modeling based on Monte Carlo simulations are in good agreement and reveal the anticorrelated behavior of the mode hopping, which results in intensity ﬂuctuations and superthermal values of the autocorrelation function. Our investigations clearly demonstrate the great beneﬁt of using photon-number-resolving detectors in nanophotonics to explore the rich physics of multimode micro- and nanolasers.


I. INTRODUCTION
The enormous progress in semiconductor laser miniaturization has opened up new research opportunities at the crossroads between optoelectronics, nanophotonics, quantum optics, and nonlinear dynamics [1][2][3]. Interesting concepts are, e.g., coupled bimodal photonic crystal cavities suitable for studying nonequilibrium thermodynamics [4] and their applications [5].
Lasers based on microcavities with high quality factors (Q factors) and small mode volumes, which operate in the cavity quantum electrodynamics (cQED) regime [6], are of particular interest. Such high-β micro-and nanolasers enable thresholdless lasing [7][8][9] and even lasing of single quantumdot (QD) devices [10][11][12]. QD micropillar lasers, which have been investigated in miscellaneous studies [13][14][15], often show a bimodal behavior with two modes with orthogonal linear polarizations competing for the common gain medium [16,17]. This peculiarity can lead to intriguing effects such as unconventional normal mode coupling [18] and stochastic temporal hopping between the fundamental mode components above the laser threshold [19], which influences the intrinsic mode dynamics and photon statistics [16]. The photon statistics of emission can be accessed by determining the second-order photon autocorrelation function which provides information about the nature of emission and temporal intensity fluctuations. Important additional information can be obtained by examining the photon-number distribution (PND) of the emission, which not only provides insight into the second-order photon correlation function, but also accesses higher-order correlations.
Quantum optical measurements using the Hanbury Brown and Twiss (HBT) configuration [20] have become an important tool for studying the photon statistics of microlasers in terms of the second-order autocorrelation function [21][22][23]. However, the typically used click detectors based on avalanche photodiodes or superconducting nanowire singlephoton detectors do not provide access to the underlying PND. In fact, these detectors cannot distinguish between one or more photons in a light field due to their Geiger-modelike operation principle. To overcome this issue, the input signal was demultiplexed to multiple detectors in a cascaded arrangement to mimic photon-number-resolving (PNR) properties in Ref. [24]. More suitable for this task, however, are superconducting transition-edge sensors (TES) with true PNR capability which enables the direct photon-number measurement in a light pulse with a dynamic range of up to 20-30 photons [25][26][27].
So far, TES detectors with sensitivity in the optical to near-infrared spectral range have mainly been used to study light sources based on parametric down-conversion (PDC) processes [28,29]. More recently, their application has been extended to the characterization of nanophotonic devices including semiconductor QDs as photon emitters [30,31]. Furthermore, a single-channel TES detector has been used for the PNR-resolved study of bimodal micropillar lasers. Despite the restriction that only one emission mode could be investigated at a time, this experiment provided important insights into the physics of such lasers and revealed a double-peak characteristic of the individual PNDs as a strong indication of a temporal mode hopping process [32].
To achieve a deeper understanding of the underlying emission properties and to go beyond the simple study of individual PNDs, in this article we use two TESs to directly measure the joint PND of two modes in an electrically driven bimodal microlaser. Measuring the joint PND enables us to distinguish with unprecedented clarity situations and bias points where the two modes emit independently with different (thermal-like and coherent) photon statistics or jointly with mixed photon statistics. Moreover, it also allows us to demonstrate the power of the advanced TES detection scheme which may stimulate further experimental and theoretical studies of coupled photonic and quantum-photonic systems. Interestingly, for the studied bimodal microlaser the gain competition between the two modes leads not only to a mode hopping but also to a mode crossing (i.e., an exchange of mode intensities) at high injection strength well above the laser threshold. The latter can only fully be revealed by the PNR technique and is explained by an extended two-state model that describes the hopping between emission dynamics that are either associated with Poissonian or with thermal-like statistics. The experimental results are supported by Monte Carlo simulations.
We describe our sample technology and the experimental setup in Sec. II, while in Sec. III the results are presented in terms of input-output characteristics, correlation functions, and PNDs. The two-state model is explained in detail here. In Sec. IV we introduce Monte Carlo simulations that allow for a direct comparison with the experimental results and a conclusion is given in Sec. V.

II. SAMPLE TECHNOLOGY AND EXPERIMENTAL SETUP
The experiments were carried out on an electrically driven micropillar laser which is based on a planar AlAs/GaAs heterostructure with a single layer of In 0.3 Ga 0.7 As QDs as the active medium. The QD layer with an areal density of 5 × 10 9 /cm 2 is placed in the middle of the central one-λ GaAs cavity which is sandwiched between an upper and lower distributed Bragg reflector (DBR). The lower (upper) DBR is composed of 30 (26) n-doped (p-doped) λ/4-thick AlAs/GaAs mirror pairs, where the doping profile was optimized to achieve a good balance between high conductivity and low free-carrier absorption [33]. After the epitaxial growth of the semiconductor heterostructure by molecular beam epitaxy, electron beam lithography (EBL) and plasma etching were used to fabricate micropillars with a diameter of 4 μm. The sample was then planarized with benzocyclobutene, and ring-shaped gold contacts that form the top (p-type) electrode of each micropillar were realized via a second EBL processing step [33].
To study the emission properties and the joint PND of the bimodal QD-micropillar lasers, the sample was mounted onto the cold finger of a continuous-flow helium cryostat set FIG. 1. Sketch of the used microelectroluminescence (μEL) setup: The bimodal micropillar is operated at T = 20 K inside a continuous-flow helium cryostat. The two linearly polarized fundamental mode components are separated with a half-wave plate (HWP) and a polarizing beam splitter (PBS). The spectral properties are measured with two grating spectrometers each equipped with a charge-coupled device (CCD) camera. The joint PND is measured using two fiber-coupled transition-edge sensors (TESs) read out by a superconducting quantum interference device (SQUID) current sensor. The TES detectors are placed in an adiabatic demagnetization refrigerator (ADR) and are operated at a temperature of 100 mK. The laser excitation is synchronized with the data acquisition (DAQ), which digitizes the analog TES output signal for further processing.
to a temperature of T = 20 K as illustrated in Fig. 1. The microlasers were driven by a pulse generator with a maximum voltage amplitude of V pulse,max = 5.1 V. For all measurements a pulse length of t pulse = 2 ns was used. Additional, a DC voltage of 1.7 V was applied to operate the microlasers close to the onset of electroluminescence for zero pulse amplitude. Overall, this defines the bias voltage as V bias = V dc + V pulse .
The light emitted by individual micropillar lasers was collected by an aspherical lens (numerical aperture NA = 0.5) in front of the cryostat. With a motorized half-wave plate and a polarizing beam splitter (PBS), the two orthogonal modes are selected and their emission is directed via different paths to two separate monochromators, each equipped with a chargecoupled device camera. In order to obtain information about the joint PND of the orthogonal emission modes, two highly sensitive single-mode fiber-coupled TES detectors were used. They were attached to the output slit of each of the two monochromators and enabled us to measure the number of photons in each incoming optical pulse [34].
The absorption of individual photons results in an increase of temperature in the TES absorber material (tungsten). Simultaneously, the absorber is also a highly sensitive superconducting phase transition thermometer, which is able to detect photons with a resolution at the single-photon level. Each TES is voltage biased to set the working point in the phase transition between superconducting and normal conducting state at a temperature of ∼150 mK. Ultimately, the photon absorption results in a change of current which is measured by an inductively coupled two-stage direct-current superconducting quantum interference device (DC-SQUID) current sensor [35]. The two fiber-coupled tungsten TESs are placed in an adiabatic demagnetization refrigerator. The temperature is regulated at 100 mK and the temperature fluctuation is 15 μK (root mean square). The detection efficiency of both channels was determined to be 87% [36]. Because of the long thermal recovery time (about 1 μs) of the signal response after photon absorption the TES cannot be used in continuous-wave mode. Therefore, the repetition frequency of the pulse generator used to excite the micropillar lasers is set to f = 5 kHz and its synchronization output is used to generate the time stamp for each joint TES detection event.
As discussed in detail in Sec. III we determined the joint PND of the microlaser in dependence of the pulse amplitude. Each corresponding data point discussed in the next section in Fig. 4 is based on the measurement and analysis of 200 800 16-μs-long TES time traces which are the digitized (125 Megasample/s) SQUID output signal (Fig. 2). For each TES time trace the pulse area can be extracted, which is assigned to a photon number in order to finally obtain the joint PND.

A. Spectral characteristics and input-output dependence of the bimodal micropillar laser
In order to explore the joint PND of a bimodal QDmicropillar laser, we first selected a device which shows suitable energy splitting of the fundamental emission mode and pronounced laser action. The fundamental mode (HE 11 ) components of the selected 4-μm-diameter pillar laser emit at an energy of about E 0 = 1.4595 eV and show an energy splitting of 30 μeV as can be seen in the μEL emission spectrum above laser threshold presented in Fig. 3. Here, the mode splitting results from a slight (unintentional) ellipticity of the pillar's cross section [11]. The low-energy mode component features a high emission intensity at intermediate pulse amplitudes between 3.0 and 4.5 V. Henceforth, we refer to this mode as the strong mode (s) and to the orthogonal mode as the weak mode (w). The two modes feature Q factors of Q = E s,w /γ s,w = 40 112 and 30 173 close to the onset of laser emission at 2.9 V, where E s,w (γ s,w ) denote the emission energy (linewidth) of the strong and weak mode. The linewidth was determined by fitting the experimental data using a Voigt line shape with a fixed Gaussian contribution with a width of 80 μeV representing the spectral resolution of the used setup. The same procedure was applied to determine the excitation- power-dependent data presented in Fig. 4. The input-output dependence of the bimodal microlaser is presented in Fig. 4(a) and shows the typical smooth s-shape expected for a highβ microlaser [6]. The two modes have a similar threshold voltage V pulse ≈ 2.7 V. Increasing the pulse amplitude beyond 4.5 V leads to a change in the emission characteristics with the weak mode becoming the superior mode above the intensity crossing point at about 4.85 V. The present input-output characteristics of the strong mode and the weak mode as well as the occurrence of the intensity crossing point are attributed to pronounced gain coupling via the joint QD gain medium. This coupling behavior is typical for bimodal micropillar lasers with a QD gain medium and has already been reported in a number of studies [16,17,37,38]. We observe a pronounced excitation-power-dependent linewidth reduction for both the strong and weak mode, which starts in the laser threshold range, as shown in Fig. 4(b). This behavior suggests the beginning of stimulated emission associated with increased temporal coherence for both modes. However, only measurements of the photon statistics can unambiguously reveal the nature of emission.

B. Photon-number distribution of the bimodal micropillar laser
Before obtaining more detailed insights into the underlying emission statistics and switching of the bimodal QD micropillar laser, let us introduce the necessary theoretical framework in terms of the joint PND. We consider a matrix P i j , in which the index pair (i, j) denotes the photon number i of the strong mode and the photon number j of the weak mode with the appropriate sets {i} = {0, . . . , m} and { j} = {0, . . . , n}, where m (n) represents the highest number of photons measured for the mode s (w). The m × n matrix P i j can be understood as a probability distribution if it has the following properties: (i) All elements are non-negative P i j 0 and (ii) the sum over all elements is unity i j P i j = 1. The corresponding singlemode PNDs are formed by summing up over the respective columns and rows: P i = j P i j and Q j = i P i j . In Fig. 5, the joint PND and the single PNDs for the strong mode P i and the weak mode Q j are shown for a pulse amplitude of 5.0 V just above the intensity crossing point (cf. Fig. 4 and cf. Sec. IV for more PNDs). Here, a noticeable feature of the measured joint PND is the concentration of probability values indicated by the darker red color. This feature is also evident in the double-peak structure of the individual PNDs in the histogram presentation of Fig. 5. The observed behavior will be explained and discussed in Sec. III D.

C. Correlation functions of joint photon-number distributions
The availability of the joint PND measured by our twochannel TES detection system allows us to extract the normalized moments of emission according to [39][40][41] wheren s =b † sb s andn w =b † wb w are the photon-number operators withb s (b † s ) andb w (b † w ) being the bosonic annihilation (creation) operators of the strong and the weak mode. The pair of non-negative integers (u, v) denotes the correlation order and :: indicate the normal ordering of operators. We only refer to the zero-delay functions g (u,v) s w (τ = 0), and from now on the time argument is omitted. The associated factorial moments are given by In Eq. (1), the expressions :n 1 sn 0 w : and :n 0 sn 1 w : are the mean photon numbers of the strong n s and the weak mode n w . Furthermore, if the correlation order of one mode is set to zero, the autocorrelation of the other single mode is re- w . The correlation functions defined in Eq. (1) are non-negative and they probe the (u + v)particle correlations. For the convenience of the reader, we connect the general correlation functions to the well-known second-order autocorrelation function, for example, for the strong mode, by inserting Eq. (2) into Eq. (1), as well as the second-order cross-correlation function Both correlation functions are in principle accessible via a classical HBT arrangement with two click detectors, although this would still hide the full photon statistics. Figure 4(c) presents excitation-power-dependent g (2) values of the strong mode and the weak mode for the bimodal QD-microlaser under study. Below threshold the coherence time of such lasers is on the order of a few tens of ps [42], which is much smaller than the applied minimum pulse length (2 ns) of the voltage source. For this reason it is not possible to resolve the expected bunching behavior with g (2) = 2 below threshold [22] and we focus our attention on the bias region at and above threshold. At threshold both modes show g (2) close to unity, indicating the emission of coherent light. By increasing the voltage amplitude, the autocorrelation signal of the strong mode increases slightly. Up to higher excitations a strong increase starts (V pulse ≈ 4.0 V) until a value just below g (2) = 2 is reached at maximal voltage amplitude. In the case of the weak mode, g (2) increases much stronger with the excitation strength and exceeds the value of thermal emission in the voltage range 3.6 V V pulse 4.5 V. In contrast to the strong mode's behavior, the g (2) value of the weak mode decreases at high excitation strength and reaches values below 2, which might indicate the transition towards coherent emission at very high excitation. Unfortunately, the very high excitation range is not accessible because it lies beyond the maximum pulse amplitude available in our experiment. Interestingly, the crossing of g (2) occurs at the same pulse amplitude as for the intensities and shows that the two modes exchange their predominantly coherent and thermal character, respectively, at this bias point. The superthermal two-photon bunching at intermediate excitation strengths of the weak mode is a consequence of mode competition [16,19,43] and will be discussed within a two-state model in Sec. III D. Superthermal light can also be generated in bimodal microlasers prepared far from equilibrium through a quench induced by a short pulse [4] and in superradiant microlasers [44,45].
The same conclusions can be drawn for the third-order autocorrelation function g (3) shown in Fig. 4(d). The maximum value for the weak mode reaches values well above the value of 6 for thermal light. Interestingly, the fourth-order autocorrelation function g (4) for the weak mode does not exceed the value for thermal light (here 24) as can be clearly seen in Fig. 4(e). Hence, the light emitted from our bimodal microlaser exhibits a two-photon and three-photon superthermal bunching but the bunching of four photons is below that of thermal light. This result clearly shows that the property "su-perthermal bunching" always has to be considered in relation to the considered number of photons.
The key strength of our experimental configuration using a two-channel TES detection system lies in the ability to measure not only the individual PNDs of the two modes, leading, for instance, to g (2) , but also the joint PND which gives access to the cross-correlation functions g (u,v) .
Let us first consider the limiting case of uncorrelated modes A and B. From B . If A and B exhibit thermal characteristics, this leads to g (1,1) A In fact, any uncorrelated distribution leads to g (1,1) A B = 1. This is in stark contrast to highly correlated sources such as those based on PDC [29] for which values of g (1,1) = 100 000 have been demonstrated experimentally [46]. Figure 4(f) shows the symmetric cross correlations g (1,1) s w and g (2,2) s w . At threshold the two-particle correlation function g (1,1) s w is close to unity corresponding to uncorrelated mode intensities. With increasing pulse excitation the cross correlation decreases below unity, indicating an anticorrelation of the two mode intensities. The cross correlation reaches a minimum value of g (1,1) s w ≈ 0.27 at V pulse ≈ 4.7 V, which is near the crossing point V pulse ≈ 4.85 V where g (1,1) s w ≈ 0.29. These values are significantly lower than the minimum value of 2/3 derived in the mesoscopic limit for the cross correlation in coupled nanolasers performing limit cycles [47]. However, the physics is slightly different and our system, in the terminology of Ref. [47], is in the nanoscopic regime, not in the mesoscopic one. In the Appendix we derive for our system an approximation of the cross correlation at the crossing point n s = n w , where we have assumed that the total photon number (sum of photons in mode s and w) have the same statistical properties as in a single-mode laser. With the experimental values g (2) s ≈ 1.74 and g (2) w ≈ 1.69 we get from Eq. (5) g (1,1) s w ≈ 0.28, which is in very good agreement with the measured value.
The four-particle correlation function g (2,2) s w measures the cross correlations of the intensities squared. This correlation function shows a similar behavior, with the difference that g (1,1) s w > g (2,2) s w over the entire excitation range. Here, a minimum value of g (2,2) s w ≈ 0.15 at V pulse ≈ 4.6 V is obtained. In order to better understand the results in relation to the circumstance of cross correlation in the photon-number representation different correlations can be considered. For example, from Eq. (4), g (1,1) can be regarded as the "cross mean value" divided by the two individual mean values. It can be seen that the probabilities near the diagonal i = j cause large contributions in the numerator of Eq. (4). If we look at the joint PND in Fig. 5, these probabilities are particularly low and validate the anticorrelating behavior. Figure 4(g) shows the asymmetric cross correlations g (2,1) s w and g (1,2) s w . The values also start close to 1 and decrease with increasing excitation. The crossover behavior is also visible at V pulse ≈ 4.85 V which implies an exchange of the mode properties. Similar to Eq. (5) we derive in the Appendix an approximation of the asymmetric cross correlation at the crossing point with g (1,2) s w ≈ g (2,1) s w , g (1,2) s w = 1 With the experimental values g (3) s ≈ 3.08 and g (3) w ≈ 3.25 [ Fig. 4(d)] follows g (1,2) s w ≈ 0.28, which is very close to the experimental result in Fig. 4(g).

D. Two-state model for the description of joint photon-number distributions
As the individual PNDs of the strong and the weak mode depicted in Fig. 5 seem to feature characteristics of thermallike as well as Poissonian-like distributions, they suggest that both of the modes are in a state of superposition of different statistics. To gain a deeper understanding of this superposition and identify its individual parts, we will make use of a twostate model as well as a decomposition of the joint PND P i j .
The intuitive idea of the two-state model for bimodal lasers is that in the time domain the two modes frequently hop between two states (l = 1, 2). In state 1 the strong mode is lasing and the weak mode is nonlasing. In state 2 it is the other way round. This model has been deduced from quantummechanical calculations of second-order correlation functions in the steady state [16,48] and, more directly, from semiclassical calculations of second-order correlation functions in the time domain [19].
The two-state model was lifted to the joint PND in Ref. [49] by assuming the incoherent superposition The non-negative numbers a l with a 1 + a 2 = 1 are the mixing parameters. a l can be considered as the probability to find the two modes in state l. P (l ) i (Q (l ) j ) is a probability distribution for the strong (weak) mode in the state l. An interpretation of each P (l ) i as a probability distribution requires (i) i P (l ) i = 1 and (ii) P (l ) i 0 for all i. The same is required for Q (l ) j . The resulting individual PND P i (Q j ) is an incoherent superposition of P (1) i and P (2) i (Q (1) j and Q (2) j ). In Ref. [49] P (2) i and Q (1) j were assumed to be thermal distributions and P (1) i and Q (2) j to be Gaussian distributions (close to Poisson distributions), which, at first glance, is also suggested by Fig. 5. We find, however, that this assumption does not provide satisfactory results over the complete range of parameters. For a more accurate decomposition of the hopping process, we employ the non-negative matrix factorization (NMF) [50] to numerically determine P (l ) i and Q (l ) j without further assumptions. Consider an m × n matrix V which is non-negative, i.e., all matrix elements are non-negative. For a given positive integer K < m, n we approximate the matrix V by a low-rank factorization with rank K, non-negative m × K matrix W , and non-negative K × n matrix H. The matrices W and H are chosen to minimize the root-mean-squared residual between V and W H.
Note that the solution is in general not unique as the optimization procedure yields only a local minimum. The relation  (7) and the non-negative matrix factorization in Eq. (8). In (e) the single distributions of strong (P (1) i ) and weak (Q (1) j ) modes for l = 1 and in (f) the single distributions of strong (P (2) i ) and weak (Q (2) j ) modes for l = 2 are shown.
to the decomposition of the joint PND in Eq. (7) becomes obvious by the identifications P i j = V i j , K = 2, and Using standard software packages we are now able to decompose the joint PND P i j and thus obtain the individual distributions of the strong and the weak mode for the respective state. As an example, the data from well attributed to a respective state. In state 1 the strong mode (P (1) i ) is close to a Poissonian distribution with a corresponding g (2) value of 1.01 [ Fig. 6(e)] and the weak mode (Q (1) j ) associated with a thermal-like distribution with g (2) = 2.39. In state 2 it is the other way around. The strong mode (P (2) i ) is associated with a thermal-like distribution with a g (2) value of 2.47 [ Fig. 6(f)] and the weak mode (Q (2) j ) is close to a Poisson distribution with a corresponding g (2) = 1.01. The mixing parameter a 1 (a 2 ) is 0.448 (0.539). Hence, here the state 2 is slightly more likely than state 1. The sum a 1 + a 2 = 0.987 is close to unity, which is a necessary condition for the applicability of the rank-2 approximation (7) and (8). If we compare Figs. 6(c) and 6(d), we identify small differences between the reconstructed distribution and the original distribution in the region between the maxima. This intermediate region cannot be attributed uniquely to state 1 or 2. It is therefore consistent that this fine structure of the original PND is not captured by the two-state model.
Next, we explore how the properties of states 1 and 2 depend on the voltage pulse height. At laser threshold and below, the NMF is not applicable, because one gets into the region of K ≈ m, n. Figure 7 therefore only shows voltages V pulse 3.5 V. In Fig. 7(a) one can see that slightly above threshold, the mixing parameter a 1 is considerably larger than a 2 . Hence, state 1 is dominant, i.e., with a high probability the strong mode is lasing and the weak mode is nonlasing. Far above threshold there is a crossing point a 1 = a 2 at V pulse ≈ 4.85 V exactly where the crossing point in Figs. 4(a), 4(c), 4(d), 4(e), and 4(g) is. Above this crossing point, the mixing parameter a 2 is larger than a 1 , i.e., the state 2 overcomes state 1, which implies that with higher probability the strong mode is nonlasing and the weak mode is lasing. Note that throughout the shown range of voltage pulse heights, the sum of the two mixing parameters is close to unity, which is an indication that the rank-2 approximation is valid. Figure 7(b) shows the second-order autocorrelation functions of the extracted PNDs P (1) i , P (2) i , Q (1) j , Q (2) j . As mentioned in the previous section, both modes reach the lasing regime for different states l (state 1 for the strong mode and state 2 for the weak mode) in which they have a g (2) value of unity above the laser threshold. The g (2) values of the strong mode in state 2 and the weak mode in state 1, on the other hand, describe their nonlasing components. The component s2 (strong mode in state 2) rises more slowly from a g (2) value of 1 to a value of 2 than component w1 (weak mode in state 1). From V pulse = 4.6 V there is a convergence of the two cases and from 5.0 V they decrease together to g (2) = 2. For V pulse ≈ 4.6-5.1 V both modes are above a value of g (2) = 2, the value for thermal light. It is remarkable here that even after the decomposition two combinations can have a g (2) value greater than 2, indicating superthermal light emission. This is consistent with the above-mentioned fact that an adaptation of the single distribution by a Poissonian and a pure thermal distribution does not provide satisfactory agreement.

A. Birth-death model for the bimodal laser
The theory is based on the birth-death model for the bimodal laser [16,17,32]. The phenomenological model deals with the diagonal elements of the quantum mechanical density matrix n, N|ρ|n, N = ρ n N , giving the probability to find the system with N excited emitters and n = (n s , n w ) photons in the strong and weak mode, respectively. The master equation is The vector ±e i denotes an added or removed photon in mode i, i.e., e s = (1, 0), e w = (0, 1). Matrix elements ρ n N with negative indices are zero. The first term describes the pump process which increases N by one; P is the pump rate. The second term represents the spontaneous emission into leaky and nonresonant modes reducing N by one; τ −1 is the spontaneous emission rate into leaky modes. The third term takes care of the emission into the mode i, where n i is increased by one and N is reduced by one; g i is the emission rate into the mode i. The fourth term corresponds to the loss of photons from mode i reducing n i by one; κ i is the optical loss rate of mode i. Finally, the last term accounts for the mode coupling where n i is reduced by one and n j is increased by one. Note that the transition rate R i→ j can be different from R j→i . This asymmetry can result from stimulated scattering due to carrier population oscillations [51] and geometric asymmetry of the cavity [52]. We mention that our model ignores the free parameter s that expresses a possible modification of the spontaneous emission between the two modes by the gain-medium-induced mode interaction [32].
The master equation (11) can be numerically solved, for instance, by direct time integration of the system of differential equations [16,32], by the fixed-point iteration [53], or by Monte Carlo simulations [17]. The latter is usually the most efficient scheme and allows us to compute the steady state and the full dynamics. The Monte Carlo approach for bimodal lasers in Ref. [17] employs the Gillespie algorithm [54,55]. The state space is given by N and n = (n s , n w ) with zero as initial condition. Each process in the master equation (11) is realized by probabilistic jumps in the state space. The corresponding jump probabilities and the variable dwell time are determined by the weights in front of ρ n N in Eq. (11): P, τ −1 N, g i N (n i + 1), κ i n i , or R i→ j n i n j . The statistics of N and n can be used to approximate the time evolution of ρ n N .

B. The extended Monte Carlo simulation
In all previous approaches to describe the emission features of microlasers based on the birth-death model [16,17,32,53] the steady-state probability ρ n N was computed. The steadystate solution was used to determine photon numbers, correlation functions, and single-mode PNDs. The comparison of these internal quantities (e.g., number of photons inside the cavity) to the externally measured experimental data (e.g., number of photons finally being collected by the detector system) requires in principle a detection model.
For photon numbers the detection model is trivial; only multiplication with the corresponding optical loss rate κ i and the total collection efficiency η are required. For normalized correlation functions such as g (2) (0), it seems that further treatment is not necessary, as the multiplied factors κ i and η cancel out. However, we have seen already in Sec. III C that below the laser threshold, the theoretical g (2) (0) (treated as an internal quantity) does not agree with the (externally) measured g (2) (0). The difference between internal and external quantities disappears only if the coherence time is larger than the detector's temporal resolution. In the opposite case, the theory has to include the finite time resolution of the detector, for instance, by a convolution with an appropriate apparatus function [22], which corresponds to a detection model.
We note that the need for a detection model is most obvious for PNDs. It can be implemented by simulating the leakage of photons in a second calculation after the steady state ρ n N is reached [17]. This, however, does not remove the detrimental influence of a short coherence time.
Our approach avoids all of these problems by (i) simulating the transient dynamics rather than the steady state and (ii) by including the external degrees of freedom from the very beginning. The latter are the numbers of photons (n collected,s , n collected,w ) collected by the detector system. Hence, we have extended the state space used in Ref. [17] to (N, n s , n w , n collected,s , n collected,w ). All numbers are set to zero initially. As a conditional process in the Gillespie algorithm, the number of collected photons n collected,i is increased by one with probability η whenever n i is decreased by one due to optical losses. The statistics of the accumulated (n collected,s , n collected,w ) immediately after the pulse directly approximates the joint PND we are looking for. Hence, the detection model is already included in the system dynamics. There is no need to model the finite collection efficiency or the finite measurement time afterwards. And, most importantly, there is no need to include the correlation time from other sources. The correlation times (actually several correlation times of various orders and for both modes) enter here exactly as in the experiment, namely by accumulating photons collected at different times.

C. Comparison to the experimental data
For a comparison with the experimental data we use the following parameters: A total collection efficiency η = 3.4 × 10 −4 , a rectangular pulse of length t pulse = 2 ns, a spontaneous emission time into the leaky modes τ = 1.33 ns, Q factors of 44 500 and 37 600, and β factors of 0.25 and 0.23 for the strong and the weak mode, respectively. The optical loss rates are then determined by κ i = γ i /h = E 0 /(Q ih ) and the rates of emission into the modes by The parameters for the mode coupling are fitted to get the best agreement with the input-output curves and the correlation functions, R 21 = 5 ns −1 and R 12 = R 21 + 0.06 ns −1 . The relation between the pump rate P in the birth-death model and the applied voltage V pulse in the experiment is not known. For simplicity, we assume a quadratic dependence P = a(V pulse − V 0 ) 2 , and a fit results in a = 3364 ns −1 V −2 and V 0 = 2.5 V. Figure 8 shows the emission intensities and the correlation functions of the collected photons as a function of the applied voltage. While there are deviations for the intensities at low voltages [ Fig. 8(a)], at intermediate and large voltages the theoretical results clearly reproduce the behavior of the strong mode and the weak mode. In particular, the exchange of the modes' properties at the crossing point V pulse = 4.85 V is accurately described. The autocorrelation functions in Figs. 8(b)-8(d) show a very good agreement between the experimental and theoretical results. It has to be emphasized that the theory is capable to reproduce the values near unity for low voltages. As explained in Sec. IV B, our theory contains inherently a detection model by incorporating external degrees of freedom. There is no need for a convolution with an apparatus function. The cross correlations in Fig. 8(e) are also reproduced albeit with a less impressive agreement. Figure 9 compares the experimental (left) and theoretical (right) joint PNDs using the same parameters for various voltages. The top row shows V pulse = 5.0 V (the experimental data have already been shown in Fig. 5) which is the case just above the crossing point of the two modes. A remarkable agreement between theoretical and experimental results can be observed. The next two rows show V pulse = 4.85 V close to the crossing point and V pulse = 4.7 V just below the crossing point. While the two-state structure with two pronounced maxima is preserved during this transition, the strong mode gains intensity if compared to the weak mode. The two-state structure gradually disappears when going to even smaller voltages. This is demonstrated in the two lower rows with V pulse = 3.8 V in the middle between the laser threshold and crossing point, and V pulse = 3.2 V slightly above the threshold. Also in this regime, the experimental and theoretical results agree remarkably well.

V. CONCLUSION
We developed a two-channel transition-edge sensor detection system which provides unprecedented access to the joint photon-number distribution of two correlated optical modes. To illustrate the capability of the detection scheme, we analyzed the light emitted from a bimodal quantum-dot micropillar laser. We experimentally demonstrated that our transition-edge sensor detection system allows for a profound insight into this bimodal quantum system as the measured photon-number distribution directly reveals a two-state structure. This two-state structure corresponds in the time domain to a stochastic hopping between two states. In state 1 the first mode is lasing while the second mode is nonlasing. In state 2 it is the other way around.
From the joint photon-number distribution one can extract various kinds of emission properties of the two modes. Here, one is not restricted to intensities, autocorrelation functions, and the cross-correlation function of second order as in conventional configurations using, e.g., simple click detectors, but one has access also to auto-and cross correlations of higher order. These higher-order moments are out of reach for conventional schemes.
We introduced a two-state model based on the non-negative matrix factorization. This factorization decomposes the measured joint photon-number distribution into two contributions, each representing one of the two states. By analyzing the second-order correlation functions of these two contributions separately we demonstrated unequivocally that, for sufficiently large excitation voltages, in a given state one mode is indeed lasing while the other one exhibits a superthermal two-photon bunching.
The experimental data and their interpretation is supported by theoretical modeling based on Monte Carlo simulations of a birth-death model. The presented approach goes well beyond known approaches by explicitly taking into account external degrees of freedom. This allows for a direct comparison of theoretical and experimental results. Good agreement between theory and experiment was observed for intensities, correlation functions, and (joint) photon-number distributions throughout the available excitation range.
We expect that in the future our photon-number-resolving detectors will prove to be a valuable tool for the investigation of various phenomena in the fields of optoelectronics, nanophotonics, and quantum optics. They may, for example, provide access to the photon-number distributions of superradiant systems which are not yet well understood [44,45].