The Overlap Factor Model of Spin-Polarised Coupled Lasers

: A general model for the dynamics of arrays of coupled spin-polarised lasers is derived. The general model is able to deal with waveguides of any geometry with any number of supported normal modes. A unique feature of the model is that it allows for independent polarisation of the pumping in each laser. The particular geometry is shown to be introduced via ’overlap factors’, which are a generalisation of the optical conﬁnement factor. These factors play an important role in determining the laser dynamics. The model is specialised to the case of a general double-guided structure, which is shown to reduce to both the spin ﬂip model in a single cavity and the coupled mode model for a pair of guides in the appropriate limit. This is applied to the particular case of a circular-guide laser pair, which is analysed and simulated numerically. It is found that increasing the ellipticity of the pumping tends to reduce the region of instability in the plane of pumping strength versus guide separation.


Introduction
Detailed control of the polarization dynamics of vertical cavity surface emitting laser (VCSEL) arrays opens up many technological opportunities. In addition to being highpower sources, potential applications may lie in emerging areas, such as reservoir computing based on polarization dynamics [1] and secure key distribution based on chaos synchronization [2,3], as well as an enhanced understanding of anticipation in polarisation chaos synchronisation of coupled VCSELs [4] for secure communications.
Key to this development is the understanding of both polarisation dynamics in a given laser and the evanescent coupling to nearby elements in the laser array. The dynamics in a single laser are governed by the interplay between spin-polarised carriers and the circularly polarised components of the optical field. The spin flip model (SFM) [5] is now wellestablished as providing an excellent quantitative description of these effects. However, despite extensive developments of the SFM, it remains restricted to the description of single lasers.
The dynamics of coupled lasers are typically studied using a coupled mode model. However, this approach has the shortcoming that it does not allow for independent polarisation of the pump in each laser, whether this is optical or via spin-polarised injection current. Nor does it allow us to investigate the spatial variation in the optical polarisation. What is required to deal properly with this scenario is a supermode (also known as normal mode) model, in which one finds the optical solutions of the entire array. This also allows for the spatial distribution of the polarisation to be studied. Whilst supermode studies of laser arrays have been undertaken [6], to our knowledge none have incorporated the optical polarisation or carrier spin-dynamics . This current work is intended to remedy this deficit, in which we extend the SFM to apply to a general supermode model for laser arrays.
The basic SFM consists of four coupled rate equations (two for spin-polarised carriers and two for polarised field components) and it includes rates of carrier recombination, photon field decay, and electron spin relaxation (spin relaxation of holes is usually assumed to be instantaneous). The nonlinear dispersion that couples the carrier concentrations to the phases of the optical fields is described by the linewidth enhancement factor, and the field interactions due to nonlinear anisotropies are included via rates of birefringence and dichroism. For conventional VCSELs, the SFM has been applied to explain the experimental results of polarisation switching (PS) [7,8]. An "extended SFM" [9] that accounts for thermal effects and includes a realistic spectral dependence of the gain and the index of refraction of the QWs has been used [10] to explain the experimental results on elliptically polarised dynamical states that occur in the polarisation dynamics of VCSELs in the vicinity of one type of PS. For a more complete discussion of polarisation dynamics in VCSELs, the reader is referred to [11].
A further development of the extended SFM [12] includes a description of the spatial variation of the electromagnetic modes and the carrier densities. The variation in the longitudinal direction is dealt with by integration over the length of the VCSEL cavity, whilst the radial and azimuthal variation is described by accurate solutions of the wave equation. The model assumes a given functional dependence of the guiding mechanisms (built-in refractive index and thermal lensing) as well as the spatial dependence of the current density. The transverse mode behaviour of gain-guided, bottom, and top-emitter VCSELs were studied, and it was shown that the stronger the thermal lens, the stronger the tendency toward multimode operation, which indicates that high lateral uniformity of the temperature is required in order to maintain single mode operation in gain-guided VCSELs. Additionally, close-to-threshold numerical simulations showed that, depending on the current profile, thermal lensing strength and relative detuning, different transverse modes could be selected.
Another version of the extended SFM [13] includes a rate equation for the temperature of the active region, which takes decay to a fixed substrate temperature, Joule heating, and heating due to non-radiative recombination into account. The temperature dependence of the PS point is characterised in terms of various model parameters, such as the roomtemperature gain-cavity offset, the substrate temperature, and the size of the active region.
The SFM has also been widely applied to describe the behaviour of spin-VCSELs whose output polarisation can be controlled by the injection of spin-polarised electrons using either electrical or optical pumping (for a review with more details, see [14]). In the latter case the polarisation of the optical pump is included [15] to reveal its effect on the output polarisation [16][17][18]. The SFM has also been used [19][20][21] to explain the experimental results on high-speed polarisation oscillations that result from competition between the spin-flip processes, dichroism, and birefringence.
From this brief summary of the SFM and its applications, it is clear that the structures studied have been limited to single lasers, either conventional electrically driven VCSELs or spin-VCSELs that may be pumped electrically or optically. In the present contribution, we seek to extend the range of application to include structures, where two or more evanescently-coupled lasers are arranged in parallel to form arrays with the possibility of different lasers having differing pumping polarisation. To the best of our knowledge, this configuration has not been analysed previously, although there is one report [22] of an experiment where optical pumping with orthogonally polarised beams was used to study the interaction between two VCSELs as a function of their separation. There is considerable literature on laser arrays because of their important practical applications as high-power sources (including, most recently, for three-dimensional (3D) sensing in smartphones [23]) and very sophisticated models of VCSEL arrays have been developed [24]. Arrays of coupled lasers are also of fundamental interest in view of the range of nonlinear dynamics that they can exhibit (see, for example, [25] and references cited therein). Although there is sometimes a need to stabilise the polarisation of such arrays of VCSELs, the possibility of manipulating the output polarisation of an array by means of independent pumping polarisations has not yet been considered. Moreover new potential applications of the theory developed here may lie in emerging areas, such as reservoir computing based on polarization dynamics [1] and secure key distribution based on chaos synchronization [2,3], as well as an enhanced understanding of anticipation in polarisation chaos synchronisation of coupled VCSELs [4] for secure communications. Hitherto, these topics have been studied by adding delayed optical injection terms into the SFM equations [26], whereas the present theory permits a more general analysis which is applicable to a much wider range of VCSEL array applications. Thus, these developments, together with the issue of how the array dynamics is affected by spin-polarised pumping, provide the motivation for the present study (the reader is also referred to complementary work reported in [27], which focusses on numerical results in terms of new dynamics and regions of bistability obtained with this theory).
In Section 2, we derive a model for guided mode lasers of general geometry with any number of real-index guides and any number of normal modes. An important aspect of this model is the introduction of the 'overlap factors', as discussed in detail in Section 3. These are calculated by integrating products of the spatial mode solutions of the Helmholtz equation over the active regions. As such, they represent a generalisation of the optical confinement factor. It is through these factors that the particular geometry of the waveguide is introduced and their effect on the laser dynamics can be quite significant, as indicated in Ref. [28] in comparison to the coupled mode model [29].
Familiarity with the overlap factors should give the necessary physical intuition into their properties and limiting behaviour that we frequently exploit in the derivation of the double-guided model in Section 4. Here, we specialise to the case of just two guides and only consider the lowest two normal modes. This model is still quite general in regards to the waveguide geometry that may be simulated, although it is particularly appropriate for the case of symmetric waveguides. In this paper, we look at two particular cases: equal slab guides and equal circular guides, both with real, stepped refractive index profiles, as shown schematically in Figure 1. The application to coupled VCSELs with circular guides is indicated in the schematic of Figure 2, omitting the Bragg mirrors, substrate, and other structural details. Note that, in Figures 1 and 2, a resonant cavity is assumed with propagation in the z-direction, i.e., normal to the plane of optical confinement. No further account is taken of the z-direction in what follows and the values of parameters appearing in the analysis are assumed to be averaged over the cavity length. For widely separated guides, we show, in Section 4.2, that the model reduces to both the SFM [5,15,30] and coupled mode model [29] in the appropriate limits.
Having established the mathematical model, we investigate the effect of varying the optical pump polarisation in each guide via numerical simulation in Section 6. A novel feature of this model is that it allows for us to examine the spatial variation of the circularly polarised components of the optical intensity and the optical ellipticity throughout the waveguide structure. Section 6.2 provides examples of this. In Section 6.3, we give some introductory examples of stability boundaries in the plane of total pump power and normalised guide separation. This illustrates how we can use this model to investigate the effect of independently varying the pump polarisation in each guide. More generally, we may also vary the overall pump power or adjust the relative sizes of each guide, thereby introducing an effective frequency detuning. Such investigations are deferred for future study.  Here, the distance between the guides is given as 2d, whilst a is used both for the half-width of a slab guide and the radius of a circular guide. Elsewhere in this work, we use n 1 = n core and n 2 = n clad for brevity. (a) Slab waveguide. (b) Circular waveguide.

