6.3 Dynamics of QCLs

The linear stability analysis, however, predicts only the instability threshold and does not describe the dynamics of the laser. To investigate the dynamics of QCLs above the instability threshold, the Maxwell-Bloch equations are solved numerically. The effect of the SA is modeled as the intensity-modulated optical field amplitude in the standard Maxwell-Bloch equations [128]:

        c        c  iμP       c
∂tE  = - --∂zE -  -------- - ---(l0 - ¯γ|E |2)E,
        n        nℏl0Dth    2n
      -iμ-       P--
∂tP  = 2 ℏDE  -  T2,
       D  - D    iμ
∂tD  =  -p------+ ---(E*P  - c.c.).
         T1       ℏ
(6.8)
E and P are the envelopes of the normalized electric field and polarization, respectively, D represents the normalized average population inversion, Dp is the normalized steady-state population inversion, Dth is the lasing threshold value of Dp for γ = 0, proportional to the pumping factor pf (pf 1), l0 is the linear cavity loss, and γ = 2γ∕μ2.

We employ a finite-difference discretization scheme to find the evolution of electric field, polarization, and population inversion in the spatial and time domain. A periodic boundary condition is applied to model a ring cavity. The equations are discretized to calculate the state of the system in each step from the state in the previous step. This procedure is shown in Fig. 6.7.


PIC

Figure 6.7: The m × n mesh points for the finite difference approximation. Here m = 8 and n = 6. The boundary conditions are shown by green dots and initial condition by pink dots. E(m,n+1), a sample step in time domain, is shown with blue dot and its value is dependent on the previous steps, the orange dots.

The number of mesh points play an important role in preventing unreasonable results and achieving consistency. The number should be large enough to make sure there is not any un-inspected variation. The x - t plane is divided into a m × n mesh points as shown in Fig. 6.7. The number of mesh points are m = 500 and n = 105 which correspond to the grid sizes of Δx = L∕m = 12 μm and Δt = Δx∕c = 120 fs. The periodic boundary conditions for the dynamical variables of the QCL (electric field E, polarization P, and population inversion D) are E(0,n) = E(L,n), P(0,n) = P(L,n), and D(0,n) = D(L,n).

The derivatives of variables in time and space are achieved by Taylor expansion about the grid point (m,n).

tE(m,n) = -E-(m,-n-+--1) --E-(m,-n)
          Δt -1-
2Δt 2
∂-E-(m,n-)-
   ∂t2 (6.9)
zE(m,n) = -E (m,n ) - E (m, n - 1)
----------Δz----------- + 1
2-Δx∂2E (m, n)
---∂z2----- (6.10)
Similar expressions are written for tP and tD. Starting with the left hand-side of Eq. 6.9 and 6.10 for the electric field, we have

Ė + c∂E--
∂z = -E(m,-n-+-1-) --E-(m,-n)
          Δt -1-
2Δt 2
∂-E-(m,-n-)-
   ∂t2
-E(m, n ) - E (m, n - 1)
-----------------------
         Δz + 1
--
2c2Δt2∂2E (m, n)
------2----
    ∂z =
E-(m,-n +-1) --E-(m---1,n)-
            Δt + 1-
2Δt[                        ]
 c2∂E-(m,-n)-- ∂E--(m,-n-)
      ∂z2          ∂t2
= k(P - E). (6.11)
The expression Ė denotes tE. Selecting the second equality helps to eliminate the terms with first derivative
k(P - E) = E-(m,-n-+--1) --E-(m---1,n)-
            Δt + 1-
2Δt[   2             2        ]
 c2∂-E-(m,-n)-- ∂-E-(m,-n)-
      ∂x2           ∂t2 (6.12)
Solving Eq.6.12 for E(m,n + 1) gives
E(m,n + 1) = kΔt(P - E) (6.13)
+ 1
--
2Δt2[                          ]
 ∂2E (m, n )   2 ∂2E (m, n )
 ------2----- c  -----2-----
     ∂t             ∂x + E(m - 1,n)
In this way, we have approximated the derivatives up to the order Δt2.

From Eq. 6.13 we have

  2
∂--E(m,-n-)
    ∂t2 = Ë = k( -Ė) - c  2
-∂-E-
∂t∂z (6.14)
 2
