Dynamics of Evanescently-Coupled Laser Pairs With Unequal Pumping: Analysis Using a Three-Variable Reduction of the Coupled Rate Equations

—The ﬁve coupled rate equations used to describe laterally-coupled pairs of lasers with weak coupling and unequal pumping are reduced to a new system of three equations. This enables approximate closed-form steady-state solutions and explicit expressions for the boundaries between regions of stable and unstable dynamics to be found. The results of applying these approximations to speciﬁc cases of coupled laser pairs are shown to be in good agreement with results obtained from numerical solutions of the original set of ﬁve equations as well as earlier resultsfromtheliterature.Inadditiontheapproximationsbasedon thereducedsetofequationsallowasystematicinvestigationoftheeffectsofmaterial,deviceandoperatingconditionsontrendsand novelfeaturesinthedynamicsoflaterally-coupledlaserpairs.Thealgebraicresultsgiveinsightintotrendswithparameterswithout theneedforextensivenumericalcomputationandshouldthereforebeofuseinmodellingtwo-elementVCSELarraysfornumerous potentialapplications.

shown to exhibit non-Hermiticity [12], [13], exceptional points [13]- [15] and parity-time symmetry breaking [1].Theoretical understanding and modelling of these effects has been largely based on coupled rate equations [1], [2], [10]- [12], [14]- [16].These apply in the general case of unequal pumping and have also been extended to include asymmetry in coupling coefficient and photon and carrier lifetimes [17].Since there is a dearth of analytical solutions, numerical methods have largely been necessary to explore the importance of key physical design and material parameters.While numerical methods have proved effective in many cases, analytic solutions can often provide more accessible routes and clearer insights towards understanding the system.For the particular case of weakly-coupled 2-element arrays with equal pumping, approximate steady-state solutions as well as analytic expressions for the boundaries of regions of stability have been presented [16].Thus the objective of the present contribution is to extend that work so as to provide approximations for these quantities in the case of unequal pumping in weakly-coupled laser pairs with symmetric real coupling.
This paper is organised as follows.Section II presents the main theoretical results in terms of a reduction of the five coupled rate equations to a new system of three equations, as well as the approximate algebraic solutions of these equations and three stability conditions for their dynamics.Details of the derivations of these results are given in Appendices A and B. Section III presents first some tests of these results against published work and then proceeds to give some examples of the use of the approximate stability boundaries for a specific example.The accuracy of these approximations is demonstrated by comparison with stability maps calculated using the method of Langrangian descriptors [18] as well as by numerical integration of the rate equations combined with a bifurcation method [19].It is shown that the approximate solutions permit novel features of the dynamic maps to be identified and their dependence on coupling coefficient and pumping asymmetry to be explored in a systematic manner.These features are summarised in the concluding section IV.

II. THEORY
We begin with the rate equations for a pair of coupled lasers whose coupling rate η is real [16]: This work is licensed under a Creative Commons Attribution 4.0 License.For more information, see https://creativecommons.org/licenses/by/4.0/ where Y A , Y B are the normalised fields and M A , M B are the normalised carrier densities in lasers A, B, respectively, φ is the phase difference between the fields in B and A, ΔΩ is the detuning between the cavity resonances of lasers B and A, τ N is the carrier lifetime, τ p is the photon lifetime, α H is the linewidth enhancement factor and Q A , Q B are the normalised pumping rates.We assume weak coupling, hence ητ p << 1.Now, following earlier approaches for VCSELs [20], [21], ring lasers [22] and spin-VCSELS [23], assume that τ p << τ N and 1) -( 4) can be reduced to three rate equations and one conservation relation, as follows: The new variables are defined by q ), and the details of the derivation of these equations are given in Appendix A. Equation ( 8) is an energy conservation law that holds on the timescale of the carrier recombination time, whilst shorter timescale dynamics are included in ( 5)- (7).The steady state solutions of ( 5) -( 7) can be found explicitly by using the approximation sin ψ s ∼ = −q which is consistent with ητ p << 1, m s << 1, where the subscript 's' denotes the steady-state value.The results can then be written as In the limit of equal pumping Q A = Q B , q = 0, (9) and ( 11) reduce to the forms in (25), (26) and (29) of [16], whilst (10) reduces to the forms of ( 27) and (28) of [16] but without the terms of order ητ p which appear in the latter.Since we assume ητ p << 1, this omission should only result in a very small error in the accuracy of (10).
By performing a small-signal analysis of ( 5) -( 7) (see Appendix B), approximate expressions for the stability boundaries of the system can be found as: Equation (15), in combination with (11), corresponds to the Hopf bifurcation.Equation ( 14) describes the saddle-node (SN) bifurcation.For this upper limit of detuning, (12a) and (12b) show that the phase values for the in-phase and out-of-phase solutions are equal: In the limit of equal pumping the results of ( 14) and ( 15) can be reduced to the forms of (30) and (31) of [16] for real coupling rate, i.e., (17) where Equations ( 5) - (15) are the main results of this paper.In the following we test them against results in the literature and give some numerical examples of their use.The accuracy of results for these examples is verified using different numerical methods of solution.

