Nonadiabatic phenomena in molecular vibrational polaritons

The hamiltonian of rotating-vibrating molecules interacting with a single lossless cavity mode can be written in the dipole approximation as

where $\hat{H}_{{\rm m}}^{(i)}$ is the ith field-free molecular rovibrational Hamiltonian, $\hat{a}_{c}^{\dagger}$ and $\hat{a}_{c}$ are photon creation and annihilation operators, respectively, $\omega_{c}$ is the frequency of the cavity mode, $\hslash$ is Planck’s constant divided by $2\pi$, $\varepsilon_{0}$ is the electric constant, $V$ is the volume of the electromagnetic mode, $\mathbf{e}$ is the polarization vector of the cavity mode, and $\hat{\mathbf{\upmu}}^{(i)}$ is the dipole moment of the ith molecule.

We assume the cavity electric field to be polarized along the lab-fixed z-axis, $\mathbf{e}=(0,0,1)$, and we consider an IR cavity, which allows restricting the molecular model to a single electronic state. Then, the

direct-product functions, also called diabatic states, provide a complete basis, where the $|\Psi_{{\rm rovib}}^{i,nJM}\rangle$ field-free rovibrational eigenstates satisfy

$J$ and $M$ are the rotational angular momentum and its projection onto the space-fixed z-axis, respectively, $n$ is all other quantum numbers uniquely defining the rovibrational states, and $|N\rangle$ is a photon number state of the cavity radiation. For a diatomic molecule $n$ can be identified with the single vibrational quantum number $v$, and using the notation $|\Psi_{{\rm rovib}}^{i,nJM}\rangle\equiv|v_{i}J_{i}M_{i}\rangle$, the matrix representation of Eq. (1) using the basis of Eq. (2) gives matrix elements of the following form:

where $\mu^{(i)}(R_{i})$ is the permanent dipole moment as a function of the internuclear distance $R_{i}$, $\theta_{i}$ is the angle between the space-fixed z axis and the ith molecule, and the energy scale was set so that the zero point energy of the photonic mode, $\hbar\omega_{{\rm c}}/2$, is omitted.

The rovibrational polaritonic states are the eigenstates of the Hamiltonian of Eq. (1). The kth polariton can be written as

where $C_{N,v_{1}J_{1}M_{1},...,v_{N_{{\rm mol}}}J_{N_{{\rm mol}}}M_{N_{{\rm mol}}}}^{k}$ are the coefficients of the kth eigenvector of the matrix having the elements given in Eq. (4).

The field-free rovibrational energies of the H³⁵Cl and H³⁷Cl molecules were computed as $E_{{\rm rovib}}^{(i),v_{i}J_{i}}=E_{{\rm vib}}^{(i),v_{i}}+B_{v_{i}}^{(i)}J(J+1)$, where the $E_{{\rm vib}}^{(i),v_{i}}$ vibrational energies were obtained from converged variational computations utilizing the spherical-DVR discrete variable representation (DVR) Szidarovszky et al. (2010) and the potential energy curve (PEC) of Ref. Magnier et al. (1993). The $B_{v}$ rotational constants (given in energy units) were evaluated for each vibrational state as $B_{v}^{(i)}=\langle v_{i}|\hbar^{2}/(2m_{{\rm red}}^{(i)}R^{2})|v_{i}\rangle$, where $m_{{\rm red}}^{(i)}$ is the reduced vibrational mass of the ith molecule. The nuclear masses used in the computations were $m_{{\rm H}}=1.00784$ u, $m_{{\rm{}^{35}Cl}}=34.9689$ u, and $m_{{\rm{}^{37}Cl}}=36.9659$ u.

For computing the $\langle v_{i}|\mu^{(i)}(R)|v^{\prime}_{i}\rangle\equiv\mu^{(i)}_{v_{i}v^{\prime}_{i}}$ vibrational (transition) dipole elements, the dipole moment curve of Ref. Zemke et al. (1981) was used. The values obtained are $\mu^{\rm(H^{35}Cl)}_{00}=0.42858$, $\mu^{\rm(H^{35}Cl)}_{10}=\mu^{\rm(H^{35}Cl)}_{01}=0.027055$, $\mu^{\rm(H^{37}Cl)}_{00}=0.42857$, and $\mu^{\rm(H^{37}Cl)}_{10}=\mu^{\rm(H^{37}Cl)}_{01}=0.027045$, all in atomic units. The $\langle J_{i}M_{i}|{\rm cos}(\theta)|J^{\prime}_{i}M_{i}\rangle$ terms were evaluated using an expression involving 3-j symbols,

