Ab Initio Linear-Response Approach to Vibro-Polaritons in the Cavity Born-Oppenheimer Approximation.

Recent years have seen significant developments in the study of strong light-matter coupling including the control of chemical reactions by altering the vibrational normal modes of molecules. In the vibrational strong coupling regime, the normal modes of the system become hybrid modes which mix nuclear, electronic, and photonic degrees of freedom. First-principles methods capable of treating light and matter degrees of freedom on the same level of theory are an important tool in understanding such systems. In this work, we develop and apply a generalized force constant matrix approach to the study of mixed vibration-photon (vibro-polariton) states of molecules based on the cavity Born-Oppenheimer approximation and quantum-electrodynamical density-functional theory. With this method, vibro-polariton modes and infrared spectra can be computed via linear-response techniques analogous to those widely used for conventional vibrations and phonons. We also develop an accurate model that highlights the consistent treatment of cavity-coupled electrons in the vibrational strong coupling regime. These electronic effects appear as new terms previously disregarded by simpler models. This effective model also allows for an accurate extrapolation of single and two molecule calculations to the collective strong coupling limit of hundreds of molecules. We benchmark these approaches for single and many CO2 molecules coupled to a single photon mode and the iron pentacarbonyl Fe(CO)5 molecule coupled to a few photon modes. Our results are the first ab initio results for collective vibrational strong coupling effects. This framework for efficient computations of vibro-polaritons paves the way to a systematic description and improved understanding of the behavior of chemical systems in vibrational strong coupling.

Introduction

Recent experimental progress in the field of polaritonic chemistry has demonstrated the possibilities of altering chemical and material properties with the strong coupling of electromagnetic fields and vibrational degrees of freedom. In this vibrational strong coupling regime, light and matter degrees of freedom hybridize forming vibro-polaritons.[1] It has been demonstrated that in this regime coupled cavity photons can be tuned to influence chemical reactivity,[2] vibrational energy redistribution,[3] optical spectra,[4,5] Raman spectra,[6] two-dimensional spectroscopy,[7] relaxation dynamics,[8] ultrafast thermal modification,[9] and even superconductivity,[10] among others. These experimental works have been complemented by various theoretical efforts,[11−19] one development, in particular, to describe these experiments is the introduction of effective vibro-polariton Hamiltonians[4,5,20,21] that include the vibrational degree of freedom via normal modes (vibrations in molecular or phonons in solid-state systems). These normal modes can be obtained, for example, experimentally from infrared spectroscopy[4,5] or numerically from first principles using electronic structure theory methods.[22,23] Although conventional electronic structure methods are not directly applicable to the light–matter strong coupling regime due to their negligence of the quantum electromagnetic field, here they can be used to calculate the vibrational normal modes of the matter system based on the force constant matrix. These vibrational normal modes are then coupled to the photon modes of the electromagnetic field. Such Hamiltonians have been applied successfully to describe various experimental findings.[4,5,21] One limitation of these vibro-polariton Hamiltonians that only include vibrational modes and photon modes explicitly is that self-consistent effects of the electron–photon interaction are neglected. In addition, this description usually aims at including only the relevant degrees of freedom of the system explicitly. While for simpler systems the relevant degrees of freedom can be known beforehand, in general, and for more complex situations, these variables are not always known.

An alternative route to simulate vibrational strong coupling is offered by first-principles methods that treat the full matter–photon Hamiltonian explicitly. Examples include the generalization of Hartree–Fock,[24,25] QED coupled-cluster (QED-CC) theory,[25−27] and quantum-electrodynamical density-functional theory (QEDFT).[28,29] In the QEDFT framework, vibrational strong coupling has been simulated in the time domain capturing the dynamics of the system to analyze optical spectra[30] or chemical reactivity,[31] but the full framework to describe vibrational strong light–matter coupling within linear-response theory has not yet been developed. While explicit calculations in the time domain have their advantages for simulating complex and anharmonic dynamics, information about vibro-polaritonic modes can be obtained from linear-response calculations more efficiently. One limitation of these first-principles methods is their relatively high computational cost, which effectively limits calculations to the single or few molecule limit, which is the opposite limit of experiments in the collective strong coupling regime.