A. Comparison With Results From the Literature
The accuracy of the approximations used in deriving ( 9) -( 15) can be tested by comparing some results with those in the literature.First, for the steady-state solutions we compare with the results of Erneux and Lenstra [24] in their special case of zero time delay of mutually delay-coupled quantum cascade lasers.These authors use similar assumptions of weak coupling to reduce the six rate equations for that system to two coupled equations for the phase of the electric fields; for the case of zero time delay these can be combined into a single equation.It is straightforward to show that our (10) and ( 11) can be combined in the form where θ o = tan -1 (α H ) -π/2.This result matches the steadystate solution of the corresponding phase equation ( 12) of [24], allowing for the differences in notation.
Contours of constant phase φ s+ , φ s-and field ratio Y B /Y A in the plane of Q B versus ΔΩ for fixed Q A , calculated from a numerical root searching technique, have been presented by Gao et al in [12].The steady-state results ( 10) and ( 12) can be used to allow comparison with these; similarly the boundaries of stable operation can be found from ( 14) and (15).Using the parameter values corresponding to the weakly-coupled "array 1" (ητ p = 0.002) of [12] with Q A = 3.2, the results of applying our equations to find tilted in-phase and out-of-phase solutions show good agreement for contours of constant phase.Encouraging agreement is also found for the limits of stability using the SN bifurcation approximation from (14).In addition we calculated the Hopf bifurcation using ( 15) and (11).For this we needed the value of τ N which is not given in [12]; since in other publications, e.g., [8], [17], [25], these authors use τ N = 2 ns, that value was assumed here.Our results indicate that only the out-of-phase case shows stability and moreover the range of stable solutions, as bounded these bifurcations, is very limited for this array.Stability is restricted to two narrow regions: one at negative detuning for Q B > Q A with phase close to 3π/2 and the other at positive detuning for Contours of constant phase as well as boundaries of stable operation in the plane of normalised detuning versus normalised pumping difference are presented by Kominis et al in [14] for very weak coupling (ητ p = 0.000315).Hence for a further test of our approximate results ( 13) -( 15) we used the same parameters as [14].The plot that appears as Figure 9(b) of [14], calculated by utilising a numerical continuation algorithm, applies to the range (0 ≤ q ≤ 1, 0 ≤ -ΔΩτ p /0.05 ≤ 0.15) in the plane of normalised detuning -ΔΩτ p /0.05 (ࣕ 'Δ' [14]) versus q (ࣕ 2 x'ΔP' [14]).For this range our results for the tilted out-of-phase solutions of this system are in very good agreement with those presented in [14].Note that comparison of our ( 3) with ( 6) of [14] indicates that the minus sign in the detuning is needed for comparison of results.Comparison of the results discussed above using parameters from [12] and [14], leaving aside the difference of axes related to pumping rate, reveals a significant difference between the topology of the regions of stability for each case.As mentioned above, in the case of parameters from [12] there are separate distinct regions of stability for positive q (Q A > Q B ) and negative q (Q A < Q B ), each bounded by a pair of SN and Hopf bifurcations.In contrast, we find for parameters from [14] there is a single stable region that includes both negative and positive q and is bounded by two SN and two Hopf bifurcations.This difference can be traced to the difference in relative sets of parameters and can be quantified in terms of a critical value η c of the coupling rate given by the condition that the square root in (17) vanishes for the case of equal pumping, i.e., For values of coupling rate that are greater than or equal to η c the argument of the square root in (17) is greater than zero and the regions of stability are similar to those for parameters of [12] where η = 1 and η c = 0.2 with our assumption of τ N = 2 ns.For values of coupling rate that are less than η c the regions of stability are similar to those for parameters of [14] where ητ p = 0.000315 and η c τ p = 0.0005.

