Mapping the Stability and Dynamics of Optically Injected Dual State Quantum Dot Lasers

: Optical injection is a key nonlinear laser conﬁguration both for applications and fundamental studies. An important ﬁgure for understanding the optically injected laser system is the two parameter stability mapping of the dynamics found by examining the output of the injected laser under different combinations of the injection strength and detuning. We experimentally and theoretically generate this map for an optically injected quantum dot laser, biased to emit from the ﬁrst excited state and optically injected near the ground state. Regions of different dynamical behaviours including phase-locking, excitability, and bursting regimes are identiﬁed. At the negatively detuned locking boundary, ground state dropouts and excited state pulses are observed near a hysteresis cycle for low injection strengths. Higher injection strengths reveal µ s duration square wave trains where the intensities of the ground state and excited state operate in antiphase. A narrow region of extremely slow oscillations with periods of several tens of milliseconds is observed at the positively detuned boundary. Two competing optothermal couplings are introduced and are shown to reproduce the experimental results extremely well. In fact, the dynamics of the system are dominated by these optothermal effects and their interplay is central to reproducing detailed features of the stability map.


Introduction
Optical injection is a key technique in many modern photonic systems and particularly unidirectional optical injection, where light from a primary laser is injected into the cavity of a secondary laser. It is used in high sensitivity signal sensing, to improve coherence of high power lasers, for spectral density enhancement in optical communications, and for arbitrary pulse shaping, among many other applications. Mapping the behaviour of an optically injected laser over control parameter space is thus of central importance for many applications. A stability map where one varies the injection level and the detuning (the frequency of the primary laser minus that of the secondary laser) and generates a diagram characterising the output of the injected device is one of the most effective tools in this regard [1]. Stability maps can take several forms and can be as simple as displaying the region where phase locking is obtained. Indeed, for many applications this is the most important regime, and this simple map is all that is needed. On the other hand, optically injected semiconductor lasers are also rich sources of fundamental non-linear dynamics, and the stability map can also be used to classify regions of different behaviour including, but not limited to, excitability, multistability, oscillatory behaviour, and chaos [1][2][3]. Such regions can be of interest both for fundamental studies and for applications. For example, device can be made to emit from the GS while free-running. Several dynamical regimes have been reported and analysed in this configuration, such as fast state switching [22], all-optical gating [23], dual state excitability [24], and neuromorphic bursting [15], with excellent agreement between experiment and theory.
We present a comprehensive combined experimental and theoretical analysis of the stability map for an optically injected dual state QD device. As well as mapping out previously discovered dynamical regions, we uncover new features such as slow oscillations with a period of 10s of milliseconds near the positively unlocked boundary. Optothermal effects are shown to play an important role in the system. Superb agreement is obtained between experiment and theory. We also note a marked absence of chaos in the system.