In this work, we introduce an efficient framework to calculate properties of systems under vibrational strong coupling from first principles. We introduce the generalized force constant matrix, where eigenvectors and eigenvalues give rise to vibro-polaritonic normal modes of the correlated matter–photon system and the frequencies of the vibro-polaritons. In addition, we develop an accurate effective model that includes light–matter feedback terms that have been previously disregarded. We show that this effective model allows for extrapolation of first-principles calculations to the collective strong coupling regime. We exemplify these methods by calculating optical spectra for single and many CO2 molecules in optical cavities, as well as for the iron pentacarbonyl Fe(CO)5 coupled to a multiphoton mode setup.

Theory of Vibro-Polaritons

In the following section, we develop the framework to describe vibro-polaritons in the linear-response regime from first principles. We start by discussing the Hamiltonian for light–matter coupled systems in the length gauge and in the dipole approximation.[30,32,33] For the vibrational strong coupling regime, it has been shown that the cavity Born–Oppenheimer approximation (CBOA) can yield an accurate description of the system.[11,13,14,16,33] This method is based on the adiabatic approximation that allows separation of the electronic degrees of freedom from the nuclear and photonic degrees of freedom. As a consequence, the photonic degrees of freedom are described as analogous to the nuclear degrees of freedom in the conventional Born–Oppenheimer approximation.[34] With this framework, the nuclear–photon dynamics of a set of Nnuc nuclei with coordinates and photon modes with photon displacement coordinates = (q1, q2, ... qα ...) are obtained by solvingwith nuclear and photonic kinetic energies T̂nuc and T̂pt, respectively, and E denotes the cavity Born–Oppenheimer (CBO) potential-energy surface for the lowest energy electronic state of the system. We emphasize here that the nuclear and photonic degrees of freedom in eq are still treated fully quantum mechanically, thus capable of describing nonclassical states of light, for example, Fock or squeezed states. In practice, we can obtain the CBO potential-energy surfaces by diagonalizing the electronic Hamiltonian of N electrons that now parametrically depends on the nuclear and the photonic coordinates withwhere

Here, T̂ describes the electronic kinetic energy, V̂ the electron–electron interaction, V̂ the electron–nuclear interaction, and V̂nuc–nuc the nuclear–nuclear interactions, respectively. In eq , we include the matter–photon coupling aswhere the α runs over photon modes. The photon displacement coordinate qα couples to the electronic and nuclear dipole moment operator, which is given by , where is the electronic position operator, and e describes the elementary charge and Z the charge of the Ith nucleus. The frequency ωα and the coupling strength λα define the parameters of the individual photon modes. Here, the coupling strength λα is related to the amplitude of the electric field at the center of charge of the molecule for the cavity photon mode α via . In general, strong coupling is defined to occur when the coupling strength is larger than decoherence rates in the system.[35] Since our simulations are of isolated systems, then as long as λ is nonzero, the system is technically in the strong coupling regime, and values where the coupling becomes “strong” depend on the decoherence rates in a given experiment.

In this work we treat the μ̂2 term in the electronic potential using a mean field approximation as described in the mean field electronic potential section of the Supporting Information.

Having setup the Hamiltonian of the matter–photon system, we can proceed to determine the vibro-polaritonic normal modes. In the first step, we define the effective nuclear and photonic forces and calculate the equilibrium configuration of the system. To derive the forces acting on nuclear and photonic degrees of freedom in the presence of matter–photon coupling, we apply the Hellman–Feynman theorem. The forces on nucleus I along the κ direction are given bywhere E is the ground-state CBO energy of the system governed by the Hamiltonian in eq , and R indicates the κ direction component of the position of nucleus I and ⟨..⟩ an expectation value evaluated using the electronic states at particular values of and . There is also an effective force on the photon displacement coordinate, which is given byThe equilibrium position with the ground-state energy E0 with respect to and is now defined by minimization of energy, as defined by eq , and thus vanishing forces, i.e., with the electronic Hamiltonian in eq .

