Obtaining an Operating Point Solution of a Traveling Wave Laser Model

This paper presents a method of obtaining an operating point configuration for a laser model based on a traveling wave model (TWM), which can then be used in a circuit-level simulator. The method first finds an approximate distributed single-mode stationary solution, this solution is then iterated using the traveling wave equations to an accurate single-mode solution, and finally a short pre-simulation is used to add harmonic content to create a multi-mode configuration of the laser approximating its behavior at an operating point. The effectiveness of this approximation is tested by initiating transient simulations from this operating point and comparing them to the output of the model started from an off state. The stochastic variation in the operating point for a particular configuration is also well predicted. Included in the formulation are gain compression and dispersion effects, laser chirp due to variation in the effective index of the laser mode, and spontaneous emission. Finally, the use of the three-stage process of finding the operating point in a circuit-level simulator is discussed. Not only does the three-stage method provide a quick, accurate operating point for the circuit simulator, but the ability to provide an orders of magnitude faster estimate for the initial circuit-level operating point is critical to the practicality of its use in the simulator. The first stage of the three-stage method does just this.


I. INTRODUCTION
The integration of optical and electronic devices to achieve lower costs and higher functionality is attractive for a wide The associate editor coordinating the review of this manuscript and approving it for publication was Su Yan . variety of information technology applications. A key to innovation at the design level for these circuits is the availability of sophisticated computer aided design (CAD) tools such as circuit-level simulators. Work has been ongoing into the development of a wide variety of approaches for optical circuit-level simulation. This progress has involved the creation of frequency and time domain methods for both linear and nonlinear devices and circuits as well as initial work into co-simulation of electronics and optical devices [1]- [5]. The work has produced a number of commercial products [6]- [8]. The research described in [9] presents a fully integrated optical/electrical simulator within the SPICE-like framework commonly used in electrical circuit simulation.
For devices such as integrated lasers an important component of circuit-level simulators is the compact model. Compact models are physically-based, fast and flexible. They take ''inputs'' and determine ''outputs'' using just enough detail to keep them extremely fast. For example the input for a laser might be an applied voltage and the outputs would be the current flow through the device and the optical output power. VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ An important criteria for such compact models is that they have representations in both the time and frequency domains and that a direct current (DC) or steady state (SS) solution is available for determining an operating point (OP) of the circuit. When determining the operating point of a circuit, a circuit/system-level simulation tool will often require the operating point of the compact model to be evaluated many times (perhaps hundreds) as the simulator iterates the nonlinear circuit equations to convergence. Therefore the requirement for a quick determination of an operating point is crucial for effective CAD of circuits and systems. A typical laser model for a circuit simulator is the rate equation model (REM) [10]. Such a model describes the laser with a set of rate equations that predict the average laser carrier and photon densities over time. The use of a lumped model, however, ignores key aspects of the laser operation due to its distributed nature. Standalone device simulation of lasers which includes distributed effects is a mature and prolific field. A number of these models are well suited to be transformed into compact models for circuit simulation.
Two such examples are models based on transmission line representations [11] and models based on traveling waves, Traveling Wave Model (TWM) [12]- [14].
A significant amount of work (see [12]- [15] and references found in [16] (chapter 31)) has established the ability of the TWM to predict a wide range of behaviour accurately and is confirmed by experimental work. These models are formulated to perform standalone device-level transient laser simulation and analysis. However, recently a methodology for creating a compact model for a traveling-wave-based device for incorporation into an MNA-based solver has been published [17].
If such models are to be used in a circuit-level simulator the model must be able to predict an operating point configuration for the laser with a requirement that this prediction is not computationally costly. This is a significant challenge for a laser as the physics is complex and nonlinear. Even a simple model of the laser leads to coupled Lorenz-like equations for the carrier and photon densities that can exhibit deterministic chaotic behaviour [18]. As such, laser operation at an operating point should not simply be understood as a stationary state where all variables are invariant in time. A laser is an oscillator, potentially resonating at a number of frequencies (modes), and although the average output power may be stable in time the internal electromagnetic fields will not be. In addition there is a random component to the fields due to spontaneous emission (SPE).
A characteristic specific to a compact model representation of an optical device in a circuit is that the output optical field will generally be referenced to a carrier frequency (ω c ). This greatly reduces the computational cost of the model. The actual optical output of the laser at a particular operating point or current bias (I d ) will be offset (detuned) from this frequency causing a linear variation in the phase of the optical output (chirp) and therefore will also not be stationary. The goal of a laser compact model must be to provide a configuration of the laser that approximates its behavior at an operating point and produces a minimal disturbance for the initial part of a transient simulation. Figure 1 illustrates a long (hundreds of ns) transient response for a laser simulation from an off state to an operating point above threshold. It is important to note that the modulation of the signal due to the optical carrier has been removed and the response is that of a complex envelope. Although at the end of the transient the laser state has evolved to a stable configuration, a significant amount of harmonic content is present due to detuning, the presence of multiple longitudinal modes and SPE. The large sinusoidal modulation of the field results from a chirp of the nominal carrier frequency due to either model detuning or a carrier dependent optical index (Henry effect). The two insets show the other sources of harmonic content-higher order longitudinal modes and SPE. Due to the presence of these other sources of harmonic content and the constant exchange of energy between the modes, the final operating point is not truly harmonic and so we refer to it as a quasi-harmonic steady state (QHSS). It should be noted that the length of the transient is considerable and precludes simply simulating mixed electrical-optical circuits by initializing them in an off state. There is a small amount of detuning to illustrate chirp. The two insets present the real part of the output at various time scales.
In this paper we will present a methodology to determine an operating point configuration of a traveling wave laser compact model [17] suitable to be incorporated into a time-domain optoelectronic circuit/system simulator such as OptiSPICE [9], [19]. The three-stage methodology will be shown in detail and then evaluated by performing transient simulations from an operating point and showing that the initial transient disturbance is not significant and that the simulation matches a simulation transitioned to steady state from an off state-which we will refer to as from-off. The ability of the method to predict the natural stochastic variation in the operating point will also be addressed. In the last section the implications of the three-stage method will be evaluated for incorporation into a circuit-level simulator.

II. TRAVELING WAVE MODEL
Traveling wave models (TWMs) of lasers are physically-based and provide a significant improvement in sophistication on REMs [12]- [14]. The model is derived directly from Maxwell's equations and produces a set of 1D first-order wave equations for the forward and backward propagating waves (see Fig. 2).

