Theory of frequency modulated combs in lasers with spatial hole burning, dispersion and Kerr

Frequency modulated (FM) frequency combs constitute an exciting alternative to generate equidistant spectra. The full set of Maxwell-Bloch equations is reduced to a single master equation for FM combs with fast dynamics to provide insight into the governing mechanisms behind phase-locking. It reveals that the recently observed linear frequency chirp is caused by the combined effects of spatial hole burning, group velocity dispersion and Kerr due to asymmetric gain. The comparison to observation in various semiconductor lasers suggests that the linear chirp is general to self-starting FM combs.

Optical frequency combs [1,2] are lasers whose spectrum consists of a set of evenly spaced modes that obey a defined phase relation. In the time domain, these lasers emit a signal, which, despite an eventual constant phase drift due to a non-zero carrier envelope offset frequency, is periodic. In literature, frequency combs are mostly linked to ultra-fast lasers that emit short pulses. However, the Fourier theorem states that a comb spectrum is generated by any periodic signal, regardless of its shape. A periodic frequency modulated (FM) signal is another example that fulfills this criterion. The first studies to generate such an FM laser output trace back to the 1960s, only a few years after the demonstration of the first laser [3]. An active intracavity phase modulator was used to generate FM oscillations in He-Ne [4] and later in ND:YAG [5] lasers. Just from the similarity of the optical spectra to the Bessel amplitudes, it was concluded that FM lasers obey a sinusoidal modulation of the output frequency [6,7] and this picture remained dominant for over 50 years.
Today, FM combs experience a renaissance. In 2012, it was shown that quantum cascade lasers (QCLs) can be used to generate combs, whose intensity remains approximately constant [8]. In contrast to the work from the 1960s, the generated FM comb in QCLs is self-starting. The possibility of generating self-starting combs using the nonlinearity provided by the gain medium is particularly appealing for fundamental laser science and the study of self-organization in complex nonlinear systems, but also of great interest for many applications. FM combs can be generated in fast gain media, e.g. QCLs that do not satisfy the conditions for passive modelocking [9]. They are self-starting, requiring no additional components e.g. saturable absorbers, which is interesting for comb generation in interband cascade lasers (ICLs) [10,11]. Both QCLs and ICLs emit in the midinfrared region that is particularly appealing for dualcomb spectroscopy [12,13].
In this work, we provide a rigorous theoretical and numerical study of FM combs that is driven by recent experimental results [14,15]. A highly optimized simu- * benedikt.schwarz@tuwien.ac.at lation tool was developed to reproduce the experimental results, to identify trends and to isolate the most relevant terms in the full set of nonlinear coupled differential equations. With this knowledge, we derive a simplified master equation for FM combs. It provides the eagerly awaited intuitive explanation of the phase-locking and answers the following questions: -What triggers self-organization of the phases in FM combs to overcome chaos?
-Why does the linear frequency chirp emerge from this competition, overcoming other solutions?
-Do FM combs require a fast gain medium?
The possibility to generate FM combs with QCLs is mostly explained through their fast gain dynamics and four-wave mixing via the occurrence of population oscillations that respond in anti-phase to oscillations of the light intensity [17,18]. A modulated intensity saturates the gain more than a constant intensity. Following the maximum emission principle [19], amplitude modulations will be suppressed to maximize the output. In fast gain media, this effect is particularly strong. This is also the reason why a slow gain medium is required for pulse generation. There, the suppression is compensated and reversed by fast saturable absorption.
The main issue with this concept is that any phase arrangement that minimizes amplitude modulations is equal in energy, which should result in a chaotic phase modulation [20]. Figure 1a shows the corresponding numerical simulation result, reproducing the expected pseudo-random behavior. Experiments, however, clearly show the formation of a distinct periodic pattern and the generation of a frequency comb [14,15]. Figure 1b shows the experimental results of a QCL frequency comb with the characteristic linear phase pattern that covers the range from −π to π. Note that we plot the intermodal phases, i.e. the phase difference between adjacent modes. This linear pattern corresponds to parabolic modal phases and a chirped instantaneous frequency.
The suppression of amplitude modulations can be interpreted as repulsive coupling of the intermodal phases.  The occurrence of self-organization in repulsively coupled systems can also be found in other fields of research, ranging from splay states in Josephson junctions [21] to cluster states in complex networks [22]. Such phenomena can be explained by additional contributions that induce an imbalance to favor one among many other solutions. Such an imbalance can for example be triggered by a finite GVD, as shown by our numerical results in figure 1c. This is particularly surprising, as the experimental observations of FM combs in QCLs were found in dispersion compensated cavities.
In the following, we explain why the GVD plays such a crucial role and which other effects are required to explain the experimental observations. The starting point is a system of eight coupled nonlinear differential equations that describes the manifold physics of a laser. The system is based on the spatio-temporally resolved Maxwell-Bloch equations in the slowly varying envelope approximation [23,24], which includes the effects of GVD and Kerr nonlinearity that have been mostly omitted previously. The details of the model and all derivations can be found in the supplementary material. While this full model is capable of a quantitative analysis, it cannot provide an intuitive understanding of the underlying physics.
Laser models with reduced complexity can yield an intuitive interpretation. An example is the Haus master equation for mode-locking with saturable absorbers [9]. Such models utilize the adiabatic approximation to eliminate variables, e.g. the induced polarization and carrier populations. However, with its application, the physics behind FM combs disappear. In the adiabatic elimination the response of these variables is assumed to be instantaneous, which is equivalent to approximating their transfer function by a constant, e.g. H(ω) = a/(1+iωT ) ≈ a. With this, the information on a remaining small phase delay that can accumulate over hundreds of round-trips is entirely lost. In order to recover sufficient information, we use a Taylor  we rewrite eq. 1 in terms of power and phase E ± = P ± exp(iφ ± ) and neglect several minor contributions, e.g. terms with ∂ t A are smaller than terms with ∂ t φ. The reduced master equation for FM combs reads: Now it will become clear why we introduced the Taylor expansion (highlighted terms) instead of the adiabatic elimination. As a frequency modulation is essentially a modulation of the phase, FM comb operation is mostly governed by eq. 3. Without the highlighted terms the phase dynamics would be lost and with that also the physics of FM combs. The first highlighted term in eq. 3 dampens phase oscillations, tending towards single mode operation. The second highlighted term in eq. 3 is due to spatial hole burning and facilitates multi-mode operation. The dispersion term proportional to k determines the evolution of the cumulative phase shape and favors a convex or concave parabola, depending on the sign of k . This results in a chirp in one direction or another. The Magnitude of k is directly related to the chirp of the intermodal phases. If it is large enough such that the intermodal phases cover the full range of 2π, the laser can produce a stable periodic output. A further increase makes the chirp larger which requires a narrowing of the spectrum (figure 2b). The residual amplitude modulation present in both the experiments and simulations (figure 1) are mostly due to the second highlighted term in eq. 2.
The situation becomes slightly more complex, when considering an additional Kerr nonlinearity β. In that case, eq. 2 and 3 are dynamically coupled. The second highlighted term in eq. eq. 2 represents a source term, forcing P ± to oscillate in a similar manner. Coupling P ± back in eq. 3 through the Kerr term, one sees a similar effect as the GVD, as approximately both are influencing the phase through a term proportional to (∂ t φ ± ) 2 . Figure 2a shows the numerical simulation of an FM comb generated by a Kerr nonlinearity with zero GVD. One can observe a slight bending of the intermodal phases, also present in the experimental data (figure 1b).
The Kerr term mostly originates from a change of the real part of the refractive index with the population distribution. This is connected to an asymmetric spectral gain profile and commonly expressed through a non-zero linewidth enhancement factor (LEF) [17]. A non-zero Kerr term or LEF strongly shift the range of GVD required for FM comb operation. This explains why QCL frequency combs have been found close to zero GVD. Figure 2c shows the range of GVD required to obtain FM comb operation for three different values of β that correspond to realistic values of LEF for QCLs at roomtemperature [25]. Altering the shape or width of the gain changes the required GVD and thus suggests con-sidering this effect in the design of broadband FM combs. It will be interesting to see if the observed behavior at very high values of GVD with the continuous narrowing of the spectrum can also be reproduced in the experiment, or if above a certain threshold the laser becomes unlocked again due to effects neglected in the reduced model. First attempts of solving the full model for high dispersion yielded again an unlocked state, but the results remain inconclusive due to numerical issues that are known from convection-diffusion problems. A detailed investigation will be part of future work.
We did not find a reason why FM comb operation should be strictly limited to fast gain media. In QCLs, with their fast dynamics, even a small gain asymmetry leads to a considerable Kerr nonlinearity. In interband semiconductor lasers, this contribution is attenuated due to the slower dynamics, but the asymmetry is much more pronounced. Moreover, dispersion driven FM comb operation appears to be independent of the gain dynamics. Hence, we believe that FM comb operation with a linear chirp is a general phenomena. As an example the presented theory also explains the experimental observations of self-starting FM oscillations in a InGaAsP laser diode [7]. The simulated intensity spectrum with the linear chirp fits the measured modal amplitudes much better than previously assumed Bessel amplitudes (fig-ure 2d). Our theory can further explain observations in numerical simulations of quantum dot and quantum well lasers [26,27] and recent results on the demonstration of FM comb operation in interband cascade lasers [10].
In conclusion, we provided detailed insights into the formation of frequency combs in single section lasers without saturable loss. Going beyond the adiabatic approximation, we derived a master equation to explain the physics and identified the most relevant contributions. In accordance to this, an FM comb requires: spatial hole burning to trigger multi-mode operation, gain saturation to suppress amplitude modulation and a minimum, but finite contribution from GVD or Kerr due to gain asymmetry that gives rise to a chirped output. A minimum amount is required such that the intermodal phases can cover a range of 2π over the spectral span to suppress amplitude modulations. Further increase will enforce a narrowing of the spectrum. The presented theory is capable of explaining experimental observations in various types of semiconductor lasers, indicating that the linear chirp is a general phenomenon behind the nature of FM combs.
This work was supported by the Austrian Science Fund (FWF) within the projects "NanoPlas" (P28914-N27) and "Building Solids for Function" (Project W1243). We begin with the density matrix formalism for a twolevel system [23]. The total Hamiltonian of the system can be represented asĤ =Ĥ 0 +Ĥ . Interaction of the system with the electric field is described witĥ H = −μE(t), whereμ is the dipole moment operator. We can write: where indices l and u denote lower and upper states and W l,u are the corresponding energies. The time evolution is governed with the von Neumann equation, with the time flow in direction ∼ e +iωt : bearing in mind that ρ ul = ρ * lu . Let us also denote the transition frequency ω 0 = ω u − ω l . Since the matrix elements ρ ll,uu represent occupation probabilities for the levels, by multiplying them with the total sheet density n tot , we obtain equations for the surface densities of the levels n l and n u . We also normalize in the same way n ul = n tot ρ ul . Furthermore, we add the dephasing processes and transitions, modeled with the polarization dephasing time T 2 and different transition lifetimes. Carrier diffusion is introduced via the coefficient D and a pumping current J to the upper level is also included. The pumping current is normalized to the elementary charge. The new set of equations are: where the index g stands for the ground level and T ij represents the transition lifetime between levels i and j.
Pumping current J is assumed to be constant over the entire laser cavity. In the next step, we find the macroscopic polarization P , where L is the thickness of the doped region: Employing the last relation, we can write the wave equation, where Γ is the confinement factor and n is the refractive index: We can now express the electric field E(z, t) as a sum of backward and forward propagating components, which allows us to do the same with n ul . Furthermore, all populations contain a grating component besides a spatial independent average component in order to account for the spatial hole burning (SHB). This ansatz is summed up in: Substitution of expressions (A6) in the density matrix equations (A3) and the wave equation (A5) yields the equations for the envelope functions, after applying the slowly varying envelope and rotating wave approximations and adding the intensity loss coefficient of the waveguide α w to the field envelope equations. Additionally, the Kerr effect is included and modeled as a phase variation through the Kerr coefficient β, so we can write: Equations (A7)-(A14) make a complete set of coupled spatio-temporal density matrix and Maxwell equations similar to the ones used in [24].

