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.

Motivated by recent experimental advancements,1–3 state-of-the-art computational tools are presently employed to consider the interaction between quantized electromagnetic radiation and a diverse range of matter systems in models reminiscent of the early Jaynes–Cummings model.4 Examples of such matter systems include quantum dots,5 excitons,6 ensemble systems,7–9 Bose–Einstein condensates,10 superconductors,11 and molecules of different sizes and complexity.12 

Central to these studies, be it experimental or theoretical, is the macro- or nanoscopic material structure giving rise to electromagnetic field-modes in a confined space, henceforth called a cavity structure, though they are not always reminiscent of the widely used Fabry–Pérot cavity. The material and geometry of these cavity structures determine the characteristics of the light–matter interaction, and precise tailoring enables a unique possibility to probe or control phenomena that take place on the scale of single particles. Examples of such phenomena include chemical reactions,1 photo-dissociation,9 electron transport,13,14 isomerization,15 and heat-transfer between molecules.16 

In theoretical works such as this, the properties of the cavity structure are typically translated into parametric values, such as the mode frequency and the vacuum electric field strength,17,18 entering the model of the system. The influence of these parameters is well studied and generally well understood.

However, not all parameters describing the cavity structure have gained equal attention. Effects arising from a finite lifetime for field excitations, also known as the cavity Q-factor or photon decay rate, have only recently been the main focus in a few ground-breaking studies,19–21 and dissipative effects, in general, have been described as a challenge for computational methods.1 

Theoretical studies often assume an infinite lifetime,15 which can be appropriate but introduces inaccuracies when used to model conditions where finite photon lifetimes are de facto non-negligible. For instance, this is typically the case for plasmonic nano-cavity structures,11,22 such as ultra-thin metallic gaps23 or nano-gap antennas,24,25 where the desire to concentrate a strong electromagnetic field in a small volume inevitably leads to short-lived excitations of the field.1 

Another previously employed strategy is to incorporate a Hamiltonian coupling between affected states and a reservoir of simple systems.26 However, this allows the energy to return from the reservoir and affect the evolution of the system, which therefore makes an approximative model for electromagnetic energy being lost into the environment. Additionally, photon decay processes inevitably introduce decoherence into the system, which makes an approach with a pure state wave function insufficient.27 Instead, a density matrix formalism is required for the time evolution of these statistically mixed states, and the Schrödinger equation is generalized to some variety of a master equation. Here, we use the Markovian Lindblad master equation, where the decay rate of field excitations will enter through the cavity decay rate κ,27,28

(1)

However, using the density matrix, ρ^, increases the computational cost considerably compared to a wave function based approach.

Recent publications19–21,29,30 on cavity decay do not deal directly with the Lindblad equation formalism; instead, the underlying formalism is a non-Hermitian evolution of a pure state, which is shown to be appropriate for both relaxation dynamics19,21 and isomerization.20 However, for general applications, the method omits some phenomenological processes, such as decoherence and time evolution of statistically mixed states, and the impact of this exclusion depends on the type of system and is observable under investigation.

In this paper, we investigate the molecular dissociation of the MgH+ molecule after excitation to an unbound electronic state and in the presence of a lossy cavity structure. Here, the photon decay will populate intermediate states that also contribute to the dissociation. Thus, a wave function based approach is not suitable for this type of problem (see the discussion at the end of Sec. II). Instead, the Lindblad master equation is used to capture the behavior of the system and deliver accurate results. We investigate how the photon lifetime and cavity vacuum field strength affect the photostability of the MgH+ molecule, and the photochemical reaction mechanisms of the coupled light–matter system are analyzed.

The Hamiltonian in Eq. (1) is a molecular Jaynes–Cummings type Hamiltonian9 modeling vibrationally and electronically excited states in the presence of a lossy cavity mode. It is composed of Ĥm for the molecule, Ĥc for the cavity mode, and the light–matter interaction Ĥcm,

(2)

We assume the cavity Born–Oppenheimer approximation,31 the rotating wave approximation, and the dipole approximation for a spatially fixed molecule. Excited electronic and vibrational states of the molecule are described on one-dimensional potential energy surfaces. The four lowest electronically excited states of MgH+ are considered (see Fig. 1). Wave-packets approaching the dissociation limit are absorbed by an imaginary potential. The cavity mode is modeled as a single mode with a photon energy of 4.3 eV (285 nm). Comprehensive details about the Hamiltonian are found in  Appendix A.

FIG. 1.

The four lowest electronic states in MgH+. Each potential energy surface is implemented with an absorbing potential, whose Gaussian shape and relative size are shown with the dotted line on the black curve.

FIG. 1.

The four lowest electronic states in MgH+. Each potential energy surface is implemented with an absorbing potential, whose Gaussian shape and relative size are shown with the dotted line on the black curve.

Close modal