A. FIELD PROPAGATION
The derivation of the field equations results from two basic assumptions: 1) propagation of the electromagnetic field in a homogeneous waveguide which is described by a set of transverse modes each with a specific group velocity v g and effective index n g and 2) identification of a reference (or carrier) wavelength associated with the field. Such a model is a good candidate for a compact model as it models the evolution of the envelope of the optical signal and not the fundamental frequency.
If we assume for simplicity that only one transverse mode is excited then the field in the waveguide can be given bŷ whereÊ f andÊ r are complex envelopes for the forward and reverse propagating fields. These capture both the local magnitude and phase of the signals. (Note a^over a symbol indicates a complex value). The carrier frequency and wavenumber for the signal are given by ω c and k c . The equation for the propagation of the envelopes of the two waves (assuming they are slowly varying) can be shown to be [12]: whereβ is a propagation constant that captures effects such as a local perturbation in the effective index or gain/loss of the field.

B. LASER MODEL
One of the advantages of the TWM is its physical nature and the ease with which phenomena can be added to the model. To complete the model of a semiconductor laser a number of effects need to be added. First, the addition of a simple distributed rate equation of the excess carrier density (N (z)) is needed. Using N (z) the propagation constantβ can be used to model the local material gain and index variation due to the applied current. In addition, if desired, coupling between the forward and reverse traveling waves can be incorporated for distributed feedback based lasers. Finally, spontaneous emission can be modeled by adding stochastic terms to the equations (F f andF r ). Given these additions the two first-order wave equations describing the evolution of these two signals are [12]: whereκ f andκ r describe coupling of the two waves. We can define a photon density S(z) = |Ê f (z)| 2 + |Ê r (z)| 2 where it is assumed that the two signals are normalized correctly.
The equation describing the evolution of the carrier density is given by a distributed first-order rate equation: where we have a differential gain G 0 , η is the quantum efficiency, V l the volume of the active region, I d is the laser current and τ n the spontaneous emission coefficient. If desired, other effects such as electron diffusion or a non-uniform current distribution can be added.
To complete a laser model it should be noted that the presence of mirrors places a boundary condition (BC) on the two fields such that: whereR is the reflectivity of the two surfaces at the end of the laser. The output of the laser will then be given by: The three first-order rate equations (1, 2) and the BC relationship (3) comprise a basic distributed model of the laser, describing both the temporal and spatial evolution of the laser operation. A key distributed phenomenon captured by the model is the presence of longitudinal modes and gain dispersion.

C. LONGITUDINAL MODES AND GAIN DISPERSION
As is well known a laser will support longitudinal modes at frequencies given by: where m is the mode number and ω = πv g /L. For the TWM it is convenient to reference frequencies with respect to the carrier frequency ω c . We are then free 1 to choose ω c to be the dominant laser mode m and we have: where l = m − m . The TWM laser model will naturally produce signals (Ê f andÊ r ) dominated by these modal frequencies. The field component associated with l = 0 is essentially a stationary field with respect to the complex envelopes. Components associated with negative and positive l values represent modes below and above ω c respectively. There are also an associated set of spatial modes for which k 0 = 0 and k m = mk where k = ω /V g . Due to the frequency dependence of the gain in the material this harmonic content will have a significant effect on the behavior of the laser. The gain of the material is captured in the TWM by the propagation constantβ. In particular the imaginary part ofβ (which represents the material gain) will be a strong function of the frequency.
A commonly used expression to capture the frequency dependence of the gain is a Lorentzian. Using a shifted frequency ω = ω − ω c the gain due to the forward traveling wave in the material can be described by: where ω L = ω L − ω c and ω L is the frequency of maximum gain. In the time domain this relationship is represented as: The physical interpretation of this relationship is that of a band-pass filter relatingP f toÊ f and that at SS (dP f /dt = 0) P f will be linearly related toÊ f . As ω diverges from ω L the band-pass nature of the Lorentzian will causeP f to be small relative toÊ f . To incorporate the effect of this polarization on the field an additional term of the form of −k p (Ê −P) is added to the propagation equations [20]: where k p represents the strength of the coupling between the polarization and the field. At frequencies near ω L the polarization term will be close to zero and the fields will be unaffected. For frequency components away from the center of the Lorentzian, however,P f will be small making the term large and negative producing a suppression of the field and capturing gain dispersion in the laser. The frequency dependence of the gain is the primary determinant of the distribution of energy between the allowed longitudinal modes for the laser. The position and width of the Lorentzian in the frequency domain will determine the dominant mode and the partition of power between modes present during transient operation and at a harmonic steady state.

D. DETUNING, HENRY FACTOR AND GAIN COMPRESSION
A second source of harmonic content in the fields is the presence of a real part ofβ. This real part has two sources; 1) a constant part δ which typically represents a static detuning and 2) a Henry factor (α H ) capturing the carrier dependence of the real part the effective optical index. The static detuning is used to align the natural (modal) frequencies of the laser with the desired carrier frequency ω c [17]. For example if δ = −1/τ p where τ p is the average lifetime of a photon in the laser then ω c is defined to be the frequency of operation when the laser is at threshold. The Henry factor is commonly modeled using a linear relationship between the effective index and the carrier concentration and we can write for the complex propagation constant, where we have defined a differential gain G 0 , a gain compression factor g f and α l specifying losses due to scattering and absorption. The real and imaginary parts of the complex propagation constant are denoted by β r and β i . The physical effect of the β r term is to induce a ''rotation'' of the field as it propagates through the laser producing a phase shift on the output and chirp in the modal frequencies.
The gain compression factor is g f defined by [10], where is the gain compression parameter. This parameter is used to represent a number of strong nonlinear effects in lasers such as spatial and spectral hole burning.

E. DISCRETIZATION OF THE TRANSIENT EQUATIONS
A numerical model of these equations is simply formed by discretizing both space and time with z → z i and t → t j . The only care that needs to be taken is to use an up-wind formulation for the spatial derivatives. This leads to the following equations [17], where the equation forP f uses a second-order (trapezoidal) time derivative to minimize integration error. These equations plus corresponding ones forÊ r andP r and the BCs represent a complete discrete transient model of the laser.

