Frequency combs induced by optical feedback and harmonic order tunability in quantum cascade lasers: supplementary materials

We study the quantum cascade laser (QCL) dynamics by using a full set of effective semiconductor Maxwell–Bloch equations (ESMBEs) for the Fabry–Perot (FP) configuration.Silvestri20; Silvestri22; SilvestriCLEO This model is based on a phenomenological expression for the optical susceptibility of a QCL active medium, and it encompasses the main properties of semiconductor materials such as a non null $\alpha$ factor, a dependence of the susceptibility from the density of carriers, asymmetric gain and refractive index spectra, and nonlinearities at each order, such as Kerr effect and four-wave mixing. In the case of FP cavity it also accounts for SHB through the inclusion of a carrier grating. The ESMBEs read:

	$\displaystyle\frac{\partial E^{+}}{\partial z}+\frac{1}{v}\frac{\partial E^{+}}{\partial t}$	$\displaystyle=$	$\displaystyle-\frac{\alpha_{\mathrm{L}}}{2}E^{+}+gP_{0}^{+},$		(S.1)
	$\displaystyle-\frac{\partial E^{-}}{\partial z}+\frac{1}{v}\frac{\partial E^{-}}{\partial t}$	$\displaystyle=$	$\displaystyle-\frac{\alpha_{\mathrm{L}}}{2}E^{-}+gP_{0}^{-},$		(S.2)
	$\displaystyle\frac{\partial P_{0}^{+}}{\partial t}$	$\displaystyle=$	$\displaystyle\mathrm{\pi\delta_{hom}}(1+i\alpha)\left[-P_{0}^{+}+if_{0}\varepsilon_{0}\varepsilon_{\mathrm{b}}\left(1+i\alpha\right)\left(N_{0}E^{+}+N_{1}^{+}E^{-}\right)\right],$		(S.3)
	$\displaystyle\frac{\partial P_{0}^{-}}{\partial t}$	$\displaystyle=$	$\displaystyle\mathrm{\pi\delta_{hom}}(1+i\alpha)\left[-P_{0}^{-}+if_{0}\varepsilon_{0}\varepsilon_{\mathrm{b}}\left(1+i\alpha\right)\left(N_{0}E^{-}+N_{1}^{-}E^{+}\right)\right],$		(S.4)
	$\displaystyle\frac{\partial N_{0}}{\partial t}$	$\displaystyle=$	$\displaystyle\frac{I}{eV}-\frac{N_{0}}{\tau_{\mathrm{e}}}+\frac{i}{4\hbar}\left[E^{+*}P_{0}^{+}+E^{-*}P_{0}^{-}-E^{+}P_{0}^{+*}-E^{-}P_{0}^{-*}\right],$		(S.5)
	$\displaystyle\frac{\partial N_{1}^{+}}{\partial t}$	$\displaystyle=$	$\displaystyle-\frac{N_{1}^{+}}{\tau_{\mathrm{e}}}+\frac{i}{4\hbar}\left[E^{-*}P_{0}^{+}-E^{+}P_{0}^{-*}\right].$		(S.6)

where $E^{+}(z,t)$, $E^{-}(z,t)$ are the forward and backward electric fields, $P_{0}^{+}$, $P_{0}^{-}$ are the forward and backward polarization terms, $N_{0}$ is the zero-order density of carriers, and $N_{1}^{+}$, $N_{1}^{-}$ are the first-order terms of the density of carriers, reproducing the carrier grating due to SHB; $v$ is the group velocity, $\alpha_{\mathrm{L}}$ is the loss coefficient, $\alpha$ is the linewidth enhancement factor, $\delta_{\mathrm{{hom}}}$ is the homogeneous part of the gain bandwidth, $f_{0}$ is the differential gain, $\varepsilon_{0}$ is the dielectric permittivity of the vacuum, $\varepsilon_{\mathrm{b}}$ is the relative dielectric constant of the QCL active medium, $I$ is the driving current of the laser, $V$ is the volume of the active region, $\tau_{\mathrm{e}}$ is the carrier lifetime, and the coefficient $g$ is given by:

$$g=\frac{-i\omega_{0}N_{\mathrm{p}}\Gamma_{\mathrm{c}}}{2\varepsilon_{0}n_{\mathrm{r}}c}\,,$$

(S.7)

