Density matrix superoperator for periodic quantum systems and its application to quantum cascade laser structures

In this work we present a generalization of the Liouvillian superoperator for periodic quantum systems that can be formulated through partitioned Hamiltonians. We formulate a compact algebraic form of the superoperator that allows efficient numerical implementation along with the possibility of further generalization and the inclusion of the system’s boundary effects (i.e. device contacts). We apply this formalism to Quantum Cascade Laser structure where we compare the second nearest and the nearest on approximation, and present the laser dynamics that is independent from the number of states considered.


I. INTRODUCTION
The Density matrix (DM) formalism was first introduced by J. Von Neumann 1-4 and its applications span through a variety of fields in Quantum Mechanics. The time evolution of the density matrix is, generally, described by the master equation in Lindblad form [5][6][7][8] and one of the key aspects of this equation is that it results in a nonlinear system of equations, if written in the matrix form. However the system can be linearized by introducing the Liouvillian superoperator and the goal of this paper is to present a mathematical formalism that can generalize and simplify the superoperator of the systems with high symmetry, in particular periodic quantum systems, such as Quantum Cascade Lasers (QCLs) and graphene. 9 QCLs 10 are powerful semiconductor sources of coherent radiation in the mid-infrared (MIR) 11 and Terahertz (THz) band 12 with potential applications in free-space communications, medical diagnostics, and chemical sensing. [13][14][15][16][17] These devices use sequential tunnelling and usually comprise a large number of semiconductor heterostructure periods (typically GaAs/AlGaAs for THz QCLs).
The DM model has been successfully applied to QCLs [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] and it represents a quantum model 35 that can describe key aspects of the underlying physics while retaining low numerical complexity when compared to the more extensive models, such as the Non-Equilibrium-Green Function (NEGF) approach. [36][37][38] The most common applications of the DM model are for systems with a few states, which usually yield analytic expressions. 4 The Liouvillian superoperator enables the generalization to systems of any size, however, treating periodic systems usually allows elimination of a large number of equations due to the system symmetry, and to our knowledge, compact form of the Liouvillian has not been considered in detail.
In this work we develop a generalization of the Liouvillian for systems that can be represented by partitioned Hamiltonian, which is commonly employed in the tight-binding models. In section II we formulate the superoperator for symmetric (2M + 1) block diagonal Hamiltonian and the density matrix that corresponds to M-neighbour approximation and also discuss the common modelling approaches and approximations that are usually introduced in DM models. Section III focuses on QCL devices and introduces dynamical aspects to the model we presented in Ref. 33. Section V presents concluding remarks, and in the Appendix we discuss algebraic derivations and present the most general formulation of the superoperator for any system that can be partitioned.