B. Numerical Results for a Specific Example
Consider the case of real index slab wave-guiding [16] with α H = 2, τ N = 1 ns, τ p = 1.53 ps and η = 83.6 exp(-2.52d/a) ns -1 where 2d is the edge-to-edge separation of the laser waveguides and 2a is the width of each.The normalised pumping rate Q A,B in each laser is related to the physical pumping rate P A,B by Q j = 11.4(Pj /P jth -1) + P j /P jth where the subscript 'th' denotes the value at lasing threshold.We will consider only the tilted out-of-phase solutions since these are the only ones allowed for these parameters at q = 0 [16].First we consider the case of laser separation given by d = 1.5a, which yields η = 1.908 ns -1 (ητ p = 0.00292), with Q A + Q B = 26.8 which corresponds to both lasers at twice threshold when q = 0. Fig. 1 shows a stability map in the plane of linear frequency detuning ΔΩ/2π Fig. 2. Stability maps in the plane of normalised detuning versus positive q for the same parameters as in Fig. 1, calculated using the method of Lagrangian descriptors: (top) for the 3 reduced equations ( 5) - (7), and (bottom) for the full set of equations ( 1) -( 4).Blue colour denotes the region of stable operation.The stability boundaries from ( 13) -( 15) are also shown.
versus q, where the dotted red, solid black and broken black lines correspond to the stability boundaries defined by ( 13), ( 14) and (15), respectively.The region of stable operation is bounded by two SN and two Hopf bifurcations, as defined by ( 14) and (15).This is because the critical value of coupling η c for these parameters is 3.35 ns -1 which is greater than the value of 1.908 ns -1 used to calculate Fig. 1.There are two points in Fig. 1 where the 3 boundary lines touch tangentially; these SN-Hopf points are at (0.7072, -1.822 GHz) and (-0.7072, 1.822 GHz).It is easy to show that these points are given by tanφ s = α H /q.
In order to verify the accuracy of our approximations, stability maps have been calculated using the method of Lagrangian descriptors (LDs) [18].These are shown in Fig. 2 for the reduced equations ( 5) -( 7), and the full set of equations ( 1) -( 4).Only the positive q half-plane is shown in each case.The colour shading indicates the contours of the Lagrangian descriptor for the system.Bifurcations are indicated by abrupt changes in these contours.Fig. 2(top) thus tests the accuracy of the approximation sin ψ s ∼ = −q which is used in deriving the stability boundaries Fig. 3. Stability map with boundaries from ( 13) -( 15) in the plane of normalised detuning versus q with η = 3.35 ns -1 and other parameters as for Fig. 1.Blue colour denotes the regions of stable operation.( 13) -( 15) (see Appendix B).These boundaries are seen to be in good agreement with the LD result in defining the region of stable operation (denoted by blue colour).Fig. 2(bottom) verifies the accuracy of the reduced set of equations ( 5) -( 7) since there is good agreement with Fig. 2(top), the only small difference occurring in the position of the Hopf bifurcation at its region of lowest q (around 0.5).
For the next example, we take the coupling rate to be the critical value of 3.35 ns -1 (ητ p = 0.00513), keeping the other parameter values the same as those for Fig. 1.This value of η corresponds to a spacing of d = 1.2766a in the model of coupled slab waveguides [16].Fig. 3 shows the boundaries defined by ( 13) -( 15) superimposed on an LD stability map using (4) - (7).Only the positive q half-plane is shown since the negative half-plane can be found by simply reversing the signs of the axes.There are two regions of stable behaviour (in blue) each of which is enclosed by saddle-node and Hopf bifurcations, and the latter boundaries intersect at q = 0.In this case the points of contact of lines defined by ( 13) -( 15) are at (0.5818, -2.7308 GHz) and (-0.5818, 2.7308 GHz).There is a good level of agreement between the approximate boundaries and those found from the LD method with an error increasing at values of q close to 1.
To complete this set of results we consider a higher rate of coupling given by η = 6.726 ns -1 (ητ p = 0.0103), corresponding to a spacing of d = a, i.e., the edge-to-edge spacing is equal to the full width of each waveguide.Fig. 4 shows two regions of stable behaviour each of which is enclosed by SN and Hopf bifurcations which meet at the points (0.3977, -4.7585 GHz) and (-0.3977, 4.7585 GHz).Again the agreement between approximation and LD boundaries is good except at higher values of q.
Another way to present these results is in terms of detuning versus normalised laser spacing d/a for various values of the normalised pumping difference q.Fig. 5 shows plots of this type for q values of 0.4 and 0.8, using the same parameters as those used previously.In Fig. 5(top) the condition for equality of the boundaries occurs for d/a = 1 at q = 0.3977 as mentioned in the  13) -( 15) in the plane of normalised detuning versus q with η = 6.726 ns -1 and other parameters as for Fig. 1.Blue colour denotes the regions of stable operation.Fig. 5. Stability maps and boundaries from equations ( 13) -( 15) in the plane of normalised detuning versus normalised spacing d/a for q = 0.4 (top) and 0.8 (bottom) and other parameters as for Fig. 1.The scale on the colour bar shows the numbers of extrema in the time series.discussion of Fig. 4; hence for the value of q = 0.4 used here the curve for ( 13) is barely observable.For q = 0.8 in Fig. 5(bottom) the change of all boundaries is clearly seen, as expected from the discussion of Figs. 1 -4.
In order to test the accuracy of the approximations used to produce Fig. 5, a numerical (Runge-Kutta) solution of the rate equations (1) -( 4) was used.From the extrema of the time series of Y A 2 + Y B 2 + 2Y A Y B cos φ at each value of d/a and detuning, one-parameter bifurcation diagrams are used to construct the stability maps in Fig. 5. Examples of one-parameter bifurcation diagrams are given in Figs 5 and 9 of [16] and details of how these are used to construct stability maps are given in [19].White colour is used to denote the regions of stability; other colours denoting regions of period 1, 2 and 3 oscillation and chaos are marked in Fig. 5(top).It is clear that there is a very good level of agreement for the larger values of d/a, whilst at smaller values the approximate results tend to underestimate the frequency range of the stable region.The SN bifurcations agree well and the underestimate is associated with error in the Hopf curves at higher coupling rates.This is thought to be associated with the small error (of order ητ p ) which was noted in the discussion of (10) in section II and implies an error of the same order in the approximation sin ψ s ∼ = −q.Thus, based on this evidence and that from the comparisons of results in Fig. 4, we estimate that the approximation of ( 15) for the Hopf bifurcation is very accurate provided ητ p is less than about 0.005.No such limit appears to apply to the approximation of ( 14) for the SN bifurcation which retains its accuracy over the entire range tested here.

