# Modeling the Lamb mode-coupling constant of quantum well semiconductor lasers

Arthur Vallet, Laurent Chusseau, Fabrice Philippe, Alain Jean-Marie, Gaëlle

Brévalle, Mehdi Alouini, Cyril Paranthoen

► To cite this version:

Arthur Vallet, Laurent Chusseau, Fabrice Philippe, Alain Jean-Marie, Gaëlle Brévalle, et al., Modeling the Lamb mode-coupling constant of quantum well semiconductor lasers. Optics Express, 2020, 28 (15), pp.21407-21419. 10.1364/OE.395790. hal-02893111

HAL Id: hal-02893111

https://hal.science/hal-02893111

Submitted on 8 Jul 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modeling the Lamb mode-coupling constant of quantum well semiconductor lasers

ARTHUR VALLET,<fnref n="1" /> LAURENT CHUSSEAU,<fnref n="1" />,\* FABRICE PHILIPPE,<fnref n="2" /> ALAIN JEAN-MARIE,<fnref n="3" /> GAËLLE BRÉVALLE,<fnref n="4" /> MEHDI ALOUINI,<fnref n="4" /> AND CYRIL PARANTHOËN<fnref n="4" />

**Abstract:** We theoretically compute the coupling constant *C* between two emission modes of an extended cavity laser with a multiple quantum-well active layer. We use an optimized Monte Carlo model based on the Markov chain that describes the elementary events of carriers and photons over time. This model allows us to evaluate the influence on *C* of the transition from a class A laser to a class B laser and illustrates that the best stability of dual-mode lasers is obtained with the former. In addition, an extension of the model makes it possible to evaluate the influence of different mode profiles in the cavity as well as the spatial diffusion of the carriers and/or the inhomogeneity of the temperature. These results are in very good agreement with previous experimental results, showing the independence of *C* with respect to the beating frequency and its evolution versus the spatial mode splitting in the gain medium.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Microwave photonics [1,2] requires compact and stable sources to generate frequencies from a few GHz to the millimeter wave range. The direction today is towards integration [3] and in the future this trend will be accentuated with increasingly smaller elementary components.

Optical beating using a dual-mode laser in combination with a high speed detector has shown to be very attractive for the generation of microwave frequencies, possibly up to the THz frequency range as has already been done for telecommunications [4]. Integrated dual-mode lasers have many advantages over two separate lasers, mainly due to their single optical cavity which allows compactness and insensitivity of the frequency difference to any first-order drift in the cavity length [5]. As a result high-purity microwave sources using dual-mode lasers have already been demonstrated [6]. The road to more compact devices uses semiconductor gain materials operating at 1.55 µm to take advantage of the enormous technological effort already made for optical communications.

These sources are our case study, dual wavelength operation is achieved by designing two distinct resonant frequencies on the semiconductor chip. Because they share the same gain material, these two sources are subject to mode competition and instabilities [7] whose counterpart is the intrinsic phase-to-noise correlation or antiphase noise [8,9] which can be turned to advantage since it reduces the RF beatnote noise [10].

<footnote n="1">¹IES, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="2">²LIRMM, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="3">³INRIA, LIRMM, Université de Montpellier, CNRS, Montpellier, France</footnote>

<footnote n="4">⁴FOTON, Université de Rennes, INSA, CNRS, Rennes, France</footnote>

<footnote>*laurent.chusseau@ies.univ-montp2.fr</footnote>
Since the earliest theory of gas lasers [7], the dual-mode emission is well understood starting from gas lasers rate equations. If mode intensities, $I_1$ and $I_2$, follow the time evolution

$$ \frac{dI_1}{dt} = (\alpha_1 - \beta_1 I_1 - \theta_{12} I_2) I_1, \qquad (1a) $$

$$ \frac{dI_2}{dt} = (\alpha_2 - \beta_2 I_2 - \theta_{21} I_1) I_2, \qquad (1b) $$

where the $\alpha_i$ are unsaturated gains and $\beta_i$ and $\theta_{ij}$ are self-saturation and cross-saturation coefficients, then the stability analysis is governed by the Lamb constant

$$ C = \frac{\theta_{12}\theta_{21}}{\beta_1\beta_2}. \qquad (2) $$

Depending on whether this constant is less or greater than 1, the two modes oscillate simultaneously or not. This model has been established within the restricted framework of “class A laser”, a system that classifies lasers according to their respective dynamics between photons, polarization and population inversion. A modern discussion of these classes and their relation to high-order photon correlations is given in [11].

Although dual-mode lasers using a solid-state active region are well known [12–14], their transposition to a semiconductor gain material has only appeared in recent years [10,15,16], the latter case involving InAs/GaAs quantum dots (QD). The procedure applied to experimentally extract C on such solid state lasers [13,17], has recently been transposed only twice to semiconductor based dual-mode lasers, consisting mainly of quantum wells (QW) active layers in Vertical External Cavity Surface Emitting Laser (VECSEL) architecture [18,19].

