The role of dephasing for dark state coupling in a molecular Tavis–Cummings model

The collective coupling of an ensemble of molecules to a light field is commonly described by the Tavis–Cummings model. This model includes numerous eigenstates that are optically decoupled from the optically bright polariton states. Accessing these dark states requires breaking the symmetry in the corresponding Hamiltonian. In this paper, we investigate the influence of non-unitary processes on the dark state dynamics in the molecular Tavis–Cummings model. The system is modeled with a Lindblad equation that includes pure dephasing, as it would be caused by weak interactions with an environment, and photon decay. Our simulations show that the rate of pure dephasing, as well as the number of two-level systems, has a significant influence on the dark state population.

One of many challenges for this developing field is to build models that incorporate the details relevant for chemical reactivity (e.g.molecular vibrations along reaction coordinates) together with collective ensemble effects (several molecules coupled to the cavity field) [27,28].Thus, in the last decade, several studies have turned their attention to such models of emitter ensembles [24,27,29,30].
In these ensemble models, a common feature is the presence of dark states, which do not couple to the field of the cavity or an external field.These states are most straightforward to consider in a Tavis-Cummings model limited to a single excitation [31]: With N identical emitters, there are two bright polaritonic states, |UP , and |LP , which are symmetrical under emitter permutation.The remaining N − 1 polaritonic states are dark.See Fig. 1.In the Tavis-Cummings model, dark states are: degenerate [32], asymmetrical under emitter permutation [33,34], have no transition dipole moment, and correspond to collective emitter excitations, which gain no population during Hamiltonian evolution [24].Note that these dark states emerge in the ensemble system, which is different from disallowed transitions in individual emitters.The absence of population in the dark states makes intuitive sense; assuming an initial symmetrical state, the Hamiltonian treats all emitters identically, so unitary evolution cannot induce asymmetry.The existence of the dark states, as well as their properties, are also quite robust; for example, they persist even if the field couplings vary among emitters [35].In this work, the model will be an extended version of Tavis-Cummings model, but properties of the dark states still hold (degeneracy, asymmetry, non-existent transition dipole moment, and zero population under Hamiltonian evolution).
The role of dephasing processes and its influence on darkstate dynamics is intertwined with identifying a mechanism for their population [48].To lift the decoupling of the dark states from the bright polariton states, an asymmetry in the Hamiltonian (or the Liouvilian) is required.To find such a process, we can consider processes that operate at the level of individual emitters, such as the interactions with each emitter's environment.We model this with dephasing operators, which can be included in the equations of motion by means of the Lindblad equation: The emitter dephasing can be motivated by fluctuations in the relative energies between the states in the emitter systems [60], or by a build-up of system-environment correlations [61].
In the context of strongly coupled molecular systems, there has been a recent focus on the utility of modelling decay processes with the Lindblad equation [23,62,63].A consequence of using a density matrix based time evolution, however, is an increase in computational complexity.For the situation where only decay processes out of the Hilbert space under consideration are involved, an effective non-Hermitian Hamiltonian together with a Schrödinger equation can be used [64][65][66].However, when dephasing becomes relevant, such simplifications will be harder to make.
The quantum trajectories method (also called the stochastic Schrödinger equation [67] or Monte Carlo wave function [68] method) can be used to bypass the explicit evaluation of the density matrix.Here, we use this method to timeevolve the Lindblad equation using an ensemble of stochastic pure-state evolutions [68,69].Quantum trajectories are not yet widely used for polaritonic chemistry problems, but some studies have started to emerge [70].
In this work, we investigate the role of dephasing for molecular strong coupling and its effects on the dark state dynamics.We include the effects of both dephasing and photon decay.The model is based on the Lindblad equation and implemented with quantum trajectories.

