Terahertz radiation generation process in the medium based on the array of the elongated nanoparticles

We conduct a theoretical study and numerical simulations of terahertz radiation generation in the medium based on armchair-edge nanoribbons and zigzag nanotubes with metallic conductivity. The multicascade mechanism of radiation generation is considered in the task of terahertz radiation generation. The level of the injection current in nanoparticle arrays has been estimated. The task expands to the similar medium where radiation current is generated with the use of infrared radiation stimulated absorption, for example, radiation of a CO2 laser. For the effective dielectric function three basic models are employed: effective media approximation (Bruggeman’s theory) [nanotube bundles], Maxwell-Gatnett theory [ordered array], and Maxwell-Garnett theory with the Clausius–Mossotti geometric correction [ordered array of nanotubes].


Introduction
A quantum cascade laser (QCL) is a device, which combines aspects from different fields such as nanoelectronics and quantum engineering, plasmonics, as well as subwavelength and nonlinear photonics. Since their invention, QCLs ) have been rapidly developed into an important class of lasers, covering a wide wavelength range 3-300 μm. Such a light amplification is possible in intersubband transitions, i.e. transitions between quantised energy states within one energy band of a semiconductor, which was first predicted by Kazarinov and Suris (1971) in a seminal paper. We agree that the fast nonradiative longitudinal optical (LO) phonon emission rate can be harnessed for the laser process [the first QCL intersubbaned laser, designed to emit at 4.3 μm wavelength, was demonstrated in Faist et al. (1994) and Capasso et al. 1999)], rather than being an hindrance for intersubbaned laser action. Zheng et al. [5,Zhang] considered, that infrared lasers had been developed on the basis of and intersubband transition in semiconductor heterostructures with quantum wells and super lattices of type-II. After appearance of works (Faist et al. 1994;Capasso et al. 1999), based on the work (Kazarinov and Suris 1971), the significant progress in the pulse lasers development is reached.
For instance, a Continuous Wave laser, functioning at cryogenic temperatures was invented (Faist et al. 1995a(Faist et al. , 1995b; a laser, functioning in a pulsed mode up to room temperature was also invented by Faist et al. (Faist et al. 1996). Considerable progress in the power output (Bai et al. 2011) and the available frequency range (Cathabard et al. 2010;Scalari et al. 2010) was made. The first far-infrared QS lasers with wavelengths larger than 20 μm (21.5 and 24 μm) were shown by Colombelli et al. (Colombelli et al. 2001). Generation of ultra-short pulses in mid-infrared (Wojcik et al. 2013) and terahertz QCLs (Barbieri et al. 2011;Freeman et al. 2012) is a promising direction in Physics.
Development of innovative QCL types and subsequent design optimization, systematic improvements of the quantum cascade lasers operating temperature, efficiency and a spectral range require detailed simulation of the underlying physical processes in these structures (Jirauschek et al. 2014). For instance, THz generation at the room temperature was obtained at the power levels up to 120 μW (Lu et al. 2013;Vijayraghavan et al. 2013). It is necessary for most technical applications to push the available room temperature output power to a few mW. Simulation of such terahertz QCL sources is especially demanding, because it requires the coupled simulation of two mid-infrared QCLs and careful simulations of the nonlinear susceptibility and the associated frequency conversion process in the heterostructure (Jirauschek et al. 2013). Thus, a quantitative simulation, reliable numerical design optimization and explorations must also include the optical cavity (Vijayraghavan et al. 2013(Vijayraghavan et al. , 2012Fathololoumi et al. 2012;Belkin et al. 2008). A further example is surface emission schemes based on one-and two-dimensional photonic crystal structures (Demichel et al. 2006;Chassagneux et al. 2009;Colombelli et al. 2003) offering tailorable emission properties and improved beam quality. The simulation of such subwavelengthstructured cavities requires advanced electromagnetic simulation, e.g., based on the coupled mode theory (Schubert and Rana 2006) or even full finite-difference time-domain simulations of Maxwell's equations (Colombelli et al. 2003;Taflove and Hagness 2000).
We can see the unique situation in the case of small-sized nanostructures. Evolution equations, written in such structures, represent the difference approximation of differential equations. Studying the processes of Low-dimensional nanostructures permits the application of well-established numerical methods and properties of different schemes. For example, the evolution of charge carriers (electrons and holes) in the terms of a tight binding model for two-atomic unit cells is described by the Dirac equation as particles and antiparticles with zero-mass (Katsnelson 2010;Morozov et al. 2008;Saito et al. 1992;Vozmediano et al. 2010). One of the prominent results is the usage of nanotubes of the continuum model (kp-type model) for description of charge carries in graphene (Onipko and Malysheva 2017;Malysheva and Onipko 2009). In Sadykov and Jolnirov (2021), based on the kp-type model, an approach for the approximate calculation of the characteristics of electromagnetic guided waves on almost circular, closely packed bundles of parallel, identical, and metallic carbon nanotubes is proposed. A kp-type model yields the polarizability scalar and the antenna efficiency of a finite-length SWCNT bundle in the long-wavelength regime over a wide frequency range spanning the terahertz and the near-infrared regimes. The approach used in Sadykov and Jolnirov (2021) is similar to the approach from Yan et al. (2016), in which the Fraunhofer's multiple slits diffraction model is used to explain the far-field pattern of a phase-locked quantum cascade laser array.
In Kibis et al. (2005) and Saito et al. (2010) a two-point unit cell is considered for zigzag SWCNTs. In White et al. (2007) and Saroka et al. (2018Saroka et al. ( , 2017 every unit cell consists of two zigzag chains, bounded by hydrogen atoms. Here as in Onipko and Malysheva (2017) the unit cell of 2 N atoms, that consists of two zigzag chains ; Fig. 4], where sublattices A, B, C and D have carbon atoms, is considered for N staggered armchair-edge GNRs in contrast to Kibis et al. (2005), Saito et al. (2010), Saroka et al. (2018) and Saroka et al. (2017).
In Portnoi et al. (2015) and Saroka et al. (2019) it is shown, that the matrix element of the velocity operator for interband dipole transitions in the vicinity of the Dirac point of an armchair graphene nanoribbon and a zigzag carbon nanotube has a universal value, which is equal to the Fermi velocity of an electron in graphene. This result is discussed in Shahnazaryan et al. (2019) and Hartmann et al. (2019).
In Saroka et al. (2019) it is also shown, that optical matrix elements in zigzag nanotubes for incident radiation, polarized parallel to the translational symmetry axis, are compared with the corresponding matrix elements in armchair graphene nanoribbons.
In Saroka et al. (2019) and Son et al. (2006) it is shown, that the curvature effects for tubes and the edge effects for ribbons result not only in a small bandgap opening, corresponding to terahertz (THz) frequencies, but also in a significant enhancement of the transition probability rate across the bandgap.
In  the author's approach combines the tight-binding model with the continuum model of kp-type to investigate finite length graphene nanoribbon with armchair edges subjected to the external electric field longitudinal to the edges and derive approximate solutions for the electron energies and wave functions around the Dirac point. In the work on the simulation of a cascade laser similar to Katsnelson (2010) and Sadykov et al. (2016) it is shown that generation of terahertz radiation is possible at reasonable experimental parameters. It is also discussed how generation of radiation by Cherenkov mechanism is possible (Sadykov and Aporoski 2019). In addition it is claimed that equations obtained for electrons in armchair graphene nanoribbons can be reduced to the Majorana equation (Sadykov and Aporoski 2017).
The paper focuses on a theoretical study and numerical simulation of the task of amplification of terahertz radiation by analogy with QC lasers in medium, based on non-interacting armchair edge nanoribbons with metallic conductivity, and on the estimation of the value of generating radiation and the value of radiative current density. The case of noninteracting metallic nanotubes for non-densely and randomly distributed composites of particles is discussed in Sadykov et al. (2016). The effective dielectric function in the geometrical model arrays is defined by the individual particle polarizability with consideration to the geometry of the array (geometric factor). It is shown that both the medium based on randomly distributed composites of particles (the effective medium-approximation, the Maxwell-Garnett's Theory [MGT]) and the array of particles (the arrays geometrical model) from parallel oriented CNTs and armchair-edge nanoribbons.
The structure of the paper is the following. Section 1 is an introduction that presents the literature review. Section 2 presents the interaction of a classical optical field with an equidistant multilevel system using the density matrix formalism. The wave functions, the eigenvalues of the energy states of charge carriers, and the matrix elements of the dipole moments in armchair nanoribbons and zigzag nanotubes of finite length are written out. In Sect. 3, in the medium based on the array of noninteracting zigzag nanotubes or noninteracting armchair-edge nanoribbons, the process of generating terahertz radiation is considered on the basis of the obtained differential-difference equation for the diagonal elements of the density matrix. Section 4 is devoted to a theoretical study of the permittivity function for a crystal lattice of scattering nanoparticles with a metallic type of conduction (CNTs and GNLs) with nanoparticle centers located at the same distance on the basis of the mean field theory and the Maxwell-Garnett approach. Section 5 provides the results of numerical simulation of the task of terahertz radiation generation in the medium based on the parallel oriented graphene armchair edge nanoribbons and zigzag nanotubes. Section 5 also presents the calculation results obtained for the effective dielectric functions in medium based on the array of non-interacting zigzag nanotubes or non-interacting nanoribbons with the edges of the chair. Section 6 is devoted to the discussion of obtained results. Conclusion is the last section.

Nanotubes
First we consider the array of chiral nanotubes as the nonlinear media and then generalise the obtained results for the array of armchair-edge nanoribbons. The theoretical part of the task for the zigzag nanotubes is viewed in Sadykov et al. (2016), in which the analytically obtained results in Sadykov and Aporoski (2017) are used. In Sadykov and Aporoski (2017) on account of an elementary phase cell, which contains sublattice A and sublattice B, the wave functions, equidistant energy levels and dipole matrix elements of laser transition in the presence of a longitudinal electric field are analytically obtained. In Sadykov and Skorkin (2013) as distinct from Sadykov and Aporoski (2017) the array of chiral nanotubes is viewed. We write the Liouville equation for the considered multilevel system for the density matrix where Ĥ ,̂ is the commutator of Hamilton operator Ĥ . We use Eq. (2.1) (where T 1 is the longitudinal time relaxation; T 2 is the transverse relaxation time) and then write the interaction of a classical optical field with a multilevel system with the use of density matrix formalism as (Sadykov et al. , 2016 where Ĥ = −d z E z is the Hamilton operator; E z = E∕2 + E * ∕2 , E and E * are the complex conjugate quantities; E ∼ exp(−i t) , (d z ) fn = e(z) fn = −2L∕ (f − n) to final state | | | k ⟩ associated with stimulated absorption, the current will have the dependence Such a conformity takes into account the changes in quantites k k ∕ t and k k ∕ t , which are caused because of the transition rate between these states. The first summand means the electrons' injection to finale state | | | k ⟩ , the second one means the reduction of the electrons' transition rate at initial state � � k ⟩ . In the case of the multicascade mechanism of radiation generation current will have a dependence

Narrow nanoribbons
We generalise the results obtained in 2a for the array of the narrow armchair-edge nanoribbons ). In  basing on the four-point unit cell model (four zigzag chains A, B, C and D; see Fig. 4) for the armchair-edge nanoribbons, we consider the staggered configuration. This model almost coincides with the model used in Onipko and Malysheva (2017) and Boyd (2003). In Sadykov and Skorkin (2013) as distinct from Sadykov and Aporoski (2017) the array of chiral nanotubes is viewed. In  A, B chains correspond to gA and gB chains from Sasaki et al. (2011) when g is odd; if g is even, C and D from our paper chains correspond to gA and gB. In  the obtained wave functions, equidistant energy levels and dipole matrix element values in the presence of a longitudinal electric field for the armchair-edge nanoribbons were absolutely the same as analogous functions from Sadykov et al. (2016) and Sadykov and Aporoski (2017) for the zigzag nanotubes. Thus, later on we consider Eqs. (2.5) and (2.6) to be valid also for the armchair-edge nanoribbons with the staggered configuration.

Approximation of differential-difference equation for the diagonal elements of the density matrix with a parabolic equation
Scientists focus on ZGNRs and armchair SWCNTs (Saroka et al. 2017(Saroka et al. , 2018Wakabayashi et al. 2010). For example, a hidden correlation between the optical absorption resonances of several zigzagr-edge nanoribbons and armchair nanotubes has been found in the nearestneighbor tight-binding model (Saroka et al. 2017(Saroka et al. , 2018. In Wakabayashi et al. (2010) we study the interrelations between optical absorption resonances in ZGNRs and armchair SWCNTs for the parallel polarization of the incident light. There is a correlation between the energy levels in ZGNR and armchair SWCNTs Sadykov and Skorkin 2013).
We consider the terahertz radiation generation process in the medium based on the array of the noninteracting zigzag nanotubes and armchair-edge nanoribbons. We are going to get energy distribution on the levels of CNTs and nanoribbons arrays. So, we write differential-difference Eq. (2.6) [see (Sadykov et al. 2016)] in the differential form (Onipko and Malysheva 2017;Malysheva and Onipko 2009; where t is time, current function J(V, T, ) in the nanotube or nanoribbon on the right of Eq. (3.1) consists of electrons with energy . In the further analysis, for convenience, we pass from discrete values to continuous function values ( ) and J(V, T, ) Here ( − k ) is the Dirac delta functions, the domain of energy 0 < < R . Lets consider L = 10 −7 m to be the length of nanoparticles. Equation (3.1) is sufficient for basic state (t = 0, ) = 0 ( ) . Taking into account, that on the left border of the energy space there is the basic state (electrons from other states pass to this state with the photons emission), Eq. (3.1 2.1) is sufficient the first boundary condition (t, = 0) = 0 ( = 0) . On its right border the equation in (3.1) is sufficient for the second boundary condition ∕ | The fact that Eq. (3.1 2.1) satisfies the second boundary condition can be explained by the fact that the spectrum becomes inequidistant far from Fermi energy, i.e. the frequency of the transition Ω between the neighbouring levels will depend on the number of the state (from quantity ). That leads to the fact, that the "diffusion coefficient" significantly reduces so, that value flow reduces on the right border [on the right of Eq. (3.1) we should take into account the dependence of the diffusion coefficient on energy ]. This means, that value ∕ on the right boundary of the domain tends to zero. Therefore we consider great, but finite number of levels and on the right border we will formulate the second boundary condition for Eq. (3.1) (Let us consider that quantity g is independent from ). In the derivation of Eq. (3.1) the equality is taken into account One can obtain the expression for current J(V, T, ) in Eq. (3.1), if the left and right parts in Eq. (2.6) are timed the Dirac Delta function ( − k ) , the obtained equation is summed over index k. As the result, we obtain the expression for current from the congruence (3.2). In this work, the Dirac Delta function is approximated by a Gaussian function or a function with a rectangular profile in energy interval k − ℏΩ ≤ ≤ k + ℏΩ where Δ k and 0 are constants. At t → ∞ , taking into account boundary conditions, the solution (3.1) will take the form of a trapeze where, in the derivation, it was taken into account (3.6) that, firstly, two electrons can be in one state in accordance with the Pauli principle (Payod et al. 2020) and, secondly, there are two types of charge carriers (electrons and holes), which contribute to the considered effect. In the derivation of the asymptotic Solution of Eq. (3.6) we also account for the max value of a function to be max ( ) = 4∕(ℏΩ) in = k o . We integrate (3.1) in the energy space, as a result, we obtain an expression for the flux of quantity in the energy space, taking into account function ( ) (3.6).

Points in the array of elongated nanoparticles for the terahertz radiation generation
Consider two mechanisms for generating current J(V, T, ) on the right side of (3.1). The first mechanism of generating current, considered in Sadykov et al. (2016), is the injection current in the nanotube in a cascade laser, where an array of nanotubes plays the role of an active region [see Fig. 2; the figure is borrowed from ]. The second mechanism for obtaining current is realized through a transition from initial state , associated with stimulated absorption. A CO 2 laser can be used to pump the medium. Next, there is a process of sequential transition from state to initial state � � k ⟩ by analogy with cascade lasers. Transition rate, associated with stimulated absorbtion of photon In the first case, in our model, in Eq. (2.6), total current I(V, T) in a nanoparticle is determined in accordance with the Landauer formula (Berestetskii et al. 1971;Landauer 1972) (3.8) V . Current J(V, T, ) is alone determined in expressions (3.2), i.e. total current equals I(V, T) = ∑ k J k . In the second mechanism for obtaining the current in Eq. (2.6), we use the transition from initial state � � k ⟩ to final state | | | k ⟩ due to the stimulated absorption of infrared radiation, for example, using the radiation of a CO 2 laser. In such a case the transition rate, associated with stimulated emission of photon W (e) and stimulated absorption of photon W (a) with polarization into solid angle do per unit time at the dipole stimulated transition from initial state � � k ⟩ to final state | | | k ⟩ is associated with the total spontaneous transition rates from an initial state to a final subband (Jirauschek and Kubis 2014;Berestetskii et al. 1971) where do is a solid angle element; N is the number of photons of a given polarization, that fall on one oscillator of the field; the number of oscillators d in frequency range d = cdk and in solid angle do is determined by expression d = V 0 k 2 dkdo∕(2 ) 3 ; the sign " × " means the cross product of two vectors; I is spectral intensity in the frequency of electromagnetic radiation with the field polarization vector in the direction of wave vector ; fn = ez fn z is the dipole matrix element of the electron in nanopaticles, = fn = ( f − n )∕ℏ ; transition rate W from initial state �n⟩ to final state �f ⟩ is obtained from Fermi's golden rule, Buttiker (1992), given by The computation of the total transition rates from an initial state to a final subband involves the summation over wave vectors from W function. Integrating (3.9) over all directions of the solid angle, the transition rate, associated with stimulated absorbtion of photon W (a) can be written as when radiation is isotopic and non-polarized ( I doesn't depend on and directions), as a result of expression (3.10)intergrating over do and summation over vector we write the transition rate with the use of radiation spectral intensity I = 2 ⋅ 4 I Transition rate W (a) from initial state �n⟩ to final state �f ⟩ helps us calculate the value of current J k ∕e = W (a) in Eq. (2.6). Later on for calculating total spectral intensity of the radiation we will use the results of Davies (1997), in which the pulse length of a CO 2 laser and the repetition rate of these impulses equal ≈ 1 s and ≈ 10 kHz , density of radiation N ≈ 10 J∕cm 2 .

Effective dielectric function for the crystal structure of nanoparticles
To investigate the THz radiation in 3D two-level resonance structure based on crystal lattices of scattering nanoparticles, one should know the effective dielectric function (EDF) eff of such medium. To define the effective permittivity of a medium, consisting of dielectric particles with dielectric constant 1 , which is also characterized by the depolarization factor n χ , it is necessary to know the dipole moment of the particles [for an ellipsoid see Grachev et al. (2014)] where χ is the individual particle polarizability, 0 is a dielectric constant, χ = (s, p) means that the polarization is parallel (s-polarization) or perpendicular (p-polarization) to the nanotube's axis.
Enter in Eq. (4.1) by analogy with (Landau and Lifshitz 1960) equivalent isotropic dielectric function 1 ( ) for anistropic nanotube, 2 is the dielectric permittivity of external medium. Suppose that the radius of solid cylinder with isotropic function 1 ( ) coincides with the outer radius of nanotube. If this isotropic solid cylinder and the anisotropic nanotube have the same polarizability under the same applied voltage, we call 1 ( ) the anisotropic nanotube's equivalent.
The effective medium-approximation (EMA) is valid for densely and randomly distributed particle composites (Landau and Lifshitz 1960). Maxwell-Garnett's theory (MGT) (Lü et al. 2000;Maxwell-Garnett 1904) is valid for non-densely and randomly distributed composites of particles. In the geometric model of arrays the effective dielectric functions is defined by the individual particle polarizability, taking into account the geometry of the array (the geometrical factor). The morphologies of both the Maxwell-Garnett theory (MGT) and the array model are not consistent with the dense and random character of this nanotube system. The Maxwell-Garnett approximation is used to analyze the propagation of electromagnetic waves in heterogeneous media to determine the effective dielectric constant of different composite structures, for example, in periodic dielectric systems (Wood and Ashcroft 1977;Datta et al. 1993). Comparing with the MGT and the array model, the EMA can describe better the experimental data. The effective medium approximation (EMA), also called the Bruggeman's theory (Tao et al. 1990;Bruggeman 1935), is another method to describe the effective dielectric properties of composites in which the particles of all components are randomly mixed together. This theory has been widely used to explain the dielectric and optical properties of composite materials (Landauer 1978;Niklasson and Granqvist 1984) and proved to be valid at all concentrations. In the case of the array model this model is used for an ensemble of separated quantum dots, which are stacked vertically with the finite period (Brouers and Berthier 1986).
Since volume fraction f is very small, one, basing on the mean field theories, can write effective dielectric function ,MG eff ( ) with the use of the Maxwell-Garnett formula (Landau and Lifshitz 1960; Savinskii and Khokhryakov 1997; where according to the mean field theories average electric field E av and electric displacement D av can be written as , Page 11 of 28 36 The average inner electric field in a particle can be written with the use of the Laplace's equation Maxwell-Garnett formula (4.2) mainly used for dilute systems. If volume fraction f is not small, one can get the effective dielectric function on the basis effective-medium approximation (EMA) (Landau and Lifshitz )1960 In terms of electrodynamics, the discrete structure of array mount of dispersive particles is equal to a continuous dielectric material, characterized by the effective dielectric function (magnetic properties of particles are neglected) (Sadykov and Aporoski 2019;Markel 2016) (see also the section "Appendix" in the article under consideration) where B is a geometrical factor; = s, p ; B s and B p are the geometrical factors, such, that the polarization is aligned along (s-polarization) or perpendicular (p-polarization) to the nanotube's axis). There is B s = 1∕3 when the nanotube centers of a tetragonal lattice are located on the same distance. Ref. (4.6) implies that the polarizability of an individual object is modified by the presence of another via its depolarizing contribution of the electric field acting on this object). Effective dielectric function of the array of nanoparticles depends on the polarizability of individual nanoparticle (4.1) ̃= ∕ V 0 = 4 ∕V and geometrical factor B , which, in turn, depends on polarization of incident light (Markel 2016;Khizhnyak 1986;García-Vidal et al. 1997;Tasaki et al. 1998).
The effective dielectric function formula (4.6) is the Maxwell-Garnett model with the Clausius-Mossotti correction (the arrays geometrical model), since it is based on a rigorous solution of Maxwell's equations under the assumption of a low density of inclusions (Markel 2016;Khizhnyak 1986;García-Vidal et al. 1997;Tasaki et al. 1998;Lakhtakia et al. 1998). Note, that the denominator in the second term of Eq. (4.6) is the Clausius-Mossotti correction responsible for the electromagnetic interaction of the quantum dots).
In order to determine the effective dielectric function using (4.6), one should use the theoretical dependence of the polarizability of an individual nanotube obtained in Slepyan et al. (1999) (CGS-system) where V is the volume of spheroid; at the determination of individual nanotube's polarization it is supposed that SWCNT is the prolate spheroid made of an orthorhombic material); 0 = p √ n ; in correspondence with (Slepyan et al. 1999) for the plasma frequency, there 1∕2 . There is B s = 1∕3 when tetragonal SWCNT the nanotube centers are located on the same distance.

Axisymmetric elongated spheroids
Elongated nanotubes can be approximated by the axisymmetric elongated spheroids, while greater semiaxis equals to the half of height nanotube c = L∕2 , transversal sizes of spheroids a = b is taken due to the equality condition of the volumes of particle and spheroid. In this case the depolarization factor along the semi-major axis have form (Markel 2016) [see Grachev et al. (2014) or Sadykov and Aporoski (2019)] where J 3 is the Newtonian potential. The polarizability of spheroid is determined in agreement with (4.1).

Elongated spheroids
Elongated nanoribbons can be approximated by the axisymmetric elongated spheroids, while greater semiaxis equals to the half of height nanotube c = L∕2 , transversal sizes of spheroids a and b is taken due to the equality condition of the volumes of particle and spheroid). The depolarization factor along the semi-major axis have form The polarizability of spheroid is determined in agreement with (4.1).

Calculation results
Based on Eq. (3.1) for the diagonal elements of the density matrix, which can be written in the differential form, we numerically get a solution for function (n, ) . Suppose in the numerical solution of Eq. (3.1) the energy definition area as 0 ≤ ≤ R , where R = R 0 , 0 = 6ℏΩ = 0.108 eV ; constant R has two values R = 2 and R = 4 . On the left border of the energy space for the quantity ( ) in (3.1) there is the first boundary condition ( = 0) . On the right border of the energy space for the quantity ( ) there is the second boundary condition ∕ | | = R = 0 . Suppose the hopping integral to be equal to 0 = 2.7eV ; L = 10 −7 m , | | = 3.24 × 10 −27 Q × m , | | | E | | | = 2.2 × 10 5 V∕m , T 1 = T 2 = 3 × 10 −12 s , kT = 0.0258 eV ( T = 300 K ). The transition frequency between the adjacent levels is Ω = 2.7 × 10 13 s −1 . The Lorentzian line shape function g( ) is defined in Eq. (2.6). In Eq. (3.11) and (3.2) function 0 ( ) and the eigenvalues of energy states of charge k carriers are in form of (4.8) Page 13 of 28 36 where F is Fermi velocity. Write the initial condition of the differential Eq. (3.1) is in form of (t = 0, ) = 0 ( ).
In the case of the multicascade mechanism we suppose the magnitude in nanotubes in Eq. (2.8) to be equal to J k 0 = 8 × 10 −6 A . The magnitude in nanotubes, which can be used as a non-linear amplification region of radiation in cascade lasers, was borrowed from Sadykov et al. (2016). In Sadykov et al. (2016) current for the ballistic transport is determined in accordance with the Landauer formula (3.8) (Berestetskii et al. 1971;Landauer 1972).
In Figs. 1 and 2 there is the dependence of function ̃(t, ) on electron energy and time t in R = 2 and R = 4 , where ̃(t, ) = (t, ) − 0 ( ) , where function ̃(t, ) is valid for the equation with the the initial condition and boundary conditions (5.1) where delta function ( − k ) in (5.2) and (5.1) are approximated by the Gauss function from (3.5); magnitude Δ k in the Gauss function is equal to Δ k = ℏΩ∕5.
From the calculation results shown in Fig. 1, it follows that in the case of the multicascade mechanism of radiation generation at R = 2 the length of the unsteady mode region is equal to 0 ≤ t ≤ 0.3 × 10 −11 s . The solution of quantity ̃(t, ) in the length of steady mode region t ≥ 0.3 × 10 −11 s has a form of plateau, i.e. the solution doesn't depend on time. From the calculation results shown in the inset of Fig. 1, a good match of the calculation results at moment t = 10 −11 s at different temperatures T = 150, 300, 600 K with analytical solution (3.6) follows. From the calculation results shown in Fig. 2, it follows that at R = 4 the length of the unsteady mode region is equal to 0 ≤ t ≤ 0.8 × 10 −11 s . The solution of quantity ̃(t, ) in the region of the length of steady mode region) t ≥ 0.8 × 10 −11 s has a ridge-shaped form.
In the mechanism of generation of radiation, when the current is obtained during the transition from the initial state to the final state, associated with stimulated absorption, we will assume that the current determined in (2.7) will have the dependence where k 0 = 6 , k 1 = 0 . In (3.11) we wrote transition rate W (a) with the help of radiation spectral intensity I = 2 ⋅ 4 I . The same dependence for radiation spectral intensity is obtained in Jirauschek et al. (2014) for the two-level system, in which dependence of spectral intensity I is expressed with intensity I 0 of radiation I ≈ I 0 ∕Δ ∼ I 0 T 2 .
For calculating the radiation spectral intensity we use the results from Davies (1997), in which the pulse-length of a CO 2 -laser and the repetition rate of these impulses are equal to ≈ 1 s and ≈ 10 kHz , density of the radiation N ≈ 10 J∕cm 2 . Thus I 0 = N∕ ≈ 10 11 J∕(m 2 × s) . Taking into account the expression for the dipole matrix element of laser transition   In calculations we consider the magnitude hoping integral to be 0 = 2.7 eV , L = 10 −7 m , as in the case of the multicascade mechanism of radiation generation. In contrast to the case of the multicascade mechanism of radiation generation, the magnitude of current and field intensity in Eq. (2.8) are equal to J k 0 = 1.48 × 10 −5 , In Figs. 3 and 4 there is dependence of function (t, ) on electron energy and time t at R = 2 and R = 4 . In calculations Eq. (5.2) was being numerically solved, taking into account given boundary conditions and initial condition ̃(t = 0, ) = 0 , where (t, ) =̃(t, ) + 0 ( ) . The inset of Fig. 3 shows dependence (t, ) at time moments t = 0.004 × 10 −11 s (curve 2) and t = 10 −11 s in the steady mode region (curve 1). Curve 3 corresponds to the analytical solution (3.6). From the calculation results shown in the inset of Fig. 3, it follows that a quantity at initial moment of time t = 0.004 × 10 −11 s (curve 2) reduces because of a transition from the initial state � � k 1 ⟩ to the final state � � k 0 ⟩ due to the stimulated absorption of infrared radiation.
Curves 1, 4 and 7 (in green) in Fig. 7 correspond to the case of = 3 × 10 11 s −1 ; curves 2, 5 and 8 (in blue) correspond to the case of = 3 ⋅ 10 12 s −1 s −1 ; curves 3, 6 and 9 (in red) correspond to the case of = 3 × 10 13 s −1 . Solid-line curve 10 in Figs. 7 and 8 corresponds to frequency 0 ∕2 ≈ 2.56 × 10 13 s −1 ; vertical dashed-line curves 11 and 12 in Fig. 7 correspond to two radiation frequencies of a CO 2 laser: c∕ 1 and c∕ 2 , where 1 = 10.6 mμ and 2 = 9.4 mμ. Depicted in the inset of Fig. 7 curves are the continuation of curves 1-9 in the frequency interval. The results of calculations have been obtained with the use of Fig. 8 Dependence of a real part of effective dielectric function eff on frequency, based on the MGT at volume fraction f = 1.44 × 10 −3 ; 3.9 × 10 −5 ; 1.44 × 10 −6 . The length and diameter of zigzag nanotubes (9, 0) are correspondingly L = 10 −7 m and d = 0.7 nm Fig. 9 Dependence of an imaginary part of effective dielectric function eff on the CNT array frequency, based on the Maxwell-Garnett theory with the Clausius-Mossotti geometric correction [ordered array of nanotubes]. The length and diameter of zigzag nanotubes (9, 0) are correspondingly L = 10 −7 m and d = 0.7 nm the Maxwell-Garnett formula (4.2). In Fig. 8 the value of function Re eff at ∕2 = 10 12 correspond to the dotted-, solid-and dashed-line curves and are equal to 5.38 , 1.12 , 1.004.

Discussion of results
The results of our research are of interest in view of the recent development of a dense array of carbon nanotubes (Slepyan et al. 2010;He et al. 2016;Komatsu et al. 2017;Gao and Kono 2019). Recently, a method of controlled vacuum filtration (CVF) has been developed. This method makes it possible to obtain crystalline films from single-walled CNTs. In 1 cm wafers the nanotubes are almost perfectly aligned and maximally packed (∼1 nanotube per cross-sectional area of 1 nm 2 (Slepyan et al. 2010;He et al. 2016;Komatsu et al. 2017). In Gao and Kono (2019) we see the results of a series of systematic experiments. They demonstrate that the CNT alignment direction can be controlled by the surface morphology of the filter membrane used in the vacuum filtration process. The results, obtained with the CVF method, allow us to consider the problem with a dense array of carbon nanotubes. In our research the value of the used density of nanotubes in the array (∼1 nanotube per cross-sectional area of 100 nm 2 ) is by three orders greater than in Sadykov et al. (2016)), but by two orders lower than Slepyan et al. (2010) and Gao and Kono (2019). We choose of such density, because high density leads to the strong interinfluence of nanotubes in the array. As a result, the parameters of nanotubes in the array greatly change compared to a single nanotube (Sadykov and Jolnirov 2021;Shuba et al. 2007;Wang et al. 2004;Kempa et al. 2007). This leads to the grand challenges in nanoscience and nanotechnology. For example: the creation of macroscopic devices by assembling nanotubes, maintaining unique one-dimensional properties of nanotubes while producing large-scale architectures of aligned SWCNTs (Kempa et al. 2007;Ma et al. 2011). According to the results of Sadykov and Jolnirov (2021) (see Fig. 5; Curve 2 [in blue]) the ratio of the nearest distance between CNTs in the array to the nanotube length in the array is L = 0.1 μm. It means, that if we have density of nanotubes in the array, when one nanotube per cross-sectional area of 100 nm 2 , the delay coefficient and the frequency of geometric resonance change slightly.
For the effective dielectric function three basic models are employed: effective media approximation (Bruggeman's theory) [nanotube bundles], Maxwell-Gatnett theory [ordered array], and Maxwell-Garnett theory with the Clausius-Mossotti geometric correction [ordered array of nanotubes]. The calculations show, that if the radiation wave length is = 2 c∕Ω, Ω = 2.7 × 10 13 s −1 , the dielectric permittivity eff of the medium, based on noninteracting parallel-oriented elongated SWNTs, is approximately equal to eff ∼ (1.5 + 0.1i) . In the frequency range, that is of interest to us, the number of the imaginary and real parts of the effective dielectric constant has reasonable values. To calculate the energy of radiation with stimulated emission of photon or with stimulated absorption of radiation, one should know the function dielectric permittivity, the longitudinal time relaxation T 1 and the transverse relaxation time T 2 in Eqs. (2.2) or (3.1). The dielectric permittivity determines the relaxation time of electromagnetic energy. Values T 1 , T 2 , c depend on the relaxation time = 1∕ in the kinetic Boltzmann equation. Generally, the changing range of varies from 3.3 × 10 −14 to 3.3 × 10 −12 s. For example, in Liu et al. (2013) the relaxation time ≈ 1.4 × 10 −12 s was obtained at the room temperature (T = 300 K).
In Jishi et al. (1993) the value ≈ 3 × 10 −12 s was experimentally obtained, that is in a good agreement with Liu et al. (2013). It was experimentally discovered in Tans et al. (1997) [see also Slepyan et al. (1999)], that the relaxation time is = 3.3 × 10 −13 s. This value is calculated using the formula [see also Sadykov and Aporoski (2017) where m c is the mass of the carbon atom, ℏ L ≈ 0.113 eV is the energy of the longitudinal acoustic phonon, q 0 ≈ 2.5 A −1 = 2.5 ⋅ 10 10 m −1 , J 0 = 2.2 eV. For k B T = 0.025 eV these values give. sinh In Borondics et al. (2006) as in Liu et al. (2013) the value of absorption (emission) probability per time unit of the longitudinal (transverse) phonon per unit of time and the length of the helical trajectory was calculated. In Borondics et al. (2006) the value of τ is approximately equal to ≈ 2.2 × 10 −12 s. If the temperature of the medium, where the nanotubes are situated, reduces k B T = 0.0125 eV, according to (6.1) we will get ≈ 1.3 × 10 −10 s. With such values of the actual task of terahertz emission generation in the array of elongate particles with the help of short high-intensity nanosecond electric pulses with a subnanosecond leading-edge width (longitudinal plasmon-polaritons) (Kibis 2001;Ginzburg et al. 2015). In this paper we, taking into account the results of Liu et al. (2013) and Jishi et al. (1993), consider that T 1 = T 2 = 3 × 10 −12 s.
The value of the dipole moment is Which is given after Eq. (2.2). For f = , n = k − 1 the value of the matrix element of the dipole moment = (d z ) k,k−1 = e(z) k,k−1 = −e2|e|L∕ 2 is commensurate with the value of the value of the matrix element of the dipole moment which is expressed through the velocity matrix element in case of interband THz transitions in narrow-gap carbon nanotubes and GNRs, where d is the electric dipole moment. Bandgaps themselves in nanotubes are induced by the magnetic field through the curvature in quasimetallic nanotubes through the edge effects in armchair nanoribbons (Portnoi et al. 2015;Shahnazaryan et al. 2019). In fact, the peak value of the matrix element of the velocity operator, because of the transitions between the top valence subband and the lowest conduction subband of the order of the Fermi velocity, does not depend on bandgaps widths. It can be seen that in this case the matrix element of the velocity operator on the transitions �v⟩ → �c⟩ in the vicinity of the Dirac point has a universal value equal to the Fermi velocity of the electron zigzag carbon nanotube. We see analogic results in our research, if we consider transitions �k⟩ → �k − 1⟩ . For f ≠ k + 1 the matrix element of the velocity operator is transformed into the matrix element of the electric dipole moment (6.2) of zigzag nanotubes and armchair nanoribbons .
For matrix elements in the case of optical transitions between the lowest conduction subband and the highest valence subbands of zigzag quasimetal CNTs, there are simple analytical expressions in the presence of curvature and magnetic field applied along the axis of the nanotube Mesyats et al. 2020) where z = −2 cos(2 ∕3 ±f ∕n) , f = Φ∕Φ 0 , Φ 0 = 2 ℏ∕|e| is the magnetic flux flowing through the polygons cross section sin(f ∕n) ≈f ∕n = √ 3b∕(2L) ≪ 1 . The magnetic field creates a magnetic flux 3bn is the diameter of the zigzag CNT, d∕L ≪ 1 . From (6.4) we definitively get the expression for the matrix element of the velocity operator It is commensurate with the matrix element of the velocity operator from (6.3).
To calculate the energy of radiation with stimulated emission of photon or with the stimulated absorption of radiation one should know the function of dielectric permittivity. For wave length = 65 μm, ≈ Ω = 2.7 × 10 13 s −1 terahertz emission fulfils condition ∕(2 ) ≈ 0.43 × 10 13 s −1 . From the calculation results shown in Figs. 5,6,7,8,9 and 10, it follows that terahertz radiation is far from maximum point 0 = 1.6 × 10 14 s −1 for function Im eff ( ) . Obtained in the work parameters of the medium consisting of the extended nanoparticles, based on the arrays geometrical model are interesting from the position of the short terahertz impulses generation (Kibis et al. 2007;Belonenko et al. 2006;Zhukov et al. 2016).

Conclusion
The results obtained in this paper show that in the case of the multicascade mechanism of radiation generation in the medium, based on noninteracting parallel aligned CNTs can be used as a non-linear amplification region of radiation. The process of the terahertz amplification with a QC-laser in the medium, based on the armchair-edge nanoribbons and the zigzag nanotubes with the metallic conductivity is theoretically studied and numerically simulated in the paper. The quantity of injection current in the nanopaticle points is rated. The task is generalized for analogous mediums in the case where radiation current can be generated with the help of the process of infrared radiation stimulated absorption, for example, using the radiation of a CO 2 laser.
In the case of the multi cascade mechanism of radiation generation in a continuous wave mode the field is equal to | | | E | | | = 2.4 × 10 5 V/m (see Sect. 3.7), where taking into account (Liu et al. 2013;Jishi et al. 1993) we consider that T 1 = T 2 = 3 × 10 −12 s. The quantity of the field intensity is the result of the stream (3.7) quantity conservation in a continuous wave mode (consequence of the charge conservation law). If concentration of nanotubes is equal to n 0 = 10 10 cm −2 = 10 14 m −2 the distance between nanotubes is 100 nm, then the output power per unit area will be about Q ∼ n 0 k 0 ℏΩJ k 0 ∕e ≈ 860 W/ cm 2 , where J k 0 = 8 × 10 −6 A, k 0 = 6 , 0 = k 0 ℏΩ = 0.11 eV. In the task of terahertz radiation generation with the use of a QC-laser in the medium, based on the armchair-edge nanoribbons and the zigzag nanotubes with the metallic conductivity, the value for transition rate W (a) ≈ 9.26 × 10 13 s −1 with stimulated absorption of radiation of a CO 2 laser (it corresponds to the injection current J k 0 = −J k 1 ≈ 1.48 × 10 −5 A in nanoparticles) has been obtained. In the discussed case the current quantity is 1.85 times larger then in the case of the multicascade mechanism of radiation generation. Field intensity in a continuous wave mode the field is equal to | | | E | | | = 3.26 × 10 5 V/m. From given in the paper numerical results for the effective dielectric function it follows, that for the terahertz radiation with wave length = 65 μm ( ≈ Ω = 2.7 × 10 13 s −1 ) condition ≪ 0 , where 0 = 1.6 × 10 14 s −1 , is executed). It means, that frequency of terahertz radiation is far from maximum point 0 of function Im eff ( ) . In this region of frequency the quantity of imaginary and real parts of the effective dielectric function has reasonable values.

Appendix
In terms of electrodynamics, the discrete structure of massive amount of dispersive particles is equal to a continuous dielectric material, characterized by the effective dielectric function (magnetic properties of particles are neglected) [see Markel (2016)] hereinafter, the axis of nanotube is directed along the z axis, where Ω = ΔxΔyΔz , Δx, Δy, Δz are the distances between nanotube centers along the axes x , y and z, correspondently; the production of operators ̂⋅ĝ paper is equal to ̂ĝ = ii ⋅ g ii , where i = 1, 2, 3, where B s and B p are the geometrical factors, such, that the polarization is aligned along ( s -polarization) or perpendicular ( p-polarization) to the nanotube's axis).
In the case where the tensor of dielectric permittivity of ellipsoid is diagonal and main direction of this tensor coincide with the main axes of ellipsoid, the matrix g ik is also diagonal ( g ij = 0, i ≠ j ) and have form [see Markel (2016)] where J 1 = J 100 , J 2 = J 010 , J 3 = J 001 , n i -is the depolarization factor (Grachev et al. 2014).
From (4.1) and (A3) it follows that g ii coincide with the polarizability i From (A1), taking into account (A2) and (A3), we obtain the expression for the effective dielectric function where f = 4 abc∕(3Ω) is the volume fraction; B is a geometrical factor, = s, p ; B s is such a geometrical factor, that the polarization is aligned along (s-polarization), B p is such (A1) eff = 2 1 + , R( ) = (a 2 + )(b 2 + )(c 2 + ),