Numerical time propagation of quantum systems in radiation fields

Atoms, molecules or excitonic quasiparticles, for which excitations are induced by external radiation fields and energy is dissipated through radiative decay, are examples of driven open quantum systems. We explain the use of commutator-free exponential time-propagators for the numerical solution of the associated Schr\"odinger or master equations with a time-dependent Hamilton operator. These time-propagators are based on the Magnus series but avoid the computation of commutators, which makes them suitable for the efficient propagation of systems with a large number of degrees of freedom. We present an optimized fourth order propagator and demonstrate its efficiency in comparison to the direct Runge-Kutta computation. As an illustrative example we consider the parametrically driven dissipative Dicke model, for which we calculate the periodic steady state and the optical emission spectrum.


Introduction
The outcome of many experimental measurements is well described by linear response theory for situations close to thermal equilibrium. Other experiments, predominantly those dealing with small quantum systems in strong external fields, require a full non-equilibrium description. One example is cavity quantum electrodynamics [1], and generally finite quantum systems in radiation fields. While the interaction of a single atom or an atomic ensemble with the quantized cavity field is weak, transitions between atomic levels can be induced with strong, classical laser fields. Through cavity losses and spontaneous emission the energy input from the external pumping dissipates. Atoms in a cavity are open quantum systems far from equilibrium.
Such situations are described either by the Schrödinger equation with a time-dependent Hamilton operator H(t) if we neglect dissipation, or more generally by a master equation with a time-dependent Liouville operator L(t), e.g. one of Lindblad-type which includes dissipation in the Markovian approximation [2]. In addition to single-time expectation values, which provide the basic information from time propagation of the wave function |ψ(t) or density matrix ρ(t), one is interested in many-time correlations functions that yield optical spectra or information about the coherence or statistical properties of the emitted light [3,4].
Since explicit solutions of linear differential equations with variable coefficients do not exist apart from simple situations, the above equations fall into the domain of numerical time-propagation. The topic of the present paper is the application of commutator-free propagators based on the Magnus series [5]. The Magnus series arises in the context of differential equations on Lie groups, where it allows, among many other things, for the systematic construction of high-order approximations to the propagator [6]. Commutator-free exponential time-propagators (CFETs) avoid the use of commutators that appear in the Magnus series. They provide an efficient and accurate algorithm for numerical time propagation [7], which we discussed for the Schrödinger equation in reference [8]. Here, we concentrate instead on master equations for open quantum systems. Although the present application lies outside of the principal Lie group setting, we feel that the numerical results presented here are promising enough to warrant closer inspection. The application to the parametrically driven dissipative Dicke model in section 5 gives an indication of the potential of this approach in non-trivial situations.
The paper is organized as follows. In sections 2 and 3 we discuss the basic numerical problem and its principal solution through the Magnus series. The commutator-free exponential time-propagators are introduced as a more practical solution in section 4, and an optimized 4th-order propagator is given in 4.2. After a demonstration of their usage with the example of a spin in a magnetic field in section 4.3 we turn to a discussion of the parametrically driven dissipative Dicke model in section 5, before we summarize in section 6.