which can be derived using the relations $\langle\phi,\theta|J\,M\rangle=Y_{JM}(\theta,\phi)=\sqrt{\frac{2J+1}{4\pi}}{D_{M0}^{J}}^{*}(\phi,\theta,\chi)$, ${\rm cos}(\theta)={D_{00}^{1}}(\phi,\theta,\chi)$, and $\int D_{K_{1}M_{1}}^{J_{1}}(\mathbf{\Omega})D_{K_{2}M_{2}}^{J_{2}}(\mathbf{\Omega})D_{K_{3}M_{3}}^{J_{3}}(\mathbf{\Omega})d\mathbf{\Omega}=8\pi^{2}\footnotesize{\begin{pmatrix}J_{1}&J_{2}&J_{3}\\ K_{1}&K_{2}&K_{3}\end{pmatrix}}\footnotesize{\begin{pmatrix}J_{1}&J_{2}&J_{3}\\ M_{1}&M_{2}&M_{3}\end{pmatrix}}$ Zare (1988), where $D_{KM}^{J}(\mathbf{\Omega})$ is the Wigner D-matrix and $\mathbf{\Omega}$ represents the three Euler angles $\phi$, $\theta$, and $\chi$.

For understanding the underlying physics, it is useful to construct the matrix representation of the Hamiltonian in Eq. (1) using vibrational-photonic basis functions $|N\rangle\prod_{i=1}^{N_{{\rm mol}}}|\Psi_{{\rm vib}}^{i,v_{i}}\rangle$. The matrix elements of the resulting vibro-photonic Hamiltonian depend on the $\theta$ rotational coordinates. After removing the rotational kinetic energy terms, the eigenvalues of the vibro-photonic Hamiltonian form vibrational polaritonic energy surfaces (VPES), which are functions of the $N_{{\rm mol}}$ rotational coordinates. Note the analogy that the VPESs provide effective potential energy surfaces (PES) for the rotational motion the same way as the usual light-dressed electronic potentials provide effective PESs for vibrational motion in polyatomics.

For two diatomic molecules (denoted as molecule 1 and 2, respectively) coupled to a single IR cavity mode, the first 4$\times$4 block in the vibro-photonic Hamiltonian reads

where

where $\mu_{kl}^{(i)}=\langle v_{i}=k|\mu^{(i)}(R)|v_{i}=l\rangle$, $\hat{T}_{{\rm rot}}^{(i)}$ is the squared rotational angular momentum of the ith molecule, and $B_{v}^{(i)}$ is the rotational constant of the ith molecule in its vth vibrational state. For the specific case of a H³⁵Cl and a H³⁷Cl molecule, Figure 1 show the three largest eigenvalues of the $\mathbf{V}_{{\rm rot}}$ matrix in Eq. (9) as a function of the $\theta_{i}$ variables for a coupling strength of $\mu^{\rm(H^{35}Cl)}_{00}\sqrt{\hbar\omega_{{\rm c}}/(\varepsilon_{0}V)}=33.26$ cm^-1. The same is shown in Fig. 2 for the case of two identical H³⁵Cl molecules. Fig. 1 demonstrates that in the case of the two different isotopologues, three polaritons are formed, as expected Houdré et al. (1996); Szidarovszky et al. (2020), and conical intersections (CI) can be formed between the VPESs. In the case of two identical H³⁵Cl molecules, Fig. 2 shows that the middle polaritonic surface is flat, which one would usually identify as a dark polaritonic state (although it is not dark in this case, see below). The upper (lower) VPES touches the dark state for red (blue) shifted cavity wave lengths.

A well-known consequence of CIs or LICIs between adiabatic electronic PESs is that the electronic wave functions of the connected surfaces become double valued, i.e., propagating the electronic wave function of a single surface along a closed curve around the CI leads to the accumulation of a topological phase of $\pi$ Longuet-Higgins et al. (1958); Longuet-Higgins (1961); Mead and Truhlar (1982); Berry (1984); Halász et al. (2011, 2018). In order to verify that the LICIs between VPESs are indeed true LICIs (and to confirm the existence of LICIs between VPESs even when the diagonal dipole terms are not omitted from the VPESs), the numerical verification of the topological phases was carried out. Now, we borrow a well-established framework developed for natural nonadiabatic phenomena, which was succesfully applied to study the topological properties of “natural” electronic degeneracies Baer (1975, 2002); Vértesi et al. (2005); Baer (2006).

Let us start from the adiabatic-to-diabatic transformation (ADT)

where $\mathbf{W}$ and $\mathbf{V}$ denote the diabatic and adiabatic PES, respectively. $\mathbf{A}$ is a unitary matrix called ADT matrix, which can be obtained by solving the equation