The states of the combined molecule-cavity system are expressed as a product state, n,M=nM, where n{0,1,2,} are the Fock states of the cavity mode and M{X, A, B, C} are the electronic states of MgH+. All product states that have distinctly higher energy than the initial state, 0,C, will never be populated under the rotating wave approximation and are removed from the description. The full Hamiltonian in Eq. (2) is then expressed in the basis of eight states {0,X, 1,X, 0,A, 0,B, 2,X, 1,A, 1,B, 0,C} covering an energy range of ∼10 eV.

The product states are grouped into three subspaces: the ground state subspace, the single excitation subspace, and the double excitation subspace. The corresponding potential energy curves are shown in Fig. 2. This partition will later be used for analyzing the dynamics of the system. Under the rotating wave approximation, states from different excitation subspaces are not coupled, which simplifies the Hamiltonian.32 Only states within the same excitation subspaces are coupled and exhibit curve crossings, as shown in Fig. 2.

FIG. 2.

Potential energy surfaces associated with the basis states for electronic and cavity degrees of freedom. The initial state is the pure molecular excitation, 0,C. Arrows indicate which states are coupled by Lindblad decay operators. States shown as dashed curves are referred to as the double excitation subspace. States shown as solid curves are referred to as the single excitation subspace. The state shown as a dotted curve is referred to as the ground-state subspace. The two lower subspaces (ground state and single excitation) are occasionally considered as a group and referred to as the decohering subspaces.

FIG. 2.

Potential energy surfaces associated with the basis states for electronic and cavity degrees of freedom. The initial state is the pure molecular excitation, 0,C. Arrows indicate which states are coupled by Lindblad decay operators. States shown as dashed curves are referred to as the double excitation subspace. States shown as solid curves are referred to as the single excitation subspace. The state shown as a dotted curve is referred to as the ground-state subspace. The two lower subspaces (ground state and single excitation) are occasionally considered as a group and referred to as the decohering subspaces.

Close modal

The photon decay, as modeled by the Markovian Lindblad master equation, assumes that the photon is lost irrevocably from the system28 (either into free space or by absorption at the surrounding bulk material).

In the product basis, these unidirectional interactions are identified by, wherever possible, decreasing the photon number by one. This yields four transitions: 1,B0,B, 1,A0,A, 2,X1,X, and 1,X0,X, which are indicated by vertical arrows in Fig. 2. Each of these unidirectional interactions forms a non-Hermitian Lindblad jump operator. The decay rates {κn} from the general case in Eq. (1) are here the same value, κ, for all photons in the mode, and κ is determined for any particular cavity structure. Since the state 2,X has two photons, the associated annihilation operator will introduce a factor of 2 in the Lindblad operator, effectively increasing the decay rate proportional to the energy increase in the field,

(3)

This model includes no Lindblad de-phasing operators, such as the effects derived from weak interactions between the molecule and its environment, which are assumed to be negligible on the timescale for molecular dissociation of MgH+.

The initial state for the simulation, ρ^0, is formed from the pure state 0,C,ψ(q), where the component ψ(q) is the nuclear vibrational wave function, created by weighting the ground-state wave function of X with the transition dipole moment between XC, which corresponds to a vertical excitation,

(4)

This method provides a consistent initial condition of a fully excited molecule without the introduction of additional parameters. With Ĥ, {L^n}, κ, and ρ^0, the model is fully defined, and the time evolution of Eq. (1) can be simulated as described in Sec. III.

As a general remark, to minimize computational cost, the Lindblad master equation can be reduced to a non-Hermitian Schrödinger equation employing an absorbing potential, given the following three conditions: The Hamiltonian does not couple any state in the subspace to other states that are not in the subspace, the initial state projected onto that subspace is a pure state, and the subspace is not the recipient of decaying states. This motivates the method in previous studies,19–21 but in this study, these conditions only apply to the doubly excited subspace (see Fig. 2 for definitions of subspaces). The main issue is that a straightforward reduction to a Schrödinger equation only accounts for the removal of population, and not the re-population that happens in the single excitation subspace, and ground-state subspace. The re-population is necessary for this study since a dominant contribution to the observable occurs after it. Additionally, when the decay from the Lindblad operators is transferring population between subspaces, the state is decohering and re-populating the receiving subspaces as a statistical mixture of states. These mixed states display a suppressed interference when evolving on the potential energy surfaces, an effect that is further discussed in relation to our data in Sec. IV. For these reasons, time evolution of a non-Hermitian Schrödinger equation, which only accommodates pure states, is not well suited to our system. Instead, the Lindblad equation will describe these influential effects.

The potential energy curves (see Fig. 1) are calculated with the program package Molpro33 at the CASSCF(12/10)/MRCI/ROOS level of theory.34,35 The time evolution of the Lindblad equation is done numerically with a Runge–Kutta scheme as it is implemented in the differential equation solver ode4536 in Octave.37 The density operator ρ^ is represented on a numerical grid for the nuclear coordinate q with 96 grid points for each of the included states n,M. The density matrix is propagated for 500 fs, a duration selected to reflect the relevant timescale of the MgH+ dissociation. Due to the increased computational cost associated with the Lindblad equation, a strategy was developed where Lindblad operators could be summed ahead of time evolution, reducing both memory requirements and computational cost (see the derivation and details about this strategy in  Appendix B).