In the first case, InGaAs/GaAsP QWs provide optical gain at 1 µm and a variable spatial separation is introduced between the two mode gain zones to reduce competition between them. A coupling constant of value C = 0.8 is extrapolated for a perfect overlap. In the second case, the dual-mode laser has been built with a perfect overlap between modes and includes InP/InAlGaAs QWs operating at 1.54 µm. A coupling constant of C = 0.84 ± 0.02 was measured. Although high, these C-values indicate a stable behavior that was not expected with models of QW lasers based on population rate equations but was strongly anticipated for QDs [20,21]. As in these models the carrier population is not adiabatically eliminated to account for the laser class B behaviour, these models do not really take into account real experimental situations. A refinement has recently been proposed using Monte Carlo simulation of the dual-mode regime taking into account the light-matter interaction directly by quantum jumps between the populations [22] and the possibility of a stable regime has been suggested. Nevertheless, the calculations in this first study were too long to investigate the question in more detail with respect to the laser parameters.

Due to the recent improvement of the Monte Carlo model [23], we are able to derive in the following a numerical strategy to illustrate the stability of dual-mode lasers with active QW regions. We achieve this by evaluating the constant C with a numerical procedure that replicates the experimental procedure [13,18,19].

## 2. Material and methods

The numerical evaluation of C using simulations is based on the variations in losses inserted selectively for each mode, thereby mimicking the experience [13]. Although complicated experimentally, it is conceptually easy to introduce in calculations since the cavity lifetime of each mode is a separate parameter. Within the framework of (1), the steady-state optical intensities are $I_i$ and the optical losses are given by $\alpha_i$. Assuming that the numerical conditions for a stable steady-state dual-mode regime are obtained, the intensity variation of mode 1 while modulating the losses of mode 2 is thus given by $\partial I_1 / \partial \alpha_2$ and should be divided by the modulated intensity
of mode 2 in the same condition $\partial I_2/\partial \alpha_2$. The converse applies when modulating the linear gain of mode 1 and yields

$$C \equiv K_{12}K_{21}, \quad \text{with} \quad K_{ij} \equiv -\frac{\partial I_i}{\partial \alpha_j} \left(\frac{\partial I_j}{\partial \alpha_i}\right)^{-1}. \qquad (3)$$

Before the previous procedure can be applied, it is necessary to obtain a steady-state operation exhibiting average photon number approximatively equals for both modes. This choice of iso-intensity operation is dictated by the experimental need for maximum beatnote contrast yielding a maximum emission. This point is called the standard dual-mode operation point and its search is sometimes cumbersome as discussed in a previous work [22]. It is only around this point that we can reproduce the experimental procedure to obtain the $K_{ij}$, which are the cross-to self-saturation ratio, by varying one by one the two cavity losses in a tiny range with aim to extract the slopes and then C.

## 2.1. Compact Monte Carlo model

Monte Carlo and stochastic laser models calculate the time-dependent traces of the intensity evolution. They have been developed most of the time from the laser master equation [24,25] and were often applied to nanolasers with aim to understand their noise performance [26,27] or specific transition to lasing [23,28,29] because fluctuations are intrinsically taken into account. Recently we studied the mode competition within a dual-mode QW laser using our specific model that includes the carrier competition within bands via their thermal relaxation and Pauli's exclusion principle [22]. In spite of extremely long numerical computations that barely allowed us to touch the modal competition in QW lasers, we were able to show that contrary to rate equation models that predict a bistable regime [20], a stable two-mode regime could be possible by adjusting the cavity losses separately on the two modes, even if this adjustment is extremely sensitive.

To go further, improvements were needed to significantly reduce computing requirements. This was done in another work and applied to the emission properties of nanolasers during the transition to lasing. We extend it here to QW dual-mode semiconductor lasers and recall only its main features which are detailed in [23]. Our picture of a QW dual-mode semiconductor laser consists of two quantized electromagnetic modes resonant with a dual-mode cavity containing $m_1$ and $m_2$ photons and a finite number of electrons occupying the valence band and the conduction band. A finite number of energy values is available for each electron, only two of such levels in each band may allow radiative transitions. Owing to level discretization and introducing thermalization constants it is possible to follow each particle (electron or photon) with a huge Markov chain whose transition rates are completely defined by the present state [20,26]. Because typical time constants are far shorter for thermalization than for photonic processes, we split the Markov chain in a rapid part, the thermalization, that is treated only once analytically using the framework of the canonical ensemble of statistical mechanics, and a slower one, the photonic and pumping events, that is calculated using the exact algorithm of Gillespie [30].