F. HARMONIC STEADY STATE EQUATIONS
As a naive approach to finding a steady state operating point of the laser one can find a set of time-independent equations by setting all time derivatives to zero and dropping the stochastic SPE terms. If, for simplicity, we assume no detuning, a Henry factor of 0, no wave coupling and ω L = 0, with no SPE, the time-independent propagation equations become decoupled and straightforward to solve. For example the propagation equation for the forward traveling wave becomes, which can be directly solved to give: This trivial solution, however, only provides information on the primary longitudinal mode (l = 0 see (4)) and does not yield any information regarding the amplitudes of the other modes that we know will be present. In order to obtain such information we need to have the frequency information implicit in the full field equations.
The obvious way to proceed is to assume we have a harmonic steady state solution with all modes present. To find the harmonic steady state equations we begin with the time domain equations and move to the frequency domain by expanding the fields and polarization response using a Fourier series of the complex envelope at a particular position z with the laser modes frequencies as a basis. Ignoring modes beyond M the fields have the form of For a particular z of the forward going field (remembering that theê coefficients are functions of z) we have from (1a) (dropping the SPE term) by equating the i th coefficients, (6). Allowing us to write, Using a spatial discretization z = j z for j = 0 . . . m witĥ e fij denoting the i th mode at position j we then assume that K i is constant over the spatial step and we have: e fi(j+1) =ê fij eK ij z andê ri(j−1) =ê rij e −K ij z (11) The carrier density N (z, t) can be assumed to be static with respect to the optical frequencies and we have N (z, t) = N (z) and discretizing gives, To obtain a harmonic steady state solution the above three sets of equations (Eq.s (11) and (12) for 2 × n × m + m unknowns) must be solved subject to the boundary conditions present at the mirrors, e f ,n,1 =ê r,n,1 eK n,1 z andê r,n,K =ê f ,n,K eK n,K z (13) The equations (11)(12)(13) using the definition of S above can be written as a set of coupled nonlinear equations where the coefficientsê fi ,ê ri and N j are the unknowns. As the equations are nonlinear they need to be solved by an iterative technique such as Newton-Raphson. Unfortunately, these equations are not self-consistent. This is because of the requirement that each mode creates a magnitude profile from an appropriate gain that matches the reflectivity at the laser ends and forms a consistent power flow over one pass through the laser (once forward and once backward). For a single mode this is possible as an appropriate N (z) can be found that creates the needed gain. However, assuming the reflectivities are constant for all modes but the total gain is different due to the Lorentzian gain dispersion this requirement is impossible to meet.
In Fig. 3 this is illustrated for the simple case of two modes and a constant but differing gain for each mode. The gain of the laser is effectively determined by the dominant high gain mode and the fields will be self-consistent over a single pass. A higher order (lower gain) mode will, after a single pass, not return to the initial field strength resulting in a FIGURE 3. Illustration of the power distribution in two modes within the laser assuming two constant but different gains and a single mirror reflectivity. Note that the ends of the fundamental (high gain and larger magnitude) meet exactly while the ends of the higher order mode (smaller amplitude and gain) do not. VOLUME 9, 2021 constantly changing distribution for the mode. This can be confirmed by running a transient to steady state and extracting the N (z) distribution. It is found that the gain distribution in the laser does indeed produce a self-consistent solution for the dominant mode of the laser. However, the other modes are not self-consistent resulting in time varying harmonic content. This consequence is a reflection of the fact that the laser does not reach a true harmonic steady state, but rather a QHSS which is tightly bound to an attractor. Due to this situation the system of equations does not converge during iteration and this approach is not fruitful.

III. OBTAINING AN OPERATING POINT
As the previous section has shown that obtaining a true harmonic steady state solution to form an operating point configuration is unworkable we proceed by attempting to find an approximate solution. We use the REM to demonstrate that there will always be a dominant mode with the other modes being almost negligible. A dominant mode solution will then be found providing an excellent approximate solution-although with no other mode harmonic content. This motivates the overall approach; to first obtain a dominant mode for steady state operation and then to fine-tune it to reduce the startup transient. Harmonic content in the form of longitudinal modes is then introduced into this solution to create an appropriate operating point configuration. This three-stage method will be well suited to the task of providing an operating point for a circuit-level simulator.

A. MULTI-MODE POWER DISTRIBUTION
The multi-mode REM [10] is used to estimate the mode amplitudes in the laser: whereN is average carrier density,S i the average photon density of the i'th mode,S T = k i=−kS i is the total photon density and we are solving for 2k +1 modes in total. The gain associated with each mode is again given by a Lorentzian and we have: The rate equations forS i (14b) determine the threshold value (N th ) at which the photon density exceeds losses and lasing occurs. This equation shows that each mode will have its own threshold as determined by that mode's gain (15). The mode with the largest gain will always have the lowest threshold. Initially, when the laser is turned on there is a rapid increase in N and it moves above the threshold of many modes and all these modes are strongly excited as shown in Fig. 4. As the transient evolves the mode closest to the peak gain always becomes dominant and N decreases until it is slightly above the threshold value for this mode. The other modes will be suppressed because the carriers are consumed by the mode that is above their threshold. This is shown in Fig 4b where we see all modes initially have similar photon densities until N stabilizes with only the dominant mode above threshold containing the vast majority of the photons. Note that the harmonic content in the higher order modes can take tens of nanoseconds to come to equilibrium.

B. SIMPLIFIED STEADY STATE TRAVELING WAVE EQUATIONS
The previous section shows that it is expected that the configuration of the operating point of the laser will be dominated by the single mode nearest the peak of the gain distribution. Assuming single-mode operation we can proceed to find a solution quite simply as we avoid the singular matrix created by the reflection of multiple modes discussed in section II-F. This dominant mode could be chirped (shifted from ω c ) producing a modulation of the laser output. The other modes and the SPE will contribute to the configuration by providing a harmonic perturbation to this mode.
As we are only dealing with a single mode we can use the trivial solution from (8) and can write for the SS state equations: These equations have been written in terms of the magnitude of the envelopes and the phase has been removed. This can be done as S is only dependent on the field magnitudes and allows a solution to be found regardless of mode chosen to be dominant. Equations (16a) and (16b) describe the growth of the two fields as they pass through each section of length z with a gain determined by N k . The total gain along the laser must exactly equal the loss at the two laser ends for a self-consistent solution. Under these conditions we obtain a set of coupled nonlinear algebraic equations for the single-mode solution subject to the boundary conditions prescribed by (13). Again, it is only the imaginary part of β k that is used for (16a)-(16b).

C. FINDING A SINGLE-MODE STEADY STATE SOLUTION
Equations (16a)-(16d) are nonlinear but with an appropriate initial guess a solution can efficiently be found by iteration. Initial values forN andS are obtained by using the REM (14) with the assumption of the existence of a single mode.
To solve (16a)-(16d) we initially set N k =N for all k and the fields to zero. Using (16a) and (16b) in combination with the mirror's reflectivities, we can update the forward and reverse traveling E fields. Repeated iterations (alternately updating the distributed values ofÊ f ,Ê r and N ) continue until the values converge. The number of iterations is small and typically of the order of 10. We will refer to this approximate operating point as the Fast Steady State (FSS) solution. To illustrate this process Fig. 5 presents the initial REM value forN and the final FSS distributions for an example of a ridge laser with a length of 300 µm operating at 1550 nm.
The fields calculated from (16) and presented in Fig. 5b are the magnitudes of the dominant mode. If the dominant mode is the fundamental then these values can be used as is. However, for a mode that is not the fundamental (and therefore not stationary) we simply add the real part of β k for whichever mode is desired. This provides the number of 2π rotations of the phase appropriate to the mode.
Even though the FSS is a good approximation to an operating point solution, when used to start the TWM it will produce a transient. This is because the TWM is sensitive to matching N with the total photon density determined byÊ f andÊ r and even a close approximation may have a significant transient. An example is shown in Fig. 6. The two figures present the average carrier density and the laser output for a full transient from-off. Also shown in these figures are laser simulations started from the FSS operating point. For the vertical scales used in these plots the FSS appears to be a very good approximate starting point. However, the inset figures in Fig. 6 present the simulation from-off after it has run to QHSS and the initial 2 ns of simulation from the FSS. As can be seen in these figures there is a small but noticeable transient in both the average carrier density and the output field. It should be noted that the transient arises not from an inappropriate value forN or E (which lie within the natural variation of the OP) but because each implementation will be subtly different. To match exactly requires an exact match of the algorithms.

