Dynamic signatures of mode competition in optically injected high-β lasers

In this paper we theoretically investigate the statistical light-emission properties of an optically injected bimodal quantum-dot micropillar laser with high spontaneous emission rates. The nanostructured device is described in terms of a stochastic, semiclassically derived rate equation model. We focus on the stochastic switching dynamics between the two fundamental modes and correlate the results with an in-depth bifurcation analysis of the underlying deterministic dynamics. By analyzing different statistical measures, e.g. average intensity, auto- and cross-correlation functions, as well as dwell-time distributions, we give a road map on how to unravel the different dynamic regimes in the presence of large noise from experimentally accessible quantities.


Introduction
Micropillar lasers with an active medium formed by quantum dots are promising candidates for nanophotonic applications. They can be implemented on a semiconductor chip and due to their small footprint can be integrated in high densities [1]. The unique properties of micropillar lasers could enable numerous innovative applications, e.g., in quantum optics [2], reservoir computing [3,4], generation of broadband chaos [5] or optical switches [6]. The emission of micropillar lasers is generally dominated by strong spontaneous emission due to the low mode volume. Additionally, the geometry of the laser cavity leads to the existence of two polarization modes, giving rise to stochastic mode switching [7].
The emission characteristics of solitary micropillar devices have been characterized experimentally [8,9] as well as theoretically with both fully quantum mechanical models [10,11] and with semi-classical approaches [7,12,13]. The emission states can also be controlled via external optical control schemes, e.g., optical injection [6,13], optical coupling of two lasers [14] or optical self feedback [15,16]. Due to the low optical power and the high relative amount of spontaneous emission noise in the small cavity lasers, their emission properties are usually characterized in terms of statistical measures. Such measures commonly include correlation functions of the photon emission, such as the second-order photon auto-correlation, which is a common tool in characterizing the emission photon statistics. A direct measurement of the dynamical laser state as in macroscopic lasers is usually not possible. In this paper we want to deepen the analysis of the noisy emission and clarify the relation between measurements of the stochastic micropillar output and the underlying deterministic bifurcation structure.
In particular, we focus on a bimodal micropillar laser subjected to optical injection. Two-mode lasers under optical injection have been the topic of previous research, e.g., for quantum well devices with two longitudinal modes [17,18], for VCSEL [19][20][21][22][23][24][25][26] with low spontaneous emission, or for two-mode QD lasers [27][28][29] where the two modes do not share the same gain medium. Generally, optical injection can lead to a multitude of different dynamics in the emitted light of a single mode laser, including stable locking, periodic oscillations, and complex or chaotic emission [30]. Depending on the damping of the relaxation oscillations of the optically injected laser, the regions of complex dynamics can be more or less pronounced [31][32][33]. In a micropillar laser the dynamics under optical injection are more involved, due to the existence of the two modes and the high amount of noise, and a detailed understanding is still missing. In this paper, we focus on the interrelation between the stochastic properties and the underlying deterministic dynamics. We investigate the deterministic bifurcation structure of the bimodal micropillar laser under optical injection and highlight the role of the mode competition in the context of multistability between the emission in the two laser modes. Subsequently we perform a deep stochastic analysis of the micropillar laser in terms of the photon auto-correlation in the individual modes as well as the cross-correlation between the modes. The combination of these two techniques allows us to predict the dependence of these experimentally accessible statistical quantities on the underlying dynamical structure.
The paper is organized as follows. In section 2 we present the numerical model used for the investigation of the bimodal micropillar laser in terms of its deterministic bifurcation structure and its stochastic properties. The deterministic and stochastic dynamics under optical injection is presented in section 3, focusing on the correlation properties of the emitted photons. We present results on the mode-switching dynamics in section 4, which reveals further insight into the underlying dynamics. A summary is given in section 5.