Experiment
The device under investigation is a 300 µm QD laser composed of InAs quantum dots on a GaAs. It has the same epitaxial structure as the device used in [14,21] but is significantly shorter (300 µm here compared to 600-900 µm in [14,21]). As a result, the device under test here never lases from the GS when free-running, and instead, lases only from the ES. It is pumped at 75 mA (1.3 times threshold at 20.5 • C). A schematic of the experimental setup is shown in Figure 1. The device is mounted and placed on an xyz stage. The primary laser (PL) is a commercial tunable laser source with minimum step size of 0.1 pm (∼0.0178 GHz). Light from the PL is injected into the secondary laser (SL) -the QD laser -via an optical circulator with an isolation greater than 30 dB. A polarisation controller is used to set the polarisation of the injected light and maximise coupling. The light from the SL enters the second port of the circulator and is directed to a filter where the ES and GS are separated. The ES emission is sent directly to a 12 GHz detector connected to a high speed, real-time oscilloscope. 10% of the GS power is sent to a power meter (PM) to monitor alignment and 90% goes directly to another 12 GHz detector, again connected to the oscilloscope. Figure 1. Schematic experimental setup. SL is the secondary (QD) laser. PL is the tunable primary laser. Light from the PL is is sent to a circulator and is then injected into the SL. A polarisation controller (PC) is used to maximise coupling. Light from the SL is sent to the circulator and then to a filter where the ES and GS light are separated. The ES light is sent directly to a 12 GHz detector and the GS light goes to a 90/10 splitter. The 10% is sent to a power meter (PM) for alignment control, and the 90% goes to another 12 GHz detector. Both detectors are connected to a high speed real-time oscilloscope. Red lines are GS light; blue lines are ES light; purple lines contain both GS and ES light; and black lines are high speed electrical cables.
Typically, one can define the experimental injection level as the ratio of the magnitude of the electric field of the injected light to the magnitude of the electric field of the free running laser. However, in this work, when the device is free-running, the GS is never lasing and so this measurement cannot be used. Instead, we define injection strength as the square root of the power of the PL reaching the facet of the SL normalised to the maximum power of the PL reaching the facet of the SL in the experiment, K = P facet P facet max . Thus, the maximum value of K is 1. Similarly, as the GS mode being injected is always subthreshold, the detuning is also difficult to define. We thus pragmatically define zero detuning to be where the ES is at a minimum power after injection.
The experimental technique to build the stability maps is similar to that of [2], but here we fix the detuning and vary the injection strength for each slice of our map. We record and analyse 20 µs-long time series of both the GS and ES intensities at regular intervals. Initially the frequency of the PL is set and the power of injected light reaching the facet of the secondary laser from the PL is swept from 0.25 mW to 2.47 mW in 120 equal steps of K. Then, a down sweep is performed where the power reaching the facet is swept from to 2.47 mW to 0.25 mW with the same 120 steps. The wavelength is then decreased in steps of 0.2 pm (∼0.0356 GHz) and the power sweeps are repeated after each step.

Stability Maps
The average power from the GS is plotted in Figure 2a,d and that of the ES is plotted in Figure 2b,e. The variance of the ES output is also calculated and plotted in Figure 2c,f. The variance plots are particularly useful for immediate identification of the constant output regions (dark blue) and areas where the intensity is oscillating (green/yellow). When the ES intensity is constant, that of the GS is also constant, and when the ES intensity is oscillating, that of the GS is also oscillating. The upper panels in all stability map figures in this section show the results generated when the injection sweep is upwards (from low strengths to high strengths) and the lower panels show the corresponding maps for the downward sweep (from high strengths to low strengths).  In Figure 2a, the bright yellow region in the centre is the phase-locked region. Here, only the GS is lasing and the ES is completely quenched as indicated by the corresponding dark blue region in Figure 2b and the zero variance in Figure 2c. To the negative detuning side of the phase locked region there is a broad area of dynamics as indicated by the green and yellow region in the variance in Figure 2c. On the positive boundary of the phase locked region there is a region of very slowly varying GS power, indicated by the non-smooth colouring/speckled pixels in Figure 2a. We defer discussion of this for now and return to it below.
To the positive side of the slowly varying speckled region, the device behaves as a dual state emitter, with both states emitting constant intensities. As the detuning is increased there is a smooth, continuous evolution of both, with the GS power decreasing and the ES power increasing. The same dual state behaviour is mirrored in the negative detuning region to the left of the dynamic oscillatory regime. We interpret this as follows: in the GS only system (or indeed any conventional optical injection system), there are unlocking boundaries outside of which there are only oscillating unlocked solutions. In the limit of large detuning, the oscillations can be physically interpreted as a beating between the injected light and the emission of the secondary laser. However, in our dual state system there is no free-running GS. Thus, this beating cannot arise. The output of our QD SL can then be thought of as a mixture of regenerated injected light [4,25] and ES emission from the secondary laser. We thus distinguish between the phase-locked output (dark blue) in Figure 2b, where the ES is completely off, and the regenerated injected light output of somewhat large detunings where there is dual state emission. Figure 3 shows just some of traces taken during a single upsweep of the injection strength K for a fixed detuning of −6.9 GHz. The coloured labels correspond to the coloured dots in Figure 4, marking their location on the maps. (The maps in Figure 4 are the same as in Figure 2, with some additional markings to allow for comparison with Figure 3). Figure 3a shows a periodic train of GS dropouts and the corresponding ES pulses as reported in [23,24]. As the injection strength is increased these periodic trains disappear, and in Figure 3b a bursting oscillation is obtained while various mixed mode oscillation (MMO) and bursting MMO regions [26,27] are observed in Figure 3c-e. In Figure 3b there is a switching between a quiescent phase, where the GS is on and the ES is off, and an active phase with oscillations in both states (similar to the dropouts and pulses of Figure 3a). A further increase to K = 0.61 leads to a different bursting dynamic as shown in Figure 3c. In this region, the switching is between a quiescent phase with the GS on and the ES off and an evolving bursting phase. The switch to the bursting phase is via an initial period of decreasing amplitude oscillations, followed by a long series of growing oscillations before a switch back to the quiescent phase. Such evolving bursts were previously reported in [15] and shown to arise via an optothermal coupling. The evolution of the bursting state arises via a determinstic thermal sweep of the detuning leading to slow passages through several bifurcations, breaking underlying bistabilities similar to the generation of excitable square waves in the GS only system in [14]. In Figure 3d K = 0.68, the trace is qualitatively similar to that in Figure 3c, but in Figure 3d the oscillations following the GS dropout die away quickly and the bursting part of the trace is shorter. Finally, when K is increased to 0.99 as shown in Figure 3e, there are no oscillations following the GS dropout and only a very short region of oscillations preceding the return to the quiescent, GS on phase.