D. FINE-TUNING FAST STEADY STATE SOLUTION
The FSS method given above is very quick and efficient but produces a small but significant transient. To eliminate this transient we run the exact TWM algorithm on the FSS without the stochastic forcing of the SPE-producing a deterministic evolution that can be run to a convergence. With no SPE present the solution will stay in the dominant mode as there is no physical mechanism to create power in the other modes. This procedure is performed with all parameters (Henry factor, detuning, etc.) at their correct value and the resultant operating point will produce very little transient. Although this step (which we will refer to as Fine Tuned Steady State (FTSS)) is significantly more computationally intensive than the FSS (typically in the order of tens of thousands of iterations) it obtains a very clean single-mode operating point solution with no significant transient as can also be seen in Fig. 6 and its insets. Note that when started from the FSS OP it is still much faster than running a simulation from-off.

E. ADDING HARMONIC CONTENT
The previous two sections determine a distribution for N and the field magnitudes of an FTSS solution. To form a complete operating point solution, however, energy must be added into the other harmonic modes of this stationary laser solution. Let us quickly review all sources of harmonic content in the laser simulation. In total there are three sources: 1) The dominant mode not being the fundamental mode.
2) Field rotation (chirp) due to a non-zero detuning (δ) or Henry factor (α H ). 3) Laser modes other than dominant mode. The first two types of harmonic content can be seen in Fig 7  which presents the real part of the laser output for three simulations. In the first two the dominant mode is the fundamental with and without detuning and the last is the case where the dominant mode is the 1 st harmonic. If the laser is chirped (either statically or by a non-zero Henry factor) the laser output is sinusoidally modulated at a frequency determined by the difference between the dominant mode and ω c . If the dominant mode is not the stationary fundamental then the output is modulated by ω l where l is the order of the dominant mode. The chirp of the laser fields by detuning or a non-zero α H causes a ''rotation'' of the internal fields of the laser and complicates the determination of which mode is dominant. This is further discussed in appendix.
The formation of the FTSS has taken care of the first two types of harmonic content. The dominant mode was selected in the approximate FSS solution and the fine-tuning has introduced the effect of the Henry factor and static detuning (chirping the laser) and found a solution completely consistent with the TWM equations (without SPE present).
The final source of harmonic content is the non-dominant harmonic modes. Fig. 8 compares the last 1 ns of the simulation from-off with a fundamental dominant mode and no detuning from Fig. 7 with a 1 ns simulation starting at an FTSS OP. Thus we compare the FTSS with the expected harmonic content of the fields. In Fig. 8a the red trace is the stationary FTSS and the blue is the full transient. Note the stochastic nature of the operating point obtained from a full transient, with harmonics present, and the constant amplitude for the single-mode stationary state. The plot of the magnitude of the E-field over time shows both a high frequency modulation and random variation present in the full transient.
The presence of this harmonic content can be seen in Fig. 8b where the Fourier transform of the laser output is presented. It can be seen that even when the dominant mode is the stationary fundamental harmonic content is present in the higher order modes. The noise floor created by the SPE is present as are the distinctive peaks of the higher order modes. The diminishing of the modes as we move away from the dominant mode is expected due to the Lorentzian.
The random variation present in Fig. 8a arises from the SPE. In the TWM the SPE is implemented by adding small random components (F f andF r ) to the fields. See (1) [15]. The broad spectrum SPE injects energy into all laser modes which is then amplified by stimulated emission and sustained only for the cavity modes. The SPE thus seeds the higher order modes and provides an intrinsic background random field.
To add the appropriate energy to the other modes of the operating point we use a short pre-simulation (PS) using the full TWM transient model including SPE to naturally add harmonic content from other modes. We will refer to such an operating point as a PS configuration with the understanding that it is normally run after an FSS followed by an FTSS. The length of the PS is, of course, of concern as it will, in part, determine the computational cost of obtaining the operating point of the laser and circuit in which it is placed. However, as we shall find, it will be quite short (approximately 1 ns) as we are simply perturbing the FTSS solution.
To illustrate this procedure Fig. 9a presents a 10 ns operating point simulation from FTSS. For this run the Henry factor was set to 2.7 which detunes the laser such that the dominant mode is the 2 nd positive mode.
To analyse the harmonic content the field was captured in a sliding time window of 1 ns length. As this window was moved through the 10 ns simulation every 20 ps an FFT was performed to decompose the signal into modal powers. The one captured at 6 ns is shown in Fig. 9b. From the FFT the power present in each mode can be extracted and plotted as function of time. This is shown in Fig. 9c for the fundamental (mode 0) and the 4 lowest positive modes. The dominant mode (second harmonic) is well above the others and remains constant while the others start at what is essentially the noise floor for our spectrum and rise quickly to their steady state levels.
Inspection of these plots shows that the length of the PS sufficient to create the higher order mode excitation, introduce SPE, and establish an appropriate operating point configuration, is roughly 0.5-2 ns. After the creation of the dominant mode using the FSS and FTSS the introduction of other harmonic content is a fairly small perturbation to create the final operating point.

F. SUMMARY
To summarize the steps in obtaining an operating point we have: 1) Obtain an FSS solution with the harmonic parameters set to zero. Any mode may be selected to be dominant. 2) Fine-tune this operating point with the harmonic parameters and iterate using the TWM to create the FTSS. 3) Run a short PS to add in the other modes and produce the PS operating point.
The first step produces an approximate single-mode solution to the operating point, the second an exact single-mode operating point and the final stage creates a multi-mode operating point. It should be emphasized that all three operating points have utility within the context of either initiating a simulation from an operating point or within the context of a larger system/circuit-level simulation.
If one wishes to simply start from an approximate operating point and can tolerate a small transient during which the laser comes to harmonic equilibrium, the FSS provides a very quick solution. The use of the FTSS eliminates the small transient but the solution will still be single-mode. Finally, the use of the PS will recreate an exact multi-mode model of the operating point of the laser from-off within the limitations of the stochastic nature of the laser. When embedded as a compact model in a circuit-level simulator all of these configurations will have a use. This will be discussed more fully in Sec. V.
All of the configurations above assume a known dominant mode configuration for the operating point to be obtained. As has been discussed the dominant mode of the laser is typically the mode nearest to the peak of the gain curve. However, as will be seen in the next section complications arise due to the precise alignment of this curve with the mode frequencies and the stochastic behavior of the laser. Details regarding the impact of the dominant mode choice are presented in Appendix.