II. SYSTEM AND MODEL
The model system is based on an optical cavity with a single carbon monoxide molecule, where the fundamental cavity mode is resonant with a transition between the 1 Σ + ground-state and an 1 Π electronically excited state in the CO molecule.To simulate a many particle system with dark states, which is also resonant with the cavity mode, N twolevel systems are added to the Hamiltonian.The Hamiltonian is thus an extension of the Tavis-Cummings model [4] (an optical cavity with N two-level systems under the rotating wave approximation), where the extension is the the CO molecule with its internuclear coordinate.
The Hamiltonian consists of the following five parts, where Ĥc corresponds to a single fundamental mode of linearly polarized light in an optical cavity, Ĥa are the N twolevel emitter systems (or atoms), Ĥca are all the resonant where â † and â are the photon creation and annihilation operators respectively, and ω c is the cavity mode frequency.The cavity photon energy is fixed at hω c = 8.27 eV (λ = 150 nm), which makes the transition between electronic states in the CO molecule resonant with the cavity.This choice puts the minima of the two potentials at similar energies, and prevents population from accumulating in either well, see Fig. 3.
where |g n and |e n are the ground and excited state of the nth two-level system, respectively.The operators σn = |g n e n | and σ † n = |e n g n | de-excite and excite the n-th two-level system, respectively.The energy of the transition is chosen to be resonant to the cavity's photon energy, hω a = hω c = 8.27 eV.
The interaction between emitters and cavity mode is derived under the dipole and the rotating wave approximations: The transition dipole moments of all two-level systems, µ a , are identical, and set to 1.5 Debye, which is on the same order of magnitude as the CO molecule (see Fig. 2

(b)).
The CO molecule is modelled by including its electronic 1 Σ ground-state potential energy curve, and one of the two electronically excited 1 Π states, which has a non-zero transition dipole µ(q) moment perpendicular to the molecular axis, and is assumed to be aligned with the field polarisation.The molecular Hamiltonian reads: Here, the first term is the kinetic energy of the nuclei, with the reduced mass of the CO molecule (M = 12 498 m e ).The potential energy curves for V Σ (q) and V Π (q) are shown in Fig. 2(a).
The previous choice of photon energy makes the molecular transition resonant with the cavity mode at an internuclear distance of q = 1.17 Å.
The transition dipole moment µ m ( q) varies between about 0 and 3 Debye according to Fig. 2(b).The vacuum electric field strength, E c (N) is the same as the atoms are experiencing (see Eq. ( 5)).
The vacuum electric field strength, E c (N), in Eqs. ( 5) and (7), is fixed at 3.00 V/nm for a single CO molecule model with no atoms (i.e.N = 0).To compensate for an increase in collective coupling strength, as we increase the number of atoms in the model, the vacuum field strength is scaled by 1/ √ N + 1 [48].
Note that the effects of single particle coupling strength can not be directly compared to a collective coupling strength of equal size.Thus, the results are expected to show a different behavior [24].The highest single particle cavity coupling strength, g avg = E c µ avg , occurs for N = 0. We compare this to the typical energy scale, g avg /hω c = 0.011.The system operates well below the ultra-strong coupling regime.Note that we assumed an average transition dipole moment of µ avg = 1.5 Debye.
In addition to the unitary part of the Hamiltonian, there are two non-unitary physical processes, which are introduced through the Lindblad equation (1).The Lindblad framework assumes that, that the system-bath correlations times are short enough to allow for the Markovian approximation [49,71].Thus, our results address interactions that can be modelled as non-Markovian.
The first non-unitary process is single-emitter dephasing of the CO molecule and all N two-level systems.The corresponding rate of dephasing, γ, is the key parameter whose effects we aim to investigate.It is the same for all emitters, which have their individual dephasing operator where, σ (n) z is the third Pauli matrix for the n-th emitter: Other choices of this operator yield the same time-evolution.However, this choice is beneficial for the quantum trajectory method, since it does not transfer population between states [72].
When constructing the Lindblad dephasing operators, we neglect the fact that the sub-systems are coupled (by the cavity mode) and build them phenomenologically.This is generally considered acceptable outside the ultra-strong coupling regime [73,74].In the appendix, section VIII A, we discuss the potential problems with this approach, and argue that the phenomenological operators are appropriate.

