Stability Boundaries in Laterally-Coupled Pairs of Semiconductor Lasers. Photonics (2), 74.

: The dynamic behaviour of coupled pairs of semiconductor lasers is studied using normal-mode theory, applied to one-dimensional (slab) and two-dimensional (circular cylindrical) real index conﬁned structures. It is shown that regions of stable behaviour depend not only on pumping rate and laser separation, but also on the degree of guidance in the structures. Comparison of results between normal-mode and coupled-mode theories for these structures leads to the tentative conclusion that the accuracy of the latter is determined by the strength of self-overlap and cross-overlap of the symmetric and antisymmetric normal modes in the two lasers.


Introduction
Arrays of laterally-coupled edge-emitting lasers (EELs) and vertical-cavity surface-emitting lasers (VCSELs) are used for many high-power or high-brightness applications [1,2]. Other potential applications include high-frequency modulation [3,4], beam-steering [5] and parity-time symmetry breaking [6]. The theoretical description of physical phenomena in laser arrays is usually based either on coupled-mode theory or normal-mode analysis. The former has the advantages of simplicity and physical insight, although it is not accurate for asymmetric structures [7], and can give erroneous results for anti-guided structures, as discussed in [8,9]. However, for symmetric structures with weak coupling between lasers, coupled-mode theory has been applied successfully to analyse the locking behaviour and dynamics of arrays with two or more elements [3,5,[10][11][12][13][14][15][16][17][18][19][20][21][22]. The normal-mode approach, which can provide a more accurate description for a wider range of structures, has been used in conjunction with a more sophisticated treatment of the electron-photon interaction to describe the beam switching and ultrafast pulsations in VCSEL arrays [23][24][25].
Comparisons of the coupled-mode and normal-mode methods have been made in terms of formalism [26] and of numerical results for bifurcations [17] for two-element arrays with real index guiding. In this situation, the coupling rate η, which in coupled-mode theory is calculated from the overlap integral of the lateral fields of the individual lasers [27], is related to the frequencies, ν s , ν a , of the symmetric and antisymmetric normal modes by η = (ν s − ν a )/2. Expressions for η in the case of purely real index guidance (i.e., ignoring any effects of gain or loss) have been derived for one-dimensional step-index (slab) waveguides [27] and for circular optical fibres that are weakly guiding (i.e., the difference between core and cladding refractive indices is much less than either index) [28]. A simplified expression for the latter that is valid for multimode fibre couplers has been given by Ogawa [29]. A more general expression for the coupling coefficient between circular cross-section VCSELs with real guidance is also available [30]. Comparison of the dynamics of a coupled pair of lasers modelled by slab waveguides indicated generally good agreement between coupled-mode and normal-mode treatments for edge-to-edge spacings greater than the waveguide width [17]. The only experimental test (to the best of our knowledge) of coupled-mode limits for laser pairs was performed by comparing predicted and measured far-field visibility of optically-pumped VCSELs [31]; the results indicated that coupled-mode theory was inaccurate for a spacing between the two pump spots of less than 13 µm when the modal radius of a solitary VCSEL was estimated as 3.5 ± 0.5 µm. In general, a clear definition of the ranges of parameters where coupled-mode theory is sufficiently accurate has not been given as yet.
The present contribution is intended to contribute to the discussion of dynamics in laterally-coupled pairs of EELs and VCSELs by presenting normal-mode theory and modelling for cases where the laser waveguides are one-dimensional (slab) or two-dimensional (circular cylindrical) structures. In each case, the guidance is purely real, so that no account is taken of guidance due to the effects of gain. We give results for the boundaries between regions of stable and unstable behaviour of these coupled pairs in the plane of normalised pump rate versus normalised spacing between the lasers. The normal-mode results are compared with those from coupled-mode theory in each case, as well as comparing between results for the slab and circular guides.

Model
The notation used for pairs of identical one-dimensional (slab) and two-dimensional (circular) waveguides is shown in Figure 1a,b, respectively. In each case, the core (cladding) refractive index is n core (n clad ) and the edge-to-edge separation is 2d. The slab half-width and cylinder radius are each a. With this notation, the conventional definitions of the normalised frequency v and the normalised decay constant of the fields in the cladding w for single solitary guides are: where λ is the free-space wavelength and n eff is the effective index of the single solitary guide. Ogawa [29]. A more general expression for the coupling coefficient between circular cross-section VCSELs with real guidance is also available [30]. Comparison of the dynamics of a coupled pair of lasers modelled by slab waveguides indicated generally good agreement between coupled-mode and normal-mode treatments for edge-to-edge spacings greater than the waveguide width [17]. The only experimental test (to the best of our knowledge) of coupled-mode limits for laser pairs was performed by comparing predicted and measured far-field visibility of optically-pumped VCSELs [31]; the results indicated that coupled-mode theory was inaccurate for a spacing between the two pump spots of less than 13 μm when the modal radius of a solitary VCSEL was estimated as 3.5 ± 0.5 μm. In general, a clear definition of the ranges of parameters where coupled-mode theory is sufficiently accurate has not been given as yet.
The present contribution is intended to contribute to the discussion of dynamics in laterallycoupled pairs of EELs and VCSELs by presenting normal-mode theory and modelling for cases where the laser waveguides are one-dimensional (slab) or two-dimensional (circular cylindrical) structures. In each case, the guidance is purely real, so that no account is taken of guidance due to the effects of gain. We give results for the boundaries between regions of stable and unstable behaviour of these coupled pairs in the plane of normalised pump rate versus normalised spacing between the lasers. The normal-mode results are compared with those from coupled-mode theory in each case, as well as comparing between results for the slab and circular guides.

Model
The notation used for pairs of identical one-dimensional (slab) and two-dimensional (circular) waveguides is shown in Figure 1a,b, respectively. In each case, the core (cladding) refractive index is ncore (nclad) and the edge-to-edge separation is 2d. The slab half-width and cylinder radius are each a. With this notation, the conventional definitions of the normalised frequency v and the normalised decay constant of the fields in the cladding w for single solitary guides are: where λ is the free-space wavelength and neff is the effective index of the single solitary guide. In what follows, we confine attention to the lowest-order normal modes in each of the structures of Figure 1, ignoring polarisation effects. Due to the symmetry of the waveguides, these modes have even parity, for the lowest order mode, and odd parity for the next highest. We refer to these as 'symmetric' and 'anti-symmetric', respectively, and are further distinguished by the fact that the antisymmetric mode goes through zero at the origin, whilst the symmetric does not. With the z-axis as the propagation direction, the total transverse electric field E(x,y,t) can then be expressed as: In what follows, we confine attention to the lowest-order normal modes in each of the structures of Figure 1, ignoring polarisation effects. Due to the symmetry of the waveguides, these modes have even parity, for the lowest order mode, and odd parity for the next highest. We refer to these as 'symmetric' and 'anti-symmetric', respectively, and are further distinguished by the fact that the anti-symmetric Photonics 2019, 6, 74 3 of 10 mode goes through zero at the origin, whilst the symmetric does not. With the z-axis as the propagation direction, the total transverse electric field E(x,y,t) can then be expressed as: where the suffix k denotes the symmetric (s) and antisymmetric (a) modes, ν k is the frequency of mode k, E k , Φ k are the time-dependent and spatial-dependent field components, respectively, and t is time.
For the one-dimensional slab waveguide of Figure 1a, the confinement is in the x direction, so that the dependence on y can be neglected.
In order to work with rate equations that are autonomous, we introduce the new field variable: The rate equation for Ẽ k is then: where κ is the cavity decay rate (=1/(2τ p ) where τ p is the photon lifetime), c is the speed of light, n g is the group index of the cavity, α is the linewidth enhancement factor, g 1 , g 2 are the mean gains per unit length in guides 1 and 2, and Γ 1kk' , Γ 2kk' are overlap factors in guides 1 and 2, defined as: where j = (1,2) labels the guide and the k,k' label the modes (s,a). The spatial components of the field are normalised as: The rate equation for carrier concentration N j in guide j is: where P j is the pumping rate in the jth guide and γ is the carrier recombination rate. The conventional linear relationship between gain, g j , and carrier concentration is assumed: where a diff is the differential gain and N 0 is the transparency concentration.

Results
For numerical calculations, we use a wavelength of 1.3 µm with core radius (half-width) 4 µm and refractive indices n core ≈ n clad = 3.4. Three values of the difference (n core − n clad ) are chosen as 0.000971, 0.002 and 0.0055. From Equation (1), these correspond to normalised frequency v = 1.571, 2.255 and 3.740, respectively. These values have been chosen to explore different regions of operation of the coupled guides in terms of the transverse modes supported by each solitary guide: the cut-offs for the first and second higher-order modes of the slab are v = π/2 and π, and those for the circular cylinder are 2.405 and 3.832, respectively. Hence, these values include operating regions, where 1, 2 or 3 transverse modes of the slab and 1 or 2 modes of the circular cylinder are present. The values of the First, we compare the coupling coefficient η for coupled circular cylindrical guides with the three values of v. We used industry standard software to numerically compute the normal modes and their complex two-dimensional field distributions. The results were benchmarked against published data wherever possible to confirm computational accuracy. For the case of symmetric one-dimensional slab waveguides, analytic solutions were used as the benchmark. Figure 2 shows the variation of η with the ratio d/a of edge-to-edge spacing to diameter for these three cases. The discrete points (indicated by circles, triangles and diamonds) are calculated from the difference of lowest-order normal-mode frequencies, η = (ν s − ν a )/2. The dashed lines are fitted to these results using the approximation due to Ogawa [29] of the form: where w is given by Equation (1) for the solitary guide. Also shown in Figure 2 (by square symbols) are the corresponding calculated results for a slab guide with v = 1.571. In this case, the dotted line is fitted to these results using the analytical form for the coupling rate [27], which yields η ∝ exp(−2wd/a). In all cases, the constant of proportionality is found by setting the functions equal to the normal-mode result at d/a = 1. There is a very good level of agreement between the numerical results and the analytic forms.
Photonics 2019, 6, x FOR PEER REVIEW 4 of 9 wherever possible to confirm computational accuracy. For the case of symmetric one-dimensional slab waveguides, analytic solutions were used as the benchmark. Figure 2 shows the variation of η with the ratio d/a of edge-to-edge spacing to diameter for these three cases. The discrete points (indicated by circles, triangles and diamonds) are calculated from the difference of lowest-order normal-mode frequencies, η = (νs − νa)/2. The dashed lines are fitted to these results using the approximation due to Ogawa [29] of the form: where w is given by Equation (1) for the solitary guide. Also shown in Figure 2 (by square symbols) are the corresponding calculated results for a slab guide with v = 1.571. In this case, the dotted line is fitted to these results using the analytical form for the coupling rate [27], which yields In all cases, the constant of proportionality is found by setting the functions equal to the normal-mode result at d/a = 1. There is a very good level of agreement between the numerical results and the analytic forms. Next, we turn our attention to the behaviour of the overlap factors defined in Equation (5) and evaluated for the lowest-order symmetric and antisymmetric normal modes. Figures 3-5 show plots of these versus d/a for coupled circular cylindrical guides with v = 1.571, 2.255 and 3.740, respectively. Each overlap factor is shown as a ratio to the conventional optical confinement factor ГS for the lowestorder mode (LP01) of the corresponding solitary single guide, given by [32]: where K0, K1 are modified Bessel functions. This ratio occurs when rate Equations (4) and (7) are cast into normalised form for computational purposes. In these figures, the subscript denoting the guide has been dropped and the values for subscripts ' k k ≠ are given as |Гsa|, since this quantity has a different sign in each guide. The symbols (circles, triangles, diamonds) in these figures refer to numerically calculated points, whilst the broken lines are empirical fits used later in the numerical solutions of the rate equations. It is noteworthy that there is significant structure in the variation of these overlap factors for low values of d/a in the operating regions considered here, and that the region of d/a where this structure occurs reduces with increasing v. Plots of the overlap factors for a specific slab waveguide have been presented in [17], including the detuning between resonant Next, we turn our attention to the behaviour of the overlap factors defined in Equation (5) and evaluated for the lowest-order symmetric and antisymmetric normal modes. Figures 3-5 show plots of these versus d/a for coupled circular cylindrical guides with v = 1.571, 2.255 and 3.740, respectively. Each overlap factor is shown as a ratio to the conventional optical confinement factor Γ S for the lowest-order mode (LP01) of the corresponding solitary single guide, given by [32]: where K 0 , K 1 are modified Bessel functions. This ratio occurs when rate Equations (4) and (7) are cast into normalised form for computational purposes. In these figures, the subscript denoting the guide has been dropped and the values for subscripts k k are given as |Γ sa |, since this quantity has a different sign in each guide. The symbols (circles, triangles, diamonds) in these figures refer to numerically calculated points, whilst the broken lines are empirical fits used later in the numerical solutions of the rate equations. It is noteworthy that there is significant structure in the variation of these overlap factors for low values of d/a in the operating regions considered here, and that the region of d/a where this structure occurs reduces with increasing v. Plots of the overlap factors for a specific slab waveguide have been presented in [17], including the detuning between resonant frequencies of the two lasers, and show somewhat less variation at corresponding d/a values for the case of zero detuning. In the limit of large d/a, all the ratios of overlap to confinement factors shown in Figures 3-5 tend to a value of 0.5. In this limit, it can be shown that rate Equations (4) and (7) reduce to the corresponding equations for the coupled-mode approximation [21]; details of this reduction as well as other aspects of the normal-mode treatment will be presented elsewhere.
frequencies of the two lasers, and show somewhat less variation at corresponding d/a values for the case of zero detuning. In the limit of large d/a, all the ratios of overlap to confinement factors shown in Figures 3-5 tend to a value of 0.5. In this limit, it can be shown that rate Equations (4) and (7) reduce to the corresponding equations for the coupled-mode approximation [21]; details of this reduction as well as other aspects of the normal-mode treatment will be presented elsewhere.
Using the results of Figures 2-5, together with similar results for slab guides, rate Equations (4) and (7) can be solved numerically and the regions of stable and unstable behaviour determined. Figures 6-8 show the results of this in terms of Hopf bifurcations in the plane of P/Pth vs. d/a, where Pth is the threshold value of P for a solitary laser.
At a given point in the plane, the steady state solutions were found using a nonlinear solver, started as close as possible to the expected solutions using approximate analytical expressions. These solutions yielded a Jacobian matrix, which was then diagonalised to find the eigenvalues. At a Hopf bifurcation, a pair of complex conjugate eigenvalues crosses the imaginary axis and the real part goes from negative (stable solution) to positive (unstable solution). A bisection routine was then written to trace the boundary between stable and unstable solutions using this criterion.   frequencies of the two lasers, and show somewhat less variation at corresponding d/a values for the case of zero detuning. In the limit of large d/a, all the ratios of overlap to confinement factors shown in Figures 3-5 tend to a value of 0.5. In this limit, it can be shown that rate Equations (4) and (7) reduce to the corresponding equations for the coupled-mode approximation [21]; details of this reduction as well as other aspects of the normal-mode treatment will be presented elsewhere.
Using the results of Figures 2-5, together with similar results for slab guides, rate Equations (4) and (7) can be solved numerically and the regions of stable and unstable behaviour determined. Figures 6-8 show the results of this in terms of Hopf bifurcations in the plane of P/Pth vs. d/a, where Pth is the threshold value of P for a solitary laser.
At a given point in the plane, the steady state solutions were found using a nonlinear solver, started as close as possible to the expected solutions using approximate analytical expressions. These solutions yielded a Jacobian matrix, which was then diagonalised to find the eigenvalues. At a Hopf bifurcation, a pair of complex conjugate eigenvalues crosses the imaginary axis and the real part goes from negative (stable solution) to positive (unstable solution). A bisection routine was then written to trace the boundary between stable and unstable solutions using this criterion.     Using the results of Figures 2-5, together with similar results for slab guides, rate Equations (4) and (7) can be solved numerically and the regions of stable and unstable behaviour determined. Figures 6-8 show the results of this in terms of Hopf bifurcations in the plane of P/P th vs. d/a, where P th is the threshold value of P for a solitary laser.  Figure 6 shows results for slab and circular cylindrical guides with v = 1.571, as well as the corresponding results obtained using the values of overlap factors in the limit of large d/a (labelled 'inf'). Figures 7 and 8 show the corresponding results for v = 2.255 and 3.740, respectively. The general form of the curves in each case resembles those reported in earlier literature on laterally-coupled lasers [11,13,21], although in some cases the branches of the curves at the lower left of each    Comparing the results in Figure 6-8, it is clear that with increasing values of v, the regions of instability shrink with their limits moving to lower ranges of d/a. This trend follows that of the overlap factors in Figure 3-5 in their deviation from the value of 0.5 at small d/a. Also, with increasing v, the curves for the slab and circular cylindrical guides become closer, illustrating the power of this normalisation. Attempts to achieve a closer match between these curves by choosing differing v values (in terms of core-cladding index difference) for each, based on matching either the w value or the coupling rate slope with d/a for the slab and circular guides, met with limited success.
Comparing the curves computed from the normal-mode equations with those from the coupledmode approximation (equivalent to the curves labelled 'inf') in Figures 6-8 indicates that in most cases the latter are at lower d/a ranges than the former. For the slab waveguide at the largest value of v (in Figure 8), the agreement between coupled-mode and normal-mode results is rather closer than for the other values, and this may result from the limited range over which the normalised overlap factors deviate from 0.5. It is worth bearing in mind also that higher-order modes (1 for the circle, 2 for the slab) are present in the guides of Figure 8, but we have considered only coupling between the fundamental modes of each guide (as represented by the lowest-order symmetric and antisymmetric normal modes of the coupled structure). The theory of coupling between higher-order modes in coupled optical fibres, including the effects of spin-orbit interaction, has been presented in [33].

Discussion
The normal-mode theory for laterally-coupled pairs of identical semiconductor lasers has been given in terms of rate equations for the lowest-order symmetric and antisymmetric modes in structures with real index guidance. Both one-dimensional (slab) and two-dimensional (circular At a given point in the plane, the steady state solutions were found using a nonlinear solver, started as close as possible to the expected solutions using approximate analytical expressions. These solutions yielded a Jacobian matrix, which was then diagonalised to find the eigenvalues. At a Hopf bifurcation, a pair of complex conjugate eigenvalues crosses the imaginary axis and the real part goes from negative (stable solution) to positive (unstable solution). A bisection routine was then written to trace the boundary between stable and unstable solutions using this criterion. Figure 6 shows results for slab and circular cylindrical guides with v = 1.571, as well as the corresponding results obtained using the values of overlap factors in the limit of large d/a (labelled 'inf'). Figures 7 and 8 show the corresponding results for v = 2.255 and 3.740, respectively. The general form of the curves in each case resembles those reported in earlier literature on laterally-coupled lasers [11,13,21], although in some cases the branches of the curves at the lower left of each Figure are not reported. Stable solutions for the antiphase mode are found to the upper right of the upper branches (away from the origin) and to the lower left of the lower branches of the curves (near the origin); in other regions, unstable steady state solutions are found. The results labelled 'inf' in Figures 6-8 are in perfect agreement with those from coupled-mode theory, using the approximation for the Hopf bifurcation [21]: Comparing the results in Figures 6-8, it is clear that with increasing values of v, the regions of instability shrink with their limits moving to lower ranges of d/a. This trend follows that of the overlap factors in Figures 3-5 in their deviation from the value of 0.5 at small d/a. Also, with increasing v, the curves for the slab and circular cylindrical guides become closer, illustrating the power of this normalisation. Attempts to achieve a closer match between these curves by choosing differing v values (in terms of core-cladding index difference) for each, based on matching either the w value or the coupling rate slope with d/a for the slab and circular guides, met with limited success.
Comparing the curves computed from the normal-mode equations with those from the coupled-mode approximation (equivalent to the curves labelled 'inf') in Figures 6-8 indicates that in most cases the latter are at lower d/a ranges than the former. For the slab waveguide at the largest value of v (in Figure 8), the agreement between coupled-mode and normal-mode results is rather closer than for the other values, and this may result from the limited range over which the normalised overlap factors deviate from 0.5. It is worth bearing in mind also that higher-order modes (1 for the circle, 2 for the slab) are present in the guides of Figure 8, but we have considered only coupling between the fundamental modes of each guide (as represented by the lowest-order symmetric and antisymmetric normal modes of the coupled structure). The theory of coupling between higher-order modes in coupled optical fibres, including the effects of spin-orbit interaction, has been presented in [33].

Discussion
The normal-mode theory for laterally-coupled pairs of identical semiconductor lasers has been given in terms of rate equations for the lowest-order symmetric and antisymmetric modes in structures with real index guidance. Both one-dimensional (slab) and two-dimensional (circular cylindrical) guides have been considered. Attention has been drawn to the important role played by the overlap factors between the modes in each guide and their variation with guide separation. In the limit of large separation, these overlap factors, when normalised by the confinement factor of the fundamental mode in a solitary laser guide, tend to the value 0.5. In this limit, the normal-mode rate equations reduce to those for the coupled-mode approximation with the corresponding result that the coupling rate is given by half the difference in frequencies of the symmetric and antisymmetric modes.
Results for the boundaries between regions of stability and instability in the plane of normalised pump rate versus normalised guide spacing have been presented for slab and circular guides with three different values of (weakly-guiding) core-cladding refractive index difference and all other parameter values the same. These results indicate that the ranges of instability shrink with increasing index difference and their boundaries move to lower values of guide spacing. The coupled-mode results become closer to those from the normal-mode theory with increasing index difference and increasing guide spacing. This behaviour is attributed to the corresponding changes of the normalised overlap factors with these parameters, in particular their deviation from the limiting value of 0.5. It is, therefore, postulated that the behaviour of the normalised overlap factors determines the level of accuracy of the coupled-mode theory, a property that has hitherto not been properly defined.
Future work in this subject will include extension of the normal-mode theory for coupled lasers to include the influence of polarisation and the effects of higher order modes.