A. Liouvillian superoperator
The DM formalism describes a quantum system as a statistical ensemble of quantum states through an operatorρ = |ψ⟩⟨ψ| which, provided with the corresponding basis of wavefunctions |ψi⟩, results in a matrix where each element ρij represents coherence of states i and j. The time evolution of the density matrix is then given by the Liouville equation: which represents the quantum mechanical analogue of the equation of motion in the classical systems. It can be further generalized by adding effects that cause decoherence of the system in a form −( dρ dt ) relax . This depends on the problem of interest and we will discuss these additional terms in the next section.
The formulation in Eq. (1) employs interaction of two operators, and outcome of such interaction is usually referred in the literature as a superoperator. The actual mathematical need for the introduction of the linear superoperator in Eq. (1) arises from the fact that the unknown is a matrix, and that in algebraic sense Eq. (1) is not linear. In mathematics, systems that take the form AXB = Y are linearized as A ⊗ B T X ′ = Y ′ where ⊗ represents the Kronecker tensor product and X ′ and Y ′ are vectorised forms of the original matrices X and Y, respectively, unpacked row by row in a column vector (for column by column unpacking, the linearization reads A T ⊗ B). We will keep the notation with an apostrophe' to refer to these vectorised forms throughout this paper. In our case the commutator linearizes in a form L = H ⊗ I − I ⊗ H T where I is an identity matrix of the same size as the Hamiltonian H. The linear operator L is called Liouvillian superoperator and the linearization of Eq. (1) in the form dρ ′ dt = − i ̵ h Lρ ′ is a well known formulation of the density matrix superoperator in the literature. 3 We also note that if we had P states in the basis of wavefunctions of the overall quantum system, the superoperator would be a matrix of P 2 × P 2 size and ρ ′ would be a column vector of length P 2 .
In many practical cases of interest, we need to deal with Hamiltonians which are partitioned in the block form, and may have high degree of symmetry, especially in cases which deal with periodic systems, where the wavefunction basis can be taken only on one period and not from the entire quantum system. In Appendix B we show a simplified algebraic way of forming a superoperator if the Hamiltonian and the corresponding density matrix are partitioned in block form without any particular symmetry. In further considerations we will focus on Hamiltonians with high level of symmetry, often employed in tight binding 39 and the nearest first or second neighbour approximations.
Consider a periodic quantum system in which we can formulate a wavefunction basis from one period and describe the Hamiltonian of that period as H 0 and its corresponding density matrix as ρ 0 . Note that if the wavefunction basis has N states, H 0 and ρ 0 are N × N matrices. The Hamiltonian H 0 can interact with the adjacent periods (as depicted in Fig. 1) given by block Hamiltonians H 1 , H 2 , . . . HM (to the right) and H− 1 , H− 2 , . . . H−M (to the left) which have the corresponding block density matrices ρ± 1 , ρ± 2 , . . . ρ±M, where M is the number of neighbours that central (H 0 ) period has in either direction. Take a system consisting of Q periods (Q ≥ 2M + 1). The Hamiltonian of the entire system can then be written as (2M + 1)-diagonal block matrix of Q × Q block size. In order to apply periodic boundary conditions, Q → ∞ and Q × Q, (2M + 1) -diagonal H and ρ need to be substituted into the Liouville equation. This yields (4M + 1) -diagonal block matrix from the commutator term in (1). This matrix will be equal to the derivative of the density matrix (times i ̵ h) and it is clear that there will be some "extra" equations that equal to zero block matrices in the overall density matrix. This occurs due to the algebraic properties of banded matrices. Generally, product of S -diagonal matrix with T -diagonal matrix is a S + T − 1 -diagonal matrix. In our case, if both H and ρ are (2M + 1)diagonal, the result of their commutator will be a (4M + 1)-diagonal matrix, and there will be 2M equations that target zero blocks in the density matrix. We can therefore group equations that Eq. (1) yields into two groups: Group 1: Group 2: Group 1 has (2M + 1) block equations, while in group 2 there are 2M additional equations which contain the same unknown density matrix blocks as group 1. The equations in Eq. (3) are the consequence of the duplication of information in the system. We assumed the form of the resulting density matrix in advance, and therefore forced the system to yield (2M + 1) equations in Eq. (2). The reason why additional 2M equations emerged is that even though we are forcing the system to have no terms for interaction with (M + 1) − th neighbour, the included neighbours (2M of them) still give their contribution to the next M + M neighbours and therefore we end up with additional 2M equations. They originate from the duplication of information from initial symmetry of the Hamiltonian and therefore can be discarded.
The number of neighbours that need to be included in the system depends on the underlying physics of the system. It is also possible to include more density matrix blocks then the number of partitions in the Hamiltonian. However this is only relevant in strongly coupled systems and in systems where the dissipator term may coherently couple periods that are not coupled by the corresponding Hamiltonian block.
Let us consider a system with the second neighbour approximation given by the Hamiltonian and the corresponding density matrix: and Let us apply general expressions in Eq. (2) and neglect Eq. (3), for M = 2 which corresponds to the Hamiltonian and the density matrix given by Eqs. (4) and (5): Formulating the system of equations in Eq. (6) can be done intuitively: for a given block ρi, equation for that block will consist of all combinations of the Hamiltonian and the density matrix blocks whose indices add to i, provided that the blocks are labelled as in this work.
The primary goal of this work is to simplify Eq. (6) (and Eq. (2)) and provide linear superoperator that can be easily numerically implemented. Equation (6) represents a system of commutator equations where we can define sub-Liouvillian operators which linearize each commutator as Li = (Hi ⊗ I − I ⊗ H T i )ρ ′ i , i = −M, . . . M. This then turns Eq. (6) into a linear system which can be written as i ̵ h dρ ′′ dt = LQρ ′′ , where LQ is a (2M + 1)-diagonal block matrix whose block elements are Li and it has a similar form as the original block form of the Hamiltonian. The unknowns will represent vectorised forms of the corresponding density matrix blocks, labelled as ρ ′ i , and we can pack these vectorised vectors in reverse order i = M, . . . − M into one vector, labelled as ρ ′′ . It is interesting that we can form linear operator for system in Eq. (6) directly from the Hamiltonian by using Khatri-Rao type of matrix product (definition and derivation is given in the Appendix A) denoted by symbol ⊠.
Equation (6) linearizes as: where I N U Q is a block matrix partitioned in the same way as the Hamiltonian, where each block is an identity matrix. Note that the second term in Eq. (7) is "dot" transpose operation which transposes only the partitions within the Hamiltonian, not the Hamiltonian itself. The important advantage of Eq. (7) is its mathematical simplicity and the similarity to the general form of density matrix superoperator, since the only difference is in the type of algebraic product (and I N U Q matrix and the "dot" transpose operation). This formulation is general for periodic quantum systems where the Hamiltonian and the corresponding density matrix can be partitioned in block form (it applies to Eq. (2) as well).
The main advantage of the formulation in Eq. (7) is that it neatly applies periodic boundary conditions and linearizes the system in the compact algebraic form. An additional difference from the Liouvillian superoperator L = H ⊗ I − I ⊗ H T is that the superoperator in Eq. (7) packs the unknowns of the system differently and requires the formulation of the system with partitioned Hamiltonian. Generally, superopetaor for any partitioned Hamiltonian can be formed by using Khatri-Rao product and this is further explained in Appendix B.
The linearization in Eq. (7) uses the Hamiltonian of (2M + 1) × (2M + 1) block size, however once Eq. (7) is solved, it is incorrect to use such Hamiltonian and the corresponding density matrix for finding the expectation values of operators. The expectation value of any operator Ô in DM formalism can be found as Tr(Ôρ). Trace operation for the partitioned matrix results in a sum of traces of each submatrix, and the product Oρ needs to be dound as a limit value when Q → ∞ with an infinite Q × Q, (2M + 1) -diagonal density matrix and the Hamiltonian as discussed in Ref. 33 in order to ensure the implementation of periodic boundary conditions.