III. METHODS
Our studied observable is energy retention, i.e. we consider the fraction of the initial excitation that, despite the photon decay, remains in the system after 500 fs.The duration is chosen in relation to the timescale of the nuclear dynamics and mirrors previous investigations [24,62].The initial excitation consists of the cavity mode being is in its first excited (single photon) state, while all emitters are in their ground-states.This limits the total number of excitations to one, which allows us to truncate the basis.
A direct solution of the Lindblad equation ( 1) with a density matrix scales quadratically with the number of states involved.Such an approach becomes prohibitive as the number of atoms in our system increases.Instead of modelling a statistical state as a single density operator, we use quantum trajectories [68,69] to obtain the statistics from an ensemble of pure state wave functions.Thus, the cost is shifted from a single memory consuming density matrix, to running multiple, but lighter, pure-state calculations with wave functions, {ψ i (t)}.Each wave function evolves stochastically, with "quantum jumps" occurring randomly in proportion to physical parameters.One can prove that in the limit of an infinite ensemble of wave functions, the state from evolution with quantum trajectories approaches the state from the Lindblad equation [68].
From the ensemble of N T trajectories, {ψ i (t) : i ∈ 1 • • • N T }, a density matrix can be recovered by summing outer products of each wave function with a uniform weight: Expectation values { Âi } can also be obtained directly from the weighted sum of all trajectories, without having to construct ρ explicitly: The quantum trajectory method requires two modifications to the time-dependent Schrödinger equation.The first one adds the last term from Eq. ( 1), i.e. ∑ n −1/2 [ L † n Ln , ρ] + , to the Hamiltonian in the form of a norm-decaying term.However, in our choice of implementation, the wave function is continuously renormalized (at each discrete time-step).Using the Lindblad operators from Eqs. ( 9) and ( 11) leads to the non-Hermitian Hamiltonian: The second term in Eq. ( 14), is responsible for the photon decay, while the third term only affects the norm of the total wave function, since σ (n) The algorithm used for the propagation renormalizes the wave function in each step, thus the last term from Eq. ( 14) has no effect and can be removed.
The second part that is required for in the quantum trajectories method are discrete, stochastic jumps.These jumps originate from the first term in Eq. (1), i.e. ∑ n Ln ρ L † n .The probability of a jump occurring during a time-interval ∆t depends on the rates κ and γ as well as the population in the subspace from which the jump occurs: Note here that the Lindblad operators L k include the rates, κ and γ, as shown in Eqs. ( 8) and (11).At each time-step, random numbers are generated to determine if a jump occurs, in which case the Lindblad operator Lk , is applied to the wave function and the wave function is normalized.
We have implemented the quantum trajectories approach with our in-house software package QDng, which allows for the time evolution of wave functions.Implementations into existing methods are recurring in the literature [70,72,78].For each statistical state (see Eq. ( 12)), N T = 2500 wave function trajectories were run.Estimates of the resulting errors are given in the caption of Fig. 4.Each wave function was timeevolved with the Arnoldi propagation method [79], at order 10 and a time-step of to 0.5 au (0.012 fs).
The time-evolution is carried out in a product basis composed of the field-free molecular states, the field free two-level system states, and the Fock states of the cavity mode.In the following, we will refer to this basis as product basis.
For the purpose of interpreting the data (section IV) we diagonalize the potential curves, to obtain polaritonic curves FIG. 3. Potential energy surfaces of the polaritonic basis.A single CO molecule and N = 2 two-level systems in an optical cavity.With N = 2 the first dark state is introduced.The upper polaritonic state is yellow, the middle polaritonic state is green, the lower polaritonic state is pink, the ground-state is blue, and dark states are black.For higher values of N > 2 the curves are essentially the same [24], but additional degenerate dark states are introduced.
(see Fig. 3).This can be achieved by a pointwise diagonalization of the Hamiltonian along the reaction coordinate, q, while omitting the kinetic energy operator.The number of polaritonic energy surfaces depends on N, but for N > 2 all additional energy surfaces are the degenerate dark states.Thus, plotting the result for N = 0 and N = 2 gives an overview for all other values of N (see Fig. 3(b)).Note that dark states are a feature of the polaritonic basis, and does not appear in the product basis.We ensure that uncoupled polaritonic surfaces are allowed to cross, by following the eigenvector tracking method described in [24].
The potential energy curves and dipole functions of the CO molecule, shown in Fig. 2, are the relevant results of a quantum chemistry calculation that originally included eight non-degenerate states and their respective transition dipole moments.The electronic structure calculations were carried out with the program package Molpro [80][81][82] at the CASSCF(10/14)/MRCI/aug-cc-pVQZ level of theory, with a state average over a total of twelve electronic states.Energies, dipole moments, and transition dipole moments are calculated at 50 internuclear distances, between 0.926 Å and 6.35 Å.The two electronic states included in this work are 1 Σ + and one of the two doubly degenerate 1 Π states.The data is in good agreement with previous calculations [83,84].The molecular potentials and transitions moments are interpolated to a spatial grid with 96 grid points in the interval 0.90 ≤ q ≤ 2.12 Å (see vertical lines in Fig. 2).The result is used to construct the molecular Hamiltonians, Ĥm and Ĥcm in Eqs. ( 6) and (7), for the numerical calculations.3(a).Colors match the ones used for the surfaces in Fig. 3.In the left column (a) the dephasing is slow (γ ≈ 0.002 fs −1 ), in the right column (b) dephasing is fast (γ ≈ 0.09 fs −1 ).The left grey arrows in Fig. 4 shows the data point of (a), and the right grey arrow is the data point for (b).The two cases, (a) and (b), behave very similar, the most apparent difference is the dampening of interference for fast dephasing, which makes the curves in (b) smoother.