Figure 2.
A three-dimensional (3D) schematic of two coupled circular waveguides encapsulating the essence of the application to a pair of vertical cavity surface emitting laser (VCSEL) cavities. Shown are the cylindrical waveguide regions incorporating the active areas. Pumping is assumed to be confined to these regions. Note that we have omitted the Bragg stack mirrors and substrate from this figure.

Summary Overview of the Paper
Briefly summarising the structure of the paper: • Section 2: derivation of the general model for the dynamics of arrays of coupled, spin-polarised lasers; • Section 3: discussion of the 'overlap factors', which incorporate the details of the geometry of the arrays; • Section 4: specialisation to a double-guided structure and demonstration that the model reduces to both the spin-flip model and coupled mode model in the appropriate limit; • Section 5: discussion of steady state solutions, giving exact and approximate expressions; • Section 6: numerical solutions, illustrating spatial profiles for the optical ellipticity and stability boundaries; and, • Section 7: conclusions.

The Optical Rate Equations
In any waveguiding structure that is defined by a spatially-dependent relative permittivity (r), we will have optical mode solutions Φ k (r) satisfying the Helmholtz equation where Ω is a reference frequency that is taken to be the average of the modal frequencies and k is the transverse mode index. The Φ k (r) are known as the normal modes or, sometimes, supermodes of the waveguide. For modes of a given order, the orthogonal polarisations have almost exactly the same spatial profile (having checked this numerically for cases of interest), and we shall assume this to be precisely true. Thus, each mode Φ k (r) may be associated with two polarisations. Later, we shall explicitly formulate this in terms of left and right-circularly polarised light. After cancelling a phase factor e iβz , where β is the propagation constant along the cavity length (see Appendix A of Ref. [31]), the total optical field may then be written as a superposition of the normal modes, as where E k (t) is a time dependent Jones vector that incorporates the polarisation and ν k is the modal frequency, determined via a solution of the Helmholtz equation for the mode. Hereafter, we shall use r = xe x = ye y for brevity wherever we need to denote spatial coordinates, but it should be remembered that r is confined to the x − y plane. Starting from the general form of Maxwell's wave equation and applying the slowly varying envelope approximation (SVEA), as described in Appendix A of Ref. [31], we may obtain a set of optical rate equations for the complex amplitudes of the normal modes Here, the ± subscripts denote the right (+) and left (−) circularly polarised components, τ p is photon lifetime, c is the speed of light, n g is the group refractive index, and α is the linewidth enhancement factor, as defined in terms of the change in the real and imaginary components of the electric susceptibility, ∆χ ± and ∆χ ± , respectively, by α = −∆χ ± /∆χ ± . Note that we have adopted this sign convention for consistency with the SFM model [5,15,30] and it is opposite to that used in Ref. [29]. Hence, the values of α that are used in this work take the opposite sign to that in the latter reference.
Each polarisation component is coupled to the other via the birefringence rate γ p and dichroism rate γ a . The time-dependent exponential factor involves the difference between modal frequencies ∆ν kk = ν k − ν k .
The summation over i in the last term of (3) is over the optically confining guides.
Here, we have defined optical overlap factors Γ (i) kk for the ith guide via ± is the average gain for each polarisation in guide (i) and the integral is over the ith guide. In practice, we take the gain to be spatially constant over a guide and zero outside it. Hence, in this paper, the overlap factors are simply defined by and we normalise the spatial profiles, so that where the integral is over all space.