where $N_{\mathrm{p}}$ is the number of stages of structure, $\omega_{0}$ is the cold cavity angular frequency closest to the gain peak and it is used as a reference frequency, $\Gamma_{\mathrm{c}}$ is the optical confinement factor, $c$ is the speed of light in the vacuum, and $n_{\mathrm{r}}$ is the effective background refractive index of the medium. We remark that $\delta_{\mathrm{hom}}$ is related to the polarization dephasing time $\tau_{\mathrm{d}}$ through the equation $\delta_{\mathrm{hom}}=\frac{1}{\pi\tau_{\mathrm{d}}}$.Silvestri20
We want to include the effect of the optical feedback in the boundary conditions of the model. Therefore, we start by adding a term $E_{\mathrm{fb}}(t)$ to the free-running boundary conditions, corresponding to an external target placed at distance $L_{\mathrm{ext}}$ from right facet of the QCL:

	$\displaystyle E^{-}(L,t)$	$\displaystyle=$	$\displaystyle\sqrt{R}E^{+}(L,t)+E_{\mathrm{fb}}(t),$		(S.8)
	$\displaystyle E^{+}(0,t)$	$\displaystyle=$	$\displaystyle\sqrt{R}E^{-}(0,t),$		(S.9)

where $L$ is the QCL cavity length and $R$ is the reflectivity of each facet of the QCL. We want to determine $E_{\mathrm{fb}}(t)$. In the frequency domain we can write the reflection coefficient due to the external feedback in our configuration as:

$\displaystyle r(\omega)=r_{\mathrm{ext}}t_{\mathrm{L}}^{2}e^{-i\Phi_{\mathrm{ext}}}e^{-i(\omega-\omega_{0})\tau_{\mathrm{ext}}}$

(S.10)

where $\omega$ is the angular frequency, $t_{\mathrm{L}}$ is the field transmissivity for both outgoing and incoming field, $\Phi_{\mathrm{ext}}=2L_{\mathrm{ext}}\omega_{0}/c$, $\tau_{\mathrm{ext}}=2L_{\mathrm{ext}}/c$, $r_{\mathrm{ext}}$ is the frequency-independent reflectivity of the external target. If we calculate the Fourier transform of the forward field $E^{+}(L,t)$, we multiply it by $r(\omega)$ given by Eq. (S.10), and then we anti-transform the obtained product, we retrieve the following expression for the feedback term.Rimoldi22

$\displaystyle E_{\mathrm{fb}}(t)=r_{\mathrm{ext}}t_{\mathrm{L}}^{2}E^{+}\left(L,t-\tau_{\mathrm{ext}}\right)$

(S.11)

Since Eqs. (S.1)–(S.6) have been derived by assuming the frequencies of the multimode fields are referenced with respect to the central frequency $\omega_{0}$ (see ref.Silvestri20), we write Eq. (S.12) consistently with this hypothesis:

$\displaystyle E_{\mathrm{fb}}(t)=r_{\mathrm{ext}}t_{\mathrm{L}}^{2}E^{+}\left(L,t-\tau_{\mathrm{ext}}\right)e^{-i\omega_{0}\tau_{\mathrm{ext}}}$

(S.12)

If we also take into account the losses in the external cavity $\epsilon_{\mathrm{L}}$, and the ones due to the reinjection $\epsilon_{\mathrm{S}}$, we obtain the boundary conditions Eqs. (1)–(2) presented in the main manuscript.

The values of the ESMBEs parameters exploited in this article are presented in Table I, and they are typical values for a THz QCL.SilvestriReview; PiccardoReview They have been used to study the free-running case in ref. Silvestri22, where two comb regions characterized by fundamental and second order harmonic combs were found.