The numerical problem
To describe the basic numerical problem we consider the Schrödinger equation (1). The standard approach obtains the wave function |ψ(t) through time stepping. The middlepoint approximation allows for propagation of |ψ(t) over a short time interval [t, t+δt]. Repeated application of equation (3) gives |ψ(t + T ) starting from |ψ(t) with N = T /δt time steps. As detailed later, straightforward expansion of the exponential shows that the error of one time step with equation (3) is ∝ (δt) 3 , such that the total error ∝ N(δt) 3 = T (δt) 2 for propagation time T scales as (δt) 2 . Conversely, the achieved accuracy scales as (δt) −2 , which we can write symbolically as error ∝ effort −2 .
As a second-order method the middle-point approximation is not efficient and requires small δt even for low accuracy demands. If we ask for a better scheme we should note that the approximation (3) has two independent sources of error. The genuine error in the situation of a time-dependent Hamilton operator arises from the replacement of H(t) by the constant H(t + δt/2), and depends mainly on the rate of change of H(t). In addition, the numerical computation of an operator or matrix exponential exp[A] involves an error determined by the spread of eigenvalues of A. In equation (3) it is roughly proportionally to δt and the norm H(t + δt/2) .
Often the rate of change of H(t), e.g. set by an external field frequency, is smaller than the largest eigenvalues of the Hamilton operator corresponding to highly excited states. Then the total error is dominated by the computation of the exponential in equation (3). Very small time steps δt and correspondingly large effort are required even if H(t) changes slowly.
This observation explains why the use of general algorithms for the solution of ordinary differential equations (ODEs), e.g. the standard 4th-order Runge-Kutta (RK4) procedure [9], cannot be recommended unreservedly for the Schrödinger equation. The problem of such (explicit) ODE solvers is that they provide only a poor approximation of the exponential in equation (3) and are inefficient already for constant H.
For example, the RK4 procedure approximates the exponential by the 5 terms exp[A] ≈ 1 + A + A 2 /2 + A 3 /6 + A 4 /24 of the Taylor series of e x . The problem is that the Taylor series is not a good approximation unless |x| is very small. This effect is shown in panel (a) in figure 1, where it is compared to an approximation using Chebyshev polynomials (of the first kind) [10]. In this example, the five term Chebyshev approximation is 16 times more accurate than the five term Taylor series. For the 4th  (4)) for the calculation of exp[−iHt] with the diagonal 10 × 10 matrix with eigenvalues H nn = n and t = π/5. We compare the 4th-order Runge-Kutta procedure (RK4) with the use of the 4th-order Chebyshev approximation in a time-stepping scheme (Cheb 4 ), and with a single propagation step using N terms of the Taylor series (Tay) or of the Chebyshev approximation (Cheb). In time-stepping the error decays as a power (here ∝ N −1/4 ) of the effort, while full computation of the exponential in a single step achieves much quicker error reduction. order approximation error ∝ effort −4 , this implies that the efficiency is increased by a factor 16 1/4 = 2.
In panel (b) in figure 1 we compare the error-effort relation of the RK4 procedure to that of the Chebyshev approximation for the calculation of exp[−iδtH], where H is a diagonal matrix with entries H nn = n as for a harmonic oscillator. Here, and also in later examples (figures 2 and 3), we give the error between an exact and numerical matrix A e/n as the maximal difference of matrix elements The Chebyshev approximation is clearly superior already for a small δt. It becomes even better with increasing eigenvalue spread or time-step |δt|. This reasoning motivates the replacement of ODE solvers by techniques which take advantage of the linearity of the Schrödinger or master equations. The calculation of a matrix exponential is better accomplished with specialized algorithms, such as splitoperator methods [11], or the Krylov [12,13] or Chebyshev technique [14] in the case of large sparse matrices. They allow for efficient propagation with time-independent Hamilton operators. Equipped with such algorithms it remains to improve on the genuine error ∝ (δt) 2 involved in equation (3) when turning to time-dependent Hamilton operators.