IV. RESULTS AND DISCUSSION
To construct a molecular Tavis-Cummings model, which includes dark states, N two-level systems were introduced along with the CO molecule.We included the values N ∈ {0, 1, 2, 3, 5, 8, 13, 22, 36, 60} and for each value the dephasing rate is varied between 0 and 10% of the cavity frequency ω c , which corresponds to the dephasing rates: 0 ≤ γ ≤ 1.26 fs −1 .After 500 fs of time evolution, the energy retention is recorded.The remaining energy in the system is then plotted as a fraction of the initial energy, in relation to the dephasing rate for each N.This constitutes the main result in this study and is shown in Fig. 4. Note that the obtained curves are not perfectly smooth, which is an expected result of the stochastic sampling of the density matrix.
Only 10% to 15% of the initial energy is retained in the system for the dephasing free case (γ = 0, see crosses on the vertical axis in Fig. 4).With no dephasing, the N −1 dark states are not populated and does not impact the behavior of the system.However, even though the number of bright states are fixed, and with a constant collective coupling strength (see Eq. ( 8)), the energy retention for γ = 0 varies.This can explained as follows: the collective Rabi splitting appears to behave similar to the single particle Rabi splitting between the upper and lower polariton.However, the splitting between middle polariton state and the upper and lower polariton state follows the single particle Rabi splitting and not the collective Rabi splitting [24].
With no two-level systems, i.e.N = 0, the impact of dephasing on the energy retention is small (black curve in Fig. 4).Fig. 5(a) and (b), compare the populations in the polaritonic basis (from Fig. 3(a)) for a slow and fast dephasing rate (and N = 0).The most obvious difference in the populations is the amount of oscillations caused by interference and by the avoided crossing in polaritonic states.We can understand the dampened oscillations as dephasing canceling the otherwise coherent population transfer between different trajectories.
Introducing dark states (N ≥ 2) will make the energy retention increasingly sensitive to changes in the dephasing rate.We consider first the case with a single dark state (N = 2).
Here, dephasing has an observable effect on energy retention.Fig. 3(b) shows the polaritonic energy surfaces for N = 2, which includes a single dark state (black curve) and middle polariton state (green curve).The corresponding population evolution is shown in Fig. 6, for slow dephasing (a) and fast dephasing (b).Note, how a fast dephasing rate imposes a rapid build-up of population in the dark state, which peaks at about 20 fs and slowly decays thereafter.This in contrast to the slow FIG. 6. Population time-evolution of the polaritonic potential energy surfaces shown in Fig. 3(b).Colors match the ones used for the surfaces in Fig. 3.In the left column (a) the dephasing is slow (γ ≈ 0.002 fs −1 ), in the right column (b) dephasing is fast (γ ≈ 0.09 fs −1 ).The left blue arrows in Fig. 4 shows the data point of (a), and the blue arrow is the data point for (b).
dephasing rate, where the dark state population occurs slowly over a timescale of ≈ 300 fs.
The highest energy retention is observed for the largest number of two-level system studied (N = 60) and a dephasing rate on the order of γ ≈ 10 −1 fs −1 .Here, the solid yellow curve in Fig. 4 shows a sharp increase in energy retention with increasing dephasing, from around 10 % to almost 90 %.The time-evolution of the populations for N = 60 is shown in Fig. 7, for both slow dephasing in (a) and fast dephasing in (b).The populations of all 59 dark states are shown as a sum (black curve).For the fast dephasing rate, the dark states population builds up rapidly, and reaches its maximum value of 90 % in about 40 fs.Thereafter, the dark states population is very slowly released back to the bright states, with a retention of 80 % at 500 fs.For N = 60 and a slow dephasing rate, the build up to about 27 % is slow, and does not reach a maximum within 500 fs.For both slow and fast dephasing, the population in the dark states has significantly increased when compared to the case N = 2 (compare Figs. 6 and 7).
The mechanism for the energy retention can be explained by the transition dipole moments of the collective system.The main energy loss mechanism is the decay of photons due to imperfect cavity mirrors.However, the dark states have no transition dipole moment in the absence of dephasing, as they represent asymmetric superpositions of matter excitations.Thus, the dark states do not couple to the cavity mode in an ideal system and can not scatter photons into the cavity mode.Population in the dark states are thus protected from photon decay.However, the dephasing breaks the symmetry in the Hamiltonian that decouples the dark states from the upper, middle, and lower, polariton states and enables population transfer.With an increasing number of particles N, the number of available dark states increases, and thus the effective rate increases with which these states are populated increases.The dark state subspace can thus serve as a reservoir which protects the systems from photon decay [49,85].
The population transfer between dark states and bright states goes both ways.As can be seen in Figs.6(b) and 7(b), dephasing will also slowly return the population to the bright states.However, this process is slower and than the population transfer into the dark states, resulting in an overall slow down of the energy loss.
The energy retention in Fig. 4 has a local maximum with respect to the dephasing rate at about γ ≈ 10 −1 fs −1 and N = 60.In this regime, the Rabi oscillations (Ω R = 2.85 × 10 −1 fs −1 ) are dampened, but not yet over dampened.The line width of the dark states and the polariton states are sufficiently narrow to be spectrally well separated.
However, if the dephasing rate the is further increased, the polariton states begin to overlap with the dark states.In this over dampened regime, the dark states are no longer sufficiently decoupled from the bright states, and there is a significant leakage from the dark states.This effect causes a decreased energy retention for large dephasing rates.Note that the maximum of the energy retention shifts to larger dephasing rates with decreasing N. Note that the decrease in energy retention, after the maximum, depends on the choice of our initial state (excited cavity), through the momentary rates populating and depopulating the dark state reservoir.

