Non-equilibrium dynamics in the dual-wavelength operation of Vertical external-cavity surface-emitting lasers

Microscopic many-body theory coupled to Maxwell's equation is used to investigate dual-wavelength operation in vertical external-cavity surface-emitting lasers. The intrinsically dynamic nature of coexisting emission wavelengths in semiconductor lasers is associated with characteristic non-equilibrium carrier dynamics which causes significant deformations of the quasi-equilibrium gain and carrier inversion. Extended numerical simulations are employed to efficiently investigate the parameter space to identify the regime for two-wavelength operation. Using a frequency selective intracavity etalon, two families of modes are stabilized with dynamical interchange of the strongest emission peaks. For this operation mode, anti-correlated intensity noise is observed in agreement with the experiment. A method using effective frequency selective filtering is suggested for stabilization genuine dual-wavelength output.


Introduction
Vertical external-cavity surface-emitting lasers (VECSELs) are wavelength tunable solid-state lasers that have been configured to produce either high-power/ultrashort mode-locked pulses or single-/multi-wavelength output [1][2][3][4][5][6][7][8][9]. Currently, the peak power produced in multi-mode operation is 106 W [5,6] and 23 W for single-mode continuous wave output [7]. VECSELs have been able to generate intracavity mode-locked single pulses of temporal duration around 100 fs [8][9][10], with the shortest intracavity pulse length at 95 fs [10]. A further unique aspect of these surface emitting lasers is that the single pass gain is low due to light amplification through a periodic or quasi-periodic stacking arrangement of thin (typically 8-10 nm thick) quantum wells (QWs). Consequently, one needs high mirror reflectivities to initiate lasing, meaning that they need to be strongly optically pumped and can sustain very strong intra-cavity circulating fields. As a result these are ideal semiconductor laser systems in which to explore strongly non-equilibrium many-body effects where carrier distributions (electrons, holes, phonons) are strongly perturbed from the standard quasi-equilibrium Fermi distributions.
For a wide range of applications, it is desirable to operate VECSELs in a regime that allows for stable two-wavelength emission. Relevant examples include spectroscopy [11,12], differential absorption LIDAR with multiple applications from detection of geological fault lines to measuring the atmospheric CO 2 concentration [13,14], semiconductor laser parameter measurement [15,16], THz-imaging in contexts ranging from security applications to detection of brain cancer cells [17][18][19][20], and THz signal generation [21,22]. In VECSELs THz output has been produced using both intracavity difference frequency generation (DFG) and mode-locked pulses centered at different frequencies with THz spacing [23,24].
In the past, two-wavelength operation in semiconductor lasers has been studied both experi-mentally and theoretically [24][25][26][27][28][29][30][31][32][33][34][35][36]. Stable dual-wavelength operation of VECSELs has been experimentally realized by multiple groups [24][25][26][27][28][29][30]. In particular, experiments show that there is persistent anti-correlated intensity dependent noise in the dual-wavelength output spectra [29,30]. Moustafa Ahmed and Minoru Yamada investigated multi-mode operation in AlGaAsâĂŞGaAs and InGaAsPâĂŞInP laser systems using multi-mode rate equations coupled to Langevin noise sources [33]. They found multiple different regimes with stable and unstable single-and multi-mode operation, which they characterize according to injection current and the linewidth enhancement factor. In Ref. [34], Yacomotti et al. have investigated multi-longitudinal-mode operation in semiconductor lasers experimentally and numerically using rate equations for the optically active gain medium. They observed strong anti-correlated mode dynamics between two modes, while the laser emitted constant total intensity. They concluded that four-wave mixing was the dominant cause of the observed phenomena.
Most of the theoretical analysis is based on rate equation approaches [31][32][33][34][35], where the optical polarization dynamics has been adiabatically eliminated and the material excitation is purely modeled by dynamically varying carrier densities. As a consequence, the induced electron-hole polarizations in the QWs are not resolved with the consequence that the ultrafast non-equilibrium light-matter coupling dynamics in VECSEL systems cannot be treated in detail. Moreover, these rate equation approaches rely on added parameterization (although well motivated physically) and directly treat the laser gain -the latter by integrating over the non-equilibrium carrier distributions is insensitive to details such as deformations due to carrier scattering and interaction with the lasing field.
To improve on the rate-equation analysis, Bäumner et al. investigated non-equilibrium carrier dynamics in a simple micro-resonator system under dual-wavelength operation [36]. In their paper, the authors simulate the Maxwell semiconductor Bloch equations (MSBE) at the Hartree-Fock level with simple rate approximations for carrier scattering and found that including non-equilibrium carrier dynamics was essential for accurate modeling of single-and dual-wavelength operation.
In this paper, we will expand the dual-mode operation analysis to VECSELs and numerically solve the MSBE with carrier scattering computed on the level of second Born-Markov approximation [37,38]. In the past, we have employed this model to systematically analyze the role of microscopic carrier dynamics in mode-locked pulse operation of VECSELs, see Refs. [39][40][41][42]. In mode-locked operation, pulses with temporal duration in the pico-to femtosecond range are common, which significantly limits the light-matter interaction time within the semiconductor quantum wells during which the detailed carrier dynamics has to be numerically evaluated. Ironically in contrast, a dual-wavelength intracavity field is continuously extracting QW carriers from the inverted states that have to be replenished by external pumping into absorbing states. Furthermore, the coexistence of several wavelength fields inside the cavity leads to dynamic sum-and difference frequency generation such that one never reaches a truly stationary regime. Hence, the high dimensional carrier scattering integrals have to be solved for the entire duration of the simulation, which poses a significant computational challenge but is unavoidable due to the strong non-equilibrium dynamics that influences the round-trip amplification of the intracavity field.