The numerical study which leads to Fig. 2 presented in the main manuscript, is conducted according to the following points. The laser is initially configured in free running operation, i.e. we solve the ESMBEs by adopting the boundary conditions with $\epsilon=0$, and the value of the bias current is chosen in order to have single mode emission; in particular, we set $I=1.08I_{\mathrm{thr}}$. Once the QCL output reaches a stable single mode emission, at $t=40$ ns the external target is introduced (Fig. S.1) by integrating the ESMBEs with the boundary conditions Eqs. (1)–(2), for fixed values of $L_{\mathrm{ext}}$ and $\epsilon$. Each simulation has total duration of $1\mu s$ and is performed for a pair ($\epsilon$, $L_{\mathrm{ext}}$); the set of values of $\epsilon$ are chosen between 0.01 and 0.09 with step 0.01, and between 0.1 and 1 with step 0.1; the values of $L_{\mathrm{ext}}$ are selected in four regions with a different step for each region: region 1 corresponds to the short cavity regime, and includes the values between 1 mm and 9 mm, with step 1 mm; region 2 includes the values between 1 cm and 9 cm, with step 1 cm; region 3 corresponds to $L_{\mathrm{ext}}$ between 10 cm and 90 cm with step 10 cm; in region 4 we consider the long cavity regime, with the values 1 m, 2 m, and 3 m. This choice allows us to study the laser dynamics in different ranges of $L_{\mathrm{ext}}$ with a limited number of simulations, providing a general overview of the laser behaviour under optical feedback.
The set of ESMBEs (S.1)-(S.6) with the boundary conditions (S.8)-(S.9) have been integrated using an optimized finite difference algorithm discretizing in both time and space, described in detail in the Appendix of Ref. Bardella17. The time step used in the simulations is $dt=20~{}$fs, which corresponds to a space step $dz=1.8~{}\mu$m. This leads to a duration for a single typical $1~{}\mu$s long simulation of 5 hours on our Intel$\copyright$ Xeon$\copyright$ W-2145 CPU 3.70 GHz processor, with 64 GB RAM. However, we specify that the duration of the transients preceding the attainment of the steady state condition varies depending on the feedback parameters. This implies that simulations have slightly different durations from each other, and the values of 1 $\mu$s and 5 hours are therefore indicative (average) values.

We would like to highlight that the study described in Section IV of the main manuscript follows the same procedure, with one difference: the bias current value is set to $I=1.5I_{\mathrm{thr}}$, which corresponds to a self-starting optical frequency comb regime. Therefore, we first simulate the free-running laser ($\epsilon=0$), allowing it to reach a steady state comb emission. Subsequently, we switch on the feedback and reproduce the map depicted in Fig. 5(a).

Finally, after solving the ESMBEs for each pair ($\epsilon$, $L_{\mathrm{ext}}$), we classify the obtained dynamical regimes following a rigorous procedure described in the following subsection.

In this subsection we show that an OFC generated in regime of short cavity starting from a free-running CW, is sensitive to the phase accumulated in the external cavity. We shift external cavity by a quantity $\Delta L_{\mathrm{ext}}$ in the same order of the QCL central wavelength $\lambda=100~{}\mu$m. In particular we consider integer multiples of $\lambda/4$. As example, in Fig. S.2 we show how the fundamental comb of Fig. 3(d) in the main manuscript evolves if we introduce a shift $\Delta L_{\mathrm{ext}}=n\lambda/4$, with n$=1,2,3,4$, reporting alternance CW-OFC with period given by $\lambda/2$. We understand, therefore, that a $\pi$ phase accumulated in the external cavity (corresponding to $\lambda/4$ and $3\lambda/4$) leads to a single-mode emission, while a $2\pi$ phase shift does not affect the comb emission. Therefore, the system is phase-sensitive, and this can be exploited for applications such as the detection of microscopic vibrations, or displacements.

Furthermore, we want to investigate how the comb emission is affected by a fine tuning of the external cavity length in the long cavity regime. If we keep $\epsilon^{2}=36~{}\%$ as in the case of Fig. S.2, but we consider $L_{\mathrm{ext,0}}=70~{}$cm, we report a mixed state emission. If also in this case we introduce a shift $\Delta L_{\mathrm{ext}}=n\lambda/4$, with n$=1,2,3,4$ of the external cavity length, we report a mixed state for all the considered values of $\Delta L_{\mathrm{ext}}$, as shown in Fig. S.3.

At this point we want to examine how a feedback-induced comb is affected by a fine tuning of the external cavity length in the long cavity regime. For this reason, we consider $L_{\mathrm{ext,0}}=30~{}$cm and $\epsilon^{2}=1~{}\%$, which corresponds to a frequency comb emission in the map of Fig. 2(a) in the main manuscript. By repeating the study performed for the two previous cases, we find that the QCL emits a comb for all the considered values of $\Delta L_{\mathrm{ext}}$, as depicted in Fig. S.4.