The accuracy of the method is tested against time evolution with our in-house software package QDng using the Chebyshev propagation method,38 and good agreement between the two was found (see  Appendix C for benchmarking results and further details about the implementation of the numerical method).

Initially, the system is vertically excited in full to the state 0,C, which corresponds to a dissociative state of the MgH+ molecule. Since this state has no field excitations, the initial state will not itself decay to any other state. The vacuum electric field strength Ec—which scales the light–matter interaction strength according to Eq. (A4)—is sampled for a range of values: Ec[0,6] GV/m. After 500 fs of time evolution, the remaining population in the system (the trace of ρ^) is recorded. This is the portion that has not been removed by absorbing potentials, and it measures the stability of the MgH+ molecule over the range of light–matter interaction strengths. The mean lifetime of the field excitation, τ = 1/κ, is then varied over eight orders of magnitude.

Figure 3 shows the population data from a batch of calculations, plotted on a two-dimensional grid. The photon lifetime τ is divided into eight sectors, one for each order of magnitude, from 105 fs in sector (a) to 10−2 fs in sector (h). For comparison, two relevant times are marked with one white solid line and one white dashed line. The solid line identifies the total duration of time evolution, which corresponds to 500 fs. The dashed line at 0.5 fs identifies the time it takes light to travel across a Fabry–Pérot type cavity with length λ/2 ≈ 140 nm. For such cavities, shorter lifetimes are not physical.

FIG. 3.

The remaining population on the vertical axis is the trace of ρ^, 500 fs of time evolution after the initial excitation to a dissociative state, and it measures the stability of the MgH+ molecule. The electric field strength, Ec, determines the strength of the light–matter coupling in the interaction Hamiltonian (A3). Mean lifetime, τ, refers to that of the field excitation from the cavity structure, and τ = 1/κ, where κ is the decay rate from the Lindblad equation (1). Black thicker lines partition the lifetime parameter into eight sectors [(a)–(h)], each corresponding to one order of magnitude. For reference, the white solid line in sector (d) marks the duration for time evolution (500 fs). The white dashed line in sector (g) marks the time duration for light to cross the length of a Fabry–Pérot type cavity (0.5 fs). The white dotted line crossing several sectors marks the points where the Rabi splitting is roughly equal to the full width at half maximum broadening of the field excitation, ΩR(q) ≈ Γ. On the left-hand side of this line resides the strong coupling regime where ΩR(q) ≫ Γ.

FIG. 3.

The remaining population on the vertical axis is the trace of ρ^, 500 fs of time evolution after the initial excitation to a dissociative state, and it measures the stability of the MgH+ molecule. The electric field strength, Ec, determines the strength of the light–matter coupling in the interaction Hamiltonian (A3). Mean lifetime, τ, refers to that of the field excitation from the cavity structure, and τ = 1/κ, where κ is the decay rate from the Lindblad equation (1). Black thicker lines partition the lifetime parameter into eight sectors [(a)–(h)], each corresponding to one order of magnitude. For reference, the white solid line in sector (d) marks the duration for time evolution (500 fs). The white dashed line in sector (g) marks the time duration for light to cross the length of a Fabry–Pérot type cavity (0.5 fs). The white dotted line crossing several sectors marks the points where the Rabi splitting is roughly equal to the full width at half maximum broadening of the field excitation, ΩR(q) ≈ Γ. On the left-hand side of this line resides the strong coupling regime where ΩR(q) ≫ Γ.

Close modal

The features calling for explanations in Fig. 3 occur exclusively at higher electric field strengths. For small values, on the other hand, Ec1 GV/m, the initial molecular excitation into the dissociative, but non-decaying, state 0,C does not exchange enough population with other dipole-coupled states, and essentially all the population is absorbed during the first 50 fs. This is the expected behavior of a free, dissociating molecule. The following discussion will therefore only consider the behavior of a molecule coupled to the field mode, i.e., Ec2 GV/m.

In sectors (a) and (b) (Fig. 3), where lifetimes are long (τ on the order of 105 fs or 104 fs), the molecule is significantly stabilized. This parameter regime can be safely considered as strong coupling. The mean lifetime here is long enough for decay processes to be negligible, and the system behaves as if the lifetime was infinite (or κ = 0). In these sectors, a sharp rise in molecular stability can be observed as the electric field strengths go beyond 2 GV/m. The stability then plateaus and oscillates around ∼0.7. The cause for the rise in stability can be attributed to the growing Rabi-splitting due to the stronger light–matter coupling.39,40 This increases the energy difference between the polaritonic states, which suppresses the transfer of population between them, and stabilizes the molecule. The observed small-scale oscillations are understood as a consequence of interference effects between nuclear wave packets in the crossing regions.41 The low rate of decay in these sectors retains the population in the double excitation subspace and the state does not decohere. Time evolution data from sector (b) is shown in Fig. 4.