Microscopic Theory
In Fig. 1, we present a schematic outline of a VECSEL structure including the QWs, the output coupling mirror and an additional intra-cavity etalon that might be introduced in order to select a specific modal output. The microscopic model of such a VECSEL must fully resolve the time dynamics of the electromagnetic cavity field and the microscopic QW carrier dynamics, where the majority of computational complexity originates from the latter. In order to simplify the numerical challenge, our analysis considers only the propagating light field along the z-axis perpendicular to the QW planes. The electromagnetic field, E(z, t), inside the VECSEL cavity is modeled using Maxwell's equation, Here, c 0 is the speed of light in vacuum, µ 0 is the vacuum permeability, n(z) is the material refractive index, P(z, t) is the QW macroscopic polarization. A standard VECSEL cavity cross-section consists of multiple material layers such as: GaAs, AlAs, and an external air cavity, as shown in Fig. 1. We assume a constant index of refraction for each material layer, which allows us to simulate the propagating field by solving Eq. (1) inside each material layer and couple the layers together using standard Maxwell boundary conditions [36]. At the boundary of two different materials, the propagating field will experience reflection and transmission and, most importantly, the field can be amplified or absorbed by interacting with the dynamically changing QW carriers through the macroscopic polarization. An optically active QW can be modeled using the multiband semiconductor Bloch equations (SBE) [37], which describe the time evolution of the microscopic electron (hole) occupation numbers, n e(h) λ(ν),k , and the microscopic polarizations, p λ,ν,k , for a given crystal momentum k in the conduction (valence) band denoted by λ (ν) The macroscopic polarization, that feeds back into Eq. (1), is defined by P(z, t) = λ,ν,k d k λν p λ,ν,k , where d k λν is the dipole matrix element. The single particle electron (hole) energies are given by e(h),λ(ν) k and couple through the Coulomb potential, V λ,ν 1 ,ν,λ 1 |k−q | . The Hartree-Fock renormalized single particle energies are e e λ,λ 1 ,k = e,λ k δ λ,λ 1 − and the propagating electromagnetic field is a source in Eq. (2) through the effective Rabi frequency, Ω λ,ν,k , defined by In Eq. (2), the higher order correlation contributions give rise to the polarization dephasing (Γ λ,ν,deph ), the carrier-carrier and carrier-phonon scattering ( d dt n e(h) Here,V λ,ν 1 ,ν 2 ,λ 1 |q | is the screened multiband Coulomb matrix element and D(E) is an energy conserving function that evaluates to unity for allowed transitions and zero otherwise. The carrier-phonon scattering is evaluated from Where g λλ 1 q is the Fröhlich matrix element [43,44], s q is the phonon distribution, ω q is the LO-phonon energy. The phonon scattering is defined by specifying the lattice parameters: the LO-phonon energy ( ω LO = 33.95 meV), lattice temperature (300 K) and the dielectric constants at high-( ∞ = 11.24) and low-( = 13.47) frequencies. These values are determined from interpolation of table values found in Ref. [45].
In order to simplify the numerical challenge, we assume that the relevant portion of the detailed semiconductor bandstructure can be approximated by a two-band effective mass model with transition energy given by ω k = E g + population inversion, the polarization dephasing is approximated using a characteristic timescale, τ deph , such that Γ cv,deph = −(1/τ deph )p cv,k . In the dephasing time approximation, carrier screening will be treated using the static Lindhard formula [37]. With two parabolic bands, Eqs. (5-6) can be further simplified by using the analytic band structure to explicitly evaluate the innermost integral over the energy conserving delta function. Employing a delta function is particularly useful for carrier-carrier scattering because it reduces the number of integrals by one and allows for a considerable computational speedup.
Our simulation domain consists of a gain chip with ten QWs, an air gap, and an output coupling mirror. This gain chip consists of a GaAs/AlGaAs distributed Bragg reflector, a barrier region with 8 nm InGaAs QWs arranged on an anti-node of the DBR standing wave, a cap layer, and an anti-reflection coating. The gain chip is optimized for ultrashort pulse generation with a peak gain at 1030 nm, similar to the one used for Ref. [39]. However, mode-locked considerations turn out to be of less importance in the context of dual-wavelength generation.
The QW carrier replenishing is modeled with a relaxation to a background Fermi distribution, f e(h) c(v),k , at a fixed density (3.0 · 10 16 m −2 ) and temperature (300 The characteristic timescale for gain recovery is on the order of tens of picoseconds or longer, which is orders of magnitudes slower than the carrier scattering [46]. Numerical tests show that the details of the results are not critically sensitive on the exact choice of τ scatt . For the studies in this paper, we use τ scatt = 30 ps which, for appropriately chosen output coupling losses, results in reasonably high intracavity field intensities that are, however, still well below the CW damage threshold for GaAs [47][48][49].
The spontaneous emission of photons into the lasing cavity is approximated using n e k n h k , as used by Bäumner et al. in Ref. [36]. Here, , β is a complex number with a random phase factor, and is the permittivity. The coupling of light from the QW into the cavity is given by | β|, and we will for simplicity assume that | β| = 1. The presence of noise in the cavity does not significantly change results and will automatically eliminate any dual-wavelength solutions that are only marginally stable.

Numerical results and discussion
In order to identify situations where a semiconductor laser operates under stable two-wavelength emission, we have to find conditions where (i) the desired two modes experience comparable round-trip amplification and (ii) where all the other possible modes, in particular those that see more gain, are effectively suppressed.

Dual-wavelength amplification
To realize the equal amplification configuration, one is tempted to look at the quasi-equilibrium carrier gain of the active QWs. However, as it turns out, this is not very useful because of the strong non-equilibrium effects introduced by the self-consistent coupling of the propagating intracavity field with the electronic QW excitations. When a laser is switched on, we have to provide a sufficiently high carrier density to allow for QW inversion, to have the gain needed for round-trip amplification. The intracavity field that builds up due to spontaneous emissions from the inverted QWs will be constructively amplified at the resonant cavity modes and the laser output stabilizes once the intracavity field intensity saturates the gain. The intensity of each spectral peak depends on the transient dynamics during buildup from noise. Under quasi-stationary operation conditions, the carrier density and the carrier distributions are dynamically established by the balancing process of competing momentum selective carrier removal through the intracavity field and the carrier replenishment by the pump process. As a consequence, the carrier density and distribution and thus the instantaneous gain significantly deviate from the quasi-equilibrium conditions. An example for such a situation in shown in Fig. 2b). Fig. 2b) shows low-momentum kinetic holes burned in the carrier distributions, where dual-wavelength lasing photons are being extracted, and high-momentum phonon scattering induced distortions near LO-phonon resonances (see more in section 3.2.1).
Using the standard gain approach, for dual-wavelength operation to occur, the instantaneous gain for each running mode has to compensate the respective losses, thus clamping the nonequilibrium gain at these frequencies. An example for such a stable dual-wavelength situation is shown Fig. 2a) where the non-equilibrium gain (blue) has clamped to the loss line (solid black). The vertical gray lines depict the dual-wavelengths. As discussed above, Fig. 2b) depicts the corresponding non-equilibrium inversion in the individual QWs that is responsible for the round-trip amplification. Comparing the actual and the quasi-equilibrium configurations, we see that that the quasi-equilibrium gain and inversion have limited use in predicting the conditions for dual-wavelength operation.
In order to evolve the full field within the VECSEL cavity we have to solve the full coupled MSBE (1-2) numerically where we account for the optical pumping that continuously generates carriers in high-momentum states whereas the lasing emission selectively removes carriers from specific regions of low-momentum states. This carrier replenishing dynamics is governed by the carrier Coulomb and phonon scattering such that the time-resolved ultrafast QW carrier dynamics forces the timesteps on the order of a fraction of a femtosecond. On a longer timescale, the convergence of the intracavity field to a stable solution takes several thousand cavity round-trips. In contrast to mode-locked pulse solutions with 100's of femtosecond duration (e.g., see Refs. [40] and [42]), the intracavity dual-wavelength field is always interacting with the optically active QWs. This continually drives the microscopic QW carriers into dynamic non-equilibrium distributions, which are never truly stationary due to intrinsic multi-mode field dynamics e.g. in the beating with the sum-and difference-frequency field components.
The need to include the full microscopic carrier scattering via Eqs. (5)(6) to arrive at an asymptotic non-equilibrium solution of the coupled MSBE poses a significant numerical challenge. The multiple dimensional integrals contain millions of multiplications, which effectively increases the computational demand of each time step by orders of magnitude relative to scattering-rate approximations. Luckily, by taking advantage of modern parallel supercomputers, and the implementation of efficient parallelization strategies we can make headway on these problems. We run simulations on an SGI UV2000 computer with 384 cores and 4TB of shared memory, which enables us to run multiple simulations in parallel. Each simulation can take from 12 hours up to 30 days to determine the intracavity field stability.
Given the complexity of these full simulations it is imperative that we devise a robust and efficient method to a priori find promising conditions for initializing the simulation with the goal of finding potentially stable dual-wavelength emission. In order to efficiently determine the regime of potentially stable two-wavelength operation, we have devised a mapping method to initialize the system close to a stable solution with two equal spectral intensities. For this purpose, we specify a gain chip with QWs initialized at a chosen background density. For each pair of frequencies (ω 1 , ω 2 ), we inject an ideal dual-wavelength field of the form, E(t) = E in (e −iω 1 t + e −iω 2 t ), until the asymptotic carrier distributions stabilize, i.e., the gain chip gives constant spectral amplification. At this point, the injection amplification of both frequencies are recorded, which are not necessarily equal. We then repeat this process for multiple input field amplitudes, E in , until an amplitude that results in equal spectral amplification of both frequencies is found. However, we note that there are frequency pairs where no field amplitude will give equal spectral amplification of modes. Fig. 3. Overview of asymptotic dual-wavelength solutions (color surface) for frequencies (ω 1 , ω 2 ) relative to peak gain computed with the injection mapping scheme for a specific gain chip. The surface amplitude gives the non-equilibrium round-trip amplification that a dual-wavelength solution with the given frequencies would experience from the gain chip. The solid red line indicates the two intersections of a constant loss and the background round-trip amplification, shown as an inset. Single frequency operation (ω 1 = ω 2 ) is indicated with a dashed black line.
The resulting dual-wavelength field with equal amplification for each pair of frequencies is then an idealized approximation of the real situation because there is no feedback, spectral broadening, phase difference between the peaks, nonlinear wave mixing, or frequency drift included in the solution. To obtain a self-consistent solution of the full problem, we use these approximate solutions as input fields. In this way, we can very efficiently limit the parameter space where we find stable dual-wavelength output and the reduced computational time allows us to study the full non-equilibrium carrier dynamics. Fig. 3, summarizes the outcome of applying the mapping technique for various choices of injected dual-wavelengths. The color surface shows the amplification that potentially stable dualwavelength solutions experience during a single pass through the gain chip. Points outside of this surface does not have a dual-wavelength solution with equal amplification of the spectral peaks. However, this would not rule out multiple families of stable unequal amplitude dual-wavelength solutions.
The red line in Fig. 3 is obtained by recording the intersections of the background amplification with a constant loss. The background amplification, shown as an inset in Fig. 3, is the undisturbed QW gain modulated by the DBR and other material layers. This red line represents the theoretical limit of dual-wavelength operation where the QW carriers would be equal to the background, i.e., undisturbed pump density. Indeed, the background round-trip amplification approximates the amplification of dual-wavelength solutions along the nearby surface edge.

