*Article* **Synchronization of Mutually Delay-Coupled Quantum Cascade Lasers with Distinct Pump Strengths**

#### **Thomas Erneux 1,\*,† and Daan Lenstra 2,†**


Received: 25 October 2019; Accepted: 5 December 2019; Published: 10 December 2019

**Abstract:** The rate equations for two delay-coupled quantum cascade lasers are investigated analytically in the limit of weak coupling and small frequency detuning. We mathematically derive two coupled Adler delay differential equations for the phases of the two electrical fields and show that these equations are no longer valid if the ratio of the two pump parameters is below a critical power of the coupling constant. We analyze this particular case and derive new equations for a single optically injected laser where the delay is no longer present in the arguments of the dependent variables. Our analysis is motivated by the observations of Bogris et al. (IEEE J. Sel. Top. Quant. El. 23, 1500107 (2017)), who found better sensing performance using two coupled quantum cascade lasers when one laser was operating close to the threshold.

**Keywords:** two delay-coupled lasers; weak coupling limit; optically injected laser

#### **1. Introduction**

Compact quantum cascade lasers (QCLs) emitting in the midwave infrared (mid-IR) are the leading semiconductor laser sources for such applications as absorption spectroscopy in the molecular fingerprint region [1,2]. Mid-IR spectroscopy has led to new applications in biology and medicine such as breath analysis, the investigation of serum, noninvasive glucose monitoring in bulk tissue, and the combination of spectroscopy and microscopy of tissue thin sections for rapid histopathology [3]. Other applications include environmental sensing and pollution monitoring, industrial process control, and security [4,5].

Recently, a novel gas sensor relying on a pair of mutually injecting QCLs has been analyzed both experimentally and numerically [6–8]. The sensing performances of the coupled QCLs have been examined in terms of the injection power, bias currents of the lasers, and their spectral detuning. High sensitivity is observed if one of the two lasers is biased around the threshold. The main objective of this paper is to explain these observations by analyzing the rate equations appropriate for two coupled QCLs. As we shall demonstrate, allowing one laser to operate close to its threshold contributes to larger domains of stable phase locked states. Physically, the transient response of the intensity of one laser slows down near its threshold, while the intensity of the second laser is keeping its fast time scale. Consequently, the fast laser quickly approaches a quasi-steady state regime, and the long time dynamics of the laser system is controlled by the slow laser. In other words, the coupled QCLs is becoming an injected laser problem where the fast and slow lasers are acting as master and slave, respectively. In a different setting, two coupled QCL cavities separated by a gap of 3 μm were studied

as a monolithic integrated photodetector [9]. An integrated detector is used for spatial sensing of the light intensity, and its control is again achieved by changing the applied bias.

The dynamics of two mutually delay-coupled semiconductor lasers (SLs), in a face-to-face configuration, has been a topic of active research [10–14]. The time delay results from the finite propagation time of the light from one laser to the other one. Of primary interest are the conditions for stable locking, and systematic studies have been undertaken in order to explore the effects of key parameters. These led to striking comparisons between experimental and numerical bifurcation diagrams [15–21]. Most often, the time delay is relatively large compared to the photon lifetime (21 to a 51 mm gap between the lasers) [16–18]. However, systems of two coupled lasers in photonic integrated circuits have recently been investigated (1–2 mm gap) [11,12]. They revive previous theoretical investigations of the short coupling regime [22].

In a different optical setting, two laterally coupled semiconductor lasers (no delay) also raised the interest of researchers for the presence of exceptional points (EPs) in parameter space [23–25]. An EP is a point where two (or more) eigenvalues simultaneously coalesce. One key difference between EPs and conventional degeneracies is their higher sensitivity to perturbations. This particular property of EPs has been proposed for use in sensor applications [26,27].

QCLs, based on intersubband transitions in semiconductor quantum wells, are characterized by ultrafast (picosecond) carrier lifetimes. An important consequence of this unique property is the absence of relaxation oscillations (RO) in the transient response of these devices. For conventional interband diode lasers (IDLs), the ROs are generating undesirable intensity oscillations for quite low feedback amplitudes. By contrast, dynamical instabilities for QCLs are only possible if the delayed feedback is strong enough [28,29].

For two coupled lasers operating at close, but distinct optical frequencies, the desired regime is when the lasers operate in a continuous wave (CW) with their frequency and phase mutually locked. They are called one color states [30] or compound laser modes (CLMs) [19,22]. To the best of our knowledge, phase locked states of two delay-coupled SLs have been investigated theoretically with equal or nearly equal pump parameters. However, the individual laser pump rates are experimentally controlled variables, and the effects of unequal pumps have been studied for two SLs without delays [23–27,31].