IV. RESULTS
This section will investigate various aspects of the work described above. It will first establish the need for an operating point solution, describe the wide variety of operating point configurations that can arise and establish the ability of the methodology of Sec: III to accurately capture these configurations. Finally, it will compare the statistical variation across harmonic steady state configurations of the method with respect to full transient solutions.   To characterize the effectiveness and reliability of the method given above a variety of simulations were performed. The basic approach is to run a given test configuration from-off in the TWM, and run it sufficiently long enough so as to reach steady state-this is the desired QHSS operating point. The calculated operating point is then used to start another TWM simulation which is compared to the desired one. The response of the output will, of course, be unsteady due to SPE and harmonic content but the variation should be within the same limits for both results. Each should be stable, agreeing on field levels and frequency content and as well, the operating point simulation should not have an initial transient in the output.
The results in this section pertain to an example III-V based laser [10] with parameters summarized in Table 1.

A. SIMULATIONS FROM-OFF FOR A VARIETY OF CONFIGURATIONS
In this section we will first present some full laser transients from-off to illustrate the complexity of dynamics in forming a QHSS. Figure 10 presents a TWM from-off transient for a simple default configuration with I d = 68 mA and the detuning, Henry factor and gain compression all set to zero. The Lorentzian gain curve is centered on the carrier frequency (ω L = 0). The first figure presents the response of the average laser carrier density (N ), the second the output electric field magnitude (|E|), and the final figure the mode intensity for the fundamental mode as well as the first four positive modes. For this case it can be seen that starting the laser from-off produces a large initial transient which manifests in the output field for over 10 ns.
The mode intensity plot (Fig. 10c) shows the evolution of the power distribution in the modes during the transient. Initially, there is quick growth in the fundamental (0 th ) mode as it is positioned at the peak of the gain curve. It quickly becomes dominant and as the transient progresses, power in the non-dominant modes drops and the power in the 0 th mode increases with the resulting power distribution determined by the relative gain of each mode. Although 15 ns is not enough time to transition to the final QHSS (the power in the 1 st mode is still dropping) the laser has come to a stable state that will continue to evolve and is very unlikely to change dominant mode. This transient is typical for a simple evolution to a QHSS with a clear dominant mode.
In Fig. 11 the output field is presented for a number of different bias currents with the same conditions for the other laser parameters. Although similar behaviour is noted for all bias conditions it can be seen that as the bias current approaches the threshold current (I th ≈ 33 mA) the transient response takes longer to settle to a QHSS. The transients presented in Fig. 10 and 11 are all roughly less than 25 ns. 2 However, much longer transients manifest themselves-either simply for stochastic reasons or due to the configuration that the laser is placed in. Some examples are presented in Fig. 12.
The first example (Fig. 12a) presents a laser with the gain curve centered on the fundamental mode. However, through sheer happenstance the mode that has become dominant after start-up is the 1 st negative mode (despite having a lower gain than the 0 th mode). This results in a long transient as this mode acquires power at the expense of the other modes. The QHSS that is formed using a dominant 1 st negative mode is potentially unstable and a longer simulation may show a transition to a QHSS formed using the fundamental mode.
In the second example (Fig. 12b) a laser biased at 38 mA with the gain curve placed at ω L = −2.5ω is shown. The placement of the gain curve midway between the 2 nd and 2 rd negative mode results in the gain of these two modes being essentially identical and both are likely to become the dominant mode of the resultant QHSS. Typically this results in a somewhat longer transient with one mode dominating 2 It is difficult to define a precise length of the transient due to the stochastic nature of the QHSS, but to quantify the time a settling time is defined by a peak to peak |E| < 10 5 being maintained for more than 10 ns. from start-up. In some cases, however, a mode will become dominant and then after a long time the other possible dominant mode will grow in intensity. The laser will then undergo a chaotic transition as the QHSS switches dominant modes. This can be seen to be the case in the mode intensity plot of this simulation. Initially, the laser evolves to a QHSS with the 2 rd mode being dominant. Although the second mode initially decreases in intensity after 400 ns it then increases in power until a chaotic transition occurs between 500 and 700 ns resulting in a QHSS with the 2 nd mode dominant. Detuning the laser or introducing a non-zero Henry factor can have similar results due to shifting the mode spectrum with respect to the gain curve.
The final case shown in Fig. 12c introduces the effect of gain compression on the laser transient in the presence of both static detuning and Henry factor. The combination of additional nonlinearity and the gain curve being misaligned with the mode spectrum produces a long region of initial instability. In this region no mode is able to become dominant. Finally after 150 ns the highest gain mode (2 rd ) becomes dominant and a stable QHSS is formed.
The simulations discussed above show a wide range of settling times from tens to hundreds of nanoseconds. It should be noted that the laser continues to settle for times longer than presented in these plots (this is not visible on the plot scale necessary to include the large transient). It can also be seen that the stochastic nature of the laser creates a distribution of settling times for a specific set of parameters. To investigate this sets of 100 runs were performed for a few laser configurations and then histograms of the settling times were compiled.
Three sets of these distributions are presented in Fig. 13 which expand on the simulation configurations used in Fig. 12. Figure 13a shows data for the default laser configuration with three differing bias currents of 38, 58 and 78 mA. As can be seen the bulk of the settling times are distributed with a peak at approximately 20 ns for all three currents and with some broadening as the current is increased. The higher current biases also exhibit a longer tail with long runs of the order of hundreds of ns being generated (such as shown in Fig. 12a for the case of I d = 78 mA which has a settling time of 110 ns.).
The second plot presents data for a laser biased at 38 mA and shifted gain curves. As we saw in Fig. 12b long settling times can occur if the gain curve is placed between two laser modes. This is confirmed in this plot with tight short distributions for ω L = −2ω and ω L = −3ω . However, for ω L = −2.5ω we have a broad distribution with a long tail. This is due to transients such as shown in Fig. 12a where neither the 2 nd or the 2 rd mode achieves dominance quickly.
The last plot (Fig. 13c) shows the effect of non-zero gain compression ( ) with increasing producing broader distributions and in particular creating a long tail even for a gain curve coincident with the fundamental laser mode.
The results presented in this section show the wide diversity of from-off transients resulting from variation in the laser configuration. They also demonstrate the chaotic nature of the VOLUME 9, 2021  Even if only a few long settling times occur, when they do it would dramatically impact the total circuit simulation time. So this presence of long transients to form the QHSS provides a strong motivation to develop the method of Sec. III for forming a QHSS operating point from which the laser can then be simulated.