Eventually the calculation requirements were reduced by a factor of more than 5 000, without loss of accuracy [23]. Thanks to this gain, the first results proposed in [22] can now be extended by adding a recursive procedure to automatically find the balance between the modes by adjusting the cavity losses and finding an appropriate couple $\alpha_i^{\text{bal}}$. As a result, a redundancy in the calculation is used to extract valid statistical results and not simple averages over a few trajectories as was done in [22] due to the heavy numerical calculations. The Fig. 1 illustrates what is generally obtained for the optical intensities by varying the optical losses $\alpha_i$ around these equilibrium values once the equilibrium between the modes $I_1 = I_2$ inside the cavity has been obtained at $\alpha_i^{\text{bal}}$. The mode competition in these lasers is so strong that we had to calculate a large number of Monte Carlo trajectories to extract the slopes from their average intensities with sufficient
precision, then the $K_{ij}$ and finally C. For each $\alpha_i$ the point reported in Fig. 1 is the statistical result of 40 runs and thus exhibits an uncertainty that sometimes includes negative intensities, which highlights the coexistence of On and Off laser operations within each trajectory of the set of 40 runs. It is interesting to note that this uncertainty is even greater at the vicinity of the mode balance, reflecting the modal competition and the more frequent mode switching around this point than when the laser is completely unbalanced, for instance at the extrema of the $\alpha_i$ considered. In practice we performed several tests before finding the acceptable solution reported in Fig. 1 for which a total of more than 10 000 runs has to be performed to extract the following values for the Lamb parameters: $K_{12} = 0.863 \pm 0.024$, $K_{21} = 1.121 \pm 0.031$ and $C = 0.97 \pm 0.04$ with specified 95% confidence intervals. Numerically, each trajectory covers more than 3 $10^5$ pump events whose statistics is poissonian to model optical pumping. An equivalent number of photonic events occurs in the cavity for each trajectory, and the total cost in calculation time including the mode-balance search and the 10 000 runs with their respective initializations was 75 minutes by requesting 40 cores on a cluster.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Evolutions of the optical intensities $I_1$, red, low energy mode, and $I_2$, blue, high energy mode of a nanolaser while varying the optical losses $\alpha_i$ once at a time. Intensity unit is the average number of photons in the cavity for each mode. Optical loss unit is the inverse of the time unit which is fixed in our model to the lifetime of the population in the excited state, generally 1 ns for a semiconductor laser. Laser gain is thus unity and assumed identical on both modes. Dots are averages over 40 runs at each $\alpha$ value and the error bars are 95% confidence intervals. Pump rate is 200 per unit time. The fraction of spontaneous emission coupled into the lasing mode $\beta$ is 0.4, which corresponds to a nanolaser with low threshold. (a) Variation of $\alpha_1$ at fixed $\alpha_2^{\text{bal}} = 0.439$. (b) Variation of $\alpha_2$ at fixed $\alpha_1^{\text{bal}} = 0.5$. Lines are fits whose slopes give the $K_{ij}$ of (3). Other laser parameters are $B = 800$ levels per band, $T = 300$ K, $\epsilon = 1$ meV level spacing within bands and a mode spacing of 4 meV, i.e. $\approx 1$ THz. Laser model details and details on how parameters are defined may be found in [23].</figcaption></figure>

## 2.2. Two-dimensionally coupled Monte Carlo model

To get closer to experiments we developed a more complex model coupling many elementary emitters as described by the previous compact model. The objective is to address any spatial effects, whether due to competition between adjacent laser zones, inhomogeneous temperature or transverse carrier scattering in the VECSEL, all of which can influence bimodal stability. Figure 2 gives a schematic representation of the 2D laser model. It is built on a meshed arrangement of elementary emitters individually described in §2.1 whose parameters may depend on the
spatial position. Each of the elementary emitters can exchange carriers with its closest neighbors, to account for electronic scattering, but all share the same photon bath since the latter is not localized in a conventional laser operating well above threshold. Eventually, the gain of each of the two modes can also be differentiated for each laser, allowing to take into account the imperfect overlap of the two modes, either desired [18] or resulting from waist differences due to possible different optical paths in the cavity [19].

<figure><img src="image_4.png" /><figcaption>Fig. 2. (a) Schematic representation of the 2D laser model. (b) Simple chain of carrier reservoirs exchanging particles.</figcaption></figure>

Aside from a simple juxtaposition of emitters, which can all exhibit different intrinsic parameters, the electron transfer requires to introduce novel rates in the Markov chain governing that transition for each elementary emitter with its neighbors. Let’s consider a simple chain of three elements A, B and C as depicted in Fig. 2(b), this rate is set proportional to the number of carrier in each element, and for B it is given by

$$n_B(t + \delta t) = n_B(t) + \mathcal{D}\delta t(n_A(t) + n_C(t) - 2n_B(t)), \quad (4)$$

where $\mathcal{D}$ is the diffusion constant in the system. A simple justification of (4) is given in Appendix A and its generalization in 2D is straightforward considering only the four closest neighbors.

Inhomogeneous temperature also influences carriers with an obvious transfer from hot devices toward cold ones. We account it as a modification of $\mathcal{D}$ in (4) by means of the Metropolis criterion for Monte Carlo. Considering two adjacent reservoirs A and B of temperatures $T_A$ and $T_B$, the $\mathcal{D}$ factor of transition from A to B is multiplied by the probability min $[\exp(1 - T_B/T_A), 1]$, thus it is favored if $T_A \ge T_B$ but still allowed with a reduced rate if $T_B > T_A$.

While we had gained a very important factor exceeding 5 000 in computation time with our ad-hoc processing of the Markov chain describing the laser, the reintroduction of multiple active elements and their coupling results in single trajectory computation times of several minutes on our multi-core computer, thus requiring many hours for averaging and producing the numerous traces required to extract the value of C. Although carrier diffusion and temperature inhomogeneity were included in the final open source software [31], their computational time sometimes made it impossible to evaluate them all.