Hysteresis and Bistabilty
In order to identify regions of bistability in the maps, we can compare the upsweep and downsweep figures in Figure 4. The upper panel shows the upsweep with a red line superimposed, marking the boundaries from the downsweep. The lower panel shows the downsweep with a white line superimposed, corresponding to the boundaries of oscillatory behaviour from the upsweep. The boundaries are largely the same, and reveal bistabilities where they differ. The negatively detuned side of the phase-locked region (the red/white line in the centre of each figure) is the same in both cases. At high injection strengths and negative detunings, the two directions match at the boundary between the dynamical region and the regeneration region. However, at moderate injection strengths from approximately K = 0.5 to K = 0.65, and between detunings of approximately −8.5 GHz and −6 GHz, there is a clear region of hysteresis. This stands out clearly in all of the plots, but most obviously in Figure 4f.  In the case of the upsweep in Figure 5a-g, dual state emission is observed (this is in the regenerated region). As K is increased to 0.645, shown in Figure 5h, the GS appears to oscillate with a very small amplitude suggestive of a Hopf bifurcation, before the onset of bursting and MMOs at K = 0.651 shown in Figure 5i. The lifetime of the bursting state decreases as the injection strength is further increased.
Similar MMOs are seen for high injection strengths during the downsweep until K = 0.645. Then, the presence of hysteresis can be seen by comparing the corresponding subplots of Figures 5 and 6. As K is further decreased below 0.645, we do not find continuous wave dual state emission until K = 0.605. Between these two values, the lifetime of the bursting state continues to increase and the switches to the quiescent phase become more rare. Eventually, trains of GS dropouts and ES pulses are observed, as seen in Figure 6b. These trains were observed for the entire 20 µs acquisition time of the oscilloscope. As there are fewer quiescent phases, the overall variance is large and is therefore represented as bright yellow in the variance map ( Figure 4f). Clearly, the yellow region extends beyond the white line marking the upsweep dynamical boundary, highlighting the hysteresis cycle. In Figure 6a, where K = 0.605, the oscillations disappear, and continuous wave dual state emission is obtained, as is also observed in the upsweep in Figure 5a-g. Thus, we find a bistability between the constant dual state emission (constant ES and regenerated injected light) and dynamic, pulsing behaviour with MMOs.