Its solution can be written in the following form

where $\gamma$ is a path ordering operator and $\bm{\tau}(s)$ is the nonadiabatic coupling matrix which contains the nonadibatic coupling terms (NAC)s

Here $\Psi_{i}$ and $\Psi_{j}$ are $ith$ and $jth$ adiabatic vibrational wave functions, respectively and $\nabla$ is the gradient operator corresponding to the rotational coordinates.

The necessary condition for the $\mathbf{A}$-matrix to yield single-valued diabatic wave function is that the $\mathbf{D}$-matrix,

needs to be diagonal and has, in its diagonal, numbers of ($+1$)s or ($-1$)s. The number of ($-1$)s in the $\mathbf{D}$ is designated as $K$ which is the number of the adiabatic vibrational eigenfunctions that flip their sign. The $\mathbf{D}$-matrix is the topological matrix which contains all topological features of a vibrational manifold in a region surrounded by its contour $\varGamma$ and $K$ is the topological number. The positions of the ($-1$)s in the $\mathbf{D}$-matrix correspond with the vibrational eigenfunctions that flip their sign. The equation for $\mathbf{D}(\Gamma)$ is a kind of quantization condition to be fulfilled by the nonadiabatic coupling matrix $\bm{\tau}$ along any closed contour $\varGamma$ in configuration space (CP).

In case of two quasi-isolated adiabatic states ($j,j+1$) the quantization condition reads:

where $\alpha_{j,j+1}(\Gamma)$ is the topological phase. When the closed contour $\Gamma$ is a circle with radii q:

The adiabatic-to-diabatic transformation angle can be defined as:

In the related numerical calculations of this study we concentrate on LICIs between the lower and middle as well as the middle and upper vibrational polaritonic states (see Figs. 1 and 2). In Figure 1 these are the two (1,2) and two (2,3) LICIs, corresponding to the three bright states, which are formed for the “H³⁵Cl + H³⁷Cl + IR cavity mode” system. These are “true” LICIs, i.e., first-order degeneracy points, as the two HCl molecules are distinguishable. However, for the two indistinguishable confined H³⁵Cl molecules, only a second-order degeneracy point is formed, as demonstrated in Fig. 2.

In Fig. 3 we report on results for the “H³⁵Cl + H³⁷Cl + IR cavity mode” system. At first we calculated the $\bm{\tau}_{1,2}$ along a circle located at the (1,2)LICI point. In Fig. 3. panel A the related ADT angle $\gamma_{1,2}$ and the corresponding topological phase $\alpha_{1,2}=-1.0169\pi$ are shown. Similar quantities are presented for the (2,3)LICI in panel B. Here we obtain $\alpha_{2,3}=1.0077\pi$. It is well noticed that both $\alpha_{1,2}$ and $\alpha_{2,3}$ are close to $\pi$ as expected. This means that there is only one (1,2)LICI and one (2,3)LICI in the studied CS region of panel A and B, respectively, and the two-state approximation works well. In panel C the closed contour, in contrast to the two previous ones, surrounds two LICI points. The calculated $\alpha_{1,2}=-1.1227\pi$ and $\alpha_{2,3}=-0.3436\pi$ topological phases are far from being close to $\pi$, indicating that the two-state approximation fails to describe accurately the topological behavior. The third state is needed to be involved in the calculations for the proper topological description. In order to derive the three-state ADT matrix elements and the corresponding $\mathbf{D}$ matrix one has to calculate all off-diagonal matrix elements of the $\bm{\tau}(s)$ NAC matrix and then apply Eqs. 12 and 14. On panel D of Fig. 3 the values for the $\mathbf{D}$ matrix elements are presented for the closed curve indicated as an inset in the panel. The values of the diagonal elements ($-1$,$+1$,$-1$) clearly imply that the distribution contains an odd number (presumably one) of LICI for each type (1,2)LICI and (2,3)LICI. Therefore, the first and third vibrational wave functions change their sign while the sign of the second vibrational wave function remains unchanged after flipping its sign twice.

In Fig. 4 topological results are presented for the “2$\times$H³⁵Cl + IR cavity mode” system in which, due to the symmetry, only a second-order degeneracy point is formed (see on Fig. 2). Consequently, if the two-state approximation worked, one would have to obtain $2\pi$ for the value of the topological phase $\alpha_{1,2}$. Instead, we obtained values of $\alpha_{1,2}=1.0878\pi$ and $\alpha_{2,3}=1.6789\pi$, which clearly show that the two-state approximation does not work in the studied region of the CS. At this point, one has to derive again the three-state ADT matrix elements and the corresponding $\mathbf{D}$ matrix by employing the appropriate nonadiabatic coupling terms. The obtained ($1$,$1$,$1$) values in the diagonal of the $\mathbf{D}$ matrix clearly demonstrate that the three-state approximation works perfectly and none of the vibrational wave functions change sign.

