Terahertz absorption-saturation and emission from electron-doped germanium quantum wells

We study radiative relaxation at terahertz frequencies in n-type Ge/SiGe quantum wells, optically pumped with a terahertz free electron laser. Two wells coupled through a tunneling barrier are designed to operate as a three-level laser system with non-equilibrium population generated by optical pumping around the 1→3 intersubband transition at 10 THz. The non-equilibrium subband population dynamics are studied by absorption-saturationmeasurements and compared to a numerical model. In the emission spectroscopy experiment, we observed a photoluminescence peak at 4 THz, which can be attributed to the 3→2 intersubband transition with possible contribution from the 2→1 intersubband transition. These results represent a step towards silicon-based integrated terahertz emitters. © 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement


Introduction
The development of optical emitters made of materials that can be directly integrated into the silicon foundry processes is a major goal of materials science and photonics [1][2][3][4]. Germanium and the SiGe alloys, which are fully foundry-compatible materials, can be employed to epitaxially grow heterostructures and quantum wells (QWs). If the QWs are doped by electrons or holes, unipolar devices can be realized exploiting their two-dimensional subbands, separated by quantized energy values, as discrete energy levels [5,6]. In particular, radiative intersubband transitions (ISBTs) can be engineered to develop optical emitters with high dipole strengths and photon energy set by the heterostructure design. The photon-emission properties of ISBTs can be exploited to circumvent the poor interband emission properties featured by indirect-bandgap group-IV materials (Si, Ge and their alloys). ISBT optical emitters, however, have been developed up to now only using III-V compound semiconductor heterostructures. Indeed, the electricallypumped quantum cascade lasers (QCL) [7], the optically-pumped quantum fountain lasers (QFL) [8,9], and the hybrid ISB Raman laser [10] have been demonstrated both in the mid-infrared (IR) and in the terahertz (THz) range using III-V compounds, but not yet using group-IV materials.
Preliminary evidence of mid-IR emission from hole-doped (p-type) Si/SiGe QWs has been previously reported [11,12], but design difficulties due to parasitic relaxation channels associated to the multi-band nature of the valence states (light-hole (LH) band, heavy-hole (HH) band, split-off band) have prevented further development. Designs for THz emission from p-type Si/SiGe circumvented some of these problems by using very narrow QWs and high strain to move the LH1 state well above the HH2-HH1 radiative transition to remove parasitic channels, but the high confinement effective mass m* = 0.3 m 0 (where m 0 is the free electron mass) still results in insufficient gain to overcome the waveguide losses [13,14].
Electron-doped (n-type) Ge/SiGe QWs have been regarded as a promising alternative due to the small effective mass m* ∼ 0.13m 0 associated to the fourfold degenerate L-valley minima of the conduction band. From this perspective the n-Ge/SiGe material system differs from the III-V compound semiconductor counterparts, such as n-type GaAs/AlGaAs, used in previous successful demonstrations of QCLs [15,16] and QFLs [8,17], which feature instead a single Γ-valley minimum. The band offset in n-type Ge/SiGe heterostructures is of the order of 100 meV [18], a value suitable for the development of THz emitters. Multi-µm-thick QW stacks are required for waveguide mode confinement in the THz range, but the growth of high-quality multi-µm-thick SiGe heterostructures is challenging due to lattice parameter mismatch between Ge and Si. Previous works [13,[18][19][20][21][22][23][24][25][26][27][28][29] have demonstrated that it is possible to strain-symmetrize the SiGe stacks up to hundreds of repetitions of the heterostructure period, and dope the QWs in a controlled way. As a consequence, reliable THz absorption spectra of ISBTs from straincompensated n-type Ge/SiGe QWs have been obtained, in good agreement with numerical simulations [19,30,31].
The scope of this work is to directly measure THz photon emission from ISBTs of n-type Ge/SiGe asymmetric coupled-quantum-well (ACQW) samples under optical pumping with a tunable THz Free Electron Laser (FEL). This achievement can be regarded as an important intermediate milestone towards the realization of an electrically pumped THz-QCL based on group-IV materials. In fact, if properly designed, ACQWs can operate as three-level laser systems (1: fundamental level, 2: lower laser transition level, 3: upper laser transition level). In previous works on GaAs-based ACQWs, ISBT photoluminescence (PL) spectra could be recorded at the 3→2 and 2→1 transition energies with optical pumping at the 1→3 transition energy [32]. To the best of our knowledge, ISBT-PL experiments have not been performed up to now on Ge/SiGe ACQWs. In our previous time-resolved pump-probe experiments at the FEL, we determined the non-radiative relaxation times of ISB populations in decoupled rectangular QWs [31,33,34] and in ACQWs [35], which are in agreement with a numerical model of ISB population dynamics under optical pumping [35]. In this work, we crucially obtained ISBT-PL emission spectra in the 2-8 THz range, featuring an emission peak centered around 4 THz, which is attributed to the 3→2 ISBT, with also a possible contribution from the 2→1 radiative relaxation. In addition, absorption-saturation by ISBTs at the 1→3 transition with picoseconds recovery times has been clearly observed, which could be relevant for future ultrafast operation of QCLs [36].

Sample description
The heterostructure design targeted a realistic compromise among: long non-radiative relaxation times of the 3→2 transition for maximum PL efficiency; a strong dipole of the 1→3 transition for efficient pump-induced population of level 3; and sufficient wavefunction overlap of levels 1 and 2 for efficient depopulation of level 2. In our ACQW samples, levels 2 and 3 originate from the coupling of the fundamental and first-excited levels of the two isolated wells [30]. Notice that the breaking of the spatial inversion symmetry in ACQWs allows optical transitions between subbands 1 and 3, otherwise prohibited in decoupled QWs by dipole selection rules. Differently from our previous optical pumping experiments with ACQWs [35], here we achieved comparable 1→2 and 1→3 oscillator strengths due to a larger degree of hybridization of levels 2 and 3. This condition enabled us to measure clear absorption signatures of both 1→2 and 1→3 ISBTs by Fourier-Transform IR (FTIR) spectroscopy. In this way we were able to precisely tune the FEL photon energy at the 1→3 absorption resonance at ∼10 THz in order to maximize the pump efficiency. Samples (sketched in Figs. 1(b) and 1(c)) were grown epitaxially on low-impurity concentration Si(100) wafers using ultra-high vacuum chemical vapor deposition (UHV-CVD). The epitaxial layer structure of the active region is made of 20 periods of identical Ge/Si 0.2 Ge 0.8 ACQWs in which a wide Ge well of width w L = 12 nm (only sample S2 has w L = 13 nm) is coupled to a narrow well 5 nm wide through a thin tunnel barrier having a thickness of 2.3 nm. The ACQW modules are separated by a 21 nm-thick Si 0.2 Ge 0.8 spacer (see Fig. 1(e)). Two samples are heavily n-doped in the wide well to N tot = 7 × 10 11 cm −2 for sample S1, and N tot = 9 × 10 11 cm −2 for sample S2. Sample S6 is n-doped to 1 × 10 11 cm −2 and one reference sample was left undoped (S7). The ratio between TM-and TE-polarized transmittance spectra (dichroic transmittance) recorded by FTIR spectroscopy at 10 K (dark-color curves) and 300 K (light-color curves) for the two doped samples (S1 and S2), the lightly doped sample (S6) and the undoped reference (S7). Intersubband transitions at energies E 12,FTIR and E 13,FTIR are clearly visible in S1 and S2, becoming less intense in less doped S6. The undoped sample S7 shows no evidence of ISBTs. The shaded areas highlight the FEL pulse energy range (blue) and the region where the strongest impurity-like transitions occur in the Si substrate (gray, see text). A sketch of the processes involved in (b) saturable absorption and (c) ISBT-photoluminescence in Ge/Si 0.2 Ge 0.8 ACQWs, with the potential profile and the relevant wavefunctions. (d) A sketch of the single-pass 70°slab waveguide used for the optical measurements. The vertical arrow indicates the radiation electric field direction. (e) A representative transmission electron microscopy image of the ACQW structures (sample S7). (f) Numerical simulations of the temporal dynamics of the normalized population n i =N i /N tot for sample S1 under optical pumping. The time envelope of the pump pulse used in the simulations is reported as dotted lines in both panels. Lattice temperature was set to 6 K.
Although integrated waveguide approaches do exist for the on-chip propagation of the longwavelength IR pump beam [37], in this work we employed a standard laboratory configuration, typically used to propagate both IR and THz beams with the electric field perpendicular to the growth plane hence almost parallel to the ISBT dipole vectors [38]. This configuration is a slab-waveguide geometry capable of guiding both pump and emission beams (see Fig. 1(d)), with the sample side facets cut and polished at 70°to the growth plane. The surface on top of the QWs was coated by metal (Ti/Au 10/80 nm), so as to excite a standing wave with transverse magnetic (TM) field and electric field exactly perpendicular to the growth plane at the QW/metal interface (see Fig. 1(d)). The length of the waveguide is about 2.5 mm for the heavily doped samples (S1, S2) and for the undoped reference (S7), in order to realize single-pass propagation of the pump pulse (no multiple reflections), while it is 6.0 mm for the lower doped sample S6 in order to increase absolute IR absorption therein. The waveguide thickness is set by the silicon wafer thickness of approximately 0.5 mm. The 5 mm clear aperture in the cryostat mount sets the actual waveguide width. Due to the lack of infrared-active phonons in both the wafer and the Ge/SiGe epitaxial layers, the background absorption is dominated by hydrogenoid impurity transitions in the much thicker Si wafer, which however do not depend on the polarization, therefore do not appear in linear dichroic spectra. FTIR dichroic transmission measurements are shown in Fig. 1(a) and, in the case of doped samples, they demonstrate two ISBTs of comparable strength, as targeted in the design, centered at the photon energies E 12,FTIR and E 13,FTIR . The ISBT energies correspond well to the values calculated with a Poisson-Schrödinger solver E ij,abs (see Table 1) after the depolarization shift has been taken into account to correct the bare transition energy E ij [30].

Optical pumping
Optical pumping of the 1→3 ISBT was achieved with quasi-monochromatic THz pulses from the FEL tuned at pump photon energies ω p in the 41.9 -48.1 meV interval centered around E 13,FTIR (338 -387 cm −1 , 10.2 -11.7 THz, or 25.8 -29.6 µm). The pulse duration was ∆t p ≈ 10 ps with a Gaussian envelope and the duty cycle was 1.3 × 10 −4 thanks to a high FEL repetition rate of 13 MHz. The maximum peak power at the cryostat position was P in,max = 15 kW. The combination of diffraction losses, electromagnetic mode mismatch at the waveguide input-facet and reflections at the cryostat windows lead to an estimated loss factor for optical coupling into the slab waveguide ∼ 6 × 10 −3 (see Appendix). Therefore, the actual pump intensity in the sample was I p,max = P in,max /πr 2 = 1.1 × 10 4 W/cm 2 , a value still sufficient for pumping a significant fraction of the electrons to level 3 as discussed in the following section. For comparison, using a QCL with peak power around 1 W for the pump beam, the achievable pump intensity would have been three to four orders of magnitude lower than I p,max with the FEL. The ISB population dynamics induced by the pump pulse has been calculated using a selfconsistent energy-balance model based on a rate equation approach [21,35] (see Fig. 1(f) for S1, similar results for the other samples not shown). The electron scattering mechanisms implemented in the model to account for inter-and intra-subband particle and energy fluxes are: electron-phonon (optical and acoustic branches); interface roughness (0.2 nm, with correlation length of 7.0 nm [39]); and ionized impurity scattering. From Fig. 1(f) one sees that the maximum values of the normalized population in level 3 and 2 are comparable, being both around 0.14, although achieved at different delay times. The shape of n 3 (t) closely tracks the shape of the pump pulse, while n 2 (t) is apparently delayed (the maximum of n 2 (t) occurs at 7 ps time delay) and shows a slower decay than n 3 (t). An exponential fit of the two decay tails, performed in a time interval located just after the maximum, gives τ relax,2 ∼ 17 ps, approximating the lifetime of level 2 in the three-level system, and τ relax,3 ∼ 3 ps (τ relax,3 ∼ 6 ps for S2) approximating the lifetime of level 3. Because τ relax,2 > τ relax,3 , population inversion between levels 3 and 2 in our samples can be achieved only during the pump pulse, thanks to the delayed maximum of n 2 (t).
Notice that 3→1 nonradiative relaxation via longitudinal optical (LO) phonon emission is possible at all temperatures in our samples, because E 31 is close to the LO phonon energy of Ge (37 meV) [34]. This fact produces the very short τ relax,3 ∼ 3 ps, which is detrimental to population inversion. To overcome this limitation, longer τ relax,3 could be envisaged in future designs by reducing the wavefunction overlap between levels 1 and 3. This improvement, however, would come at the cost of a lower optical pump efficiency, which is also proportional to the wavefunction overlap between levels 1 and 3. The sustained rise of n 2 (t), also hampering population inversion, is due to the arrival of electrons to subband 2 through elastic scattering from subband 3 via the 1→3→2 sequence. The electrons arriving in subband 2 through this process have high kinetic energies (resulting from their high initial potential energy in 3), and thus produceexcess electron temperature in subband 2 of about 70 K, a value significantly higher than the lattice temperature of 6 K. At short delay times, such a high electronic temperature enables thermally activated phonon-mediated electron scattering events from subband 2 to 1. Since the involved electrons occupy an initial state with large kinetic energy, this intersubband relaxation channel is responsible for a rapid cooling of the subband 2 electrons. In turn, this rapid cooling slower the n 2 (t) population decay rate at larger delay time. . As a result, the maximum of n 2 (t) is similar to the maximum of n 3 (t), and this prevents sustained population inversion. At delay times approaching 15 ps, n 2 (t) decays with characteristic times of few tens of picoseconds, in line with pump-probe data measured in two-level n-type Ge/SiGe QWs featuring E 12 smaller than the LO phonon energy [31,[33][34][35].
In order to check the pumping efficiency in our experimental conditions, we measured the saturation of the 1→3 absorption by increasing step-wise the pump fluence F p = ∆t p P in /πr 2 from zero to the maximum value provided by the FEL with metal mesh attenuators. We simultaneously recorded the incident peak power P in and the transmitted peak power P T (see Fig. 2(a) for details) and we calculated the fluence-dependent transmittance T(F p ) = P T / P in (the fluence being proportional to P in ) for different samples, pump photon energies ω p and temperatures, shown in the plots of Figs. 2(b)-(e). The zero-fluence transmittance T 0 = T(F p ∼0) is calibrated with a Lorentz oscillator model of the transmittance based on FTIR data for each sample, temperature and pump photon energy (see an example in Fig. 2(f)). Two doped samples (S1, S6) and one undoped reference (S7) were investigated, mounted in the sample holder so as to have the FEL radiation in TM-polarization inside the slab waveguide. A nonlinear increase of T(F p ) is clearly seen in Figs. 2(b)-(e) for the doped samples S1 and S6 at all the investigated pump photon energies and temperatures, maybe apart from one case in S6 (red curve in Fig. 2(c), possibly affected by pump depletion due to a longer waveguide length of sample S6). In contrast, the undoped reference sample S7 is the only sample that shows a non-monotonic behavior (see Fig. 2(b)). This is a clear demonstration of the main role played by the bleaching of the ISBTs in the increase of T(F p ) of the doped samples. The lightly doped sample S6 in Fig. 2(c) demonstrates a transmission increase of 10% at F p,max at all ω p , whilst the heavily doped sample S1 in Figs. 2(d) and 2(e) reaches an increase of transmission above 25% for resonant pumping at ω p = 41.9 meV. This is consistent with the effectiveness of heavy doping in increasing the amount of saturable ISBT losses with respect to non-saturable optical transmission losses.
The overall scenario is quantitatively confirmed by our numerical model of the non-equilibrium population dynamics. From the simulated saturated absorption coefficient α 13 (F p ) we calculate the fluence-dependent transmittance T(F p )=exp(-α 13 (F p )L). Here α 13 (F p )=C(N tot f 13 /m* √ ε r )(n 1 -  The relative transmittance change ∆T/T 0 for sample S1 at 6 K (g) directly compared to simulations (h). The dashed line in (g) is the ratio between S1 and S7 data in (b). The zero-fluence FTIR transmittance value T 0 has been subtracted from the data in panel (d) to calculate ∆T. Note that the abscissa of (b-e, g-h) represents the fluence inside the sample corrected for the optical loss factor . n 3 are taken at t = 0 from Fig. 1(f), ∆ω ∼ 10 13 s −1 is the ISBT linewidth, and L is the equivalent optical path length in the QWs. L depends on the number of repetitions of the ACQW structure (here equal to 20), on the internal reflection angle θ in =75°, and on the physical waveguide length [19]. From the simulated transmittance data, we then evaluate, for three different pump photon energies close to the E 13 resonance, the differential transmittance ∆T/T 0 = (T(F p ) -T 0 )/T 0 as a function of F p and we find fair agreement with the same quantity determined experimentally for all the pump photon energies employed (compare Figs. 2(g) and 2(h)). Notice that the numerical model data in Fig. 2(h) does not include any effect extrinsic to the ISBT, therefore it does not include pump-induced absorption producing the apparent "flattening" of experimental curves at high fluence in Fig. 2(g). Indeed, an even better agreement between data in Fig. 2(g) and numerical simulations in Fig. 2(h) is obtained by dividing the ∆T/T 0 data for S1 at ω p = 41.9 meV by the corresponding data for S7 (empty symbols in Fig. 2(g)), thereby normalizing for the pump-dependent silicon wafer transimssion. From the ∆T/T 0 calculated from numerical simulations we obtain a saturation fluence F p,sat,num = 17 × 10 −8 J/cm 2 [40], which we use as a rough estimate of the optimal pump intensity in a ISBT-PL emission experiment.

Terahertz intersubband emission
The analysis of absorption-saturation indicates that a fraction 0.14 of the electrons in the fundamental level 1 is excited to the level 3 under optical pumping with the maximum available FEL pulse fluence F p,max = 11 × 10 −8 J/cm 2 ∼ 0.6 F p,sat,num . We have then used the same fluence for measuring the radiative THz-PL emission of ISBTs. To this aim, we replaced the second pyroelectric detector in Fig. 2(a) with a Michelson interferometer and a liquid-He bolometer ( Fig. 3(a)). The collection mirror after the cryostat was a 50-mm diameter parabolic reflector with a focal length of 100 mm, therefore emission angles up to ± 14°outside the waveguide were collected, which corresponds to ± 4°inside the waveguide according to Snell's law. For this measurement, the most heavily doped sample S2 was selected to maximize the output emission power. The FEL pump photon energy was 44.4 meV with a detuning of ∆ = +4.8 meV from E 13,FTIR = 39.6 meV (see Fig. 1(a)). The detuning was meant to limit pump absorption in the silicon wafer associated to the tail of the impurity transition line at 39 meV. The effect of this small detuning on the pump efficiency is not critical, as observed in the analysis of Fig. 2.
A long-wavelength pass filter (2 mm thick z-cut quartz plate, whose polarization-dependent transmittance is shown in Fig. 3(c)) was used in front of the bolometer to screen it from the direct FEL pump beam. The attenuation of the quartz plate at the pump photon energy of 44 -45 meV is estimated in the range of 10 −8 , making the residual FEL power smaller than the blackbody background. A wire-grid polarizer was inserted before the bolometer to switch between TM-or TE-polarized emission, while the pump was always kept in TM-polarization. The TM and TE emission spectra of sample S2 are shown in Fig. 3(c) on a logarithmic intensity scale and zoomed in Fig. 3(d) on a linear intensity scale. The horizontal scale of the plot of Fig. 3(d) is cut at a photon energy of 10 meV (2.4 THz) because the features appearing in the TM spectrum of Fig. 3(c) below that frequency are tentatively attributed either to coherent synchrotron radiation emission from electron microbunches in the FEL cavity, or to greybody emission from the lattice and/or the electron gas, and are therefore not related to ISBT emission. In Fig. 3(c), an emission peak centered around 16 meV is clearly observed in the TM emission spectrum (green curve) but not in the TE one (purple curve). This value is in good agreement with the estimated 3→2 ISBT energy (E 32 = 14.6 meV). However, since the 2→1 ISBT energy E 21 = 20.3 meV is close to E 32 , also photon emission from this transition may contribute to the observed PL emission peak. Indeed, our numerical simulations indicate that the time-integrated number of photons associated to the 2→1 transition should be higher by a factor of 1.6 than that the one due to the 3→2 transitions (black dashed curve in Fig. 3(c)). To provide an explanation for the fact that the observed PL peak position is closer to the 3→2 ISBT energy, we observe that reabsorption events involving photons emitted across the 2→1 transition, are more likely to occur with respect to those related to 3→2 photons, due to the large carrier density in the fundamental subband. Since in our model we have not taken into account energy dependent reabsorption effects, this mechanism can be invoked to address the observed discrepancy between the experimental and calculated peak energy of the PL spectrum. As an alternative explanation, we notice that in our numerical model, also an increase of the intrasubband cooling rates moves the PL spectral weight towards E 32 . As a matter of fact, to explore the consequences of a possible underestimation of the intrasubband cooling rates, we have simulated the carrier dynamics with an externally fixed and subband independent electron temperature. In this case we obtain that upon lowering Te, the radiative relaxation of E 2 electrons into the fundamental subband slows down, due to the occurrence of Fermi blocking effects. In fact, when considering vertical transition between cold carriers, the final state in the fundamental subband is not empty, because of the large carrier population in the ground state. As a consequence, PL spectrum dominated by photons emitted across the 3→2 transition are predicted for low electron temperatures. To provide evidence in favor of this scenario, we notice that Fermi blocking suppressing the emission of 2→1 photons has been reported also for the quantum fountain structure studied in [32]. Finally, one should consider that, due to the relative small energy difference between E 32 and E 21 it could also happen that the observed PL feature is indeed mainly related to 2→1 radiative transitions as predicted by the model, which however might not sufficiently accurate in reproducing the exact value of the E 21 emission energy. This could be attributed to uncertainties related to the lack of a precise knowledge of the dopant distribution profile (which has a non negligible impact on the subband energy spacing) or to corrections associated to many body effects influencing the emission energy, which have not been included in our calculations.
The emitted ISBT-PL peak power has been estimated from the comparison of the TM emission spectrum (green curve in Fig. 3(c)) with the (attenuated) pump line measured with the same detector (blue curve in Fig. 3(c)). The energy-integrated ISBT-PL emission peak power is found in the range of 10 to 100 µW, where the uncertainty is due to non-linearity of the detector response. Assuming isotropic in-plane emission of the vertical ISB dipole (red shaded circular area in Fig. 3(b)) [41], we obtain a factor 4°/180°= 0.022 for the collected emission and then a calculated peak power emitted in all directions of 0.45 to 4.5 mW, corresponding to a ISBT-PL efficiency of the order of 10 −7 , which is in rough agreement with the ratio between typical radiative and non-radiative lifetimes of THz ISBTs (of the order of 10 µs and 1 ps, respectively) [33,42,43]. The time-averaged CW power of our device (60 to 600 nW) compares well with previous ISBT-PL experiments [42], which were also carried out with pulsed lasers (470 nW emission with 7.8 W pump from CO 2 laser), and with the power emitted by THz-ISBT electroluminescent devices (3 nW in Ref. 1 and 6 nW in Ref. [44]). Our optically-pumped THz emission peak power of 0.45 to 4.5 mW is still significant if compared to optimized electrically pumped THz-QCLs based on AlGaAs/GaAs QWs, which emit up to 100 mW peak power [16], however with a much higher power density (typical QCL ridge width is 0.1 mm, compared to 5 mm of the present device). Note that the typical input electrical peak power of a THz-QCL around 20 W [16] is comparable to the input optical peak power inside our device of 90 W (calculated from the optical peak power of the FEL of 15 kW corrected by the optical coupling loss factor ). In summary, the energy conversion efficiency of our emitter is approximately lower by a factor 10 2 -10 3 than that of a THz-QCL, mostly due to the lack of optical cavity enhancement and amplified stimulated emission in our devices.

Conclusions
In conclusion, we have observed terahertz photoluminescence emission from germanium quantum wells. To this aim, Ge/SiGe heterostructures have been epitaxially grown on silicon wafers realizing a three-level laser system exploiting the so-called "quantum fountain" scheme, in which two asymmetric-coupled quantum wells are doped and pumped by an external optical beam (here, a free electron laser) promoting electrons to the upper laser transition level 3 through the 1→3 inter-subband transition. The required experimental fluence level of the optical pump pulses has been determined by 1→3 absorption-saturation measurements in the range 10-12 THz. Numerical simulations indicate that approximately 14% of the electrons were promoted to the upper laser transition level in the experiment. Photon emission took place around 4 THz at cryogenic temperatures, with a photoluminescence efficiency of the order of 10 −7 . A spectral peak appears only in the TM polarization, in which the electric field direction is normal to the quantum well plane. The peak position at 4 THz matches the inter-subband transition energy between levels 3 and 2, as determined by both band-structure calculations and spectroscopy at equilibrium, thus supporting the assignment to inter-subband photoluminescence emission. These results represent a significant step towards silicon-based integrated terahertz emitters, which have been explored only by numerical simulations up to now [45][46][47].

Appendix A: evaluation of the waveguide optical-coupling-loss factor
The surface-plasmon slab waveguide geometry with input facet at 70 degrees with the quantum well (QW) plane, purposely realized for terahertz transmission and emission experiments (see Fig. 1(d) and Fig. 2(a)), guarantees pure TM polarization of the incoming radiation in the QW region. The deposition of a metal film on top of the QWs ensures that the TM polarization has a nonzero field in the QW region, while the TE polarization serves as the reference with null electric field in the QW region. In this way, by taking linear dichroic spectra with a FTIR spectrometer, one obtains self-referenced absorption spectra of intersubband transitions (ISBTs), where all other absorption contributions are canceled out (see Fig. 1(d)). Similarly, the difference between TM and TE emission spectra in Figs. 3(c)-(d) can be entirely assigned to emission from QWs. As a downside, the edge-coupling of terahertz radiation into the sample is not efficient because of the mismatch between the focal spot diameter of 8.0 mm (FTIR) or 1.0 mm (FEL) and the wafer thickness (0.5 mm). Furthermore, we used a rectangular aperture of 0.3×5.0 mm 2 located in front of the sample edge to eliminate stray light (see Fig. 2(a)). To quantify the terahertz radiation intensity that is lost in edge-coupling, we have measured the absolute transmittance at 30 µm wavelength by FTIR, using the intensity measured with a focal spot diameter of 8.0 (1.0) mm as reference for the transmittance of the slit-sample mounting in FTIR (FEL) experiments. In the FTIR experiment, the absolute transmittance is in the range of 0.1%. In the FEL experiment, the transmittance is estimated to be around 0.9%. A further loss factor comes from the diamond window of the cryostat. At THz frequencies, the diamond/air interface has a reflectance of 0.17.
Considering that the window is slightly wedged to avoid coherent spectral oscillations of the transmittance, the reflection loss factor has to be considered twice because of the two facets of the window, T diam ∼ (1-R) 2 . Therefore, the total loss factor is quantified as follows: (1)

Appendix B: fitting equation for absorption-saturation data
In Fig. 2 we show a fit to the absorption-saturation data in order to extract a direct experimental estimate of the saturation fluence F p,sat,exp . Widely employed analytical models for saturable absorption [48][49][50] do not consider the finite length of the absorber and thus are valid for ultrathin absorbers only. Here, instead, we have to take into account the finite length of the sample and the non-saturable loss mechanisms, therefore we employ the optical model for the fluence-dependent transmittance of a saturable absorber proposed in [51] and recently used in [52]. The fitting equation, reported as continuous curves in Figs. 2(b)-(e) is: The fitting formula accounts for the fact that, even in the absence of ISBT absorption, the transmittance does not reach unity due to reflection losses and absorption in the silicon substrate wafer, through the non-saturable transmittance parameter T ns . The value of T ns = 0.39 has been obtained from FTIR analysis of the undoped sample S7. The finite effective length L appears in the calculation of the linear (i.e. zero-fluence) transmittance T lin , which ideally corresponds to the quantity T 0 measured by FTIR spectroscopy on samples S1, S2 and S6. As reported in Table 2, T lin depends on the doping level, the photon energy and the sample length and varies in the range from 0.12 (for S1 at the ISBT resonance) to 0.35 (for S6 out-of-resonance). The maximum theoretical level of saturable losses for a given sample and pump photon energy is given by the difference T ns -T lin , which is larger for heavier doping and/or at the ISBT resonance. a For S7 the value of the saturation fluence is related to the saturation of impurity transitions in the silicon wafer, and not the saturation of intersubband transitions b The best-fit values F p,abs > 180 × 10 −8 J/cm 2 have large uncertainties because they are more than ten times higher than the fitting range of fluencies F p < 12 × 10 −8 J/cm 2 .
The non-monotonic behavior of the undoped sample S7 in Fig. 2(b) is reproduced by the exponential term exp(-F p /F p,abs ) that accounts for pump-induced increase of absorption [51].
In our experiment, we attribute pump-induced absorption mainly to the weak interaction of the FEL pump beam with shallow impurities in the silicon wafer. The interaction is weak because the pump photon energy is far from the impurity photoionization lines which clearly appear in the terahertz absorption spectra at 25-39 meV, however the wafer has a much larger interaction volume with the pump than the ACQWs, because the wafer is thicker by a factor 10 3 . The interaction of the FEL pump with silicon wafers with low density of impurities can be summarized as follows [53][54][55]: weak bleaching of transitions between bound impurity states for low and intermediate fluences, dominant pump-induced absorption for high fluence, and this behavior is actually observed in sample S7 (Fig. 2(b)). In the fitting of Eq. (1) to the data of sample S1 of Fig. 2(b), we have therefore fixed the value for F p,abs (within a 25% range of allowed variation) to the value obtained for S7 at the same pump photon energy ( ω p = 42.3 meV). We have then performed a fitting with F p,sat,exp as the only free parameter and we obtain a best-fit value F p,sat,exp = 7.2 × 10 −8 J/cm 2 for sample S1 at ω p = 42.3 meV. This value has then been kept fixed (within a 10% range of allowed variation) in the fitting procedure for all doping levels, temperatures and pump photon energies of Figs. 2(c)-(e). This procedure is motivated by the fact that F p,sat is weakly doping-and temperature-dependent as it is determined only by the ISBT absorption cross-section and by the non-radiative lifetimes [40]. While fitting all the data of Figs. 2(c)-(e), we have then left F p,abs as free fitting parameter. F p,abs decreases for each given sample and temperature when the photon energy approaches the Si-impurity transition at 39 meV (see Table 2). This confirms that the further is the pump photon energy from the Si impurity transition region, the less important is pump-induced absorption increase.
Experimental difficulties in the alignment of the sample precisely at the focal plane along the Gaussian beam waist resulted in the addition of another free parameter to the fitting of the data in Figs. 2(c)-(e), introducing a sample-dependent spot size. It has been demonstrated [44], also at THz frequencies [45], that the experimentally-determined saturation fluence is minimum when the sample is exactly in the focal plane, because the same input pump power P in can result in different fluences F p if the sample is positioned at an axial distance from the focal plane. This effect has to be taken into account in the calculation of the fluence inside the sample by multiplying P in by a scaling factor c p : F p = c p P in ∆t p πr 2 = P in ∆t p πr 2 where r' = r/ √ c p in our experiment ranges from 0.56 to 0.28 mm, and it is to be compared with the value of r = 0.50 mm measured with a THz camera without any sample. Reduction of the spot size in the presence of the sample (r' < r, or c p > 1, see Table 2) could be due to the high refractive index of silicon √ ε r,Si = 3.4. Indeed, the ideal value predicted by geometrical optics is r' = r/ √ ε r,Si = 0.50/3.4 = 0.15 mm. The maximum ideal value of c p is then ε r,Si = 11.6, so the values shown in Table 2 ranging from 0.8 to 3.1 indicate that, inside the sample, we are far from the ideal case of a diffraction-limited spot size due to poor waveguide input coupling. In summary, the fit of Eq. (1) to the data of Fig. 2 gives F p,sat,exp = (7.2 ± 3.0) × 10 −8 J/cm 2 . The large error on F p,sat,exp is dominated by the standard deviation of around 50% of the misalignment coefficient c p (see Table 2). The fact that c p obtained for different measurement runs is distributed symmetrically around its mean value confirms that c p is indeed a random experimental error variable. The value of F p,sat,exp is in disagreement with F p,sat,num = 17 × 10 −8 J/cm 2 , however the order of magnitude is the same and the two values are off by less than a factor two, therefore we believe that both the limitations of the saturable absorber model behind Eq. (1) [44] and the assumptions made in fixing the fitting parameters to specific values are presumably the cause of disagreement, rather than a failure of the numerical model.

Disclosures
The Authors declare no conflict of interest.