FIG. 4.

Time evolution data for Ec=3.0 GV/m and τ = 3.6 × 104 fs [sector (b) in Fig. 3]. States not shown have a peak population of less than 0.002. (a) Populated states in the double excitation subspace. (b) Solid line shows the total population, and the dashed line shows the renormalized purity.

FIG. 4.

Time evolution data for Ec=3.0 GV/m and τ = 3.6 × 104 fs [sector (b) in Fig. 3]. States not shown have a peak population of less than 0.002. (a) Populated states in the double excitation subspace. (b) Solid line shows the total population, and the dashed line shows the renormalized purity.

Close modal

In sector (c), where the mean lifetime is measured in thousands of femtoseconds (τ on the order of 103 fs), the impact of photon decay starts to become noticeable. However, when compared to the infinite lifetime case, the sector is still qualitatively similar to a varying Ec. The error introduced by an infinite lifetime approximation, in sector (c), can be quantified by population deviations from the case of τ = . To get the worst-case error in the sector for the population, the lifetime is fixed at its shortest, τ = 1000 fs, and the deviation from τ = is calculated for all such points, which gives an average deviation of ∼0.07.

The analysis is repeated for a lifetime one order of magnitude larger than the timescale of the studied phenomena. Here, 500 fs of time evolution gives τ = 5000 fs. The average deviation is then calculated to ∼0.02, which can be an acceptable deviation for many applications.

In sector (d), with mean lifetimes on the order of hundreds of femtoseconds (τ on the order of 100 fs), the nuclear dynamics is now on the same timescale as the photon decay and cannot be neglected anymore. A small change in the lifetime now changes the stability of the molecule. In sectors (a)–(c), the long lifetime implies that dissociation is primarily happening in the double excitation subspace but this changes in sector (d). The decay rate is now sufficient for the system to also dissociate from the single excitation and ground-state subspaces via the states 0,A, 0,B, 1,X, and 0,X. For Ec2.5 GV/m, this type of dissociation dominates in sector (d). The loss of coherence in the state ρ^(t) results in a suppression of the interference effects and the finer oscillations are dampened. The trapping in polariton states that has contributed to the stabilization now becomes less efficient. The photon decay projects the wave packet partially onto unbound vibrational eigenstates in the single excitation subspace, as well as the ground state. The population remaining in these unbound states for long enough will then contribute to dissociation.

In sector (e), where the lifetime is on the order of tens of femtoseconds (τ on the order of 101 fs), the local minimum in Fig. 3 can be understood as an optimum for dissociation via both the single excitation subspace and the ground state subspace. All subspaces are populated just long enough such that each subspace can contribute to the dissociation. In this sector, the behavior of the system is dominated by effects from decay, and we are now firmly in the dissipative regime (see Fig. 5 for population data from time evolution in the product states).

FIG. 5.

Time evolution data for Ec=3.0 GV/m and τ = 48 f [sector (e) in Fig. 3]. All states are noticeably populated. (a) Populated states in the double excitation subspace. (b) Populated states in the single excitation subspace. (c) Populated states in the ground-state subspace. (d) Solid line shows the total population, and the dashed line shows the renormalized purity.

FIG. 5.

Time evolution data for Ec=3.0 GV/m and τ = 48 f [sector (e) in Fig. 3]. All states are noticeably populated. (a) Populated states in the double excitation subspace. (b) Populated states in the single excitation subspace. (c) Populated states in the ground-state subspace. (d) Solid line shows the total population, and the dashed line shows the renormalized purity.

Close modal

In sectors (f) and (g), where the photon lifetime is on the order of femtoseconds or tens of femtoseconds (τ on the order of 100 fs or 10−1 fs), the lifetime is short enough that the population decays into the ground state before there is enough time to contribute to dissociation from the other subspaces. The time evolution is shown in Fig. 6, where the initial population in 0,C decays almost instantaneously to the mostly bound vibrational states in 0,X (via 0,A).

FIG. 6.

Time evolution data for Ec=3.0 GV/m and τ = 0.58 fs [sector (g) in Fig. 3]. States not shown have a peak population of less than 0.08. (a) Populated states in the double excitation subspace. (b) Populated states in the single excitation subspace. (c) Populated states in the ground-state subspace. (d) Solid line shows the total population, and the dashed line shows the renormalized purity.

FIG. 6.

Time evolution data for Ec=3.0 GV/m and τ = 0.58 fs [sector (g) in Fig. 3]. States not shown have a peak population of less than 0.08. (a) Populated states in the double excitation subspace. (b) Populated states in the single excitation subspace. (c) Populated states in the ground-state subspace. (d) Solid line shows the total population, and the dashed line shows the renormalized purity.

Close modal

