## Abstract

Multimode dynamics in bidirectional laser cavities can be accurately described by folding space into time delay. This results in a set of delayed algebraic equations that preserve the dynamics of all cavity modes while drastically reducing number of degrees of freedom. This reduction allows for both linear stability analysis and bifurcation diagram reconstruction, as well as integration times reduced by orders of magnitude.

© 2012 Optical Society of America

A remarkably rich area in laser physics is that of multi-longitudinal mode dynamics, which result from the nonlinear interaction of the electromagnetic field with the active material. Multimode regimes are more prominent in systems where the gain spectrum is broad as compared to the longitudinal mode spacing and they usually involve widely different time scales. Depending on the particular system considered they can manifest as mode-partition phenomena [1], mode hopping [2], mode itinerancy [3], or mode-locking [4]. Mode-locked operation has important technological applications since it allows to generate short optical pulses at high repetition rates [5]. Yet, our understanding of these multi-mode regimes is limited due to the inherent difficulties that exist for studying them with the standard tools of theoretical nonlinear dynamics such as Linear Stability Analysis (LSA) or bifurcation diagrams using software like [6, 7].

Multimode dynamics are naturally described in terms of a bidirectional Travelling Wave Model (TWM) [8] that requires considering a huge number of Degrees of Freedom (DOF) which impedes LSA and bifurcation studies, hence limiting TWMs to numerical analyses [9–14]. The reason is that to properly account for the broad gain spectrum requires an appropriately small time step *δt*, while the Courant-Friedrichs-Lewy (CFL) condition for stability imposes an accordingly fine spatial discretization *δz*.

This has stimulated the search for more tractable models that allow for analytical or semi-analytical study, but these have been developed only in quite specific situations. For an almost conservative cavity, the Uniform Field Limit (UFL) [15] can be invoked, which permits to derive a system of coupled ordinary differential equations (ODE)s —the so-called multimode Rate Equations (RE); however, which modes participate in the dynamics and how they are coupled still require phenomenological considerations.

When the length of the active material is small as compared with the optical wavelength, i.e. *L _{a}* ≪

*λ*, the dynamics can be described by a set of ODEs complemented by Delayed-Differential Equations (DDEs) or Delayed-Algebraic Equations (DAEs) [16, 17] that represent propagation in passive sections, as e.g. an external cavity; when coupled to a RE approximation for the laser dynamics, this leads to the well-known Lang and Kobayashi DDE model for lasers within the weak optical feedback approximation [18].

In the case of an extended active medium with *L _{a}* ≫

*λ*, a tractable model has been developed only for unidirectional cavities under strong approximations; linear gain and recombination, no internal losses and the neglect of gain dispersion due to the presence of a external frequency filter [19]. The resulting DDE model has been applied to the study of semiconductor passive-mode-locking (PML) with an intracavity saturable absorber (SA) section, and it has recently been extended by including internal losses in a lumped approximation as well as the carrier dynamics of multi-level quantum-dot materials [20].

In this letter we present a DAE description for bi-directional laser cavities containing an extended medium that permits to drastically reduce the number of DOFs as compared to a TWM while accurately preserving the dynamics. This large reduction of the number of DOF allows for a direct linear stability analysis (LSA) and bifurcation diagram reconstruction and improves time integration speed by two orders of magnitude; this enables extensive stochastic analysis required to evaluate e.g. jitter in a pulse train. The method, based on folding space into time delay, comprises all of the aforementioned alternative approaches and it can also be adapted to convective problems in other fields.

The equations in a bidirectional TWM for the forward and backward amplitudes *E*_{±} read

*z,t*) are normalized to the cavity length

*L*and to the cavity time of flight

*τ*=

_{c}*L/υ*, respectively. The source terms contain both the internal losses

_{g}*α*and the active medium polarization

*P*

_{±}(

*z,t*) whose dynamical evolution is given by a local equation determined by the type of active material. The solution of Eq. (1) reads

_{±}is the advected integral of the source. Notice that the results in [17] correspond to Eq. (2) for a cavity without internal loss and with a localized active medium. If the source terms

*S*

_{±}do not strongly vary along the characteristic lines

*ξ*

_{±}=

*z*∓

*t*in the integration interval, one can approximate

