Terahertz radiation generation from metallic electronic structure manipulated by inhomogeneous DC-fields

We propose a feasible, high-efficiency scheme of primary terahertz (THz) radiation source through manipulating electronic structure (ES) of a metallic film by targeted-designed DC-fields configuration. The DC magnetic field is designed to be of a spatially inhomogeneous strength profile, and its direction is designed to be normal to the film, and the direction of the DC electric field is parallel to the film. Strict quantum theory and numerical results indicate that the ES under such a field configuration will change from a 3D Fermi sphere into a highly-degenerate structure whose density-of-state curve has pseudogap near Fermi surface. Wavefunctions’ shapes in this new ES are space-asymmetric, and the width of pseudogap near Fermi surface, as well as magnitudes of transition matrix element, can be handily controlled by adjusting parameter values of DC fields. Under available parameter values, the width of the pseudogap can be at milli-electron-volt level (corresponding to THz radiation frequency), and the magnitude of oscillating dipole can be at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-9} C*m$$\end{document}10-9C∗m-level. In room-temperature environment, phonon in metal can pump the ES to achieve population inversion.

The consensus that electronic structure (ES) of matter determines its optical property guides people to tailor/design suitable ES for achieving various optical applications. Compared with those secondary radiation sources  , primary radiation sources based on suitable ES are more efficient  . Here, the phrase "secondary radiation sources" refers to that radiation generation arises from nonlinear optics effect or nonlinear interaction between driving radiation beam (at other wavelength) and driven emitter materials. The phrase "primary radiation sources" refers to that the driving is not radiation or electromagnetic wave. Impurity doping, alloying, growing/synthesizing various heterostructures such as semiconductor quantum well and superlattice are popular methods of preparing novel ESs  . These methods can be categorized into chemical manipulation of ES, which refers to change matter's chemical constituent. Clearly, the chemical manipulation can lead to "permanent" change in the ES of matter.
In contrast, fields manipulation is to modify the ES having a "volatile" change, which refers to the ES recovering to its primitive form when removing the manipulating fields. For example, De Hass-Van Alphen effect is an example of DC-(magnetic)-fields-manipulation of metal's ES, and field-effect-transistor (FET) is an example of DC-(electric)-fields-manipulation of semiconductor carriers' ES 44,45 . Such a "volatile" change is more appropriate for many applications because of its simplicity, efficiency and economy. The fields-manipulation is more flexible in achieving tunable solid-state radiation source because it does not need to change materials for achieving another center-frequency output.
Obvious application prospect of radiations at terahertz (THz) band 46 promotes intensive investigations on achieving THz primary radiation sources. Scarcity of materials with ES characterized by THz-level gap leads to so-called "THz gap" 47 . Main limit of above-mentioned band-structure-engineering  is that the designed ES comprises too many levels which often form narrow bands or so-called mini-bands. In room-temperature environment, radiation generation from such an ES in which two mini-bands spaced by a meV-level gap will have a too broad spectrum undesirable for many applications. Moreover, the dipole matrix element in such an ES is of finite enhancement than those in atomic and molecular ES.
The advantage of the fields-manipulating metallic ES over band-structure-engineering of semiconductor carriers  is that larger charge density, as well as larger dipole magnitude, can be achieved in the former. Usually, the density of free electrons in a metal is at 10 22∼23 cm −3 -level, higher several order of magnitude than that of carriers in semiconductor (about 10 10 cm −3 ). Therefore, manipulating free electrons in a metal is more hopeful to achieve higher output power than manipulating carriers in a semiconductor.
In this work, for more practical purpose, we investigate manipulating the ES of a metallic film by a special DC-fields configuration in which the magnetic field is space-inhomogeneous. It is well-known that a homogeneous DC magnetic field B 0 can lead to an ES in which the general form of 3D wavefunctions is a product of a 2D plane-wave (in terms of (z, x) ) and a 1D harmonic-oscillator wavefunction (in terms of y). As shown latter, orthogonality among 1D harmonic-oscillator wavefunctions enables the dipole matrix elements for the transition between two Landau levels in such an ES to be confined at 0 and hence prohibits dipole radiation generation from the transition between two Landau levels. In contrast, targeted using inhomogeneous B 0 e z ∼ x or y can lead to an ES whose 3D wavefunctions have 2D localized parts losing orthogonality. In such an ES, the above-mentioned prohibition on dipole radiation generation can be removed, and hence a practical radiation source is possible.