Slow Oscillations
As mentioned above, on the positive detuning side close to the transition between the phase-locked and regeneration regions, there is an area of dynamics indicated by a speckled, non-uniform colour along the diagonal boundary in Figure 2a,d (and also, of course, in Figure 4a,d). In fact, there are essentially two colours, indicating that there are two distinct average GS intensities. The binomial nature of the colouring suggests a slow dynamic of which we only sample a short, nearly constant part within the 20 µs measurement interval. The sampling window is increased to 1 second and new maps are created and shown in Figure 7. We note that there are also very narrow regions of hysteresis on both the left and right of this region. Figure 8 shows timetraces from an upsweep at −0.6 GHz marked by white lines in the upper panels of Figure 7. The period of the oscillations is extremely long: approximately 40 ms. The output is strongly dominated by the GS, but very small corresponding oscillations are also observed in the ES as can be seen in Figure 8. The cycle resembles a periodic oscillation between the phase-locked GS output and the dual state regenerated light output.
At the lower injection strength boundary of this region (i.e. at the bottom of the white line), the lower intensity section of the timetrace is longer lived than the upper intensity section, as is clear in Figure 8b. Increasing K, the duty cycle changes, eventually reaching 0.5, as seen in Figure 8h. Continuing to increase the injection strength, the higher GS intensity section becomes longer lived than the lower GS intensity section, as is most clearly seen in Figure 8k. Qualitatively, this behaviour is extremely similar to the optothermal square waves on the negative detuning side in the GS only system [14,16] and to the square wave/bursting phenomena reported here and originally in [15], where optothermal coupling breaks bistabilities. We thus interpret this region as another broken bistability between a low power GS output such as that shown in Figure 8a and a high power GS output such as that shown in Figure 8l. However, and importantly, the oscillations here are ∼100,000 times slower.
The qualitative similarity of the oscillations with the previously investigated optothermal effect is very suggestive. However, even apart from the vast difference in timescales, there is a significant difference that prevents the explanation of the oscillations via the original thermal coupling. Consider the phase locked bistability near the negative detuning boundary of [14]. In that case, the low power solution is on the negatively detuned side of the region and the high power solution on the positive side. The physical source of the coupling in that case is non-radiative carrier recombination. The low-power solution results in high carrier density and an ensuing higher temperature via increased non-radiative recombination. In the high power solution, the carrier density is lower and so the temperature is lower. Thus, higher power leads to lower temperature and vice versa; the deterministic cycle is driven by the low power solution being pushed towards positive detuning and the high power solution pushed towards negative detuning, yielding an anticlockwise phasor cycle as shown in [14,16]. However, in the new region identified in this work, the high power solution is the more negatively detuned of the two. Thus, the high power solution would push the system further towards the negative detuning direction and vice versa for the low power solution, and so the cycle would not arise. Thus, heating due to non-radiative recombinations cannot account for the phenomenon, and a coupling with the opposite sign is required. Such a coupling does exist, as can be seen by considering several important ways in which heating can arise in the device. In particular, reabsorption of light in the device can lead to such heating. In [28,29] both recombinative and re-absorptive heating was included. It is also known that self-absorption arises in QD lasers in particular, due to the wide distribution of states and dot sizes [30]. Thus, for re-absorption, a higher optical power leads to an increased temperature, and so the effect results in an effective detuning sweep in the opposite direction to that of the first optothermal effect, as required. Furthermore, the timescale for the second effect should be significantly longer as the time over which re-absorptive heating is distributed through the device can be of the order of tens of milliseconds, or even longer [28,31]. This matches the observed timescale extremely well. We see below that including such an effect in the model allows for excellent agreement with the experiment. The slow oscillations shown in panels (b-k) are between these two intensities, suggesting there is an underlying broken bistability. As K is increased the duty cycle moves from longer lived lower sections as most clearly seen in (b), through 50:50 in (f,g), and on to longer lived upper sections as most clearly seen in (k).