Theoretical model
The model used in this paper is based on semiclassical rate equations [34] taking into account the electron scattering mechanisms into the quantum dots as derived in our previous works [7,35,36]. Due to structural asymmetries in micropillar lasers the fundamental cavity mode splits into two orthogonal linearly polarized modes [37]. In our model this is accounted for by two separate complex electric field equations, denoted by E j , with j ∈ {s, w} for the strong and weak mode, respectively. Weak and strong corresponds to the mode with lower and higher lasing intensity beyond threshold. We model the electrical fields E j of the two modes, the occupation probability of the active and inactive quantum dots ρ (in)act , and the reservoir carrier density n r . As the QD transitions are inhomogeneously broadened, only optical transitions which are sufficiently close to the lasing mode in frequency contribute to the stimulated emission. Thus, parts of the QDs present in the device are able to contribute to the lasing process and are therefore called active, whereas the remaining QDs at the flanks of the inhomogeneous distribution are called inactive and increase the optical losses inside the cavity [38,39].
We implement the following model equations: The laser is pumped by applying an electric current J to the reservoir n r from where electrons may either recombine without contributing to the lasing mode, denoted by the recombination rate τ r , or scatter into quantum dots with the rate S in . We account for experimental details in the pumping process by assuming an injection efficiency η, and a parasitic current J p , determined by fitting the experimental input-output curves [6]. A and V are the effective mode area and volume, respectively, e 0 is the elementary charge, and Z QD is the number of quantum dots which actively contribute to lasing. The spontaneous recombination lifetime in the quantum dots is given by τ sp . The central lasing wavelength is given by ω, κ j denotes the cavity loss rate of the respective laser modes, and bg = n 2 bg is the background permittivity. The linewidth enhancement factor is given by α.
The occupation probability of the inactive quantum dots is calculated from the steady-state value of an equation equivalent to equation (2) but without the stimulated emission term and reads: The electric fields of weak and strong mode both interact with the active QDs by stimulated emission. Since the frequencies of the two modes differ by only a few tens of μeV, we consider only one carrier population that is interacting with both optical modes. Thus, both modes compete for the same gain medium and the gain coefficient g j in equations (1) and (2) is given by where c 0 is the vacuum speed of light. The two laser modes are considered to only interact via the common gain medium, without a direct interaction of the optical fields. Our results therefore hold for a large range of frequency spacing between the modes, provided a common gain medium exists. The gain compression coefficients ε ij with i, j ∈ {s, w} have been chosen to best fit the experimental results, however, the coefficients can also be derived from the spin-flip model (or San-Miguel-Feng-Moloney model) [40], which is commonly used to describe polarization dynamics in vertical-cavity surface emitting lasers. The derivation of our model equations from the spin-flip model is detailed in appendix A. Spontaneous emission into the lasing modes is modeled via a complex Gaussian white noise source ξ j (t) ∈ C, where ξ j (t) j = 0 and ξ j (t)ξ j (t + t ) = δ(t )δ j,j , where · denotes the time average, such that The factor β describes the fraction of spontaneously emitted photons coupled into the laser modes. Without the contribution of the spontaneous emission, i.e. for β = 0, we have a deterministic differential equation system, which we will also analyze in terms of its bifurcation structure in section 3. We do not consider charge-carrier noise in the model equations as the spontaneous emission noise in the optical field dominates the stochastic properties of the laser. The field equation (1) is formulated within the rotating frame of the cold-cavity frequency ω = 2πν. The free-running frequency of each micropillar mode, ν (0) j , deviates from ν due to the amplitude-phase coupling. We define this shift as δν (0) j . When we consider optical injection into the strong mode of the micropillar, the detuning Δν inj = ν inj − ν (0) s measures the distance between the injecting master laser frequency ν inj and the free-running strong mode frequency ν (0) s . Within the rotating frame of ω, the optical injection term in equation (1) for strong mode injection reads: The injection strength is normalized to E 0 = |E s | 2 + |E w | 2 , the averaged free-running amplitude of the electric field inside the micropillar laser device. The parameter K is related to the optical power of the free-running laser P (0) and the injected power P inj via K = P inj /P (0) . The weak mode is assumed to be perfectly isolated from the injected signal, and the injection term is only added in the strong mode equation. The parameters used in the simulations are listed in table D1.