Dual-wavelength stabilization
Classically, gain clamping to the cavity loss will dictate single-wavelength output from a laser. Hence, in order to establish dual-wavelength operation, it is necessary to suppress emission at the peak gain wavelength by restricting the possible lasing modes. This suppression can be accomplished using either frequency selective gratings, etalons, or a sufficiently short optical cavity with properly placed longitudinal modes.
The first method involves a carefully selected external cavity length, which will ensure that the longitudinal cavity modes overlap with the two chosen wavelengths and no other mode experiences higher gain. The possible cavity lengths are usually short (less than a hundred micrometers) and result in strict spectral filtering because only the resonant cavity modes will be amplified. The second method involves inserting an additional etalon inside the cavity. In the third method, a finite impulse response (FIR) frequency filter is inserted into the cavity.

Short cavity
To illustrate the dual-mode dynamics, we start with a very short-cavity configuration since this provides the most straightforward mode selection method. Fig. 4 gives an overview of the cavity field and microscopic dynamics for a cavity of length 31 µm resulting in frequencies at (ω 1 , ω 2 ) = (-26.7 meV, -8.8 meV), i.e., 4.3 THz spacing, in Fig. 3. Fig. 4(a) presents the stable output spectrum whereas (b) shows a snapshot of the inversion in all ten QWs (blue) during stable operation. We notice the pronounced k-dependent differences between the actual and the background inversion (black dashed line). The strongest deformations result from the k-selective carrier removal by the running modes. In addition, we also notice modifications at higher k-values away from the main spectral peaks that originate from the LO-phonon resonances in the carrier-phonon scattering sums, Eq. 6.
In Fig. 4(c), we depict the dynamic inversion change for the example of the 6 th QW away from the DBR. The dual mode intracavity field is continually extracting carriers. The dynamic mode beating results in a dynamic k-dependent inversion oscillation around an equilibrium value. Even though these oscillation amplitudes are reduced once the intracavity field has stabilized, the carrier distributions never become truly stationary for dual-wavelength output. Fig. 4(e-f) shows the time dynamics of the electron replenishing and carrier scattering. The carrier pumping relaxes the carriers to the background distribution, shown as a black dashed line in Fig. 4(b). In this particular example, hot carriers are injected into the QWs with a peak right below k = 5. Carrier scattering is equilibrating the hot high momentum carriers into the low momentum region, which the cavity field can then extract, see Fig. 4(f).