The Carrier Rate Equations
The rate equations for spin-polarised populations of carriers may be derived from the optical Bloch equations. The general result for the spatially dependent carrier concentrations N ± (r, t) and circularly polarised optical fields E ± (r, t), including spin relaxation, may be found to be given by where τ N is the carrier lifetime, Λ ± is the pumping rate, and γ J is the spin relaxation rate. The ± subscripts on N refer to spin up (+) and spin-down (−) carriers, which directly couple with right (+) and left (−) circularly polarised photons respectively. Here, we assume that all carrier pumping, whether that be optical or electrical, is confined to the active region, and, hence, the effects of lateral diffusion are neglected at this time. Note that we take |E ± (r, t)| 2 = S(r, t) to be the photon density and, hence, E ± does not have dimensions of electric field in (5). This is unproblematic, since the optical rate Equations (3) may be multiplied by any arbitrary factor to match the dimensions of E ± in (5) without changing the dynamics.
The earlier assumption that the gain is spatially constant over a given guide and zero between guides requires a similar assumption for the carrier concentrations. We shall assume a linear gain model of the form g(N) = g 0 (N − N 0 ), where N 0 is the transparency concentration, so, if g(N) is a step function, N ≤ N 0 outside the active regions. With the pumping confined to the active regions and no spatial diffusion, we may take any optical loss in the cladding regions to have been absorbed into the cavity loss rate κ. Therefore, we may take the carrier concentration in this region to be exactly N 0 .
Taking the spatial dependence to be in the x − y plane only, we may then put where N (i) (t) is the time dependent carrier concentration in the ith guide and ξ (i) (x, y) is a step function. Because we may subtract N 0 from either side (which we do on normalisation), we may take this as effectively giving zero outside the active regions. Applying this assumption, we find that the rate equations for the spin-polarised concentrations in the (i)th guide are given by where the (i) superscripts on a quantity label the values of that quantity in each guide. Appendix B of Ref. provide the details of the derivation [31].
In this study, we assume that |∆ν k,k | ν k , ν k . This assumption is physically relevant, provided that the coupled waveguides are well separated (relative to a characteristic length). In that case, the modal frequency of the symmetric and anti-symmetric modes will become very similar. However, the assumption may not be so accurate for the frequency difference between different orders of transverse modes. Our main interest lies only in the symmetric and anti-symmetric versions of the lowest order mode of a two-guide structure, in which case the assumption is well-justified. Under this assumption, the non-autonomous Equations (3) and (7) may have their explicit time dependence removed and simplified to and where we have used the tilde notation to denote the transformed optical field. Equations (8) and (9) then represent the general model for any number of normal modes and any number of confining guides. In the following, instead of the model (3) and (7), we will analyse (8) and (9). Nevertheless, the analysis of the latter equations will still be valid in recognising unstable solutions of the former ones with a critical eigenvalue λ that is much larger than |∆ν k,k |. which is because, before the factor exp(i∆ν k,k t), which is slowly varying, starts to have any effect in the system, the unstable solution will already show its instability.
Stable solutions of (8) and (9) will also correspond to stable solutions of the model (3) and (7) if the time frame is of order O(1/|∆ν k,k |). In this way, we also conjecture that, if all the eigenvalues of a solution are far away from the imaginary axis, then the presence of the slowly varying phase exp(i∆ν k,k t) should not change the eigenvalues much.