### 3. Results

First of all, we used the compact model described in §2 to evaluate the sensitivity of C to some typical laser parameters. Let’s consider $\beta$ which governs the emission ratio in the laser mode and
therefore the laser threshold. The typical case already illustrated in Fig. 1 has been reproduced for $\beta$ values ranging from $10^{-2}$ to 1, which corresponds to the value for a well confined VCSEL and that of a perfect nanolaser. The cross- to self-saturation coefficients $K_{ij}$ and the Lamb coupling constant $C$ are reported with respect to $\beta$ in Fig. 3 following the statistical analysis on more than 10 000 trajectories per point discussed in Fig. 1.

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Influence of $\beta$ on (a) $K_{12}$ (red) and $K_{21}$ (blue) and (b) $C$ calculated using the Monte Carlo model. The span considered covers the transition from a VCSEL to a perfect nanolaser.</figcaption></figure>

At first glance all the mean values remain insensitive to $\beta$, which is not the case for the error bars produced by the statistical analysis. In any case the dynamical operation of the laser is calculated as stable dual-mode even if $C$ is very close to unity. Close examination of the trajectories shows in this case very noisy intensities for each mode but with noise clearly in antiphase and therefore, as expected, the beat intensity noise is much lower than that of the individual modes. High $C$ values are thus of no penalty for applications as long as only the beat signal is used (See Appendix B). Anyway, as $\beta$ decreases, the error bars increase more and more to exceed by far the critical value $C = 1$, which is a sign of increasing instability. It has been shown in [23] that low values of $\beta$ promote a mode shutdown whose duration is inversely proportional to its value. This results in large emission peaks when this mode is re-ignited, causing a large increase in photon number variance and in turn the observed increasing error bars. We thus admit that such behavior is intimately provoked by the very small number of photons and electrons considered in our Monte Carlo simulation. If we consider that a significant increase in this number of electrons and photons is impossible in the simulations without making the computation times prohibitive, and if we refer to the computed mean values that do not vary significantly over such a large scale, we must conclude that $\beta$ has a very small influence on the effective stability, so we will assume in the following that the computed behavior for a typical $\beta$ of 0.5 (or 0.4 in Fig. 1) is representative of the experimental VECSEL.

The dynamics of semiconductor lasers is very dependent on their class, and therefore on their optical cavity. A VECSEL is intrinsically class A because of its large cavity, which considerably increases the photon lifetime, whereas an integrated VCSEL laser is clearly class B with a slower carrier dynamic than that of photons. It can therefore be expected that the dual-mode behaviour may also depend on the class considered for the study laser. We have set up a simulation to illustrate this transition. Keeping the parameters of Fig. 1 with $\beta = 0.5$ and increasing step-by-step $\alpha_1$, the cavity losses of the low energy mode, the class of the laser gradually changes from A to B. Results are given in Fig. 4.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Calculated (a) $K_{12}$ (red) and $K_{21}$ (blue) and (b) *C* constant as a function of $\alpha_1$ the ratio between cavity lifetime to population lifetime. From left to right the evolution here mimics the transition from a class A to a class B laser.</figcaption></figure>

For numerical convergence reasons, $\alpha_1$ was chosen in the narrow range $[10^{-1}, 1]$ and the solution $\alpha_2$ satisfying the mode balance was determined each time before computing the $K_{ij}$ and Lamb constant *C* using statistics over 10 000 trajectories. Actually, $\alpha_2$ is proportional to $\alpha_1$ since it can be adjusted according to the linear law $\alpha_2 = 1.06\alpha_1 - 0.089$. The negative constant in the fit law thus prohibits dual mode behavior if $\alpha_1$ becomes too small, and it is also responsible for the growing gap between the powers of the two modes when $\alpha_1$ decreases.

A regular transition with rapidly evolving $K_{ij}$ values is obtained (see Fig. 4(a)). At the smallest $\alpha_1$, which corresponds to a class A laser with a carrier dynamic ten times faster than that of photons, $K_{12}$ is already very close to zero. This explains why it was not possible to consider lower $\alpha_1$ without making the laser bistable. Conversely, this is the point where $K_{21}$ is the highest, showing an extreme sensitivity of the higher energy mode intensity to the losses of the low energy one. At the opposite of the curve where $\alpha_1 = 0.79$, the two constants $K_{12}$ and $K_{21}$ tend towards the same value very close to unity, showing an equivalent sensitivity for one or the other of the two modes towards modifications made on its own losses or on those of his companion. At that point the carrier and photon dynamics are now almost equivalent, a situation close to that of a class B laser. Again, it was not numerically possible to push $\alpha_1$ to values far beyond without resulting in unreasonable pumping.

As seen in Fig. 4(b), in the same time the Lamb constant *C* increases regularly. For class A lasers, we obtain our lowest value of *C* = $0.85 \pm 0.12$ which is in very good agreement with published experiments [18,19]. On the other hand, moving towards the class B laser, *C* increases steadily and becomes so close to unity that it excludes any practical stable dual-mode regime, according to rate-equation models [20]. In our view, the fact that bulk semiconductor or QW lasers are intrinsically class B makes them unsuitable for dual mode operation. Only the addition of an external cavity to make them class A and eventually the advent of spatial hole burning makes this possible. This is not the case if quantum dots are used instead because then intrinsic mode decoupling occurs through the inhomogeneous gain broadening, and we have high expectations for the realization of compact and stable dual-mode lasers with this type of gain material [16,21].