The organization of the paper is as follows. Section 2 introduces the rate equations for two coupled QCLs, as well as their asymptotic approximation, valid in the limit of weak coupling, weak frequency detuning, and arbitrary pump parameters. It consists of two delay coupled Adler equations for the phases of the fields. The CLMs are then investigated in Section 3 in terms of their frequencies. As functions of the detuning, these frequencies appear as close orbits in the bifurcation diagram. As the pump parameter of one laser comes close to threshold (*P*<sup>2</sup> → 0), these orbits overlap progressively larger domains of detuning. In Section 4, the limit *P*2/*P*<sup>1</sup> → 0 is analyzed in detail taking into account that our problem now depends on two small parameters, namely *P*2/*P*<sup>1</sup> and the small coupling rate. A new asymptotic analysis of the original laser equations is performed and leads to equations for an optically injected single mode laser where the delay no longer appears in the arguments of the dependent variables. The stability of the locked states is then analyzed. If one laser is operating slightly below the threshold, as in the experiments in [7], the locked state is always stable. Last, we discuss in Section 5 the impact of our results for conventional IDL lasers.

#### **2. Dimensionless Equations**

The response of a QCL subject to a delayed feedback is analyzed using rate equations formulated in [32–34] on the basis of a three level model. In [28], it was shown that these equations for a QCL subject to delayed optical feedback can be reduced to the classical Lang and Kobayashi (LK) equations derived for IDLs. The LK equations consist of the rate equations for a conventional SL supplemented by a term describing the optical feedback of the electrical field. Two key parameters control the dynamical stability of the laser, namely the ratio of the carrier to photon lifetimes *T* and the linewidth enhancement factor *α*. For a QCL, *T* is typically an *O*(1) quantity compared to the large *O*(103) value of an IDL. Moreover, *α* = 0 − 1 for a QCL is relatively small compared to the IDL *α* = 2.5 − 3.5. These two essential properties of a QCL explain the observed high tolerance with respect to optical feedback. Mathematically, we therefore consider two coupled LK equations as the rate equations for two QCLs coupled face-to-face. We use the formulation detailed in [19]. Specifically, the evolution equations are in terms of the optical fields *Ejopt* = *Rj* exp(*iφ<sup>j</sup>* + *iωjt*) where *ω<sup>j</sup>* is the optical angular frequency of laser *j* and the carrier densities *Nj* (*j* = 1, 2). Introducing the frequency detuning Δ = *ω*<sup>2</sup> − *ω*<sup>1</sup> and the averaged frequency *ω* = (*ω*<sup>1</sup> + *ω*2)/2, *φ*<sup>1</sup> = Δ*t*/2 + Φ1, and *φ*<sup>2</sup> = −Δ*t*/2 + Φ2, it is mathematically convenient to reformulate the two fields as:

$$E\_{1opt} \quad = \, \, R\_1 \exp(i\frac{\Delta t}{2} + i\Phi\_1 + i\omega\_1 t) = R\_1 \exp(i\Phi\_1 + i\overline{\omega}t),\tag{1}$$

$$E\_{2opt} = -R\_2 \exp(-i\frac{\Delta t}{2} + i\Phi\_2 + i\omega\_2 t) = R\_2 \exp(i\Phi\_2 + i\overline{\omega}t). \tag{2}$$

The rate equations for the amplitudes *Rj*, phases Φ*j*, and densities *Nj* are then given by [31]

$$R\_1^{\prime} = \underset{\pi \quad \text{ $N\_1 R\_1 + \varepsilon R\_2(t - \tau)$ cos}{\text{ $R\_1(t - \pi)$ cos}}}{\text{ $N\_2(t - \pi)$ cos}} (\theta + \Phi\_2(t - \tau) - \Phi\_1 - \mathcal{C}),\tag{3}$$

$$\Phi\_1^{\prime} = \begin{array}{c} -\frac{\Delta}{2} + aN\_1 + \varepsilon \frac{R\_2(t-\tau)}{R\_1} \sin(\theta + \Phi\_2(t-\tau) - \Phi\_1 - \mathbb{C}), \tag{4}$$

$$TN\_1^{\prime} = \,^1P\_1 - N\_1 - (1 + 2N\_1)R\_{1\prime}^2 \tag{5}$$