B. Dissipator
Liouville equation (1) describes only the interactions included in the Hamiltonian, however, it is quite common to model some interactions separately by adding phenomenological relaxation terms that damp the equation of motion in the form −( ρ τ ) relax or, in a more detailed form −γmn(ρmn −ρ (eq) mn ), m, n = 1, 2, . . . , N, where γmn are decay rates that damp density matrix terms ρmn to their equilibrium value ρ (eq) mn , where 4 ρ eq mn = 0, m ≠ n. Furthermore, the state populations (described by the diagonal elements of the density matrix) decay can be also introduced in the form ∑ E n >E m ρ nn τ mn − ∑ E n <E m ρ mm τ nm and it applies γmn = 1 2 ( 1 τ m − 1 τ n ) + γ col mn , where τm and τn are total decay rates of populations out of levels m and n, while γ col mn are the dephasing rates due to the processes (i.e. elastic collisions) that do not affect the state populations, and influence only the phase, these rates are often referred to as proper dephasing rates. 4 We will refer to all these additional terms in Liouville equation as a dissipator D = −( ρ τ ) relax of the system; this term follows the symmetry of the Hamiltonian described in the previous subsection, but does not follow a neat algebraic formulation as Eq. (7), even though these terms originate from various interaction Hamiltonians that could have been included in the total Hamiltonian H. This term is linear and we can transform it as D → D ′′ ρ ′′ , however the form of D ′′ can be derived by mathematical induction and it depends on the system considered. In the Appendix C we show how D ′′ is derived for QCLs. 25,26,33 Equation (1) now acquires an additional term D and Eq. (7) can be written as: The reason why D is often excluded from H is because of possible simplifications and physical interpretation of the quantum system under consideration. The Hamiltonian of the system usually consists of several different interaction Hamiltonians that are of interest, however some of them can be considered as perturbations and their effect can then be handled by perturbation theory (i.e using Fermi's golden rule). For example, in laser systems, the rates τmn can be directly expressed via Fermi's golden rule, and τm would simply represent state lifetimes. It is important to note that D must not break the positivity of the density matrix since the main requirement for a physical solution of Eq. (1) is that ρ is a positive semi-definite matrix which ensures that all diagonal elements of the density matrix are positive (which must physically be satisfied). This condition will be satisfied if various scattering processes are included in D in Lindblad 5 form, and Fermi's golden rule satisfies this condition. Note that carious forms of dissipators can be found in the literature 9,30-32 that may also be applied.

C. Dipole approximation
The density matrix formalism is capable of describing optical macroscopic effects of the quantum system along with the effects of external electrical perturbation.
We can assume that the Hamiltonian can be split as H = H free + V where H free is the Hamiltonian of the unperturbed system, and V is the energy of interaction. The dipole approximation physically models the system as a dipole antenna that will resonate at a specific transition frequency, in this case the energy difference of the levels, and assumes that the interaction energy can be written as V = −μA, where μ = −er is the electric dipole moment operator 4 and A(r, t) is the optical electric field.
The most common implementation of the DM formalism describes H free as a diagonal tight binding matrix (filled with the corresponding energy states Em) and the substitution into the general Liouvillian yields a system of equations in the form: 4 This equation is actually a closed form of the Liouvillian (L = H ⊗ I − I ⊗ H T + D), where H = H free + V and it represents the most common DM formulation 4 that is often combined with Maxwell equations. The formulation we presented in Eq. (8) is not equal to the one in Eq. (9). It is somewhat valid to claim that Eq. (8) can be derived from Eq. (9), but it is important to point that the Hamiltonian is different, because it includes interactions between the periods and not just the tight-binding Hamiltonian. Additionally, the order in which both systems of equations are written is different. In Eq. (9) the density matrix ρ is unpacked row-wise, and the first two terms simply correspond to the algebraic transformation of the Kronecker product (L = H ⊗ I − I ⊗ H T ) and periodic boundary conditions are not applied. In Eq. (8) the density matrix was partitioned in blocks as in Eq. (5) and the overall system stores the unknowns by unpacking each block row-wise. For example, in Eq. (6) the overall density matrix first unpacked block ρ 2 into N × N column vector, then ρ 1 , ρ 0 . . ., and then joined them into one 5N 2 column vector which is referred to as ρ ′′ in Eq. (7). This somewhat cumbersome repacking in Eq. (8) resolved the problem of deciding which equations in the periodic system can be ignored, how to properly set the periodic boundary conditions and, most importantly, offered the convenience in physical interpretation of the system that is described by a partitioned Hamiltonian and the corresponding density matrix.
Furthermore, the approach we present in Eqs. (4,5) can be further generalized, in order to include the finiteness of system considered (i.e. contacts of the laser structure), which is one of the main advantages of the more general models. 36

D. Non-rotating wave approximation
In most cases of interest the optical electric field A has a physically expected form which can greatly simplify the DM model. Generally, A(r, t) should be obtained from Maxwell equations coupled to the DM model via polarization, which is given as P = ND Tr(μρ), where ND is the density of electrons in the system. However, one of the most common approaches in handling the Maxwell equations is to look for a solution in the particular harmonic form that is expected from a system. In case of a laser (or absorber) structure, it is expected that the optical field has a dominant term at a particular frequency. A variety of approaches that take different forms of A(r, t) were used, 26,29,33 and we will discuss some of them. A very common approach is to assume that the interaction V is weak (its matrix elements are much smaller then elements of H free ). In this case Eq. (9) can be treated by perturbation theory, and a detailed study can be found in ARTICLE scitation.org/journal/adv literature, in particular for two-level systems. 4 Another common approach treats Eq. (9) by rotating-wave approximation (RWA), however for compact usage in Eq. (8) we will focus on a more general approximation from which RWA was originally derived. Assume that the optical electric field has a general plane-wave form: where β is the propagation constant. Without loss of generality, we will assume that the optical field has only the time dependence, as given by the second equation in Eq. (10). Both of these expressions are used in treatment of DM and Maxwell equations, but the second form can be taken in cases where individual quantum systems are much smaller them the wavelemgth of light. Non-rotating wave approximation (NRWA) assumes that the system's response will have identical harmonic form, while RWA discards some terms based on the underlying physics of the system. Note that implementing the RWA is actually more complicated because we would need to manually select which terms of the density matrix would be allowed to oscillate at +w and which at −w, depending on the energy states that correspond to the particular density matrix elements. NRWA does increase the numerical complexity, however it significantly simplifies further formulation of the model.
The drive terms in differential equations (8,9) have a harmonic form due to Eq. (10), and it would appear natural to assume that the solution will follow this structure, i.e. ρ(t) = ∑ k ρ + k (t)e iω k t + ρ − k (t)e −iω k t . However, this is possible in some special cases when the underlying physics allows it. The generalization of the system in Eq. (8) is far from trivial.
Interestingly, the generalization of the effect of Eq. (10) on DM formulation in Eq. (8) is possible when only one frequency component is considered, and it can be expressed by a simple algebraic formulation.
To illustrate that, we consider H = H dc + H ac e −iωt + H ac e iωt and ρ = ρ dc + ρ ac − e −iωt + ρ ac + e iωt , where each term is partitioned and has the (2M + 1) diagonal form as in Eqs. (4,5). Therefore, each block Hj and ρj has three harmonic terms: Hj = H dc j + H ac j e −iωt + H ac j e iωt and ρ k = ρ dc k + ρ ac − k e −iωt + ρ ac + k e iωt . Consider the equations resulting from the commutator of Hj and ρ k : Each commutator in Eq. (11) can be linearized with the corresponding sub-Liouvillian, and we will only have two of them: and the linear form of Eq. (11) is: The consequence of Eq. (12) is that it affects each equation in Eqs. (2,3), and it breaks each of them into three equations. This causes the system size to increase three times, however the approach we presented in Eq. (8) needs to be just slightly re-modified. In Eq. (12) we can notice that linearization of Eq. (11) can be written in a similar Khatri-Rao notation as used in Eq. (8). We need to apply Eq. (12) to each equation in Eqs (2,3), however we can notice that the structure of Eq. (12) is tridiagonal and we can transform the initial Hamiltonian by defining an expansion rule for non-rotating wave approximation (NRWA) in the form: (13) is applied to each block of our initial Hamiltonian, the expanded Hamiltonian enlarges 3 times, but the formulation of Eq. (8) is still correct and slightly modified by adding NRWA superscript in each term in Eq. (8), and also adding a term which is a consequence of the time derivative: NRWA expansion affects Eq. (8) as follows: • H NRWA is a 3Q × 3Q Hamiltonian obtained by applying Eq. (13) on each block where Q = 2M + 1. Depending on the problem, some blocks may have only some of the frequency terms and this can simplify the problem of interest.
is enlarged 3 times and represents 3Q × 3Q matrix filled with identity matrices only in positions where the corresponding Hamiltonian has non-zero blocks. For algebraic convenience, we could also define I N NRWA U Q as block matrix with 9Q 2 submatrix identity blocks of N × N size, which has direct algebraic formulation from a unity matrix U as is banded, and appropriate numerical algorithms can be applied, in contrast to the system in Eq. (9) whose sparse properties are generally broken by the dissipator. usually follows the rule in Eq. (13) and needs to comply with the order of storage of the unknown ρ ′′ NRWA . As discussed above, the dissipator of the system depends on the system under study, and therefore its linearization and incorporation into the system cannot be defined generally.
• Frequency terms populate the main diagonal of the system, and this is formulated by Ω NRWA . These terms originate from the time derivative in Eq. (14). The algebraic formulation is Note that the generalization presented so far is only valid for one frequency in Eq. (10). The expansion rule in Eq. (13) can clearly be generalized further if we wanted to include integer multiples of frequency ω. Generally, if we include F multiples of ω, the expansion rule in Eq. (13) will then be (2F + 1)-diagonal (2F + 1) × (2F + 1) matrix.
The difficulty for further generalization lies in frequency interactions. For example, consider a case where we have two frequency components in Eq. (10), ω 1 and ω 2 . The potential terms are ω 1 , ω 2 , ω 1 − ω 2 , ω 1 + ω 2 , 2ω 1 , 2ω 2 which in NRWA expansion would expand the system up to 12 times, however, depending on the values of ω 1 and ω 2 and the state energies, some terms may be neglected due to the underlying physics. A good example can be found in Ref. 26, where three different frequencies were used for the study of nonlinear effects in QCLs.
In general, any linear system of commutators can be formulated through Khatri-Rao product, and different expansion rules would be needed, depending on the problem of interest.

E. Boundary conditions
Boundary conditions represent the crucial part of the system formulation. In this work we focused on periodic systems described by the banded partitioned Hamiltonian, but the effects of device boundaries will generally appear at all four corners of the Hamiltonian and therefore break its banded structure. A vast literature is available for boundary conditions of banded matrices 40,41 and one of the main advantages of the formulation in this work is that one can include the effects of device contacts, similarly to NEGF methods. 36 Expectation value of any operator Ô in DM formalism can be found as Tr(Ôρ). For infinite periodic systems Q → ∞, and the common approach is to find the limit value of Tr(Ôρ) as in Ref. 33, however the result will strongly depend on the chosen boundary conditions, which depend on the problem considered.

III. QUANTUM CASCADE LASER
The derivation in the previous section is general for any periodic system and the differences in implementation of the model lie in the specific forms of the coupling off-diagonal Hamiltonian blocks and the dissipator. In Refs. 24, 26, and 33 we presented one of the simplest nearest neighbour DM implementation that is, in essence, a generalization of approach used in Ref. 22. The simplicity of the model comes from the fact that blocks H 1 and H− 1 which describe the interaction between adjacent periods, only have the dc term which consists of Rabi coupling strengths given in Ref. 42. Note that this approximation is ambiguous for QCL structures since formulation in Ref. 42 considers tunnelling coupling rate of a two well system. This may be circumvented by using generalized scattering approach 32,34 which also offers generalization of the dissipator beyond the Fermigolden rule. We will focus on the simpler model and formulate it with Khatri-Rao notation introduced in Eq. (8) and also consider the dynamic effects by coupling the model with Maxwell wave equation.
Consider a 3 × 3 tridiagonal Hamiltonian and the corresponding tridiagonal density matrix and the dissipator: where H 0 is an N × N matrix involving its dc and two ac terms in the form H 0 = H dc + H ac e iwt + H ac e −iwt . H dc is the tight-binding Hamiltonian of one QCL period, which has the main diagonal filled with bound state energies. Matrix Υ describes the effect of the applied bias K and it will slightly influence Eq. (14). The dissipator blocks can all be expressed in the form ρ j τ , and this was reviewed in Ref. 33. Interestingly the linearization of these blocks satisfies D ′ 1 = D ′ −1 , more details are given in the Appendix. H ac = eZA 1 (t), where Z is the dipole matrix that corresponds to the single period basis wavefunctions and A 1 (t) is the optical electric field at radiation frequency ω.
In Ref. 33 we performed steady-state analysis by taking A 1 = const and ran a minimization algorithm that varied A 1 until the gain of the device was clamped to the loss due to the saturation effect. In this work, we couple the DM model with Maxwell wave equation in order to explore the dynamics of the system. Wave equation can be implemented through travelling wave 29 approach or standing wave approach which correspond to the forms in Eq. (10), respectively. In this work we will couple the DM model with the time-dependent-only wave equation which is derived when Maxwell's equation is approximated with standing waves and we will assume that only one mode in the structure is lasing. Additionally, we will apply slow-varying envelope approximation which is common in laser considerations 43 and which turns the wave equation into the first order differential equation.
In Eq. (15) we have introduced notation different from that in Ref. 33, more similar to that in Ref. 26 because this notation allows the generalization given by Eqs. (2,3). We have also added the effect of external bias. Our main assumption is that the applied bias is weak and that the wavefunction basis can be taken on one period ARTICLE scitation.org/journal/adv and then translated to the neighbouring periods, while the effect of the bias K over the period length L would only be noticeable in the Hamiltonian blocks as the energy shift of ±eKL. Note that this approximation allows the density matrix in Eq. (15) to have equal diagonal blocks, even though its Hamiltonian broke that symmetry. This will have a significant influence on the derivation of the output parameters.
Generally, the addition of the effect of bias influences Eqs. (2,3) by the addition of the term keKLH 0 ρ k , however substitution in Eq. (15) does not change significantly; the term eKLΥ⊠I N U Q should be added in Eq. (8), its NRWA expansion is eKLΥ → eKLΥNRWA = eKLΥ ⊗ I 3N ×3N which needs to be added in Eq. (14), as is further discussed in the Appendix C. The system of equations that Liouville's equation yields is: where γ = c n αL represents the total loss αL multiplied by the group velocity, n is the refractive index and P is the expectation value of the polarization operator. The system under study is represented by an infinite periodic matrix and the form in Eq. (16) focuses on minimal number of equations that need to be solved, however the expected values of the operators need to be derived through usage of infinitely sized banded matrices, as discussed in Ref. 33 L Tr(Zρ + 0 ) where n 2D is the sheet electron density and Z is the dipole matrix that corresponds to the central period basis.
The linearization of the DM part of Eq. (16) is given by Eq. (14) (with addition of +eKLΥ NRWA ρ NRWA ′′ term) and numerical implementation is straightforward after the Khatri-Rao formulation is introduced. Full matrix form of all terms in Eq. (14) is given in the Appendix C, together with discussion on the linearization of the dissipator part. Note that the DM part of the system in Eq. (16) has three N × N blocks of the density matrix, and NRWA expands the system three times, which, after the linearization, has 9N 2 differential equations. For the steady-state analysis we can set the time derivatives in Eq. (14) to zero, and the system would be linear and homogeneous. For the study of dynamics, the system is nonlinear, however the system is of the first order and various ordinary differential equations numerical procedures are available in numerical packages. Normalization condition of the density matrix requires that its trace has unity value. This property only affects the blocks that directly correspond to the basis (Trρ dc 0 = 1) and we can replace one of the equations in the 9N 2 system with this condition in order to solve the system most efficiently or introduce a substitute that one of the elements is one minus all the others. Note that it is not arbitrary which equation is replaced. Only an equation that contains diagonal elements of the block ρ 0 can be replaced.
The algorithm, runs self-self-consistent method for solving Schroedinger-Poisson equation under equithermal subband approximation 33 with A 1 = 1Vm −1 (very low optical interaction) which yields the highest possible value of the optical gain. If that value is larger then the total loss, the system then uses the steady-state wavefunction basis obtained for A 1 = 1V m −1 and couples it with the wave equation as in Eq. (16). Figure 2 shows the dynamic behaviour of 2THz QCL structure 33  Physically, Fig. 2 illustrates the essential property of QCL structures. The relaxation oscillations are barely present as the system reaches the steady-state much faster then its round trip time (≈60ps), even at the peak power. This dynamic behavior is independent on the number of states considered, and it can sustain a system with any number of states, in contrast to the common Max-Bloch approaches. In this case, we included complex interactions between twelve states per QCL period. Numerical implementation of the model was made straightforward due to the compact form of the superoperator, Eq. (14), and apart from the tunneling strengths, input for the presented DM model is identical as in RE approach, 33 and the main difference is in the mathematical aspect of the system construction.
Note that further numerical simplification is possible if we separate real and imaginary parts of the equations in Eq. (16) and use Hermitian symmetry of the system. Additionally, it may be possible to determine the Jacobian of Eq. (16), use the steady-state point value for A 1 = const as stationary point and solve the system

IV. THE SECOND NEAREST NEIGHBOUR INTERACTION
The nearest neighbour interaction is valid for QCL devices due to the fact that we expect a rapid decay of wavefunctions inside the injection barriers, which makes the interaction of the first and the third period very weak. Here we consider a simple two-well periodic structure whose wavefunctions extend across three periods, which requires including higher order Hamiltonian and the corresponding density matrix blocks, labelled as H± 2 and ρ± 2 . The matrix form of the Hamiltonian, density matrix and dissipator have similar form as in Eq. (15), but they are now pentadiagonal matrices, where each matrix has additional super-and sub-diagonal containing block H± 2 , ρ± 2 and D± 2 . The system of equations that Liouville's equation yields is: This equation is then linearized as Eq. (14) (with the addition of the term due to the applied electric field). As an example, consider a GaAs/Al 0.1 Ga 0.9 As two well periodic structure that has three states, their wavefunctions extending even into the second neighbour periods (Fig. 3), and therefore solving system in Eq. (17) is necessary. The calculated current density dependence on the applied electric field is shown in Fig. 4.
In Figure 4 we can notice that the second nearest neighbour approximation did improve the initial result. Note however, that numerical complexity of pentadiagonal model requires solving 15N 2 equations when NRWA is employed. This structure is not a QCL device, but rather an exemplary structure where including the second nearest neighbour interaction is necessary. In Table I  the model difference is most notable in Fig. 4. We calculated the overlap integrals of central period wavefunctions across the first neighbouring Fi = ∫F|ψi| 2 dz and the second neighbouring Si = ∫S|ψi| 2 dz period. The third wavefunction is strongly coupled in both periods which can also be detected in Fig. 3, and in this bias range, there is also a strong coupling of all three wavefunctions in the first neighbouring period. Outlining the necessity for the second neighbour approximation is nominally not straightforward. One of the wavefunctions may strongly couple with the second neighbour, however it might not have a significant influence on the transport. The values in Table I give a rough estimate of the effect. We can additionally illustrate optical field effects by investigating absorption saturation, displayed in Fig. 5.   FIG. 4. Current density versus applied external electric field for different nearest neighbour DM models. The nearest neighbour approximation is shown in red and results in tridiagonal DM problem, while the second nearest neighbour approximation is shown in blue and results in pentadiagonal DM problem.

ARTICLE
scitation.org/journal/adv  6. Current density versus applied electric field obtained by a tridiagonal density matrix 33 and pentadiagonal density matrix for QCL structure. Inset shows optical power dependence on current.
In Figure 5 we can observe how absorption coefficient gets saturated with the increase of optical electric field. Note that here we are not coupling the DM model with Maxwell equation as in the previous considerations, but rather directly changing optical electric field value within the Hamiltonian. Similarly to Fig. 4 there is a difference in results depending on number of neighbours considered in the DM model, shown in the inset.
Implementing the second nearest neighbour approximation for QCL device uses identical expressions as Eq. (17) and comparison between the first and the second nearest neighbour approximations for QCL structure considered in Ref. 33 is presented in Fig. 6. Figure 6 shows nearly identical results as 33 and confirms the expectation that expanding the model does not affect the results significantly. All other results (for material gain and optical power) presented in Ref. 33 are identical as well.

V. CONCLUSION
We presented a generalization of the Liouvillian for periodic quantum systems that can be represented by (2M + 1)-diagonal partitioned Hamiltonian and the corresponding density matrix. We showed that Khatri-Rao product for the system of commutator equations can be applied directly to the initial form of the Hamiltonian, and that it can be written with a similar algebraic notation as the commonly known Liouvillian of non-partitioned Hamiltonian.
In Appendix A we give further generalization for partitioned approaches which can be applied to a variety of quantum systems. Our formulation allows intuitive physical interpretation of the equations in the system and it can extend the DM model to include the effects of the device contacts, which should make this model competitive with more complex models, such as NEGF approach. We illustrated how NRWA affects the periodic system and discussed key approximations that are commonly introduced in the DM formalism. We applied our model to an exemplary periodic structure and QCL structure considered in our previous work, where we illustrated that number of partitions in the density matrix and the Hamiltonian depends on localization and span of wavefunctions over neighbouring periods. We showed that, when applied to a QCL structure, the superoperator given in Eq. (14) requires nearly identical input and algorithm approach as rate equation models, and the only difference lies in algebra needed for model construction. Note that the sparse-banded nature of the equations that this model yields also enables efficient numerical optimization.
We analyzed the absorption saturation with optical field on periodic exemplary structure and presented dynamical response of the bound to continuum THz QCL structure with twelve states considered. The model can be further developed to investigate dynamical behavior of QCLs.

ACKNOWLEDGMENTS
This work was supported by COST action BM1205. Numerical part of this work was undertaken on ARC2, part of the High Performance Computing facilities at the University of Leeds, UK.

APPENDIX A: ALGEBRAIC SIMPLIFICATION FOR PERIODIC PARTITIONED SYSTEMS
Kronecker tensor product is a well known type of product in algebra 44 and for matrices A and B, A ⊗ B creates a block matrix where a block in ijth place is given as aijB. Kronecker product with the identity matrix Im×m has a very intuitive form, when performed from the left as Im×m ⊗ A, the resulting matrix is block diagonal with A blocks on the main diagonal; when performed from the right as A × I the resulting matrix takes every element in A and multiplies it by identity matrix which visually looks like every element in A was diagonally stretched. One of the great advantages of the Kronecker product is that it enables linearization of the system. The system in the form AXB = Y can be linearized as (A ⊗ B T )X ′ = Y ′ for row-wise vectorisation of X and Y.
Khatri-Rao product 45-47 is defined as a "dot" product of partitioned matrices, where the "dot" is the Kronecker product, however alternative definitions exist in the literature 48 and we will focus on the definition as in Refs. 49 and 50, which requires A and B to be partitioned in the similar form. The Khatri Rao product would then behave as originally stated. This product is useful, as we have shown in this work, to describe linear systems given by Kronecker products. We made significant algebraic simplifications which allow efficient numerical implementation along with the intuitive physical interpretation of equations given in the block form.
Let us linearize each commutator in Eq. (6) as Li = (Hi ⊗ I − I ⊗ H T i )ρ ′ i and separate the terms that multiply the unknowns from the left and from the right into two separate groups as: The first term on the right in Eq. (A1) results in a block matrix that has the same form as the Hamiltonian in Eq. (4), where each block goes through Kronecker ⊗ product with the identity matrix I. Similarly, the second term results in a matrix that has a form of block transpose (or "dot" transpose) of the Hamiltonian, where blocks go through Kronecker tensor product with I from the left. It is clear that we can extract the Hamiltonian and its block transpose in Eq. (A1) and let them undergo "dot" product like operation with a (2M + 1) diagonal matrix filled with identity matrices, where "dot" would represent Kronecker product. This type of product is known as Khatri-Rao product ⊠. Let us therefore define a matrix I N U Q as: I I I I I  I I I I I  I I I I I  I I I I I  I I I I where I is N × N identity matrix and UQ×Q is Q × Q unity matrix. The I N U Q matrix is then a Q × Q block matrix filled with identity matrices. Note that minimal value for Q is (2M + 1), but if the system is strongly coupled, Eq. (3) cannot be discarded and the density matrix would need more blocks. It is clear that Eq. (7) follows directly from Eq. (A1) when Khatri-Rao operation is used. Note that the formulation in Eq. (7) can be performed with I N U Q matrix defined as (2M + 1) block diagonal matrix (filled with identity sub-matrices), however the formulation in Eq. (A2) has direct algebraic form.

APPENDIX B: GENERALIZATION OF THE BLOCK FORMULATION
Khatri-Rao notation can be used for constructing a superoperator for general forms of the Hamiltonian and the density matrix which are partitioned. Without loss of generality, we can consider the 2 × 2 case: Note that the symmetry in Eq. (B2) resembles the Kronecker product with the identity matrix and Liouvillian of the commutator L = H ⊗ I − I ⊗ H T . Equation (B2) indeed folds into Liouvillian if all Hij and ρij were scalars, however for submatrix blocks Eq. (B2) is a slightly generalized formulation of the Kronecker product. In Eq. (B2) we have four equations with the terms Hijρij and ρijHij and all these terms can be linearized by a well known identity Hijρij → (Hij ⊗ I)ρ ′ ij , ρijHij → (I ⊗ H T ij )ρ ′ ij , and by using Khatri Rao product Eq. (B2) can be linearized as: This linearization can now be expressed in compact algebraic manner. The first matrix in Eq. (B3) can be obtained from the initial Hamiltonian as (U 2×2 ⊗ I 2×2 ) ⊠ H, because if we had 2 × 2 matrix filled with 2 × 2 identity matrices, the Khatri-Rao product with H would create the first term, and similarly the last matrix needs to repeat transpose of H two times and this can be done as I 2×2 ⊗ H. Without loss of generality, the third and the fourth term could be replaced by one matrix filled with identity matrices, which has trivial algebraic form (IU 4 = U 4×4 ⊗ I), however, if it is of interest to the reader, the third term can be formulated as (U 2×2 ⊗ I 2×2 ) ⊠ (U 2×2 ⊗ I) and the forth as I 2×2 ⊗ (U 2×2 ⊗ I). Note that there are various algebraic formulations of these matrices, and our formulation is not unique. It is clear that generally, for fully partitioned Hamiltonian H and the corresponding density matrix ρ of size QN × QN, whose partitions have the size N × N, the commutator would be given as: where ρ ′′ is first unpacked row-wise in the block form, and then each block is unpacked row-wise, which at the end forms one column vector and Liouville equation becomes linear. This formulation is general, for a variety of problems and in this work we focused on a special case with high symmetry. Note that we focused on square matrices with square blocks, however this is not a strict requirement for Eq. (B3), and further generalization with blocks of different size is possible. Note that C in Eq. (B4) describes the same system as Liouvillian L = HQN×QN ⊗ IQN×QN − IQN×QN ⊗ H T QN×QN would, and the main difference is in the order in which the equations are written. Liouvillian unpacks the unknown density matrix row by row, while C unpacks the blocks first, and then the blocks themselves. Formulation of superoperator in Eq. (B4) offers an intuitive physical interpretation, because in many cases partitions of the Hamiltonian have physical interpretation as well. An additional advantage of Eq. (B4) is that it is also a sparse matrix, as is Liouvillian, and for a periodic system it remains sparse when the dissipator is added to the Liouville equation, which opens the possibility for numerical optimizations. If H and ρ are given as in Eqs. (4) and (5) their commutator yields a multidiagonal matrix with high symmetry as well, and the formulation in Eq. (B4) is not necessary because a large portion of the equations can be discarded, and algebraic generalization needs to be derived separately, as we presented through Eqs. (4)- (14).

APPENDIX C: PARTICULAR TERMS FOR THE QCL MODEL
In Ref. 33 we discussed the dissipator which is present in the most common DM implementation in Eq. (9). Adding the dissipator to the Liouville equation corresponds to modelling interactions in the system by perturbation theory. In QCL, modelling various transport mechanisms is usually performed by introducing relaxation times which are obtained from Fermi golden rule. Along with the simplicity of these terms, Fermi golden rule does not break the positivity of the density matrix and it is always a safe method of modelling the transport. These terms, unfortunately, do not have a general algebraic formulation as the commutator, however, for the QCL, we can represent them by a banded matrix in a similar form as the density matrix. Let us first clarify the physical interpretation of these terms. In a laser system, we expect that transport occurs between various energy levels, and this causes relaxation of the state populations. State populations are described only by the diagonal elements in the overall density matrix, and this relaxation is modelled as which, for the QCL and the banded density matrix we introduced, influences the diagonal elements. These occur only with the blocks ρ 0 and relaxation of the state populations is: 33 τ ij is the total decay rate out of ith level, commonly known as state lifetime. Note that D relax 0 only depends on diagonal elements in ρ 0 . The additional effect that needs to be ARTICLE scitation.org/journal/adv included is the dephasing. Dephasing can be interpreted from simple dipole example, 43 if a wavefunction is given as a mixture of two states described by their wavefunctions ψi = ai(t)e −iE i t/ ̵ h ψi(r) and ψj = aj(t)e −iE j t/ ̵ h ψj(r), the probability will have three terms, individual occupancies = |ai(t)ψi(r)| 2 , |aj(t)ψj| 2 and the third term which will have sinusoidal component at frequency equal to the energy spacing between states i and j. This term resembles the dipole oscillation in the classical mechanics, and directly corresponds to the off-diagonal elements in the density matrix. This term decays as ai(t)aj(t)ψi(r)ψj(r), and this decay will be proportional to the decay of |ai(t)aj(t)|, however these terms describe level occupancy and they can have a phase ai(t) = |ai(t)|e iϕ i which creates additional term e i(ϕ j −ϕ i ) . This phase can be randomized by dephasing processes that do not change occupancies (amplitudes of ai, aj) which is referred as the pure dephasing. Note that the dephasing affects every off-diagonal element in the density matrix, regardless of the block it is located in, and it can be described as: The first two terms in τ ∥ describe the decay of state occupancies, while the next three terms represent an approximation for pure dephasing. 33 For the notation in Eq. (9), γmn = 1 where τm can be easily compared with Eq. (C2). Equations (C1) and (C2) now target all the elements in the overall density matrix, however, we want to partition it in the form similar to the way density matrix is partitioned. We can first formulate the dissipator blocks as D k = ρ k τ k , k = −M, . . . , M where 1 τ k is a tensor. Notice the symmetry of Eq. (C2): it is irrelevant whether states i and j are in the same period or not, if the basis was chosen on one period. This means that the dissipator blocks on off-diagonal of the block formulation of the dissipator have entirely equivalent tensor forms 1 τ ∥ = 1 τ k , k ≠ 0 and the central block will be missing the main diagonal of this tensor 1  33, and then generally the dissipator block form follows the symmetry of the block form of the density matrix, and transport can then be symbolically described by a tensor matrix given in Eq. (5) in Ref. 33.
We are not particularly interested in tensor forms of the dissipator, but rather in their linearization that needs to be inserted in Eq. (14). We only need to make three matrices which we will conveniently name τ −1 , τ −1 ∥ and τ −1 ′′ ∥ , and it applies τ −1 ′′ ∥ = τ −1 ∥ − diag(τ −1 ∥ ) (the main diagonal is set to 0). The first matrix is an N × N matrix filled with scattering rates that come from various scattering mechanisms included in the model and are obtained by Fermi golden rule on the off-diagonal positions, while the terms −τ −1 i on the main diagonal represent the state lifetimes. The second and the third matrix are also N × N matrices with terms given by Eq. (C2), and the diagonal terms in the third matrix are set to zero.
Let us assume that NRWA has not been applied yet, and that we want to write linear form of the dissipator. Each block of D k of the dissipator is N × N matrix, and the linearization requires us to make an N 2 × N 2 matrix that would target the corresponding vectorised form of the density matrix block ρ ′ k . The off-diagonal blocks D k , k ≠ 0, linearize trivially. Matrix 1 τ ∥ is vectorised row-wise and this is then placed on the main diagonal of a N 2 × N 2 matrix which we will refer to as τ −1 lin∥ , and tensor formulation The central block has two tensor components, the dephasing part τ −1 ′′ lin∥ linearizes in the same way as τ ∥ , however some elements will be missing because of the slight difference between these matrices. Linearization that comes from τ −1 is not intuitive, and it was originally derived by mathematical induction in Ref. 33. Visually it looks like stretching of τ −1 matrix, and mathematically this can be formulated as in Ref. 33: matrix τ −1 is written as a sum of N 2 matrices of N × N size, where each matrix in the aforementioned sum is obtained from τ −1 by setting all elements to zero except the ij-th element (any matrix can be written in this form). This sum needs to be written with respect to the element position and to correspond to the row-wise packing. The next step is to place terms of this sum in a N 2 × N 2 matrix row-wise, which we will refer to as τ −1 lin . Alternative definition would be to expand every element in τ −1 by multiplying it with a matrix of N × N size which has unit value in ij-th place and zeros elsewhere. Mathematical formulation can be written as τ −1 lin = τ −1 ⊗ (δijUN×N) (Kronecker delta δij activates ij-th element in the unit matrix and Kronecker tensor product would repeat the operation for every i and j, however this is not fully valid mathematically because δij depends on i and j). Equation (C3) visually illustrates the linearization on the 2 × 2 example. We can now linearize central blocks of the dissipator as D ′ 0 ρ ′ 0 , D ′ 0 = τ −1 ′′ lin∥ + τ −1 lin . Note that this linearization did not consider the NRWA in which every density matrix block is split into three terms, but fortunately the inclusion of NRWA is trivial and, in essence, follows the expansion rule in Eq. (13), the difference being that none of the scattering times are frequency dependent and NRWA will cause D ′ k to simply expand into 3 × 3 diagonal block matrices. This can also be formulated as D NRWA ′ k = I 3×3 ⊗ D ′ k . Full linear form of Eq. (14) that corresponds to the numerical example in section III is:

ARTICLE
scitation.org/journal/adv In Eq. (C4) all blocks in H NRWA are N × N in size, and once this matrix (and its "dot" transpose) gets into Khatri-Rao product with I N U 15 matrix (which is 9 × 9 block matrix filled with IN×N submatrices) we would obtain 9N 2 × 9N 2 matrix where each block would be a corresponding Liouvillian of every block in H NRWA with size N 2 × N 2 . Note that we introduced expansion rule in Eq. (13) to act on Hamiltonian blocks, however, this expansion could have been defined for submatrix Liouvillians as well (similarly to Eq. (12)). Blocks H± 1 do not have ac terms and this simplifies H NRWA . We should underline again that I N U 9 does not need to be completely filled with identity matrices, and that the only real requirement is that it has IN×N blocks at the same places as H NRWA does. As we discussed, D NRWA ′′ , Ω NRWA and Υ NRWA are block diagonal with very high symmetry (Ω NRWA is nearly the identity matrix, while D NRWA ′′ has equal linearization blocks in every position, except the positions that contain ρ dc,ac± 0 ). Note that the overall system is sparse, and that this property can be used in numerical implementation. For the nearest (first) neighbour approximation the system above is reduced to 9N 2 equations by deleting the first and the last three rows and columns in each matrix in Eq. (B3) and setting H± 2 to zero in the remaining equations, which then results in the equivalent model as in Ref. 33.
The derivation of Ω NRWA is simple in the sense that this term is a consequence of the time derivative of i ̵ hρ ′′ NRWA which has 3(2M + 1) submatrices in the form ρ ′ac− k , ρ ′dc k , ρ ′ac+ k , k = −M, . . . , M, and these terms originally have exponentials e ±iwt . After the ARTICLE scitation.org/journal/adv differentiation and multiplication with i ̵ h, the terms with frequency would be present only in front of ac terms, as ∓ ̵ hwρ ′ac± k . Linearization of these terms would simply represent identity matrix I N 2 ×N 2 when this term is moved to the left hand side of Eq. (14) it would alternate as depicted in Eq. (C4). Algebraic formulation of such alternation is Ω NRWA = (I 2M+1×2M+1 ⊗ G 3×3 ) ⊗ I N 2 ×N 2 where G is 3 × 3 diagonal matrix with elements − 1, 0, 1 on the main diagonal.