This spatio-temporal DAE —together with the corresponding boundary conditions and the evolution equations for the active material— provide an alternative formulation of bidirectional multimode laser dynamics that generalizes previous approaches for extended media:

- In the case of a quasi-monochromatic field and under the UFL, it is natural to expand
*E*(*z,t*−*τ*) =*E*(*z,t*) −*τĖ*(*z,t*) + ··· which leads to the RE approximation, whose particular form depends on both the boundary conditions and the material description.

When Eq. (4) is applied to cavities where the UFL holds, the correct information for all cavity modes is preserved even if only the endpoints of the cavity, *z* = 0 and *z* = 1 are considered, i.e. *τ* = 1. The reason is that the intermediate “missing” spatial points have been folded into the time delays. As an example, we consider a FP cavity containing an active medium composed of homogeneously broadened two-level atoms, whose evolution equations are those given in [10]. The boundary conditions for the fields are *E*_{+}(0,*t*) = *R*_{1}*E*_{−}(0*,t*) and *E*_{−}(1*,t*) = *R*_{2}*E*_{+}(0*,t*) where *R*_{1,2} denote the complex reflection coefficients of the facets. The LSA of the “off” solution *E*_{±} = 0 yields the threshold current *J _{m}* and lasing frequency

*ω*of the cavity modes through the implicit equation

_{m}*g*(

*ω*) =

*γJ/*(

*γ*

*− iω*) −

*α*. The solutions for different values of

*R*

_{1}=

*R*

_{2}=

*R*are plotted in Fig. 1 (dots) together with the exact result (lines) in the TWM [10], given by Clearly, Eq. (5) corresponds to the [1/1] Padé approximant of Eq. (6) given by

*e*≈ (1 +

^{g}*g*/2)/(1−

*g*/2), which is accurate provided that

*g*< 0.5 [21]. Indeed Eq. (5) gives an excellent approximation to

*all*cavity modes for high reflectivities

*R*≳ 0.75, getting worse as

*R*reduces, the net threhsold gain increases and the field amplitudes become non uniform in space. We emphasize that in obtaining Fig. 1 we have considered only the endpoints of the cavity,

*z*= 0 and

*z*= 1, but we recovered the correct information for all cavity modes in the UFL.

When the UFL does not hold, the same idea can still be applied by dividing the cavity into *N* sections short enough for the UFL to apply within each of them. The fields in section *j* ∈ [1..*N*] are then coupled to their nearest neighbors by a temporal DAE of the form

This is at variance with a TWM approach where the temporal and spatial increments are defined one from the other according to the CFL condition. This allows to largely reduce the number of DOF in the DAE approach as compared to an equivalent TWM and it also opens the way to non uniform spatial grids. In the general case of a local material response given by a convolution integral [11] that involves past values of the local fields *E*_{±} (*z*,·) over a temporal extent related to the decoherence time of the polarization. The stencils of the TWM and the DAE approach are represented in Fig. (2), with *N _{TWM}* and

*N*spatial points, respectively. Accordingly, the reduction in DOF is given by

_{DAE}*d*=

*N*. This reduction can be different within each section of a compound devices as exemplified in Fig. (2)b where

_{TWM}/N_{DAE}*d*= 8 and

*d*= 2 in the first and second section, respectively. Finally, notice that the DAE approach contains the TWM as a particular case.

The most obvious advantage of the DAE approach is that it substantially shortens computation time as compared to TWM. This is illustrated with a most demanding example, namely a two section PML cleaved-facet FP semiconductor laser. The evolution equations for the active medium and parameter values for III-V materials are given in [12], and the simulations are performed with Freetwm [22]. The total length is *L* = 1mm, with a Saturable Absorber length of 3.2%, and the roundtrip is *τ _{c}* =

*L/υ*= 12.5ps. In order to correctly describe the polarization of the active medium, the time step is

_{g}*δt*= 23.5fs. The CFL condition thus imposes

*δz*=

*υ*= 1.89

_{g}δt*μ*m for the full TWM, hence the gain and SA sections are discretized over

*N*

_{1}= 513 and

*N*