Statistical evaluation of the lasing dynamics
The solitary bimodal micropillar laser without optical injection (K = 0) emits light at both cavity modes, which are typically separated by only a few tens of μeV in energy. A sketch of the laser is shown in figure 1(a). A typical input-output curve is shown in figure 1(b). Here, the averaged output intensity of the two modes I j = 2κ j ε 0 ε bg V |E j (t)| 2 is plotted as a function of the pump current J. For injection currents J below threshold (J th ≈ 17 μA), the averaged mode intensities are mainly determined by spontaneous emission and coincide. A statistical representation of the stochastic emission can be seen in figures 1(c) and (d), where the histograms of the intensity distribution are projected onto the (|E s | 2 , |E w | 2 )-plane (darker colors correspond to higher probability).
For the following investigation of the dynamics under optical injection (strong mode injection), we choose a pump current of J = 30 μA (thick vertical dashed line in figure 1). At this operation point the average intensities of the two modes are already clearly separated. Nevertheless, even well above threshold stochastic switching between weak and strong modes occurs [compare figures 1(c) and (d)]. A switching event denotes a momentary, usually very short, exchange of intensities between weak and strong mode. For the solitary laser, it is due to the coexistence of two stable solutions [7].  To characterize the dynamics of the optically injected bimodal emitting high-β laser, usually the mean intensities of both weak and strong mode as a function of the injection parameters (injection strength K and frequency detuning Δν) are evaluated. This has been done experimentally and numerically for this kind of lasers in [6,13]. Figure 2(a) shows the ratio between the intensity of both modes color-coded as a function of the injection parameters. Within the locking region (light colors) the strong mode is locked to the injected signal and the weak mode is suppressed, i.e. the spontaneous switching is much less likely to occur. Outside the locking cone we observe regions where the strong mode is suppressed in favor of the weak mode (dark blue colors).
In order to explain the occurrence of the different dynamic regimes, we investigate the underlying deterministic dynamics, i.e., the bifurcation structure of the injected two-mode laser without noise (β = 0). The bifurcation lines define the boundaries between different solutions in the injection-parameter plane. We obtain the deterministic bifurcation lines via numerical path continuation of equations (1)-(3) using the AUTO-07p package [41]. In the deterministic system, for β = 0, the dynamics are characterized by the coexistence of single-mode solutions, in which only the strong mode is lasing, and two-mode solutions with simultaneous strong and weak-mode emission. The bifurcation lines are shown in figure 2(a), with the color code showing the intensity ratio of the strong and weak mode, obtained from stochastic simulations. The orange lines mark transcritical bifurcations (TC) where two-mode solutions intersect single-mode solutions (TCLC mark TC bifurcations of limit cycles). Within these regions the strong mode is suppressed in the stochastic simulations [dark blue color in figure 2(a)], indicating that the weak mode has turned on. In the stochastic simulations these transitions are smoothed out, but nevertheless can be attributed to the nearby deterministic bifurcations. In the regions where only the deterministic strong-mode solution is stable and the weak mode is not lasing, a small weak mode intensity is found due to the spontaneous emission. The locking region [white in figure 2(a)] is bounded by an SNIPER bifurcation (saddle-node-infinite-period bifurcation) for negative detuning and by a TC bifurcation line at positive detuning, beyond which the weak mode slowly turns on. It is noted that only the physically relevant bifurcations, i.e. the ones that occur on stable and thus observable solutions are plotted in order to aid clarity of the bifurcation diagram.
A more detailed picture of the solutions and their bifurcations is shown in figure 2(b) along a vertical cut through the bifurcation diagram at K = 0.75. It displays the average intensity of the deterministically found stable and unstable solutions of the strong mode emission. Black lines indicate single-mode solutions where the weak mode is off, while red lines indicate solutions with two-mode emission (we show only the strong mode intensity for clarity). The nonlinear oscillations outside the locking region (limit cycle dynamics, LC) are indicated with thick lines (the lines indicate the mean intensity of the oscillations, the time dependence for the single mode oscillations can be seen in figure 5(a)). Dashed lines represent unstable solutions. Starting from negative detuning [upward along the vertical dashed line in figure 2(a)] two solutions coexist; a single-mode limit cycle and a stable two-mode emission which is stable below the Hopf bifurcation (green cross in figure 2(b)). At the Hopf bifurcation, a two mode LC is born (thick red line) which loops back and changes stability at the saddle node of limit cycles (SNLC) around -10GHz. At the transcritical bifurcation of limit cycles (TCLC), the strong-mode limit cycle intersects the two-mode LC solution under a change of stability. Thus, in between the transcritical bifurcation curves [within the dark blue region in figure 2(a)], only the stable two-mode fixed-point solution remains. In the region between the SNLC and SNIPER bifurcations (open circles) only the oscillating single-mode solution is stable, which disappears in the SNIPER bifurcation, in which the LC connects to the phase-locked solution. In the single-mode phase-locked solution the weak mode is suppressed until the TC bifurcation is reached at positive detuning (orange square at 5GHz), in which again a stable two-mode fixed-point solution becomes stable. From the TCLC bifurcation at 15 GHz on, again a bistability between a single-mode oscillation and a two-mode steady-state emission is found, which exists for all positive values of detuning beyond this point.
Comparing the results from the bifurcation analysis (the solution structure) with the information gained from the mean intensities of the micropillar emission, it becomes obvious that an experimental characterization via the mean intensities does not reveal the full dynamics. Especially the identification of oscillating dynamics or multistability is not possible. Numerically, switching between strong and weak mode emission, can be detected directly from the intensity time series. Experimentally, the fast time scale of the switching and the low intensities in micropillar devices requires other statistical measures. We will show that information about the underlying stability structure can be inferred from measurements of the photon statistics. A well established and experimentally accessible statistical method is the measurement of autoand cross-correlation functions g (2) ij (τ ) [42,43], where i, j indicate either strong or weak mode. These quantities characterize the emitted light of the laser and can be measured, e.g. by a Hanbury-Brown and Twiss setup [44]. For our purposes, we use the classical representation of this quantum-mechanically defined statistical measure: This allows us to calculate the correlation functions from the time series. Usually, for single mode emission, g (2) ii (τ = 0) is of interest, because it allows to differentiate among thermal, coherent and (when using the quantum-mechanical definition) non-classical light. For thermal, i.e., spontaneous, emission the auto-correlation function takes the value g (2) ii (0) = 2 while g (2) ii (0) = 1 is obtained when the laser emits coherent light. Values of g (2) ii (0) > 2 are an indicator for super-radiance [45] or super-thermal bunching [7]. Since we are also interested in the correlation between strong and weak mode, we need to determine the cross-correlation function g (2) ij (τ ) [6]. Here, anti-correlation of the two competing modes is the strongest, where the cross-correlation function g (2) ij (0) reaches its lowest value. A value of g (2) ij (0) = 1 indicates uncorrelated emission between the strong and weak mode.
To illustrate the signature of deterministic features in the stochastic dynamics, we again look at a line scan at K = 0.75 as a function of the detuning and plot the mean intensities I j of the emission for the deterministic simulation and the stochastic simulation in figures 3(a) and (b), respectively. Additionally, the auto-and cross-correlation of the two modes, g (2) ww (0), g (2) ss (0), and g (2) sw (0), are plotted in figure 3(c). We obtain the deterministic results from direct numerical integration of the model equations using a fourth-order Runge-Kutta integration scheme with a time step of 1 ps. The stochastic time series is produced using a stochastic Euler integration scheme with a time step of 0.1 ps. Going from negative to positive detunings Δν inj along the line scan (compare figure 2(b)) one comes across four distinct regions  figure 3(b), the change between region I and II shows up as an exchange of the role of strong and weak mode, however, not at the bifurcation (vertical dashes lines in figure (3)) but well beyond at about -27GHz. The auto-correlation in region I, figure 3(c), shows a similar behavior and low stochastic mean intensities [ figure 3(b)] are accompanied by high g (2) ii (0) values. The locking region IV is clearly identified by g (2) ww (0) = 2 and g (2) ss (0) = 1 while we can identify the two-mode region II by g (2) ww (0) = 1. In region II, a change in detuning induces a redistribution of the lasing power between the strong and weak modes, while the total output power remains constant [46,47]. The bifurcation scenario observed here for the micropillar laser is similar to that observed in conventional VCSELs with parallel optical injection [25] or comparable two-mode lasers [48] which exhibit comparable polarization switching dynamics. As the strong-mode power is reduced in region II, its emission becomes dominated by spontaneous emission, smoothly increasing g (2) ss (0) toward 2 (see appendix B). We note that  (2) ss (0), g (2) ww (0) and (c) cross-correlation functions g (2) sw (0) color-coded within the (K, Δν inj )-parameter space. The cross-correlation decay time τ sw shown in (d) is obtained as full-width-at-half-minimum of g (2) sw (τ ) (cf figure 6). The bifurcation lines of the deterministic system without noise are shown as in figure 2(a). The dynamical regions I-IV are labeled as in figure 3. in region II there is destructive interference between the injected signal and the laser mode, reducing the strong-mode power. This region extends beyond the locked region obtained for a single-mode laser for which the weak mode does not lase. The laser remains phase-locked only due to the presence of the weak mode.
The bifurcations are not immediately visible in the stochastic data due to the smoothing of dynamical transitions. Nevertheless, their statistical signature is visible as dips in the cross-correlation function g (2) sw (0) [black line in figure 3(c)], that show an increased anti-correlation between the strong and weak mode close to the bifurcations. Especially the SNLC, SNIPER and TC bifurcations show up as narrow dips in g (2) sw (0) [figure 3(c), black line]. The two pronounced regions of anti-correlation at large absolute detunings occur within the region of bistability, region I. They can be understood qualitatively by the evolution of the stability of the solutions beyond the outermost TCLC bifurcation. In figure 3(d) we plot the real part of the largest eigenvalue λ of the solutions. For the limit-cycle solution we evaluate λ = (log M)/T from the largest non-trivial Floquet multiplier M and the period T. While the strong-mode limit-cycle solution is least stable close to the outermost TCLC bifurcations, its stability increases (Reλ becomes more negative) when moving further away from zero detuning. Simultaneously, the stability of the dual-mode fixed-point solution decreases as Reλ approaches zero. Consequently, the strong-mode solution becomes dominant far away from the locked region, leading to the near-complete suppression of the weak mode as seen in figure 3(b). In the transition region between the outermost TCLC bifurcation and the edges of the shown detuning range, frequent and long-lived switching events take place in the regions where both bistable solutions are still sufficiently stable. In the stochastic simulations, the bistability born in the TCLC bifurcation thus becomes apparent only in a limited region of pronounced anti-correlation. Across the whole detuning range, the cross-correlation function is closest to unity (uncorrelated strong and weak modes) in regions where no bistability exists (II, III, IV), as the switching events can be expected to be least frequent in these regions.
In figure 4 we plot numerical results for the auto-correlation (a), (b) and cross-correlation (c) g (2) ij (0) within the full (K, Δν inj )-plane. Additionally, we show the cross-correlation time τ sw in figure 4(d). It is calculated from the full-width-at-half-minimum of g (2) sw (τ ) [see figure 6(b) for selected plots of g (2) sw (τ )]. From the auto-correlation data the locking region IV with stable strong mode lasing and thermal weak mode characteristics is easily identified by the auto-correlation values g (2) ss (0) = 1 and g (2) ww (0) = 2 (region in between SNIPER and TC). The stable two-mode lasing regimes II are found in between the transcritical bifurcation lines, characterized by g (2) ss (0) ≈ 1, g (2) ww (0) = 1. The very high value of g (2) ww (0) at the edge of the locking region IV at negative detuning is an indication for strong but rare spikes, which is commonly found close to an SNIPER bifurcation [49]. For the strong mode, regions with g (2) ss (0) > 2 are found within the bistability region I close to the border of the two-mode regions II where the weak mode dominates. The reason is the continuously increasing stability of the strong mode emission in I which eventually leads to rare on-switching events of the strong mode and thus to g (2) ss (0) > 2. Considering the cross-correlation g (2) sw (0) in figure 4(c) of the stochastic dynamics under optical injection, we observe a similar structure in the (K, Δν inj )-plane, however this time we gain different information. As expected, within region IV where the strong mode is locked, there is no correlation between the lasing strong and the thermal weak mode (g (2) sw (0) ≈ 1; yellow colors). Interestingly, the same holds for region II with stable two-mode lasing (the regions enclosed within the TC and TCLC curves). Thus, when the laser operates in the two-mode lasing regime II, there is only weak anti-correlation between the strong and weak mode emission.
Regions with the strongest anti-correlation [blue colors in figure 4(c)] are located within the bistability region I, where the two-mode solution and the strong-mode limit cycle oscillations exist and are equally stable. Their position is found well beyond the TCLC line, and next to the region where the strong mode shows g (2) ss ( . Values up to ≈ 6.5 ns are found which suggest a strongly enhanced and long-lived anti-correlation between the modes. We want to clarify this further by analyzing the switching dynamics in the next section.

Dwell-time statistics of mode switching events
In the following, we will investigate the interplay of the mode-switching dynamics and the underlying deterministic dynamics. So far we have seen that depending on the injection parameters, the auto-and cross-correlation functions of the emission show distinct features that allow us to characterize the emission. From the numerically obtained time series, we can additionally extract the dwell-time distribution of the occurring switching events. Within the stochastic time series of the two-mode emission, the time between consecutive switching events can yield additional information about the underlying dynamical system. Both the stability and the distance in phase space of coexisting solutions will influence the switching frequency and the dwell-time statistics. Or, conversely, the measurement of the switching statistics along with correlation functions will yield information about the underlying deterministic bifurcation structure.
The dwell time t dw of a given switching event is defined by the time interval over which the intensity of one mode becomes dominant over the other. We mark the beginning of a switching event when the mode intensity exceeds the other mode by one third of its value, I i > 4 3 I j , and the switching event ends when the inverse relation is fulfilled, I j > 4 3 I i . An example of such switching events and the corresponding dwell time is sketched in figure 5(c) for a section of the simulated stochastic time series. Here, at first a switching event of the weak mode is seen, with a dwell time t weak dw , followed by a strong-mode switching event with a dwell time t strong dw . The time series thus consists of a series of alternating switching events of the two modes. Depending on the injection parameters, the switching dynamics can change, with varying statistics of the switching frequency and the dwell times. For comparison, we plot the two stable deterministic solutions that coexist for parameters in region I without noise and thus without switching events in figures 5(a) and (b), (Δν inj = 10 GHz, K = 0.2). The intensities perform either stable beating oscillations in the strong mode close to the detuning frequency single mode limit cycle oscillations shown in ( figure 5(a)) or stable two-mode emission occurs [ figure 5(b)]. The beating frequency of the LC oscillation is still observed in the noisy time series in figure 5(c). It is, however, interrupted by the switching events. We note that the two-mode solution, while being a fixed-point solution, possesses a pair of eigenvalues with imaginary parts close to the beating frequency found in figure 5(a). In the stochastic time series, noise will therefore excite oscillations at that frequency during the switching events also when the weak mode is dominant.
In the following, we evaluate the length of individual switching events, yielding a statistic distribution of dwell times t dw . To calculate the dwell-time distribution, we record all dwell times during a 10 ms long time series. The results for Δν inj = ±10 GHz (within region I at K = 0.2) are plotted as black and gray lines in figures 6(a) and (b) for weak-mode and strong-mode dwell times, respectively. They show a multi-peaked behavior, i.e. different intervals between switching events occur with higher probability while others are improbable. We find minima in the probability distribution at integer multiples of the beating frequency, i.e., close to the detuning |Δν inj |. To understand those multi-peak distributions, one has to consider the limit cycle solution of the underlying beating oscillations in region I. The strong and weak mode oscillate  , for different detuning Δν inj generated from time series of 10 ms at K = 0.2 (solid lines) and for K = 0 (dashed line, multiplied by 10 for better visibility). Black and gray lines show multi-peaked distributions obtained in region I where switching is induced by the underlying beating oscillations, red lines correspond to the dwell-time distributions found within the two-mode region II. (c) Cross-correlation function g (2) sw (τ ) as a function of τ for the detuning values chosen in (a) and (b). The arrows show the correlation time τ sw for Δν inj = 2.1 GHz.
periodically and in anti-phase (see figure 5). At time instances where both intensities are close, it is more likely that a noise perturbation excites a switch than at time instances where they are far apart in phase space. Thus, a higher probability of switching times repeats every full period of the oscillations and the maxima can be described by where the deterministic beating period T osc satisfies At the detuning of Δν inj = ±10 GHz (region I) both the strong-mode and weak-mode dwell-time distributions show a similar histogram structure [black lines in figures 6(a) and (b)]. On the contrary, in region II at Δν inj = 2.1 GHz (red lines), a multi-peaked structure is only seen for the weak-mode dwell-time distribution, whereas the strong mode exhibits a single exponential decay toward longer dwell times (a semi-logarithmic plot of the distribution in shown in appendix C). The detuning of Δν inj = Long-time correlations that are induced by the underlying deterministic dynamics are directly visible in the time-resolved cross-correlation between the two modes. This is shown in figure 6(c) where we plot the corresponding g (2) sw (τ ) for the previously discussed detuning values. In the case of pronounced deterministic oscillations (region I), the cross-correlation function has long exponential tails with a modulation at the oscillation period. A time-resolved measurement of the second-order cross-correlation function can therefore be used to gain further insight into the underlying deterministic dynamics.
The cross-correlation function and the dwell-time statistics can change significantly depending on the injection parameters. We thus show the corresponding quantities along a detuning scan at an injection strength of K = 0.2 in figure 7. Panels (a) and (b) show the cross-correlation depth g (2) sw (0) and cross-correlation time τ sw , respectively. In panel (c) we show the corresponding dwell-time distribution for weak-mode switching events t weak dw . The horizontal axis represents the different times between switches (as in figure 6) while the color code indicates the probability of a switch occurring with dwell time t weak dw . We show the corresponding plots for K = 0.75 in appendix C.
The multi-peaked dwell-time distributions can be seen in the whole bistable region I. They form a curved rippling structure in the (Δν inj , t dw )-plane as the underlying deterministic beating frequency depends on the detuning Δν inj . The white lines in figure 7(c) are calculated from equation (9) and coincide with the peaks in the dwell-time distribution. The ripple structure is most pronounced in the regions of strong anti-correlation. This is accompanied by a significant increase in the correlation time, which can now be understood from the deterministic oscillations introducing long-term correlations in the switching dynamics. Conversely, close to the TC bifurcation in region II, a significant anti-correlation is observed, while the correlation time remains relatively short. This is accompanied by the disappearance of the multipeaked structure in the dwell-time distribution as the imaginary part of the leading eigenvalue vanishes when approaching the TC bifurcation, thus suppressing the long-term correlation between the modes. In contrast, the existence of a bistability greatly enhances the correlation time, as seen in region I [cf figure 4(d)]. In region III, behind the Hopf bifurcation that destabilizes the strong-mode oscillations, remnants of the riple structure can be seen in figure 7(c) due to noise induced excitations of the former stable solution, however with a correlation time close to zero. Within the locking region IV hardly any switches occur [white in figure 7(c)] and the cross-correlation g (2) sw (0) reaches 1, indicating no correlation. The dwell times in optically injected micropillar lasers have not been directly measured and are accessible only via indirect measures such as the cross-correlation time. In such measurements, correlation times of a few hundred picoseconds to a few nanoseconds have been observed, which matches the results found here [6].

Conclusion
We characterized the stochastic emission of a two-mode micropillar laser under external optical injection into its strong mode. The deterministic dynamics were investigated in terms of the bifurcation structure in dependence of the injection strength and relative detuning frequency, revealing the interaction between the two laser modes. Around the phase-locked region we found transitions between single-mode and two-mode solutions, with significant regions of bistability which are delimited mainly by transcritical bifurcations of fixed points and limit cycles. By comparison of the deterministic bifurcation analysis with the stochastic simulations, we revealed the signatures of the underlying deterministic bifurcation structure in stochastic measures. In particular, we found an enhanced anti-correlation between the two laser modes at these bifurcations. This anti-correlation is found to be the result of enhanced stochastic mode-switching, i.e., intermittent emission on the respective weaker laser mode. Experimentally, the cross-correlation can be measured even for the low output power encountered in these micropillar lasers, which could be used to recreate the underlying bifurcation structure from the statistical analysis of high-β lasers.
Furthermore, we evaluated the stochastic dwell-time distribution of individual mode-switching events. A multi-peaked dwell-time distribution is found in regions where the deterministic solutions possess an intrinsic period, i.e., in foci or around periodic orbits. In regions of multistability between single-mode beating oscillations and steady-state two-mode solution this multi-peaked distribution is especially pronounced, and the deterministic dynamics introduce long-time correlations leading to a slowly decaying cross-correlation between the laser modes. Instead, if no multistability exists, the dwell-time distribution exhibits a single peak, as known from standard Kramers escape processes, and the correlation time is low. These observations will aid the characterization of the dynamics of micropillar lasers and other low-power devices by means of experimentally accessible quantities such as the second-order photon cross-correlation.  With these definitions, we can rewrite the field and gain equations in terms of the x and y polarization, corresponding to strong and weak mode: We identify 2γ p as the frequency difference between the two modes, and 2γ a the difference in the cavity decay rate. The two linearly polarized modes are coupled via the population difference δρ. In the limit of fast spin relaxation rate, we can adiabatically eliminate the last dynamic equation by setting d dt δρ = 0, which yields the quasi-static relation   We write E y = |E y |e iθ and set E x = |E x |, such that θ gives the phase difference between the two linearly polarized modes. Inserting this into equation (A.8) yields: For the noise-dominated and anti-correlated dynamics observed in the studied micropillar laser, we can assume θ to be randomly and uniformly distributed. We thus calculate the phase-average of the above equation using e iθ cos θ = 1 2 : The second term can be neglected for |E x (t)| 2 2 2 γ s |μ| 2 T 2 , which, for our parameters, corresponds to output powers below ≈ 10 μW. The remaining term yields the cross-gain compression factor, with an equivalent expression for d dt E y (t).