The most complete model, with its 2D coupled Monte Carlo emitters, is expected to provide further insights, especially as it supports the spatial dimension of VECSEL lasers. First of all, coupling parameters were introduced one by one to see their influence on $K_{ij}$ and *C*. In practice,
it was not possible to go as far as the situation described in the Fig. 2, because then the Markov chain to be simulated explodes in size. In fact the number of transitions to be taken into account is multiplied by the number of individual elements included. It grows even bigger if coupling is introduced since the rate in (4) is proportional to the total number of carriers considered and thus becomes very high if $\mathcal{D}$ is high. In practice, we limited calculations to a few emitters that we considered regularly spaced rather than spread out on a grid, because our maximum allowable computer time has thus been reached. Our typical situation includes 4 of such individual Monte Carlo lasers for which we consider the influence of carrier coupling, pump intensity, lasing level spacing and mode overlap. We will see that temperature does not matter in our framework.

The Fig. 5(a) illustrates the influence of carrier coupling between emitters. Although the value of this coupling is varied over more than 8 decades, the $K_{ij}$ and $C$ are fairly constant over the whole range of variation and identical to that of no-coupling, and only differences in the uncertainties are noticeable. As concerns the final stability of the dual mode regime defined by the Lamb constant $C$, it is almost unchanged, even in the case of strong spatial carrier exchanges between emitters. We therefore conclude that carrier diffusion has no influence in our calculations, which is beneficial because it is time consuming to take into account. In the same way and in order to lighten the calculations we have ignored in the following the effects of temperature as previously discussed as a modification of the diffusion coefficient. However, it should not be ignored that temperature also acts on many material parameters including the gain and loss of optical modes as well as the position of gain maximum. Intense optical intensity in the cavity can also influence the gain through non-linear saturation or through heating due to higher absorption. The first of these effects was demonstrated in our first Monte Carlo simulations using the microscopic model which demonstrated a shift in mode balance with optical intensity [22]. This is no longer possible with the canonical formalism used here [23], which is the only way to carry out the amount of computation required by this study. However, if the effect of temperature is ignored in the following, it would still be possible to add it as a phenomenological gain model.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Calculated dual-mode coupling constants (a) versus $\mathcal{D}$, the carrier coupling parameter and (b) the pumping intensity $J$ in carrier per unit time. The laser simulated includes 4 emitters like the one of Fig. 1 with $\beta = 0.5$. In (a) the value determined without coupling is indicated by dashed colored horizontal lines. For both, the dotted gray line figures the stability limit for $C$.</figcaption></figure>
The Fig. 5(b) accounts for the influence of optical pumping level, *J*, in the same 4 emitter system. Again we varied *J* over a wide range with only tiny changes either on the *K*<sub>ij</sub> or on *C*. With the exception of side effects of pumping such as temperature increase which influences the gain and refractive index of the active material, two effects that we do not take into account in our model, we can therefore neglect the influence of pumping on the stability of the dual mode regime.

Closer to the experiment [19] is the evaluation of *K*<sub>ij</sub> and *C* when the energy spacing between the two lasing modes is varied. This is reported in Fig. 6 together with the experimental results. The parameter values are the same as in Fig. 1 including 4 emitters in the same cavity with identical gains on both modes. In our model, the energy spacing between modes is an integer multiple of the discretization of energy levels in the conduction band or valence band. This discretization is usually 1 meV but we have separately evaluated in a few cases what happens with half this value and found exactly the same results. This is the reason why calculations are given in Fig. 6 starting from 1 meV up to 5 meV. Although measured *K*<sub>ij</sub> values exhibit variations with mode spacing, the variation of calculated *K*<sub>ij</sub> is larger. However, they are very similar in average, at least between 2 and 4 meV of mode separation. With a mode separation of 5 meV the calculated *K*<sub>12</sub> is very low. This situation already occured in Fig. 4 with the lowest values of *α*<sub>1</sub>. Again, this comes along with a large *K*<sub>21</sub> and an extreme sensitivity to parameters for the simulated laser. Nevertheless, like experiments, the resulting Lamb constant *C* is almost independent of mode spacing, although its calculated value always exceeds that of the measured value.

<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Calculated (lines) and measured (filled dots and squares) [19] (a) *K*<sub>12</sub> (red) and *K*<sub>21</sub> (blue) and (b) *C* constant as a function of energy spacing between the two modes. Open symbols are calculations given with surrounding colored bands for 95% confidence levels.</figcaption></figure>