IV. CONCLUSION
In this paper we report on how the five coupled rate equations which describe laterally coupled laser pairs with weak real coupling (ητ p << 1) can be reduced to a system of three equations The variables in these equations are the phase difference between the fields in the two lasers, the ratio of the field amplitudes (expressed as an angle) and the difference in normalised excess carrier densities in the two lasers.The underlying physical parameters are the coupling rate η, the linewidth enhancement factor α H , the carrier lifetime τ N , the photon lifetime τ p , the detuning ΔΩ between the cavity resonances of the lasers and the normalised pumping rates expressed in terms of their sum ).The approximate steady-state solutions of this reduced set of equations have been tested against published results which were found using numerical methods.Additionally, the use of a small-signal analysis has revealed closed-form expressions which predict the boundaries of stable operation.These have also been compared against published results and our own numerical solutions for the case of weakly-coupled laser pairs, and good agreement between the numerical methods and those from the approximations are found.Finally, in the limit of equal pumping the algebraic approximations are shown to reduce to those already known from earlier work [16], [26].
The new closed-form expressions allow a systematic investigation of the dependence of the dynamics on parameters such as pumping rates, coupling rate and linewidth factor.Results have been presented for a specific example showing the effect of varying η and q whilst keeping (Q A + Q B ) fixed.These results show novel effects in terms of the asymmetric behaviour of the stability boundaries and how this develops as the parameter values are changed.In all cases analysed it is found that pairs of saddle-node and Hopf bifurcations touch tangentially at points whose co-ordinates can be found algebraically.In addition the new expressions offer a route for future exploration of the various forms of nonlinear dynamics that exist in a system of coupled laser pairs in terms of physical designs and materials parameters, and via external control using differential pumping.Such results should be useful in modelling 2-element VCSEL arrays for applications such as beam steering, enhanced modulation response, etc., since trends with various parameters can easily be tracked without the need for extensive numerical simulation.
The approach used here follows a method applied previously to VCSELs [20], [21], [23] and ring lasers [22] in that a reduction in the number of rate equations is achieved by assuming a power conservation law and a new variable that is related to the ratio of the field amplitudes.This approach would also be applicable to coupled nanowire lasers which have been predicted to have potential for ultra-high frequency modulation [27].Whilst the present contribution has been limited to the case of real coupling rates, it is straightforward to extend this to take account also of the phase of the coupling in order to give a more general description of the effects of gain-guiding and index antiguiding [16].Further generalisation to deal with larger arrays of weaklycoupled lasers, including two-dimensional arrays of VCSELs, is also possible but in all cases it is only possible to reduce the number of rate equations by two.