Group velocity dispersion
In reality, a laser cavity always possesses a non-zero dispersion which has a profound impact on the intermode dynamics of the laser. Hence, if one aims for modeling the delicate process of phase-locking of the laser, it is of interest to obtain a time-domain wave equation that could fully describe the evolution of the electromagnetic field inside the cavity. We start from a one-dimensional wave equation that considers the possibility of dispersion: where P stands for the nonlinear macroscopic polarization and the dielectric permittivity ε(t) is calculated as: The convolution integral in equation (A15) is not convenient for calculation, so we will try to obtain a more appropriate form. After applying the Fourier transform to equation (A15), knowing the identity F( ∂ n ∂t n x(t)) = (iω) n X(ω), we obtain the exact form of the equation in the frequency domain, where ω stands for the instantaneous frequency: We can now introduce the wave number k(ω) = ωn(ω)/c, where ε(ω) = n 2 (ω). Also, let ω 0 be the carrier frequency, which is introduced in the ansatz (A6) and k 0 = k(ω 0 ). From (A17) we then obtain: Knowing that the instantaneous frequency ω is in the vicinity of ω 0 , we can write the Taylor expansion of k(ω): where k (m) 0 = ∂ m k/∂ω m and the group velocity can be defined as v g = 1/(dk/dω)| ω0 ≈ c/n. Now it is a good idea to insert ansatzes (A6) in equation (A17) and after some calculation, equations for the envelope functions can be derived: which can also be written as: Keeping in mind that E, n u l ∼ e iωt , we can conclude that E ± , σ ± ∼ e i(ω−ω0)t . This means that (ω−ω 0 ){E ± , σ ± } = −i∂/∂t{E ± , σ ± }. After inserting this in equation (A19), we get: Combining expressions (A20) and (A21) one can obtain: After some derivation and disregarding all derivatives of the same or higher order than O ∂ 3 ∂t 3 in the above equation, one gets: Looking at the last expression, we can neglect the second order derivatives in terms of both time and space in the first bracket for two reasons. Namely, they come with opposite signs and practically cancel out in addition to being minor in value, ∂ 2 ∂z 2 << k 0 ∂ ∂z and ∂ 2 ∂t 2 << ω 0 ∂ ∂t . Furthermore, both the first and second order time derivative in the bracket on the right-hand side of the equation can be neglected since the induced polarization σ is a perturbation and is small compared to ε 0 E. Lastly, we can single out the term i k (2) 0 2 ∂ 2 E± ∂t 2 as the term which describes the group velocity dispersion (GVD) in the cavity. This is a straightforward step, keeping in mind the fact that we have obtained this term from the Taylor expansion of the wavevector around the central frequency. One should note that the actual value of the GVD is given with k Equations (A7) -(A12) combined with the wave equations (A22) make a complete system that is used in this work.