Finally we used the most complete model with multiple emitters to evaluate the dependence of decoupling on mode overlap in the active semiconductor region. This point was experimentally investigated in [18]. Numerically, we chose 4 different emitters in the same cavity with mode gains according to a prescribed spatial splitting *d* within the active medium. Again the laser parameters of Fig. 1 with *β* = 0.5 were used. To represent the overlapping modes we set these modal gains according to gaussian profiles for the two modes. This is represented in the insets of Fig. 7 for selected values of splitting distances *d*. Furthermore and in order to keep similar overall gain for the two modes, their total gains were equalized. Nevertheless, the method accounts
carefully of the experimental condition used in [18] using only a restricted number of emitters. Our simulated points and the corresponding fitting curve $C_0 \exp - (\frac{d}{w_0})^2$ expected with gaussian profiles [18] are plotted in Fig. 7 together with the experimental data.

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Lamb C constant of two non-overlapping gaussian modes as a function of their effective separation *d* in gain medium. Four emitters are set in the cavity with different mode gains adjusted to emulate the behavior of two non-overlapping gaussian modes. Monte Carlo results are the brown dots with error bars and the brown curve is a gaussian fit. Green stars are measurements from Ref. [18]. The insets for different separation distances are mode profile drawings where the dots indicate the exact position and relative intensity for each emitter while the curves are eye guides corresponding to the simulated continuous mode profiles. Total gain was equalized for each mode and *d*.</figcaption></figure>

As anticipated *C* decreases with *d* following a gaussian fit whose starting point at *d* = 0 µm is *C*₀ = 0.89 and radii *w*₀ = 34 µm. When the spacing between modes increases, the Fig. 7 exhibits an excellent overall agreement with the experimental results even if *C*₀ is slightly higher than that of the two reported experimental values [18,19]. Our complete model is thus an accurate tool for predicting laser dynamics, including mode competition in a device as complex as a dual-mode VECSEL with adjustable mode separation in the gain medium. Notice that we obtained in Fig. 6(b) larger Lamb C factors in what seems to be the same situation in Fig. 7 with four emitters at perfect overlap. The difference here comes from the gain values which are no more identical in Fig. 7 for all four emitters but decrease as we move away from the beam center, according to the first inset of Fig. 7. Therefore, taking into account a realistic mode profile would probably slightly alter the results in Fig. 6 in favor of a better agreement with the experiment.

**4. Conclusions**

Owing to drastic improvements in the Markov chain which represent a dual-mode laser, we have set up very efficient Monte Carlo simulations of dual-mode semiconductor lasers and have considerably extended the explored range of parameters compared to previous work [22]. The significant increase in speed obtained on the program makes it possible to multiply the simulations, allowing for the determination of the Lamb factor *C* for quantum well lasers by mimicking the experimental modulation of one mode at a time and observing all intensity modulations. A complete QW semiconductor laser simulator has been built, able to take into account multiple modes with different spatial and spectral profiles and ultimately lateral carrier diffusion and inhomogeneous temperature.

Numerical results show that some parameters do not influence the dual-mode stability, for instance the spontaneous emission factor $\beta$ or lateral carrier diffusion, whereas the intrinsic
nature of the laser (class A or class B) is of tremendous importance. Up until now, experimental evidence on the stability of dual-mode QW lasers has apparently contradicted model predictions. This originated from the discrepancy in dynamical properties since the measurements were taken in class A VECSELs while models were tailored to compact class B semiconductor lasers.

The application of the model to previously published experimental situations is done with success and good precision, both concerning the independence of C with the spectral separation of the two modes and its decrease when optical modes are spatially separated in the gain medium. Nevertheless, all this required a significant computational effort, and as the possibility of extending the computation to more complex lasers has been included in the software, this work will deserve to be expanded in the future for specific applications.

## A. Justification of (4)

Carrier diffusion is usually expressed in one dimension using

$$ \frac{\partial n(x, t)}{\partial t} = D \frac{\partial^2 n(x, t)}{\partial x^2}, \qquad (5) $$

where D is the diffusion coefficient.

The finite difference schemes for derivatives in the left-hand and right-hand sides yields

$$ \frac{\partial n(x, t)}{\partial t} = \frac{n(x, t + \delta t) - n(x, t)}{\delta t} + o(\delta t), \qquad (6a) $$

$$ \frac{\partial^2 n(x, t)}{\partial x^2} = \frac{n(x - \delta x, t) - 2n(x, t) + n(x + \delta x, t)}{\delta x^2} + o(\delta x^2), \qquad (6b) $$

which, inserted in (5) yields

$$ n(x, t + \delta t) = n(x, t) + \delta t D \left( \frac{n(x - \delta x, t) - 2n(x, t) + n(x + \delta x, t)}{\delta x^2} \right) + o(\delta t + \delta x^2). \qquad (7) $$

This result is the exact analogous of (4) for the discrete chain of Fig. 2(b) provided that $D = D/\delta x^2$.

## B. Estimate of the beating noise from Monte Carlo trajectories

Since the objective of beating two modes of a dual-mode laser is to obtain a pure frequency, the question of beatnote noise is a major concern. This appendix proposes a noise estimate starting from the intensity fluctuations of the two modes as obtained from Monte Carlo trajectories. In essence, this result is a lower bound for the beating noise since it neglects any phase noise. However, it illustrates how the intrinsic noise of each mode may partly cancel in the beating because of antiphase coupling. Additional improvements would, for instance, require keeping track of the temporal carrier evolution in order to add phase fluctuations via the Henry amplitude-phase coupling factor in semiconductor lasers. These improvements could be introduced into Monte Carlo models in future work, following the procedure introduced in [32].