Intracavity etalon
A common method to generate dual-wavelength operation in experiments is to use an intracavity etalon. To study this scenario, we present in Fig. 5 the results of a simulation with a glass (n = 1.5) etalon of length 99.1 µm which has been placed inside an external cavity of length 3.3 mm. The etalon is configured to induce dual-wavelength output with a frequency spacing of In Fig. 5b), we see the spectral intensity of the two highest output modes competing for amplification. In particular, the black dashed circles show where an adjacent mode overtakes the momentarily dominant wavelength. The total energy in the mode families surrounding the spectral peaks at -19.4 meV and -10.9 meV are compared in Fig. 5c), by integrating over a ±4 meV region. The time dynamics clearly shows that the spectral peaks are strongly anti-correlated. Similar behavior has been observed in experiments [29,30].
In panels d-f), we present snapshots of the output spectra at 250 ns, 336 ns, and 418 ns after initialization. The three snapshots show the spectral dynamics surrounding the two main families of modes located around -19.4 meV and -10.9 meV. Starting from Fig. 5d) at 250 ns, we see that the dominant mode near -19.4 meV is overtaken by an adjacent mode around 336 ns, shown in Fig. 5e-f). These snapshots are marked in Fig. 5b) with thick vertical grey bars and show a typical situation. Indeed, there are similar events indicated by the dashed black circles in Fig. 5b). Altogether, even though we run the simulations for up to 1600 ns, i.e., for a total simulation time corresponding to about 30 days, the emission spectra never stabilize due to the ongoing competition between the different modes. This clearly emphasizes the relatively poor performance of the intracavity etalon which not only provides frequency selective filtering but also a certain amount of intra-cavity reflected light that leads to destabilization of genuine two-mode operation. However, in a variety of practical settings (e.g., intracavity DFG THz generation) a clean dual-wavelength THz signal might not be required [30].