B. OPERATING POINT DETERMINATION AND PRE-SIMULATION
The simulations in the previous section demonstrate that a wide array of physical effects determine the final QHSS present after a transient initiated from-off. In addition to the full transients from-off the plots in Fig. 10 and 11 present simulations started from an operating point constructed using this paper's methodology. As can be seen, at the scale used in these plots, the operating point appears to be well constructed and reproduces the QHSS of the transitions from-off. In this section a more detailed analysis and comparison of the operating point configuration is given with respect to the average values of the optical output and the harmonic content present at the operating point.
As discussed in Sec. III the basic method of constructing the operating point configuration is to first create an approximate dominant mode configuration (FSS), fine-tune it (FTSS), and then use a short PS to introduce harmonic content to the configuration. The result of this procedure for a laser biased at 68 mA with a chirp from a Henry factor of α H = 2.7 is shown in Fig. 14. Note that these plots are different from the full from-off plots in Fig. 10 in that the TWM simulations are not shown from-off or t = 0. These simulations plot the tail end of the TWM from-off simulation so they can be used as a reference for comparison with the calculated operating point configuration simulation.
The figure shows four plots for each of three simulations: carrier density, E-field magnitude, modes intensity over time and temporal spectrum. Because the FSS is a quick approximation we see that the transient in both the carriers and E-field is significant albeit short. The spectrum compares well, although the mode plot shows that there was a transient in the spectrum, including the dominant mode.
The FTSS compares well with the TWM. There is now only a transient in the mode plot for the non-dominant modes. The carriers and E-field vary in the same manner as the actual TWM output. The transient in the spectrum is only in the non-dominant modes because we have fine-tuned this configuration to the dominant mode eliminating the other modes. In these examples the initial mode power distribution VOLUME 9, 2021 of the FSS appears better than the FTSS. This is due to the approximate FSS containing power in the non-dominant modes. However, the initial power distribution in the FTSS is not determined by the model equations and is not strictly correct.
The PS simulation is an excellent approximation to the TWM steady state. There are no transients on any outputs and the power in the non-dominant modes matches that present in the simulation from-off.
These results show that for the case of a chirped fundamental as the dominant mode all three methods of creating an OP result in a viable configuration-with the understanding that the FSS will produce a small transient and the FTSS will produce a small spectrum transient if used alone. The methodology provides for a choice of dominant mode to form the QHSS and that will now be discussed.

C. EFFECT OF DOMINANT MODE
Choosing a valid dominant mode is critical for forming a stable operating point. As was seen in the simulations from-off given above, the mode nearest (or even at) the Lorentzian peak may not be the mode which forms the QHSS to which the laser evolves. For a situation where the peak is near a mode it will most likely be the dominant mode in the QHSS. However, the peak can be located anywhere in the modal spectrum and, in particular, the Henry factor will move the modes through the spectrum as a function of laser bias.
Simulations from-off show that a QHSS close to the peak of the Lorentzian can be quasi-stable even if not the nearest to the peak. For from-off simulations such quasi-stable operating points will produce long transients with chaotic transitions from one QHSS to another (see Fig. 12). If such a fundamentally non-dominant mode is chosen to form an operating point the TWM will at some point undergo a chaotic evolution to a higher gain mode-a process that can take longer than starting the laser from-off. Fig. 15 illustrates the situation where the choice of the initial mode used for the starting operating point is important. For this simulation the Henry factor was set to α H = 2 and the choice of parameters places the gain peak roughly half way between the 1 st and 2 nd mode. In Fig. 15a we see a simulation of 50 ns for an OP initiated in the 1 st mode-which is stable. The second set of figures (Fig. 15b) presents the same simulation but with the 2 nd mode used to form the OP and it is again stable. However, in the final set of figures the OP is formed from the fundamental (0 th ) mode and the simulation displays a chaotic transition to a QHSS with the 1 st mode dominant.
The transient simulation from the OP formed using the fundamental mode is quite stable lasting approximately 10 ns, but the laser state is noticeably different from a typical simulation from-off with this set of laser parameters. At the point of transition the TWM envelope changes the direction of rotation, settles on the first harmonic as the dominant mode and a stable laser state forms. The transition is very clearly seen in the modal evolution with the 1 st mode rising in magnitude until at 20 ns it becomes the dominant mode. This plot also makes obvious the chaotic period from 10 ns to 30 ns where there is a large amount of instability in the mode powers-also apparent in the plot of the field magnitude.
The methodology for the creation of the OP is agnostic to the choice of the dominant mode and can be used to determine the behaviour of the laser when placed in a particular configuration. However, care must be used in choosing the mode. If a mode is chosen that is not adjacent to the gain peak with respect to frequency, the OP may well be unstable and the simulation will reflect that instability.

D. VARIATION OF PARAMETERS
A requirement of the method is that it provides an effective operating point solution for a full range of parameter values. This should include bias current ranges from threshold to realistic maximum values, Henry factors typical for the laser in question, variation in the position and width of the Lorentzian gain profile, appropriate values for the gain compression constant and static detunings for accommodating system-level requirements on the specification of the carrier frequency. The basic requirements of a successful operating point construction are the absence of a significant initial transient, an appropriate set of modal powers and long-term stability.
To investigate the robustness of the method a large number of simulations were performed varying the parameters of interest. For all cases the method was found to work well-subject to the selection of an appropriate initial dominant mode. A few characteristic examples are presented in Fig. 16 where the E-field magnitude, mode evolution and spectrum are shown for six different parameter configurations. The bias current, chirp from detuning and chirp from Henry factor are varied. For the output E-field the final 5 ns of a TWM run from-off is compared with a run from the operating point solution. The spectra are calculated on the last 1 ns of E-field data. The spectrum is shown first as a plot of the mode strengths over time and then a normal spectrum frequency plot comparing the reference TWM with the operating point simulation.
For all the configurations the operating point simulations agree well with the TWM QHSS obtained from simulations from-off. As expected for the QHSS there is a continuous fluctuation in the output due to SPE and mode competition. In the time domain the variation due to the stochastic SPE means we cannot expect an identical output, but we do expect the fields to stay within the limits of field variation. The E-field plots show that the natural amplitude variation is in good agreement and that the operating point has no noticeable transient at the start.
The mode plots over time show that the modal energies have the same average behaviour and variability as the from-off TWM reference-subject to the continuous fluctuation due to SPE. These plots also show that the operating point creates no transients with respect to the power in the non-dominant modes. The spectrum of the final 1 ns shows excellent agreement. Although these simulations present a transient simulation of 5 ns in length, simulations of these configurations (and others) were run for hundreds of ns to confirm the long-term stability of the operating point created. We conclude that the OP provides a very accurate operating point solution for a wide range of the TWM parameters with negligible transient.