Let us start from the electromagnetic fields of the two modes in the laser cavity

$$ E_1(t) = \mathcal{E}_1(t) \exp(j\Omega_1 t), $$

$$ E_2(t) = \mathcal{E}_2(t) \exp(j\Omega_2 t). $$

If we define $\Omega_0 = (\Omega_1 + \Omega_2)/2$ the average optical pulsation and $\Omega = (\Omega_1 - \Omega_2)/2$ the half of the beating pulsation, the total field is

$$ E(t) = E_1(t) + E_2(t) = \mathcal{E}_1(t) \exp(j(\Omega_0 + \Omega)t) + \mathcal{E}_2(t) \exp(j(\Omega_0 - \Omega)t). \qquad (8) $$
Since the electromagnetic energy is $E(t)E^*(t)$, we have

$$E(t)E^*(t) = \mathcal{E}_1(t)\mathcal{E}_1^*(t) + \mathcal{E}_2(t)\mathcal{E}_2^*(t) + \mathcal{E}_1(t)\mathcal{E}_2^*(t) \exp(j2\Omega t) + \mathcal{E}_1^*(t)\mathcal{E}_2(t) \exp(-j2\Omega t),$$

which simplifies by regrouping the last two terms

$$E(t)E^*(t) = \mathcal{E}_1^2(t) + \mathcal{E}_2^2(t) + 2\mathcal{E}_1(t)\mathcal{E}_2(t) \cos(2\Omega t). \quad (9)$$

In essence the $\mathcal{E}_i(t) \propto \sqrt{m_i(t)}$ are real stochastic process, therefore the number of "photons" at the beating frequency $2\Omega$ inside the cavity is

$$u_{2\Omega} \propto \mathcal{E}_1(t)\mathcal{E}_2(t) = \sqrt{m_1(t)m_2(t)}. \quad (10)$$

The result of (10) is illustrated using the particular trajectory of Fig. 8(a) where we have plotted the stochastic process $u_{2\Omega}$. As visible from the figure, the geometric average representing the beat intensity is considerably smoother than the highly fluctuating individual modes. This can be attributed to the predominantly antiphase nature of the large individual fluctuations which considerably reduce the variations in total photon number.

<figure><img src="image_10.png" /><figcaption>Fig. 8. (a) Example of calculated trajectories – photon number versus time – in dual-mode regime and (b) corresponding $g^{(2)}(\tau)$ versus time for the three processes. Red: low-energy mode $m_1$. Blue: high-energy mode $m_2$. Orange: beating process $u_{2\Omega}$.</figcaption></figure>

As a result the corresponding second-order delayed autocorrelation $g^{(2)}(\tau)$ – the most adequate statistical indicator [11] in the presence of a few intracavity photons per mode – is very close to unity for the beatnote, proving good coherence, as opposed to the larger noise figures of the individual modes (Fig. 8(b)).

## Funding

Agence Nationale de la Recherche (Swiss-French IDYLIC contract ANR-15-CE24-0034).

## Disclosures

The authors declare no conflicts of interest.

## References

1. A. J. Seeds and K. J. Williams, “Microwave photonics,” J. Lightwave Technol. 24(12), 4628–4641 (2006).

2. J. Yao, “Microwave photonics,” J. Lightwave Technol. 27(3), 314–335 (2009).

3. D. Marpaung, J. Yao, and J. Capmany, “Integrated microwave photonics,” Nat. Photonics 13(2), 80–90 (2019).

4. T. Nagatsuma, G. Ducournau, and C. C. Renaud, “Advances in terahertz communications accelerated by photonics,” Nat. Photonics 10(6), 371–379 (2016).
5. L. Chusseau, G. Almuneau, L. A. Coldren, A. Huntington, and D. Gasquet, "Coupled-cavity vertical-emitting semiconductor-laser for continuous-wave terahertz emission," *IEE Proc. J* **149**, 88–92 (2002).

6. G. Pillet, L. Morvan, M. Brunel, F. Bretenaker, D. Dolfi, M. Vallet, J.-P. Huignard, and A. Le Floch, "Dual-frequency laser at 1.5 µm for optical distribution and generation of high-purity microwave signals," *J. Lightwave Technol.* **26**(15), 2764–2773 (2008).

7. M. Sargent, I. I. I. M. O. Scully, and W. Lamb, *Laser Physics* (Addison-Wesley, 1974).

8. K. Otsuka, P. Mandel, S. Bielawski, D. Derozier, and P. Glorieux, "Alternate time scale in multimode lasers," *Phys. Rev. A* **46**(3), 1692–1695 (1992).

9. M. Wichmann, G. Town, J. Quante, M. Gaafar, A. Rahimi-Iman, W. Stolz, S. W. Koch, and M. Koch, "Antiphase noise dynamics in a dual-wavelength vertical-external-cavity surface-emitting laser," *IEEE Photonics Technol. Lett.* **27**(19), 2039–2042 (2015).