Theory and method
Electronic structure manipulated by homogeneous DC-fields, revisiting De Hass-Van Alphen effect. According to many textbooks 44,45 , free electrons under homogeneous DC magnetic field can be described by a class of 3D eigenstate wavefunctions whose general form is a product of a 2D plane-wave exp(ik z z + ik x x) and a 1D harmonic-oscillator wavefunction ψ 1D n y . Here, if 3D eigenstate wavefunctions of a quantum system are of a general form which is a product of a 1D extended-state wavefunction and a 2D localized-state one, i.e. ψ 3D k z ,n (x, y, z) = exp(ik z z)ψ 2D n (x, y) , it is not certain for these 3D wavefunctions being orthogonal mutually. For example, if the 2D wavefunctions ψ 2D n (x, y) are governed by a 2D Hamiltonian operator as follows: where α , β x and β y are constants and the "potential-energy" function V does not contain differential operators such as ∂ x , it is easy to prove that orthogonality exists among them. It is well-known that orthogonality exists among harmonic-oscillator wavefunctions, and same conclusion is also hold for plane-waves.
In such a ES, dipole matrix elements ψ 3D k z ,n * x, y, z ψ 3D k z +q,m dxdydz to be confined at 0 if q = 0 and m = n . For generating a photon with momentum q and energy qc , momentum conservation and energy conservation demand two levels having a difference in their wavevectors. Thus, ψ 3D k z ,n * x, y ψ 3D k z +q,m dxdydz are = 0 because of orthogonality among plane-waves, and ψ 3D k z ,n * z ψ 3D k z +q,m dxdydz is = 0 because of orthogonality among harmonic-oscillator wavefunctions. In short, if 2D, as well as 1D, localized-state wavefunctions have orthogonality, these dipole matrix elements are 0, which prohibits dipole radiation generation.
In contrast, 2D localized-state wavefunctions can loose orthogonality for 2D Hamiltonian operator in suitable forms. For example, if the 2D Hamiltonian operator reads where f x and f y are functions of variables x, y rather than two constants, its eigenstate 2D wavefunctions can be not orthogonal mutually. Once 2D localized-state wavefunctions loose orthogonality, these dipole matrix elements, such as ψ 3D k z ,n * z ψ 3D k z +q,m dxdydz , can be ∼ q −2 � = 0 for cases with q = 0 and m = n in which ψ 3D k z ,n * x, y ψ 3D k z +q,m dxdydz = 0 still exist. This removes a prohibition on dipole radiation generation.
Electronic structure manipulated by inhomogeneous DC-fields configuration. Above discussion enlightens us to design/tailor the Hamiltonian governing a quantum system. Hamiltonian in following form has a 2D-part meeting the general form Eq. (2) where inhomogeneous DC magnetic field and its vector potential read B 0 , β and E 0 are constants. Although Hall current, which is common for metal under E × B configuration, can exist (because B is along z-axis and E is along y-axis in this configuration and hence a Hall current is along x-axis), this Hall current represents a secondary effect and is too weak compared with controlling current in solenoid responsible for producing B. Usually, perturbation Hamiltonian associated with such a Hall current (which is a consequence of another stronger perturbation Hamiltonian), is taken as too weak to being added into total Hamiltonian. Corresponding experimental setup is illustrated as Fig. 1. A pair of Helmholtz square-cross-section coils is on purpose mis-aligned to be not co-axial and hence produces a magnetic "slope", which refers to B 2 z is the function of y. The metallic film is separated from two electrode plates by two gaps, which might be filled with highpermittivity insulator or not. According to detailed description elsewhere 48 , it is feasible to attain a dimensionless r.h.s = α * ∂ xx lnφ + (∂ x lnφ) 2 + ∂ y ∂ y lnφ + iwxy + ∂ y lnφ + iwxy 2 − sy + α * k 2 z ; x; (12) y = s + 1 s 2 + 1 su + s − 1 s 2 + 1 sv; Figure 1. Sketch of field configuration and experimental setup. Note that two gaps might be filled with highpermittivity insulator, which can also act as an "holder" holding the film. If the "holder" is undertaken by thin dielectric substrate, suitably adjusting charge current in coils/solenoids can still warrant the metallic film to feel a magnetic slope. Note that cross-sections of solenoids are of square shape, which is not reflected in the scheme (due to the limit of plot software) and hence have to be stated clearly here. where coefficients f i,j and g i,j are complex-valued constants. These expansion coefficients can be solved from above equations through following strict procedure: for Eqs. (19,20), order-by-order comparing their terms proportional to constant u 1 and v 1 , we can obtain 6 algebra equations. For Eq. (21), 2 algebra equations can be obtained by order-by-order comparing its terms proportional to u 1 and v 1 . Moreover, ∂ x Eq.(21) = 0 and ∂ y Eq.(21) = 0 also lead to 2 algebra equations. Note that 2 formulas g 0,0 = i k y + k x and f 0,0 = i k y − k x can be derived from definitions Eq. (19). Thus, 12 algebra equations and formulas (see "Appendix") can determine solutions of these these low-order expansion coefficients g 1,0 , g 0,1 , f 1,0 , f 0,1 , f 1,1 , g 1,1 . In this strict procedure, g 0,2 , g 2,0 , f 2,0 , f 0,2 act as intermediate variables.
Finally, E k − E 0 k is a function of g 1,0 and f 0,1 and hence can be known after these low-order expansion coefficients are solved according to this procedure (whose detailed formulas are presented at "Appendix" section).
According to formulas in "Appendix" Eq. (83), E k is still ∼ k 2 z but has a more complicated dependence on This will lead to corresponding variation in density-of-state (DOS) curves. Moreover, all wavefunctions shape |ψ k | are space-asymmetric and hence correspond to dipole moment Density-of-state of the ES. When external fields are absent, conduction-band of a metal corresponds to a Fermi sphere ES in 3D k-space and has a top at Fermi surface |k| = k F and a bottom at |k| = 0 . Usually for light metals, their k F are at 10 4 µm −1 -level 44,45 . In the De Hass-Van Alphen effect, the ES becomes a set of 2D cylindric surfaces in 3D k-space and each cylindric surface corresponds to a sub-band 44,45 . Above result implies that, under inhomogeneous DC-fields configuration, the ES changes from a 3D Fermi sphere in which E 0 k ∼ k 2 x + k 2 y + k 2 z and 0 ≤ |k| < k F exist, to a corrugated 3D Fermi sphere. As shown in Fig. 2, each 2D cross-section (on the k x − k y plane) of the 3D Fermi sphere can be divided into 6 segments, 3 corresponding to positive energy shift and 3 for negative. For θ = π 4 , 3π 4 , 5π 4 , 7π 4 , the energy shift E can be ±∞ and correspond to nearly zero DOS, or dθ d�E ∼ 0 if θ = π 4 , 3π 4 , 5π 4 , 7π 4 . Thus, for minimizing total energy, some regions on the k x − k y plane, which are "unoccupied" or outside 2D Fermi circle, will be occupied because of negative energy shift over them, and some regions, which are "occupied" or inside 2D Fermi circle, will be unoccupied because of positive energy shift over them. For each |k � |-value, www.nature.com/scientificreports/ the dependence of the energy shift on θ , according to Eq. (83), will lead to a DOS curve with pseudogap (see Fig. 2c,b) because the DOS can arrive at its maximum at two θ-values: θ = 0 and θ = π . If each |k � |-value has such a DOS curve with pseudogap, their superpositions, as shown in Fig. 2d, will be possible to remain a pseudogap. The pseudogap can be tuned, under available values of controllable parameters, to be at meV-level. According to Eq. (83), the fact s 2 < 0 (see Eq. 9) enables the width of the pseudogap at each k -value, or 8ws 4 [s 2 +1] 2 * 1 |k � | , to be large enough. Moreover, the maximum of allowed k -value corresponds to the narrowest among pseudogaps. Thus, as shown in Figs Space-asymmetric wavefunction shape |φ k x ,k y | implies each k-state corresponding to a dipole moment d x,y ≡ |ψ k | 2 x, ydxdy . According to "Appendix", the lowest order expansion terms of Relnφ k x ,k y is ∼ xy . Thus, values of d x,y are dependent of a parameter h 1,1 ≡ ∂ x ∂ y lnφ | x=0,y=0 . Larger |h 1,1 |-value will correspond to larger dipole moment. As reflected by Fig. 5, larger k -values, which correspond to states near Fermi surface, correspond to larger |h 1,1 |-values.
According to familiar formula for radiation generation from quantum transition 44,45 , the radiation strength I, or transition matrix element, is proportional to both the product of DOS of initial state and that of final state and the product of dipole moment of initial state and that of final state: Note that high dipole moment and high DOS can be simultaneously achieved at some θ-regimes. For example, Fig. 5a indicates that h 1,1 can be large enough at θ = 0.5π , and Fig. 2c indicates that θ = 0.5π corresponds to a very flat segment of the θ − �E curve and hence to a large DOS-value. Therefore, the corrugated ES described above, which has a tunable width of pseudogap and sufficiently high dipole moments at states near peaks around pseudogap in DOS curves, is very favorable to powerful radiation output.
The partition function of the ES, which is taken as a canonical ensemble, reads where K B is Boltzmann's constant. The ground state of such an ES is to fill the complicated 3D body and macroscopically corresponds to a static dipole (because all k z , k , θ -states are of "polarized" shapes).