E. STATISTICAL VARIATION
The results presented above show that the method of Sec. III produces individual operating points that reproduce the QHSS created by from-off simulations with essentially no initial transient and that simulations run from these operating points have the appropriate modal energies and behaviour. However, due to the stochastic nature of the laser model, when run from-off the QHSS obtained from a set of runs will exhibit a natural variation in average output intensity and internal variables such as the average carrier density. We now wish to investigate if the method proposed will also produce this natural variability.
The first two stages of the method of Sec. III (the FSS and FTSS) are deterministic so the source of the variability in a QHSS created using this method is the third stage PS. It was shown above that a 1 ns PS is sufficient to produce the modal energy distribution appropriate for an operating point. In this section we show that it is also sufficient to create the expected variability in the operating point. This is done by running 100 simulations of each configuration and comparing the time averaged magnitude of the E-field and average carrier density for the final 1 ns of the TWM and its corresponding operating point. These values will be plotted in a histogram to show the statistical variation of the two simulations and how well they agree. Fig. 17 shows histograms for five different bias currents ranging from 38 mA (slightly above the threshold current I th = 33.2 mA) to 78mA. Fig. 17a shows the carrier density histograms and Fig. 17b the E-field histograms. The carrier density varies inversely with the E-field magnitude as expected. The spread of the E-field also increases as the bias current decreases and approaches threshold, at which point it will cease to lase. Note that the E-field magnitude scales are not all identical to facilitate comparisons with respect to the width of the distributions, but are shown approximately relative to each other. A comparison of the histograms illustrates that the operating point distributions match the TWM from-off distributions.    18 shows histograms for three different detunings: δ = 0, δ = 0.5ω , δ = ω . The introduction of a detuning creates a situation where the from-off simulations can create differing operating point configurations dependent on the dominant mode of the QHSS. We see that with nonzero detunings carrier density and E-field cluster around two values-a different value for each dominant mode. It should be noted that for the simulations run from an OP the dominant mode is a parameter and runs were initiated so that the number of runs in each configuration matched the from-off simulation.
A similar set of histograms is shown for the Henry Factor in Fig. 19. The Henry Factor has a similar effect to detuning in that it produces a chirp and so we expect similar results. The plot shows three values of Henry Factor: α = 0, α = 2.7, α = 5.4. And again we see the nonzero Henry Factor producing two high frequency values corresponding to the two steady state dominant modes. The distributions match well with the from-off simulations.
The last set of histograms shown in Fig. 20 show the statistical results while varying the gain compression parameter . The values are = 0, = 2, and = 5. A nonzero gain compression results in a distribution of dominant modes, each with its characteristic carrier values and E-field magnitudes. Note that the gain compression of = 5 results in three different dominant modes although the −1 mode occurs only once in 100 simulations. This set of histograms also shows the excellent agreement of the operating point distributions with the TWM from-off simulations.
From this set of results we can conclude that a pre-simulation of 1 ns is sufficient to create the expected variability in the operating points.

V. DISCUSSION AND COMPUTATIONAL CONSIDERATIONS
Although the ability to start a stand-alone laser simulation from an operating point is a useful capability, the VOLUME 9, 2021 primary motivation for the methodology described above is to facilitate the use of the TWM in a circuit-level simulator. As described in Smy and Rasmussen [17] the TWM is well suited for incorporation into a circuit/system-level simulator due to the common use of a complex envelope modulating an underlying carrier. It is also a very sophisticated physical description of the laser capturing an array of phenomena not modeled by simpler descriptions such as the REM. Generally, a circuit-level designer will desire to simulate a transient simulation of a circuit from an operating point and the initial task of the simulator will be to obtain an operating point for the circuit and therefore, of course, of all the elements in the circuit. Although the type of simulator being used will determine the methodology for obtaining this operating point, in essence a SS solution of a coupled set of nonlinear algebraic-differential equations will need to be found. This will typically be done by using a Newton-Raphson iteration procedure during which the element compact models will be called upon to determine an operating point configuration for a given set of inputs. Even in moderately sized circuits 3 this procedure can require many iterations and an efficient simulation engine will require quick compact model updates.
The iterative procedure used by a simulator requires the outputs of each compact model to converge to a specified criteria. For a compact model of a laser such as the TWM this presents an issue as the model is inherently stochastic due to SPE and the operating point itself is actually a quasi-harmonic steady-state (QHSS) due to chirping of the envelope and higher order modes. However, the three-stage methodology for finding this QHSS given above is well suited for use in obtaining an operating point. Reviewing the three stages:

A. FAST STEADY STATE (FSS)
This initial stage provides a very quick deterministic approximate solution to the QHSS. Because it is single-mode the internal quantities are not in precise agreement with the fundamental TWM equations, leading to a small transient when used as the starting OP. However, the outputs of the compact model 4 fall within the natural range of variation due to SPE and harmonic variation. The method itself requires typically only ten internal iterations to converge to a new OP and thus requires on the order of a thousand floating point calculations. This computational intensity is not larger than one would find using an array of complex transistors modeled with BSIM 4 The optical fields produced as the facets are typically accurate to within 1% of the value of the final QHSS mosfet compact models. Each update takes roughly 50 ms for a 51 element model in Matlab. 5

B. FINE TUNED STEADY STATE (FTSS)
The second stage in finding the QHSS involves fine tuning the FSS by iteration using the TWM equations with SPE set to zero. This provides a deterministic single-mode solution that is consistent with the TWM and produces no transient when used as an OP. However, it is computationally relatively expensive taking tens of thousands of iterations to converge-roughly equivalent to running the TWM equations for 1 ns. Each update takes roughly 10 s for a 51 element model in Matlab.

