Calculation of the steady states in dynamic semiconductor laser models

We discuss numerical challenges in calculating stable and unstable steady states of widely used dynamic semiconductor laser models. Knowledge of these states is valuable when analyzing laser dynamics and different properties of the lasing states. The example simulations and analysis mainly rely on 1(time)+1(space)-dimensional traveling-wave models, where the steady state defining conditions are formulated as a system of nonlinear algebraic equations. The performed steady state calculations reveal limitations of the Lang-Kobayashi model, explain nontrivial bias threshold relations in lasers with several electrical contacts, or predict and explain transient dynamics when simulating such lasers.


Introduction
Semiconductor lasers (SLs) and coupled SL systems are beneficial in many modern applications requiring specific characteristics of dynamic or stationary emission. A variety of models are used for simulations of well-above-threshold dynamics of SLs. Advanced models, defined by the complex systems of 1(time)+3(space)-dimensional (1+3-D) or 1+2-D PDEs (Hess and Kuhn 1996;Inoue et al. 2019), can give a deep insight into the spatiotemporal dynamics. These models, however, typically rely on a large number of not very well-known parameters, require advanced numerical tools, are computationally expensive, and allow only limited analysis. On the other hand, in this work considered 1+1-D PDE (Javaloyes and Balle 2009;Radziunas 2017) and even simpler DDE/ODE models (Lang and Kobayashi 1980;Yamada 1989;Danckaert et al. 2002) may lack quantitative precision but can be quickly solved on standard computers and admit a variety of analytic and semianalytic methods for their analysis.
Calculating and analyzing stable and unstable steady states is crucial for understanding SL dynamics (Schatz 1992;Sieber 2002;Bauer et al. 2004;, estimating different lasing characteristics, or designing SLs for specific real-world E parts ("sections") composed SL device cavity, coupled to the real carrier density N in m a "active" sections S a of this cavity, see, e.g., pink-colored boxes in Fig. 1. Within the remaining m p = m − m a "passive" sections S p , the carriers are not present (e.g., air gaps between different coupled lasers, light-blue boxes in Fig. 1), or the field and carrier rate equations are decoupled (passive waveguides, yellow boxes in the same figure). and N are vector functions of time and, in the PDE model case, space. Typically, dynamic SL models are given by the systems of equations governing the evolution of fast fields and relatively slow carriers N: Vector P(| | 2 ) , defined by a set of squared moduli of field vector components, represents the local or global photon number or power within the whole optical field, separate optical modes, or different polarization components. Typically, P(| | 2 ) enters functions and N with small (nonlinear gain compression) prefactors and is linear with respect to its arguments. Complex propagation factor consists of the static and dynamically varying parts 0 and . 0 in all m laser parts is defined by field losses and built-in or temperatureinduced nonuniformities of the refractive index. In the m a active sections determined is mainly defined by the carrier-dependent refractive index change and gain (vector) functions n R (N) and g(N), which can be extended by field power-depending nonlinear corrections. F sp is a Langevin spontaneous emission term. Since the carrier dynamic is usually slow, and the gain compression impact is small, changes of the propagation factor within short and even moderate time intervals are also small. Thus, in more complex models, freezing = ⋆ (obtained during time integration of the model equations at the fixed time instant, for example) and analyzing by operator H( ⋆ ) determined linear with respect to field equation (1), one can extract interesting information on the structure of the optical fields, which remains almost preserved within these short or moderate time intervals (Radziunas and Wünsche 2005).
In some ODE models (standard or multimode rate equations (Yamada 1989;Danckaert et al. 2002)), instead of complex vector the real photon number function is used. N and H , in this case, are scalar-, vector-, or matrix-functions of carrier and photon numbers. In DDE models, (t) is defined by a single or several ( s ′ ) complex components, while m ′ -component real vector-function N(t) represents possibly multilevel carrier dynamics in active sections. N is a real vector-function, whereas H , besides the simple scalar-or matrix-type functional part, also includes a time-delay operator D . A most prominent representative of DDE models for SLs with s � = m � = 1 is a Lang-Kobayashi (LK) system for lasers with delayed optical feedback (Lang and Kobayashi 1980), see Fig. 1c. A normalized version of this model is determined by (4) H = (N) + D , where g = N , n R = H N . H , , I are the linewidth enhancement factor, feedback rate, and bias current; is the delay time, whereas the small parameter , the ratio of photon and carrier lifetimes, indicates a slow-fast nature of the SL model. In 1+1-D PDE models, typically used for describing narrow-waveguide edge-emitting lasers, is a (vector-) function determining the field distribution along the whole SL device cavity. The field operator H accounts for the spatial derivative z along the longitudinal coordinate z and the boundary-interface-conditions (BIC), given by a set of field reflection-transmission relations at the edges of all m SL device sections, see thick vertical bars (edges) and corresponding thin red and blue arrows (transmission and reflection of the counterpropagating fields) in Fig. 1a. The vector function N can also depend on z within the active sections S a , indicating a local carrier density in this way. In simpler PDE models, one neglects spatial hole burning of carriers (Schatz 1992) and considers sectionallyuniform components of N, sectionally-averaged N and P , and, thus, section-wise constant propagation factor in Eq. (3). In this case, the carrier rate equations (2) for corresponding ODE, DDE, and PDE models are nearly identical. Different versions of the 1+1-D PDE traveling wave (TW) model (Javaloyes and Balle 2009;Radziunas 2017) are used to simulate and analyze the spatiotemporal dynamics in various multisection SLs, including those schematically represented in Fig. 1. A simplest scaled TW model (Radziunas et al. 2006) of a solitary SL (lower right scheme in Fig. 1a) for (z, t) = ( + E , − E ) T and a singlecomponent spatially-uniform N(t) is defined by where l, , 0 are scaled diode length, coupling factor of the counterpropagating fields + E and − E in the presence of Bragg grating (BG), and complex scaling-induced factor. r 0 and r l are complex field reflection factors at the laser facets, whereas the remaining parameters are the same as in the LK model (4) case.

Spectral problem and characteristic equation
The field equations (1) with neglected F sp for optically uninjected SL usually are rotationally-invariant. Let us assume that the propagation factor or, more precisely, its correction is fixed, = ⋆ , and the optical fields are given by Here is a time-independent normalized complex vector with the same s ′ components as the field function , e.g., ≡ 1 in the LK model (4). The real and imaginary parts of the complex frequency stand for the optical frequency (relative to the central frequency 0 ) and damping of the field. When squared, a nonnegative real scalar factor f represents the photon number or field intensity. Finally, is a real phase factor that can be selected arbitrarily in rotationally invariant systems. Substitution of the Ansatz (6) into the field equations (1) implies the spectral problem for the complex eigenfrequencies and eigenvectors , formulated for the fixed ⋆ . Typically, this problem implies a complex algebraic equation whose roots are eigenfrequencies of the spectral problem (7). To find the (simple) roots , we use the homotopy method, Newton's iterative procedure, and exploit our software LDSL-tool (Radziunas xxx). For each solving Eq. (8), a corresponding eigenvector is a scaled nontrivial solution of Eq. (7). Thus, can be interpreted as a vector-function of the provided ⋆ and 1 .
In DDE models, (7) is an -dependent system of s ′ linear equations w.r.t. , whereas in (8) is the determinant of this system. In 1+1-D PDE models, (7) still is a system of ODEs along the longitudinal coordinate z, possibly supplemented with several z-depending algebraic equations. For TW models (Radziunas 2017) relying on the first-order derivatives z only, this system within each laser section can be reduced to a couple of linear 1-st order ODEs for counterpropagating optical fields + E and − E , determining ⋆ -and -dependent C 0 -semigroup. At least for the piece-wise constant ⋆ (a natural approximation during the numerical treatment of model equations) these semigroups can be represented by analytic transfer matrices M j (z � , z �� ; ⋆ , ) , translating the vector (z) from z ′′ to z ′ within any device section S j . Combining sectional matrices M j with all BIC allows deriving Eq. (8) for simple linear (Radziunas and Wünsche 2005) and a variety of nontrivial (e.g., ring) SL device configurations , whereas the vector at any z can be written as an analytic function of ⋆ and : Here z ⋆ is an appropriately selected edge of one of the SL device sections, whereas the fixed vector ⋆ satisfies the BIC at this edge and imposes the same scaling for all vector functions (z) . For example, z ⋆ and ⋆ in linear SL cavities can be the front facet of the SL device and the fixed vector satisfying reflecting boundary conditions at this facet, respectively (Radziunas and Wünsche 2005). For more complex cavities, the selection of z ⋆ and ⋆ is less trivial and should be made with care.
We also note that T defined in Eq. (9) is an analytic function of arbitrary ⋆ and . Only those vector functions defined by ( ⋆ , ) satisfying the characteristic Eq. (8) can fulfill all BIC imposed on the optical fields.

Steady state equations
Each steady state of the dynamic SL model (1-3) is the solution having the form Here ̃ and Ñ are time-independent complex and real vector functions, f is a nonnegative real constant, whereas is an arbitrary phase factor. In contrast to typical stationary states of the dynamical systems, complex optical fields of the steady states rotate in time with a constant real relative optical frequency ̃.
For DDE models, is a preselected vector with s ′ components determining the scaling of . Thus, P in (12) is known, and the steady state defining system (11) consists of one complex characteristic equation and m ′ real equations, relating two real constants f 2 , and m ′ -component real vector N. The roots of this nonlinear system define the steady states of the DDE models. In LK-type models, these steady states are best known as external cavity modes, ECMs.
In the case of 1+1-D PDE models (or, more precisely, the TW model), P(z) and N(z) are real vector functions, representing local photon and, possibly, multilevel carrier densities in m a active sections of the SL device. Eqs. (11,12) is a system of algebraic and functional (z-dependent) equations, which can hardly be resolved on the functional level. Thus, we treat it on the discrete level, induced, e.g., by the numerical method during model simulations with LDSL-tool (Radziunas xxx). Namely, we assume that • all active sections are subdivided into m D small subintervals s a = [z � a , z �� a ]; • N(z) is constant within each s a and is fully defined by m ′ D real numbers; • for any z from s a spatial average ⟨P As a consequence, ⋆ = (N, f 2 P) and N are piece-wise constant w.r.t. z, and evaluation of T(z, ⋆ , ) from Eq. (9) for z = z � a and z ′′ a is sufficient for finding an exact value or numerical approximation of P(⟨� � 2 ⟩ s a ) . Thus, P can be expressed as a function of ⋆ (i.e., N and f 2 P ) and , whereas the steady-state defining system (11,12) can be rewritten as a system of nonlinear algebraic equations Here, one complex equation = 0 (equivalent to two real equations), m ′ D real equations N = 0 , and m ′′ D real equations involving P relate two real constants f 2 , , and m � D + m �� D real constants defining real vector-functions N and P. Once the system (13) is solved, the reconstruction of ̃ in Eq. (10) is performed using transfer matrices and function T (9). To solve numerically (Radziunas xxx) this system, we use the homotopy method, Newton's iterations, and exploit analytic expressions for partial derivatives of all algebraic functions in Eqs. (13). Figure 2 represents calculations of ten steady states in monolithically integrated external cavity (EC) diode laser (DL), schematically represented in Fig. 1b. The TW model and the basic set of parameters within the active and dispersive reflector sections are identical to those used for simulations of Fig. 4 in Ref. (Wenzel et al. 2021), except for in the present case used linear gain and refractive index functions, related by n R (N) = H g(N) (Radziunas 2017). In contrast to Ref. (Wenzel et al. 2021), here we used an additional 2mm-long passive waveguide section, which induced a coexistence of multiple stable steady states. We assume a single carrier level (i.e., m � D = m D ) and account only for a combined photon density in both counterpropagating fields, For m D = 100 subintervals in the single 1mmlong active section, we have 202 real algebraic equations in (13) relating the same number of variables defining the steady states of the (discretized) TW model.
Curves of different colors in panels (a) and (b) of Fig. 2 show distributions of the carrier density Ñ and field intensity f 2P within the active section for ten calculated steady states. The reconstructed counterpropagating field intensities |f̃+ E | 2 and |f̃+ E | 2 along the whole cavity of two stable steady states are given in panel (c). Panel (d) shows the intensity of the complex EC reflection spectrum R( ) , which is the response of the passive reflector, consisting of the passive waveguide and BG, to the incident monochromatic wave e i t . A similar (instantaneous) response G( , ) of the active part of the ECDL (including possible reflection and transmission of the fields at the interface to the EC) depends additionally on the instant value of the factor , see Ref.  and schematics in Fig. 1b and (c). Finally, large empty bullets in Fig. 2e show the relative frequency ̃ and corresponding threshold gain 1+f 2P ⟩ S a of these states. Only two of the calculated steady states with the lowest gain values were stable in the considered case, which was confirmed by transient simulations of the dynamic TW model.
The algorithm for the location of the steady states in the general TW model was developed during the recent study of SL emission's linewidth (Wenzel et al. 2021) in ECDLs Krüger et al. 2019). In many applications, however, simplified TW with sectionally-averaged N, N , and P (where a set of m D subintervals s a is reduced to m a sections S a ) or even more simple DDE (such as LK-) models are  N and (N, ) , respectively), the vector variable P is defined by (N, )-dependent function P and can be eliminated from Eq. (13). The steady-state defining system is reduced to one complex -and m a real N -equations, together defining real m a -component vector Ñ and real factors ̃ and f 2 .
Study of these low dimensional algebraic systems, including calculation of the steady states and their continuation with the change of model parameters, can be performed using standard numerical path-following techniques. For example, small black bullets in Fig. 2e are the steady states of the corresponding basic TW model for considered ECDL. Dashed and dotted black curves in the same figure are branches of these states calculated for tuned field phase shift and losses in the passive waveguide of the device, respectively. Even though the basic model neglects the nonuniformity of carrier densities (see Fig. 2a) and the gain compression, in the considered case it was perfectly suited for the identification of frequency and gain of the steady states of more general TW model (cf. large colored and small black bullets in Fig. 2e). In the remainder of this work, we shall exploit the basic TW model for location and analysis of the steady states in SLs with one or two active sections.

SL device with a single active section
Once the basic TW model describes dynamics in SLs with a single active section ( m a = 1 ), the steady state-defining system (13) can be written as where , , and I are model parameters actively used for location of the steady states. is a function of only two real unknowns and N, and each steady state is fully defined by a real-valued triple (̃,f 2 ,Ñ) . By resolving complex characteristic equation in (14) we find a set of pairs (̃,Ñ) . Substitution of each such pair into the remaining equation allows finding the corresponding f 2 , which can be tuned by changing bias current parameter I. Only nonnegative f 2 can represent steady states. Bias current I th , at which f 2 is zero, is a threshold current of this state.
For a variety of SLs (such as the ECDL from which for any pair ( , N) uniquely define the feedback factors and . The first and the second equations are independent on and , such that they determine two-dimensional manifolds in three-dimensional (N, , ) -or (N, , )-space, respectively. The intersection of these manifolds at fixed feedback level (FFL) or fixed feedback phase (FFP) (modulus 2 ) define curves in ( , N) plane. The FFL curves for several values of are shown by solid colored lines in Fig. 3. Variables and N on the axes are nonscaled frequency and carrier density relative to those of one of the solitary SL resonances. For small and moderate , we see closed FFL loops around each resonance of the solitary laser. At the critical feedback level c , which is the field amplitude reflection �r c � = √ 0.05 at the diode edge facing the EC, the two-dimensional manifold defined by the first equation in (16) has multiple saddle points resulting the collision of the adjacent loops at infinitely large N . For even larger feedback, a single periodically modulated FFL curve is formed. For comparison, thin dashed curves in Fig. 3 show corresponding ellipses of ECMs in the properly normalized Lang-Kobayashi model (Radziunas et al. 2006). Obviously, the standard LK model can only account for a single solitary laser resonance. One can see that closed loops of two models around the selected resonance are in good agreement for small feedback levels but rapidly diverge when > 0.1 . Behind c , the ECM ellipses and FFL curves have different topologies.
The intersection of multiple fixed feedback phase (FFP) curves determined by the second equation in (16) with FFL curves for given and define all steady states of the model. A subset of such states calculated for the TW model with included material gain dispersion (Radziunas 2017) is shown by black dots on the orange = 0.1 level curves in Fig. 3b. All but one FFP curve in this diagram either bypass the selected FFL, or cross it twice, defining node-type (lower crossing) and saddle-type (upper crossing) steady states ("modes" and "antimodes" in LK-model language), or touch the FFL curve, indicating the saddle-node bifurcation positions (dash-dotted curve in the same diagram). The remaining single FFP curve (solid black line in this diagram) terminates at the solitary laser resonance, where vanishes. It is noteworthy that the LK and TW models in the close vicinity of the solitary laser resonance in the considered case with = 4.5 ns have nearly the same minimal feedback levels enabling the saddle-node bifurcation ( ≈ 1.063 ⋅ 10 −4 and ≈ 1.064 ⋅ 10 −4 ) and the feedback level ( ≈ 3.87 ⋅ 10 −4 and ≈ 3.86 ⋅ 10 −4 ) admitting a steady state at the mode degeneracy (Sieber 2002), also known as the exceptional point. Due to the simplicity of the field equations in the considered laser, nearly all curves in Fig. 3 have analytic parametric representations. In more general cases, given, e.g., by a distributed feedback (DFB) laser with a passive dispersive reflector (Radziunas 2017), or an ECDL considered above, see black curves and bullets in Fig. 3e, we still exploit the splitting of the resonance equation (16) ( e i = e 2 0,ph l ph with l ph denoting the length of the middle phase-tuning section in this case) but rely on more involved formulas and numerical path-following techniques (Radziunas xxx).

SL device with two active sections
Let us consider the basic TW model, describing the dynamics of SLs with two active sections ( m a = 2 ). An example of such a device is the ECDL, see Fig. 1b, with the active Bragg grating section. The parameter set exploited below was already used to calculate the steady states in Fig. 2e. The missing material parameters of now active BG are identical to those of the amplifying section. The steady-state defining system (13) in this case can be rewritten as Besides relating elements of the steady state-defining quadruplet (f 2 ,̃,Ñ 1 ,Ñ 2 ) , the first, second, and third equations of this system depend on the bias currents I 1 , I 2 in both active sections and the field phase tuning factor ∝ ℑ 0,ph in the middle passive waveguide section, respectively. These parameters are explored when constructing the steady-state branches and analyzing the dynamics of the whole TW model. Being linear with respect to f 2 , I 1 , and I 2 in the basic TW model, the first two equations in (17) can be resolved with respect to these factors: The complex characteristic equation does not depend on f 2 , and for considered ECDL can be split into two real equations, similar to those of Eq. (16): where G 1 and G 2 are responses of both active sections to through the internal section side injected monochromatic field. Like in the previously considered case, the first of these new equations does not depend on and can be used for the definition of 2-D steady-state manifolds (SSM) in ( , N 1 , N 2 )-space, see the semitransparent surface in Fig. 4a. The second real equation for each regular point of this surface uniquely determines the phase . For fixed , the characteristic equation defines multiple mode threshold density (MTD) (17) N 1 (N 1 , N 2 , f 2 , ;I 1 ) = 0, N 2 (N 1 , N 2 , f 2 , ;I 2 ) = 0, ( , N 1 , N 2 ; ) = 0.
curves on this surface, each determined by pairs (N 1 , N 2 ) at which the imaginary part of the complex frequency of optical mode vanishes, ≡ ℜ = . These MTD curves, given by solid black lines in Fig. 4a, can be found by standard path-following algorithms, provided one or several points on each such curve are known 2 . Moreover, MTD curves are essential in explaining the slow-fast SL model dynamics. In the case of ultimately slow carrier dynamics, several segments of these curves define low-dimensional exponentially attracting center manifolds containing all regular and irregular attractors of the considered dynamical system (Sieber 2002). The carrier densities (N 1 , N 2 ) on these specially chosen MTD curve segments imply vanishing ℑ of the corresponding mode and damping ( ℑ > 0 ) of all remaining modes. In the present example, the center manifolds are defined by the "outer" (lower-left) border of the (N 1 , N 2 )-plane projection of the MTD curve set; see thick MTD curve segments for modes I, II, and III in Fig. 4b. For (N 1 , N 2 ) belonging to the remaining "inner" segments of the MTD branches, at least one mode has negative damping, ℑ < 0 . Since the carrier dynamic of our system is relatively slow, we expect that (N 1 , N 2 )-plane projections of calculated transients should be accumulated around the outer borders of such MTD-curve representation. Each triplet (̃,Ñ 1 ,Ñ 2 ) on the 2-D SSM and on MTD curves in particular, together with an arbitrarily chosen non-negative factor f 2 , defines a steady state of the basic TW model with properly chosen parameters I 1 , I 2 , and . Whereas is determined by the last equation in (19) and the triplet alone (e.g., = 0 at black solid MTD curves in Fig. 4a), bias currents I 1 and I 2 are defined by the quadruplets and the functions J 1 and J 2 from Eq. (18). For vanishing f 2 , these are pairs of current thresholds, allowing to reach the threshold densities (Ñ 1 ,Ñ 2 ) of the selected mode but not big enough to excite this mode. A corresponding lasing threshold representation of MTD curves is shown in Fig. 4c. The laser is off for small injections, which are left and beneath all curves in this diagram. By increasing one or both currents, we can reach and cross one of the outer lasing threshold curves, which causes switching on the lasing. In our example, depending on the relation of the currents, we can switch on the lasing by crossing the branches of modes I, II, or III, see thick black  Fig. 4c, which become excited in the first part of the switching-on process. A switching-on of the laser by setting (I 1 , I 2 ) = (20, 2) mA (which is just slightly above threshold, see Fig. 4c) is represented by the lowest yellow-orange-brown curve in Fig. 4b. This curve is the projection of the calculated 40-ns transient trajectory onto (N 1 , N 2 )-plane. At the initial stage, the system reaches the threshold of Mode I and, shortly afterward, Modes IV and II, such that all these modes are excited. Initially excited Mode I dominates, and the calculated trajectory follows its MTD curve, relaxing to by the bias-currents determined stable steady-state on it (green bullet in Fig. 4b).
In general, each triplet on the 2-D SSM alone defines the steady states for multiple sets of currents I 1 and I 2 , related to each other by a real nonlinear equation F 1 = F 2 from (18). This equation is fulfilled along 1-D curves on the 2-D SSM for fixed (but not arbitrary) pairs (I 1 , I 2 ) . Green, red, and blue thin solid curves in Fig. 4a are examples of these by the phase -parametrized steady state branches calculated for three different sets (I 1 , I 2 ) . Solid bullets at the intersections of these thin curves with the solid thick black MTD curves represent the steady states of the TW model for a corresponding couple of bias currents and = 0 . Projection of all these steady state branches and MTD curves onto the bottom ( , N 1 )-plane implies the steady state representations reminding those of Fig. 2e. Another projection of the steady states and the MTD curves onto the (N 1 , N 2 )-plane is shown in Fig. 4b.
For fixed , a monotonic change of a single bias current typically causes a gradual shift of the steady-state triplets along the MTD curves but can also imply the creation or annihilation of the steady state pair on the MTD curve in the saddle-node bifurcation. Tuning of the bias current can also lead to the Hopf-bifurcations of the steady states, which are particularly expected close to the intersection of MTD curve projections in (N 1 , N 2 )-plane, see Fig. 4b. Separation of frequencies ̃ of the involved MTD branches defines the periodicity of the generated orbit. Such mode-beating pulsations available for large bias injection regions are of practical interest and were studied, for example, in Ref. (Bauer et al. 2004). In the present example, starting from the Mode I-defined stable steady-state at (20, 2) mA currents, green bullet in Fig. 4b, we alternated a slight step-wise increase of I 1 with transient simulations. In the beginning, we could observe a slow upwards shift of the still stable steady state along the MTD branch of Mode I in Fig. 4b. At I 1 ≈ 24.6mA, already slightly behind the crossing of MTD branches of Modes I and II, the steady state underwent Hopf bifurcation and could not be traced along the MTD curve of Mode I anymore. Instead, a newly stabilized steady state on the MTD branch of Mode II could be observed and traced.
A middle yellow-orange-brown curve in Fig. 4b represents a calculated 3-ns long trajectory after a large-step switch of the bias currents from (20, 2) mA (green bullets) to (50, 2) mA (red bullets). At the initial stage after the bias switch, nearly all field power is contained in Mode I, and the trajectory follows the MTD branch of Mode I, trying to converge to the steady state (red bullet) on this branch. This state is unstable, Modes II and III start to contribute to the total field power, and the trajectory first turns towards the steady state on the MTD branch of Mode III and, later, of Mode II, exhibiting mode-beating oscillations on the way.
A next switch of the bias currents to (50, 50)mA, see the upper yellow-orange-brown curve in Fig. 4b, however, was not leading to the stable steady-state (blue bullet) on the outer MTD branch of Mode III, but to the "inner" state still determined by the same Mode II, even though Mode III at the achieved carrier densities has negative damping, ℑ < 0 . In this case, Hopf bifurcation by a further small-step increase of I 2 could not be found, and two stable steady-states determined by Modes II and III were coexisting. We explain this contradiction to our original assumption of the close vicinity of the attractors to outer MTD branches in (N 1 , N 2 )-plane projections by insufficient slowness of the carrier dynamics in the considered system. Indeed, halving the spontaneous Page 13 of 14 121 recombination parameters (doubling carrier lifetime) and bias currents does not affect the MTD branches or the steady state locations but destabilizes the previously stable "inner" state determined by Mode II. A sequence of transient simulations combined with a small-step increase of I 2 starting from the Mode II-defined stable steady-state, the red bullet in Fig. 4b, has confirmed the existence of Hopf bifurcation before reaching the blue-bullet-defined position on the MTD branch of Mode II.

Conclusions
In conclusion, we discussed algorithms and challenges in calculating (stable and unstable) steady states in the TW model of the SLs. Approximation of the continuous functions by the piece-wise constant ones allows replacing the steady-state-defining functional equalities with a system of algebraic equations, resolvable using the homotopy method and Newton's iterative procedure. We have also shown that the steady states of the basic simplified TW model relying on the sectional averages of the carrier and photon densities can provide good approximations of the states in more general models. We have shown how stable and unstable steady states of the basic TW model governing the dynamics of the SL device with m a active sections can be defined by m a + 2 real nonlinear equations relating the same number of real steady state defining variables. We compared the steady states of the TW and LK models in the numerical example of the FP laser with optical feedback (case of m a = 1 ). Constructed using a first-order approximation of the characteristic equation ( , ) at the solitary SL resonance (Radziunas et al. 2006), the LK model is in good agreement with the TW model when feedback is small, and a single-mode operation of the solitary SL is pronounced. For multimode SLs or SLs with a feedback ratio exceeding the field amplitude reflection at the diode's facet facing EC, this good agreement is lost, and the usage of the LK model becomes questionable. For the ECDL with active Bragg grating section (case of m a = 2 ), we constructed the mode threshold density branches. These branches can be used to define all steady states of the system and derive nontrivial relations of bias currents at both active sections needed to switch on the laser device. Several transient simulations have shown how calculated trajectories approach and follow predominantly outer MTD branches in (N 1 , N 2 ) projection plane. The existence of the attractors on the inner MTD branches and, thus, multistability are possible due to the insufficient slowness of the carrier dynamics. Finally, we note that steady-state calculations do not provide direct information about their stability and other properties, such as the linewidth or small signal modulation response, but they still are extremely useful for the understanding of laser dynamics.
Funding Not applicable. Open Access funding enabled and organized by Projekt DEAL.
Data availibility can be requested by sending an email to the author. Note that the software kit LDSLtool (Radziunas xxx) explored in this work is not an open-source code.

Declarations
Conflict of interest there are no known competing financial interests or personal relationships that could influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.