∂-E-(m,-n-)-
   ∂x2 = k-
c(          )
  ∂P--- ∂E--
  ∂z    ∂z-1-
c 2
∂-E--
∂z∂t (6.15)
If we multiply Eq. 6.15 by c2 and then subtract it from Eq. 6.14
∂2E (m, n)
---∂t2----- - c2∂2E (m, n)
---∂x2----- = k( -Ė) - ck(           )
  ∂P    ∂E
  ∂z--- -∂z-. (6.16)
As
= (1∕T2)(DE - P),
and
Ė = k(P - E) - c∂E--
∂x,
the right side of the Eq. 6.16 can be simplified to
∂2E (m, n )
------2----
    ∂t - c2∂2E (m,n )
-----2-----
   ∂x =
k[                                    ]
  (1 ∕T )(DE  - P ) - k(P - E ) + c∂E--
      2                           ∂x + ck(          )
  ∂E--  ∂P--
  ∂x    ∂x =
k[(1∕T2)(DE  -  P) - k(P -  E)] + 2ck∂E
----
∂x - ck∂P
----
∂x =
k[(1∕T )(DE  -  P) - k(P -  E)]
     2 + 2kE-(m,-n)---E-(m---1,n-)
          Δt
- kP (m, n) - P(m  - 1,n)
-----------------------
          Δt =
k[(1∕T )(DE  -  P) - (1∕T  + k)P +  kE ]
     2                  2 + 2kE(m,-n-) --E-(m---1,n)-
          Δt
- kP (m, n) - P(m  - 1,n)
----------Δt----------- (6.17)
where we used the following relationship for the speed of wave propagation in vacuum
c = Δx
----
Δt,
Therefore, evolution of electric field in time and space is given by
E(m,n + 1) = E(m - 1,n) + kΔt(P - E) (6.18)
+ 1
--
2Δt2k[(1∕T2)DE  - (1∕T2 )P - kP  + kE ]
+ kΔt[E(m,n) - E(m - 1,n)]
-1-
2kΔt[P(m,n) - P(m - 1,n)],
rearranging the equation results in:
E(m,n + 1) = (1 - kΔt)E(m - 1,n) (6.19)
+ (                          )
  1k Δt - 1-kΔt2(1∕T  +  k)
  2       2          2P(m,n)
+ ( 1                  1           )
  -(kΔt )2E (m, n) + -k(1∕T2 )Δt2
  2                  2D(m,n)E(m,n)
+ (      )
  1
  -k Δt
  2P(m - 1,n).
Following a similar procedure, the equations for P and D are obtained. The details and the full derivations of this method can be found in  [224].

Because we have neglected noise in our treatment, the solutions would remain |E| = |P| = 0 for all times. In order to get the laser started, one must therefore assume a small initial disturbance of the electric field, e.g., a Gaussian  [224]

E(z, 0) = 0.1e-100(z
L-1
2)2 . (6.20)
The initial inversion D is set to the threshold pumping pf which is the minimum injection current needed to start the lasing action.
P(z, 0) = 0, D(z, 0) = pf. (6.21)
Figure 6.8 indicates the transient buildup of the intensity for a ring cavity without the saturable absorber effect. After a few round trips, E,P, and D reach the approximate periodic steady state. To investigate the saturable absorber effect, we recalculate the variables including the saturable absorber term in the Maxwell-Bloch equations, see Eq. 6.8.

PIC

Figure 6.8: Intensity as a function of time. In this graph, pf = 16,L = 6×10-3m,k = 1(10T 2), and T1 = 2T2. After a few round trips, the variables saturate to the approximate continuous wave solution in which the intensity is transformed from the initial Gaussian to a the stable operation mode.

For the parameters corresponding to the optimized 3QW QCL including the SA, the lasing instability appears as the rise of the side modes with the increase of the SA coefficient. The energy in the Rabi sidebands can change either discontinuously or continuously at the RNGH instability threshold [126]. In lasers with slow gain recovery time, the transition in the Rabi sidebands is discontinuous [225], however, because of the fast gain recovery time in QCLs, Rabi sidebands continuously grow around the central cw mode, see Fig. 6.9. Furthermore, more Rabi side modes appear with the increase of the pumping strength.