Finite Impulse Response Filter
In experimental setups, the external cavity can include a frequency selective grating which will spatially separate the intracavity field spectrum [50][51][52]. The separated frequencies are reflected by a frequency selective mirror and returned to the cavity. This method allows for a high degree of spectral control in the filtering of the intracavity field. Methodically, this procedure first applies a Fourier transform to the intracavity field and then a spectral filter to the back reflected signal which is returned to the cavity [50]. In our simulations, an idealized version of this method can be implemented through a FIR filter with a customizable filter function. The FIR filter is placed in the external air cavity and provides a frequency dependent loss to the intracavity field through a convolution with the impulse response of a given filter function.
An example the simulation results is shown in Fig. 6 where we include a FIR frequency filter in a 3.3 mm external cavity with spectral peaks centered at -19.4 meV and -10.9 meV. Here, we have constructed the filter function as two Gaussian peaks with 1 THz full width half maximum and 2 THz peak-to-peak spacing. The resulting output spectra is shown in Fig. 6a) and it is clearly very different from the one produced by the etalon, in Fig. 5a). As we can see in Fig. 6, the FIR filter successfully stabilizes 2 THz dual-wavelength operation at (-19.4 meV, -10.9 meV). The FIR filter is also useful to realize other solutions on the surface of Fig. 3 with mode separation less than 2 THz, which are unreachable with a short external cavity.