_{2}= 17 points, respectively. Figure 3 depicts the bifurcation diagrams for the output intensity (pulse height, duration and skewness) as the current in the gain section is scanned in the full TWM (left column) and for two different DAE descriptions (center and right columns). For the full TWM we observe that immediately after threshold, the laser emits in pulsed mode. The pulse train is unstable due to a Q-switching instability which remains up to

*J*≈ 3

*J*; above this value, the pulse train becomes stable and regular. As the current increases, the pulses shorten down to ≃ 1 ps, and at the same time they become more asymmetrical. The results in the second and third columns —obtained with our DAE approach with

_{th}*N*

_{1}= 33 and

*N*

_{2}= 5, and

*N*

_{1}= 9 and

*N*

_{2}= 2, respectively— completely agree with those for the full TWM for all the variables. However, they have been obtained much more easily: in the third column, the number of DOF and the calculation time have been divided by ∼ 50 and while the temporal increment is still 23.5fs, the distance between points in the gain section is close to Δ

*z*∼ 100

*μ*m. In the three cases, the number of modes is identical (∼ 1000) as determined by 2/

*δt*.

The second great advantage of the reduced number of DOF is that it permits to apply the standard tools of non linear dynamics which are based on the LSA of the different solutions. This permits to determine the different regimes of operation and how they are connected. Different LSA strategies are implemented in [6], [7] or [23], but they all require performing the QR decomposition [21] of a characteristic matrix, a process whose complexity grows cubically with the number of DOFs. This restricts this kind of analysis to relatively small systems and makes it impossible for the full TWM. Instead, the aforementioned huge reduction of the number DOF allows us to perform the LSA of the coupled DAE system easily. As a test example, we compute the eigenvalue spectrum around the off solution for a simple FP semiconductor laser for different number of sections *N* = 9, 17 and 65. We use a method similar to the one in [23], and we remark that while the results with *N* = 9 are obtained in a few seconds, those for *N* = 65 required several hours as well as several Gigabytes of memory; the LSA of the full TWM (with *N* = 513) could not be performed. The results are shown in Fig. 4, where we plot the real part of the eigenvalues versus its corresponding imaginary part. Each eigenvalue corresponds to a laser mode, with the imaginary part giving the modal frequency and the sign of the real part determining whether the mode is stable (ℜ(*λ*) < 0) or unstable. One recognizes the typical asymmetric gain spectrum of an active semiconductor as given by the susceptibility of [24]. The insets show the eigenvalues in the regions of maximum amplification (top) and maximum absorption (bottom). The small differences that can be observed for the eigenvalues obtained with different section number arise from the residual non uniformity of the field within each section. As expected, these vanish as the number of section is increased and the obtained eigenvalues converge to their exact values. Hence, as in the previous example there is a very good agreement between the results obtained with *N* = 9 and those obtained with *N* = 65.

In summary, we have presented an alternative approach for the modelling of spatially extended laser cavities where the non linear travelling wave equations are recast into an ensemble of DAEs. The folding of space into time delay drastically reduces the number of DOF, which allows both for a direct LSA of the underlying system as well as one to two orders of magnitude improvements for direct time integration. The method here proposed should also enable the Floquet analysis of periodic regimes as required to assess the stability of a PML.

## Acknowledgments

We acknowledge financial support from the Ramón y Cajal program and the project ALAS ( TEC2009-14581-C02-01).

## References and links

**1. **R. Linke, B. Kasper, C. Burrus, I. Kaminow, J.-S. Ko, and T. Lee, “Mode power partition events in nearly single-frequency lasers,” J. Lightwave Technol. **3**, 706–712 (1985). [CrossRef]

**2. **M. Ohtsu, Y. Teramachi, Y. Otsuka, and A. Osaki, “Analyses of mode-hopping phenomena in an AlGaAs laser,” IEEE J. Quantum. Elect. **22**, 535–543 (1986). [CrossRef]

**3. **L. Furfaro, F. Pedaci, M. Giudici, X. Hachair, J. Tredicce, and S. Balle, “Mode-switching in semiconductor lasers,” IEEE J. Quantum. Elect. **40**, 1365–1376 (2004). [CrossRef]

**4. **L. E. Hargrove, R. L. Fork, and M. A. Pollack, “Locking of he-ne laser modes induced by synchronous intracavity modulation,” Appl. Phys. Lett. **5**, 4–5 (1964). [CrossRef]