We can notice that also in this case the phase shift affects the waveform, but the type of dynamical regime (a comb in this case) is still the same for all the simulations. We have repeated this numerical experiment and verified the occurrence of this behaviour for other combinations of parameters, such as $L_{\mathrm{ext,0}}=20~{}$cm, $\epsilon^{2}=4~{}\%$ and $L_{\mathrm{ext,0}}=20~{}$cm, $\epsilon^{2}=1~{}\%$, where comb emission is obtained for all the considered values of external cavity length shift $\Delta L_{\mathrm{ext}}=n\lambda/4$, with n$=1,2,3,4$; for the combinations of$L_{\mathrm{ext,0}}=40~{}$cm, $\epsilon^{2}=4~{}\%$, $L_{\mathrm{ext,0}}=1~{}$m, $\epsilon^{2}=36~{}\%$, mixed state emission is obtained for all the considered values of external cavity length $\Delta L_{\mathrm{ext}}=n\lambda/4$. Therefore, for the case $I=1.08I_{\mathrm{thr}}$ (free-running CW as initial condition) these results show that in the short cavity regime an alternance between comb states and single mode emission occurs by performing a fine tuning of the external cavity on the $\lambda$ scale. On the contrary, if we consider the long cavity regime, a variation in the length of the external cavity on the wavelength scale does not alter the type of dynamic regime emitted by the laser, although it implies a slight variation in the obtained waveform.

We estimated the critical value of $\alpha$ for which irregular dynamics is observed, considering the case $L_{\mathrm{ext}}=1$ cm, $\epsilon^{2}=36~{}\%$, $\tau_{\mathrm{e}}=5$ ps, and simulating the QCL dynamics for the values of $\alpha$ between -0.1 and 0.7, with a step 0.1, and the other parameters fixed as in Table I. We chose this case because for $\alpha=-0.1$ we have a comb regime, while for $\alpha=0.7$ we report irregular dynamics, so that we are interested to estimate the value where the transitions between these two types of regimes occur. For a rigorous estimation, we calculate the comb indicators $M_{\sigma_{P}}$ and $M_{\Delta\Phi}$. In Fig. S.5 we show the plot of the OFC indicators for different values of $\alpha$ for the mentioned case, and we observe a clear transition from comb to unlocked dynamics for $\alpha=0.5$, a value which is compatible with THz QCLs.SilvestriReview; PiccardoReview We verified that the unlocked states reported for $\alpha\geq 0.5$ are actual irregular regimes and not mixed states (which are also unlocked according to the definition of the indicators, even if they present regular amplitude modulations with the periodicity given by the EC roundtrip). We remark that we referred to these states as "irregular" rather than "chaotic", because we did not assess their level of chaos using methods such as analyzing Lyapunov exponents, which would be necessary for accurately labeling them as "chaotic".
An example of these regimes, corresponding to $\alpha=0.5$, is reported in Fig. S.6. It can be observed that the output power exhibits irregular modulations (Fig. S.6(a)), which do no present the periodicity given by the roundtrip of the external cavity. For this reason, they can not be classified as mixed states. If we look at the zoom on the first beatnote in the power spectrum (Fig. S.6(b)), we notice that this is split in several equidistant peaks spaced by the difference between FSR and $\mathrm{FSR_{ext}}$. In this case the nominal value of FSR is $20.8~{}$GHz, $\mathrm{FSR_{ext}}=15~{}$GHz, and the found spacing in Fig. S.6(b) is $5.8~{}$GHz. This shows that these irregular states arise from the beating between the external and laser cavity modes. By looking at the maps of Figs. 4(c)-(d) in the manuscript, we notice that they occur in a region where the FSR and $\mathrm{FSR_{ext}}$ have the same order of magnitude. We understand, therefore, that in presence of high value of both $\alpha$ and feedback strength, if FSR and $\mathrm{FSR_{ext}}$ have the same order magnitude we do not obtain a comb emission but an unlocked dynamics as described.