Appendix B: Frequency modulated comb -master equation
Based on the recent experimental data that can be found in the main body of the paper, self phase-locked lasers give rise to frequency modulated (FM) frequency combs with a particular linearly chirped frequency output. The aim of this work is to try to shed some light on the mechanisms that are responsible. As already mentioned, equations (A7) -(A12) and (A22)) give an accurate quantitative model that is valid in the general case. However, they represent a coupled system of differential equations. Such a system is not very helpful if one aims for acquiring an intuitive understanding of the laser intermode dynamics that lead to a formation of an FM frequency comb. Furthermore, changing one parameter influences several equations and it is not clear what is the dominant effect that emerges as a consequence. Hence, it is preferable if it is possible to eliminate some of the equations and simplify the system. We will do so by considering that the gain medium possesses fast dynamics. This will allow us to obtain a model based only on the wave equation containing additional terms -a master equation. For a fast gain medium, it serves as a quantitative tool almost as good as the full model provided previously. However, concerning the qualitative abilities, it is far more superior, giving insight into how different terms affect the process of self phase-locking of the laser.
We can start with neglecting the lower level population, assuming that the extraction from the lower level is efficient, as well as the pumping to the upper level. Lifetimes T 1 and T g can be introduced for the static population n u0 and the grating part of the population n u2 , respectively, as: Upon dividing the polarization into contributions from n u0 and n u2 as σ ± = σ 0± +σ 2± ,the system can be written as: n c Let us first eliminate n u0 . Its value is strictly real and has no contribution to the the phases of the field envelopes, so it does not affect the process of phase locking directly. Hence, one could expect that it is enough to calculate it from equation (B3) by simply employing the adiabatic approximation, e.g. setting the time derivative to zero. We could conclude since in the FM comb, dynamics of the modal phases is much more significant than the dynamics of the modal amplitudes. We can eliminate the polarization σ ± by just considering the σ 0± contribution which we calculate after making the adiabatic approximation in equation (B5). Hence, after replacing σ ± = i µT2 2h n u0 E ± in the equation (B3), we have, after some calculation: where the saturation field E sat has been introduced as E 2 sat = 2h 2 /(µ 2 T 1 T 2 ). Let us now analyze equation (B4). After applying the Fourier transform, bearing in mind that F( ∂ ∂t ) = i(ω − ω 0 ) since we are dealing with the envelope functions, we can write: so that after applying the inverse Fourier transform: We have utilized the assumption that the gain medium is fast in the derivation of the last relation when we made the approximate step in the Taylor expansion of (1+i(ω− ω 0 )T nu2 ) −1 . We can again eliminate the polarization σ ± by replacing σ ± = i µT2 2h n u0 E ± in the equation (B9): In the last equation it was possible to take n u0 outside of the brackets since the value is strictly real and hence its derivative would not affect the phases of the field envelopes. A similar procedure can be done concerning σ 0± . After applying the Fourier and inverse Fourier transform and also keeping terms up to O ∂ 2 ∂t 2 in the Taylor expansion, one can obtain: The reason behind keeping the second derivative in the Taylor expansion will be clear afterwards. After combining equation (B8) with the previous relation, the expression for σ 0± is obtained: Lastly, it is possible to calculate σ 2± from equation (B6) in a similar way: ∂t .
We have kept the term ∂nu2 ∂t since n u2 is a complex value and has impact on both the amplitude and the phase of the field envelope. Combining equations (B8) and (B10) with the previous relation, one can obtain the final expression for σ 2± : The final step is to insert relations (B11) and (B12) into the wave equation (B7). Upon introducing the unsaturated intensity gain factor g 0 in units of m −1 as: we can finally write the master equation for the forward and backward propagating components E ± of the com-plex field envelope: In this way, the system of eight coupled differential equations was reduced down to a system of just two coupled differential equations. They provide a good quantitative tool for a laser which possesses fast gain dynamics. However, even if that is not the case, equation (B14) is still useful as a qualitative tool to gain insight into the underlying physics. Concerning the quantitative analysis for the slow gain medium, one could utilize the Taylor expansion approach only to the polarization equations and implement the starting differential equations for the populations and complex field envelopes.
From here it is possible to simplify the equation even more. Let us represent the envelopes of the electric field as: where A ± are the strictly real amplitudes and φ ± are the phases of the forward and backward propagating components of the complex field. For simplicity, we can also introduce the saturated intensity gain g(P ) as: where P = P + + P − = |E + | 2 + |E − | 2 is the normalized power and P sat = E 2 sat . Inserting equations (B15) and (B16) into (B14) and grouping real and imaginary terms together yields equations for the amplitudes and phases separately: We can give some comments now about the influence of different terms. From the equation for the amplitude (B17) it is clear why we needed to keep the second derivative in the Taylor expansion in the derivation of σ 0± in the equation (B11). It can be easily seen that it is the only term that gives feedback from the phases φ ± to the amplitudes A ± (if we neglect the dispersion) through , so leaving it out results in chaotic behavior, since the amplitudes become independent of the phases. Furthermore, it is the largest in value, and hence has the largest impact, compared to the other second derivative terms that would emerge from Taylor expansions of σ 2± or n u2 .
Moreover, SHB has a crucial role in the intermode laser dynamics, since it is the main effect that induces multimode lasing. It decreases the gain competition between cavity modes and allows a large number of side modes to overcome the lasing threshold. We verified this easily by setting the SHB grating term in the population n u2 = 0, which consequently makes σ 2± = 0. This alteration results always in single mode laser operation. Concerning the master equation (B14), the term that dominantly describes SHB effect is (T 2 + T nu2 )E ± E ∓ ∂E * ± ∂t given in the second square bracket. Its exclusion causes a single mode solution to arise. This can be understood better by analyzing the equations for the amplitudes and the phases. In equation (B17), this term is represented by (2T 2 +T nu2 ) ∂A∓ ∂t A ± A ∓ which gives the contribution from the amplitude derivative of the opposite propagating envelope. It is even more clear when looking at the phase equation (B18). There, the SHB term is T nu2 A 2 ∓ ∂φ∓ ∂t and it is the only term that results in the crosstalk between the phases φ + and φ − in the corresponding equations.
Equations (B17) and (B18) are valid both for amplitude and phase modulated laser. However, in an FM comb, one can assume that the laser dynamics are dominantly described with the phase equation (B18), since it is the phase that is being modulated. The amplitude is approximately constant and is less important. From this follows that the derivatives of the phases φ ± are more significant than the derivatives of the amplitudes A ± . That means that we can write simplified equations in the case of an FM frequency comb, disregarding multiple terms from equations (B17) and (B18), which we have confirmed via numerical simulations: n c Alternatively, we can express the last two equations using the power P ± instead of the amplitude if we write E ± = P ± e iφ± : We will first comment the power equation (B21). One sees that the entire second square bracket in the equation (B17) is left out. Terms in that bracket originate from SHB and their influence on the power (amplitude) can be neglected, bearing in mind that we are analyzing the dynamics of an FM comb which is governed by the phase evolution. The term proportional to ( ∂φ± ∂t ) 2 is the only one that gives a feedback from the the phase to the power (amplitude) and is necessary in the case of phase-locking in the presence of the Kerr nonlinearity. Intriguingly, it turns out that it can be left out in the case of phase-locking with GVD. Furthermore, the power can be calculated from the simple equation ±∂P ± /∂z = (g(P )−α w )P ± and then plugged back in the phase equations. This yields a time independent power with a completely constant output.
We will now turn to the equation (B22). The first term in the square brackets leads only to a single mode solution since it dampens the oscillations. The second term in the square brackets is due to SHB and it links together the phases of the opposite propagating fields. This leads to multi-mode operation, proving again that SHB is an absolute necessity. Next, it is seen that both the Kerr nonlinearity and the group velocity dispersion dominantly influence the phase, and not the amplitude of the field. The GVD term proportional to k shapes the phases into a parabola, which corresponds to linear intermodal phase differences and a chirped frequency. The sign of k determines whether the parabola is convex or concave, giving rise to a chirp in one direction or another.
Concerning the Kerr influence, we can give a following explanation. The second term in the square brackets in the equation (B21) is a source term, which drives the power to oscillate in a similar manner as ( ∂φ± ∂t ) 2 . Then it is clear that the Kerr term, which proportional to β in the equation (B22), has a similar impact as the GVD. original state. As a conclusion, carrier grating induced by SHB is necessary (but not sufficient) for FM phaselocking with a linear chirp.