The above studies undoubtedly proved that the topological features of VPES and electronic PES are very similar. Earlier results Baer (2002); Vértesi et al. (2005); Baer (2006) obtained for the latter situation can almost be transferred one by one so as to understand the topological features of the VPESs. The topological properties of the VPESs imply that direct observables should exist, which prove the existence of the nonadiabatic effects among the VPESs. Some of these are presented in the following sections.

Another physical property which can be effected by CIs or LICIs is the spectrum. The spectroscopy of polaritonic systems is a valuable tool for providing insight into their various physiochemical properties Zeb et al. (2018); Ribeiro et al. (2018b); Szidarovszky et al. (2018a, 2020), and nonadiabatic effects resulting from the strong coupling of the electronic, vibrational, rotational, and photonic degrees of freedom can cause a significant modification on the molecular absorption spectra. In one possible theoretical formulation of light-dressed spectroscopy Szidarovszky et al. (2019), transitions between light-dressed states are evaluated using first-order time-dependent perturbation theory, that is, for the $T_{l\leftarrow k}$ transition amplitude between the kth and lth light-dressed state $T_{l\leftarrow k}\propto\langle\Psi_{l}|\hat{\mathbf{d}}\mathbf{e}^{(p)}|\Psi_{k}\rangle\delta(E_{l}-E_{k}\pm\hbar\omega_{p})$ is assumed, where $\mathbf{e}^{(p)}$ and $\hbar\omega_{p}$ are the polarization vector and photon energy of the probe pulse, respectively, while $|\Psi_{i}\rangle$ and $E_{i}$ are the light-dressed wave functions and energies, respectively. In this work the light-dressed states (polaritonic states), $|\Psi_{i}\rangle$, have the form shown in Eq. (5), and the $T_{l\leftarrow k}$ transition amplitudes can be computed with formulae similar to that in Eq. (4). Note that within this framework the computed spectra represent the absorption of light running parallel to the cavity mirrors, in the medium, not that of which is transmitted through the mirrors.

Figure 5 shows the IR absorption spectra of H³⁵Cl molecules and the mixture of H³⁵Cl and H³⁷Cl molecules, confined in IR cavities of different coupling strengths. Blue, red, and green colors in the transition lines represent vibrational excitation in H³⁵Cl, vibrational excitation in the other H³⁵Cl (left panels) or H³⁷Cl (right panels), and photonic excitation, respectively. The color of each transition line is obtained by mixing the above three colors according to the weight of the different types of basis functions in the final polaritonic states of the transition. The panel (A1) in Fig. 5 demonstrates that at the weakest coupling strength presented in this work, the light-dressed spectrum of two confined H³⁵Cl molecules shows one dominant and several minor peaks, arising from the splitting of the single $|v=1,\,J=1,\,M=0\rangle\leftarrow|v=0,\,J=0,\,M=0\rangle$ field-free transition. The strongest peak is split into two in panel (B1) of Fig. 5. The small energy splitting reflects the slight differences in the rovibrational energies of H³⁵Cl and H³⁷Cl, while the significant difference in the peak heights reflects the different VPES landscape and related nonadiabatic effects (under field-free conditions the height of the split peak is nearly identical). The small green peaks represent excitation in the photonic mode (photon number increased by one), which appears in the spectrum due to the contamination of the photonic mode with molecular excited (ro)vibrational states, having non-zero transition amplitude. No significant difference can be seen in the low-resolution envelopes of the two panels in the first row.

The light dressed spectra at larger cavity coupling strengths demonstrate that the individual peak positions and peak heights can be slightly different for the confined 2$\times$H³⁵Cl and H³⁵Cl+H³⁷Cl systems, but no significant differences arise in their spectrum envelope. Due to the strong mixing between the different molecular excitations and the photonic mode, assigning the peaks at larger cavity coupling strengths becomes complicated and less informative. Nonetheless, the peak colors represent the character of the different transitions to some extent. Panels (B2),(B3) and (B4) of Fig. 5 show that some peaks are dominated by a H³⁵Cl, H³⁷Cl or photonic transition (peaks with mostly red, blue or green color), but other peaks show significant mixing between the rovibrational modes of the HCl molecules and the photonic mode (purple, dark green or brown peaks). Despite the subtle differences in the spectra of the “2$\times$H³⁵Cl + cavity mode” and “H³⁵Cl + H³⁷Cl cavity mode” systems, some of which can be attributed to the nonadiabatic couplings between the VPESs, the low-resolution spectrum envelopes remain nearly identical, thus the IR spectrum doesn’t seem to be the most efficient tool to highlight the nonadiabatic effects in this case.