Theory
We use a rate equation model to numerically investigate the system. Rate equation models have been used extensively to great effect in studies of optically injected QD lasers, such as in [3,9,32,33] and further references within. Superficially, there are several different models. Of course, all arise from the same starting point: the Maxwell-Bloch equations. The differences then arise via choices about how to describe carrier capture and escape between the carrier reservoir, the excited state(s), and the ground state of the QDs and even whether to treat the electrons and holes separately or together in an excitonic model. Further, one can also consider whether to include effects such as phase-amplitude coupling and inhomogeneous broadening microscopically or phenomenologically, or indeed at all. We necessarily need to extend any model to include two optothermal effects with extremely slow timescales and so, for purposes of computational efficiency, we use an excitonic model for the carrier dynamics, while for the light dynamics we consider the electric field of the GS and the intensity of the ES. Our model is, where E GS is the electric field of the GS, I ES is the power of the ES, ρ GS and ρ ES are the occupation probabilities of the GS and ES, respectively, and N is the normalised chargecarrier number in the quantum well reservoir. The rate equation model in Equations (1)- (5) is thus a two state extension of that used in [3,9] where it was used to model GS only QD lasers under optical injection. The scattering terms in Equations (6)- (8), on the other hand, are of the form used in [32,33] and are given by There are three ∆ terms controlling the detuning in the system. ∆ 0 allows us to choose our reference frame. We define ∆ 0 ≡ δ 0 + δ with δ 0 chosen so that the GS frequency of the free-running laser is at zero. Thus, δ 0 compensates the frequency shift due to the phaseamplitude coupling and the optothermal effects in the free-running laser. Then, δ = 0 is at the central GS frequency of the free-running quantum dot laser and it is δ that defines the x-axis in each of our maps below. The other two terms, ∆ 1 and ∆ 2 , are used to implement the aforementioned optothermal coupling mechanisms. ∆ 1 accounts for heating of the active region due to recombinative heating, while ∆ 2 models overall device heating via re-absorption. We implement these two different effects with individual coupling strengths c 1,2 , and characteristic time scales γ 1,2 in the following two equations, where the subscript "th" denotes the value of the corresponding charge-carrier variable at the laser threshold. (We note that in [14,16] the thermal effect was coupled to the power even though the physical explanation is via the carriers. In the Class A case, this is unavoidable, as the carriers are adiabatically eliminated.) The meaning of all the other parameters used in the equations and their values are given in Table 1. It turns out to be quite instructive to first consider the system in the absence of any optothermal effects, then with only one effect, and finally with both effects included, thereby allowing for a comprehensive description of how the observed features emerge. There are many features that are common to all three cases. In particular, there is a large, central, phase-locked region; there are large dual state, regenerated, injected light regions for large magnitude detunings; and for low injection strengths at negative detuning, an antiphase dropout and pulsing dynamic is found in each case.
In the absence of any optothermal coupling, there are no square wave phenomena for either negative or positive detuning. Instead, we find multiple bistabilities near the negative detuning boundary. One of these is very clear from Figure 9 in the top left corner of each subplot. Consider the leftmost edge of Figures 9a,d. In the upsweep, the regenerated, two state output persists right up to the top of the figure. Figure 9c shows that the regenerated two state output eventually destabilises via a Hopf bifurcation along the bright line in the top left corner. There is only a very narrow region over which the resulting oscillatory solution exists, after which the system moves to the phase locked output. In the downsweep case, however, the bright yellow phase-locked region extends down to K∼4.8 at the extreme left of Figure 9d. Thus, there is a clear bistability between the phase-locked output and the regenerated two state output over a large region and, indeed, another bistability between the phase-locked output and the oscillating two state output over a small region. (d-f) show the corresponding downsweep. Bistabilty 1, the bistability between two phaselocked states, can be seen by comparing the top left corners of (a,d) at high injection strengths. It can also be seen clearly by comparing the top left corners of (b,e) at high injection strengths. Bistability 2, the bistability between a phase locked solution and the dual state output at moderate injection strengths, is clear when looking at the top left corners of the variance plots (c,f). Bistability 3, the low injection strength bistability between a constant (dual state) output and an oscillating (dual state) output, is visible at the bottom of the oscillating regime seen in (c,f). It is even clearer in the supplemental gif.
There is another region of hysteresis around K∼4. This is not as easy to see in the figure as it is over a much smaller area. We include a gif in the supplementary material showing the transition from Figure 9c,f, in which the hysteresis is clear. At this injection strength, the phase-locked solution first destabilises into the oscillating GS dropout/ES pulse regime reported in [24] as the injection strength is decreased. This phase-locked/pulsing regime is bistable with the regenerated two state output and its destabilising Hopf induced cycle discussed above (again, as is clear in the supplemental gif). Notably, though, there is not a bistability for positive detuning (where we find the extremely slow dynamic in the experiment).
One might reasonably expect the negative detuning bistabilities from previous GS only analyses [9,17], and some have even already been described in [15,21,32]. In [15], it was shown that the introduction of an optothermal coupling leads to the destruction of high injection strength bistability and that a periodic bursting dynamic is instead obtained. We repeat this analysis here with our model.