Parameters in laser rate equations.
The new ES has a main feature: wavefunction |ψ k x ,k y ,k z >= |ψ k x ,k y x, y |exp ik x x + ik y y exp(ik z z)exp(iω k t) have a space-varying shape |ψ k x ,k y x, y | . The where ω 0 ≡ E 2 − E 1 > 0 , and other parameters are of respective customary definitions 48 . The value of A 12 depends on not only energies of two levels, E 2 and E 1 , but also space shapes of their wavefunctions, |ψ 2 | and |ψ 1 |. It is well-known that the dipole matrix element of an ES is proportional to the characteristic scale of wavefunctions of the ES, | < ψ 1 |r|ψ 2 > | ≡ ξ −3 |ψ 1 (r/ξ )|r|ψ 2 (r/ξ )|d 3 r ∼ ξ . For small quantum systems such as atoms and molecules, their wavefunctions have small space extension or highly localized. For atomic wavefunctions, these exist |ψ atom | ∼ exp(−r/ξ ) and ξ ∼ 0.053 nm or Bohr radius (see many textbooks). As a consequence, the value of dipole matrix element of an atom | < ψ 1 |x, y, z|ψ 2 > | ∼ ξ will be about at nm-level and hence corresponding t spon,12 -values, for |E 2 − E 1 | at eV-level, is ∼ 10 −9 s (see Eq. 8.3.12 at page 168 of Ref. 50 ). Similar results hold for small molecules. Moreover, in semiconductor quantum well and superlattice, wavefunctions of carriers are of larger ξ-value, usually tens nm, than atomic and molecular wavefunctions. Impurity levels also have larger ξ-value because of larger effective electronic mass in Hydrogen-like model on impurity. In contrast, extended-state wavefunctions in the ES of a metal is of larger ξ-values, as well as larger values of dipole matrix element, than those in quantum well and superlattice.
The wavefunctions' shape also affects electron-phonon coupling constant. Electron-phonon interaction can be represented by a perturbation Hamiltonian density 45 where d n = A q cos q · R n − ωt is the displacement of the n-th atom and V is atomic potential. Corresponding scattering probability, or electron-phonon coupling constant g 1,2 , reads 45  www.nature.com/scientificreports/  www.nature.com/scientificreports/ where phonon number A 2 q can be obtained from Bose-Einstein distribution function: g BE ω q ; T ≡ 1 exp( ω q /K B T)−1 . Clearly, M 1,2 depends on phonon number A 2 q and space shapes of |ψ 1 > and |ψ 2 > . Although both the electron-phonon coupling constant and the spontaneous emission lifetime are affected by the space shapes of wavefunctions, their dependencies on the shapes are different. This implies that with the ES changing (i.e., shapes of wavefunctions changing), the ratio among these coefficients in rate equations of radiation generation (see below) also change and hence so do its solutions 51,52 . Moreover, there is difference between the scattering by longitudinal-wave acoustic phonon u lo = exp(iq x x + iω lo q t)e x and that by transverse-wave one u tr = exp(iq x x + iω tr q t)e y . In the former case, u · ∇V ∼ x * x while in the latter case, u · ∇V ∼ x * y (because ionic potential V is a function of x 2 + y 2 + z 2 ). Thus, the difference between two electron-phonon coupling constants is dependent on shapes of wavefunctions Likewise, the spontaneous emission coefficient A 12 also depends on shapes of wavefunctions in a different manner The rate equations of radiation generation based on 2-level read where P 2 and P 1 are populations of two levels. Here, because conventional components for achieving amplified radiation, such as reflecting mirrors or resonant cavity 52 , are not used, there is no feedback usage of generated radiation and hence the population inversion is merely due to the phonon pump, i.e., not amplified. The output radiation is indeed an enhanced spontaneous emission from an inverted equilibrium population. Because no radiation is input, such an enhanced spontaneous emission can also be viewed as an infinite-fold amplification, i.e., the enhancement in spontaneous emission is viewed as infinitely amplifying 0 input radiations. Therefore, we borrow the phrase "stimulated emission", despite not very accurate, to refers to this enhancement in spontaneous emission.
Equations (32,33) can be derived in following way. If phonon scattering is not taken into account, populations of two levels are where N(E 2 ) and N(E 1 ) are DOS of two levels, respectively. During t the number of electrons transiting from E 2 to E 1 , due to phonon-scattering, is 45 where the factor 1 − g FD (E 1 ; T) represents un-occupation probability of the E 1 -state and hence means that the transition is prohibited if the E 1 -state has been occupied. Likewise, during t the number of electrons transiting from E 1 to E 2 is 45 Therefore, after taking into account the phonon scattering, populations of two levels read where Equations (38,39) can be re-written as Eqs. (32,33) if introducing (29) g tr 2,1 ∼ |ψ 1 |iqxy|ψ 2 |dxdy; (30) g lo 2,1 ∼ |ψ 1 |iqx 2 |ψ 2 |dxdy; (31) A 12 ∼ |ψ 1 | x 2 + y 2 |ψ 2 |dxdy.