Magnus propagators
The importance of accurate evaluation of exponentials for the Schrödinger equation is related to the fact that the exponential maps the hermitian Hamilton operator H onto the unitary propagator exp[−itH]. Many differential equations involve a Lie algebra (here: of hermitian Hamiltonians) and a Lie group (here: of unitary propagators) in this way. The idea behind "geometric numerical integration" of ODEs [15] is that also an approximate propagator should stay in the respective Lie group.
Let us consider general linear differential equationṡ where x(t) is a vector and the coefficient matrix A(t) is time-dependent. The formal solution of equation (5) is provided by the propagator U(t) that gives for all x(0). For a scalar equationẋ = a(t)x, the propagator is obtained through integration For operators or matrices [A(t 1 ), A(t 2 )] = 0 is possible, such that this expression does not generalize. We can still read this expression as an approximation For small t = δt, this is nothing else than the approximation (3), if the integral over τ is approximated (also with error (δt) 3 ) using the middle-point value A(δt/2). Although equation (7) involves a finite, maybe large, error its exponential form guarantees that the approximate propagator lies in the Lie group. The question is whether we can improve on the (δt) 3 scaling of the error and preserve the exponential form. An affirmative answer is given by the Magnus series [5,6], which gives as an exponential of Lie algebra elements Ω n (t) and provides a systematic scheme for their construction. The first term in (8) is the term known from the scalar case. The non-commutativity of A(t) is accounted for by correction terms Ω n (t), for n ≥ 2. The term Ω n (t) is given by a time-ordered integral of n-fold nested commutators of A(τ i ) and can be obtained through a recursive calculation. The first two terms are and Ω 3 (t) = 1 6 t 0 dτ 1 Explicit expressions for higher-order terms become quickly unwieldy. Importantly, by building terms from commutators of Lie algebra elements A(τ i ), every Ω n (t) stays in the Lie algebra. The Magnus series solves two problems. On the one hand, it preserves the Lie group structure of a differential equation. For the Schrödinger equation, where A(t) = −iH(t), the propagator U(t) is unitary as the exponential of an anti-hermitian matrix. Furthermore, unitarity of U(t) is preserved for any truncation of the Magnus series. On the other hand, since the term Ω n (t) involves an n-fold integration over time, its size scales as (δt) n . Working with a truncated series including terms Ω n for n ≤ N only, the error of the obtained approximate propagator itself scales as (δt) N +1 . We can thus improve systematically on the middle-point approximation (3) by including more terms from the Magnus series.
Unfortunately, the Magnus series does not solve the practical problem of finding an efficient numerical time-propagation algorithm. Computation of the nested commutators and multiple integrals is difficult to implement and consumes computational resources. Fortunately, there is a simpler and more convenient way.

Commutator-free exponential time-propagators
The use of commutator-free exponential time-propagators (CFETs) has been discussed in references [7,8,16]. They are, basically, a reformulation of the Magnus series that avoids integrals and commutators and gives the propagator as a product of exponentials of simple linear combinations of A(t).
The simplest CFET is the middle-point approximation (3) itself. A 4th-order CFET, where the error scales as (δt) 4 , was introduced in [7,16]. It gives the approximate propagator as the product of two exponentials U CFET (δt) = exp δt g 1 A (1) + g 2 A (2) exp δt g 2 A (1) + g 1 A (2) , (11) which involve a linear combination of A(t) specified by the coefficients and uses only the values of A(t) evaluated at two points in [0, δt] given by This expression is the simplest non-trivial CFET. A better, optimized, 4th-order CFET is presented below in equation (22). Higher-order CFETs can be constructed, and the freedom in the choice of coefficients can be exploited for their optimization, i.e. the minimization of the error. The construction of CFETs is rooted in the theory of abstract free Lie algebras that underlies the Magnus series. Its description is beyond the scope of this paper, and we refer the reader to reference [7] and our reference [8] for details. Here, we proceed in the opposite way and give a direct check of the validity of equation (11) that avoids most of the language of free Lie algebras.