APPENDIX A DERIVATION OF THE REDUCED SET OF COUPLED RATE EQUATIONS
Defining , adding and subtracting the two versions of (4) yields For the situation of ητ p << 1 we will assume that (A1) can be replaced by a conservation law that holds for dynamic time scales longer than the relaxation oscillation period: We have thus assumed that the time derivative in (A1) is zero and that the final two terms on the RHS can be neglected.We show below in (A5) that the penultimate term is zero.
From (A3) it follows that Using ( 1) and ( 2) this result leads to Following the example of [22], we define a new variable ψ by using again the conservation law (A3): It follows that (A5) can be written as Using this result and the definitions (A6), (A2) becomes Also, (3) becomes From the definitions (A6), we have Hence ( 1) and ( 2) can be transformed to Multiplying (A13) by tan[(ψ + π/2)/2] and (A14) by cot[(ψ + π/2)/2], and adding gives Defining new variables m = m Am B and q = (Q A -Q B )/ Π, (A10), (A11) and (A15) can be written as These are the reduced set of equations given as ( 5) - (7) in the main text.The steady state versions are where the subscript 's' denotes the steady-state value.Since we assume that ητ p << 1, m s << 1, it follows that a good approximation for the solution of (A19) is Using this approximation in (A21) and (A8) it follows that In the limit of equal pumping q = 0, these results reduce to the expressions given in [16] in the case of real coupling.In the limits of q tending to 1 or -1, which correspond, respectively, to Q B = 1 or Q A = 1, i.e., one or other laser at threshold, the result in (A23) for the laser above threshold suggests that the quantity m j tends to infinity.However, this is also the case for the steady-state solutions of (1) and (2) since at threshold the field amplitude is zero.
For zero detuning , in the tilted in-phase case the result is For zero detuning , in the tilted antiphase case the result is For equal pumping, (B16) reduces to η < Q/(2α H τ N ), also first derived in [26].Note that these and other authors use the notation P = (Q -1)/2 for the normalised excess currents.

Fig. 4 .
Fig. 4. Stability map with boundaries from (13) -(15) in the plane of normalised detuning versus q with η = 6.726 ns -1 and other parameters as for Fig.1.Blue colour denotes the regions of stable operation.