# Non-Hermitian wave packet approximation of Bloch optical equations

###### Abstract

We introduce a non-Hermitian approximation of Bloch optical equations. This approximation provides a complete description of the excitation, relaxation and decoherence dynamics of ensembles of coupled quantum systems in weak laser fields, taking into account collective effects and dephasing. In the proposed method one propagates the wave function of the system instead of a complete density matrix. Relaxation and dephasing are taken into account via automatically-adjusted time-dependent gain and decay rates. As an application, we compute the numerical wave packet solution of a time-dependent non-Hermitian Schrödinger equation describing the interaction of electromagnetic radiation with a quantum nano-structure and compare the calculated transmission, reflection, and absorption spectra with those obtained from the numerical solution of the Liouville- von-Neumann equation. It is shown that the proposed wave packet scheme is significantly faster than the propagation of the full density matrix while maintaining small error. We provide the key ingredients for easy-to-use implementation of the proposed scheme and identify the limits and error scaling of this approximation.

###### pacs:

42.50.Ct, 78.67.-n, 32.30.-r, 33.20.-t, 36.40.Vz, 03.65.Yz## I Introduction

Optics of nanoscale materials has attracted considerable attention in the past several years Maier and Atwater (2005); Lal et al. (2007); Ozbay (2006); Berini (2009) due to various important applications Stockman (2011). Exploring electrodynamics of near-fields associated with subwavelength systems, researchers are now truly dwelling into nanoscale Maier et al. (2001); Garcia de Abajo (2007); Ebbesen et al. (2008); Schuller et al. (2010); Stockman (2011). Owing to both new materials processing techniques Hutter and Fendler (2004) and continuous progress in laser physics Slusher (1999) the research in nano-optics is currently transitioning from linear systems, where materials and their relative arrangement control optical properties Barnes and Murray (2007), to the nonlinear regime Kroo et al. (2008). The latter expands optical control capabilities far beyond conventional linear optics as in the case of active plasmonic materials Wuestner et al. (2011), for instance, combining highly localized electromagnetic (EM) radiation driven by surface plasmon-polaritons (SPP) with non-linear materials Krasavin et al. (2006). Yet another promising research direction, namely optics of highly coupled exciton-polariton systems, is emerging Bellessa et al. (2004); Chang et al. (2006). It basically reincarnates a part of research in semiconductors Agranovich and La Rocca (2005); Khitrova et al. (2006), bringing it to nanoscale via deposition of ensembles of quantum emitters (molecules Dintinger et al. (2005); Ebbesen et al. (2009); Lekeufack et al. (2010); Berrier et al. (2011), quantum dots Park et al. (2007); Komarala et al. (2008); Andersen et al. (2011); Livneh et al. (2011)) directly on to plasmonic materials.

Even in the linear regime, when the external EM radiation is not significantly exciting the quantum sub-system, SPP near-fields can be strongly coupled to quantum emitters. This manifests itself as a Rabi splitting widely observed in transmission experiments Gomez et al. (2010). Moreover, a new phenomenon, namely collective molecular modes driven by SPP near-fields, has been observed Ebbesen et al. (2009) and recently explained Salomon et al. (2012). It was also shown that nanoscale clusters comprised of optically coupled quantum emitters exhibit collective scattering and absorption Sukharev and Nitzan (2011) similar to Dicke superradiance A. V. Andreev and Il’inskiĭ (1980). It is hence important to be able to account for collective effects in a self-consistent manner. We also note that a series of works by Neuhauser et al. Gupta and Neuhauser (2001); Lopata and Neuhauser (2007); Lopata and Neuhauser (2009a, b) clearly demonstrated that the presence of a single molecule nearby a plasmonic material can significantly alter the scattering spectra.

In many applications (such as optics of molecular layers coupled to plasmonic materials Dintinger et al. (2005); Ebbesen et al. (2009); Lekeufack et al. (2010); Berrier et al. (2011); Salomon et al. (2012), for instance) self-consistent modeling relies heavily on the numerical integration of the corresponding Maxwell-Bloch equations, assuming that static emitter-emitter interactions can be neglected, which is true for systems at relatively low densities. Such an approximation results in expressing the local polarization in terms of a product of the local density of quantum emitters and the local averaged single emitter’s dipole moment Allen and Eberly (1975). One of the first efficient numerical schemes for simulations of nonlinear optical phenomena of quantum media driven by external classical EM radiation was proposed by Ziolkowski et al. Ziolkowski et al. (1995). Using a one-dimensional example of ensembles of two-level atoms it was shown that the corresponding Maxwell-Bloch equations can be successfully integrated using an iterative scheme based on the predictor-corrector method (strongly coupled method). Later on this approach has been extended to two- Slavcheva et al. (2002) and three-dimensional systems Fratalocchi et al. (2008). Although such a scheme accurately captures the system’s dynamics, it can become extensively slow for multidimensional systems Sukharev and Nitzan (2011). Moreover this method is limited to two-level systems only. A more efficient technique based on the decoupling of Bloch equations from the Ampere law was proposed in 2003 by Bidégaray Bidégaray (2003). This latter method, usually referred to as a weakly coupled method, noticeably improves the efficiency of the numerical integration of Maxwell-Bloch equations, allowing to consider multilevel quantum media Sukharev and Nitzan (2011).