In sectors (g) and (h), where the lifetime of field excitations is on the order of tens or hundreds of femtoseconds (τ on the order of 10−1 fs or 10−2 fs), a new phenomenon has to be introduced to understand why the molecular stability is declining for the fastest decay. Basis states with a non-zero photon number inherit the relevant lifetime, and here, it is short enough for the energy of these states to become a non-negligible superposition of energies according to a Lorentz distribution with the full width at half maximum Γ = ℏκ. Or put differently, states with short lifetimes experience energy broadening. Which states are affected and the extent of this broadening are shown in Fig. 7. The superposition of energies means that states that were otherwise resonantly dipole-coupled acquire an increasing average detuning, and population transfer is suppressed. The outcome is a smaller spectral overlap and a diminished population transfer between the (sharp) initial state 0,C and the (broadened) decaying states. This keeps most of the population on the potential energy surface of 0,C, where it is quickly absorbed at the end of the grid, realizing the declining stability of the MgH+ molecule that is most pronounced in sector (h) of Fig. 3.

FIG. 7.

The left-hand side shows the basis states for the system (the same as in Fig. 2). States in blue have some of its energy in the field mode. With the finite photon lifetime, their energy will spread over a superposition of energies in a Lorentzian broadening. The degree of broadening is shown on the right-hand side, plotted on the same energy scale as the left-hand side (and, for clarity, renormalized to equal amplitudes). Graph (f) corresponds to the lifetime of τ = 3.2 fs, which is the logarithmic center of sector (f) in Fig. 3, thus estimating the typical energy broadening that sector. The pattern continues, with graph (g) estimating the broadening in sector (g) and graph (h) estimating the broadening in sector (h). The dashed state has two photons and thus twice the broadening of solid blue states.

FIG. 7.

The left-hand side shows the basis states for the system (the same as in Fig. 2). States in blue have some of its energy in the field mode. With the finite photon lifetime, their energy will spread over a superposition of energies in a Lorentzian broadening. The degree of broadening is shown on the right-hand side, plotted on the same energy scale as the left-hand side (and, for clarity, renormalized to equal amplitudes). Graph (f) corresponds to the lifetime of τ = 3.2 fs, which is the logarithmic center of sector (f) in Fig. 3, thus estimating the typical energy broadening that sector. The pattern continues, with graph (g) estimating the broadening in sector (g) and graph (h) estimating the broadening in sector (h). The dashed state has two photons and thus twice the broadening of solid blue states.

Close modal

The increase in energy broadening from Fig. 7 implies that additional states will overlap in energy, and thus, their interactions should be included in the model, which occurs about half-way through sector (g) and certainly in sector (h). This means that the initial assumption about independent subspaces (see Fig. 2) will start to break down. Fortunately, the decline in population seen in sector (h) from Fig. 3 is explained by the weaker coupling between the non-broadening initial state 0,C and the decaying states. This suppression will be unaffected by any missing couplings, and therefore, the qualitative behavior with a decline in stability for very short lifetimes will still hold.

Changes in chemical properties due to the interaction with lossless optical cavities have been under intense theoretical study in the last decade.1,15 However, from this study, it is clear that long lifetimes are not a prerequisite for interesting modifications of chemical properties. In Fig. 3, a white dotted line is drawn where the energy width of the field excitation Γ = ℏκ is roughly equal to the Rabi splitting ΩR(q)=2Ecμ(q). Transition dipole moments are shown in Fig. 8 as functions of q. For finding ΩR(q) ≈ Γ, the transition dipole moment is thus estimated as μ = 1 × 10−29 cm, and a white dotted line for ΩR = Γ can be drawn. Some orders of magnitude to the left of this line lie the strong coupling regime, where ΩR(q) ≫ Γ, and to its right, energy broadening is the dominant effect. Still, an equally significant modification of MgH+ stability is observed in this regime.

FIG. 8.

Four transition dipole moments are used in the calculations. XA, XB, and AC are resonant with the cavity mode frequency and thus included in the interaction Hamiltonian [Eq. (A3)]. The transition XC is used to create the initial excited state.

FIG. 8.

Four transition dipole moments are used in the calculations. XA, XB, and AC are resonant with the cavity mode frequency and thus included in the interaction Hamiltonian [Eq. (A3)]. The transition XC is used to create the initial excited state.

Close modal

Using the same estimate for the transition dipole moment (μ = 1 × 10−29 cm), the entire parameter regime shown in Fig. 3 falls just below the ultrastrong coupling regime (where Ecμ/ωc>0.1), thus motivating the approximations made in the Hamiltonian (see  Appendix A).

We have studied the photostability of MgH+ in a lossy cavity. The dynamics have been simulated by performing nuclear wave packet dynamics via the Lindblad equation. This approach includes the vibronic decoherence caused indirectly by the loss of photons from the cavity. The studied parameter range includes decay rates inside and outside the strong coupling regime.

Deep in the strong coupling regime, for cavity lifetimes longer than the nuclear dynamics (τ ≫ 100 fs), stabilization of MgH+ is achieved through the formation of well-separated polariton states. Interference effects at the curve crossing can be observed due to the nearly fully coherent time evolution. For cavity decay rates on the order of tens of femtoseconds, the stabilization effect decreases and the aforementioned interference effects disappear. Even though this regime can still be regarded as the strong coupling regime, the coherent wave packet time evolution now competes with the effects from dissipation. The separation of polaritonic states, which is responsible for the stabilization, is now affected by the photon decay.