An additional approach to highlight nonadiabatic couplings in molecular systems is to investigate their laser-induced dynamics. In this work such simulations were carried out by directly solving the time-dependent Schrödinger equation in the diabatic direct-product basis representation, assuming $T=0$ K. At room temperature, due to the non-zero populations in $J>0$ rotational levels, the numerical values for the expectation values shown below would vary, but the qualitative conclusions should not change.

As already demonstrated for the vibrational dynamics on electronic PESs Halász et al. (2011), an efficient way to pinpoint nonadiabatic effects is to monitor the populations on the adiabatic surfaces during the dynamics. Here we perform similar studies for the case of adiabatic vibrational polaritons. In the case of VPESs, the population on the $a$th adiabatic state is evaluated as

where $|\Psi(t)\rangle$ is the time-dependent wave function expanded in terms of the direct-product diabatic states, see Eq. (2), and $\hat{P}^{(a)}=|\psi^{(a)}\rangle\langle\psi^{(a)}|$ is a projector onto the $a$th adiabatic state, i.e., the $a$th eigenstate of the matrix in Eq. (9). Note that

that is, the $C_{N,v_{1},v_{2}}^{(a)}(\theta_{1},\theta_{2})$ expansion coefficients depend on the rotational coordinates, $\theta_{1}$ and $\theta_{2}$, as predictable from Figs. 1 and 2.

Fig. 6 shows the temporal evolution of the adiabatic populations for the “2$\times$H³⁵Cl + IR cavity mode” system and the “H³⁵Cl + H³⁷Cl + IR cavity mode” system, when a short IR laser pulse, tuned to the field-free vibrational fundamental of H³⁵Cl (2926.1 cm^-1), is used to excite these systems from their ground states. After 60 fs, when the laser field has subsided, the population of the ground adiabatic state is 0.52 for both systems, indicating that the degree of laser-induced molecular vibrational excitation is identical for the two systems. For both the “2$\times$H³⁵Cl + IR cavity mode” and “H³⁵Cl + H³⁷Cl + IR cavity mode” systems, shown in the the upper and lower panels of Fig. 6, respectively, the laser populates all three excited adiabatic states. However, the population is distributed among the three excited adiabatic states in a different manner for the two systems. Somewhat surprisingly, the middle polaritons, which one would intuitively address as dark states, are also populated by the light field for both systems. This is due to the permanent dipole coupling with the ground adiabatic state (see first row and column of second matrix in Eq. (9)), test simulations without this coupling show that the middle polaritons are indeed almost “dark”, and are much less populated. In either case, due to the strong nonadiabatic couplings between the VPESs, rapid, hundred-femtosecond timescale population transfer proceeds among all three adiabatic states following the excitation by the laser field. Overall, Fig. 6 demonstrates that when molecular rotations are feasable, strong nonadiabatic couplings influence the dynamics of multimolecule (ro)vibrational polaritons. Furthermore, the nonadiabatic dynamics shown in Fig. 6 verify that breaking the molecular permutational symmetry, in this case by introducing a different Cl isotope, the polariton dynamics can deviate substantially.

As an alternative approach to probe the presence of nonadiabatic couplings in the laser-induced dynamics of our rovibrational polaritonic test systems, the orientation, i.e., the $\langle{\rm cos}(\theta)\rangle$ expectation value, was evaluated for a single H³⁵Cl (H³⁷Cl) molecule, when confined in an IR cavity with another H³⁵Cl (H³⁷Cl) molecule or a H³⁷Cl (H³⁵Cl) isotopologue. The numerical results are shown for the same and mixed isotopologue cases in the upper and lower panels of Figure 7, respectively. The upper panel of Figure 7 demonstrates that the isotope effect causes only minor differences in the orientation dinamics, the curves for H³⁵Cl and H³⁷Cl and nearly identical. However, the lower panel of Figure 7 shows that on the timescale of the rotational dynamics of HCl (picoseconds), the orientation curves deviate substantially when the two isotopologues are mixed in the cavity. Although the exact value of the numerical results presented here should be treated with caution, because cavity loss was neglected in this work, such a dependence of the orientation on the partner molecule is a clear fingerprint of molecular distinguishability on the VPESs (and corresponding nonadiabatic couplings) and, therefore, on the rotational dynamics.