The approach proposed in the present paper is based on this weakly coupled method Sukharev and Nitzan (2011). It further improves the numerical efficiency of the Maxwell-Bloch integrator for ensembles of multilevel quantum emitters. By incorporating a new non-Hermitian wave packet propagation technique into the weakly coupled method, we demonstrate that our approach can be successfully applied to ensembles of multilevel atoms and diatomic molecules. We show with this method that it is sufficient to propagate a single wave function instead of the complete density matrix. Relaxation and dephasing are taken into account via empirical gain and decay rates whose time-dependence is automatically adapted for an optimal description of dephasing processes.

The paper is organized as follows. We first introduce our model in section II.1, using, as an example, an ensemble of two-level atoms. The non-Hermitian wave packet approximation is then described in section II.2. Applications of the proposed method to the case of a nano-layer comprised of two-level atoms are discussed in section II.3. We then generalize our method to the case of interacting multi-level emitters in section II.4. Finally, using this generalized approach, we calculate in section II.5 the optical properties of a nano-layer of coupled molecules, taking into account both the vibrational and rotational degrees of freedom, and we reveal several interesting features in the absorption spectra. Last section III summarizes our work.

## Ii Theoretical model and applications

As pointed out in the previous section, the importance of collective effects manifested in recent experiments and theoretical papers calls for the development of self-consistent models capable of taking into account mutual EM interaction in ensembles of quantum emitters. While direct numerical integration of corresponding Maxwell-Bloch equations can be based on either a strong coupled method Ziolkowski et al. (1995) or a weakly coupled one Bidégaray (2003), for multi-level systems in two- and especially three-dimensions such a brute-force approach becomes numerically very expensive, both in terms of computation (CPU) times and memory requirements. Indeed, the main disadvantage of such approaches is that both the CPU times and memory requirements scale generally at least as , where is the dimension of the Hilbert space of the system. This unfavorable scaling is directly related to the size of the reduced density matrix used to describe the system and to the associated number of time-dependent equations one has to solve to follow the system’s dynamics. By contrast, one has only to solve a reduced set of time-dependent equations when the system can be described by a single wave function. It would therefore be extremely useful to be able to derive a “Schrödinger-type” approximation which could describe on an equal footing the field-induced coherent dynamics of a multi-level quantum system and the associated relaxation and decoherence processes. This type of idea is not entirely new. For small excitations, the time-dependent density functional theory Runge and Gross (1984) allows for instance, in its real-time version, the extraction of density matrix dephasing without evolving the full density matrix. This goal has also been achieved in the past using different approaches, the three main contributions being the stochastic Schrödinger method used in conjunction with a Monte Carlo integrator Dalibard et al. (1992); Dum et al. (1992a); Gardiner et al. (1992); Dum et al. (1992b); Gisin and Percival (1992); Makarov and Metiu (1999), the Gadzuk jumping wave packet scheme Gadzuk et al. (1990); Finger and Saalfrank (1997) and the variational wave packet method Gerdts and Manthe (1997); Pesce et al. (1998). These methods require the propagation of different wave functions and are therefore mainly attractive when .

The approach we propose here is different since it is based on perturbation theory. It also significantly speeds up calculations since it only requires the propagation of a single wave function under the action of an easy-to-implement time-dependent effective Hamiltonian in order to reproduce accurately the full dynamics. We demonstrate the efficiency and accuracy of our method using, first, ensembles of two-level atoms, and then, interacting diatomic molecules including both the vibrational and rotational degrees of freedom. We also discuss the conditions under which our method is no longer valid. All results are compared with data obtained via direct integration of Maxwell-Bloch equations using a weakly coupled method Sukharev and Nitzan (2011).

### ii.1 Atomic two-level system

Let us first consider a two-level quantum system interacting with EM radiation. We label the two levels as and , with associated energy eigenvalues and , respectively. The corresponding density matrix satisfies the dissipative Liouville-von Neumann equation Blum (1996)

(1) |