Shortening the photon lifetime even further, down below 1 fs, increases the stabilization again, as the molecule is rapidly funneled back into its ground state of the system. Our results suggest that, for an optimal photon lifetime in this region, low Q-factor cavities may facilitate the control of photochemical reactions,21 where several competing mechanisms are responsible for the observed phenomena: The population is transferred efficiently between the polariton states, which may be interpreted as an optimized spectral overlap of the cavity mode42 with the energy width of the nuclear wave packet. The short photon lifetime contributes to the cooling of the system. With fast enough decay, the molecule can dissipate the energy stored in the electronic excitation before dissociation takes place. The combination of both effects results in an optimal stabilization of MgH+ in this particular configuration.

It is worth noting that the stability optimum is at the border of the strong coupling regime, which has also been found in other studies.21 This suggests that it may be an interplay of strong coupling and cavity cooling, which is required to explain polaritonic chemistry experiments.

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).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The molecular Hamiltonian is here expressed in terms of the four lowest electronic states {X, A, B, C},

(A1)

The reduced mass, M, of the MgH+ molecule is 1763.2 me. 1D potential energy surfaces {VX(q), …, VC(q)} are calculated under the Born–Oppenheimer approximation with the program package Molpro33 at the CASSCF(12/10)/MRCI/ROOS level of theory34,35 and interpolated using splines to a grid with 96 grid points. Wave packets approaching the right-hand side of the grid are removed by Gaussian-shaped absorbing potentials. How much norm that has been removed is stored as a function of time, separately for each product state (see Fig. 1).

The frequency of the fundamental cavity mode ωc is chosen to match a bright transition of the Mg atom (at 285 nm).9 The cavity mode is modeled in the Fock basis {0,1,2,}. With the photonic creation and annihilation operators, â and â, respectively, the cavity Hamiltonian is formed,

(A2)

Light–matter interactions are modeled under the dipole and rotating wave approximations. Molecular transitions resonant with the cavity mode frequency are XA, XB, and AC (see Fig. 1 for their potential energy surfaces),

(A3)

The transitions AB, BC, and XC are not resonant with the cavity mode frequency, but the transition XC 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,

(A4)

To do time evolution with the Lindblad equation (1), all {L^n} 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 time-independent. □

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,

(B1)

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,

(B2)

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,

(B3)

These two terms in Eq. (B3) are performing the functions of re-populating states, and depopulating states, in that order of appearance.

The first term, responsible for re-populating states, can unfortunately not be factorized,

(B4)

However, we treat it as if it could be anyway and define the left parenthesis from Eq. (B4) as Ŝ1,

(B5)

The multiplication from Eq. (B4) is then expanded,

(B6)

Terms with a single κn are identical to the desired terms from the Lindblad equation (B3), but there are also several unwanted cross-terms. Expressing ρ^ in the employed basis {ei} and using the result that all operators {L^n}{L^ij} are of the form of Eq. (B2), we simplify the operator for each type of term, beginning with the single κn term,

(B7)

These desired terms thus correspond to increasing population in diagonal matrix elements of ρ^ (due to eiei), in proportion to other diagonal matrix elements of ρ^ (due to ρjj). Then, the cross-terms are simplified as

(B8)

The {L^ij} operators are required to be orthogonal, which means that either im or jn. That is, either eiem 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 cross-terms from Eq. (B6). Let D[] be this operation, which sets off-diagonal matrix elements to zero in the employed basis. The initial sum over multiple matrix multiplications can now be rewritten using Ŝ1,

(B9)

In contrast to the re-populating term from the Lindblad equation (B3), the other term, responsible for depopulation, can be factorized,

(B10)

Defining the left parenthesis above as the operator Ŝ2, we can reduce this to two matrix multiplications,

(B11)

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 {L^n} 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,

(B12)

