Simulating photodissociation reactions in bad cavities with the Lindblad equation

Optical cavities, e.g., as used in organic polariton experiments, often employ low finesse mirrors or plasmonic structures. The photon lifetime in these setups is comparable to the timescale of the nuclear dynamics governing the photochemistry. This highlights the need for including the effect of dissipation in the molecular simulations. In this study, we perform wave packet dynamics with the Lindblad master equation to study the effect of a finite photon lifetime on the dissociation of the MgH+ molecule model system. Photon lifetimes of several different orders of magnitude are considered to encompass an ample range of effects inherent to lossy cavities.

Light-matter interactions are modeled under the dipole and rotating wave approximations. Molecular transitions resonant with the cavity mode frequency are |X⟩ ↔ |A⟩, |X⟩ ↔ |B⟩, and |A⟩ ↔ |C⟩ (see Fig. 1 for their potential energy surfaces), Hcm = EcμXA(q)(â † |X⟩⟨A|+â|A⟩⟨X|)+EcμXB(q)(â † |X⟩⟨B|+â|B⟩⟨X|) The transitions |A⟩ ↔ |B⟩, |B⟩ ↔ |C⟩, and |X⟩ ↔ |C⟩ are not resonant with the cavity mode frequency, but the transition |X⟩ ↔ |C⟩ is used for creating the initial excited state. The transition dipole moments {μ(q)} are obtained at the same level of theory as the potential energy surfaces, and Fig. 8 shows the ones that are used in the model. For a minimal model of the molecule, permanent dipole moments and self-energy are omitted. The vacuum electric field strength, Ec, determines the strength of the light-matter interaction and depends on the mode volume of the cavity structure, V,

APPENDIX B: SUMMING LINDBLAD OPERATORS AHEAD OF TIME EVOLUTION
To do time evolution with the Lindblad equation (1), all {Ln} operators must be stored in memory, and given N such operators, 6N matrix multiplications are required for each time step. A method for reducing memory requirements and computational cost is demonstrated here. Memory storage is reduced from N matrices to 2. The computational cost goes from 6N matrix multiplications to four matrix multiplications and two operations where off-diagonal elements are set to zero. Three assumptions are used, which are fulfilled in this study.
Assumption 1. In the basis employed for the computation, each Lindblad operator is decaying a single basis state to another basis state. ◽ Assumption 2. There are no Lindblad pure de-phasing operators in the model. ◽ Assumption 3. Lindblad operators and decay rates are timeindependent. ◽ If assumptions 1 and 2 do not hold, this method may still be useful in combination with a change of basis and/or treating de-phasing operators separately.
There are two useful notations for the set of all Lindblad operators. The most minimal notation is to simply enumerate them, However, according to assumption 1, each operator involves only two states from the employed basis {|ei⟩}. Thus, Lindblad operators can also be expressed in terms of the basis states that comprise them, The associated decay rates are denoted as {κn} or {κij}, respectively. Depending on whether the emphasis is on the number of operators or on the basis states that comprise them, these two notations will be used interchangeably in the following. The focus is on treating the Lindblad operators inside the sum of the Lindblad equation, These two terms in Eq. (B3) are performing the functions of re-populating states, and depopulating states, in that order of appearance.

ARTICLE scitation.org/journal/jcp
The first term, responsible for re-populating states, can unfortunately not be factorized, However, we treat it as if it could be anyway and define the left parenthesis from Eq. (B4) asŜ 1 , The multiplication from Eq. (B4) is then expanded, Terms with a single κn are identical to the desired terms from the Lindblad equation (B3), but there are also several unwanted crossterms. Expressingρ in the employed basis {|ei⟩} and using the result that all operators {Ln} ≡ {Lij} are of the form of Eq. (B2), we simplify the operator for each type of term, beginning with the single κn term,L These desired terms thus correspond to increasing population in diagonal matrix elements ofρ (due to |ei⟩⟨ei|), in proportion to other diagonal matrix elements ofρ (due to ρjj). Then, the cross-terms are simplified aŝ The {Lij} operators are required to be orthogonal, which means that either i ≠ m or j ≠ n. That is, either |ei⟩⟨em| is an off-diagonal matrix element after matrix multiplication or ρjn is an off-diagonal matrix element before matrix multiplication (or indeed both). Setting off-diagonal elements ofρ to zero both before and after matrix multiplication withŜ 1 will therefore remove only the crossterms from Eq. (B6). Let D[⋅] be this operation, which sets offdiagonal matrix elements to zero in the employed basis. The initial sum over multiple matrix multiplications can now be rewritten usingŜ 1 , In contrast to the re-populating term from the Lindblad equation (B3), the other term, responsible for depopulation, can be factorized, Defining the left parenthesis above as the operatorŜ 2 , we can reduce this to two matrix multiplications, With the definition ofŜ 1 andŜ 2 from Eqs. (B5) and (B11), they are calculated once ahead of time evolution and stored in memory, and the set {Ln} can be discarded. Together with the operation D[⋅], the initial sum from Eq. (B3) with 6N matrix multiplications is reduced to four matrix multiplications and two D-operations where off-diagonal elements are set to zero,

APPENDIX C: NUMERICAL METHOD AND VERIFICATION
The Lindblad equation is solved numerically, using the ode45 36 (Runge-Kutta) differential equation solver, as implemented in Octave, 37 hereafter referred to as the Octave method. Verification and accuracy benchmarks of the Octave method are done by comparison to time evolution with the Schrödinger equation, in the in-house QDng package, using the Chebyshev propagator, 38 hereafter referred to as the QDng method. This QDng method uses 512 grid points and a fixed duration time step set at 0.0242 fs, previously determined to have high accuracy in similar systems when compared to even shorter time steps. 9 The Octave method uses 96 grid points and a variable duration time step. Accuracy is here specified with a 1 × 10 −6 absolute tolerance and 1 × 10 −3 relative tolerance. For both methods, we calculate the remaining population after 500 fs of time evolution (i.e., the same data that make up our main result in Fig. 3). The QDng method requires there to be no photon decay; thus, κ = 0 in the Octave method. A comparison of the two methods is shown in Fig. 9.
At a standard deviation in the remaining population of 0.0086, the agreement between the two methods is satisfactory. Noticeable deviations occur for higher field strengths, where the impact of interference effects causes some disagreement. However, the overall profile of the oscillations is still accurately tracked. Details from time evolution, for three data points in the parameter range, are shown in Figs. 4, 5, and 6. The first data point is in the polaritonic strong coupling regime, with parameters Ec = 3.0 GV/m and τ = 3.6 × 10 4 fs, directly after the sudden rise in sector (b) of Fig. 3. The second data point has parameters Ec = 3.0 GV/m and τ = 48 fs, which puts it in the low stability region of sector (e). The third data point has parameters Ec = 3.0 GV/m and τ = 0.58 fs, which is a data point on top of the high stability region in sector (g) of Fig. 3. At each data point, populations are plotted individually for each product state (see Fig. 2 for the states). States that are not noticeably populated (peak population less than 0.08) are not shown. The total population remaining in the system is also shown, along with a renormalized purity. The renormalization is done to eliminate loss of purity caused directly from the loss of norm. The full time evolution is 500 fs, but plots are limited to the first 250 fs where the influential processes are taking place.