where is the total Hamiltonian and is a superoperator, taken in the Lindblad form Breuer and Petruccione (2002), describing relaxation and dephasing processes under Markov approximation. The field free Hamiltonian reads

(2) |

The interaction of the two-level system with EM radiation is taken in the form

(3) |

denotes here the instantaneous Rabi frequency associated with the coupling between the quantum system and an external field. In Eq. (1), the non-diagonal elements of the operator include a pure dephasing rate , and the diagonal elements of this operator consist of the radiationless decay rate of the excited state. Equations (1)-(3) lead to the well-known Bloch optical equations Allen and Eberly (1975) describing the quantum dynamics of a coupled two-level system

(4a) | |||||

(4b) | |||||

(4c) | |||||

(4d) |

where is the Bohr transition frequency and where the dot denotes the time derivative.

We assume in the following that the system is initially in the ground state , and we will show that, under some assumptions, the subsequent induced excitation and relaxation dynamics can be accurately described by a non-Hermitian wave packet approximation.

### ii.2 Non-Hermitian two-level wave packet approximation

Within the aforementioned simplified model, the system’s wave packet can be expanded as

(5) |

Let us now assume that the ground and excited states energies include an imaginary part that we denote as and , respectively. Inserting expansion (5) into the time-dependent Schrödinger equation

(6) |

and projecting it onto states and yields the following set of coupled equations for the coefficients

(7a) | |||||

(7b) |

We will now derive the differential equations describing the temporal dynamics of the products , where the subscript corresponds to this simplified non-Hermitian Schrödinger-type model. Our goal is to obtain the gain and decay rates and which will allow for an approximate description of the system’s excitation and relaxation dynamics. For the evolution of the populations and , one gets

(8a) | |||||

(8b) |

The conservation of the total norm

(9) |

thus results in

(10) |

It therefore appears that, since the populations and are generally time dependent, at least one of the two rates or , which are not yet fully defined, must be taken as a time-dependent function.

Taking Eq. (10) into account, one finally obtains the following set of Bloch equations for the approximate density matrix

(11a) | |||||

(11b) | |||||

(11c) | |||||

(11d) |

Two obvious choices can then be made for the empirical gain and decay parameters and :

The optimal choice for applications in linear optics of nano-materials should be as follows: in weak fields, the variation of the populations, as described by perturbation theory, is a second order term with respect to the coupling amplitude while the variation of the coherences is a first order term. For a correct description of the quantum dynamics in weak fields, it is important to describe first order terms accurately. Therefore we proceed with the second choice

(12) |

From Eqs.(10) and (12) we obtain a set of empirical time dependent gain and decay rates and that we can insert in the Schrödinger-type approximation

(13a) | |||||

(13b) |

We finally obtain the following set of time-dependent coupled equations

(14a) | |||||

(14b) |

This system is solved numerically using a fourth order Runge-Kutta algorithm Press et al. (2002). Eqs. (14a), (14b) include two non-Hermitian terms which can (as we will illustrate below) accurately reproduce the dissipative dynamics in weak fields, i.e. for . Indeed, as one can easily show that, in this limit, Eqs. (11a)-(11d) are strictly equivalent to the exact Bloch Eqs. (4a)-(4d).

At first sight, it might seem that the proposed formalism rather relies on a mathematical trick. It does however have physical insights. Indeed, the time-dependent gain and decay rates, in some cases, have a deep physical meaning, as it was shown in Ref. Klaiman et al. (2008), for instance, where the non-Hermitian formalism was applied to study EM wave propagation in so-called PT-symmetric waveguides. In this study, it was shown that the gain and loss coefficients could be used to control the beat length parameter which describes waveguides. In our approach and are gain and decay rates inherent to the studied quantum system since they only depend on the decay and dephasing rates and of this system and on the relative population of the quantum states involved. In the past, such type of non-Hermitian approaches have been proven to be very powerful tools in many branches of physics, including resonant phenomena, quantum mechanics, optics, and quantum field theory. A recent comprehensive survey of various applications of the non-Hermitian approach can be found in Ref. Moiseyev (2011).

### ii.3 Application to a uniform nano-layer of atoms

While our approximation is equally applicable to three-dimensional systems, we have chosen, for the sake of simplicity, to test it on a simplified one-dimensional system consisting of a uniform infinite layer of atoms whose thickness lies in the range of a few hundred nanometers. An incident radiation field propagating in the positive -direction is represented by a transverse-electric mode with respect to the propagation axis. It is characterized by a single in-plane electric field component and a single in-plane magnetic field component . To account for the symmetry of the atomic polarization response, the atoms in the layer are described as two-level systems with the following states: an -type ground state and a -type excited state. This model is a one-dimensional simplification of a more general approach used in Ref. Sukharev and Nitzan (2011) and we refer the reader to the body of this paper for details.

