**2. Non-equilibrium Green's function method**

The most fundamental tool in the design of THz-QCLs structures and analysis is a numerical package to calculate subband wavefunctions and energies. Because the

subband energy position is critical for the discussion of the parasitic absorption in this work, it needs to estimate the high-lying energy position more precisely. Here, two factors effecting the energy separation between subbands are considered, **a.** In THz-QCLs, the quantum structure contains the layers with thickness of only several nanometers; the non-parabolicity can largely effect on the energy of confined subbands, especially on HCS as it is lifted further away from the bottom of the conduction band; and the high-electric filed operation of THz-QCLs will worse this issue. Here, the band structures are based on three-band Hamiltonian that accounts for the conduction (*c*), light-hole (*lh*), and split-off (*so*) bands, as follows:

$$H = \begin{pmatrix} E\_c(\mathbf{z}) + \mathbf{S}(\mathbf{z}) \frac{\hbar^2 k\_x^2}{2m\_0} & i\sqrt{\frac{2}{3}} \mathbf{P}(\mathbf{z}) k\_x & -i\sqrt{\frac{1}{3}} \mathbf{P}(\mathbf{z}) k\_x \\\ -i\sqrt{\frac{2}{3}} \mathbf{P}(\mathbf{z}) k\_x & E\_{\text{fl}}(\mathbf{z}) + (\mathbf{1} + L(\mathbf{z})) \frac{\hbar^2 k\_x^2}{2m\_0} & \mathbf{0} \\\ i\sqrt{\frac{1}{3}} \mathbf{P}(\mathbf{z}) k\_x & \mathbf{0} & E\_{\text{av}}(\mathbf{z}) + (\mathbf{1} + L(\mathbf{z})) \frac{\hbar^2 k\_x^2}{2m\_0} (\mathbf{z}) \end{pmatrix} \\ (1)$$

where *P* is the interband momentum matrix element related to the Kane energy *Ep* through: *P z*ð Þ¼ ffiffiffiffiffiffiffiffiffiffiffiffi *m*0*Ep*ð Þ*z* 2 q . By comparing the three-band and one-band models for calculating the high-lying subband energy, it finds an approximate 2 meV difference. **b.** the alignment of subbands and also the energy position is also very sensitive to the conduction band offset (CBO) values, especially in short period with tall barriers (here, AlAs% of AlGaAs barrier is 0.3). We follow the latest calibration of CBO based on a machine learning method reported in Ref. [7]. The exact numbers of subbands participating in the transports were controlled by the axial cut-off energy. These subbands were transformed into localized basis modes (reduced real space basis) and used in the NEGF algorithm [26, 27]. The subband energy broadening can play significant roles for estimating the tunneling current (increase the dephasing) and the optical gain (widen the radiation linewidth), especially for THz-QCL studied in this work, as this subband energy broadening (�10 meV) is similar as the photon energy (15 meV). In THz-QCLs, this broadening originates from multiple scattering couplings. Here, the self-energy terms are calculated for all scatterings, including the optical phonons, acoustic phonons, charged impurities, interface roughness, alloy disorder, and electron-electron interactions [28–32]. The critical part of the model is a self-consistent NEGF solver that starts from an initial guess of the Green's functions, the self-energies are then presented roughly, and the Green's functions are again calculated iteratively. Simultaneously, the mean field electrostatic potential is calculated self-consistently (Poisson's equation). Such iterations are performed until convergence is reached. The current density as well as the carrier density distribution is finally displayed. The optical gain or absorption in pairs of intersubband transitions follows linear response theory.