Direct validation of the 4th-order CFET
The principle idea is to combine the two exponentials in equation (11) with the Baker- and compare the resulting expression with the original Magnus series. Let us begin with the Taylor series of A(t), in the vicinity of t = 0. For the 4th-order CFET, only terms A 1 , . . . , A 4 have to be considered. We insert the Taylor series in the Magnus series (8), and keep the first three terms to Ω 4 (t). The terms Ω n (t) for n ≥ 5 give contributions of order (δt) 5 and higher. A simple counting of indices shows that only the seven terms can contribute in fourth order. The Magnus series thus gives the propagator Note that the commutator [A 1 , [A 1 , A 2 ]] does not contribute. This expression is the exact reference for comparison. It is easy to see that the middle-point approximation (3) is correct to second order (δt) 2 : The terms δtA 1 +((δt) 2 /2)A 2 are reproduced exactly, but the commutator [A 1 , A 2 ] in the third order term is missing.
For the 4th order CFET from equation (11), it is for k = 1, 2, which inserted gives . We now use the BCH formula to combine the two exponentials. Only the three commutators shown in equation (19) need to be evaluated, the following commutators in the BCH formula contribute terms of order (δt) 5 or higher. We then obtain the expression which allows for direct comparison with the Magnus series in equation (16). We immediately recognize the seven terms and commutators from there. Comparison of the coefficients ξ i , which we evaluated with the BCH formula, to the prefactors in (16) gives the conditions To complete our check we insert x 1 , x 2 from equation (14) and g 1 , g 2 from (12) and find that all the seven conditions are satisfied. Note that condition (21g), and also (21d) and (21f), are redundant. For the construction of a CFET, this process has to be reversed. We start from an ansatz with a product of exponentials, derive the relevant conditions for a CFET of required order, and solve the resulting polynomials equations for possible coefficient values. Several adjustments of the direct calculation done here simplify bookkeeping, and reveal underlying structures which reduce the number of conditions. Still, the construction of higher-order CFETs is involved and not entirely free of brute force computations. We refer the reader to reference [8] to get an impression, where CFETs up to order 8 are presented. The usage of a given CFET, however, is plain and simple.

Optimized 4th order CFET
For the application to dissipative systems we recommend here the use of an optimized 4th-order CFET, which is equation (43) in our reference [8]. Extending the simpler expression (11) it gives the approximate propagator as the product of three exponentials where with and The error of this CFET scales again as (δt) 5 , but the prefactor in front of the error term is considerably smaller than for the CFET (11). The reduction outweighs the increase of effort using three instead of two exponentials.
In contrast to the original Magnus series usage of the CFET (22) is compellingly easy. Only linear combinations of A(t) evaluated at three different points in the interval [0, δt] need to be formed. All commutators and integrations have been removed from the expression. Of course, we assume the existence of an algorithm for the computation of matrix exponentials.
Let us stress the advantage of easy usage with the even simpler formulation that is obtained for an A(t) = B + f (t)C (it is easily generalized to include more terms). In this case equation (22) can be written as with the time steps and the coefficients Through the CFET the full propagation with a time-dependent term f (t)C is replaced (approximately) by piecewise propagation with a constant term f i C. The choice of the coefficients f 1 , f 2 , f 3 according to equation (28) guarantees that the error of this approximation scales as (δt) 5 . This is the signature of geometric integration [15]: Figuratively speaking, instead of moving along a curve in the Lie group we move repeatedly along short straight lines, the direction of which is given by the linear combinations in equation (22) or (26).
Note that the time steps δt 1 , δt 2 are positive such that propagation proceeds in the forward direction. This is important for dissipative systems where a negative δt i would push eigenvalues of L(t) into the right half complex plane, which leads to exponentially growing terms and corresponding numerical instabilities. For this reason we restrict ourselves to 4th-order CFETs here.

Exemplary application of the optimized 4th-order CFET
Let us apply the 4th-order CFET (22) to a simple example and compare with the RK4 procedure which, as we claimed, should be less efficient because it does not properly compute exponentials. We have to stress that the efficiency of CFETs depends on a good algorithm for the computation of matrix exponentials. Otherwise, when the additional computational overhead involved exceeds the savings achieved with a large time-step δt, the simple Runge-Kutta procedure is more efficient.
Good algorithm for the symmetric case (A † = ±A, e.g. for a hermitian Hamiltonian) are the Chebyshev, Krylov and split-operator techniques mentioned before. For the unsymmetric case ([A, A † ] = 0) encountered for dissipative systems, all techniques meet problems which are only partially solved. After suitable modifications of the standard procedure the Chebyshev technique behaves most favourably. The exploration of this point has to be left for a future publication, here we take the virtues of the Chebyshev technique for granted.