The time-domain Maxwell’s equations for the dynamics of the electromagnetic fields

(15a) | |||||

(15b) |

are solved using a generalized finite-difference time-domain technique where both the electric and magnetic fields are propagated in discretized time and space Taflove and Hagness (2005). In these equations, and denote the magnetic permeability and dielectric permittivity of free space. The macroscopic polarization of the atomic system

(16) |

is taken as the expectation value of the atomic transition dipole moment , where is the atomic density.

A self-consistent model is based on the numerical integration of Maxwell’s equations (15a) and (15b), coupled via Eq. (16) to the quantum dynamics. In the mean-field approximation employed here it is assumed that the density matrix of the ensemble is expressed as a product of density matrices of individual quantum emitters driven by a local EM field. In order to account for dipole-dipole interactions within a single grid cell we follow Bowden and Dowling (1993) and introduce Lorentz-Lorenz correction for a local electric field term into quantum dynamics according to

(17) |

where is the solution of Maxwell’s equations (15a), (15b) and macroscopic polarization is evaluated according to Eq. (16).

The quantum dynamics is evaluated by computing the atomic dipole moment either using the single-atom density matrix (4a)-(4d) or the single-atom wave packet (14a)-(14b). To compare two approaches in the linear regime (i.e. for , ), we calculate the transmission , reflection and absorption spectra of an atomic layer of thickness nm as a function of the incident photon energy for an atomic transition energy of eV. The results are shown in Fig. 1.

The transmission and reflection spectra are obtained from the normalized Poynting vector

(18) |

at a specific location under and above the atomic layer, respectively, where , and , are the Fourier components of the total and incident EM fields. The absorption spectrum is then simply obtained as

(19) |

In Fig. 1 the solution of Maxwell-Liouville-von Neumann equations is shown as a blue solid line while the open red squares are from the solution of our approximate non-Hermitian Schrödinger model. One can notice a perfect agreement of the two methods irrespective of the atomic density. These results clearly demonstrate that, in the linear regime where the atomic excitation probability remains small and varies linearly with the incident field intensity, a simple wave packet propagation is sufficient to mimic the excitation and relaxation dynamics of an ensemble of quantum emitters. Therefore in the weak field regime, the propagation of the full density matrix is unnecessary.

At low density, most of the incident radiation is simply transmitted and a small part of the incident energy is absorbed by the atomic ensemble at energies close to the atomic transition energy of 2 eV. The absorption spectrum shows the conventional Lorentzian profile. At higher densities, the absorption spectrum is strongly modified due to the appearance of collective excitation modes Sukharev and Nitzan (2011), and a large part of the incident energy is either absorbed or reflected from the atomic nano-layer at photon energies close to the transition energy.

The blue lines with squares in Fig. 2 show, as a function of the incident field intensity, the relative error obtained in the calculation of the absorption spectrum at the transition energy using the Schrödinger approximation when compared to the solution of the full Liouville-von Neumann equation . The solid blue line is for the lowest atomic density m and the dashed blue line is for the highest density m.

One can see that irrespective of the atomic density the relative error of the non-Hermitian wave packet approximation scales linearly with the field intensity. This is not surprising. Indeed, solving coupled Schrödinger equations (14a)-(14b) is strictly equivalent to solving the approximate Bloch equations (11a)-(11d). The latter differ from the exact Bloch equations (4a)-(4d) by a term proportional to the excited state population . The maximum excited state population is also shown in Fig. 2 as a function of the field intensity. It can be seen that it also varies linearly in the present weak field regime. We can conclude (from Fig. 2 and numerous calculations we have performed) that as long as the excited state population remains smaller than 1% our wave packet approximation can be used absolutely safely.

As shown in Fig. 3, we could verify that the quality of the calculated spectra is still rather good when the excited state population reaches 35%. This figure shows the reflection spectra calculated using the “exact” Liouville-von-Neumann equations and using our approximate Schrödinger model for a relatively high laser field amplitude, chosen such that the maximum excited state population reaches 35%. It is only when the excited state population approaches 50% that one can observe a very sudden failure of the present Schrödinger model, as could be expected from the divergence of the time-dependent gain and decay rates (13a) and (13b) when . This example shows that the present non-Hermitian Schrödinger approximation still holds in the case of relatively large couplings.