Equal Guides
The overlap factors that are defined by (4) are calculated from the spatial modal solutions of the Helmholtz equation Φ(r). Details of the solutions used in this work are given in Appendix E of Ref. [31]. Firstly, we shall just consider equal guides with the same refractive index n 1 in the core regions and n 2 elsewhere, although our treatment of the rate equations in Section 4 is general enough to deal with twin guides of any geometry. Table 1 lists the parameters used in our example calculations. The guides may be characterised by a normalised decay constants u (in the core regions) and w (in the cladding regions) and where a is either the half-width of a slab waveguide or the radius of a circular guide, as illustrated in Figure 1. These decay constants can also be combined into a conventional 'normalised frequency' v, which is defined as In practice, we shall take (10)- (12) to refer to the value for a single isolated guide. In a solitary slab guide, for values of v < π/2, only one guided TE mode is supported. For equal double slab guides with this same index profile, only two TE modes are supported. Therefore, we refer to such structures as being 'weakly-guiding'.
Because of the symmetry of equal guides, the lowest mode has even parity and the second lowest, odd parity. We refer to these as the 'symmetric' and 'anti-symmetric' modes, respectively, and then label them by suffixes s and a, respectively. Note that the antisymmetric mode Φ a (x) always goes through zero in between the guides, whereas the Φ s (x) does not. This qualitative behaviour persists, even when we break the symmetry of the guides, so that we may still use s and a as labels, although they would then be distinguished by topology rather than geometric symmetry.
The overlap factors Γ (i) k k are found by integrating the products of the spatial modes over each guide, as in (4). For closely spaced guides, the products Φ 2 s (x), Φ 2 a (x) and the modulus |Φ s (x)Φ a (x)| are noticeably different (see Figures 3 and 4 of Ref. [31]). Hence, we note that, in general, Γ However, in the case of equal guides, by symmetry, we do have aa (equal guides).
Additionally, by symmetry, the integral of |Φ s (x)Φ a (x)| will be the same in each guide, although the sign will be opposite, so As the separation between the guides gets larger, the spatial profiles become like those of isolated guides. The difference between these profiles and those of an isolated guide are: (i) the normalisation-the integral of the squared modulus of the spatial modes will be half that of the optical confinement factor--and (ii) the sign of the anti-symmetric mode is flipped in one of the guides. In fact, as the separation tends to infinity, we may obtain the wave functions of the isolated guides by adding and subtracting the modes via This is the basis for the definition of the 'composite modes' in terms of the normal modes that are defined in (18) and (19) defined in the next section. We then have where Γ S is the optical confinement factor of a single guide. In this limit, we also have  . Plot of the overlap factors relative to Γ S for a weakly-guiding (v = π/2) symmetric slab structure as a function of spatial separation between the guides. Here, 2d is the edge-to-edge separation between the guides and 2a is the guide width (8 µm in this case). Only the overlap factors for guide (1) are shown. The factors for guide (2) are the same, except that Γ  . Spatial modes for slab guides with widths w 1 = 7.9 µm and w 2 = 8.1 µm (the wider guide is on the right). In this case, we define a = (w 1 + w 2 )/4. The normalised frequency for an averaged guide is v = π/2. Note that, as the guide separation increases, the s mode becomes more greatly confined to the wider guide, whilst the a mode is confined to the narrower. (a) Spatial modes for unequal guides with d/a = 1. (b) Spatial modes for unequal guides with d/a = 2.
For such a weakly guiding structure, there is significant variation of the factors for d/a < 2. For more strongly guiding structures, v > π/2, this variation from Γ (1) sa is greatly reduced.