10. S. De, V. Pal, A. E. Amili, G. Pillet, G. Baili, M. Alouini, I. Sagnes, R. Ghosh, and F. Bretenaker, "Intensity noise correlations in a two-frequency VECSEL," *Opt. Express* **21**(3), 2538–2550 (2013).

11. N. Takemura, M. Takiguchi, and M. Notomi, "Probing the Ginzburg-Landau potential for lasers using higher-order photon correlations," *arXiv:1908.08679* (2019).

12. M. Brunel, F. Bretenaker, and A. L. Floch, "Tunable optical microwave source using spatially resolved laser eigenstates," *Opt. Lett.* **22**(6), 384–386 (1997).

13. M. Alouini, F. Bretenaker, M. Brunel, A. Le Floch, M. Vallet, and P. Thony, "Existence of two coupling constants in microchip lasers," *Opt. Lett.* **25**(12), 896–898 (2000).

14. M. Brunel, A. Amon, and M. Vallet, "Dual-polarization microchip laser at 1.53 µm," *Opt. Lett.* **30**(18), 2418–2420 (2005).

15. G. Baili, L. Morvan, M. Alouini, D. Dolfi, F. Bretenaker, I. Sagnes, and A. Garnache, "Experimental demonstration of a tunable dual-frequency semiconductor laser free of relaxation oscillations," *Opt. Lett.* **34**(21), 3421 (2009).

16. K. A. Fedorova, A. Gorodetsky, and E. U. Rafailov, "Compact all-quantum-dot-based tunable THz laser source," *IEEE J. Sel. Top. Quantum Electron.* **23**(4), 1–5 (2017).

17. A. McKay, J. M. Dawes, and J.-D. Park, "Polarisation-mode coupling in (100)-cut Nd:YAG," *Opt. Express* **15**(25), 16342–16347 (2007).

18. V. Pal, P. Trofimoff, B.-X. Miranda, G. Baili, M. Alouini, L. Morvan, D. Dolfi, F. Goldfarb, I. Sagnes, R. Ghosh, and F. Bretenaker, "Measurement of the coupling constant in a two-frequency VECSEL," *Opt. Express* **18**(5), 5008–5014 (2010).

19. G. Brévalle, S. Pes, C. Paranthoën, M. Perrin, C. Vallois, C. Hamel, A. Mereuta, A. Caliman, E. Kapon, A. Vallet, L. Chusseau, H. Folliot, and M. Alouini, "Direct measurement of the spectral dependence of Lamb coupling constant in a dual frequency Quantum Well-based VECSEL," *Opt. Express* **27**(15), 21083–21091 (2019).

20. L. Chusseau, F. Philippe, P. Viktorovitch, and X. Letartre, "Mode competition in dual-mode quantum dots semiconductor microlaser," *Phys. Rev. A* **88**(1), 015803 (2013).

21. L. Chusseau, A. Vallet, M. Perrin, C. Paranthoën, and M. Alouini, "Lamb mode-coupling constant in quantum dot semiconductor lasers," *Phys. Rev. B* **98**(15), 155306 (2018).

22. L. Chusseau, F. Philippe, and F. Disanto, "Monte Carlo modeling of the dual-mode regime in quantum-well and quantum-dot semiconductor lasers," *Opt. Express* **22**(5), 5312–5324 (2014).

23. A. Vallet, L. Chusseau, F. Philippe, and A. Jean-Marie, "Markov model of quantum fluctuations at the transition to lasing of semiconductor nanolasers," *Phys. E* **105**, 97–104 (2019).

24. R. Loudon, *The Quantum Theory of Light* (Clarendon, 1983).

25. P. R. Rice and H. Carmichael, "Photon statistics of a cavity-QED laser: A comment on the laser–phase-transition analogy," *Phys. Rev. A* **50**(5), 4318–4329 (1994).

26. L. Chusseau and J. Arnaud, "Monte-Carlo simulation of laser diodes sub-poissonian light generation," *Opt. Quantum Electron.* **34**(10), 1007–1023 (2002).

27. A. Lebreton, I. Abram, N. Takemura, M. Kuwata-Gonokami, I. Robert-Philip, and A. Beveratos, "Stochastically sustained population oscillations in high-β nanolasers," *New J. Phys.* **15**(3), 033039 (2013).

28. K. Roy-Choudhury, S. Haas, and A. Levi, "Quantum fluctuations in small lasers," *Phys. Rev. Lett.* **102**(5), 053902 (2009).

29. G. P. Puccioni and G. L. Lippi, "Stochastic simulator for modeling the transition to lasing," *Opt. Express* **23**(3), 2369–2374 (2015).

30. D. T. Gillespie, *Markov Processes: An Introduction for Physical Scientists* (Academic, 1992).

31. A. Vallet, "Laser2D: A two-band 'semiconductor-like' markovian multimode laser using the Monte-Carlo method," https://github.com/crazythutur/Laser2D (2019–2020).

32. S. De, V. Potapchuk, and F. Bretenaker, "Influence of spin-dependent carrier dynamics on the properties of a dual-frequency vertical-external-cavity surface-emitting laser," *Phys. Rev. A* **90**(1), 013841 (2014).