We consider $I=1.5I_{\mathrm{thr}}$, corresponding to a free running comb emission. We add the feedback and, for some pairs ($L_{\mathrm{ext}}$, $\epsilon^{2}$) we observe still emission of an OFC, as depicted in the map of Fig. 5(a) in the main manuscript (see the light blue region). Let us consider a pair in the light blue region, i.e. $L_{\mathrm{ext}}=7~{}$mm, $\epsilon^{2}=1~{}\%$, which is also one of the cases displayed in Fig. 5(b) of the main manuscript. We perform a fine tuning of the EC length on the wavelength scale for this comb solution, reporting some differences with respect to the $I=1.08I_{\mathrm{thr}}$ case discussed in Sec. S.3. In fact, we observe a passage from a comb regime (Fig. S.7(a)) to irregular dynamics (Fig. S.7(b)) when the EC length shift $\Delta L_{\mathrm{ext}}$ increases from 0 to $\lambda/4$, and then a return to comb for $\Delta L_{\mathrm{ext}}=\lambda/2$ (Fig. S.7(c)). Then, we again obtain irregular dynamics for $\Delta L_{\mathrm{ext}}=3\lambda/4$, and a comb for $\Delta L_{\mathrm{ext}}=\lambda$. Therefore, our numerical simulations show that the comb solutions are obtained with periodicity $\lambda/2$.
If we consider a a case in the long cavity regime ($L_{\mathrm{ext}}=60~{}$cm, $\epsilon^{2}=1~{}\%$), we observe that the comb states alternate with mixed states, characterized by a multi-peaked beatnote, as shown in Fig. S.8

This scenario is the result of the nonlinear competition between the EC and laser cavity modes in the QCL medium. These results reproduce the experimental findings presented in ref.Liao22, testifying the reliability of our simulator.

In this section we replicate the study of Sec. III in the main manuscript for the case of a mid-IR QCL. We consider values of $\delta_{\mathrm{hom}}$, $\alpha$ factor, central emission frequency $\nu_{0}$, and carrier lifetime $\tau_{\mathrm{e}}$ that are typical of mid-IR QCLs, as indicated in Table 2. For the ESMBE parameters not specified in Table 2 we utilize the values provided in Table 1. The value bias current value is $I=1.08~{}I_{\mathrm{thr}}$, and corresponds to single mode emission in free-running operation, in analogy with the study presented in Fig. 2 in the main manuscript.
The feedback diagram is shown in Fig. S.9.

The higher value of gain bandwidth and $\alpha$ with respect to the THz QCL case of Fig. 2 leads to a larger number of modes involved in the competition inside the active medium. In fact, the phase-amplitude coupling provided by $\alpha$ favors the multimode dynamics of the laser,SilvestriReview; Opacak2019; Silvestri20 which can also manifest itself more easily in presence of a larger gain profile. For this reason, the dynamical scenario presented in the feedback map of Fig. S.9 is characterized by a large number of multimode states, both locked and unlocked. A comb region is observed for $4~{}\%<\epsilon^{2}<10~{}\%$ in the short cavity range ($L_{\mathrm{ext}}<1~{}$cm). However, for $\epsilon^{2}>10\%$ we report a majority of irregular regimes if $L_{\mathrm{ext}}<10~{}$cm. The high level of feedback, in combination with a larger $\alpha$ and more modes involved, implies an enhanced competition between external and laser cavity modes, and prevents the locking of the involved optical lines. Then, for $L_{\mathrm{ext}}>10~{}$cm we observe mixed states, in analogy with the THz QCL diagram of Fig. 2.
Furthermore, we notice that the parameter configuration discussed in this section differs from the case of Fig. 4(c) in the main manuscript for the larger value of $\delta_{\mathrm{hom}}$. This allows us to highlight the role of the gain bandwidth in the comb/multimode state formation in presence of an external target. In fact by comparing the two maps of Fig. S.9 and Fig. 4(c), we clearly observe that a larger gain bandwidth is crucial in order to have multimode regimes and comb formation in the short cavity range, and compensates the effect of the short carrier lifetime, which tends to keep stable the single-mode emission.mezzapesa2013; columbo2014

Finally, we observe that in the mid-IR QCL case, we report HFCs with order 5 (Fig. S.10(a)), 6 (Fig. S.10(b)), and 7 (Fig. S.10(c)), that were not found with the parameters of a THz QCL (in that case the maximum reported order was 4). This shows that in correspondence of lower values of both polarization dephasing time $\tau_{\mathrm{d}}$ and carrier lifetime $\tau_{\mathrm{e}}$ (and therefore faster carrier and polarization dynamics), HFCs characterized by a faster time scale can be generated in presence of feedback.