Derivatives of the forces in eqs and 6 provide sufficient information on the electronic energy surface to construct the effective Hamiltonian for the nuclei and photons of eq . The CBO energy of the coupled light–matter system with small perturbations around the equilibrium configuration can be expressed aswhere ΔR are displacements of atom I along direction κ and Δqα perturbations of photon displacement qα; E0 is the energy of the equilibrium configuration. We have defined the matricesThe vibro-polariton eigendisplacements η of the light–matter coupled system and the vibro-polariton eigenfrequencies ω can be obtained by solving the generalized eigenvalue problemwhere M = Mδδ, and M is the mass of nucleus I and a identity matrix.a The matrices C and M̃ as well as generalized eigendisplacements η can be used to rewrite eq in a more compact formwhere now C acts as a generalized force constant matrix which includes both nuclear and photon degrees of freedom. The analogous generalized dynamical matrix can then be defined aswith eigenvalues ω2 and vibro-polariton eigenvectors U. For a normalized set of U, the eigendisplacements are related by η = U/(M̃)1/2, where the eigendisplacements are normalized to obey ηM̃η = 1. The eigenvalues of this generalized dynamical matrix give the frequencies ω of eq which in the harmonic limit are equivalent to the energies ϵ of eq .

Analyzing the structure of the force constant matrix, we find a 2 × 2 block structure of C (left side of eq ) reminiscent of the electron–photon linear-response polaritonic Casida equation.[36] We find the matter block C( and photon block C( on the diagonal coupled by an off-diagonal block C(, which introduces the matter–photon coupling. For the case of λ = 0, the off-diagonal blocks vanish, and the matter block reduces the standard force constant matrix.[23] We further note that while in the polaritonic Casida equation the photon block is strictly diagonal, since there is no explicit photon–photon interaction present, the same is not true for the generalized force constant matrix here. In this case, the photon block C( is not diagonal due to an effective photon–photon interaction between individual photon modes that originates from the electron–photon interaction. The manifestations of this effective photon–photon interaction is discussed in Section . A schematic representation for a single vibrational mode coupled to a single photon mode is given in Figure .

Figure 1

Schematic representation of vibro-polaritonic excitations in an optical cavity. On the left, an infrared-active vibrational excitation of CO2 is depicted at a particular energy level. On the right, a particular photon mode of the cavity is depicted as an excitation of photon displacement coordinate qα within a harmonic potential. Under strong coupling, these vibration and photon modes hybridize leading to upper and lower polaritons as depicted by the two states in the center. Note that the eigenvectors of these hybrid states have opposite signed qα components as depicted with the blue arrows.

We can now obtain the infrared spectrum from the eigenvectors and eigenfrequencies of the generalized dynamical matrix and the vibro-polariton mode effective charges. Both quantities are defined analogously to the case of conventional linear-response theory of vibrations/phonons.[22] The mode effective charge of vibro-polariton normal mode m along direction κ is given byUsing these effective charges, the corresponding infrared spectrum I can be constructed aswhere the peaks at the frequencies ω with amplitudes Z are broadened by the Lorentzian L(Ω, ω, δ)[37]

This definition of the infrared spectra corresponds to the electric dipole autocorrelation function, or the Fourier transform of the dipole moment as a function of time, and is typically used in standard first-principles calculations.[38] In an optical cavity, this IR spectrum corresponds to a situation where the IR spectrum is measured perpendicular to the cavity axis.[39] Note that in the definition presented here we have included the electronic response of the dipole to photon coordinate qα at fixed ionic positions (the second term in eq ), a term often disregarded in other first-principles works. We note that that in the context of optical cavities other spectra are also measured, such as reflection and transmission spectra through the mirror, which can be defined in terms of a photonic autocorrelation function.[40]

Model for Vibro-Polaritons

To elucidate the various microscopic contributions to the results of the full first-principles theory, we now develop an equivalent vibro-polaritonic model. We first rewrite eq with the matter degrees of freedom rotated into a basis of uncoupled vibrational normal modes. The nuclear displacements from the equilibrium configuration can be expressed in terms of vibrational mode amplitudes N which specify the change of ionic positions. For a general ionic displacement given by a set of ΔR, the corresponding set of N are given by where are the normal vibrational mode eigendisplacements of the uncoupled problem at λ = 0. Force constant matrix elements can be written in terms of N and qα by expanding the expectation values present in eqs and 6 to linear order. Then, the dipole expectation value readsWe express the first force contribution on the left side of eq in terms mixed second derivatives given by matrices Θ and Ξ, defined explicitly in the model derivation section of the Supporting Information. With the above expansions and change of basis, we can define the following harmonic modelwhere T includes the kinetic energy of the vibrational modes (N) and the photon modes (qα). Z is the ionic contribution to the uncoupled vibrational mode effective charge of eq and given explicitly by . Additional details on the derivation of the model can be found in the Supporting Information.

We find that the matter–photon coupling strength, i.e., the term proportional to Nqα, depends on the quantity . As discussed before, the dipole moment of the system consists of two contributions, a nuclear one and an electronic one. As a consequence, the term also includes two contributions: the nuclear dipole moment, as well as the change of the electric dipole moment due to a change in nuclear configuration. Analogously, the term describes the change of the electric dipole moment due to a change in photon coordinate qα. We note that while photon modes are not explicitly coupled in eq , i.e., there is no photon–photon coupling term, the term introduces effective photon–photon coupling in the vibro-polariton model. Since the model describes a set of interacting quantum harmonic oscillators, it can also be solved analytically.[41]

For a detailed illustration, we now consider the model of eq for a single photon mode coupled to a single vibration mode with the relevant vibration only influencing the dipole moment along the direction of photon polarization. With these simplifications, we can drop the mode indices, label the vibration mode frequency with subscript N and the photon mode with subscript q, and treat the dipole moment μ and coupling strength vector λ as scalars. Then, eq reduces toHere, we find two effective frequencies: (i) the effective frequency of the vibrational normal mode that is given byand (ii) the effective frequency of the photon mode that is given by

In addition, we have the effective interaction strength that is given by

The resulting eigenvalues are then the upper and lower polaritons with frequenciesWe find that the photon frequency ω at which resonance occurs is then not that of the bare phonon mode (ω), instead resonance occurs when ω̃ = ω̃. The model of vibro-polaritons in eq contains three parameters which have a dependence on the coupling strength λ. These are the derivatives of the dipole with respect to photon displacement d⟨μ̂⟩/dq and nuclear positions d⟨μ̂⟩/dN, as well as derivatives of the Coulomb forces on nuclei expressed as Ξ, where derivatives with respect to nuclei positions are in a basis of uncoupled vibrational normal modes. For an uncoupled system (λ = 0), both Ξ and d⟨μ̂⟩/dq are zero. The λ dependence of all three of these parameters is a result of coupling strength and q dependence of the electronic state. Changes in the electronic state with λ and q change the force terms written as expectation values (i.e., within ⟨..⟩) in eqs and 6. This effect is then captured by these λ-dependent parameters of the model, written in the basis of vibrational normal modes of the uncoupled system.

An alternative approach for treating vibro-polaritons from first principles is to use a model which couples the cavity photon mode to matter vibrations. The parameters, which characterize the matter vibrations, are then obtained from standard first-principles methods.[4,5,20,21] In such an approach, the modification of the electronic potential energy due to the cavity is not taken to account consistently. Such models correspond to neglecting the coupling strength-dependent terms in eq . Setting Ξ and equal to zero and setting equal to its λ = 0 value recovers such a model which can be constructed without cavity modification of the electronic potential energy. We refer to this approximation as the “μ2 model” as it still contains quadratic dipole terms from eq . If one further neglects this term in the NN coupling that is of order λ2, one arrives at a system of bilinearly coupled vibrational and photon oscillators similar to the Hopfield model.[42] In this simplified model, the single vibration–single photon effective frequencies are simply the bare vibrational and cavity normal modes, and any λ dependence of is neglected in the effective coupling strength term.

Results and Discussion

In this section, we illustrate the developed approach on single and many CO2 molecules, as well as the iron pentacarbonyl Fe(CO5). We list the numerical details for these calculations in the numerical details section of the Supporting Information. We start by discussing the case of CO2 molecule(s).

Single CO2 in an Optical Cavity

Figure shows the computed vibro-polariton normal mode frequencies (vertical lines) and Lorentzian broadened infrared spectra (black curves) at various values of the coupling strength λ. The color of the vertical lines corresponds to the absolute value of the photon component of the corresponding vibro-polariton normal mode eigenvector. In this calculation, one photon mode is included with frequency ωα = 2430 cm–1 chosen to be near resonance with the 2436 cm–1 asymmetric stretching vibration mode of the uncoupled system. We choose this slight detuning to be consistent with the calculation in ref (30). The direction of the λα vector, which sets the photon mode polarization direction, was chosen to be aligned with the oscillating dipole moment along the C–O bonds as indicated by the blue arrow in the inset of the bottom plot in Figure . A Lorentzian broadening of δ = 25 cm–1 was used. By increasing the coupling strength from λ = 0 to λ = 0.1, we observe the Rabi splitting of the vibrational mode at 2436 cm–1 between the upper and lower vibro-polariton branches. We note that the observed values are in quantitative agreement with the fully time-dependent results of ref (30). As expected, neither the noninfrared (IR)-active symmetric stretching mode at 1363 cm–1 or the degenerate bending modes at 607 cm–1 couple to the cavity. The latter of which is only IR active along directions orthogonal to the cavity polarization. The eigenvectors of the two polariton modes are linear combinations of the asymmetric stretching mode and the photon displacement. The lower polariton has a photon displacement aligned with the vibration mode dipole and a larger photon component. While for the upper polariton eigenvector, the photon displacement is antialigned with the vibration mode dipole, and the photon component is smaller. Definitions of ultrastrong and deep strong coupling are given in terms of Rabi model parameter g in eq 1of ref (35). The closest analog of g in this work would be . If is taken to be the value from the zero coupling case, then the onset of ultrastrong coupling would be λ > 0.125 au, and deep strong coupling would be λ > 1.25 au for the case of the CO2 molecule. However, our calculations include both changes in with λ as well as . Both of these quantities influence the level of splitting and make direct comparison with definitions of ultrastrong and deep strong coupling in terms of the Rabi model parameters not straightforward. An alternative quantity to measure the strength of the light–matter coupling can be obtained in terms of the ratio of the Rabi splitting and the cavity frequency, where we find 8.5% in the case of λ = 0.05 au and over 18% for λ = 0.1 au.

Figure 2

CO2 IR spectra (black curve) for different λ (in atomic units) values for cavity frequency ωα = 2430 cm–1 and eigenvalues (vertical lines) colored by photonic character. The inset in the λ = 0 plot shows the CO2 molecule with the blue arrow indicating the polarization of the photon mode.

Notable in the results is the asymmetry in the Rabi splitting, especially in the strong coupling regime. The lower polariton is seen to have a more intense IR peak and a larger frequency shift with respect to the frequency of the bare photon mode than the upper polariton. This behavior is despite the finding that the lower polariton has a smaller matter and larger photon contribution than the upper polariton as indicated by the peak color. We find that the IR amplitudes here are dominated by the change in the electronic contribution to the dipole moment due to the change in photon displacement q, i.e., the term in eq . While the derivative of the dipole moment with respect to q is smaller in magnitude than the corresponding matter contribution (the Born effective charges), the photon component of the eigendisplacements can be much larger than the matter components as the photon components are not reduced by a factor relating to their mass (from M̃ of eqs and 12). Reference (20) found for the IR spectra of chemical systems under vibrational strong coupling in resonant setups that the lower polariton also has a smaller matter contribution than the upper polariton. This is identified as the source of the smaller IR peak of the lower polariton which is also observed in ref (43). However, in both of these works, changes in (electronic) dipole moment due to the cavity mode displacement (at fixed nuclei postions) are not accounted for. When this effect is included, the IR peak of the lower polariton increases, especially at large coupling strengths where the effect is larger and the mode has more photon character. Interestingly, in the case of strong coupling to an electronic excitation, refs (36) and (44) also show larger amplitude absorptions for peaks with more photon characters.

The asymmetry in the frequency splitting for the upper and lower polaritons can be understood by examining the two mode model presented in eqs –22. The λ-dependent parameters Ξ, , and enter in a manner which shifts the effective frequencies of both the vibrational and photon modes. Then, even when the cavity mode is tuned to the frequency of the vibration mode, these effective frequencies differ, and thus, splitting is not symmetric around the original vibration frequency. The λ dependence of each of these terms is a result of the electronic response to the cavity potential.

In the next step, we compare different effective models to the discussed first-principles results and analyze the individual terms in eq in more detail. Figure A compares the upper and lower polariton frequencies at different levels of modeling as a function of coupling strength. The results of the generalized dynamical matrix approach using the first-principles theory described in Section (shown in black) are seen to be in near perfect agreement with the full two mode model (shown in blue) of eq . In addition, we compare to two additional approximate models, which show significant differences from the full model in the ultrastrong coupling regime. The first model, shown in green dotted lines, which we term the μ2 model, corresponds to results where all of the λ dependence of the model parameters in eq have been neglected so that Ξ = d⟨μ⟩/dq = 0 and d⟨μ⟩/dN are taken as the values from the uncoupled case. The second model, shown in the orange dotted line corresponds to a Hopfield type model where in addition to the approximations made for the μ model the λ2 term from eq is also set to zero (equivalent to dropping the μ2 term in Vpt–). For a cavity mode precisely in resonance to a vibration mode, the Hopfield model maintains perfect symmetric splitting up to the extremely strong coupling regime. While the inclusion of the μ2 term does permit some asymmetry in the splitting, it is seen that when parametrized by first-principles results from the λ = 0 limit this asymmetry is relatively minimal, and results do not differ much from the Hopfield model. While some asymmetry is also present due to the small detuning of photon and vibration mode in our calculations, it is only when coupling-dependent model parameters obtained from QEDFT are included that the more dramatic asymmetric splitting is recovered. Figure B shows how various terms in the model vary with coupling strength λ. The largest λ-dependent contribution is seen to come from the d⟨μ⟩/dq term. This change in electronic dipole moment due to the photon displacement shifts the effective cavity mode frequency away from resonance with the phonon mode.

Figure 3

(A) CO2 mode splitting at various levels of modeling, see main text for definitions with cavity frequency ωα = 2430 cm–1. (B) Change of the different model parameters with coupling strength λ (in atomic units). The signs of N and q have been chosen such that all curves are positive everywhere.

Collective Strong Coupling Limit with Many CO2 Molecules

First-Principles Results

Rigorous first-principles approaches in the treatment of strong light–matter coupling have largely been applied to the problem of a single molecule strongly coupled to cavity photon modes. However, experimentally strong coupling is typically achieved via “collective coupling” where coupling strength is enhanced by increasing the number of emitters in the cavity.[45] The increased computational efficiency of the linear-response method presented in Section enables some aspects of the collective coupling regime to be accessible within QEDFT. We have simulated chains of CO2 molecules aligned along their C–O bond directions coupled to a cavity mode polarized along this same direction. Molecules are chosen to be spaced 20 Bohr apart to simulate the dilute gas limit. Figure shows comparisons between QEDFT results for a single molecule, Nmol molecules, and the results of the many molecule model presented in Section . A Lorentzian broadening of δ = 10 cm–1 was used. In each of these plots, one can see the lower and upper polaritons similar to those observed in the single molecule case but also Nmol – 1 dark modes near 2436 cm–1 with no IR amplitude. The Rabi splitting and IR spectra in the very strongly coupled single molecule case and more weakly coupled Nmol case are nearly identical with only some differences in the lower polariton frequencies at a very large number of molecules/very strong coupling.

Figure 4

Comparisons between full QEDFT results for NmolCO2 molecules at λ = 0.05 (in orange) with results for single CO2 molecules at (in blue). (λ is in atomic units) Also shown (in green) are results for a model of the form in eq for the Nmol case but constructed using parameters from QEDFT calculations with only two CO2 molecules. Vertical lines are used to indicate frequencies of the normal modes, while curves show broadened IR spectra. To facilitate comparison, the IR amplitudes have been scaled by Nmol1/2.

Similar to the case of a single coupled molecule, the lower (upper) polariton eigendisplacements consist of the original asymmetric stretching mode aligned (antialigned) with the photon displacement. However, now in the multimode case, the collective upper and lower polaritons consist of every molecule experiencing this asymmetric stretching in phase. The multimolecule setup also results in a number of dark modes which correspond to combinations of the original asymmetric stretching modes on each molecule but in such a way that the overall dipole moment when freezing has one of these collective dark modes as zero.

Modeling Larger Numbers of Molecules

The similarity between the results of a single strongly coupled molecule with multiple, more weakly coupled molecules suggests that within the level of theory applied in this work the microscopic description of one or two molecules can capture the relevant physics for many molecules coupled to the cavity in the dilute limit. To this end, we construct a model of the form presented in eq with Nmol CO2 molecules coupled to the same cavity mode as in previous sections at a coupling strength of λ(. Nearly all parameters in this model can be obtained from first-principles calculations of a single molecule with coupling strength except for certain elements of the Ξ matrix which we obtain from first-principles calculations for two molecules with coupling strength of .b The Nmol model consists of the same photon modes as single molecule case so ωα( = ωα(1) and Nmol copies of the vibration modes from a single uncoupled molecule. To simplify the notation for mapping model parameters of the Nmol system to the parameters of corresponding one or two model parameters, we have introduced the superscript indicating the number of molecules in the model a particular parameter corresponds to. Since we include copies of the original, single molecule, and vibrational modes as our starting basis, it is convenient to write our nuclear degrees of freedom with two indices: a molecular index and vibrational mode index I which corresponds to a normal mode of the uncoupled single molecule system. Together the pair of indices corresponds to an atomic displacement on molecule according to the eigendisplacement of the single molecule vibrational mode given by . So for Nions ions in each molecule in three dimensions , , and . Within the dipole approximation, a change in qα will result in a change in dipole moment for all molecules in the system, so the susceptibility must be scaled for the model as . For the choice of basis consistent with the above definitions, Ξ has a block structure where on-diagonal blocks correspond to coupling between vibration modes on the same molecule, and off-diagonal blocks correspond to coupling between vibration modes of different molecules. While the on-diagonal blocks can be obtained via ab initio calculations on a single molecule, the latter requires ab initio treatment of two molecules. The details of this construction are presented in the multimolecule delta matrix section of the Supporting Information. Since the molecules are sufficiently separated and since the long-range μ2 term is in practice handled with the mean field approximation (eq 1 in the Supporting Information), the impact of any one molecule on another is nearly independent of their distance. The impact of two molecules on a third is equivalent to a single molecule contributing the same change in dipole moment. So to harmonic order, the case of two molecules captures nearly all relevant interactions to describe Nmol molecules within the level of theory used in this work.

It is seen that within the dipole approximation there is almost no difference in the IR spectrum between a single molecule strongly coupled and a collection of molecules more weakly coupled aside from the appearance of dark modes. However, the coupling used in eq when applied to the many molecule case assumes equal coupling to all molecules in the system as there is no spatial dependence of λα. A more realistic simulation of collective coupling would facilitate better understanding of the similarities and differences between local and collective strong coupling and will be the subject of subsequent work.

Fe(CO5) in Multiple Photon Mode Setup

In the previous section, a single cavity mode was coupled to numerous degenerate vibration modes each on different molecules. In this section, we investigate a cavity coupled to multiple degenerate and nondegenerate vibration modes of a single iron pentacarbonyl molecule. Experimental data of a similar system setup have been published in ref (4).

Our system is studied with a cavity mode in resonance with several IR-active vibrations as well as with two additional photon modes to simulate additional harmonics of the cavity. As shown in the inset of the bottom panel of Figure , the coupled cavity polarization is set to be along an axis 45° from the axis of the 3-fold rotational symmetry. This setup leads to a coupling to both the vibrational mode at 2013 cm–1 which involves polar distortions along the 3-fold axis and the two degenerate vibrational modes at 1995 cm–1 which involve distortions perpendicular to the 3-fold axis. A cavity mode at 1995 cm–1 is set to couple most strongly while additional “harmonics” at frequency ratios of 3/4 (1496 cm–1) and 5/4 (2494 cm–1) are set to have a coupling strength 0.3 times that of the central mode. There are also two non-IR-active vibrational modes nearby in energy at 2016 and 2097 cm–1 which do not couple to any cavity modes. Figure depicts the normal modes of the system as vertical lines colored by photon character as well as the Lorentzian broadened IR spectra at several coupling strength magnitudes. A Lorentzian broadening of δ = 7 cm–1 was used.

Figure 5

Fe(CO)5 IR spectra for different λ values (black curve) and eigenvalues (vertical lines) colored by photon character. Here, three photon modes are present, corresponding to third, fourth, and fifth harmonics of the optical cavity. The λ values indicated in each plot are in atomic units and correspond to coupling strength of the fourth harmonic (at 1995 cm–1); the other two harmonics have been set to have one/third the coupling strength of the fourth harmonic. The inset of the λ = 0 plot shows the Fe(CO)5 molecule, with the blue vector indicating the direction of the photon mode polarization. Annotations give the values of the projections of the vibro-polariton eigenvectors to the third and fourth harmonic uncoupled photon states.

At coupling strengths with |λ| ≤ 0.02, the two outer cavity modes at 1496 and 2494 cm–1 are approximately uncoupled from the vibrational modes of the system and the central photon mode at 1995 cm–1. The IR amplitudes for the outer modes in the regime are dominated by the effect the cavity mode has on the electronic system (through the dμ/dq term). The central cavity photon mode interacts the three IR-active vibrational modes nearby in energy: the polar along the 3-fold axis (z) mode at 2013 cm–1 and the two degenerate polar modes within the plane perpendicular to the 3-fold axis (xy) modes at 1995 cm–1. The result of this cavity-induced coupling is four nondegenerate modes: a dark state which is a linear combination of the two xy modes and three polaritons which are linear combinations of the cavity photon mode, xy modes, and z mode. The dark state is still IR active but not along the cavity mode polarization direction. Similar to the case with CO2 as coupling strength is increased, the frequencies of the upper (lower) most polariton continue to grow larger (smaller), respectively, while the photon character of the polariton mode decreases (increases). The middle polariton rapidly converges to a frequency of 2006 cm–1. As coupling strength increases, the photon character of this mode decreases until there is no photon character, and the mode is made up of a linear combination of the polar z and xy vibrations. The cavity has induced a coupling between these polar vibration modes changing the eigenstate even in a regime where this eigenstate has no photon character.

At extremely strong coupling strengths with |λ| ≥ 0.05, the outer cavity mode harmonics begin to interact with other modes of the system. In the top two panels of Figure , it can be seen that even the lower frequency IR-active modes below 700 cm–1 begin to pick up some small photon character. Furthermore, as coupling strength increases to this very strong regime, the lower polariton begins to mix with this lower frequency cavity mode. The two normal modes between 1200 and 1600 cm–1 become a linear combination of vibrations, both the lowest harmonic cavity photon as well as the central cavity photon. At |λ| = 0.1, we observe that this effective photon–photon interaction has grown so strong that the photon mode components of these two modes essentially swap so that the eigenvector of the mode at 1226 cm–1 has a larger component coming from the cavity photon mode at 1995 cm–1, and the mode at 1587 cm–1 has a larger component from the cavity mode at 1496 cm–1.

Summary and Conclusion

In this work, we have introduced a first-principles framework to calculate the vibro-polaritonic normal modes of systems when light and matter are strongly coupled. Employing the cavity Born–Oppenheimer approximation to separate electronic from nuclear and photonic degrees of freedom and constructing dynamical matrices that include the photonic degree of freedom enables us to characterize these vibro-polariton states. Our approach is based on QEDFT, which makes it applicable to large system sizes while including effects of the cavity on electronic states. We demonstrate the framework on calculations for single and many CO2 molecules and iron pentacarbonyl Fe(CO5). In addition, we derive and compare to a first-principles-based model that allows for the extrapolation of first-principles calculations of few molecules to the collective strong coupling limit of molecular ensembles.

Our work opens many different avenues to explore. The techniques used here can be extended to other properties related to the system normal modes such as the low frequency Raman spectra. The vibro-polaritonic normal modes computed using the methods developed could be used as an efficient basis for exploring anharmonic couplings including interactions between polaritonic excitations.[46] The collective setup employed in this work assumes the same coupling strength for all molecules. A more realistic description where different molecular positions imply different coupling strengths due to the profile of the cavity mode could provide insight to potential differences between the collective coupling limit and small numbers of very strong coupled molecules. Such techniques can also be used to explore other related questions such as the engineering of strong coupling on single atoms[47] and local modifications of impurities due to collective coupling.[48] We have utilized the cavity Born–Oppenheimer approximation and treated the electronic portion of the two-body operator μ̂2 via a mean field potential. More sophisticated treatment of exchange-correlation effects both of electron–photon interactions and how the presence of the cavity can modify electron–electron interactions are of interest. Such more advanced treatments will be especially important when energy surfaces are sufficiently close together, and the validity of the CBOA should be carefully tested. Utilizing the methods developed in this work, potentially along with these extensions, experimentally relevant molecules can be studied to gain new insights on cavity modification of chemical reactivity. Also of interest is the extension of QEDFT approaches, including the linear-response technique presented here, to solid-state systems treated with periodic boundary conditions to study the effects of optical cavities on phonons and phonon-polaritons.[49]