Unequal Guides
For double-guided structures with unequal guiding regions, we lose the symmetric relations that were previously found. Figure 4 illustrates two examples with guide widths w 1 = 7.9 µm and w 2 = 8.1 µm, where the wider guide is on the right. In these cases, we have put a = (w 1 + w 2 )/4, whilst the same refractive index difference has been used as for the equal width guides. It can be clearly seen that the component of each mode is different in each guide, and it is now the case that aa (unequal guides).
We also note that each mode becomes more localised to a particular guide with increasing separation, with the s-labelled mode tending to the wider guide. This may be understood from basic waveguiding theory, since the s mode has the biggest effective index (or propagation constant) and the effective index (or propagation constant) of a mode of a single guide varies as the width. As the separation between the guides tends to infinity, each normal mode approaches the mode of a single isolated guide. Hence, labelling the narrow and wide guides (1) and (2), respectively, where Γ 1 and Γ 2 are the optical confinement factors of single guides of width w 1 and w 2 . Figure 5 illustrates these behaviours, which shows the variation of the overlap factors for a non-symmetric slab guide as a function of spatial separation. The calculations presented here use the same guide widths as for the modes that are shown in Figure 4. Note that, since Γ S is the optical confinement factor of the averaged isolated guide, the limiting values of ss /Γ S and Γ (1) aa /Γ S are not exactly the same. Figure 5. Plot of the overlap factors for a weakly-guiding (v = π/2) non-symmetric slab structure as a function of spatial separation between the guides. Here, 2d is the edge-to-edge separation between the guides and 2a = (w 1 + w 2 )/2 is the average guide width, using w 1 = 7.9 µm and w 2 = 8.1 µm, as in Figure 4. In contrast to the case of symmetric guides, in this case the s mode tends to occupy guide (2) (the wider guide) and the a mode occupies guide (1). As the separation increase, the corresponding overlap factors tend to the optical confinement factor, whilst all other factors tend to zero. Here, the modulus of Γ (2) sa is shown, as this factor is negative.

Circular Guides
For circular guides, we used a commercial eigensolver to find solutions of the Helmholtz equation for symmetric structures at various guide separations. These solutions are for the two polarisation components of the lowest-order symmetric and antisymmetric mode. To interpolate between these results, we found that the overlap factors could be fitted very well by the following empirical formulae: and Γ (i) Table 2 lists the parameters used. Here, Γ S is again the optical confinement factor associated with the lowest mode of a single isolated guide. For circular guides, this is HE11 (LP01) mode.
In Section 4.2.2, we discuss the reduction of the normal mode model to the coupled mode model. It is found that the coupling coefficient µ is given by the difference in normal mode frequencies µ = (ν s − ν a )/2. Using the calculated values of these frequencies, the coupling coefficient for equal circular guides is found to be well-approximated by the Ogawa [32] expression where w is given by (11). The constant of proportionality may be found by fitting this to the calculated value of (ν s − ν a )/2 at d/a = 1. Table 2. Parameters used in fitting functions for circular waveguides. Here, the cladding refractive index is n 2 = 3.4 and the core refractive index is n 1 = n 2 + ∆n.
Parameter ∆n = 0.000971 Unit In this paper, we confine our attention to structures that only involve two weaklyconfining guides supporting only two guided modes. The treatment that we shall follow in this section will be valid for the general case of unequal guides, although it is of particular use for the symmetric case, which we focus on in this paper.
For a double-guided structure, we may put reference frequency to Ω = (ν s + ν a )/2. The optical rate equations of (8) may now be written where k = s, a for the symmetric and anti-symmetric modes, respectively and where k = a, s. These are then the equations for the evolution of the normal modes. However, solutions in terms of the normal modes do not lend themselves well to physical intuition. When thinking of optical guides in close proximity, it is more natural to think of the optical intensity in each guide. To this end, it is convenient to introduce new optical field variables, as defined by and The motivation for this is that the squared modulus of these 'composite' modes becomes the optical intensity in each guide at infinite separation, which greatly aids in the visualisation of the laser dynamics.
In the following sections, we write all of the equations in their real form. Although these are often quite ungainly, there are practical advantages to rendering them in this way. In particular, any stability analysis of these equations must start with writing the equations out in real form anyway. Moreover, the real form is generally better suited, if not required, for numerical calculations, whether this be for dynamical simulation via the Runge-Kutta method or finding the steady state solutions using a nonlinear solver.
Defining the φ 21±± as the phase difference between E 2± and E 1± and φ kk+− as the phase difference between E k+ and E k− , in Appendix C of Ref. [31] we find that the optical rate equations may be written in real form as and Note that, since φ 22+− = φ 21++ − φ 21−− + φ 11+− , we do not need an equation for this last phase variable.
In these equations, we have used whilst the gain terms are defined by and The Γ terms that are introduced above are further defined in terms of the optical overlap factors via and It is worth noting a general limiting behaviour as the separation between guides tends to infinity that Γ is half the optical confinement factor in an isolated guide. Hence, in this limit (recalling that the separation between the guides is 2d), For equal guides, 2Γ (i) ∞ = Γ S , the optical confinement factor of an isolated guide. For unequal guides, we may take Γ S to be an average of the optical confinement factors. This does not undermine the generality of (20)-(23), since, in all cases, the factor of Γ S in the denominator of the gain terms cancels with the factor multiplying it. The inclusion of Γ S here is one of convenience to elucidate the limiting behaviour of the rate equations, as discussed later in Section 4.2.