First Optothermal Effect
Here, we introduce the first optothermal effect, arising from non-radiative recombinations, and plot the updated maps in Figure 10. As expected, (and again, as was already shown in [15]) this breaks the upper bistability, and yields the periodic bursting square wave regime. The bistability between the phase locked/pulsing output and the dual state/Hopf dynamic is preserved although its precise location and size has changed. However, more intriguing here is the effect that this has on the positive detuning boundary. Now, where before there was no bistability, we find one. It is a bistability between the phase-locked output and the constant dual state output as indicated in Figure 11 near 1 GHz. This arises as the first optothermal effect shears the map so that the two solutions now overlap. It arises in a relatively narrow region precisely where the extremely slow square waves are obtained in the experiment. The appearance of this bistability motivates the introduction of a second optothermal coupling. As discussed in the experimental section above, this coupling must take the opposite sign to the original optothermal effect in order to yield the desired deterministic oscillations. The region of high variance seen in (c,f) is much larger than in Figure 9 as the large phase locked bistability is broken and is replaced with a square wave regime. γ 1 = 40 µs −1 , c 1 = 3 ns −1 .

Second Optothermal Effect
We use a characteristic timescale of γ 2 = 4 µs −1 for the second optothermal effect, which is not quite as slow as the experiment but allows for a reasonable computation time while maintaining a large enough timescale separation with the existing dynamic timescales. We have verified that an even longer timescale for γ 2 reproduces qualitatively identical dynamics. With the addition of the second optothermal effect, the square wave bursting solutions and the experimentally observed bistability are preserved as shown in Figure 12. We note that the value of c 1 when the two effects are included is different to the one we chose for just one optothermal effect. This is because the second optothermal effect counteracts the first one and so, in order to maintain the same overall shape, c 1 must change when the second coupling is introduced. On the positive detuning side, the new thermal effect yields a region of very slow oscillations in the region where there had been a bistability. The region over which these arise is most clearly seen when looking at the GS variance in Figure 12c,f and in the number of maxima as shown in Figure 13c,f in the narrow lines at the positive detuning boundaries. (In fact, there is even still a very small region of hysteresis as can be seen in Figure 11b in agreement with the experiment.) Figure 11. 1D bifurcation diagrams showing the bistabilities. The many diagonal lines are brought about by the bursting phenonema, where the amplitude is growing even though the input parameter for detuning remains fixed. (a) only one optothermal effect is included. A bistability is created on the positive unlocking boundary. At the negative unlocking boundary, there is a large cycle between the remnants of a bistability, with some of the bistability remaining intact. (b) both optothermal effects are included. The bistability at the positive boundary has been replaced with a limit cycle. In both cases K = 4.5.  Figure 13 is very instructive in highlighting the evolution of dynamical features with the inclusion of optothermal effects. We plot the number of maxima in the time series when there is no optothermal coupling in Figure 13a,d, when there is recombinative optothermal coupling in Figure 13b,e, and when both optothermal couplings are included in Figure 13c,f. Without optothermal coupling, there are very few regions of dynamics and those that exist do so only over small areas. When one optothermal coupling is included as in Figure 13b,e dynamics become much more prevalent in the system and a large region of dynamics appears in the central region of the map. In Figure 13c,f we find that the two competing optothermal effects endow the system with all of the experimentally observed dynamics. Most notably, as already mentioned, the experimentally observed extremely slow oscillations are found near the positive boundary.
Apart from the dynamical features, there are several distinct features in the shape of experimental map. One of these is the near vertical boundary for negative detuning for K∼0.6 in Figure 1. This feature is recovered with the inclusion of both optothermal effects as seen in the lower panels in Figures 12 and 13c,f for K ∼ 4. A second feature is the shape of the dynamical region (and consequently the phase locked region). In the experiment, there is a large and very asymmetric region of dynamics for negative detuning and a very narrow region of dynamics along the positive detuning boundary. Figure 13 shows the evolution of the dynamics arising in the model as the optothermal effects are included. There is one very narrow region of dynamics for negative detuning when no optothermal coupling is included. There is a large, somewhat symmetric, central region of dynamics when the recombinative optothermal effect is included. Finally, when both optothermal couplings are included, the map is sheared so that the dynamical region is rendered asymmetric, folding back towards negative detuning as the injection strength is increased, just as in the experiment.