V. CONCLUSION
In summary, we have investigated a molecular Tavis-Cummings model under the influence of dephasing and photon decay.The model system consisted of a CO molecule with a varying number of resonantly coupled two-level systems, which has been solved with a quantum trajectory approach instead of a direct solution of the Lindblad equation.
This approach allowed us to use a wave function base calculation which scales more favorably with respect to the number of states than evaluating a density matrix explicitly.
In this atomistic model, we could show that the dark states become increasingly coupled to the polariton states as the dephasing rate is increased.As a result, population gets trapped in the dark states, and it is protected from photon decay processes.Our findings are in line with earlier studies [49,85].Under the influence of dephasing the dark states are no longer decoupled from the polariton states.The effective transfer rate into the dark states even increases with an increasing number of dark states, and thus providing an increasingly efficient protection against photon decay.In the investigated systems, the photonic excitation was rapidly transferred into the dark state reservoir and slowly released back into the bright polariton states.
Our results show that dephasing, as it would occur under experimental conditions in condensed phase, plays a significant role in the dynamics of such a system.Realistic models should thus not only include photon decay, which has been demonstrated to play a crucial role [24], but should also include the effects from pure dephasing in condensed phase.

VII. DATA AVAILABILITY
The data that support the findings of this study are available upon reasonable request.