Self-consistent dual-wavelength results
To obtain a general overview, we now scan the parameter space for optimal dual mode operation conditions. For this purpose, we perform series of numerical simulations solving the full microscopic equations, always initializing the runs utilizing the results summarized in Fig. 3. As we described earlier in Sec. 3.1, we choose a pair of frequencies, (ω 1 , ω 2 ), as well as the initial field amplitude, E in , which gives equal spectral amplification for both frequencies. At first, we introduce the QWs into a stable state and allow the approximate -not yet self-consistentdual-wavelength field to fill the cavity, before we set the output coupling loss equal to the spectral amplification. The full simulation is then started from this initial configuration and is run until the feedback from the cavity equilibrates the system into a self-consistent solution. Fig. 7. An overview of dual-wavelength solutions computed from initial conditions in Fig. 3. The black dashed line is an outline of the lower part of Fig. 3. Colored circles indicate stable dual-wavelength operation and black triangles indicate simulations that converge to single-wavelength operation. The color of the circles indicates the difference in relative spectral intensity between the two stable modes.
The asymptotic dual-wavelength amplification is far from the linear background amplification, as seen in Fig. 3. There is a broad range of potential dual-wavelength solutions with frequency separation ranging from around 1 THz to above 20 THz with multiple different field intensities. These solutions all distort the QW carriers to varying degrees depending on the field intensities and cavity field spectral locations.
Generally, the successful simulations have a relatively short equilibration phase after initialization followed by a convergence into stable dual-wavelength operation. It is during this phase that the self-consistent nonlinear solution is generated from our ideal initial guess in Fig. 3. The non-equilibrium QW dynamics will to some degree pull the initial spectra away from the cold cavity modes, nonlinear wave mixing creates spectral peaks that feed back into the cavity, and the non-equilibrium modifications of the inversion can modify or even destabilize the dual-wavelength operation. Altogether, these effects can change the relative spectral peak intensities from the ideal initial guess or even lead to destabilization of the dual-mode operation.
Using the short-cavity mode stabilization, we find multiple independent solutions on the surface of Fig. 3 with frequency spacings from about 4 THz and close to 20 THz. In this range, the calculations converge in a manner similar to Fig. 4, but with different transient dynamics and final peak intensities. For spectral peak separations below 4 THz, the longitudinal modes are too closely spaced and additional spectral filtering is required to stabilize dual-wavelength operation. In this frequency range, a FIR filter has been inserted into a longer external cavity to stabilize dual-wavelength solutions, with results similar to Fig. 6.
In Fig. 7, we present a summary of our full simulation results. In particular, we find that some mode-pairs near the red line, seen in Fig. 3, converge directly to single-wavelength operation. This is because of nonlinear influences from the QWs that pull the wavelengths away from the initial guess and into a region with only single-mode operation. Furthermore, we find that the individual spectral peaks are not necessarily equally strong. In particular, we note that the deviations from our initial guesses tend to increase as we approach the single-wavelength region. Nevertheless, the our initialization method turns out to be very useful for reducing the parameter space of dual-wavelength solutions and reduces simulation time for very time consuming calculations.

Conclusion
In summary, we have systematically investigated dual-wavelength operation in VECSELs using a microscopic many-body model with carrier scattering computed on the level of second Born-Markov approximation coupled to Maxwell's equation. In order to study dual-wavelength operation, an injection mapping scheme is used to approximate the asymptotic dual-wavelength operation. This approximation is used to identify ranges of stable dual-wavelength operation in phase space and we implement a few different methods to realize select dual-wavelength output.
Our numerical results show that dual-wavelength operation in VECSELs produces strong nonequilibrium QW carrier dynamics. The carriers are driven away from quasi-Fermi equilibrium dynamics and carrier scattering produces pronounced distortions in microscopic QW inversions. Even under stable dual-mode operation conditions, the carriers are not stationary but oscillate at the beating frequency.
Frequency selection via an intracavity etalon results in strong anti-correlated output intensity noise. The etalon does not produce clean dual-wavelength operation, but establishes two dominant families of modes. These modes compete for energy from energetically nearby carriers, which in turn create anti-correlated intensity fluctuations. Finally, we show that a promising method to stabilize the dual-wavelength output can be realized using an ideal frequency filter.

Funding
This material is based upon work supported by the Air Force Office of Scientific Research under award number FA9550-17-1-0246.