The Lindblad equation is solved numerically, using the ode4536 (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.

FIG. 9.

Data for method accuracy benchmark. The remaining population in the system after 500 fs of time evolution. The black curve is obtained from 1000 QDng, Chebyshev propagator, calculations (QDng method) and considered the baseline for accuracy comparisons. Pink crosses are the results from 100 Octave, ode45, calculations (Octave method).

FIG. 9.

Data for method accuracy benchmark. The remaining population in the system after 500 fs of time evolution. The black curve is obtained from 1000 QDng, Chebyshev propagator, calculations (QDng method) and considered the baseline for accuracy comparisons. Pink crosses are the results from 100 Octave, ode45, calculations (Octave method).

Close modal

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 × 104 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.

1.
R. F.
Ribeiro
,
L. A.
Martínez-Martínez
,
M.
Du
,
J.
Campos-Gonzalez-Angulo
, and
J.
Yuen-Zhou
, “
Polariton chemistry: Controlling molecular dynamics with optical cavities
,”
Chem. Sci.
9
,
6325
6339
(
2018
).
2.
T.
Schwartz
,
J. A.
Hutchison
,
C.
Genet
, and
T. W.
Ebbesen
, “
Reversible switching of ultrastrong light-molecule coupling
,”
Phys. Rev. Lett.
106
,
196405
(
2011
).
3.
A.
Thomas
,
L.
Lethuillier-Karl
,
K.
Nagarajan
,
R. M. A.
Vergauwe
,
J.
George
,
T.
Chervy
,
A.
Shalabney
,
E.
Devaux
,
C.
Genet
,
J.
Moran
, and
T. W.
Ebbesen
, “
Tilting a ground-state reactivity landscape by vibrational strong coupling
,”
Science
363
,
615
619
(
2019
).
4.
E. T.
Jaynes
and
F. W.
Cummings
, “
Comparison of quantum and semiclassical radiation theories with application to the beam maser
,”
Proc. IEEE
51
,
89
109
(
1963
).
5.
N.
Jia
,
N.
Schine
,
A.
Georgakopoulos
,
A.
Ryou
,
L. W.
Clark
,
A.
Sommer
, and
J.
Simon
, “
A strongly interacting polaritonic quantum dot
,”
Nat. Phys.
14
,
550
554
(
2018
).
6.
T.
Byrnes
,
N. Y.
Kim
, and
Y.
Yamamoto
, “
Exciton–polariton condensates
,”
Nat. Phys.
10
,
803
813
(
2014
).
7.
H. L.
Luk
,
J.
Feist
,
J. J.
Toppari
, and
G.
Groenhof
, “
Multiscale molecular dynamics simulations of polaritonic chemistry
,”
J. Chem. Theory Comput.
13
,
4324
4335
(
2017
).
8.
O.
Vendrell
, “
Coherent dynamics in cavity femtochemistry: Application of the multi-configuration time-dependent Hartree method
,”
Chem. Phys.
509
,
55
65
(
2018
), part of the Special Issue: High-dimensional quantum dynamics (on the occasion of the 70th birthday of Hans-Dieter Meyer).
9.
E.
Davidsson
and
M.
Kowalewski
, “
Atom assisted photochemistry in optical cavities
,”
J. Phys. Chem. A
124
,
4672
4677
(
2020
).
10.
F.
Brennecke
,
T.
Donner
,
S.
Ritter
,
T.
Bourdel
,
M.
Köhl
, and
T.
Esslinger
, “
Cavity QED with a Bose–Einstein condensate
,”
Nature
450
,
268
271
(
2007
).
11.
D. G.
Baranov
,
M.
Wersäll
,
J.
Cuadra
,
T. J.
Antosiewicz
, and
T.
Shegai
, “
Novel nanostructures and materials for strong light–matter interactions
,”
ACS Photonics
5
,
24
42
(
2018
).
12.
J.
Flick
and
P.
Narang
, “
Cavity-correlated electron-nuclear dynamics from first principles
,”
Phys. Rev. Lett.
121
,
113002
(
2018
).
13.
F.
Herrera
and
F. C.
Spano
, “
Cavity-controlled chemistry in molecular ensembles
,”
Phys. Rev. Lett.
116
,
238301
(
2016
).
14.
E.
Cortés
,
W.
Xie
,
J.
Cambiasso
,
A. S.
Jermyn
,
R.
Sundararaman
,
P.
Narang
,
S.
Schlücker
, and
S. A.
Maier
, “
Plasmonic hot electron transport drives nano-localized chemistry
,”
Nat. Commun.
8
,
14880
(
2017
).
15.
A.
Mandal
and
P.
Huo
, “
Investigating new reactivities enabled by polariton photochemistry
,”
J. Phys. Chem. Lett.
10
,
5519
5529
(
2019
).
16.
S. M.
Ashrafi
,
R.
Malekfar
,
A. R.
Bahrampour
, and
J.
Feist
, “
Long-distance heat transfer between molecular systems through a hybrid plasmonic-photonic nanoresonator
,” (unpublished 2020).
17.
J.
Flick
,
M.
Ruggenthaler
,
H.
Appel
, and
A.
Rubio
, “
Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (QED) chemistry
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
3026
3034
(
2017
).
18.
M.
Kowalewski
and
S.
Mukamel
, “
Manipulating molecules with quantum light
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
3278
3280
(
2017
).
19.
I. S.
Ulusoy
and
O.
Vendrell
, “
Dynamics and spectroscopy of molecular ensembles in a lossy microcavity
,”
J. Chem. Phys.
153
,
044108
(
2020
).
20.
P.
Antoniou
,
F.
Suchanek
,
J. F.
Varner
, and
J. J.
Foley
, “
Role of cavity losses on nonadiabatic couplings and dynamics in polaritonic chemistry
,”
J. Phys. Chem. Lett.
11
,
9063
9069
(
2020
).
21.
S.
Felicetti
,
J.
Fregoni
,
T.
Schnappinger
,
S.
Reiter
,
R.
de Vivie-Riedle
, and
J.
Feist
, “
Photoprotecting uracil by coupling with lossy nanocavities
,”
J. Phys. Chem. Lett.
11
,
8810
8818
(
2020
).
22.
J. T.
Hugall
,
A.
Singh
, and
N. F.
van Hulst
, “
Plasmonic cavity coupling
,”
ACS Photonics
5
,
43
53
(
2018
).
23.
J. J.
Baumberg
,
J.
Aizpurua
,
M. H.
Mikkelsen
, and
D. R.
Smith
, “
Extreme nanophotonics from ultrathin metallic gaps
,”
Nat. Mater.
18
,
668
678
(
2019
).
24.
R.-Q.
Li
,
D.
Hernángomez-Pérez
,
F. J.
García-Vidal
, and
A. I.
Fernández-Domínguez
, “
Transformation optics approach to plasmon-exciton strong coupling in nanocavities
,”
Phys. Rev. Lett.
117
,
107401
(
2016
).
25.
C.
Climent
,
J.
Galego
,
F. J.
Garcia-Vidal
, and
J.
Feist
, “
Plasmonic nanocavities enable self-induced electrostatic catalysis
,”
Angew. Chem., Int. Ed.
58
,
8698
8702
(
2019
).
26.
H.-P.
Breuer
and
F.
Petruccione
,
The Theory of Open Quantum Systems
(
Oxford University Press
,
2010
).
27.
D.
Manzano
, “
A short introduction to the Lindblad master equation
,”
AIP Adv.
10
,
025106
(
2020
).
28.
S.
Haroche
and
J.-M.
Raimond
,
Exploring the Quantum: Atoms, Cavities, and Photons (Oxford Graduate Texts)
, 1st ed. (
Oxford University Press
,
2006
).
29.
J.
Fregoni
,
G.
Granucci
,
E.
Coccia
,
M.
Persico
, and
S.
Corni
, “
Manipulating azobenzene photoisomerization through strong light–molecule coupling
,”
Nat. Commun.
9
,
4688
(
2018
).
30.
J.
Fregoni
,
G.
Granucci
,
M.
Persico
, and
S.
Corni
, “
Strong coupling with light enhances the photoisomerization quantum yield of azobenzene
,”
Chem
6
,
250
265
(
2020
).
31.
J.
Flick
,
H.
Appel
,
M.
Ruggenthaler
, and
A.
Rubio
, “
Cavity Born–Oppenheimer approximation for correlated electron–nuclear-photon systems
,”
J. Chem. Theory Comput.
13
,
1616
1625
(
2017
).
32.
M.
Kowalewski
,
K.
Bennett
, and
S.
Mukamel
, “
Cavity femtochemistry: Manipulating nonadiabatic dynamics at avoided crossings
,”
J. Phys. Chem. Lett.
7
,
2050
2054
(
2016
).
33.
H. J.
Werner
,
P. J.
Knowles
,
R.
Lindh
,
F. R.
Manby
,
M.
Schütz
 et al., MOLPRO, version 2006.1, a package of ab initio programs, Cardiff, UK,
2006
.
34.
P.-O.
Widmark
,
P.-Å.
Malmqvist
, and
B. O.
Roos
, “
Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions
,”
Theor. Chim. Acta
77
,
291
306
(
1990
).
35.
P.-O.
Widmark
,
B. J.
Persson
, and
B. O.
Roos
, “
Density matrix averaged atomic natural orbital (ANO) basis sets for correlated molecular wave functions
,”
Theor. Chim. Acta
79
,
419
432
(
1991
).
36.
See https://octave.sourceforge.io/octave/function/ode45.html for Function refecence: ode45; accessed 16 October 2020.
37.
See https://www.gnu.org/software/octave/index for Gnu octave; accessed 16 October 2020.
38.
H.
Tal-Ezer
and
R.
Kosloff
, “
An accurate and efficient scheme for propagating the time dependent Schrödinger equation
,”
J. Chem. Phys.
81
,
3967
3971
(
1984
).
39.
J.
Galego
,
F. J.
Garcia-Vidal
, and
J.
Feist
, “
Cavity-induced modifications of molecular structure in the strong-coupling regime
,”
Phys. Rev. X
5
,
041022
(
2015
).
40.
M.
Kowalewski
,
K.
Bennett
, and
S.
Mukamel
, “
Non-adiabatic dynamics of molecules in optical cavities
,”
J. Chem. Phys.
144
,
054309
(
2016
).
41.
D.
Wang
,
Å.
Larson
,
H. O.
Karlsson
, and
T.
Hansson
, “
Molecular quantum wavepacket revivals in coupled electronic states
,”
Chem. Phys. Lett.
449
,
266
271
(
2007
).
42.
M.
Motsch
,
M.
Zeppenfeld
,
P. W. H.
Pinkse
, and
G.
Rempe
, “
Cavity-enhanced Rayleigh scattering
,”
New J. Phys.
12
,
063022
(
2010
).