The time dependence of the excited state population and of the coherence is illustrated in panels (a), (b) and (d) of Fig. 4. These results were obtained with the same parameters as in Fig. 1, with an atomic density of m. One can see from the panels (b) and (d) of this figure, showing the square modulus and the real part of the coherence , that our non-Hermitian Schrödinger model reproduces quite accurately the coherence dynamics of the system. On the other hand, as shown in panel (a), this is obtained at the cost of a poor description of the excited state population dynamics. Indeed, to describe correctly the coherence of the system one is led to overestimate the excited state decay rate. However, as seen in Fig. 1, this overestimation of the decay rate does not have any impact on the accuracy of the calculated absorption, reflection and transmission spectra when the excited state population remains small compared to the ground state population. Finally, panel (c) of Fig. 4 shows the time dependence of the gain coefficient . This gain rate basically follows the evolution of the excited state population, as it could already be inferred from Eq. (13a). The decay rate is not shown in this figure since it is essentially constant and equal to (see Eq. (13b)) in the present situation with .

### ii.4 Generalization to multi-level systems

In order to generalize our Schrödinger-type approximation of the excitation and dissipation dynamics of an ensemble of quantum emitters to a multi-level system we now consider the case of neutral diatomic molecules. More specifically, we have chosen the particular case of the ground and first excited electronic states of the Li molecule and we follow the electronic dynamics and the nuclear motion by expanding the total molecular wave function using the Born-Oppenheimer expansion

(20) |

where and denote the electronic wave functions associated with the ground and first excited electronic states of Li, respectively. The electron coordinates are denoted by the vector , and the vector represents the internuclear vector.

We now separate the global electronic coordinate of all electrons into the coordinate of the core electrons and the coordinate r of the active electron Charron and Raoult (2006). The ground electronic state is considered as a 2s state, and the electronic wave function is expressed in the molecular frame (Hund’s case (b) representation) as the product

(21) |

where and are the radial and angular parts of the electronic wave function associated with the active electron. Similarly, the 2p excited state of symmetry is expressed in the molecular frame as

(22) |

Due to the symmetry of both electronic states, the ro-vibrational time-dependent wave functions and can be expanded on a limited set of normalized Wigner rotation matrices in order to take into account the rotational degree of freedom, following Charron et al. (1994)

(23a) | |||||

(23b) |

where denotes the molecular rotational quantum number while denotes its projection on the electric field polarization axis of the laboratory frame.

Introducing these expansions in the time-dependent Schrödinger equation (6) describing the molecule-field interaction and projecting onto the electronic and rotational basis functions yields, in the dipole approximation, the following set of coupled differential equations for the nuclear wave packets and

(24a) | |||||

(24b) | |||||

where is the electronic transition dipole. The ro-vibrational nuclear Hamiltonians are defined as

(25) |

where denotes the molecular reduced mass. The potential energy curves and associated with the ground and first excited electronic states of the molecule are taken from Ref. Schmidt-Mink et al. (1985) and the matrix elements which couple the nuclear wave packets evolving on these electronic potential curves can be written using -symbols as

(26) | |||||

We assume that the molecules are prepared at time in the ro-vibrational level and of the ground electronic state. In weak linearly polarized fields and except for a phase factor this ground state component of the molecular wave function remains unaffected and the excited state component is limited to and . The ground and excited nuclear wave packets are thus finally expanded in terms of ro-vibrational eigenstates as

(27a) | |||||

(27b) |

where and are time-dependent complex coefficients. denote here the bound ro-vibrational eigenstates of the electronic potentials Charron and Suzor-Weiner (1998). It is not necessary here to take into account the dissociative nuclear eigenstates associated with the excited potential since their coupling with the ground vibrational level of the ground electronic state is negligible.

We thus arrive at a multi-level system which is very similar to the atomic case described in sections II.1 and II.2, except that the initial ground state is now coupled with a large set of excited levels. For convenience and for an easy comparison with the atomic case, we will label the ground state as state number 0 and the excited states as states number . Our reference calculations will be based on the numerical solutions of the corresponding Liouville-von Neumann equations

(28a) | |||||

(28b) | |||||

(28c) | |||||

(28d) |

where is the total energy of the excited state and where denotes the instantaneous Rabi frequency associated with the molecule-field interaction as defined in Eqs. (24a) and (24b). and denote the pure dephasing rate and the relaxation rate associated with the excited states, respectively.

In comparison with this “exact” model, our Schrödinger-type approximation will be based on the numerical solutions of the coupled equations for the time-dependent expansion coefficients

(29a) | |||||

(29b) |

where the time-dependent gain and decay rates and are now defined as

(30) | |||||

(31) |

One can show that with such a definition of and , Eqs. (29a) and (29b) are strictly equivalent to Eqs. (28a)-(28d) in the limit of weak couplings. In the next section we demonstrate that our Schrödinger-type model accurately reproduces the excitation and dissipation dynamics of multi-level quantum systems in the limit of weak couplings, i.e. for .