**5. **H. A. Haus, “Mode-locking of lasers,” IEEE J. Sel. Top. Quantum Electron. **6**, 1173–1185 (2000). [CrossRef]

**6. **E. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, B. Sandstede, and X. Wang, “Auto97: Continuation and bifurcation software for ordinary differential equations,” (2011).

**7. **K. Engelborghs, T. Luzyanina, and G. Samaey, “Dde-biftool v. 2.00: a matlab package for bifurcation analysis of delay differential equations,” Tech. Rep., Department of Computer Science, K.U.Leuven, Belgium. (2001).

**8. **J. A. Fleck, “Emission of pulse trains by Q-switched lasers,” Phys. Rev. Lett. **21**, 131–133 (1968). [CrossRef]

**9. **J. Javaloyes and S. Balle, “Mode-locking in Fabry-Pérot lasers,” IEEE J. Quantum Electron. **46**, 1023–1030 (2010). [CrossRef]

**10. **A. Pérez-Serrano, J. Javaloyes, and S. Balle, “Bichromatic emission and multimode dynamics in bidirectional ring lasers,” Phys. Rev. A **81**, 043817 (2010). [CrossRef]

**11. **J. Javaloyes and S. Balle, “Quasiequilibrium time-domain susceptibility of semiconductor quantum wells,” Phys. Rev. A **81**, 062505 (2010). [CrossRef]

**12. **P. Stolarz, J. Javaloyes, G. Mezosi, L. Hou, C. Ironside, M. Sorel, A. Bryce, and S. Balle, “Spectral dynamical behavior in passively mode-locked semiconductor lasers,” IEEE Photon. J. **3**, 1067–1082 (2011). [CrossRef]

**13. **J. Javaloyes and S. Balle, “All-optical directional switching of bistable semiconductor ring lasers,” IEEE J. Quantum Electron **47**, 1078 –1085 (2011). [CrossRef]

**14. **A.G. Vladimirov, A.S. Pimenov, and D. Rachinskii, “Numerical Study of Dynamical Regimes in a Monolithic Passively Mode-Locked Semiconductor Laser,” IEEE J. Quantum Electron **45**, 462 –468 (2009). [CrossRef]

**15. **L. Narducci and N. B. Abraham, *Laser Physics and Laser Instabilities* (World Scientific, Singapore, 1988).

**16. **J. Mulet and S. Balle, “Mode locking dynamics in electrically-driven vertical-external-cavity surface-emitting lasers,” IEEE J. Quantum Electron. **41**, 1148–1156 (2005). [CrossRef]

**17. **L. A. Lugiato and F. Prati, “Difference differential equations for a resonator with a very thin nonlinear medium,” Phys. Rev. Lett. **104**, 233902 (2010). [CrossRef] [PubMed]

**18. **R. Lang and K. Kobayashi, “External optical feedback effects on semiconductor injection laser properties,” IEEE J. Quantum Electron. **16**, 347–355 (1980). [CrossRef]

**19. **A. G. Vladimirov and D. Turaev, “Model for passive mode locking in semiconductor lasers,” Phys. Rev. A **72**, 033808 (2005). [CrossRef]

**20. **M. Rossetti, P. Bardella, and I. Montrosset, “Modeling passive mode-locking in quantum dot lasers: A comparison between a finite-difference traveling-wave model and a delayed differential equation approach,” IEEE J. Quantum. Electron **47**, 569 –576 (2011). [CrossRef]

**21. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes: The Art of Scientific Computing* (Cambridge University Press, 2007.

**22. **J. Javaloyes and S. Balle, “Freetwm: a simulation tool for semiconductor lasers.” (2012). Available at http://nova.uib.es/ONL/Softwares/Softwares.html.

**23. **A. Pérez-Serrano, J. Javaloyes, and S. Balle, “Longitudinal mode multistability in ring and Fabry-Pérot lasers: the effect of spatial hole burning,” Opt. Express **19**, 3284–3289 (2011). [CrossRef] [PubMed]

**24. **S. Balle, “Simple analytical approximations for the gain and refractive index spectra in quantum well lasers,” Phys. Rev. A **57**, 1304–1312 (1998). [CrossRef]