Mode-switching induced super-thermal bunching in quantum-dot microlasers

The super-thermal photon bunching in quantum-dot (QD) micropillar lasers is investigated both experimentally and theoretically via simulations driven by dynamic considerations. Using stochastic multi-mode rate equations we obtain very good agreement between experiment and theory in terms of intensity profiles and intensity-correlation properties of the examined QD micro-laser’s emission. Further investigations of the time-dependent emission show that super-thermal photon bunching occurs due to irregular mode-switching events in the bimodal lasers. Our bifurcation analysis reveals that these switchings find their origin in an underlying bistability, such that spontaneous emission noise is able to effectively perturb the two competing modes in a small parameter region. We thus ascribe the observed high photon correlation to dynamical multistabilities rather than quantum mechanical correlations.


Introduction
Quantum-dot (QD) light sources based on high-Q optical microcavities are widely discussed to be future optoelectronic devices for data communication [1][2][3][4][5], and with regard to quantum communication, where they can act as sources of single photons or of entangled photon pairs. Micropillar lasers are a further step towards revolutionizing the field of laser devices which started with the development of the nowadays widely used vertical cavity surface emitting lasers (VCSEL), as there exist plentiful advantages over edge-emitting lasers e.g. in terms of ultra-low thresholds, circular beam profiles and on-wafer testing capability [6]. Interestingly most conventional VCSEL devices show polarization mode instabilities and polarization mode switching [7][8][9][10] which is a drawback in many applications, but also offers great opportunities regarding the observation of optical chaos in bench-top experiments [11]. The much smaller bimodal microcavity lasers discussed here are on the crossroad between conventional laser systems and single photon sources, and are singled out by their relatively large and tunable spontaneous emission coupling factors.
We explore the time dependent light emission dynamics of the bimodal micropillar lasers and provide the first experimental evidence of statistical switching reported for these kind of microscale devices. We clearly show that the nonlinear dynamics of these nanostructured devices can be described (qualitatively and quantitatively) with relatively simple rate equations adapted from the well known two-mode ring laser [12][13][14] and VCSEL [6,7,[15][16][17] models. Fully quantum mechanical theories that derive the equations of motion for photons and correlation expectation values [18,19] have their strength in describing single photon emitters [20] and have also been used to describe the micropillar lasers and their bimodal properties [21]. We show that they do not necessarily need to be applied here. Especially, the super-thermal bunching found in the micropillar lasers [21] can be understood in terms of dynamical instabilities. Thus, we pursue an extended semi-classical approach via stochastic two-mode rate equations that differs from existing models, as a common carrier reservoir and nonlinear gain-compression effects are included. We report on polarization switching in the transition region between classical nonlinear dynamics and nonlinear quantum optics. This includes a discussion of the bifurcation structure which shows that the underlying bistability is a crucial feature to explain the physical mechanism behind the intense multiphoton pulsations of the bunching laser mode. We compare our theoretical modeling with experimental results for an electrically driven bimodal QD micropillar laser.
Importantly, our description of the system includes cavity-quantum electrodynamically (cQED) enhanced light-matter coupling in the form of stochastic spontaneous emission noise. Thus, it is well suited to investigate the free running micropillar QD lasers under the influence of stochastic perturbations. This approach is used as the preferred method to achieve a simplified model for an excellent proportion between computational effort and physical insight. With our results we extend the scope of the ongoing vivid discussion on microlasers to the field of nonlinear dynamics and bifurcation studies.
In section 2 we will investigate the second order intensity correlation function ( ) ( ) t g 2 for analyzing the photon statistics of the light emission. This function can be conveniently calculated from simulations using time-series and is measured experimentally with the commonly used Hanbury-Brown and Twiss setup (HBT), as presented in section 4. In combination with experimental intensity profiles and correlation time measurements, experimental ( ) ( ) t g 2 measurements and energy spectra will be discussed in section 4. Section 5 compares these results with numerical simulations based on our rate equation described in from section 3. We discuss the findings in section 6 by singling out one possible mechanism for the super-thermal bunching found experimentally.

Photon statistics
CQED effects in low mode volume microlasers lead to an enhanced coupling of spontaneous emission into the cavity mode and result in high β values close to unity [22]. Here, the β-factor describes the fraction of spontaneous emission rate r laser coupled into the laser mode normalized by the total spontaneous emission rate [18,23,24]. Due to their Purcell enhanced β-factors, microcavities are characterized by high noise amplitudes on the one hand and lower losses into non-lasing and leaky modes [5,18] on the other hand. In the investigated microlasers, the high β-factor is of particular importance and differs significantly compared to large mode volume VCSELs. Due to its statistical character, the method of choice for investigating the light output and to reveal the underlying physics is to examine correlation functions. In contrast to ( ) ( ) t g 1 , which is the field correlation function, this paper will deal with intensity correlation functions ( ) ( ) t g 2 only. It can be written as [23,25]: where i and j denote the different cavity modes in our system. In the context of ring lasers the same quantity is discussed as normalized auto-and cross-correlation (l + 1 ij ) [14]. It is noted that above used definition is the classical representation of the quantum mechanical definition, which reads: In this equation, † b i and b i describe the photon creation and annihilation operators of one photon in mode i. Of particular importance is the case of t = 0 i.e. . The intensity correlation function determines the quantum state of the light field [26], namely it discriminates between thermal, coherent and non-classical light. Using the quantum mechanical photon distributions for Fock, coherent and thermal states, three prominent cases for

Model
The model used in this paper is based on semiclassical rate equations [27] taking into account the electron scattering mechanisms into the QDs as derived in our previous works [28,29] extended to the contribution of two optical modes. This is crucial as due to the cylindrical symmetry the microlaser exhibits two polarization modes, see sketch in figure 1(a). The mode degeneracy is lifted due to structural asymmetries and the fundamental cavity mode splits into two orthogonal linearly polarized mode components [30]. The amplitude of the complex electric field of weak and strong mode are denoted by E w and E s , respectively, where weak and strong corresponds to the modes with lower and higher lasing intensity. Both interact with the QDs via induced emission of radiation, see figure 1(b). Assuming only excitonic recombination processes the gain of both modes is proportional to the inversion of the QD system ( ) r -2 1 with the excitation probability for excitons ρ. The differential gain coefficients are labeled g s and g w for the strong and the weak mode, respectively. Because of the very low mode splitting of just a few tens of meV and the linear polarization of the modes, we assume only one carrier population that is interacting with both active optical modes. Hence, in contrast to the indirect coupling in the spin-flip model [7][8][9][10], the modes are directly connected via gain competition.
We introduce gain compression factors ε to the gain of both modes [31,32] which changes the differential gain according to , whereẽ e º n c 0 bg 0 is the conversion factor that translates | | E 2 into intensities. The cross-compression factors e ij are chosen to be non-zero, hence, high electrical fields in one mode will not only decrease induced photon emission in the same mode but in the second mode as well. The carriers are connected to an external reservoir n r , which is pumped by an injection current J. We assume an injection efficiency η and a parasitic current J p to account for the experimental details of the current injection. Carriers from the reservoir scatter into the dots with a constant scattering rate S in . It is noted that in the limit of small fields, degenerate modes and fast carrier dynamics our model can be rewritten into the two-mode rate equations used for Dye lasers with gain competition dating back to the 80s [12][13][14]. However, to capture a realistic pump current dependence the additional nonlinearities induced by the nonlinear gain and coupling terms are needed.
Leakage currents in the reservoir ( ) included. Spontaneous emission coupling into lasing modes is modeled two different ways in order to allow both bifurcation analysis and quantitative numeric simulation. For the latter we implement a Gaussian white noise For the bifurcation analysis the deterministic representation of the spontaneous emission is used, with an equal average number of photons that are emitted into the lasing modes [33]: The dynamic equations for the two slowly varying electrical field amplitudes E s , E w , the QD occupation probability ρ and the reservoir carrier density n r read: All parameters used for the simulations are given in table 1.

Experimental observations
The experimental studies were performed using electrically driven QD micropillar lasers. The devices are based on high-quality AlAs/GaAs microcavity structures with a single layer of self-assembled In Ga As QDs with an arial density of´-5 10 cm 9 2 in the active layer [34]. Using electron beam lithography and plasma etching micropillars with a diameter of 3.6 μm were patterned and electrically addressed using an upper ring-shaped contact. For more details on the layout and the fabrication process of electrically contacted micropillars we refer to [35]. We investigated the emission characteristics of the microlasers by means of micro-electroluminescence (μEL) spectroscopy using continuous excitation via an external current source. The experimental setup (see figure 2) has a spectral resolution of 25 μeV and includes a He-flow cryostat to operate the microlasers at cryogenic temperature (10 K). It is further equipped with a fiber-coupled HBT setup with Si based single photon counting modules in combination with a quTAU time-to-digital converter and a streak camera for recording single-shot timetraces.
A typical μEL emission spectrum is depicted in figure 3 at an injection current of 96 μA. The two emission lines of the micropillar correspond to the two orthogonal linearly polarized components of the fundamental HE 11 mode with a Gaussian far-field pattern. Here, the fundamental cavity more is split into two components, because of a slight structural asymmetry [30]. They are indicated by strong mode (black line) and weak mode (orange line) in the following. The mode splitting amounts to 103 μeV. The quality factors (Q) of strong and weak mode are determined by Lorentzian line shape fitting to be Q s =13 900 and Q w =13 100, respectively.

Numerical simulation
The numerical simulations of the dynamic equations (5)- (7) were performed using the realistic parameters from table 1 and the stochastic noise term of equation (3). Here, major fit parameters are the gain compression factors  ij (crucial for shaping the input-output curve), the pump efficiency η and the value for the parasitic currents J P . The adjustment of the injection current is given due to the fact that the device is placed on one sample together with 60 other microlasers that are all pumped simultaneously. Hence, the injection current into the measured device obeys some uncertainty. Also the losses and the dipole transition moment were adjusted, but within very limited borders, as both modes needed to stay almost identical.
Starting with a characterization of the steady states, figure 4(a) shows the change of the laser output intensity of both modes versus the injected pump current J (black and orange lines). The experimentally measured intensities are plotted in the same subplot, but indicated with symbols. The two polarization modes are distinguished with a color code of black and orange for strong and weak mode, respectively. The theoretical calculations of the input-output curves depicted in figure 4(a) show an excellent agreement with the experimental data. Basic characteristics of the micropillar laser are reproduced: for low injection currents, strong and weak mode show almost equal intensities. Above the threshold of approximately m 100 A, they show a separation with a typical s-shape for the strong mode and a maximum for the weak mode with continuously decreasing lasing intensity at higher injection currents. In the numerics, the shaping of this intensity profile is mainly achieved due to the phenomenologically introduced gain compression factors  ij . These factors operate in combination with the actual mode intensity and decrease the gain not only of the corresponding mode itself, but also of the competing mode. A high intensity of the strong mode leads to a lowered gain for the weak mode  and less photon emission, consequently. Our b = 0.0056 is chosen to fit the observed s-shape. The somewhat moderate β-factor is explained by the rather large mode volume of the present devices with a diameter of 3.6 μm.
Significantly larger values have been reported for micropillar lasers with adiabatic cavity design and diameters below 2 μm [36]. Thus, we focus on the mechanism of the experimentally observed super-thermal photon bunching by keeping the model as simple as possible to allow for an in-depth bifurcation analysis. 140, 180 A. This super-thermal bunching and its interpretation is the central aspect of this paper. It is accompanied by a peak observed for the measured correlation times (symbols in figure 4(c))that is slightly shifted towards higher injection currents.
Interestingly, this qualitative behavior of the auto-correlation function was also reported in Dye ring-lasers for asymmetric modes [13,37] and seems to be a basic effect of mode-interaction regardless of the systems dimensionality. For the ring laser, switching between the modes was observed and correlation times on a timescale of several microseconds were measured [14] (also interpreted as first passage times [38]). The calculations of the temporal correlation times exhibited by our micropillar laser are much faster (on the order of nanoseconds) and depicted in figure 4(c). The faster dynamics is due to the much faster internal processes. Note that no additional peaks are observed for ( ) ( ) t g ii 2 and thus a beating between the modes can be precluded. The numeric results for the correlation time t corr agree well with the experimental data, as can be seen by comparing the orange line with the symbols in figure 4(c). Again, the increasing correlation time with output intensity, i.e. with pump, is already evident in much simpler models [14,39] which supports our assumption that we can capture the emission properties of the micropillar laser with our rate-equation approach.
Outside of the super-thermal bunching regime there is a discrepancy between experiment and theory found in the auto-correlation function of the weak mode ( ) ( ) g 0 ww 2 (gray area in figure 4(b)). We explain that by the very fast decay of ( ) ( ) t g 2 , i.e. low correlation times of less than 1 ns that are in conflict with the temporal resolution (approximately 1 ns) of the experimental method. There the auto-correlation function cannot be measured accurately, hence we find an underestimation of the measured To understand the underlying nonlinear mechanism responsible for the super-thermal bunching, we calculate the solutions of our model using the deterministic expression of the spontaneous emission term, as described in equation (4). For each mode we find three solutions where two of them are stable. Figure 5(a) shows the possible lasing intensities of weak and strong mode in orange and black, respectively, as a function of the β factor for an injection current of m 160 A (thus within the region of pump currents where large ( ) g ww 2 values are observed). Stable solutions are indicated with solid, unstable solutions with dashed lines. The three solutions that exist for low β factors collide in a saddle-node (SN) bifurcation, and for higher β no bistability exists. The position of this SN bifurcation in parameter space was followed in parameter space (injection current J-and βfactor) by means of path continuation using the tool AUTO07p [40] and depicted as a line in figure 5(b). The hatched areas distinguish between regions of one (red) and three (blue) solutions. Thus, the blue hatched area indicates the parameter region of bistability in our bimodal mode laser system. Following the depicted arrow in figure 5(b) to the right corresponds to increasing the injection current as done for the investigated microlaser in figure 4. Along this line the laser enters the region of bistability which we will show is the reason for the experimentally observed high ( ) g ww 2 values. Taking again the stochastic noise of the spontaneous emission into account (using equation (3)), switching between both stable solutions (e.g., high lasing intensity of strong mode with low intensity of the weak mode and vice versa) is observed. To visualize these switching events, which turn out to be crucial for the observed superthermal bunching, we show numerically calculated time traces in figure 6.
We chose two different ways to visualize the time traces of the two modes. The left-hand side of each plot in figure 6 shows intensity-intensity histograms of the strong (x-axis) and the weak mode (y-axis). The color code indicates how often a state with a certain weak and strong mode lasing intensity is reached, i.e., orange/red indicates high and yellow/white indicates low probability. The right-hand side of each diagram in figure 6 is a snapshot of 15 ns, where black and orange lines show the time evolution of the lasing intensity of the strong and the weak mode, respectively. For low injection currents slightly higher than the threshold ( figure 6(a)) a frequent switching of the lasing modes is observable. Further, the lasing state is anti-correlated, meaning that both modes switch simultaneously. The frequent switching is always accompanied by transients, such that a transversely spread line shows up in the histogram (figure 6(a) left side). Comparing the results to Monte-Carlo simulations of a two-mode ring laser a striking similarity can be found [39].
For higher injection currents, the switching events are less often as seen in the right-hand side picture of figure 6(b). This can be explained with a lower relative noise strength, when the noise from the spontaneous emission term is compared to the distance of both possible stable solutions for the modes. This distance increases with the injection current (as evident from the input-output curves), and the almost constant noise strength can drive the switching events less frequently (see also discussion of figure 8 later on). Hence, the influence of transient states to a whole time series is less and the histogram reveals both emission states of the system with a much more pronounced spot for the strong mode (see figure 6(b) left). The injection current chosen for figure 6(b) is in the regime where we find super-thermal bunching. Hence, the experimentally observed bunching can be regarded as an effect, where most of the energy of one mode (with respect to a certain considered) is conserved in a small number of pulses. This way the mean lasing intensity is low, but its variance is comparatively high. The fraction of both is exactly the classical representation of the auto-correlation function explained in equation (1) and explains intuitively the dynamical mechanism for the photon bunching found here.
A transition back to thermal emission of the weak mode is depicted in figure 6(c), where the system is far above the threshold. For this very high injection currents, the spontaneous emission is not strong enough to kick the system from one lasing state to another. This means that the most stable (strong mode) solution survives and the other mode is almost solely driven by spontaneous emission. This also explains the big variance in the numerical results of the correlation times discussed before in figure 4(c). For long time series (calculation times of 1 ms and longer), very rare but strong pulses for the weak mode can occur, accompanied with a long dwell time due to the low relative noise strength. As these events only take place in one of several realizations, the variance of the auto-correlation function is increased significantly, leading to high values of ( ) . A streak camera image obtained for the emission of our microlaser is depicted in figure 7. It shows rare switching events between the lasing modes, which is similar to what is found for the time traces described in figure 6(b) and thus supports our interpretation. The small dots at the polarization angle of the weak mode show the usual switching events. They are some few pixels long and correspond to sub ms switching that has been demonstrated by simulations both in the time series diagrams and by correlation times on the order of »100 ns in figure 4(c).
To further underline the described mechanism, we complete our numerical results obtained for the microlaser with a systematic parameter study. In figure 8 we present stochastic results for the complete parameter space already investigated in figure 5. The numerically obtained photon auto-correlation function ( ) is color coded in figure 8(a) and plotted as a function of injection current J and β factor. Note that the  (a switching event is defined as a pulse of the weak mode that is stronger than the intensity of the strong mode.) earlier discussed figure 4(b) is a horizontal cross section of figure 8(a) but deviations in the absolute values can occur due to only one simulation being performed per pixel. The super-thermal bunching (red area) shifts for higher β to higher injection currents, which is not surprisingly, because the noise strength increases with increasing β. The black solid line indicates the saddle-node bifurcation, derived by path continuation (also depicted in figure 5). Thus we can state that strong bunching occurs within the bistability regime (right of the black line in figure 8(a)) close to the saddle-node bifurcation line.
Of course, due to both the stochastical character of the system and the transient behavior, we also find ( ) > g 2 ww 2 for injection currents below that bifurcation line. Even if the deterministic analysis shows only one solution, the system is still weakly attracted to the ghost of the second solution. Thus we still find excursions in the phase space which consequently result in photon bunching. The same argumentation applies for superthermal bunching at low injection currents in figure 4(b).
Looking at the distance between the two stable solutions (depicted in figure 8(b)) we observe an increase with J and therewith a decrease of the relative noise strength with J. The spontaneous emission manages less frequently to force the weak mode towards the strong intensity solution. This is also obvious in figure 8(c), where we evaluate the average switching frequency between the modes. The blue colored area indicates low switching frequencies and coincides with the high ( ) ( ) g 0 ww 2 areas in figure 8(a), in accordance with our interpretation of rare switching events causing the bunching. The bigger the distance between the two emission states is, the less often the system switches from one state to another, until no switching occurs and only noise is left (bottom right in the picture figure 8(c)).

Discussion
The super-thermal photon bunching observed in microlasers with QD active material is explained within a semi-classical framework focussing on the nonlinear dynamics of multi-mode interaction. We derived a model that is not only able to yield results in qualitative agreement with the experimental data but is also able to relate the observed effects to the underlying bifurcation structure, thus extending the vivid discussions on microlasers to the field of nonlinear dynamics and therewith also connecting to the Dye ring-laser community. Our bifurcation analysis revealed an underlying bistability which, combined with the cQED-enhanced spontaneous emission noise, evokes random mode-switching. Time series diagrams showed that these switching events lead to super-thermal photon bunching as soon as one mode has a very low mean intensity but can be excited by spontaneous emission to yield fast transients to an intense but short-lived light output. Our results provide a bridge between semiconductor microlasers and the dynamics of the competitive modes in lasers in general. Thus we hope to foster new interdisciplinary research at the crossroads of important and so far independent communities dealing with nonlinear laser dynamics, nanophotonics and microscopic modeling of QD devices. As such they will open up new avenues for the understanding and future applications of cQED enhanced microsystems.