{'fraction_non_alphanumeric': 0.039344879518072286, 'fraction_numerical': 0.015248493975903615, 'mean_word_length': 4.564594972067039, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Method for the Maxwell-Liouville-von Neumann Equations using Efficient Matrix Exponential Computations

We present a novel method to solve the Maxwell-Liouville-von Neumann (MLN) equations in an accurate and efficient way without invoking the rotating wave approximation (RWA). The method is a combination of two established concepts, namely the operator splitting method as well as the adjoint representation of the Lie algebra SU(N) (or pseudospin representation). The former concept ensures the accuracy of the approach, but is computationally expensive. The latter concept provides an efficient representation of the problem and two optimization possibilities. We have implemented and verified both optimization approaches and demonstrate that substantial speedup can be achieved.


Introduction
The Maxwell-Liouville-von Neumann (MLN) equations describe the interaction of electromagnetic fields with quantum mechanical systems and are an important tool in nonlinear optics, for example to model the dynamics of quantum cascade lasers (QCLs) [1] or quantum well structures [2]. Usually, numerical methods are required to solve the MLN equations. The rotating wave approximation (RWA) is commonly used in order to save computation time but is avoided in this work since it omits certain features of the solution.
Several numerical methods to solve the MLN equations without invoking the RWA have been published in the last decades, starting with the pioneering work by Ziolkowski et al. [3]. This work uses the finite-difference time-domain (FDTD) method for Maxwell's equations as well as the Crank-Nicholson scheme and a predictor-corrector (PC) method for the Liouville-von Neumann equations. In the latter, only two energy levels are considered (in this case, the equations are commonly called Maxwell-Bloch equations).
In the work by Slavcheva et al. [4,5], this approach was extended to cope with more than two energy levels. The methods in [3][4][5] use the adjoint representation of the Lie algebra SU(N) for N energy levels (also called pseudospin representation) and hence eliminate redundant computations and memory requirements. Quantum mechanical operators such as the density matrix are Hermitian, therefore it is sufficient to calculate and store only the half of the offdiagonal elements. Additionally, the trace of the density matrix must equal 1, so the diagonal elements can be expressed using only N − 1 real quantities. Regarding the numerical method, the predictor-corrector method is also computationally efficient.
However, Bidégaray et al. [6] proved that this method may produce unrealistic results when applied to multi-level systems. As an alternative that provides long-term stability in the simulations, an operator splitting approach was suggested [7,8]. The major drawback of the operator splitting technique is the required computational effort to calculate the matrix exponentials at every time step. Different ways to perform this calculation efficiently can be found in literature: The approximation used in [6][7][8], employing the Expokit [9] software package and replacing the FDTD with the pseudo-spectral time-domain (PSTD) method [10], the Magnus expansion via Sylvester's formula [11], the scaling and squaring method as well as a Krylov subspace method [12], and diagonalization of the matrix [13].
Although several research groups used the adjoint representation [3-5, 11, 14] or the operator splitting technique [6-8, 10, 13], there is -to the best of our knowledge -no approach that combines both concepts. In the work at hand we present such a combination and describe two optimization possibilities for the matrix exponential calculations.
In the following, we describe the Maxwell-Liouville-von Neumann equations (Section 2) and the Liouville-von Neumann equation in the adjoint representation (Section 3). In Section 4 we focus on the numerical treatment of the Liouville-von Neumann equation. The resulting method can be used in combination with both FDTD and PSTD, therefore the treatment of Maxwell's equations shall receive less attention in the scope of this work. The verification of the presented methods and their computational performance are discussed in Section 5 using two test cases. Finally, we summarize the performance improvements of the presented methods and give an outlook on future work.

The Maxwell-Liouville-von Neumann Equations
In the work at hand, we consider one-dimensional Maxwell's equations for the electric and magnetic field components E z (x, t) and H y (x, t), where x is the propagation direction, y and z denote the transversal coordinates, and t is time. In Maxwell's equations the conductivity σ, permittivity , and permeability µ of the active region material as well as the polarization term P z (x, t) are taken into account. The second part of the MLN equations is the Liouville-von Neumann equation. It is used to determine the behavior of the quantum mechanical multilevel systems along the propagation direction (described by the density operatorρ (x, t), which can be written as N × N matrix for N discrete energy levels). The equation reads where the right hand side consists of the relaxation superoperator G (ρ) (that may include e.g. scattering processes) and the Liouvillian L (ρ). Here, denotes the reduced Planck constant. The HamiltonianĤ (x, t) =Ĥ 0 +Ĥ I (x, t) consists of a time-independent partĤ 0 and a time-dependent interaction partĤ I (x, t) = −μE z (x, t), whereμ is the dipole moment operator. It is practical to define a corresponding pair of Liouvillian superoperators L 0 (ρ) and L I (ρ). We assume that only the electric field component E z (x, t) is relevant for the interaction process. Finally, the polarization term P z (x, t) in Eq. (1) (or its derivative, respectively) is calculated as where N a is the density of quantum mechanical particles in the system. Since one can write Eq. (3) as where L 0 (ρ) represents the time-independent part of the Liouvillian.

The Liouville-von Neumann Equation in Adjoint Representation
The adjoint representation (also called pseudospin representation) is one of the most efficient ways to describe the density matrix, since only the non-redundant N 2 − 1 real elements are considered [15]. It is practical to write them as a vector ì d. The density matrixρ can be composed asρ whereÎ is the N × N identity matrix,ŝ j are generators of the Lie algebra of SU(N), and the vector elements are defined as d j Tr ρŝ j . The generators are traceless Hermitian N × N matrices and fulfill the condition Tr ŝ jŝk = 2δ jk , where δ jk denotes the Kronecker delta. One possible choice [15] for the generatorsŝ = {û 12 , . . . ,v 12 , . . . ,ŵ 1 , . . . , where the indices satisfy 1 ≤ j < k ≤ N and 1 ≤ l ≤ N − 1 and the transition-projection operators are defined ast jk | j k |. For N = 2 and N = 3 these generators produce the Pauli matrices and the Gell-Mann matrices, respectively.
In order to transform the Liouville-von Neumann equation to the adjoint representation, Eq. (6) is inserted into Eq. (2). Subsequently, the Frobenius inner product1 ·,ŝ k F is applied to the result. The left hand side of Eq. (2) then reads Since both superoperators are linear, we can write for the right hand side. As a result, the Liouville-von Neumann equation can be transformed to where the elements of the N × N matrices L and G are respectively, and the equilibrium vector ì d eq has the elements With the separation of the Liouvillian in mind, we can split the matrix L + G into a timeindependent matrix M and a time-dependent part UE z , whose elements are 1The Frobenius inner product is defined as â,b F = Tr âb † . Note that ·,ŝ k F = Tr {·ŝ k }, since the generators are Hermitian. respectively, and write Eq. (10) as which is the Liouville-von Neumann equation in the adjoint representation. Analogously, the derivative of the polarization must be expressed as function of the vector ì d. In order to do so, we transform the dipole moment operator to a vector ì v µ = 1 2 using the generators as basis2 and insert Eqs. (6) and (15) into Eq. (5). The argument of the trace function in Eq. (5) then reads and can be simplified using the definitions in Eqs. (11) and (13) once the trace function is applied. Then, the derivative of the polarization is calculated as

Numerical Treatment
The achievements by Bidégaray et al. [6,7] serve as basis for the numerical treatment of the equations described in the previous sections. The method in [7] uses the finite-difference timedomain (FDTD) method [16] to solve Maxwell's equations and an operator splitting approach for the Liouville-von Neumann equation. Figure 1 provides a schematic overview. As already mentioned, we focus on the efficient numerical treatment of the Liouville-von Neumann equation in this study. Our efforts can be reused in an approach that uses a different method for Maxwell's equations, e.g. the pseudo-spectral time-domain (PSTD) method used in [10]. Since the adjoint representation was not used in [6,7], we revisit the derivation of the operator splitting approach in the following. Subsequently, we describe two methods to compute the required matrix exponentials efficiently.

Operator Splitting Approach
The right hand side in Eq. (14) can be split up in two operators L 1 (time-independent) and L 2 (time-dependent), yielding t x Fig. 1. Schematic of the discretization. The FDTD method uses a staggered grid for the discretization of the electric field E z (marked with crosses) and the magnetic field H y (circles). The density matrixρ (or ì d, respectively) and the resulting polarization P z is evaluated at the same location as the electric field, but shifted half a time step (squares). Subsequently, the differential equation is solved individually for each operator. For operator L 1 , the solution is where ì d in M −1 ì d eq is the inhomogeneous part of the solution. The solution for operator L 2 can be determined as3 The complete differential equation can be solved using the symmetric Strang splitting [17], i.e., by updating ì d first with the solution of operator L 1 for half a time step then with the solution of operator L 2 for a full time step

Verification and Performance Evaluation
We implemented the diagonalization technique as well as the method using the generalized Rodrigues formula based on the mbsolve project [21]. This project has recently served as base for our operator splitting approach using the Eigen library [22] and the Padé approximation for the matrix exponential calculations [18]. The Padé approximation approach serves as reference in terms of verification as well as performance in two test cases, which we describe in the following. Finally, we discuss the performance of each technique.

Two-level system (Maxwell-Bloch case)
The first test case is the self-induced transparency (SIT) simulation by Ziolkowski et al. [3]. In this simulation, a two-level system is described with the Hamiltonian where ω 12 = 4π × 10 14 s −1 is the transition frequency and µ 12 = 10 −29 A s m is the dipole moment between the two levels, and the relaxation superoperator where T 1 = 10 −10 s −1 and T 2 = 10 −10 s −1 are decay terms (upper level lifetime and dephasing time, respectively).6 These components are transformed to adjoint representation as described in Section 3. Additionally, a initial density matrix can be transformed to an initial condition ì d 0 . The simulation was set up with 32 768 spatial grid points, which corresponds to a spatial discretization size of 4.578 nm. By setting the Courant number C = 0.5, we chose a time step size of 7.635 × 10 −18 s. Then, the simulation was executed and the populations of both levels were recorded. A snapshot after 187.5 fs is depicted in Fig. 2. It shows perfect agreement between the results of all methods and the data depicted in [3].

Three-level system
The second test case incorporates the three-level setup presented by Song et al. [23]. The Hamiltonian isĤ 6Note that a different sign convention is used for µ 12 in [3]. However, the convention has no effect on the result if applied consistently. where ω 2 = 2.372 s −1 and ω 3 = 2.417 s −1 are the eigenfrequencies, µ 12 = 1.48 × 10 −29 A s m is the dipole moment between the levels 1 and 2, and γ = √ 2 is the ratio between the dipole moments. Furthermore, we define the relaxation superoperator as where g 1 = ρ 33 + ρ 22 − 2ρ 11 , and T 1 = 10 −10 s is the relaxation time. As initial condition, the density matrixρ was transformed to the adjoint representation.
Using the choice of grid point size and time step size from the two-level test case, the results from [23] could be reproduced. Figure 3 shows the agreement of the different methods for this test case.

Performance
The simulations were executed on an Intel Xeon Processor E7-4870. Although the code is designed to run efficiently in parallel, only one thread was used since we wanted to measure the quality of a single-thread optimization. In order to ensure reproducibility, the measurements were repeated five times. The performance was determined as number of grid point updates per time unit, i.e. P = N x N t /t exec , where N x and N t are the number of spatial and temporal grid points, respectively, and t exec is the measured execution time. This allows the metric to be used to compare the performance of problems with different sizes.
It should be noted that the execution time does not contain the time required to construct or delete data structures (which is negligible for the simulation setups in question). It does, however, contain the time required to store the result data, in particular the time to convert the vector ì d to the density matrixρ.
The performance values of the different methods applied to the two test cases are shown in Fig. 4. Compared to our baseline (the Padé approximation method), the diagonalization approach performs 1.26x better for the two-level test case and 1.55x better for the three-level test case. Using the generalized Rodrigues formula, we were able to achieve a speedup of 3.48x and 7.35x, respectively.

Conclusion
Our method solves the Maxwell-Liouville-von Neumann equations in an efficient and accurate way. It is a novel combination of two established concepts. The adjoint representation is most efficient in terms of eliminating the redundancy inherent to quantum mechanical operators and additionally leads to interesting properties of the resulting description. These properties are exploited in order to provide two efficient implementations of the operator splitting technique, which is accurate and stable, yet computationally expensive.
Both implementations are tested with the help of a two-level and a three-level test case and can reproduce the results found in related literature. Regarding the performance, the diagonalization approach yields a small improvement and the method using the generalized Rodrigues formula provides substantial speedup, at least at test cases with small level count.
We expect that both methods provide even larger performance improvements compared to standard methods (e.g. the Padé approximation) in test cases with many energy levels. Also, the performance ratio between the diagonalization approach and the Rodrigues method may differ. However, the behavior of the methods applied to a many-level setup has to be investigated first.
Finally, the method is able to cope with a broad range of problems. However, there are two limitations, namely the restriction to one-dimensional models as well as the strict definition of the interaction Hamiltonian. In future work we will aim to generalize our method in order to overcome those limitations.

Supplementary Material
See [21] for source code, build instructions, and basic documentation. In this paper we used the development branch riesch2017c, the Eigen library version 3.3.4, and the Intel C++ compiler 17.0.