Radiation generation from 3-level model and 2-level model. Many authors have taken into account
phonon scattering in radiation generation from semiconductor 36,43,51 . For metal film held by dielectric "holder", environment temperature can be felt by the film due to the connection between the "holder" and vacuum chamber wall. Because of the difference between g tr and g lo , radiation generation is easy to involve in three levels: one from a state below the pseudogap and two above. An electron initially at a state below-pseudogap is scattered by a high-energy longitudinal-wave phonon to a state above-pseudogap and then is scattered by a low-energy transverse-wave phonon to another state above-pseudogap, and then return to the state below-pseudogap by emitting a photon. Because both phonon and photon have linear dispersion relation and their phase velocities satisfy c tr ≡ ω tr q < c lo ≡ ω lo q < c (because of different values of elastic constants for longitudinal-wave and transversewave 53 ). As illustrated in Fig. 6, this inequality relation ω tr q < ω lo q < c can warrant a closed energy-momentum conservation process among phonon and photon. Moreover, two electronic states within a same parabolic band cannot ensure the difference between their momentum and that between their energy to satisfy a linear dispersion relation hold by photon, but two states from different parabolic bands can. These features well illustrate why the beating between longitudinal-wave phonon and transverse-wave one can allow an inter-band transition. Why this beating cannot cause a 3D Fermi sphere ES to generate photon can also be explained (because any two electronic states are in a same parabolic band E k ∼ k 2 ). Moreover, as previously mentioned, because different parabolic bands correspond to different 2D shaped-parts of 3D wavefunctions, losing orthogonality among 2D shaped-parts also favor the dipole matrix elements for the transition k z + q, m → (k z , n) to be = 0 . Namely, there are two prohibitions on the dipole radiation generation, one is the orthogonality of 3D wavefunctions and the other is parabolic dispersion relation of the 1D plane-wave part of 3D wavefunctions.
Because sound speed in many metals can be 10 3∼4 m/s-level, if c lo = 5 × 10 3 m/s and c tr = 3 × 10 3 m/s, there will be k 3 q 2 = c lo −c tr c−c tr ∼ 2 × 10 −5 and hence yielding a photon k 3 c ∼ meV will need a phonon q 2 c lo ∼ c−c tr c lo −c tr c lo c meV ∼ 2.5 meV , which can be supplied sufficiently for most metals in room-temperature environment 300 K.
Rate equations of radiation generation based on 3-level model read As shown in Eq. (28), the dependence of M 1,2 on phonon strength A 2 q ∼ g BE ω q ; T determines A 2 (because g BE ω q 13 ; T < g BE ω q 23 ; T if q 13 > q 23 ), and hence M 1,3 < M 2,3 or σ 13 < σ 23 .
As previously mentioned, eV-level gap between atomic levels can lead to A ∼ 10 9 s −150 . Thus, even if assuming meV-level gap exists in atomic ES, corresponding A will be ∼ 10 0 because of A ∼ ω 3 0 as shown in Eq. (26). Luckily, in the ES above-described, the ξ-value is ∼ 2π/ k 2 x + k 2 y ∼ µm , about 10 5 -folds of that of atomic one. Thus, the A-value in the new ES can reach ∼ 10 10 s −1 . Moreover, for many metals with a Fermi sphere ES, phonon scattering can lead to a relaxation-time about 10 −13 ∼ 10 −14 s or σ ∼ 10 13 ∼ 10 14 s −145 . The new ES enables these σ-values to be kept at this level.
Actually, above-mentioned two characteristic time scales: 1/A and 1/σ , represent respectively the strength of electron-photon interaction/coupleing and that of electron-phonon one, which are expectation values of two perturbation energies on electrons ψ * k+q H e−photon ψ k d 3 r and ψ * k+q H e−phonon ψ k d 3 r , where H e−photon = eE · r and H e−phonon = u · ∇V . Here, under room temperature environment, ionic displacement u can be estimated, according to usual values of elastic constant K of metals ∼ 10 2 GPa = 10 11 J/m 344 for yielding elastic energy 1 2 Ka 2 * u ∼ 1 meV = 1.6 × 10 −22 J , where a ∼ 100pm represents the size of a cell, to be ∼ 10 −2 a ∼ 1pm . Such a displacement, 10 −2 folds of inter-atom distance a, can correspond to a 10 meV-level energy perturbation on an electron because the inner electric field V/a in most metals can be at V/[100pm]-level. Namely, for most metals in 300K environment, ψ * k+q H e−photon ψ k d 3 r < ψ * k+q H e−phonon ψ k d 3 r , especially ψ * k+q H e−THz ψ k d 3 r << ψ * k+q H e−phonon ψ k d 3 r , can be warranted. Note that the energy increment of an excited electron is not the kinetic energy of a phonon, instead, it is the variation of Coulomb energy associated with the ionic displacement, u · ∇V c , that contribute this energy increment. Although the phonons number is not too large, electronic excitation is still very obvious. Therefore, it is not necessary to seek for applying large-amplitude THz accoustic wave on the metal or inputting large mount of phonons. This greatly lessen experiment condition requirement. Otherwise, we have to seek for a powerful THz mechanics vibration source whose difficulty is illustrated as follows: to input THz accoustic wave into the metal, we utilize the well-known formula reflecting cantilever transverse oscillation frequency ν inverse proportional to its extruded length L (see page 120-121 of ref. 53 ), i.e. ν = π 2 L 2 EI ρS , for a typical metal material with following representative values: elastic module E = 2 × 10 11 Pa = 2 × 10 11 kg × 1 m × 1 s 2 , cross-section area S = 1 × 10 −6 m 2 , density ρ = 2700 kg/m 3 , moment of inertia I = ab 3 /12 , I/S ∼ m 2 , extruding length L = 10 −3 m,there will be ν ∼ MHz. Namely, it is difficult to make an object whose material is common and size is in a macroscopic scale ( > µm ) to do a vibration with a frequency at THz. A man-made macroscopically-sized THz mechanics vibration source is nearly not available.
But the metal's "native" THz phonons are of sufficient mount at 300 K. Because the number of phonons at a specified frequency ν is the product of Bose-Einstein distribution function g BE at ν and phonon DOS at ν : ≫ 0 if T = 300 K and = 400 K or higher. Fig. 6c indicates that this ratio is close enough to 1 at some typical situations. Moreover, the frequency  Figure 6. (a) A sketch illustrates that three waves with linear dispersion relation can form a closed energymomentum conservation process. The slope of each straight-line segment represents its phase velocity. Because of c tr < c lo < c , q 2 = q 1 + k 3 and c lo q 2 = ω 2 = ω 1 + E k = c tr q 1 + ck 3 can be satisfied, where c lo and c tr are phonon's phase velocities and c is light velocity. (b) A sketch illustrates that a transition between two electronic states k 4 , k 5 within a same parabolic band cannot ensure the momentum difference k 5 − k 4 = k 3 and the energy difference E 5 − E 4 = k 2 5 − k 2 4 to satisfy a linear dispersion relation E 5 −E 4 k 5 −k 4 = const , which can be ensured by a transition between two states from different parabolic bands or E 5 = E 4 + 2k 3 * k 4 + k 2 3 + E g . (c) The dependence of the ratio between the number of phonons whose frequencies > 1 THz and total number of phonons on Debye temperature. www.nature.com/scientificreports/ corresponded by maximum phonon DOS in many metals can be up to tens meV 44,45 , therefore a metal's "native" THz phonons are in sufficient mount for pumping electrons. By choosing suitable values of controllable parameters, rate equations for 3-level model can yield solutions satisfying population inversion. This can be referred to Fig. 7. Actually, even for 2-level model, suitable values of controllable parameters can also warrant population inversion even for 2-level model (see Fig. 7).
T Hz is defined likewise. Usually, free electron density in a simple metal is at 10 22∼23 cm −3 -level, which corresponds to an inter-atom interval at nm-level. Because THz frequency is non-transparent relative to metal medium with such a 10 22∼23 cm −3 -level density , it is only a thin layer near the surface that can have contribution to THz output. For a film 1 cm × 1 cm × 0.01 cm, the thickness 0.01 cm is at the level of THz = 0.03 cm, the wavelength of THz photon, and hence it is suitable to view the layer of a 0.01 cm-thickness as the "skin" layer of a THz EM wave.
Thus, a metallic film 1 cm × 1 cm × 10 −2 cm can contain 10 20∼21 conduction-band electrons. For the new ES, states whose energies are within the range [E F − 0.05 * K B T, E F + 0.05 * K B T] have obvious contribution to various quantum transition processes, and the number of these states is estimated to be proportional to 0.1 * K B T/E F , which is ∼ 2.5 meV/10 eV ∼ 10 −4 for metals in room-temperature environment 300 K, and hence is ∼ 10 16∼17 . Of course, if the width of the range is 0.01 * K B T , the number of these states will be ∼ 10 15∼16 , still higher several orders of magnitude than that of carriers in semiconductor. According to above formula, it is feasible to an output power at Watt-level from a film 1 cm × 1 cm × 10 −2 cm which can contribute 10 15 electrons to participate oscillation with a THz frequency and hence yield 10 15 * meV = 1.6 × 10 −7 J energy over a time-scale at ps-level.
Moreover, it is easy to understand that intraband carrier acceleration by radiated THz electric field can lead to significant THz absorption and hence decrease the THz output power. Luckily, the THz-active layer is only of a thickness below the THz wavelength, ∼ 0.01 cm < THz . Before E x,y is significantly depleted by intraband acceleration, it has transmitted the thin layer. Moreover, because space-asymmetric wavefunctions' shape, intraband carrier acceleration to higher k x -state might imply larger THz dipole moment and hence does not definitely imply "significant THz absorption". A definite answer will be the task of future works.
Beside, the coupling loss to free-space THz radiation is also a factor affecting THz output power. For outstanding main purpose, this work avoids to consider various complicated secondary effects/mechanisms affecting THz output power. For example, intraband carrier acceleration above-mentioned can "heat" carriers and  www.nature.com/scientificreports/ then "heat" ions, through carrier-ion collisions, and cause phonons energy rising and then in turn cause more THz EM generation. The main purpose is on the feasibility of fields-manipulating ES for achieving a specified wavelength EM output. These secondary effects/mechanisms will be content of future works.

Summary
Fields-manipulating the ES of a metal enables high-density free electrons to form a new "volatile" ES, which can couple with acoustic phonon to achieve population inversion between states spaced by meV-level pseudogap. The fact that mobile electrons in metal is higher several orders of magnitude than that in semiconductor determines the technical scheme presented here to be of marked application prospect.