VIII. APPENDIX A. Impacts from phenomenological dephasing operators
For a single two-level system, the dephasing operator in Eq. ( 9) has the expected effect; to decay away phase information without moving population between states.However, when a system is built from many strongly coupled sub-systems (such as in our case), this operator will include some unintended amount of driving and decay.To probe how significant this effect is, we can fully diagonalise the Hamiltonian from Eq. ( 2), and apply the basis transformation to a dephasing operator.(Note that we can omit the ground-state from the diagonalisation since it does not couple to any other state.)This basis transformation into a fully diagonalised Hamiltonian will give a set of eigenstates containing some unphysical states; when the energy of the vibrational motion is enough to dissociate the molecule.However, since our time-evolution does not exhibit dissociation we do not populate such questionable states, and they could be excluded from this analysis without affecting the conclusion.
A dephasing operator (for N = 2 atoms) after transformation to the energy ordered eigenbasis is shown in Fig. 8. Driving and decay will appear as off-diagonal elements.
On close inspection we can see that there are off-diagonal elements, but that the matrix is essentially block diagonal.This means that there are some driving and decay between eigenstates that are close to degeneracy, but the population is trapped within each block of states.Thus, we cannot drive or decay population further than between a few energetically adjacent states.This is not of concern to our investigation and a phenomenological operator is considered sufficient.For further discussions, see for instance [49,86,87].

FIG. 1 .
FIG. 1. Schematic representation of the Tavis-Cummings model.The upper |UP and lower polartion state |LP are split by a collective Rabi frequency.The N − 1 dark states |DS do not couple to the electromagnetic field and are thus unaffected by the cavity.Line broadening, cause by dephasing for example, introduces a coupling of the dark states with the polariton states.

FIG. 2 .
FIG. 2. Potential energy curves and transition dipole curves for the CO molecule.Dots marks frozen internuclear distances from the quantum chemistry calculation.Vertical lines surround the region used in the time-evolution.(a) Potential energy surfaces for the 1 Σ ground-state and one 1 Π excited state.(b) Transition dipole moment between 1 Σ and 1 Π.

FIG. 4 .
FIG. 4. Energy retention in relation to dephasing rate.The vertical axis shows energy retention, i.e. the remaining energy as a fraction of initial energy, E final /E initial where E initial = hω c , after 500 fs of time-evolution.The horizontal axis is the single emitter dephasing rate, γ, which is the same for all emitters.Each curve corresponds to a particular number of two-level systems, N. The number of dark states for each curve is N − 1. Crosses on the vertical axis shows the energy retention for γ = 0. Arrows mark points where dynamics plots are supplied, see Figs. 5, 6, and 7. Standard deviations due to the Trajectories method (from a set of 25 runs) are also calculated at these points.Grey arrows both have σ ≈ 0.008, left blue has σ ≈ 0.011, right blue has σ ≈ 0.012, left yellow has σ ≈ 0.009, right yellow has σ ≈ 0.008.

FIG. 7 .
FIG. 7. Population time-evolution of the polaritonic potential energy surfaces for N = 60 (which are very similar to the ones shown in Fig. 3. Colors match the ones used for the surfaces in Fig. 3(b), but the black curve is a sum of all 59 dark states.In the left column (a) the dephasing is slow (γ ≈ 0.002 fs −1 ), in the right column (b) dephasing is fast (γ ≈ 0.09 fs −1 ).The left yellow arrows in Fig. 4 shows the data point of (a), and the right yellow arrow is the data point for (b).
VI. ACKNOWLEDGMENTS This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 852286) and European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement no.E.D. thanks Prof. Jonas Larson for valuable discussions and insightful feedback.