### ii.5 Application to a uniform nano-layer of molecules

We consider a simplified one-dimensional system similar to the one discussed in case of two-level atoms. A uniform infinite layer with a thickness of nm comprised of Li molecules is exposed to incident linearly polarized radiation. The incident field propagates in the positive -direction. We calculate the absorption spectrum of this molecular layer just as we did in section II.3 for atoms.

To ascertain the validity of the proposed Schrödinger-type approximation, absorption spectra are represented in Fig. 5 as a function of the incident photon energy , in the linear regime, for two different molecular densities: m in the left column (panels (a) and (b)) and m in the right column (panels (c) and (d)). The solutions obtained via integrating Maxwell-Liouville-von Neumann equations are shown in the upper panels (a) and (c) as blue solid lines while the ones obtained from our approximate non-Hermitian Schrödinger model are shown in red in the lower panels (b) and (d) (inverted spectra).

As in the case of two-level atoms, one can notice a perfect agreement of the two methods irrespective of the molecular density. This shows that, in the linear regime where the molecular excitation probability remains small and varies linearly with the incident field intensity, a simple wave packet propagation is sufficient to mimic the excitation, dissipation and decoherence dynamics of an ensemble of multi-level quantum emitters. The propagation of the full density matrix is then, again, unnecessary.

It is important to note that at low density, we observe a series of overlapping vibrational resonances which reflects the vibrational structure of the excited molecular potential and which follows the Franck-Condon principle Franck and Dymond (1926); Condon (1926). The green labels seen in panel (a) of Fig. 5 indicate the excited state vibrational level responsible for the observed resonance.

At higher density, the absorption spectrum is strongly distorted and one observes, just like in the atomic case Sukharev and Nitzan (2011), the appearance of collective excitation modes, the difference being that a vibrational structure may still be present in some of these collective molecular excitation modes. A detailed analysis of the physics underlying the appearance of these intriguing molecular collective modes will be presented in another paper. For high densities, we could also observe that the radiation field is not transmitted anymore through the molecular layer in the absorption window of the molecule. Within this window, the field is either absorbed or reflected. The small inset seen in panel (c) of Fig. 5 shows a magnification of the absorption spectrum in the energy range corresponding to the “normal” low-density absorption spectrum. One can see in this inset, and by comparing with panel (a) of Fig. 5, that the absorption spectrum is not strongly modified at high molecular density in this energy region.

Figure 6 finally shows the computation time necessary for the calculation of these absorption spectra as a function of the number of quantum states introduced in this multi-level model for both the Maxwell-Liouville-von Neumann (blue line with circles) and Maxwell-Schrödinger (red line with squares) approaches. One can observe a substantial difference in computation times which can prove of crucial importance when one has to deal with realistic three-dimensional systems. This origin of the observed gain in computation time relies on the necessity of propagating a single wave function instead of a full density matrix.

## Iii Summary and conclusions

We proposed a new and simple non-Hermitian approximation of Bloch optical equations where one propagates the wave function of the quantum system instead of the complete density matrix. Our method provides an accurate, complete description of the excitation, relaxation and decoherence dynamics of single as well as ensembles of coupled quantum emitters (atoms or molecules) in weak EM fields, taking into account collective effects and dephasing. We demonstrated the applicability of the method by computing optical properties of thin layers comprised of two-level atoms and diatomic molecules. It was shown that, in the limit of weak incident fields, the dynamics of interacting quantum emitters can be successfully described by our set of approximated equations, which result in a substantial gain both in computation time and computer memory requirements. These calculations also reveal some intriguing new collective molecular excitation modes which will be presented in detail in another publication. The proposed approach was demonstrated to provide a substantial increase in numerical efficiency for self-consistent simulations.

## Iv Acknowledgements

E.C. would like to acknowledge useful and stimulating discussions with O. Atabek and A. Keller from Université Paris-Sud (Orsay) and with E. Shapiro from the University of British Columbia (Canada). M.S. is grateful to the Université Paris-Sud (Orsay) for the financial support through an invited Professor position in 2011. E.C. acknowledges supports from ANR (contract Attowave ANR-09-BLAN-0031-01), and from the EU (Project ITN-2010-264951, CORINF). We also acknowledge the use of the computing facility cluster GMPCS of the LUMAT federation (FR LUMAT 2764).

## References