First example: Spin in a magnetic field
As an example for a non-dissipative system consider a spin (length j) in a rotating magnetic field, with Hamilton operator The time-evolution of the wave function can be determined exactly after transformation with exp[iωtJ z ] to the rotating frame (see below).
In figure 2 we plot the error-effort relation for the optimized CFET (22) and the RK4 procedure for a typical set of parameters (the behavior for other parameters is identical). The error ε is determined as in equation (4). For the effort N H we count the number of evaluations of matrix-vector multiplications with the Hamilton operator, which is generally the most time-consuming step. The matrix exponentials needed for the CFET are calculated with the Chebyshev technique to machine precision.
We see that the RK4 procedure requires lesser effort for low accuracy only, but use of the CFET becomes quickly advantageous as the spin length or propagation time increases. For j = 10 and t = 100 in panel (b), the CFET is more efficient for error goals less than 1%, with an a efficiency gain of a factor 2 . . . 4.

Second example: Driven dissipative two-level system
We keep the spin in the rotating magnetic field as an example and include dissipation. With dissipation, its time evolution is described by a master equation (2) for the spin density matrix ρ(t).
In the Lindblad formalism, the Liouville operator L(t) = L H (t) + L D is the sum of two or more terms with different meaning [2]. The first term contains the Hamilton operator H(t) and will be time-dependent. The second and further terms have the form They introduce eigenvalues of L(t) with finite (negative) real part and thus describe dissipation. Within the Lindblad formalism the form of D[A] guarantees that the structural properties of the density matrix -hermiticity, normalization, positive semidefiniteness -are strictly preserved. Note that the numerical time propagation scheme does not depend on the precise form of L(t), as long as the master equation remains linear and local in time.
For the driven spin, dissipation is included through the Lindblad term and the full Liouville operator is with the dissipation rate γ > 0. For j = 1/2, the exact solution of this problem is possible with a transformatioñ ρ(t) = exp[iωtJ z ]ρ(t) exp[−iωtJ z ] to the rotating frame, which gives a time-independent Hamilton operatorH = 2(∆ − ω)J z + 2V J x . Note that the transformation leaves J z invariant.
The stationary state in the rotating frame, corresponding to the eigenvalue zero of the transformed Liouville operatorL for γ > 0, is In particular, J z (t) converges for t → ∞, with constant value In panel (a) in figure 3 we plot J z (t) starting from the initial state with J z (0) = +1/2. Transient oscillations decay as a result of finite dissipation γ > 0. The value of J z (t) in the quasi-equilibrium state depends on the field frequency ω, which we increase slowly during time-propagation. J z (t) follows closely the value J z ∞ from equation (36), with a short delay, and we identify the resonance at ω(t) = ∆.
These expressions allow for comparison with numerical results.
For the numerical solution of this problem, we do not transform the problem but keep the time-dependence of H(t) explicitly. In panel (b) in figure 3 we plot the erroreffort relation as in figure 2. The relation between the CFET and the RK4 procedure is similar to the non-dissipative case, and we recognize the factor 1/2 of error reduction for j = 1/2. The advantage of the CFET is however not quite as distinct as for the dissipation-free case in figure 2, panel (a).