Absence of Chaos
As mentioned in the introduction, there is a marked absence of chaos when a QD laser emitting from the ES is optically injected near the GS. We have not observed any chaotic regimes in the experiment or in the numerical studies. The presence of chaos in optically injected QD lasers for the GS only case has been reported several times both experimentally and theoretically, such as in [3,32,34] and we have confirmed it arises using our model by adjusting the parameters to allow for GS only lasing. However, it appears to be absent in the dual state system analysed in this work. While there are multiple maxima in large regions of the map as shown in Figure 13c,f, these correspond only to the bursting dynamics, which is complex but not chaotic.

Discussion
We have presented a detailed stability map for the optically injected dual state QD system. This is the first comprehensive experimental map in the literature to the best of our knowledge. We also present the first numerical map taking optothermal effects into account. These effects have dramatic effects on the overall system, but allow for excellent agreement between our experimental and numerical results. Unlike with the optically injected conventional semiconductor laser, optothermal effects play a somewhat prominent role in the map, and there is an array of intact and broken bistabilities due to competing optothermal effects. A particularly beautiful route to an extremely slow square wave dynamic near the positive detuning boundary is obtained. As with the implementation of optothermal effects in other systems, a bistability is broken yielding a deterministic optothermal cycle. However, the bistability only exists due to a different, competing, optothermal coupling in the first place! There is no sign or even hint of the dynamic in the bare system.
The need to consider optothermal effects in this system, and indeed in the conventional GS only system, arose entirely due to experimental findings where unexpected square wave and bursting phenomena were observed. However, the inclusion of optothermal effects in the model is shown here to improve the mapping beyond simply recovering the dynamics, with characteristic shapes and features found in the experiment recovered only when the full system is analysed. This suggests that perhaps including such couplings in other injection systems might be worthwhile.
With a view to the future, QD lasers are particularly suited to photonic integrated circuits, given their potential for isolator free operation on chip due to the enhanced stability to optical feedback [35][36][37]. Novel neuromorphic information processing may be possible using dual state QD lasers, as described recently in [38]. The stability map is an indispensable tool for analysing such applications, revealing the regions and parameter ranges over which dynamics are obtained. In fact, many of the dynamic regions obtained in this system are neuromorphic. The excitable dynamic reported in [24] persists at low injection strength and is of clear relevance while the neuromorphic potential of the optothermal bursting has previously been highlighted in [15].