Carrier Rate Equations
The carrier rate equations are straight-forward to render in the notation for the double guided structure. From (9), these are where I (i) Using Γ as , we can expand (32) in terms of E 1,± and E 2,± , as Expressing this in terms of (28) and (29), we have

Normalised Rate Equations
Upon normalisation (described in detail in Appendix D of Ref. [31]), Equation (31) may be re-written as where the M (i) ± are the normalised carrier densities and Here, we have used assumed a linear gain model, having put where g 0 is the differential gain and N 0 is the transparency concentration. Hence, in the steady state, we find the threshold carrier concentration N S for the single guide to be given by The normalised carrier concentrations appearing in (33) are then of the general form Hence, M will be zero at transparency and unity at threshold.

Further defining
and and the cavity decay rate κ = 1/(2τ p ), the normalised optical rate equations are and Here, the amplitudes A are normalised according to Equations (33), (34), and (40)-(43) constitute the dynamical model of double-guided structure. This may be applied to any waveguide geometry that is restricted to two guides and the two lowest confined modes.

Reduction to the Spin-Flip Model (SFM)
In this section, we assume equal guides and, so, employ the results for the identities and limiting behaviour of the overlap factors that are found in Section 3.1. Because d → ∞, the factors that are defined earlier in (28) and (29) become The carrier terms that are defined in (37)-(39) then become ± and ∆M ± = 0 (44) and the expression for the optical intensity of (34) appearing in the carrier rate equations becomes I (i) for i = 1, 2. Additionally, note that the term given earlier in (24) as µ = (ν s − ν a )/2 approaches zero as the guide separation approaches infinity and the frequencies of the normal modes become equal. Hence, Equations (40)-(43) for the normalised optical rate equations reduce to for i = 1, 2, Meanwhile, (33) for the carrier rate equations becomes The guides are now completely uncoupled and we have two sets of equivalent equations for each. Defining new variables N = (M (i) for each guide, we may re-write (46)-(49) as and Here, η = (η + + η − )/2 (dropping the (i) superscripts), we have defined an effective spin relaxation rate γ s = γ + 2γ J and P is the pump ellipticity defined by Equations (50)-(53) are the real form of the spin-flip model (SFM) [5,15,30]. Note that Gahl et al. [15] have the factor of (1 + iα) multiplying the cavity loss term in the complex rate equations, which we take to be unphysical, so there will be a discrepancy between their expressions and those above.

Reduction to the Coupled Mode Model (CMM)
We may consider an alternative scenario in which we retain the coupling term µ between the guides, but let the overlap factors take their limiting values as the guide separation tends to infinity. Additionally, we may remove the coupling between the spin polarised components by setting γ J = γ a = γ p = 0. In this case, Equations (40)-(43) reduce to and whilst the carrier rate equations become These give us two independent sets of equations for the polarisation components, each of which reproduces the coupled mode model of Ref. [29] with real coupling coefficient µ and no frequency detuning (although with a difference in sign on the α factor, due to opposite sign definitions). It is of particular note that the coupling coefficient µ, given in terms of the difference between the mode frequencies in (24), is consistent with the analysis of Marom et al. [33], who found the same relation between the coupled mode and normal mode models.

Steady State Solutions
In general, analytical steady state solutions of the double-guided structure are not obtainable. However, exact expressions and very good approximations are both available in certain limiting cases. In this section, we continue to consider only the case of symmetric guides. Only results are given here; for additional details of the derivations, see Ref. [31].

Effect of Spin Relaxation
In the steady state, the carrier rate equations of (33) yield Alternatively, we may make the intensities the subject, giving Close to threshold, we may take the optical intensities in (60) to be negligible and put these to zero, giving Note that, in general, the different polarisation components of the intensity will not go to zero at the same overall pumping rate η (i) = η (i) (62) is not exact. However, with a large spin relaxation rate γ J γ, Equation (62) is found to be a good approximation (in practice, γ J > 10γ).

Equal Pumping
A useful simplification to make is to assume equal pumping in each guide. That is, the total pump power and pump ellipticity are both the same η ± . Under the equal pumping assumption, Equation (42) reduces to In the steady state, this is satisfied for φ 21±± = 0, π. On reduction to the coupled mode model [29], these are referred to as the 'in-phase' and 'out-of-phase' solutions, respectively.
Because |A 1,± | = |A 2,± |, Equations (40) and (41) may be written as where n = 0, 1 for the in-phase and out-of-phase solutions, respectively, and i = 1, 2 for each guide. We may simplify (64) even further by assuming that there is no direct optical coupling between the polarisation components. That is, the dichroism and birefringence rates are both zero, γ a = γ p = 0 (note that these components may still be indirectly coupled via the spin-polarised carrier concentrations).
With this simplification, we find for the in-phase solution and for the out-of-phase solution (at infinite separation, Γ and With these results for |A i,± |, M ± and φ 21±± , we see that nothing depends on φ 11+− (moreover, it may also be shown that ∂φ ii+− /∂t = 0 follows without assuming it to be so). Hence, in this simplified case, we have found exact, analytic steady state solutions for all of the variables (whilst φ 11+− may be set arbitrarily).
If we now allow for γ a and γ p to be non-zero, we may obtain the condition where we have defined the modal optical ellipticity in the (i)th guide via Equation (69) agrees with the result of taking the ratio of Equations (15) and (16) in Adams et al. [18], in the case where there is no misalignment of birefringence and dichroism.
We describe (70) as the 'modal' ellipticity, since it is terms of the composite mode amplitudes. Although this is defined for each guide, there is spatial dependence beyond this. Later, in Section 6.2, we shall define a spatially varying ellipticity, hence the reason for the specific nomenclature here. Note that, since tan(φ ii+− ) = tan(mπ + φ ii+− ), where m is an integer, we have two possible solutions for φ ii+− of φ 0 and φ 0 + π, where φ 0 is the solution of the arctangent of Equation (69) for φ 0 ∈ [−π/2, π/2]. Applying this consideration, we find where m = 0, 1. In the reduction to the spin-flip model, these solutions are also referred to as being 'in-phase' (m = 0) and 'out-of-phase' (m = 1). However, we will not use this terminology here, in order to avoid confusion with the previously defined meaning of these terms in the context of the coupled mode model. In this case, no closed form expression for the optical amplitudes can be found. However, the results of numerical simulation show that using (67) and (68) in conjunction with (69) provides a very good approximation for the pump ellipticities of |P (i) | <∼ 0.8.

Results and Discussion
In this paper, we mainly focus on the role of pump ellipticity. Further investigations into bistability and polarisation switching based on the theoretical model that was developed here were carried in [27], which may be viewed in conjunction with the present paper.

Numerical Solutions of the Rate Equations
For general solutions of the model, we employ a combination of numerical methods. For time series simulations of (33) and (40)-(43), we use an adaptive Runge-Kutta method of orders 4 and 5 [34,35]. This is very useful for finding both stable steady state solutions and simulating the temporal dynamics in regions of instability. For finding the unstable steady state solutions and computing the Jacobian, we use a nonlinear solver implementing a trust-region dogleg algorithm [36] based on the interior-reflective Newton method [37,38]. Where steady state solutions exist, this latter method is much faster than time series simulation, facilitating efficient routines for tracing stability boundaries and analysing the Jacobian eigenvalues to establish the nature of bifurcations.
In this section, we give the results for weakly guided structures with v = π/2. Table 1 lists the waveguide parameters, whilst Table 3 depicts the specific laser parameters used.

Spatial Profiles
The spatial profiles of the intensity and output ellipticity presented in this section were calculated directly from the numerical solutions of the Helmholtz equation for circular guides (i.e., no empirical interpolation was used). The spatially-dependent output ellipticity is defined as where the circularly-polarised intensities are given in terms of the normal mode amplitudes A k,± , by The normal mode amplitudes A k,± are reconstructed from the composite mode solutions, as described in Appendix C of Ref. [31]. In Figure 6, we show the intensity results with pump ellipticities of P (1) = 0 and P (2) = 1 and a total normalised pump power in each guide of η (i) = η (i) + + η (i) − = 100 for circular guides of radius a = 4 µm and an edge-to-edge separation that is given by d/a = 1 with an optical wavelength λ 0 = 1.3 µm.
Here, the mid-line joining both of the guides is in the y-direction and guide (2) is in the positive y half of the plane (to the left in the diagrams). We note a residual component of left-circularly-polarised light in guide (2). At this pumping power, this is largely due to spatial coupling between the guides. Figure 7 shows the corresponding output ellipticity for these intensities. We note a dip in between the guides where the ellipticity goes to −1, even though the pump polarisation is P (2) = 1 in guide (2) and P (1) = 0 in guide (1). The reason for this can be seen in Figure 8, which plots the modal amplitudes against the pump ellipticity in guide (2). Across the whole range, the optical polarisation is dominated by the anti-symmetric mode, which follows the ellipticity of P (2) . However, at y = 0, the anti-symmetric mode goes through zero, so the only contribution at this point comes from the much smaller symmetric component, for which the left-circularly polarised amplitude is slightly larger. Hence we see this dramatic dip. However, note that the optical intensity for both components is very small in this region.

Stability Boundaries
An initial comparison of the normal mode model to the coupled mode model [29] (neglecting polarisation) has highlighted the importance of the overlap factors on the laser dynamics [28]. This is more significant for the weakly-guided structures that are designed to support only the lowest optical modes. Here we extend this initial investigation to explore the effects of including both optical and carrier-spin polarisation. In particular, we consider the stability of the laser dynamics as a function of the ratio of the optical pump power to the threshold pump against the normalised guide separation. In all cases, we take the total pump power in both guides to be equal, which is η (1) = η (1) − , and so we may drop the guide index.
Following Ref. [29] in the case of zero pump ellipticity, we may define the pump-topump threshold ratio in terms of Q = (η + + η − )/2 (the factor of two accounts for the normalisation method used) and consider the threshold pump power Λ th at infinite guide separation. This gives us the relation where C Q = g 0 N 0 /g th and the threshold gain is g th = 2κn g /(cΓ S ). The values of the differential gain g 0 , transparency density N 0 and other laser parameters are given in Table 3. With the optical confinement factor Γ S for a circular guide with refractive index step ∆n = 0.000971 listed in Table 2, this gives a value of C Q = 43.9. For non-zero pump polarisation, we may expect the threshold pump to be affected. We retain the definition Q = (η + + η − )/2 and, again, consider the behaviours at infinite separation. With no pump polarisation, the threshold pump would be Q = 1. In the same limiting conditions, our model is equivalent to the spin-flip model.
Following the same analysis that led to the steady state threshold density expression (71) and noting that, when the laser just turns on Q th = N th = (M (i) where ε is the modal ellipticity in either guide, as given by (70) and φ 0 is the solution of (69) for φ 0 = φ ii+− ∈ [−π/2, π/2]. Note that we have taken the limiting condition Γ S /(2Γ (i) kk ) → 1. Using the parameters shown in Table 3, we find, via numerical simulation, that, even for P (i) = 1, the modal ellipticity does not exceed 0.6, and Q th − 1 is of order 0.001. Hence, we may safely neglect the effect of the pump ellipticity on (74).
In the Λ/Λ th − d/a plane, we find a Hopf bifurcation separating the stable steady state solutions in the upper right half of the plane from unstable solutions in the lower left. Regions of stability also appear at low pump powers and small separation, again being separated by a Hopf bifurcation. In both of the regions, these stable solutions are the out-of-phase solutions, which have predominantly anti-symmetric modal components.
Solutions for slab waveguides neglecting the polarisation have been initially reported in Ref. [28] and compared to the coupled mode model results. It was found that, generally, the overlap factors tended to push the boundaries up in the direction of both increasing power and increasing separation, thus enlarging the regions of instability.
In this work, we focus on the effect of optical pump ellipticity and the simulation results described here are for weakly-guiding (∆n = 0.000971, v = π/2) circular guides. Starting with equal pump ellipticities P (1) = P (2) in Figure 9. Here, we have plotted curves from P (1) = 0 to P (1) = 1 in steps of 0.2. The instability region is steadily reduced as P (1) = P (2) is increased. On the other hand, the boundary of the lower stability region remains insensitive to these changes.
In Figure 10, the ellipticity in guide (1) is kept fixed at P (1) = 0, whilst P (2) is varied from 0 to 1 in steps of 0.2. Here, we see the same qualitative behaviour, with the stability moving towards the origin with increasing P (2) , although not to the same extent as in Figure 9. Note that the lower stability region has not been plotted in this case. Figure 11 effectively shows the continuation of this set of results (but starting with P (1) = 1 and P (2) = 0) and increasing P (2) in steps of 0.2 to 1, ending with P (1) = P (2) = 1, as in Figure 9. Finally, Figure 12 shows the effect of putting P (1) = −P (2) and increasing P (1) from 0 to 1. In this case, the final curve with P (1) = 1 and P (2) = −1 does not diminish the instability region to quite the same extent as P (1) = 1 and P (2) = 1, although we still see the same qualitative reduction of the region with increased pump ellipticity.

Conclusions
We have derived a set of general rate equations for the laser dynamics in a waveguide of arbitrary geometry supporting any number of guiding regions and normal modes.
The details of the geometry are encoded into the 'overlap factors' (generalisations of the optical confinement factor), which may then have a significant effect on the laser dynamics.
We have focused on the particular case of a double-guided structure in the case of just two supported normal modes and derived a set of real rate equations in terms of 'composite modes'. This treatment is particularly useful for the consideration of symmetric guides (with reflection symmetry) and it can be shown to reduce to both the spin flip model and the coupled mode model in the appropriate limiting case.
Assuming symmetric guides, we have found both exact and approximate analytical expressions for steady state solutions in certain simplified cases. In particular, we have looked at the case of equal power pumping in both guides and investigated the spatial solutions for the circularly polarised intensities (and optical ellipticity) and stability maps for different pump ellipticities. We have found that, in general, increasing the pump ellipticity reduces the region of instability in the pump power versus the normalised separation plane.
The effect of high birefringence has been explored in subsequent work [27] based on the model that was developed here, revealing new dynamics and regions of bistability. It has been shown that optical switching of the polarisation states of the lasers may be controlled through the optical pump and that, under certain conditions, the polarisation of one laser may be switched by controlling the intensity and polarisation in the other.
Although the equations that we have derived are general enough to deal with nonsymmetric guides, where there is a large difference between guides in double guided structure, it may be more appropriate to use the general solutions, since the behaviour of the overlap factors tends to diverge from the symmetric case quite rapidly. We have not investigated the resulting dynamics that are associated with this in this paper and leave this matter to be addressed in more detail in future work.