The parametrically driven dissipative Dicke model
The previous examples serve as a benchmark for the CFET approach. We can now demonstrate its usefulness for a less academic case, the parametrically driven dissipative Dicke model. Optical properties of the Rabi case j = 1/2, corresponding to a single qubit, have been explored in reference [17], and quantum phase transitions in the parametrically driven Dicke model without dissipation are studied in [18].
The Hamilton operator of the Dicke model [19], describes an ensemble of two-level atoms (transition energy ∆) as a pseudo-spin of length j, which couples to the cavity field (frequency Ω). We assume a time-dependent interaction constant  Figure 4.
Comparison of the recommended (equation (22)) and simpler (equation (11) with loss rate κ ≥ 0. In contrast to the standard quantum optical treatment in rotating wave approximation, it is not possible to eliminate the explicit time-dependence of H(t) through a transformation to the rotating frame. Because the rotating wave approximation is not applicable, physical properties of the parametrically driven dissipative Dicke model have to be extracted from time propagation of the density matrix of the joint atom-photon system. The number of entries of the density matrix, which grows as ∝ (2j + 1) 2 , becomes large already for moderate pseudo-spin length j. In addition, highly excited states contribute to the dynamics if j grows. Therefore, as we discussed in section 2, the advantage of CFETs over general ODE solvers such as the RK4 procedure will be pronounced for this more complex example.
In figure 4 we compare the recommended CFET from equation (22) to the RK4 procedure and the naive middle-point approximation (3). We see that already for j = 5 we can easily reduce the numerical effort by a factor eight if we use CFETs. The reduction is achieved independently of the intended error ǫ. The middle-point approximation, which is only a second order scheme, is not able to compete even with the RK4 procedure. This clearly supports our recommendation for the CFET (22), and extends the positive results from [8] to dissipative systems. Note that the CFET (22) is better (by 50%) than the simpler CFET (11) although it requires computation of three instead of two matrix exponentials.
Note again that the advantage of CFETs stems from the fact that the propagator is approximated as an product of exponentials of the Hamilton or Liouville operator. We know that this strategy is favourable for non-dissipative systems because it respects the unitary geometry of the Schrödinger equation [6,15], but apparently it works well also for density matrix propagation in driven dissipative systems. Conversely, CFETs rely on good computation of matrix exponentials. As mentioned at the beginning of section 4.3 we use the Chebyshev technique here because it proved to be more efficient and reliable than a Krylov computation for non-hermitian matrices.

Steady state resonances
In panel (a) in figure 5 we show the time evolution of the initial state ρ(0) = |ψ ψ|, where |ψ = |−j/2 ⊗ |vac is the state with no atomic or field excitations. We plot the population inversion for a linear variation of the modulation frequency ω p from 1.5 to 2.5. Transient oscillations are observed for t 500, before two resonances evolve at a time corresponding to ω p ≈ 2 ± 0.12. Beyond the resonances, I z (t) decays again to a small value with weak oscillations. The energy level diagram in panel (b) in figure 5 explains the appearance of resonances in panel (a) through transitions between the states of the Jaynes-Cummings ladder [20].
The eigenstates of the zero coupling Hamiltonian are the J z , a † a eigenstates |m, n , with m + j atomic excitations and n photons. For Ω = ∆ the states |m + k, n − k , for several integer k, are degenerate with energy (m + n)Ω. At weak coupling λ 0 ≪ Ω, ∆, degenerate states are split ∝ λ 0 by the atom-field coupling. Counting energies relative to the energy of the lowest state |−j, 0 , the energies of the two singly excited states (|−j + 1, 0 ± |−j, 1 )/ √ 2 are given by E 1,± = Ω ± √ 2jλ 0 .
A weak modulation of λ(t) introduces transitions between |−j, 0 and the doubly excited states (vertical arrows in panel (b)). Note that the splitting of energy levels in the diagram is determined by the co-rotating terms J − a † , J + a in the coupling term J x (a + a † ), which preserve the number of excitations in the sense of the rotating wave approximation, but the transitions arise from the counter-rotating terms J + a † , J − a and change the number of excitations by two. In the Rabi case j = 1/2, the two doubly excited states (|1/2, 1 ± |−1/2, 2 )/ √ 2 have energy E 2,± = 2Ω ± √ 2λ 0 . Resonances are expected at these energies [17].
In the Dicke case j > 1/2, the splitting of the three degenerate states |−j + 2, 0 , |−j + 1, 1 , |−j, 2 has to be calculated. This gives the energies E 2,± = 2Ω ± √ 8j − 2λ 0 for the odd/even parity combination, reproducing the j = 1/2 result. The third state with unshifted energy E 2,0 = 2Ω is a linear combination of |−j + 2, 0 , |−j, 2 only. It does not couple to |−j, 0 through the counter-rotating terms, and the transition to this state is forbidden in leading order of perturbation theory. Therefore, we also expect only two resonances in the Dicke case. For the parameters from figure 5, they occur at E 2,± = 2 ± 0.02 √ 38 ≈ 2 ± 0.1233. To identify the resonances numerically we propagate the system with fixed ω p until the periodic steady state is reached. Then we calculate the quantities averaged over one modulation period 2π/ω p . Here, is the number of cavity bosons. In figure 6 the quantities I z , N b are shown as a function of ω p . We recognize the two resonances ω p ≈ E 2,± , which are broadened due to the cavity losses ∝ κ. For the calculation we used the optimized 4th-order CFET (22) together with a Chebyshev computation of the exponential.

Emission spectrum
To study the optical properties of this system we compute the cavity emission spectrum S(ω). It is obtained as the Fourier transform of the correlation function which we calculate with the quantum regression theorem [3] through time propagation of the operator aρ(t) (for τ ≥ 0). The correlation function involves the average over one modulation period [t, t + 2π/ω p ] for large t, i.e. in the periodic steady state. We include in figure 6 the total emission which is given by the integral over positive ω in accordance with the fact that emission of a (real) photon can only decrease the energy. We note the normalization We see that S tot ≈ N b close to resonance, when emission is strong. Away from resonance S tot drops below N b , since N b counts also bound photons that do not contribute to emission. However, S tot remains finite as a consequence of the Markovian approximation used here for the dissipative term [17].    interpretation of the emission spectrum is again possible using the energy level diagram from figure 5.
For ω p ≈ E 2,+ in figure 7, the higher of the two doubly excited states is populated. Since the operator a changes the number of excitations by one, an atom in this state can decay to the lowest state only through the intermediate singly excited states. Four transitions corresponding to the four vertical arrows in panel (a) can be identified in this situation. For the Rabi case j = 1/2, they are in order of increasing energy The numerical values correspond to the parameters from figure 7, and the transitions are marked correspondingly in panels (a), (b). In transitions 2 and 4 the doubly excited state decays. These transitions are shifted by ω p − E 2,+ away from resonance (transitions 2 ′ , 4 ′ in panels (c) and (d)). In transitions 1 and 3 the singly excited states decay, and their energy does not depend on the modulation frequency. Note that transition 4 (dashed arrow) is parity forbidden in lowest order perturbation theory and gives only a weak signal.
The opposite case ω p ≈ E 2,− in figure 8 has an analogous interpretation that follows from the energy diagram in panel (a). Now the lower of the doubly excited states is populated, and the order of transitions is reversed, with transition 1 as the parity forbidden transition.
For the Dicke case j > 1/2 in figures 9, 10 we expect the same qualitative behavior as for j = 1/2 since the additional state with energy E 2,0 does not participate in the transitions. For the higher resonance ω p ≈ E 2,+ in figure 9, the situation corresponding The numerical values correspond to the parameters from figure 9. As a difference to the Rabi case j = 1/2 we note that ω 2 ≈ ω 3 , such that the two transitions cannot be distinguished in panel (b) because of the finite linewidth acquired through cavity losses. If we change ω p transition 2 is shifted but transition 3 remains fixed. This allows for the separate identification of both peaks in panels (c), (d). In the same way we can understand the opposite case ω p ≈ E 2,− in figure 10, which is similar to figure 8.
For stronger coupling, additional peaks appear in the emission spectrum. The peak position and height is shown in Fig. 11 for the situation corresponding to figure 7. Note that we adjust the modulation frequency ω p = E 2,+ for different λ 0 to remain close to the higher resonance. For λ 0 0.1 a fifth transition peak "(5)" becomes visible, and the previous interpretation of S(ω) via the Jaynes-Cummings ladder breaks down. At least such situations require numerical time propagation because a simple perturbative interpretation is no longer possible.

Summary and Outlook
Numerical time propagation allows for the theoretical description of experimentally relevant non-equilibrium situations beyond the linear response regime. Such situations arise in particular if small systems such as atoms or molecules are manipulated by strong radiation fields. Important directions of research include the optical properties of ensembles showing collective behavior, such as polariton or exciton condensates [21,22]. A characteristic optical signature, as of the spatial shape and energy distribution of the optical emission [23,24] or the coherence properties of the emitted light, can provide the proof of existence for a condensate. A fundamental theory of optical properties of collective phases of (quasi-) particles with finite lifetime requires the description of the driven open quantum system that is realized, e.g., by the excitonic condensate in a semiconductor [25]. Different but on the fundamental technical level related questions arise in the field of non-equilibrium transport problems [26].
Increasing complexity of the physical situation under study coincides with an increase of the computational effort, which underlines the need for powerful numerical algorithms. We argued here in favor of commutator-free exponential time-propagators as a convenient alternative to the original Magnus series. Shifting the focus of previous studies, where CFETs were shown to be well suited for the time propagation of driven systems without dissipation, we here applied CFETs to open quantum systems. Using the parametrically driven Dicke model as an illustrative example we calculated the optical emission spectrum with this technique.
Conceptually, CFETs are recipes for the reduction of the original problemthe solution of the Schrödinger or master equations with a time-dependent Hamilton operator -to the computation of matrix exponentials. Because they replace the naive approximation of the time propagator by a single exponential per time step, as it is encoded in the second order middle-point approximation (3), with a more sophisticated combination of exponentials, higher-order CFETs significantly reduce the numerical effort. The particular appeal of CFETs is that they can be combined with any technique for the computation of matrix exponentials, for example the powerful Krylov (Arnoldi) or Chebyshev techniques in the context of large sparse matrix computations. CFETs do not compete with such techniques, but instead serve the complementary purpose of achieving a favorable error-effort scaling also for equations with a time-dependent H(t). Therefore, CFETs should be of interest to anyone presently using Krylov (Arnoldi) or Chebyshev techniques in studies of driven quantum systems: It is straightforward and simple enough to add the CFET computational scheme from equations (22), (26) to an existing program or implementation [27,28].
In addition to the principal research directions mentioned above many open problems arise within the more restricted context of the present work. A notorious problem is the efficient evaluation of the matrix exponential for non-symmetric large sparse matrices, which is essentially a problem of polynomial approximation in the complex plane without precise knowledge of the approximation domain. One may also question the principal usage of CFETs for dissipative systems, where dynamical semigroups replace the Lie group setting. Modifications of CFETs can be tailored to this situation and circumvent the negative time-step problem that occurs for methods beyond the 4th-order CFETs to which we restricted our present considerations.
A more physical question concerns the use of the Lindblad formalism in the description of optical emission.
The Markovian approximation is not entirely satisfactory here, since it cannot distinguish between energy increasing ("virtual") and energy decreasing ("real") transitions in the emission process. For this reason, the total emission rate remains finite (of the order of the cavity loss rate κ) even without external pumping although it should drop to zero. The use of non-Markovian master equations [2,17] might overcome these problems, which should also be relevant if we ask for the coherence properties of the emitted light [3,4]. Another possibility is the combination of CFETs with polynomial techniques for the numerical representation of open quantum systems [29][30][31].
We hope to be able to return to some of these issues soon, which we had to leave unresolved here. Presently, we can conclude that CFETs are one promising contribution to numerical time propagation of complex non-equilibrium quantum systems and warrant further exploration.