- Maier and Atwater (2005) S. A. Maier and H. A. Atwater, J. Appl. Phys. 98, 011101 (2005).
- Lal et al. (2007) S. Lal, S. Link, and N. J. Halas, Nat. Photon. 1, 641 (2007).
- Ozbay (2006) E. Ozbay, Science 311, 189 (2006).
- Berini (2009) P. Berini, Adv. Opt. Photon. 1, 484 (2009).
- Stockman (2011) M. I. Stockman, Opt. Express 19, 22029 (2011).
- Maier et al. (2001) S. A. Maier, M. L. Brongersma, P. G. Kik, S. Meltzer, A. A. G. Requicha, and H. A. Atwater, Adv. Mater. 13, 1501 (2001).
- Garcia de Abajo (2007) F. J. Garcia de Abajo, Rev. Mod. Phys. 79, 1267 (2007).
- Ebbesen et al. (2008) T. W. Ebbesen, C. Genet, and S. I. Bozhevolnyi, Phys. Today 61, 44 (2008).
- Schuller et al. (2010) J. A. Schuller, E. S. Barnard, W. S. Cai, Y. C. Jun, J. S. White, and M. L. Brongersma, Nat. Mater. 9, 193 (2010).
- Hutter and Fendler (2004) E. Hutter and J. H. Fendler, Adv. Mater. 16, 1685 (2004).
- Slusher (1999) R. E. Slusher, Rev. Mod. Phys. 71, S471 (1999).
- Barnes and Murray (2007) W. L. Barnes and W. A. Murray, Adv. Mater. 19, 3771 (2007).
- Kroo et al. (2008) N. Kroo, S. Varro, G. Farkas, P. Dombi, D. Oszetzky, A. Nagy, and A. Czitrovszky, J. Mod. Optic. 55, 3203 (2008).
- Wuestner et al. (2011) S. Wuestner, A. Pusch, K. L. Tsakmakidis, J. M. Hamm, and O. Hess, Philos. T. R. Soc. A 369, 3525 (2011).
- Krasavin et al. (2006) A. V. Krasavin, K. F. MacDonald, A. S. Schwanecke, and N. I. Zheludev, Appl. Phys. Lett. 89, 031118 (2006).
- Bellessa et al. (2004) J. Bellessa, C. Bonnand, J. C. Plenet, and J. Mugnier, Phys. Rev. Lett. 93, 036404 (2004).
- Chang et al. (2006) D. E. Chang, A. S. Sorensen, P. R. Hemmer, and M. D. Lukin, Phys. Rev. Lett. 97, 053002 (2006).
- Agranovich and La Rocca (2005) V. M. Agranovich and G. C. La Rocca, Solid State Commun. 135, 544 (2005).
- Khitrova et al. (2006) G. Khitrova, H. M. Gibbs, M. Kira, S. W. Koch, and A. Scherer, Nat. Phys. 2, 81 (2006).
- Dintinger et al. (2005) J. Dintinger, S. Klein, F. Bustos, W. L. Barnes, and T. W. Ebbesen, Phys. Rev. B 71, 035424 (2005).
- Ebbesen et al. (2009) T. W. Ebbesen, A. Salomon, and C. Genet, Angew. Chem. Int. Edit. 48, 8748 (2009).
- Lekeufack et al. (2010) D. D. Lekeufack, A. Brioude, A. W. Coleman, P. Miele, J. Bellessa, L. D. Zeng, and P. Stadelmann, Appl. Phys. Lett. 96, 253107 (2010).
- Berrier et al. (2011) A. Berrier, R. Cools, C. Arnold, P. Offermans, M. Crego-Calama, S. H. Brongersma, and J. Gomez-Rivas, ACS Nano 5, 6226 (2011).
- Park et al. (2007) H. Park, A. V. Akimov, A. Mukherjee, C. L. Yu, D. E. Chang, A. S. Zibrov, P. R. Hemmer, and M. D. Lukin, Nature 450, 402 (2007).
- Komarala et al. (2008) V. K. Komarala, A. L. Bradley, Y. P. Rakovich, S. J. Byrne, Y. K. Gun’ko, and A. L. Rogach, Appl. Phys. Lett. 93, 123102 (2008).
- Andersen et al. (2011) M. L. Andersen, S. Stobbe, A. S. Sorensen, and P. Lodahl, Nat. Phys. 7, 215 (2011).
- Livneh et al. (2011) N. Livneh, A. Strauss, I. Schwarz, I. Rosenberg, A. Zimran, S. Yochelis, G. Chen, U. Banin, Y. Paltiel, and R. Rapaport, Nano Lett. 11, 1630 (2011).
- Gomez et al. (2010) D. E. Gomez, K. C. Vernon, P. Mulvaney, and T. J. Davis, Nano Lett. 10, 274 (2010).
- Salomon et al. (2012) A. Salomon, R. J. Gordon, Y. Prior, T. Seideman, and M. Sukharev, Phys. Rev. Lett. 109, 073002 (2012).
- Sukharev and Nitzan (2011) M. Sukharev and A. Nitzan, Phys. Rev. A 84 (2011).
- A. V. Andreev and Il’inskiĭ (1980) V. I. E. A. V. Andreev and Y. A. Il’inskiĭ, Sov. Phys. Usp. 23, 493 (1980).
- Gupta and Neuhauser (2001) A. Gupta and D. Neuhauser, International J. Quant. Chem. 81, 260 (2001).
- Lopata and Neuhauser (2007) K. Lopata and D. Neuhauser, J. Chem. Phys. 127, 154715 (2007).
- Lopata and Neuhauser (2009a) K. Lopata and D. Neuhauser, J. Chem. Phys. 130 (2009a).
- Lopata and Neuhauser (2009b) K. Lopata and D. Neuhauser, J. Chem. Phys. 131, 014701 (2009b).
- Allen and Eberly (1975) L. Allen and J. H. Eberly, Optical resonance and two-level atoms (Wiley, New York, 1975).
- Ziolkowski et al. (1995) R. W. Ziolkowski, J. M. Arnold, and D. M. Gogny, Phys. Rev. A 52, 3082 (1995).
- Slavcheva et al. (2002) G. Slavcheva, J. M. Arnold, I. Wallace, and R. W. Ziolkowski, Phys. Rev. A 66, 063418 (2002).
- Fratalocchi et al. (2008) A. Fratalocchi, C. Conti, and G. Ruocco, Phys. Rev. A 78, 013806 (2008).
- Bidégaray (2003) B. Bidégaray, Numer. Meth. Part. D. E. 19, 284 (2003).
- Runge and Gross (1984) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
- Dalibard et al. (1992) J. Dalibard, Y. Castin, and K. Mølmer, Phys. Rev. Lett. 68, 580 (1992).
- Dum et al. (1992a) R. Dum, P. Zoller, and H. Ritsch, Phys. Rev. A 45, 4879 (1992a).
- Gardiner et al. (1992) C. W. Gardiner, A. S. Parkins, and P. Zoller, Phys. Rev. A 46, 4363 (1992).
- Dum et al. (1992b) R. Dum, A. S. Parkins, P. Zoller, and C. W. Gardiner, Phys. Rev. A 46, 4382 (1992b).
- Gisin and Percival (1992) N. Gisin and I. C. Percival, J. Phys. A: Math. Gen. 25, 5677 (1992).
- Makarov and Metiu (1999) D. E. Makarov and H. Metiu, J. Chem. Phys. 111, 10126 (1999).
- Gadzuk et al. (1990) J. Gadzuk, L. Richter, S. Buntin, D. King, and R. Cavanagh, Surface Science 235, 317 (1990).
- Finger and Saalfrank (1997) K. Finger and P. Saalfrank, Chem. Phys. Lett. 268, 291 (1997).
- Gerdts and Manthe (1997) T. Gerdts and U. Manthe, J. Chem. Phys. 106, 3017 (1997).
- Pesce et al. (1998) L. Pesce, T. Gerdts, U. Manthe, and P. Saalfrank, Chem. Phys. Lett. 288, 383 (1998).
- Blum (1996) K. Blum, Density matrix theory and applications (Plenum Press, New York, 1996), 2nd ed.
- Breuer and Petruccione (2002) H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, New York, 2002).
- Press et al. (2002) W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes: The Art of Scientific Computing (Cambridge University Press, Cambridge, 2002).
- Klaiman et al. (2008) S. Klaiman, U. Günther, and N. Moiseyev, Phys. Rev. Lett. 101, 080402 (2008).
- Moiseyev (2011) N. Moiseyev, Non-Hermitian quantum mechanics (Cambridge University Press, 2011).
- Taflove and Hagness (2005) A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference time-domain method (Artech House, Boston, 2005), 3rd ed.
- Bowden and Dowling (1993) C. M. Bowden and J. P. Dowling, Phys. Rev. A 47, 1247 (1993).
- Charron and Raoult (2006) E. Charron and M. Raoult, Phys. Rev. A 74, 033407 (2006).
- Charron et al. (1994) E. Charron, A. Giusti-Suzor, and F. H. Mies, Phys. Rev. A 49, R641 (1994).
- Schmidt-Mink et al. (1985) I. Schmidt-Mink, W. Müller, and W. Meyer, Chem. Phys. 92, 263 (1985).
- Charron and Suzor-Weiner (1998) E. Charron and A. Suzor-Weiner, J. Chem. Phys. 108, 3922 (1998).
- Franck and Dymond (1926) J. Franck and E. G. Dymond, Trans. Faraday Soc. 21, 536 (1926).
- Condon (1926) E. Condon, Phys. Rev. 28, 1182 (1926).