C. PRE-SIMULATION (PS)
The final stage of the creation of an QHSS OP is a short 1 ns pre-simulation using the full TWM equations with SPE. This adds energy to the higher order modes, variability to the OP and is obviously not a deterministic solution due the random nature of the SPE. A PS of 1 ns takes roughly 10 s for a 51 element model in Matlab.
The above model update times were for a Matlab implementation and comparisons to commercial simulators written in a compiled language such as 'C' can not be readily made. However, it is clear that when finding a circuit-level OP the use of the FTSS and PSS are precluded from the initial iteration due to the computational intensity of each stage. Although obviously much faster than simulating the device using the full TWM equations from-off 6 both the FTSS and PS require a large number of iterations (in the order of tens of thousands). The FSS, however, is ideally suited for this purpose as it is fast (200 times faster in the Matlab implementation) and deterministic.
Therefore a procedure can be proposed for the use of the three-stage methodology within a circuit simulator. During the initial stage of convergence the FSS is used. Once a circuit-level OP has been found that is consistent with the FSS the procedure can be continued for a small number of iterations using the FTSS. Finally one last pass is done to obtain the PS to provide a QHSS solution appropriate to start a transient solution.

VI. CONCLUSIONS
This paper presents a method of obtaining an operating point configuration for a laser model based on a traveling wave model. Intended for use in a circuit-level simulator this configuration can be used as an initial starting point in a circuit starting at a bias point. A three-stage method is presented that provides a fast approximate solution for convergence to an initial circuit-level operating point, a precisely tuned OP consistent with the traveling wave model equations and finally a harmonically complete operating point. A key aspect of the model is the addition of harmonic content to the stationary solution of the laser model, recreating its behavior at a QHSS solution.
The effectiveness of this solution is shown by initiating transient simulations from a wide variety of operating points and observing the size of the initial disturbance in the laser output and the evolution of the power in the laser modes. Included in the model formulation are gain compression and dispersion effects, laser chirp due to variation in the effective index of the laser mode, and spontaneous emission. The results presented show the effectiveness of the method in obtaining a sufficiently accurate operating point, as well as the need to include harmonic content to incorporate multi-mode behaviour, spontaneous emission, and chirp in the laser due to a Henry factor.
It must be stressed that the effectiveness of this solution is not only measured by the accuracy of the final operating point but by the practicality of the solution in a circuit simulator. Because the operating point will be calculated many times for each operating point required, the computational cost needs to be extremely small for the circuit simulator's initial stage of convergence. The three-stage methodology not only provides a reasonably fast final operating point calculation but also an extremely quick estimate for the circuit simulator's first stage.

APPENDIX DOMINANT MODE DETERMINATION
The dominant mode at which the laser operates is largely dependent on the position of the Lorentzian gain curve with respect to the mode structure. If the peak of the gain curve is close to a particular mode then this mode will generally be the dominant mode in any OP formed. As the dominant mode is often chosen to be ω c (in the absence of significant detuning) this will typically be the fundamental mode. However, if ω L lies midway between two modes the dominant mode can be either of these two modes. For example if ω L is halfway between the 2 nd and 2 rd harmonics (ω L = 2.5ω ) the dominant mode after a full transient from-off will be the either the 2 nd or 2 rd mode. The dominant mode frequency is also effected by the presence of a non-zero value of static detuning (δ) or Henry factor (α H ) which cause an internal rotation of the laser field shifting the mode frequencies with respect to the gain curve. This is illustrated in Fig. 21 which shows the forward and reverse traveling wavesÊ f andÊ r along the laser length at an operating point with an effective non-zero detuning of δ = 0.4ω for three cases. In Fig. 21a the fields are identical to the original stationary solution with the fundamental mode dominant but are rotating over time at rate of − ω. Fig. 21b shows a harmonic solution with the first harmonic mode dominant that is also rotating over time but at a rate of ω 1 − ω. The third figure shows a configuration for α H = 2 which produces an output approximately equivalent to the two other cases. At the output of the laser all three operating points will have the frequency shift (chirp) equal to the base frequency of the dominant mode plus the chirp of ω = −0.4ω .   Fig. 22a the first case shows the dominant temporal frequency of the fundamental ω 0 shifted by − ω due to a static detuning. The spatial frequency (in Fig. 22b) is not shifted and remains at k 0 . For the second case in which the dominant mode is the first harmonic the temporal frequency is seen to be ω 1 − ω while the spatial frequency is now k −1 . The gain curve of the Lorentzian is shown to illustrate that two configurations have similar gain. It should be noted that the these two configurations are at identical operating conditions and use the same laser parameters and yet the two steady state solutions are quite distinct operating points each determined solely by which mode is chosen to be dominant. The choice of which one to use is therefore arbitrary.
The overall affect of a non-zero Henry factor is quite similar to static detuning. From (17) we see that both contribute to the real part of β. Because α H is multiplied by the gain (which is N dependent) it will vary as N . Still, the effect of the Henry factor will be approximately the same as detuning with the average value of β due to the Henry factor. This is shown in Fig. 21c where the case of a dominant 1 st harmonic was chosen. For a Henry factor of α H = 2 the frequency shift is approximately −0.4ω with a β ave = 39. This agrees with the detuning constant β = δ/2 = 41.5 in Fig. 21b producing similar spectra to this case (see Fig. 22), with a temporal frequency shift from ω 1 of − ω and the dominant spatial mode being at k −1 .
If, as for the first two configurations presented in Fig. 21, the peak of the Lorentzian lies between two laser modes the random nature of the SPE implies that even with identical parameters and inputs the laser will not always end up at the same operating point after a full transient. Not only will each simulation be unique but the dominant mode will not always be the same for each. There will be a statistical distribution of mode configurations for each set of parameters.
To choose the dominant mode with which to form the operating point configuration one needs to determine the position of the laser modes in frequency with respect to the FIGURE 22. E field frequency spectra for an operating point with an effective detuning of ω ≈ −0.4ω for the three cases in Fig. 21. (a) Temporal frequency spectrum of the laser output. The gain curve of the Lorentzian is shown (not to scale) for illustrative purposes. (b) Spatial frequency spectrum which is always an integral number of k .
Lorentzian peak-including the effects of detuning and the Henry factor. Detuning simply shifts the mode structure and can easily be accommodated. To calculate the shift due to the Henry factor requires an appropriate carrier distribution N (z). As N (z) varies only slightly through the laser, the distribution appropriate with α H = 0 can be used to determine the effective detuning. Once the position of the modes are determined with respect to the Lorentzian the mode nearest to the peak can be chosen. As other nearby modes might be semi-stable (particularly if the Lorentzian peak lies close to the point midway between two modes) one can also investigate other modes nearby. TOM J. SMY received the B.Sc. and Ph.D. degrees in electrical engineering from the University of Alberta, Edmonton, AB, Canada, in 1986 and 1990, respectively. He is currently a Professor with the Department of Electronics, Carleton University, Ottawa, ON, Canada. He has published over 100 journal articles. He has coauthored the OptiSPICE, Atar, SIMBAD, and 3-D-films simulators. His current research interests include optical, thermal, and multi-physics simulation of electronic devices/packages and modules, the study and simulation of thin film growth and microstructures, and a variety of backend processing projects. VOLUME 9, 2021