$$R\_2^{'} = \underset{\text{A}}{N\_2 R\_2 + \varepsilon R\_1 (t - \tau)} \cos(\theta + \Phi\_1 (t - \tau) - \Phi\_2 - \mathbb{C}),\tag{6}$$

$$\Phi\_2^{\prime} = \frac{\Delta}{2} + aN\_2 + \varepsilon \frac{R\_1(t-\tau)}{R\_2} \sin(\theta + \Phi\_1(t-\tau) - \Phi\_2 - \mathbb{C}),\tag{7}$$

$$T N\_2^{\prime} = \,^1P\_2 - N\_2 - (1 + 2N\_2) R\_2^2. \tag{8}$$

In these equations, time *<sup>t</sup>* is measured in units of the photon lifetime *<sup>τ</sup><sup>p</sup>* ∼ <sup>10</sup>−11*s*. Prime means differentiation with respect to *t*. *P*<sup>1</sup> = *O*(1) and *P*<sup>2</sup> = *O*(1) are the pump parameters measuring the amount of electrical current used to activate the individual lasers. The complex mutual coupling is accounted for by *ε* exp(*iθ*). *τ* and *C* ≡ *ωτ* represent the delay time and the (mean) induced phase, respectively. The distance *L* between the lasers is a few centimeters, which then implies that the delay *τ* ≡ (*L*/*c*)/*τp*, where *c* is the speed of light, is around 10. From Equation (3) with *R*<sup>1</sup> = *O*(1) and *R*<sup>2</sup> = *O*(1), we note that *N*<sup>1</sup> needs to be an *O*(*ε*) small quantity in order to balance the first two terms in the right hand side of Equation (3). Similarly, balancing all three terms in the right hand side of Equation (4) requires that the detuning |Δ| is small like *ε*. The same conclusions apply for Equations (6) and (7).

We analyze Equations (3)–(8) assuming weak coupling (*ε* << 1) and small detuning (Δ = *O*(*ε*)). The analysis leads to the following two coupled Adler delay differential equations for the phase of the electrical fields (see Section 1 of the Supplementary File):

$$\frac{d\Phi\_1}{dt} = -\frac{\Delta}{2} + \varepsilon \sqrt{\frac{P\_2}{P\_1}(1+a^2)} \sin(\theta\_0 + \Phi\_2(t-\tau) - \Phi\_1),\tag{9}$$

$$\frac{d\Phi\_2}{dt} = -\frac{\Delta}{2} + \varepsilon \sqrt{\frac{P\_1}{P\_2}(1+a^2)} \sin(\theta\_0 + \Phi\_1(t-\tau) - \Phi\_2) \tag{10}$$

where

$$
\theta\_0 \equiv \theta - \mathbb{C} - \arctan(\mathfrak{a}).\tag{11}
$$

These equations were formulated in [8], using the theory developed in [31] for the zero delay case. They were also the starting point of the investigations in [35]. We mathematically rederived those equations in a more systematic way by using an asymptotic method where Δ is scaled with respect to *ε*. This analysis is necessary as we later consider the case of one laser operating close to its threshold for which Equations (9) and (10) are no longer valid.

Equations (9) and (10) with *P*<sup>1</sup> = *P*<sup>2</sup> and *θ*<sup>0</sup> = 0 have been studied in detail [15,36,37]. Particular attention has been devoted to (1) constant phase solutions (Φ1, Φ2)=(0, *σ*) [36], (2) compound laser modes (Φ1, Φ2)=(*ωt*, *ωt* + *σ*) [19,36], and (3) time-periodic unbounded solutions (Φ1, Φ2)=(*C*/2 + Φ(*t*), −*C*/2 + Φ(*t*)) with < *d*Φ/*dt* >= *cst* [15].

It is worthwhile to briefly review the case of zero delay, which was analyzed in detail [38–40] since the pioneering work of Winful and Wang [41], who considered the case of zero detuning (Δ = 0), equal pumps (*P*<sup>1</sup> = *P*2), coupling phase *θ* = *π*/2, and -*α* replacing *α* in Equations (4) and (7). If *τ* = 0, Equations (9) and (10) can be combined into a single equation for *σ* ≡ Φ<sup>2</sup> − Φ<sup>1</sup> given by:

$$\frac{d\sigma}{dt} = \Lambda + \varepsilon \sqrt{1 + a^2} \left[ \sqrt{\frac{P\_1}{P\_2}} \sin(\theta\_0 - \sigma) - \sqrt{\frac{P\_2}{P\_1}} \sin(\theta\_0 + \sigma) \right]. \tag{12}$$

The steady states are the phase locked states. They are shown in Figure 1 for equal and non-equal pump values. We observe that the size of the locking domain increases as one of the pumps comes closer to its threshold, a feature for which we again see if the delay is not zero.

**Figure 1.** Frequency-locked states in the case of zero delay. *σ* = Φ<sup>2</sup> − Φ<sup>1</sup> is shown as a function of Δ (−*π* < Φ < *π*). Full and broken lines correspond to stable and unstable branches. *θ*<sup>0</sup> = *π*/4, *ε* = 0.02, *α* = 1, *P*<sup>1</sup> = 1, and the value of *P*<sup>2</sup> is indicated in the figure.

#### **3. Compound Laser Modes**

The compound laser modes or CLMs are the basic solutions of our problem. They are the solutions of Equations (9) and (10) of the form:

$$
\Phi\_1 = \omega t,\ \Phi\_2 = \omega t + \sigma.\tag{13}
$$

From Equations (1) and (2), we understand that after coupling, the optical frequency *ωop* is given by:

$$
\omega\_{0p} = \overline{\omega} + \omega.\tag{14}
$$

Inserting Equation (13) into Equations (9) and (10) leads to two equations for *ω* and *σ* given by:

$$
\omega \quad = \quad -\frac{\Delta}{2} + \varepsilon \sqrt{\frac{P\_2}{P\_1}(1+a^2)} \sin(\theta\_0 - \omega \tau + \sigma),
\tag{15}
$$

$$
\omega\_{\parallel} = \frac{\Delta}{2} + \varepsilon \sqrt{\frac{P\_1}{P\_2} (1 + a^2)} \sin(\theta\_0 - \omega \tau - \sigma). \tag{16}
$$

Expanding the trigonometric functions, Equations (15) and (16) are rewritten as:

$$\left(\omega + \frac{\Lambda}{2}\right) \frac{1}{\varepsilon \sqrt{\frac{P\_0}{P\_1} (1 + a^2)}} = \begin{pmatrix} \sin(\theta\_0 - \omega \tau) \cos(\sigma) \\ + \cos(\theta\_0 - \omega \tau) \sin(\sigma) \end{pmatrix} \tag{17}$$

$$\left(\omega - \frac{\Delta}{2}\right) \frac{1}{\varepsilon \sqrt{\frac{P\_0}{P\_2} (1 + a^2)}} = \begin{pmatrix} \sin(\theta\_0 - \omega \tau) \cos(\sigma) \\ -\cos(\theta\_0 - \omega \tau) \sin(\sigma) \end{pmatrix}.\tag{18}$$

From Equations (17) and (18), we determine cos(*σ*) and sin(*σ*):

$$\cos(\sigma) = \frac{1}{2\sin(\theta\_0 - \omega\tau)\varepsilon\sqrt{(1+\alpha^2)}} \begin{bmatrix} \left(\omega + \frac{\Lambda}{2}\right)\sqrt{\frac{P\_1}{P\_2}}\\+ \left(\omega - \frac{\Lambda}{2}\right)\sqrt{\frac{P\_2}{P\_1}}\\- \omega \end{bmatrix},\tag{19}$$

$$\sin(\sigma) \quad = \frac{1}{2\cos(\theta\_0 - \omega\tau)\varepsilon\sqrt{(1+a^2)}} \begin{bmatrix} \left(\omega + \frac{\Lambda}{2}\right)\sqrt{\frac{P\_1}{P\_2}}\\ -\left(\omega - \frac{\Lambda}{2}\right)\sqrt{\frac{P\_2}{P\_1}} \end{bmatrix}.\tag{20}$$

#### *3.1. Equal Pumps*

Before we consider the effect of unequal pump parameters, it is worthwhile to first analyze the case of equal pumps. The expression Equations (19) and (20) are considerably simplified as:

$$\cos(\sigma) \quad = \frac{\omega}{\sin(\theta\_0 - \omega \tau) \varepsilon \sqrt{(1 + \mathfrak{a}^2)}},\tag{21}$$

$$\sin(\sigma) \quad = \frac{\Delta}{2\cos(\theta\_0 - \omega\tau)\varepsilon\sqrt{(1+a^2)}}\tag{22}$$

and provide a solution in parametric form. We first extract *σ* = *σ*(*ω*) from Equation (21):

$$\sigma = \arccos\left(\frac{\omega}{\sin(\theta\_0 - \omega \tau)\varepsilon\sqrt{(1+a^2)}}\right) \tag{23}$$

and then compute Δ = Δ(*ω*) using Equation (22):

$$
\Delta = 2\cos(\theta\_0 - \omega\tau)\varepsilon\sqrt{(1+a^2)}\sin(\sigma). \tag{24}
$$

If *τ* = 0, the expression Equation (24) tells us that the locking domain verifies the inequality:

$$|\Delta| \le 2\varepsilon \sqrt{(1+a^2)} \cos(\theta\_0). \tag{25}$$

The expression Equation (25) is in agreement with Equation (36) in [31]. Figure 2 represents *ωτ* and *σ* as functions of Δ*τ*. The values of the dimensionless parameters are based on the following values of the original parameters for the photon lifetime *τph*, the delay *τe*, the feedback rate *εe*, and the detuning <sup>|</sup>Δ*e*<sup>|</sup> : *<sup>τ</sup>ph* <sup>=</sup> 7.5 <sup>×</sup> <sup>10</sup>−12*s*, *<sup>τ</sup><sup>e</sup>* <sup>=</sup> <sup>10</sup>−10*s*, *<sup>ε</sup><sup>e</sup>* <sup>=</sup> 2.8 <sup>×</sup> 109 *<sup>s</sup>*<sup>−</sup>1, and 0 <sup>&</sup>lt; <sup>|</sup>Δ*e*<sup>|</sup> <sup>2</sup>*<sup>π</sup>* ≤ 109 *<sup>s</sup>*<sup>−</sup>1. The dimensionless parameters are then obtained as *τ* ≡ *τe*/*τph* = 13.33, *ε* ≡ *ε<sup>e</sup>* × *τph* = 0.021, and |Δ|≡|Δ*e*| × *τph* ≤ 0.047.

**Figure 2.** compound laser mode (CLM) frequencies *ωτ* and phase difference *σ* for the case *P*<sup>1</sup> = *P*2. They are determined from the parametric solution Equations (23) and (24). The fixed parameters are *θ*<sup>0</sup> = *π*/4, *τ* = 13.33, *ε* = 0.021, and *α* = 1. The extrema of *ωτ* are *ωτ*<sup>−</sup> = −0.32 and *ωτ*<sup>+</sup> = 0.22. The extrema of Δ*τ* are the limit points Δ*τ* = ±Δ*τLP* = ±0.58.The interval [−Δ*τLP*, Δ*τLP*] is the locking range, i.e., the detuning range where the lasers lock their frequencies.

Figure 3 shows *ωτ* as a function of Δ*τ* for different values of *θ*0. The different orbits are bounded by limit points located at Δ*τ* = ±Δ*τLP*. These points mark the extreme detuning values where the two coupled lasers lock to each other. Figure 4 shows Δ*τLP* > 0 as a function of *θ*<sup>0</sup> for the interval 0 ≤ *θ*<sup>0</sup> ≤ *π*. Δ*τLP* is the largest at *θ*<sup>0</sup> = 0 and *π*. It motivates examining the limit *θ*<sup>0</sup> → 0. Figure 3 with *θ*<sup>0</sup> = 0.01 suggests that the nearly flat CLM orbit is bounded by two limit points appearing close to *ωτ* = 0. Therefore, the locking condition Equation (25) evaluated at *θ*<sup>0</sup> = 0 applies for this case.

**Figure 3.** CLM frequencies for the case *P*<sup>1</sup> = *P*<sup>2</sup> and for different values of *θ*<sup>0</sup> (indicated in the figure). The fixed parameters are *τ* = 13.33, *ε* = 0.021, and *α* = 1. As we decrease *θ*<sup>0</sup> from *π*/2, the double orbits progressively change into a single orbit.

**Figure 4.** The limit point Δ*τLP* > 0 is shown as a function of *θ*0. The maximum appears at *θ*<sup>0</sup> = 0 and *<sup>π</sup>* and is given by <sup>Δ</sup>*τLP* <sup>=</sup> *ετ*<sup>√</sup> 1 + *α*2.

#### *3.2. Unequal Pumps*

We are now ready to explore the case of unequal pumps. Using Equations (19) and (20), we eliminate the trigonometric functions of *σ* and obtain the following equation for *ω*:

$$\left\{ \begin{array}{c} \frac{1}{\sin^{2}(\theta\_{0} - \omega \tau)} \left[ \begin{array}{c} \omega \left( \sqrt{\frac{P\_{1}}{P\_{2}}} + \sqrt{\frac{P\_{2}}{P\_{1}}} \right) \\ + \frac{\Delta}{2} \left( \sqrt{\frac{P\_{1}}{P\_{2}}} - \sqrt{\frac{P\_{2}}{P\_{1}}} \right) \end{array} \right]^{2} \\\ + \frac{1}{\cos^{2}(\theta\_{0} - \omega \tau)} \left[ \begin{array}{c} \omega \left( \sqrt{\frac{P\_{1}}{P\_{2}}} - \sqrt{\frac{P\_{2}}{P\_{1}}} \right) \\ + \frac{\Delta}{2} \left( \sqrt{\frac{P\_{1}}{P\_{2}}} + \sqrt{\frac{P\_{2}}{P\_{1}}} \right) \end{array} \right]^{2} \right\} = 4\epsilon^{2}(1 + a^{2}). \tag{26}$$

Equation (26) is equivalent to a quadratic equation for Δ given by:

$$
\frac{\Delta^2}{4} \left( F\_-^2 \mathcal{C}\_1 + F\_+^2 \mathcal{C}\_2 \right) + \Delta \omega \mathcal{F}\_+ \mathcal{F}\_- \left( \mathcal{C}\_1 + \mathcal{C}\_2 \right) + \omega^2 \left( F\_+^2 \mathcal{C}\_1 + F\_-^2 \mathcal{C}\_2 \right) - 4 \varepsilon^2 \left( 1 + a^2 \right) = 0 \tag{27}
$$

where:

$$F\_{\pm} \equiv \sqrt{\frac{P\_1}{P\_2}} \pm \sqrt{\frac{P\_2}{P\_1}}, \ C\_1 \equiv \frac{1}{\sin^2(\theta\_0 - \omega \tau)}, \text{ and } \mathcal{C}\_2 \equiv \frac{1}{\cos^2(\theta\_0 - \omega \tau)}. \tag{28}$$

We only need to explore the domain 0 ≤ *θ*<sup>0</sup> ≤ *π* since *C*<sup>1</sup> and *C*<sup>2</sup> remain unchanged with −*ωτ* replacing *ωτ* and *θ*<sup>1</sup> = 2*π* − *θ*<sup>0</sup> replacing *θ*0. Figure 5 illustrates the case of small values of *P*2/*P*1. The CLM orbits are now close to the line *ωτ* = −Δ*τ*/2 and are bounded by two critical values of *ωτ* = *ω*±*τ*. They delimit the domain of real solutions of the quadratic Equation (27). We note that the CLM orbit increases in size as *P*2/*P*<sup>1</sup> → 0. An analysis of the discriminant of Equation (27) allows us to determine *ω*±*τ* (see Section 2 of the Supplementary Materials File). They delimit the domain of real solutions for Δ = Δ(*ω*). Note that they are not the values of *ωτ* corresponding to the limits points ±Δ*τLP*, but are very close if *P*2/*P*<sup>1</sup> → 0. Figure 6 shows *ω*±*τ* as functions of *P*2/*P*1. In implicit form, *x* ≡ *P*2/*P*<sup>1</sup> = *x*(*ωτ*) satisfies the quadratic equation:

$$\mathbf{x}^2 + \mathbf{x}\left[-2\cos(2(\theta\_0 - \omega\tau)) - \frac{4\omega^2}{\varepsilon^2(1+a^2)}\right] + 1 = 0. \tag{29}$$

**Figure 5.** CLMs for small values of *P*2/*P*1. *P*<sup>1</sup> = 1, and the value of *P*<sup>2</sup> is indicated in the figure. The fixed parameters are *τ* = 13.33, *ε* = 0.021, *θ*<sup>0</sup> = *π*/4, and *α* = 1.

**Figure 6.** The extrema *ωτ*<sup>±</sup> as functions of *x* = *P*2/*P*<sup>1</sup> are obtained by solving the quadratic Equation (29). Their approximations as *x* → 0 are given by Equation (30) (dotted red lines). The fixed parameters are *τ* = 13.33, *ε* = 0.021, *θ*<sup>0</sup> = *π*/4, and *α* = 1. The horizontal dotted lines mark the values of *ωτ*<sup>−</sup> = −0.36 and *ωτ*<sup>+</sup> = 0.22 at *x* = 1 (*P*<sup>1</sup> = *P*2) previously documented in Figure 2.

As seen in Figure 6, <sup>|</sup>*ω*±*τ*| → <sup>∞</sup> as *<sup>x</sup>* <sup>→</sup> 0. From Equation (29) and assuming *<sup>ω</sup>*<sup>2</sup> <sup>=</sup> *<sup>O</sup>*(*x*−1), we find the limit:

$$
\omega\_{\pm} \tau \to \pm \frac{\varepsilon \tau}{2} \sqrt{(1 + a^2) \frac{P\_1}{P\_2}} \text{ as } x \to 0. \tag{30}
$$

The corresponding values of Δ*τ* are given by Δ*τ*<sup>±</sup> = −2*ω*±. Therefore, the size of the CLM orbits satisfies the inequality:

$$|\Delta| \le \varepsilon \sqrt{(1+a^2)\frac{P\_1}{P\_2}} \ (P\_2/P\_1 \to 0) \tag{31}$$

Moreover, solving the quadratic Equation (27) and then taking the limit *P*2/*P*<sup>1</sup> → 0 lead to:

$$
\Delta \to -2\omega \pm \frac{4P\_2}{P\_1} \sqrt{\frac{4}{\mathcal{C}\_1 \mathcal{C}\_2} \left(\omega\_+ - \omega\right) \left(\omega - \omega\_-\right)}.\tag{32}
$$

In summary, we have determined the CLMs and their limits as *P*2/*P*<sup>1</sup> → 0. Because of the square roots in Equation (28), the ratio *P*2/*P*<sup>1</sup> needs to be positive. Therefore, we cannot explore the case of Laser 2 slightly below the threshold, and a different asymptotic analysis is needed where the ratio *P*2/*P*<sup>1</sup> is scaled with respect to *ε*.

We did not analyze the stability of the CLMs for arbitrary values of the pump parameters. However, we know that, because of the absence of relaxation oscillations, Hopf bifurcation instabilities are possible only for large delays [28,29]. This is not the case here. In the limit *P*2/*P*<sup>1</sup> small, we note from Equations (9) and (10) that |Φ1| freely increases while Φ<sup>2</sup> satisfies a single Adler equation. In Section 3 of the Supplementary Material File, we show that a Hopf bifurcation is not possible. Branches of CLMs are either stable or unstable, and their changes of stability occur at the limit points Δ*τLP* (saddle node bifurcation points).

#### **4. The Limit of Small Ratios of the Two Pumps**

In Section 4 of the Supplementary File, we examine the limit *<sup>P</sup>*2/*P*<sup>1</sup> <sup>→</sup> <sup>0</sup><sup>+</sup> and find that our previous theory becomes invalid as soon as:

$$P\_2 = O(\varepsilon^{2/3}).\tag{33}$$

In other words, the two coupled phase Equations (9) and (10) failed to describe the correct dynamics of the mutually injected lasers if *P*<sup>2</sup> is comparable to *ε*2/3 or smaller. Section 4 of the Supplementary Material File describes a new asymptotic analysis taking into account the scaling Equation (33). We find *<sup>R</sup>*<sup>1</sup> <sup>=</sup> <sup>√</sup>*P*<sup>1</sup> and <sup>Φ</sup><sup>1</sup> <sup>=</sup> <sup>−</sup><sup>Δ</sup> <sup>2</sup> *t*, in the first approximation, while *R*<sup>2</sup> and Φ<sup>2</sup> ≡ Φ<sup>2</sup> + <sup>Δ</sup> <sup>2</sup> *t* satisfy the equations for an optically injected laser:

$$\frac{d\mathcal{R}\_2}{dt} = \left(P\_2 - R\_2^2\right)\mathcal{R}\_2 + \varepsilon\sqrt{P\_1}\cos(\theta\_1 - \overline{\Phi}\_2),\tag{34}$$

$$\frac{d\overline{\Phi}\_2}{dt} = \Delta + a(P\_2 - R\_2^2) + \frac{\varepsilon\sqrt{P\_1}}{R\_2}\sin(\theta\_1 - \overline{\Phi}\_2) \tag{35}$$

where:

$$
\theta\_1 \equiv \theta + \frac{\Delta}{2}\tau - \omega\_1 \tau. \tag{36}
$$

The delay *τ* does not explicitly appear in the arguments of the dependent variables, but its effect appears in the expression of *θ*1. Using Equations (1) and (2), the leading expressions of the optical fields are:

$$E\_{1opt} = \sqrt{P\_1} \exp(i\omega\_1 t) \text{ and } E\_{2opt} = R\_2 \exp(i\omega\_1 t + \overline{\Phi}\_2). \tag{37}$$

where *ω*<sup>1</sup> is the optical frequency of Laser 1. The expression Equation (37) clearly indicates that Laser 1 and Laser 2 are operating as master and slave lasers, respectively. Equations (34) and (35) are the equations of an optically injected Class A laser with parameter *α* [42,43].

The steady state solution for the intensity *R*<sup>2</sup> <sup>2</sup> satisfies:

$$\left[ (P\_2 - R\_2^2)^2 R\_2^2 + \left[ \Delta + a(P\_2 - R\_2^2) \right]^2 R\_2^2 = \varepsilon^2 P\_1. \tag{38}$$

From Equation (38), we extract the solution in implicit form:

$$
\Delta \pm = -\mathfrak{a} \left( P\_2 - R\_2^2 \right) \pm \sqrt{F} \tag{39}
$$

where:

$$F \equiv \frac{\varepsilon^2 P\_1}{R\_2^2} - (P\_2 - R\_2^2)^2 \ge 0. \tag{40}$$

The two branches of solution <sup>Δ</sup> <sup>=</sup> <sup>Δ</sup>±(*R*<sup>2</sup> <sup>2</sup>) are shown in Figure 7.

**Figure 7.** Steady state solution Equation (39). The broken straight line <sup>Δ</sup> <sup>=</sup> <sup>−</sup>*α*(*P*<sup>2</sup> <sup>−</sup> *<sup>R</sup>*<sup>2</sup> <sup>2</sup>) delimits the branches <sup>Δ</sup>−(*R*<sup>2</sup> <sup>2</sup>) and <sup>Δ</sup>+(*R*<sup>2</sup> <sup>2</sup>). The parameters are *ε* = 0.021, *α* = 1, *P*<sup>1</sup> = 1, and *P*<sup>2</sup> = 0.1.

From Equations (34) and (35), we determine the linearized equations for the steady state Equation (39) and obtain the following characteristic equation for the growth rate *λ*:

$$
\lambda^2 - A\lambda + B = 0\tag{41}
$$

where:

$$\begin{array}{rcl} A & \equiv & 2(\mathcal{P}\_2 - 2\mathcal{R}\_2^2), \\ \mathcal{P}\_{\mathcal{D}\_{\Delta}} & = & (\mathcal{D}\_{\Delta} - 2\mathcal{D}\_2^2)(\mathcal{D}\_{\Delta} - \mathcal{D}\_2^2) & \mathcal{P}\_{\Delta}\mathcal{D}\_{\Delta}^2 (\mathcal{A}\_{\Delta} - \mathcal{A}\_{\Delta}\mathcal{D}\_{\Delta}) \end{array} \tag{42}$$

$$\begin{array}{rcl} B & \equiv & \left( P\_2 - 3R\_2^2 \right) \left( P\_2 - R\_2^2 \right) - 2aR\_2^2 \left( \Delta + a \left( P\_2 - R\_2^2 \right) \right) \\ & + \left( \Delta + a \left( P\_2 - R\_2^2 \right) \right)^2. \end{array} \tag{43}$$

The stability conditions are thus given by:

$$B > 0 \text{ and } A < 0. \tag{44}$$

We next analyze these two conditions. Using Equation (39), we computed *<sup>d</sup>*Δ±/*dR*<sup>2</sup> <sup>2</sup> and found that:

$$B = \mp 2R\_2^2 \left[ \frac{\varepsilon^2 P\_1}{R\_2^2} - (P\_2 - R\_2^2)^2 \right]^{-1/2} \frac{d\Delta\_{\pm}}{dR\_2^2}. \tag{45}$$

The expression Equation (45) relates *B* to the slope of the steady state branches of solutions, namely *<sup>d</sup>*Δ±/*dR*<sup>2</sup> <sup>2</sup>. *<sup>B</sup>* <sup>&</sup>gt; 0 for <sup>Δ</sup> <sup>=</sup> <sup>Δ</sup><sup>−</sup> because *<sup>d</sup>*Δ−/*dR*<sup>2</sup> <sup>2</sup> > 0 (see Figure 7). On the other hand, *B* > 0 for only parts of the branch Δ = Δ<sup>+</sup> ,verifying the inequality *d*Δ+/*dR*<sup>2</sup> <sup>2</sup> < 0 (see Figure 7). The critical points for *B* = 0 correspond to saddle node bifurcation points characterized by a zero eigenvalue and a negative or positive real eigenvalue. The condition *A* < 0 requires that *R*<sup>2</sup> <sup>2</sup> > *P*2/2. The critical points *R*2 <sup>2</sup> = *P*2/2 are Hopf bifurcation points provided that *B* > 0.

Figure 8 shows typical bifurcation diagrams. Note from Equation (42) that the stability condition *A* < 0 is always satisfied if *P*<sup>2</sup> ≤ 0, meaning no Hopf bifurcation instabilities. Figure 9 illustrates this case showing a complete branch of stable steady states.

**Figure 8.** Top: Stability diagram in terms of the pump strength *P*<sup>2</sup> and detuning Δ. The domain of a stable steady states is delimited by two Hopf bifurcation lines. They verify the scaling law |Δ*H*| → ∞ as *P*<sup>2</sup> → 0. The region *c* exhibits the coexistence of three steady states. The regions denoted by *U* correspond to an unstable steady state. Bottom: Bifurcation diagram for the intensity *R*<sup>2</sup> <sup>2</sup> as a function of Δ. The parameters are *ε* = 0.021, *α* = 1, *P*<sup>1</sup> = 1, and the value of *P*<sup>2</sup> is indicated in the figure; *H* and *SN*denote Hopf bifurcation and saddle node bifurcation points, respectively.

**Figure 9.** Same values of the parameters as for Figure 8 except *P*<sup>2</sup> = 0.

#### **5. Discussion**

In summary, we presented a rigorous asymptotic derivation of two coupled Adler delay differential equations in the limit of weak coupling and low detunings. This analysis was necessary in order to evaluate their mathematical validity as the ratio *P*2/*P*<sup>1</sup> was progressively decreased. It also suggested an alternative theory when the coupled Adler equations failed to provide the correct dynamics. This was the case if *P*2/*P*<sup>1</sup> was small like *ε*2/3 or smaller, where *ε* was the coupling strength.

Particular attention was devoted to describe analytically the locking width. The latter strongly depended on both the coupling strength and delay induced phases and increased in size as *P*2/*P*<sup>1</sup> → 0. For very low values of *P*2/*P*1, a new asymptotic analysis led to the equations of an optically injected Class A laser [44]. Laser 1 and Laser 2 were acting as master and slave lasers, respectively.

The analysis developed in this paper also applies for IDLs, but requires then taking into account the relaxation oscillations exhibited by the solitary lasers. If we define the relaxation oscillation frequency as *<sup>ω</sup>* <sup>=</sup> <sup>√</sup>2*P*1/*<sup>T</sup>* where *<sup>T</sup>* <sup>∼</sup> 102–103 is the ratio of the carrier to photon lifetimes, we verify that the derivation of the two delayed Adler phase equations described in this paper remains valid provided *ω*<sup>2</sup> >> *ε*. A general theory is more complicated than for QCLs because we need to take into account different scalings between three small parameters, namely *ε*, *P*2/*P*1, and *ω*. A preliminary analysis of the limit *P*2/*P*<sup>1</sup> → 0 indicated that the coupled laser equations reduced to the rate equations for an optically injected Class B laser [44] provided that *ω* and *P*2/*P*<sup>1</sup> verified specific scalings with respect to *ε*. In future work, we plan to investigate this case in more detail.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2304-6732/6/4/125/s1. **Author Contributions:** The two authors equally contributed to the investigation of the laser problem and its mathematical analysis.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Stability Boundaries in Laterally-Coupled Pairs of Semiconductor Lasers**

#### **Martin Vaughan 1, Hadi Susanto 2, Nianqiang Li 3, Ian Henning <sup>1</sup> and Mike Adams 1,\***


Received: 30 May 2019; Accepted: 20 June 2019; Published: 25 June 2019

**Abstract:** 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 confined 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.

**Keywords:** semiconductor lasers; coupled lasers; stability; normal modes; coupled modes

#### **1. 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–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–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 η = (ν*<sup>s</sup>* − ν*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.

#### **2. 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 2*d*. 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:

$$
\sigma^2 = \left(\frac{2\pi a}{\lambda}\right)^2 \left(n\_{\rm core}\,^2 - n\_{\rm clad}\,^2\right) \text{and} \, w^2 = \left(\frac{2\pi a}{\lambda}\right)^2 \left(n\_{\rm eff}\,^2 - n\_{\rm clad}\,^2\right) \tag{1}
$$

where λ is the free-space wavelength and *ne*ff is the effective index of the single solitary guide.

**Figure 1.** (**a**) Schematic of two coupled slab waveguides. (**b**) Schematic of two coupled circular cylindrical waveguides.

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 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:

$$E(\mathbf{x}, y, t) = \sum\_{k} E\_k(t) \Phi\_k(\mathbf{x}, y) \exp(-i\nu\_k t),\tag{2}$$

where the suffix *k* denotes the symmetric (s) and antisymmetric (a) modes, ν*<sup>k</sup>* is the frequency of mode *k*, *Ek*, Φ*<sup>k</sup>* 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:

$$
\overline{E}\_k = E\_k \exp(-i\nu\_k t). \tag{3}
$$

The rate equation for *Ẽ<sup>k</sup>* is then:

$$\frac{d\mathbb{E}\_k}{dt} = \left[i\frac{\nu\_k - \nu\_{k'}}{2} - \kappa\right]\widetilde{\mathbb{E}}\_k + \frac{c}{2n\_{\mathcal{S}}}(1 - ia)\sum\_{k'} \widetilde{\mathbb{E}}\_{k'} \{\overline{\mathbb{S}}\_1 \Gamma\_{1kk'} + \overline{\mathbb{S}}\_2 \Gamma\_{2kk'}\}.\tag{4}$$

where κ is the cavity decay rate (=1/(2τ*p*) where τ*<sup>p</sup>* is the photon lifetime), *c* is the speed of light, *ng* is the group index of the cavity, α is the linewidth enhancement factor, *g*1, *g*<sup>2</sup> are the mean gains per unit length in guides 1 and 2, and Γ1*kk'*, Γ2*kk'* are overlap factors in guides 1 and 2, defined as:

$$
\Gamma\_{jkk'} = \int\_{guide;j} \Phi\_k(x,y)\Phi\_{k'}(x,y)dxdy,\tag{5}
$$

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:

$$\int\_{-\infty}^{\infty} \int\_{-\infty}^{\infty} \left| \Phi\_k(x, y) \right|^2 dx dy = 1. \tag{6}$$

The rate equation for carrier concentration *Nj* in guide *j* is:

$$\frac{dN\_j}{dt} = P\_j - \gamma N\_j - \frac{c}{n\_\mathcal{S}} \sum\_{k,k'} \widetilde{E}\_k^\* \widetilde{E}\_{k'} \overline{g}\_{j} \Gamma\_{j \text{kk'}}, \quad j = 1, 2,\\ k = s, \\ a, k' = s, \\ a\_{\prime \prime} \tag{7}$$

where *Pj* is the pumping rate in the *j*th guide and γ is the carrier recombination rate. The conventional linear relationship between gain, *gj*, and carrier concentration is assumed:

$$\overline{\mathcal{g}}\_{j} = a\_{diff} \Big( \mathcal{N}\_{j} - \mathcal{N}\_{0} \Big), \tag{8}$$

where *adi*ff is the differential gain and *N*<sup>0</sup> is the transparency concentration.

#### **3. Results**

For numerical calculations, we use a wavelength of 1.3 μm with core radius (half-width) 4 μm and refractive indices *ncore* ≈ *nclad* = 3.4. Three values of the difference (*ncore* − *nclad*) 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 other parameters appearing in the rate equations are <sup>κ</sup> <sup>=</sup> 327 ns<sup>−</sup>1, <sup>α</sup> <sup>=</sup> 2, <sup>γ</sup> <sup>=</sup> 1 ns<sup>−</sup>1, *adi*ff <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>15</sup> cm2 and *<sup>N</sup>*<sup>0</sup> <sup>=</sup> <sup>1</sup> <sup>×</sup> <sup>10</sup><sup>18</sup> cm<sup>−</sup>3.

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, η = (ν*<sup>s</sup>* − ν*a*)/2. The dashed lines are fitted to these results using the approximation due to Ogawa [29] of the form:

$$
\eta \propto \frac{1}{d^{1/2}} \exp\left(-2w\frac{d}{a}\right),
\tag{9}
$$

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(−2*wd*/*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.

**Figure 2.** Variation of coupling rate with *d*/*a* for pairs of circular cylindrical guides with three values of *v* and for a slab guide with *v* = 1.571. Symbols (circles, squares, triangles, diamonds) are calculated from normal-mode theory; dashed and dotted lines are fitted from analytic results given in the text.

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 Γ<sup>S</sup> for the lowest-order mode (LP01) of the corresponding solitary single guide, given by [32]:

$$\Gamma\_{\rm S} = 1 - \left(1 - \frac{w^2}{v^2}\right) \left(1 - \frac{K\_0^2(w)}{K\_1^2(w)}\right) \tag{10}$$

where *K*0, *K*<sup>1</sup> 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.

**Figure 3.** Overlap factors versus *d*/*a* for coupled circular cylindrical guides with *v* = 1.571.

**Figure 4.** Overlap factors versus *d*/*a* for coupled circular cylindrical guides with *v* = 2.255.

**Figure 5.** Overlap factors versus *d*/*a* for coupled circular cylindrical guides with *v* = 3.740.

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.

**Figure 6.** Stability boundaries in the plane of *P*/*Pth* versus *d*/*a* for coupled circular cylindrical guides with *v* = 1.571. Curves labelled 'inf' are obtained using the values of overlap factors in the limit of large *d*/*a*.

**Figure 7.** Stability boundaries in the plane of *P*/*Pth* versus *d*/*a* for coupled circular cylindrical guides with *v* = 2.255.

**Figure 8.** Stability boundaries in the plane of *P*/*Pth* versus *d*/*a* for coupled circular cylindrical guides with *v* = 3.740.

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]:

$$\frac{P}{P\_{\rm th}} = \frac{2\alpha \frac{\eta}{\mathcal{V}} - \frac{a\_{diff} N\_0}{\overline{\mathcal{S}}\_{\rm th}}}{1 + \frac{a\_{diff} N\_0}{\overline{\mathcal{S}}\_{\rm th}}}.\tag{11}$$

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].

#### **4. 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.

**Author Contributions:** M.A. conceived the project and the general flow of the theoretical work. M.V. carried out the theoretical derivation of the equations. M.V. and I.H. performed the numerical computations. M.J.A., I.H. and H.S. directed various aspects of the theory and modelling. All authors contributed to writing the paper.

**Funding:** This research was funded by the Engineering and Physical Sciences Research Council (EPSRC) under Grant No. EP/M024237/1.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
