# Cavity-induced Non-Adiabatic Dynamics and Spectroscopy of Molecular
Rovibrational Polaritons studied by Multi-Mode Quantum Models

Eric W. Fischer<EMAIL_ADDRESS>Theoretische Chemie, Institut für
Chemie, Universität Potsdam, Karl-Liebknecht-Straße 24-25, D-14476 Potsdam-
Golm, Germany  Peter Saalfrank<EMAIL_ADDRESS>Theoretische
Chemie, Institut für Chemie, Universität Potsdam, Karl-Liebknecht-Straße
24-25, D-14476 Potsdam-Golm, Germany  Institut für Physik und Astronomie,
Universität Potsdam, Karl-Liebknecht-Straße 24-25, D-14476 Potsdam-Golm,
Germany

###### Abstract

We study theoretically the quantum dynamics and spectroscopy of rovibrational
polaritons formed in a model system composed of a single rovibrating diatomic
molecule, which interacts with two degenerate, orthogonally polarized modes of
an optical Fabry-Pérot cavity. We employ an effective rovibrational Pauli-
Fierz Hamiltonian in length gauge representation and identify three-state
vibro-polaritonic conical intersections (VPCIs) between singly-excited vibro-
polaritonic states in a two-dimensional angular coordinate branching space.
The lower and upper vibrational polaritons are of mixed light-matter hybrid
character, whereas the intermediate state is purely photonic in nature. The
VPCIs provide effective population transfer channels between singly-excited
vibrational polaritons, which manifest in rich interference patterns in
rotational densities. Spectroscopically, three bright singly-excited states
are identified, when an external infrared laser field couples to both a
molecular and a cavity mode. The non-trivial VPCI topology manifests as
pronounced multi-peak progression in the spectral region of the upper
vibrational polariton, which is traced back to the emergence of rovibro-
polaritonic light-matter hybrid states. Experimentally ubiquitous spontaneous
emission from cavity modes induces a dissipative reduction of intensity and
peak broadening, which mainly influences the purely photonic intermediate
state peak as well as the rovibro-polaritonic progression.

##  I Introduction

The interaction of infrared-active molecular vibrations with quantized
electromagnetic field modes of optical Fabry-Pérot-type cavities has been
experimentally shown to have a peculiar impact on chemical ground state
reactivity and molecular propertiesgeorge2015 ; ebbesen2016 ; thomas2016 ;
thomas2019 ; herrera2020 . The origin of this impact is attributed to the
formation of vibrational light-matter hybrid states, known as vibrational
polaritonsebbesen2016 . Linearthomas2016 ; thomas2019 ; chervy2018  and
nonlinearxiang2018 ; zhang2019  infrared spectroscopic techniques have been
the experimental methods of choice to characterize these vibrational light-
matter hybrid states. Consequently, there has been a decent amount of
theoretical effort to address the calculation of infrared spectra for vibro-
polaritonic systemssaurabh2016 ; ribeiro2018 ; flick2019 ; li2020 ; li2021 ;
fischer2021 ; bonini2021 . There, one or more molecular vibrations interact
with a single or multiple quantized cavity modes and the composite system is
commonly described by an effective vibrational Pauli-Fierz Hamiltonian in the
length-gauge representation of non-relativistic cavity quantum electrodynamics
(cQED)flick2017a ; flick2017b ; schaefer2018 ; lihuo2021a . In the latter, the
light-matter interaction is mediated by the molecular ground state dipole
moment and its orientation relative to the cavity mode polarization
directions. Often, the dipole moment is assumed to be either exactly aligned
with the polarization direction of a cavity mode or fixed at a certain
orientation with respect to a single cavity mode polarization. The impact of
molecular rotations on vibrational polaritons, which naturally influences the
light-matter coupling due to the orientation dependent nature of the
interaction, has been addressed only recently in a few worksflick2018 ;
sidler2020 ; szidarovszky2021 ; triana2021 . There, the orientation of the
dipole moment was allowed to vary with respect to the polarization direction
of a single cavity mode and rotation related cavity-induced conical
intersections have been studied in comparison to the role of classical laser
fieldstriana2021 . However, the quantization of the cavity radiation field
naturally leads to doubly degenerate optical modes with two orthogonal
polarization directionsthiru1998 , which both should be relevant for the
light-matter interaction of rovibrating molecules.

Here, we consider a rovibrating diatomic molecule in an optical cavity, which
interacts with a pair of degenerate and orthogonally polarized cavity modes.
Inspired by recent work of Vibók and coworkersszidarovszky2021 , we treat the
rovibro-polaritonic problem from a non-adiabatic perspective by energetically
separating the rotational and vibro-polaritonic dynamics. As a result, we
observe non-adiabatic effects that manifest as three-state vibro-polaritonic
conical intersections (VPCIs) between vibro-polaritonic surfaces. We discuss
the three-state VPCIs in detail from a topological, a dynamical and a
spectroscopical perspective. Further, we examine the effect of experimentally
ubiquitous dissipation on the infrared spectra of rovibrating molecules in an
optical cavity due to spontaneous emission from excited states. We note the
model character of our study, as gas phase experiments in optical cavities
have not yet been performed due to issues in reaching the vibrational strong
coupling regime.

This paper is structured as follows: In Sec.II, we introduce the effective
rovibrational Pauli-Fierz Hamiltonian (II.1), our ab initio molecular model
(II.2) and discuss the adiabatic separation of rotational and vibro-
polaritonic degrees of freedom (II.3). In Sec.III, we present and discuss our
results. First, we examine the formation of three-state vibro-polaritonic
conical intersections and their topological properties (III.1). Second, we
discuss vibro-polaritonic population dynamics and the time-evolution of
reduced vibro-polaritonic rotational densities (III.2). Third, we characterize
the light-matter hybrid system by means of infrared spectra (III.3) and
discuss the spectroscopic manifestations of non-adiabatic features in
comparison to purely vibro-polaritonic systems (III.3.1). Finally, we study
the influence of dissipation on the infrared spectra due to spontaneous
emission effects(III.3.2). Sec.IV summarizes our work and connects our
findings to other branches of molecular polaritonics.

##  II Theory and Model

###  II.1 Rovibrational Pauli-Fierz Hamiltonian

We consider a diatomic molecule in a two-mode optical cavity (cf. Fig.1(a))
described by an effective length-gauge, rovibrational Pauli-Fierz Hamiltonian
in dipole and cavity Born-Oppenheimer type approximationfischer2021 ;
flick2017a ; flick2017b ; schaefer2018 , as given by

$\hat{H}=\hat{H}_{S}+\hat{H}_{C}+\hat{H}_{SC}+\hat{H}_{DSE}\quad.$ (1)

In the specific example chosen in this work, $\hat{H}_{S}$ is the Hamiltonian
of a single rovibrating CO molecule

$\hat{H}_{S}=\dfrac{\hat{j}^{2}}{2I}\,\underbrace{-\dfrac{\hbar^{2}}{2\mu}\dfrac{\partial^{2}}{\partial
r^{2}}+V(r)}_{=\hat{H}_{\mathrm{vib}}}\quad,$ (2)

with moment of inertia, $I=\mu r^{2}$, where $r$ is the CO-stretching
coordinate and $\mu$ the reduced mass, and angular momentum operator

$\hat{j}^{2}=-\hbar^{2}\dfrac{1}{\sin\theta}\dfrac{\partial}{\partial\theta}\,\sin\theta\dfrac{\partial}{\partial\theta}+\dfrac{1}{\sin^{2}\theta}\dfrac{\partial^{2}}{\partial\phi^{2}}\quad,$
(3)

with polar angle, $\theta\in[0,\pi]$, and azimuthal angle, $\phi\in[0,2\pi)$,
respectively. The vibrational Hamiltonian, $\hat{H}_{\mathrm{vib}}$, is
determined by the molecular Born-Oppenheimer ground state potential energy
surface, $V(r)$.

Further, we consider a Fabry-Pérot type cavity with two energetically
degenerate and orthogonally polarized, transverse cavity modes

$\hat{H}_{C}=\sum_{\lambda=z,y}\hbar\omega_{c}\left(\hat{a}^{\dagger}_{\lambda}\hat{a}_{\lambda}+\dfrac{1}{2}\right)\quad,$
(4)

with polarization index $\lambda$, harmonic cavity frequency $\omega_{c}$ and
photon creation/annihilation operators
$\hat{a}^{\dagger}_{\lambda}/\hat{a}_{\lambda}$, respectively. Here, we chose
the cavity modes to be polarized along the $z$\- and $y$-axis of the molecular
space-fixed frame (perpendicular to the wave vector, $\underline{k}$), as
shown in Fig.1(a). Further, the orientation dependent light-matter interaction
is mediated by the molecular dipole moment

$\underline{d}(r,\theta,\phi)=d(r)\begin{pmatrix}\sin\theta\cos\phi\vspace{0.1cm}\\\
\sin\theta\sin\phi\vspace{.1cm}\\\ \cos\theta\end{pmatrix}\quad,$ (5)

where $d(r)$ is the molecular dipole moment along the CO-bonding axis. The
light-matter interaction is given by the projection of
$\underline{d}(r,\theta,\phi)$ on the cavity mode polarization directions,
$\left(\underline{\epsilon}_{\lambda}\cdot\underline{d}(r,\theta,\phi)\right)$,
with $\lambda=z,y$. Accordingly, the bare cavity-molecule interaction
Hamiltonian reads

$\displaystyle\hat{H}_{SC}$
$\displaystyle=g\sum_{\lambda=z,y}\biggl{(}\underline{\epsilon}_{\lambda}\cdot\underline{d}(r,\theta,\phi)\biggr{)}\left(\hat{a}^{\dagger}_{\lambda}+\hat{a}_{\lambda}\right)\quad,\vspace{0.2cm}$
(6)
$\displaystyle=g\,d(r)\biggl{(}\left(\hat{a}^{\dagger}_{z}+\hat{a}_{z}\right)\cos\theta+\left(\hat{a}^{\dagger}_{y}+\hat{a}_{y}\right)\sin\theta\sin\phi\biggr{)}\quad,$

and the dipole self-energy, which has been thoroughly studied in
Refs.rokaj2018  andschaefer2020 , is given by

$\displaystyle\hat{H}_{DSE}$
$\displaystyle=\dfrac{g^{2}}{\hbar\omega_{c}}\sum_{\lambda=z,y}\biggl{(}\underline{\epsilon}_{\lambda}\cdot\underline{d}(r,\theta,\phi)\biggr{)}^{2}\quad,\vspace{0.2cm}$
(7)
$\displaystyle=\dfrac{g^{2}}{\hbar\omega_{c}}d^{2}(r)\biggl{(}\cos^{2}\theta+\sin^{2}\theta\sin^{2}\phi\biggr{)}\quad.$

Note, there are two independent contributions to the DSE term related to
independent z- and y-polarized cavity modes without any cross terms coupling
modes along different polarization directions.

The light-matter interaction strength, $g$, in Eqs.(6) and (7), carrying the
dimension of an electric field strength (V/m), is related to a dimensionless
parameter $\eta$ viakockum2019

$g=\dfrac{\hbar\omega_{c}}{|d_{fi}|}\,\eta\quad.$ (8)

Here, $d_{fi}$ is a vibrational transition dipole matrix element between
selected initial and final molecular states (see below), to which we choose
the cavity mode frequency $\omega_{c}$ to be resonant. Generally, the coupling
parameter $g$ (or $\eta$) is a function of the cavity volume, the dielectric
constant inside the cavity and the number of molecules in the
cavity.fischer2021  Here we treat $g$ (and $\eta$) as an adjustable parameter.
Further, $\eta$ specifies the light-matter interaction regime: Vibrational
strong coupling (VSC) for $0.0<\eta<0.1$ and vibrational ultrastrong coupling
(VUSC) for $\eta\geq 0.1$, respectively.kockum2019

Figure 1: (a) Schematic sketch of diatomic molecule with vibrational
coordinate, $r$, and angular coordinates $\theta,\phi$ in optical two-mode
cavity with polarization vectors
$\underline{\epsilon}_{z},\underline{\epsilon}_{y}$, wavevector,
$\underline{k}$, and molecular space-fixed frame with axes, $x,y,z$. A
$z$-polarized cavity mode is indicated in yellow. (b) Anharmonic potential,
$V(r)$, and dipole function, $d(r)$, of CO stretching mode (interpolated
CCSD(T)/aug-cc-pV5Z results) with vibrational ground state, $\ket{0_{r}}$, and
first excited state, $\ket{1_{r}}$. $d(r)$ is given in units of Debye (Db).

In what follows, by assuming that the center of mass of the rovibrating CO is
fixed in the cavity, the full Hamiltonian is five-dimensional (three molecular
coordinates, $r$, $\theta$ and $\phi$, and two cavity modes). Below we shall
work in a mixed basis to represent $\hat{H}$, namely a grid basis for the
angular coordinates and a state-representation for the vibro-polaritonic modes
as characterized by quantum numbers $v_{r}$ (for the CO-vibration) and
($n_{z}$, $n_{y}$) for the cavity modes. Further, we emphasize the CO
molecule’s free rotation is described by two angular coordinates
$\theta,\phi$, i.e., the CO molecule is allowed to take arbitrary orientations
with respect to the cavity polarization plane spanned by vectors,
$\underline{\epsilon}_{z}$ and $\underline{\epsilon}_{y}$. Moreover, we do not
take into account the molecule’s center of mass motion as cavity boundary
effects are assumed to play no role under the long-wavelength approximation
and the neutral CO’s center of mass motion does not directly couple to the
cavity modes, which would require a net molecular chargesidler2020 .

The dynamics as generated by the rovibrational Pauli-Fierz Hamiltonian in
Eq.(1), is governed by a time-dependent Schrödinger equation

$\mathrm{i}\hbar\,\dfrac{\partial}{\partial
t}\Psi(v_{r},n_{z},n_{y},\theta,\phi,t)=\hat{H}\,\Psi(v_{r},n_{z},n_{y},\theta,\phi,t)\quad,$
(9)

with a five-dimensional rovibro-polaritonic wave packet,
$\Psi(v_{r},n_{z},n_{y},\theta,\phi,t)$. We will discuss two different initial
states to solve Eq.(9), as specified below, given by an isolated rovibrational
excitation of the molecule and a light-matter superposition state,
respectively.

###  II.2 Molecular Model

We study a single carbon monoxide molecule with reduced mass,
$\mu=m_{\mathrm{C}}m_{\mathrm{O}}/m_{\mathrm{CO}}=12506\,m_{e}$ (electron mass
$m_{e}$). The molecular ab initio PES, $V(r)$, and molecular dipole function,
$d(r)$, were calculated as function of the CO-bond-length $r$ (cf. Fig.1(b))
on the CCSD(T)/aug-cc-pV5Z level of theory via the software package
Gaussian16g16 . We find a ground state equilibrium bond length of
$r_{e}=2.145\,a_{0}$ and a rotational constant of $B=\frac{\hbar}{4\pi c\mu
r^{2}_{e}}=1.91\,\mathrm{cm}^{-1}$, which is in close agreement with the
experimental value of $B_{\mathrm{exp}}=1.92\,\mathrm{cm}^{-1}$rank1965 .

Further, we numerically obtained the two lowest vibrational
eigenvalues/eigenstates (cf. Fig.1(b)) of the vibrational Hamiltonian,
$\hat{H}_{\mathrm{vib}}$, in Eq.(2) based on a Colbert-Miller discrete
variable representationcolbert1992  with $N_{r}=1501$ grid points. The
fundamental anharmonic transition frequency,
$\hbar\omega_{10}=2137\,\mathrm{cm}^{-1}$, between the two lowest lying
vibrational eigenstates of the CO stretching mode compares well to an
experimental value of $2143\,\mathrm{cm}^{-1}$camilde1996 . The corresponding
vibrational transition dipole moment takes $d_{10}=0.066$ ea0, which is in
agreement with Ref.huajun2015 . The molecular dipole function $d(r)$ changes
roughly linearly in the $r$-range shown in Fig.1(b) and takes an absolute
equilibrium value of $|d(r_{e})|=0.12\,\mathrm{Db}$, which is in close
agreement with literaturescuseria1991 . In the following, we set the cavity
mode frequency $\hbar\omega_{c}=\hbar\omega_{10}$ and the transition dipole
moment $d_{fi}=d_{10}$.

###  II.3 Diabatic Vibro-Polaritonic Basis

The rotational constant, $B$, and the fundamental vibrational transition
energy, $\hbar\omega_{10}$, set two different excitation energy scales, which
allow for adiabatic separation of vibrational (“fast”) and rotational (“slow”)
degrees of freedomszidarovszky2021 . Accordingly, we consider a restricted
basis of zero-order “vibro-polaritonic” states, $\ket{v_{r},n_{z},n_{y}}$,
which constitute eigenstates to $\hat{H}_{\mathrm{vib}}+\hat{H}_{C}$, and is
given by

$\begin{cases}X_{0}:&\ket{0_{r},0_{z},0_{y}}\vspace{0.2cm}\\\
X_{1}:&\ket{1_{r},0_{z},0_{y}},\ket{0_{r},1_{z},0_{y}},\ket{0_{r},0_{z},1_{y}}\end{cases}\quad,$
(10)

with one-dimensional ground state manifold, $X_{0}$, and three-dimensional
singly-excited state manifold, $X_{1}$. The $X_{1}$-manifold allows for the
description of the lowest lying excited vibro-polaritonic states in the VSC
regime. Higher-lying excited states are neglected here.

In the following, we denote the four basis states in Eq.(10) generically as
$\ket{D_{k}}$ and expand a rovibro-polaritonic wave packet as

$\Psi(v_{r},n_{z},n_{y},\theta,\phi,t)=\sum_{k}\varphi_{k}(\theta,\phi,t)\ket{D_{k}}\quad,$
(11)

with time-dependent, rotational wave packets, $\varphi_{k}(\theta,\phi,t)$.
Further, the matrix representation of the rovibrational Pauli-Fierz
Hamiltonian in the zero-order basis, which constitutes a matrix operator in
$(\theta,\phi)$-space, is given by

$\underline{\underline{H}}=\underline{\underline{T}}(\theta,\phi)+\underline{\underline{V}}(\theta,\phi)\quad,$
(12)

with ($4\times 4$)-blocked potential energy,
$\underline{\underline{V}}(\theta,\phi)$, and rotational kinetic energy
operator matrices, $\underline{\underline{T}}(\theta,\phi)$, respectively. The
latter is given by

$\displaystyle\underline{\underline{{T}}}(\theta,\phi)=\frac{\hat{j}^{2}}{2\mu}\begin{pmatrix}\langle
0_{r}|\frac{1}{r^{2}}|0_{r}\rangle&\langle
0_{r}|\frac{1}{r^{2}}|1_{r}\rangle&0&0\\\ \langle
1_{r}|\frac{1}{r^{2}}|0_{r}\rangle&\langle
1_{r}|\frac{1}{r^{2}}|1_{r}\rangle&0&0\\\ 0&0&\langle
0_{r}|\frac{1}{r^{2}}|0_{r}\rangle&0\\\ 0&0&0&\langle
0_{r}|\frac{1}{r^{2}}|0_{r}\rangle&\\\ \end{pmatrix}$ (13)

where $\langle\dots\rangle$ indicates integration with respect to the
$r$-coordinate. The zero-point energy shifted potential energy matrix is given
by

$\displaystyle\underline{\underline{V}}(\theta,\phi)=\begin{pmatrix}\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{00}&\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{10}&g\,d_{00}\cos\theta&g\,d_{00}\sin\theta\sin\phi\\\
\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{10}&\hbar\omega_{10}+\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{11}&g\,d_{10}\cos\theta&g\,d_{10}\sin\theta\sin\phi\\\
g\,d_{00}\cos\theta&g\,d_{10}\cos\theta&\hbar\omega_{c}+\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{00}&0\\\
g\,d_{00}\sin\theta\sin\phi&g\,d_{10}\sin\theta\sin\phi&0&\hbar\omega_{c}+\dfrac{g^{2}}{\hbar\omega_{c}}f_{\theta\phi}\braket{d^{2}}_{00}\\\
\end{pmatrix}\quad,$ (14)

with dipole self-energy matrix elements,
$\braket{d^{2}}_{v_{r}v^{\prime}_{r}}=\braket{v_{r}}{d^{2}(r)}{v^{\prime}_{r}}$,
and

$f_{\theta\phi}=\cos^{2}\theta+\sin^{2}\theta\sin^{2}\phi\quad.$ (15)

In $\underline{\underline{V}}(\theta,\phi)$, the first diagonal entry
corresponds to the zero-order ground state and the diagonal elements of the
lower $3\times 3$-block resemble states forming the $X_{1}$-manifold in the
order given by Eq.(10). (The same assignment can be made for the kinetic-
energy matrix elements in Eq.(13).)

The $X_{0}$\- and $X_{1}$-manifolds are energetically well separated by
$\hbar\omega_{10}=\hbar\omega_{c}$ and subject to a weak anharmonic coupling
($g^{2}\braket{d^{2}}_{10}$) induced by the dipole self-energy and
comparatively strong light-matter coupling ($g\,d_{00}$). Further, the DSE
constitutes a potential in angular coordinates on the diagonal
($g^{2}\,\braket{d^{2}}_{v_{r}v_{r}}f_{\theta\phi}$) and the singly-excited
molecular state couples to the two singly-excited cavity mode states
($g\,d_{10}$).

There is only a weak kinetic coupling in the kinetic-energy matrix, by the
term $\frac{\hat{j}^{2}}{2\mu}\braket{0_{r}}{\frac{1}{r^{2}}}{1_{r}}$. Hence,
the Pauli-Fierz Hamiltonian takes a diabatic-like representation in the zero-
order basis Eq.(10). Following arguments given in Ref.szidarovszky2021 , we
interpret the eigenstates of $\underline{\underline{V}}(\theta,\phi)$ as
vibro-polaritonic states, which are given by the vibro-polaritonic ground
state $\ket{G}$, the lower $\ket{L_{1}}$, the middle $\ket{M_{1}}$ and the
upper $\ket{U_{1}}$ vibro-polaritonic states, respectively. Note, the latter
take the role of adiabatic states here. Further, the corresponding eigenvalues
$\varepsilon_{G}(\theta,\phi),\varepsilon_{L_{1}}(\theta,\phi),\varepsilon_{M_{1}}(\theta,\phi)$
and $\varepsilon_{U_{1}}(\theta,\phi)$, are functions of the angular
coordinates, $\theta$ and $\phi$, and constitute vibro-polaritonic potential
energy surfaces (PES) for the rotational dynamics of the molecule.

##  III Results and Discussion

###  III.1 Vibro-Polaritonic Conical Intersections

Figure 2: Top: Three-state vibro-polaritonic conical intersections (VPCIs)
between lower vibro-polaritonic, $\varepsilon_{L_{1}}(\theta,\phi)$, middle
vibro-polaritonic, $\varepsilon_{M_{1}}(\theta,\phi)$, and upper vibro-
polaritonic surfaces, $\varepsilon_{U_{1}}(\theta,\phi)$, for (a) VSC regime
with $\eta=0.05$ and (b) onset of VUSC regime with $\eta=0.1$. Bottom: One
dimensional cuts through vibro-polaritonic surfaces (bold), linear (dotted-
dashed) and quadratic approximations (dotted) with
$\varepsilon_{l}(\theta,\phi=0)$ and
$\varepsilon_{l}(\theta=\frac{\pi}{2},\phi)$ for (c) $\eta=0.05$ and (d)
$\eta=0.1$ with $l=L_{1},M_{1},U_{1}$.

We first discuss cavity-mode induced non-adiabatic effects in the manifold of
singly-excited vibro-polaritonic states. Under vibrational strong coupling, we
identify three-state vibro-polaritonic conical intersections (VPCIs) between
vibro-polaritonic surfaces $\varepsilon_{L_{1}}(\theta,\phi)$,
$\varepsilon_{M_{1}}(\theta,\phi)$ and $\varepsilon_{U_{1}}(\theta,\phi)$
under the conditions

$\hbar\omega_{10}=\hbar\omega_{c},\hskip 8.5359pt\theta=\dfrac{\pi}{2},\hskip
8.5359pt\phi=0,\pi\quad.$ (16)

The three-state VPCIs, as shown in Figs.2(a) and (b) for $\eta=0.05$ and
$\eta=0.1$, respectively, are located in a two-dimensional angular coordinate
branching space. Two distinct intersections are identified with intersection
coordinates $(\theta,\phi)=(\frac{\pi}{2},0)$ and $(\frac{\pi}{2},\pi)$, due
to the periodicity of the azimuthal angle, $\phi$.

The three-state VPCIs exhibit a characteristic double-cone topology formed by
the $L_{1}$\- and $U_{1}$-surfaces. They are triply degenerate at the
intersection point due to an additional crossing with the middle vibro-
polaritonic surface, $\varepsilon_{M_{1}}(\theta,\phi)$, which exhibits a
local minimum here. In the VSC regime ($\eta=0.05$), the splitting of the
$L_{1}$\- and $U_{1}$-surfaces is slightly asymmetric with respect to the
$M_{1}$-surface as can be seen in Figs.2(a) and (c). At the onset of the VUSC
regime ($\eta=0.1$), the $L_{1}$/$U_{1}$-splitting turns out to be strongly
asymmetric (cf. Figs.2(b) and (d)), with a dominant inverse cone in
$\varepsilon_{U_{1}}(\theta,\phi)$. The $M_{1}$-surface forms more pronounced
minima at the intersection points for stronger light-matter interaction and
the $L_{1}$-surface is subject to a mild “Mexican-hat”-type topology close to
the intersection coordinate.

To reveal detailed characteristics of the VPCI close to the intersection
coordinate, we expand $\underline{\underline{V}}(\theta,\phi)$ in
$\theta,\phi$ around the intersection point at
$(\theta,\phi)=(\frac{\pi}{2},0)$ (cf. Appendix B for details) in analogy to
vibronic problems. In Figs.2(c) and (d), we show cuts through the VPCIs along
$\theta$ and $\phi$, with exact surfaces (bold) besides linear (dotted-dashed)
and quadratic (dotted) approximations. We observe three characteristic
features of the VPCI: (i) the degeneracy is lifted linearly along both angular
coordinates in $\varepsilon_{L_{1}}(\theta,\phi)$ and
$\varepsilon_{U_{1}}(\theta,\phi)$ in analogy to linear vibronic coupling
theorykoeppel1984 , (ii) $\varepsilon_{M_{1}}(\theta,\phi)$ exhibits a
harmonic character at the intersection point, which results from the potential
character of the dipole self-energy and (iii) the mild local “Mexican-
hat”-type topology of $\varepsilon_{L_{1}}(\theta,\phi)$ at $\eta=0.1$ also
stems from the DSE.

We close by noting, that the topological motive of a three-state VPCI is
similar to accidental three-state conical intersections in molecular vibronic
coupling theory.matsika2003 ; matsika2005 ; coe2005 ; matsika2011

###  III.2 Rovibro-Polaritonic Dynamics

We now turn to the dynamics of different zero-order rovibro-polaritonic wave
packets, which mimic a light-matter hybrid system initially excited by an
external classical, $z$-polarized laser field. The additional action of a
laser could be included in the Pauli-Fierz Hamiltonian by adding a term,
$-\underline{\mu}\cdot\underline{E}(t)$ (with classical laser field
$\underline{E}(t)$ and where $\underline{\mu}$ is a “dipole” function
depending on molecular and cavity degrees of freedomflick2019 ), but here we
simply assume that such an excitation already took place and appropriate
initial rovibrational(-polaritonic) states have been prepared. In order to
represent different laser-polariton coupling scenarios, or, implicitly,
molecule-cavity dipole functions, we consider two different initial states,
namely

$\displaystyle\Psi^{(r)}_{0}$
$\displaystyle=\ket{1_{r},0_{z},0_{y}}\,Y^{1}_{0}(\theta,\phi)\quad\mathrm{and}\vspace{0.2cm}$
(17) $\displaystyle\Psi^{(r,z)}_{0}$
$\displaystyle=\frac{1}{\sqrt{2}}\biggl{(}\ket{1_{r},0_{z},0_{y}}+\ket{0_{r},1_{z},0_{y}}\biggr{)}\,Y^{1}_{0}(\theta,\phi)\,,$
(18)

with first excited rotational state $Y^{1}_{0}(\theta,\phi)$, where we adopt
the notation, $Y^{j}_{m_{j}}(\theta,\phi)$, for spherical harmonics. Here,
$\Psi^{(r)}_{0}$ resembles a purely molecular rovibrational excitation and
$\Psi^{(r,z)}_{0}$ constitutes a superposition of a singly-excited molecular
and a singly-excited cavity mode state, respectively. The first case resembles
a scenario, where an external laser field had coupled exclusively to the
molecule, whereas the latter considers additionally a coupling of the laser-
field to the $z$-polarized cavity mode.rovibXite

Using these initital states, we discuss the time-evolution of the light-matter
hybrid system in terms of zero-order (diabatic) populations

$\displaystyle P^{(k)}_{\mathrm{dia}}(t)$
$\displaystyle=\int^{2\pi}_{0}\int^{\pi}_{0}|\varphi_{k}(\theta,\phi,t)|^{2}\,\sin\theta\,\mathrm{d}\theta\,\mathrm{d}\phi\quad,$
(19)

with $\varphi_{k}$ given in Eq.(11), and their vibro-polaritonic (adiabatic)
counterparts, $P^{(l)}_{\mathrm{ad}}(t)$, with $l=G,L_{1},M_{1},U_{1}$. The
latter arise from adiabatic wave packets,
$\varphi^{\mathrm{ad}}_{l}(\theta,\phi,t)$, by using an expansion analogous to
Eq.(11 but with adiabatic basis states $\ket{G}$, $\ket{L_{1}}$, $\ket{M_{1}}$
and $\ket{U_{1}}$. Further, we study the adiabatic rotational dynamics by
means of reduced rotational densities

$\rho^{\mathrm{ad}}_{l}(\theta,t)=\int^{2\pi}_{0}|\varphi^{\mathrm{ad}}_{l}(\theta,\phi,t)|^{2}\,\mathrm{d}\phi\quad,$
(20)

which resemble the reduced dynamics of rotational wave packets,
$\varphi^{\mathrm{ad}}_{l}(\theta,\phi,t)$, on the vibro-polaritonic surfaces
as depicted in Fig.2. We note, that dynamics along the $\theta$-coordinate
turns out to be particularly illustrative compared to the $\phi$-coordinate.

We solve the TDSE (9) with the multiconfigurational time-dependent Hartree
(MCTDH) method in its multi-set formulationmeyer2012  as implemented in the
Heidelberg MCTDH package mctdh2019 . Numerical details on the MCTDH method as
well as on evaluation of above mentioned populations and reduced densities are
provided in Appendix A.

####  III.2.1 Vibro-Polaritonic Population Dynamics

Before we discuss the time-evolution of populations, we examine zero-order
state contributions to stationary vibro-polaritonic states at $t=0$. For the
rovibrationally excited initial state, $\Psi^{(r)}_{0}$, we find only
$\ket{L_{1}}$ and $\ket{U_{1}}$ equally populated for both light-matter
interaction scenarios with $\eta=0.05$ and $\eta=0.1$ (cf. Figs.3(a) and (b)).

Figure 3: Population dynamics for different initial states under VSC. Zero-
order (diabatic), $P_{\mathrm{dia}}(t)$, (left) and vibro-polaritonic
(adiabatic), $P_{\mathrm{ad}}(t)$, (right) population dynamics for
rovibrational singly-excited initial state, $\Psi^{(r)}_{0}$, for (a)
$\eta=0.05$ and (b) $\eta=0.1$. Zero-order (left column) and vibro-polaritonic
(right column) population dynamics for singly-excited superposition state,
$\Psi^{(r,z)}_{0}$, for (c) $\eta=0.05$ and (d) $\eta=0.1$.

In particular, we find no population in the middle vibro-polaritonic state
$\ket{M_{1}}$, i.e., no molecular vibrationally excited state is contributing
here, which identifies $\ket{M_{1}}$ as purely photonic in character. Further,
for the superposition state, $\Psi^{(r,z)}_{0}$, we observe contributions from
all three vibro-polaritonic states, with equally populated $\ket{L_{1}}$ and
$\ket{U_{1}}$ constituting the dominant contribution opposed to $\ket{M_{1}}$
(cf. Figs.3(c) and (d)). Accordingly, $\ket{L_{1}}$ and $\ket{U_{1}}$ are
identified as “true” light-matter hybrid states containing both molecular and
photonic contributions.

Turning to the dynamics, we consider a total propagation time of
$t_{f}=3000\,\mathrm{fs}$ for zero-order (“diabatic”) populations,
$P_{\mathrm{dia}}(t)$, and their vibro-polaritonic (“adiabatic”) counterparts,
$P_{\mathrm{ad}}(t)$. We find the dynamics up to $t_{f}$ to be mainly
dominated by states spanning the $X_{1}$-manifold as shown in the left and
right columns of Fig.3, which results from the relatively large energetic
separation of $X_{0}$\- and $X_{1}$-manifolds.

For the zero-order (“diabatic”) dynamics initiated by $\Psi^{(r)}_{0}$ and
depicted in Figs.3(a) and (b) (left column), we observe a fast initial
population transfer to both singly-excited cavity mode states followed by
coherent exchange dynamics. The latter are dominated by the molecular and the
$z$-polarized excited states and subject to quantum beats with a period of
roughly $1000\,\mathrm{fs}$. For an increased light-matter interaction
strength ($\eta=0.1$), the quantum beat period shortens to roughly
$600\,\mathrm{fs}$ and the frequency of coherent zero-order population
transfer is strongly enhanced. In contrast, vibro-polaritonic (“adiabatic”)
dynamics (cf. Figs.3(a) and (b), right column) are dominated by slow
population transfer from $\ket{U_{1}}$ to $\ket{L_{1}}$, accompanied by a
gradual population of $\ket{M_{1}}$. For $\eta=0.05$ we observe a
characteristic maximum in the $\ket{L_{1}}$-population at around
$700\,\mathrm{fs}$, which is shifted to roughly $500\,\mathrm{fs}$ for
$\eta=0.1$. Later times are characterized by several decaying recurrences in
$\ket{U_{1}}$.

For the dynamics initiated by $\Psi^{(r,z)}_{0}$, we find a significantly
slower zero-order population transfer between the molecular and the
$z$-polarized excited mode states for both coupling scenarios (cf. Figs.3(c)
and (d), left column). From the vibro-polaritonic perspective, as depicted in
the right column of Figs.3(c) and (d), $\ket{U_{1}}$ is initially again
depopulated in favor of $\ket{L_{1}}$, accompanied by a slight population
increase in $\ket{M_{1}}$. The first maximum in $\ket{L_{1}}$ appears at same
times as observed above and for longer times, the population dynamics is less
structured here. We note, for both $\Psi^{(r)}_{0}$ and $\Psi^{(r,z)}_{0}$ the
zero-order ground state is weakly contributing at $\eta=0.1$ and exhibits a
strongly oscillatory dynamics, which is especially pronounced for
$\Psi^{(r,z)}_{0}$.

Finally, from this population-based perspective, we conclude that (i) only
$\ket{L_{1}}$ and $\ket{U_{1}}$ are light-matter hybrid states opposed to the
purely photonic $\ket{M_{1}}$ state and (ii) the inclusion of rotational
degrees of freedom provides a cavity-induced transfer channel between vibro-
polaritonic excited states. Notably, the latter is absent in purely
vibrational problems, where an additional bath is required to mediate
population transfer between vibro-polaritonic states.

####  III.2.2 Reduced Rotational Dynamics

We now turn to rotational dynamics of the CO molecule based on vibro-
polaritonic (“adiabatic”) reduced rotational densities,
$\rho^{\mathrm{ad}}_{l}(\theta,t)$ (cf. Eq.(20)), for different vibro-
polaritonic surfaces. In order to distinguish reduced densities for different
initial states, we introduce the notation $\rho^{(r)}_{l}(\theta,t)$ for
$\Psi^{(r)}_{0}$ and $\rho^{(r,z)}_{l}(\theta,t)$ for $\Psi^{(r,z)}_{0}$,
respectively.

In Fig.4, the time-evolution of $\rho^{(r)}_{l}(\theta,t)$ is shown for
excited-state surfaces $\varepsilon_{U_{1}}$ (top row), $\varepsilon_{M_{1}}$
(middle row) and $\varepsilon_{L_{1}}$ (bottom row), respectively, with
$\eta=0.05$ (left column) and $\eta=0.1$ (right column). The horizontal
$\theta$-axis runs from $\pi$ to $0$ (left to right) and the vertical time-
axis runs from $0\,\mathrm{fs}$ to $3000\,\mathrm{fs}$ (top to bottom). The
position of the intersection point at $\theta=\frac{\pi}{2}$ is marked by a
vertical blue line.

Figure 4: Time-evolution of adiabatic reduced rotational densities,
$\rho^{(r)}_{l}(\theta,t)$, for rovibrational singly-excited initial state,
$\Psi^{(r)}_{0}$, with vertical time-axis and horizontal $\theta$-axis on
upper, $\varepsilon_{U_{1}}(\theta,\phi)$ (top row), middle,
$\varepsilon_{M_{1}}(\theta,\phi)$ (middle row) and lower,
$\varepsilon_{L_{1}}(\theta,\phi)$ (bottom row), vibro-polaritonic potential
energy surfaces for $\eta=0.05$ (left column) and $\eta=0.1$ (right column).
Intersection at $\theta=\frac{\pi}{2}$ indicated by blue vertical line. Total
density, $\rho^{(r)}(\theta,t)=\sum_{l}\rho^{(r)}_{l}(\theta,t)$, with
$l=G,L_{1},M_{1},U_{1}$ normalized for fixed time $t$.

At $t=0$, both $L_{1}$\- and $U_{1}$-surfaces are nearly equivalently
populated (in line with Figs.3(a) and (b)), and $\rho^{(r)}_{L_{1}}(\theta,t)$
and $\rho^{(r)}_{U_{1}}(\theta,t)$ show two symmetric maxima along $\theta$
with respect to the intersection point.

In the VSC regime with $\eta=0.05$ (cf. Fig.4, left column), density on the
$U_{1}$-surface symmetrically approaches the VPCI and is transferred to the
$L_{1}$-surface with minor contributions on $M_{1}$. In parallel,
$\rho^{(r)}_{L_{1}}(\theta,t)$ tends away from the intersection point and
develops a rich interference pattern as time evolves. This pattern can be
rationalized in terms of rotational wave packet interference, where non-
adiabatic transfer between $U_{1}$\- and $L_{1}$-surfaces induces interference
between wave packets initially located on different surfaces. At later times,
several recurrences at the $U_{1}$-surface are observed due to
$\rho^{(r)}_{L_{1}}(\theta,t)$ reentering the intersection region. Notably,
significant density on $\varepsilon_{M_{1}}$ is mainly found in the harmonic
region close to the VPCI. For VUSC with $\eta=0.1$, the density transfer
between $L_{1}$\- and $U_{1}$-surfaces is enhanced but qualitatively
equivalent to the VSC regime.

Figure 5: Time-evolution of adiabatic reduced rotational densities,
$\rho^{(r,z)}_{l}(\theta,t)$, for singly-excited superposition state,
$\Psi^{(r,z)}_{0}$, with vertical time-axis and horizontal $\theta$-axis on
upper, $\varepsilon_{U_{1}}(\theta,\phi)$ (top row), middle,
$\varepsilon_{M_{1}}(\theta,\phi)$ (middle row) and lower,
$\varepsilon_{L_{1}}(\theta,\phi)$ (bottom row), vibro-polaritonic potential
energy surfaces for $\eta=0.05$ (left column) and $\eta=0.1$ (right column).
Intersection at $\theta=\frac{\pi}{2}$ indicated by blue vertical line. Total
density, $\rho^{(r,z)}(\theta,t)=\sum_{l}\rho^{(r,z)}_{l}(\theta,t)$, with
$l=G,L_{1},M_{1},U_{1}$ normalized for fixed time $t$.

Turning to the initial superposition state, $\Psi^{(r,z)}_{0}$, we find the
reduced vibro-polaritonic rotational density to be initially distributed over
all three vibro-polaritonic surfaces (cf. Fig.5) in line with the population
analysis. Further at $t=0$, densities $\rho^{(r,z)}_{L_{1}}(\theta,t)$ and
$\rho^{(r,z)}_{U_{1}}(\theta,t)$ show a highly asymmetric distribution with
respect to the VPCI, in contrast to the symmetric character of
$\rho^{(r,z)}_{M_{1}}(\theta,t)$.

As time-evolves, we again observe a non-adiabatic density transfer via the
VPCI between $U_{1}$\- and $L_{1}$-surfaces for both interaction regimes (cf.
Fig.5, left and right column). The transfer is enhanced for VUSC in line with
the population analysis given above and characterized by a rich interference
pattern. Further, the asymmetric character of the reduced densities on the
$U_{1}$\- and $L_{1}$-surfaces is clearly observable for the whole time-
interval studied and the $M_{1}$-surface is substantially stronger explored
for both coupling regimes. We attribute the latter effect to a finite initial
population of $\ket{M_{1}}$, as the non-adiabatic population transfer to
$\varepsilon_{M_{1}}$ is rather inefficient as seen before.

In summary, rotational densities allow to reveal the non-adiabatic character
of vibro-polaritonic population transfer, which is dominated by (i) a
funneling effect of the VPCI and (ii) interference of rotational wave packets
initially located on different vibro-polaritonic surfaces.

###  III.3 Rovibro-Polaritonic Infrared Spectra

We now consider the spectroscopic characterization of the rovibro-polaritonic
system with an emphasis on non-adiabatic signatures of the three-state VPCI.
Infrared (IR) spectra, $\sigma(\omega)$, are calculated as

$\sigma(\omega)=A\displaystyle\int^{\infty}_{0}C(t)\,e^{\mathrm{i}\omega\,t}\mathrm{d}t\quad,$
(21)

with an autocorrelation function, $C(t)$, and a constant here chosen as$A=1$.
The autocorrelation function $C(t)$ is defined as the overlap between the
initial states $\Psi(0)=\Psi_{0}^{(r)}$ and $\Psi_{0}^{(r,z)}$ as defined in
Eqs.(17) and (18) and the propagated state, $\Psi(t)$, under the influence of
the coupled Pauli-Fierz Hamiltonian. (Note that this is an
approximation/variant to the “usual” procedure in which an eigenstate of the
Hamiltonian is multiplied by the dipole function and then propagated in
time.heller ) We obtain two different spectra due to the two initial states,
or implicitly, molecule-cavity dipole functions $\underline{\mu}$ (see above).
Further details on the evaluation of $C(t)$ and $\sigma(\omega)$ can be found
in Appendix A.

Figure 6: Vibro-polaritonic infrared spectra for different initial states and
light-matter interaction parameters $\eta$ with single-mode limit initial
states, $\tilde{\psi}^{(r)}_{0}=\ket{1_{r},0_{z}}$, and,
$\tilde{\psi}^{(r,z)}_{0}=\frac{1}{\sqrt{2}}\left(\ket{1_{r},0_{z}}+\ket{0_{r},1_{z}}\right)$,
(grey curves), VPCI energy (grey vertical lines) and intensities,
$I^{(r)}\propto|\braket{\Psi^{(r)}_{0}}{\phi_{i}}|^{2}$, and,
$I^{(r,z)}\propto|\braket{\Psi^{(r,z)}_{0}}{\phi_{i}}|^{2}$. Top row:
Infrared-spectra for the vibrational strong coupling (VSC) regime with
$\eta=0.05$ for (a) molecular rovibrationally excited initial state,
$\Psi^{(r)}_{0}$, and (b) light-matter superposition state,
$\Psi^{(r,z)}_{0}$. Bottom row: (c) and (d) analogous to (a) and (b) in top
row for onset of vibrational ultrastrong strong coupling (VUSC) regime at
$\eta=0.1$.

All spectra shown below have been obtained for a total propagation time,
$t_{f}=6000\,\mathrm{fs}$. In order to reveal rotational effects and effects
of a second cavity mode, we consider IR spectra obtained for a purely vibro-
polaritonic system in the single-cavity-mode limit as a reference, i.e., a CO
molecule solely interacting with the $z$-polarized cavity-mode. Here, we take
into account initial states $\Psi(0)=\tilde{\psi}^{(r)}_{0}=\ket{1_{r},0_{z}}$
and
$\Psi(0)=\tilde{\psi}^{(r,z)}_{0}=\frac{1}{\sqrt{2}}\left(\ket{1_{r},0_{z}}+\ket{0_{r},1_{z}}\right)$.
In order to interpret spectroscopic signatures, we also consider contributions
of eigenstates, $\phi_{i}$, of the effective Hamiltonian
$\underline{\underline{H}}$ (cf. Eq.(12)), to initial states, as given by
intensities $I^{(r)}\propto|\braket{\Psi^{(r)}_{0}}{\phi_{i}}|^{2}$ and
$I^{(r,z)}\propto|\braket{\Psi^{(r,z)}_{0}}{\phi_{i}}|^{2}$, respectively
(details are provided in Appendix A).

####  III.3.1 Infrared Spectra

We first consider spectra in the VSC regime with $\eta=0.05$ and initial
states $\Psi^{(r)}_{0}$ and $\Psi^{(r,z)}_{0}$ as shown in Fig.6(a) and (b).
For both initial states, we observe a series of transitions with a dominant
$L_{1}$-peak at $2076\,\mathrm{cm}^{-1}$ below the intersection point energy
of $2137\,\mathrm{cm}^{-1}$. In the region of the $U_{1}$-surface, we find a
progression of five peaks between $2207\,\mathrm{cm}^{-1}$ and
$2309\,\mathrm{cm}^{-1}$ depending on the initial state with spacing
$17-33\,\mathrm{cm}^{-1}$, which decreases for increasing peak energy. The
peak intensity increases with energy with a significantly less intense high-
energy peak terminating the progression. Notably, the latter transitions do
not reflect purely rovibrational states, but as discussed below correspond to
rovibro-polaritonic hybrid states. Moreover, for $\Psi^{(r,z)}_{0}$ a
prominent peak is found at $2171\,\mathrm{cm}^{-1}$, which resembles an
excitation of the purely photonic intermediate polariton state and is absent
for $\Psi^{(r)}_{0}$.

Further, the single-cavity-mode limit reveals the well-known $L_{1}$\- and
$U_{1}$-peaks for $\tilde{\psi}^{(r)}_{0}$ with Rabi-splitting
$\Omega_{R}=214\,\mathrm{cm}^{-1}$, respectively (grey curves). For
$\tilde{\psi}^{(r,z)}_{0}$, the chosen linear combination exactly captures the
$U_{1}$-state such that only a single peak is observed here. Accordingly, the
inclusion of rotational effects leads to a much richer excitation spectrum
(colored curves) than accessible in the single-mode limit (grey curves).
Explicitly, non-adiabatic effects show up as multi-peak progression in the
spectral region of the upper vibro-polaritonic state, which is characterized
by the inverse cone topology of the $\varepsilon_{U_{1}}$-surface.

From an analysis of intensities, $I^{(r)}$ and $I^{(r,z)}$, we are able to
reveal the detailed character of observed excitations, which are hidden in
spectra due to finite peak widths resulting from finite propagation times of
rovibro-polaritonic wave functions. Both, $I^{(r)}$ and $I^{(r,z)}$ exhibit
similar dominant contributions to the light-matter hybrid $L_{1}$\- and
$U_{1}$-peaks, whereas only $I^{(r,z)}$ exhibits the purely photonic
$M_{1}$-peak as expected. Notably, $I^{(r,z)}$ shows additional contributions
from several energetically close lying states with lower intensity to
individual spectroscopically resolved peaks. We attribute a dominant
rotational light-matter hybrid character to the latter as they are absent for
the rovibrationally excited initial state, i.e., seem to be dominated by a
photonic contribution.

Turning to IR spectra at the onset of the VUSC regime with $\eta=0.1$, we
observe a broadening of the spectrum and additional peaks emerge as shown in
Fig.6(c) and (d). In the single-mode limit, we find a large Rabi-splitting of
$\Omega_{R}=433\,\mathrm{cm}^{-1}$. Further, the progression in the
$U_{1}$-spectral region related to rovibro-polaritonic states is now more
pronounced with seven peaks showing slightly increased spacing of
$29-49\,\mathrm{cm}^{-1}$, which again decreases for increasing peak energy.
The middle-polariton peak at $2171\,\mathrm{cm}^{-1}$ for the light-matter
excited initials state splits into three peaks at $2116,2239$ and
$2258\,\mathrm{cm}^{-1}$, respectively. As previously for VSC, we observe a
number of (potentially energetically close) rovibro-polaritonic hybrid states
with lower intensity, which contribute to peaks in the spectrum obtained from
$\Psi^{(r,z)}_{0}$ as indicated by $I^{(r,z)}$ and $I^{(r)}$.

In conclusion, IR spectra reveal the formation of light-matter hybrid states
with contributions from both molecular rotations and vibrations as well as
cavity mode excitations, i.e., rovibro-polaritonic states, with non-adiabatic
signatures prominently manifesting as rovibro-polaritonic progression in the
spectral region of the $U_{1}$-inverse cone.

####  III.3.2 Cavity Loss Effects

Finally, we consider the impact of experimentally ubiquitous spontaneous
emission effects from cavity modes on rovibro-polaritonic infrared spectra.
Spontaneous emission manifests in finite excited cavity mode lifetimes (Finite
lifetimes of molecular rovibrational states are neglected in what follows).
These are considered here by a phenomenological approach recently employed in
electronic strong coupling studiesulusoy2020 ; felicetti2020 ; antoniou2020 .
There, the cavity mode Hamiltonian, $\hat{H}_{C}$, is replaced by a non-
Hermitian operator

$\hat{H}^{(\kappa)}_{C}=\sum_{\lambda=z,y}\left(\hbar\omega_{c}-\mathrm{i}\,\frac{\kappa}{2}\right)\hat{a}^{\dagger}_{\lambda}\hat{a}_{\lambda}\quad.$
(22)

The imaginary contribution accounts for an effective cavity decay rate,
$\kappa$, which we here set to $\kappa=43\,\mathrm{cm}^{-1}$ to model an
infrared cavity with quality factor $Q=\frac{\hbar\omega_{c}}{\kappa}=50$,
being slightly lower as in Ref.shalabney2015 . Note that the zero-point energy
contribution is neglected in Eq.22 to avoid artificial ground state
decayulusoy2020 .

Figure 7: Vibro-polaritonic infrared spectra subject to dissipation (colored)
with cavity decay rate, $\kappa$, for different initial states and light-
matter interaction regimes with non-dissipative reference (grey) and VPCI
energy (vertical line). Top row: Infrared-spectra for the vibrational strong
coupling (VSC) regime with $\eta=0.05$ for (a) molecular rovibrationally
excited initial state $\Psi^{(r)}_{0}$, (b) light-matter superposition state
$\Psi^{(r,z)}_{0}$ and (c) comparison of both dissipative spectra. Bottom row:
(d)-(f) analogous to (a)-(c) in top row for onset of vibrational ultrastrong
strong coupling (VUSC) regime at $\eta=0.1$.

In Fig.7, we compare infrared spectra subject to spontaneous emission
(dissipation) with non-dissipative results as presented in Fig.6. Most
prominently, cavity-loss effects manifest themselves in a significant
intensity reduction and peak broadening, especially for the $M_{1}$-peak
themselves . The latter again supports the purely photonic nature of the
$\ket{M_{1}}$-state. For the coupling scenario $\eta=0.05$, spectra resulting
from $\Psi^{(r)}_{0}$ and $\Psi^{(r,z)}_{0}$ differ only slightly in their
intensity as depicted in Figs.7(a)–(c), with slightly more intense
$L_{1}$-/$U_{1}$-transitions for $\Psi^{(r)}_{0}$. Further, the rovibro-
polaritonic progression induced by the VPCI is only weakly observable and
resembles several shoulders in the $U_{1}$-spectra region. For increased
light-matter interaction strength (cf. Figs.7(d)–(f)), the progression is
clearly visible for both initial state although it suffers from the peak
broadening effects.

##  IV Summary and Outlook

We studied a model system composed of a single rovibrating diatomic molecule
(carbon monoxide), which interacts with two energetically degenerate and
orthogonally polarized optical modes of a Fabry-Pérot-type cavity. The cavity
modes were tuned resonant to the fundamental vibrational transition of the CO
stretching mode and we discussed the vibrational strong coupling regime up to
the onset of the vibrational ultrastrong coupling regime. Based on an
energetically motivated adiabatic separation of rotational and vibro-
polaritonic degrees of freedom, we identified a three dimensional single-
excitation manifold of vibro-polaritonic surfaces, which is subject to two
distinct three-state vibro-polaritonic conical intersections in a two-
dimensional angular coordinate branching space. The lower and the upper vibro-
polaritonic states are of mixed light-matter hybrid character and the middle
polaritonic state is purely photonic in nature. In the VSC regime, the lower
and upper vibro-polaritonic surfaces exhibit a symmetric double-cone topology
with respect to a harmonic middle polaritonic surface at the triply degenerate
intersection point. For increasing light-matter interaction, the symmetric
splitting of upper and lower vibro-polaritonic surfaces becomes strongly
asymmetric, with the inverse cone of the upper vibro-polaritonic surface
dominating the VPCI topology and a mild “Mexican-hat”-type topology in the
lower vibro-polaritonic surface.

From a dynamical perspective, we studied the time-evolution of vibro-
polaritonic populations and reduced vibro-polaritonic rotational densities.
The population dynamics reveal, that the presence of rotational degrees of
freedom induces an efficient transfer channel between vibro-polaritonic states
in terms of the VPCI. Moreover, the rotational density shows a rich
interference pattern throughout the dynamics, due to interfering rotational
wave packets initially located on different vibro-polaritonic surfaces.

Further, we studied infrared spectra of the rovibrating light-matter hybrid
system with respect to different initial conditions, which mimicked classical
laser-field excitations of a (i) molecular rovibrational excited state and
(ii) a light-matter superposition state. We observe bright transitions to all
three vibro-polaritonic states, only when the laser field is allowed to
initially excite a light-matter superposition state. The non-trivial topology
of the VPCI manifest as a pronounced progression in the spectral region of the
upper vibro-polaritonic surface, which we find to correspond to rovibro-
polaritonic states containing contributions of molecular rotational,
vibrational and cavity degrees of freedom. Further, we considered finite
lifetimes of cavity mode excitations due to spontaneous emission effects,
which manifest in a significant reduction of the purely photonic middle
polariton peak and broadening of the rovibro-polaritonic progression.

Possible extensions of the present study refer to increasing the reduced
diabatic “vibro-polaritonic” basis to go beyond the (effective) single-
molecule picture and aim at the study of collective effects. In particular, it
would be instructive to investigate how the orientation dependent light-matter
interaction effects the nature of “dark states” in rovibro-polaritonic
problems. Further, we did not consider the nature of the vibro-polaritonic
ground state here, which renders central for purely rotational dynamics of
molecules in infrared Fabry-Pérot cavities. Moreover, we note that for
electronic strong coupling (ESC) problems, it had already been suggested in
Refs.triana2018 ; triana2019  that (single-mode) cavity induced CIs could be
formed, how one could detect them and what is the influence of molecular
rotation on themtriana2019 . This work is complemented by Ref.farag2021 ,
where a similar problem has been studied with respect to Berry Phase effects
in a polaritonic CI setting. Here, differences to “classical” light-induced
conical intersectionssindelka2011 , which emerge from a linearly polarized,
classical laser field with a single polarization direction, could be expected
in the fully quantized cQED setting of ESC scenarios featuring (at least) two
degenerate, orthogonally polarized cavity modes. Finally, the study of
rovibro-polaritonic effects on cavity-altered chemical reactions and on
polaritonic infrared spectra might be rewarding for a deeper understanding of
molecular cQED.

## Acknowledgements

We acknowledge fruitful discussions with Shreya Sinha, Foudhil Bouakline,
David Picconi and Tillmann Klamroth (all from Potsdam). The authors thank the
Deutsche Forschungsgemeinschaft (DFG) for financial support through project Sa
547/18. E.W. Fischer also acknowledges support by the International Max Planck
Research School for Elementary Processes in Physical Chemistry.

## Data Availability Statement

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

## Conflict of Interest

The authors have no conflicts to disclose.

## Appendix

###  IV.1 Numerical Details

We employ the MCTDH approachmeyer2012 ; meyer1990 ; manthe1992 ; beck2000 ;
meyer2003  as implemented in the Heidelberg MCTDH package, version
8.6.0mctdh2019 . The state-specific rotational wave packets
$\varphi_{k}(\theta,\phi,t)$ in Eq.(11) are here expanded as

$\varphi_{k}(\theta,\phi,t)=\sum^{n}_{i=1}A_{ki}(t)\,\varphi^{(k)}_{i}(\theta,\phi,t)\quad,$
(IV.11)

with time-dependent coefficients, $A_{ki}(t)$, and state- and time-dependent
rotational single-particle functions (SPFs),
$\varphi^{(k)}_{i}(\theta,\phi,t)$, respectively. For all calculations, we
chose $n=3$ SPFs for each diabatic state and represented the states-specific
rotational wave packets in a two-dimensional Legendre discrete variable
representation (PLeg) with $N_{\theta}=51$ and $N_{\phi}=37$ grid points.

Vibro-polaritonic (adiabatic) populations $P^{(l)}_{\mathrm{ad}}(t)$ are
obtained from Eq.(19) by taking into account the adiabatic representation of
the zero-order rotational wave packets, $\varphi_{k}(\theta,\phi,t)$, given as

$\displaystyle\varphi^{\mathrm{ad}}_{l}(\theta,\phi,t)$
$\displaystyle=\sum_{k}u_{lk}(\theta,\phi)\,\varphi_{k}(\theta,\phi,t)\quad,$
(IV.12)

where $u_{lk}(\theta,\phi)$ are elements of a unitary matrix
$\underline{\underline{U}}(\theta,\phi)$, which diagonalizes the potential in
Eq.(14). The vibro-polaritonic (adiabatic) rotational wave packets,
$\varphi^{\mathrm{ad}}_{l}(\theta,\phi,t)$, evolve on the surfaces as depicted
in Fig.2. Further, infrared spectra are calculated as (setting a possible
prefactor to 1)

$\sigma(\omega)=\displaystyle\int^{T}_{0}C(t)\,f_{W}(t)\,e^{\mathrm{i}\omega\,t}\mathrm{d}t\quad,$
(IV.13)

with autocorrelation function $C(t)$ and window function

$f_{W}(t)=\left(1-\dfrac{t}{T}\right)\cos\left(\dfrac{\pi}{T}\,t\right)+\dfrac{1}{\pi}\sin\left(\frac{\pi}{T}\,t\right)\quad.$
(IV.14)

For the Hermitian, non-dissipative case, we propagate a given initial wave
packet $\Psi(0)$ under the influence of the Pauli-Fierz Hamiltonian in Eq.(1),
to a time $t_{f}=3000\,\mathrm{fs}$ and the spectrum is obtained with
$T=2\,t_{f}=6000\,\mathrm{fs}$ since the autocorrelation function is evaluated
asbeck2000

$C(t)=\braket{\Psi^{\star}(t/2)}{\Psi(t/2)}\quad.$ (IV.15)

In the non-Hermitian scenario, $\hat{H}^{(\kappa)}_{C}$ as defined in Eq.(22)
is used instead for the cavity modes, and the autocorrelation function is
computed as

$C(t)=\braket{\Psi(0)}{\Psi(t)}\quad,$ (IV.16)

where we propagate up to $T=t_{f}=6000\,\mathrm{fs}$. Finally, for analysis,
intensities (stick spectra) are obtained as

$\begin{matrix}I^{(r)}=\sum_{i}|\braket{\Psi^{(r)}_{0}}{\phi_{i}}|^{2}\,\delta(\hbar\omega-\varepsilon_{i})\vspace{0.2cm}\\\
I^{(r,z)}=\sum_{i}|\braket{\Psi^{(r,z)}_{0}}{\phi_{i}}|^{2}\,\delta(\hbar\omega-\varepsilon_{i})\end{matrix}\quad,$
(IV.17)

with eigenvalues $\varepsilon_{i}$ and eigenstates $\phi_{i}$ of
$\underline{\underline{H}}$ calculated by means of a Lanczos algorithm as
implemented in the Heidelberg MCTDH package. All converged results have been
obtained based on a two-dimensional Legendre discrete variable representation
(PLeg) with $N_{\theta}=51$ and $N_{\phi}=37$ grid points as above, with
$9000$ iteration steps.

###  IV.2 3-State-VPCI Details

We expand trigonometric functions in $\underline{\underline{V}}(\theta,\phi)$
up to first order in $\theta$ and $\phi$ around the intersection point at
$(\theta,\phi)=(\frac{\pi}{2},0)$, as

$\cos\theta\approx\theta-\dfrac{\pi}{2}\quad,\hskip
14.22636pt\sin\theta\sin\phi\approx\phi\quad,$ (IV.21)

and

$f(\theta,\phi)\approx\phi^{2}+\left(\theta-\dfrac{\pi}{2}\right)^{2}=f^{(1)}(\theta,\phi)\quad.$
(IV.22)

The resulting approximate potential energy matrix is given by

$\underline{\underline{V}}^{(1)}(\theta,\phi)=\begin{pmatrix}G_{00}\,f^{(1)}(\theta,\phi)&G_{10}\,f^{(1)}(\theta,\phi)&-g_{00}\,\left(\theta-\dfrac{\pi}{2}\right)&g_{00}\,\phi\vspace{0.2cm}\\\
G_{10}\,f^{(1)}(\theta,\phi)&\hbar\omega_{10}+G_{11}\,f^{(1)}(\theta,\phi)&-g_{10}\,\left(\theta-\dfrac{\pi}{2}\right)&g_{10}\,\phi\vspace{0.2cm}\\\
-g_{00}\,\left(\theta-\dfrac{\pi}{2}\right)&-g_{10}\,\left(\theta-\dfrac{\pi}{2}\right)&\hbar\omega_{c}+G_{00}\,f^{(1)}(\theta,\phi)&0\vspace{0.2cm}\\\
g_{00}\,\phi&g_{10}\,\phi&0&\hbar\omega_{c}+G_{00}\,f^{(1)}(\theta,\phi)\end{pmatrix}\quad.$
(IV.23)

with $\hbar\omega_{10}=\hbar\omega_{c}$,
$G_{ij}=\frac{g^{2}}{\hbar\omega_{c}}\braket{d^{2}}_{ij}$ as well as
$g_{ij}=g\,d_{ij}$ where $i,j=0,1$. The linear approximation, holds for small
deviations from the intersection coordinate where $f^{(1)}(\theta,\phi)\approx
0$. For larger values of $\theta$ and $\phi$, the quadratic approximation then
takes into account the DSE-related term, $f^{(1)}(\theta,\phi)$. The
corresponding eigenvalues of $\underline{\underline{V}}^{(1)}(\theta,\phi)$
correspond to the different approximations of vibro-polaritonic surfaces with
cuts shown in Fig.2, bottom-row.

## References

  * (1) J. George, A. Shalabney, J. A. Hutchison, C. Genet, T. W. Ebbesen; J. Phys. Chem. Lett. 6, 1027, (2015).
  * (2) T. W. Ebbesen, Acc. Chem. Res. 49, 2403, (2016).
  * (3) A. Thomas, J. George, A. Shalabney, M. Dryzhakov, S. J. Varma, J. Moran, T. Chervy, X. Zhong, E. Devaux, C. Genet, J. A. Hutchison, T. W. Ebbesen, Angew. Chem. Int. Ed. 55, 11462, (2016).
  * (4) A. Thomas, L. Lethuillier-Karl, K. Nagarajan, R. M. A. Vergauwe, J. George, T. Chervy, A. Shalabney, E. Devaux, C. Genet, J. Moran, T. W. Ebbesen, Science 363, 615, (2019).
  * (5) F. Herrera, J. Owrutsky; J. Chem. Phys. 152, 100902, (2020).
  * (6) T. Chervy, A. Thomas, E. Akiki, R. M. A. Vergauwe, A. Shalabney, J. George, E. Devaux, J. A. Hutchison, C. Genet, T. W. Ebbesen; ACS Photonics 5, 217, (2018).
  * (7) B. Xiang, R. F.Ribeiro, A. D. Dunkelberger, J. Wang, Y. Li, B. S. Simpkins, J. C. Owrutsky, J. Yuen-Zhou, W. Xiong; PNAS 115, 4845, (2018).
  * (8) Z. Zhang, K. Wang, Z. Yi, M. Suhail Zubairy, M. O. Scully, S. Mukamel; J. Phys. Chem. Lett. 10, 4448, (2019).
  * (9) P. Saurabh, S. Mukamel; J. Chem. Phys. 144, 124115, (2016).
  * (10) R. F. Ribeiro, A. D. Dunkelberger, B. Xiang, W. Xiong, B. S. Simpkins, J. C. Owrutsky, J. Yuen-Zhou; J. Phys. Chem. Lett. 13, 3766, (2018).
  * (11) J. Flick, D.M. Welakuh, M. Ruggenthaler, H. Appel, A. Rubio; ACS Photonics 6, 2757, (2019).
  * (12) T. E. Li, J. E. Subotnik, A. Nitzan; PNAS 117, 18324, (2020).
  * (13) T. E. Li, A. Nitzan, J. E. Subotnik; J. Chem. Phys. 154, 094124, (2021).
  * (14) J. Bonini, J. Flick; arXiv:2108.11564, (2021).
  * (15) E. W. Fischer, P. Saalfrank, J. Chem. Phys. 154, 104311, (2021).
  * (16) J. Flick, H. Appel, M. Ruggenthaler, A. Rubio, J. Chem. Theory Comput. 13, 1616, (2017).
  * (17) J. Flick, M. Ruggenthaler, H. Appel, A. Rubio, PNAS 114, 3026, (2017).
  * (18) C. Schäfer, M. Ruggenthaler, A. Rubio, Phys. Rev. A 98, 043801, (2018).
  * (19) X. Li, A. Mandal, P. Huo; Nat. Commun 12, 1, (2021).
  * (20) J. Flick, P. Narang; Phys. Rev. Lett. 121, 113002, (2018).
  * (21) D. Sidler, M. Ruggenthaler, H. Appel, A. Rubio; J. Phys. Chem. Lett. 11, 7525, (2020).
  * (22) T. Szidarovszky, Péter Badankó, Gábor J Halász, Ágnes Vibók; J. Chem. Phys. 154, 064305, (2021).
  * (23) J. F. Triana, J. L. Sanz-Vicario; J. Chem. Phys. 154, 094120, (2021).
  * (24) D. P. Craig, T. Thirunmachandran; Molecular Quantum Electrodynamics, Dover Ed., Mineola, N.Y., (1998).
  * (25) V. Rokaj, D. M. Welakuh, M. Ruggenthaler, A. Rubio; J. Phys. B: At. Mol. Opt. Phys. 51, 034005, (2018)
  * (26) C. Schäfer, M. Ruggenthaler, V. Rokaj, A. Rubio; ACS Photonics 7, 975, (2020).
  * (27) A. F. Kockum, A. Miranowicz, S. De Liberato, S. Savasta, F. Nori, Nat. Rev. Phys. 1, 19, (2019).
  * (28) Gaussian 16, Revision C.01, M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian, Inc., Wallingford CT, 2016.
  * (29) D. H. Rank, A. G. St Pierre, T. A. Wiggins; J. Mol. Spec. 18, 418, (1965).
  * (30) D.T. Colbert, W.H. Miller; J. Chem. Phys. 96, 1982, (1992).
  * (31) N. Mina-Camilde , Carlos Manzanares I., J. F. Caballero; J. Chem. Educ. 73, 804, (1996).
  * (32) C. Hua-Jun, W. Jie, L. Hao, C. Xin-Lu. Chin. Phys. B, 24, 083102, (2015).
  * (33) G. Scuseria, M. D. Miller; J. Chem. Phys. 94, 6660, (1991).
  * (34) H. Köppel, W. Domcke, L.S. Cederbaum; Adv. Chem. Phys. 57, 59, (1984).
  * (35) S. Matsika, D. R. Yarkony; J. Am. Chem. Soc. 125, (2003), 10672.
  * (36) S. Matsika; J. Phys. Chem. A 109, 7538, (2005).
  * (37) J. D. Coe, T. J. Martinez; J. Am. Chem. Soc. 127, 4560, (2005).
  * (38) “Three-State Intersections” in Conical Intersections: Theory, Computation and Experiment by W. Domcke, D. R. Yarkony, H. Köppel; Advanced Series in Physical Chemistry. 17. World Scientific. (2011).
  * (39) We assume, $\mu^{(r)}=d^{(r)}_{10}\ket{1_{r},0_{z},0_{y}}\bra{0_{r},0_{z},0_{y}}\cos\theta$, with, $\mu^{(r)}\ket{0_{r},0_{z},0_{y}}Y^{0}_{0}=\Psi^{(r)}_{0}$, and, $\mu^{(r,z)}=d^{(r,z)}_{10}\left(\ket{1_{r},0_{z},0_{y}}\bra{0_{r},0_{z},0_{y}}+\ket{0_{r},1_{z},0_{y}}\bra{0_{r},0_{z},0_{y}}\right)\cos\theta$, such that, $\mu^{(r,z)}\ket{0_{r},0_{z},0_{y}}Y^{0}_{0}=\Psi^{(r,z)}_{0}$. For both, we set, $d^{(r)}_{10}=d^{(r,z)}_{10}=1$.
  * (40) H.-D. Meyer, WIREs Comput. Mol. Sci. 2, 351, (2012).
  * (41) G. A. Worth, M. H. Beck, A. Jäckle, and H.-D. Meyer. The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). O. Vendrell and H.-D. Meyer Version 8.5 (2013). Version 8.5 contains the ML-MCTDH algorithm. See http://mctdh.uni-hd.de. Used versions: 8.6.0 (2022).
  * (42) E.J. Heller, Acc. Chem. Res. 14, 368, (1981).
  * (43) I. S. Ulusoy, O. Vendrell; J. Chem. Phys. 153, 044108, (2020).
  * (44) S. Felicetti, J. Fregoni, T. Schnappinger, S. Reiter, R. de Vivie-Riedle, J. Feist; J. Phys. Chem. Lett. 11, 8810, (2020).
  * (45) P. Antoniou, F. Suchanek, J. F. Varner, J. J. Foley IV; J. Phys. Chem. Lett. 11, 9063, (2020).
  * (46) S. Shalabney, J. George, J. Hutchinson, G. Pupillo, C. Genet, T.W. Ebbesen; Nat. Commun. 6, 5981, (2015).
  * (47) J. F. Triana, D. Pelàez, José Luis Sanz-Vicario; J. Phys. Chem. A 122, 2266, (2018).
  * (48) J. F. Triana, J. L. Sanz-Vicario; Phys. Rev. Lett. 122, 063603, (2019).
  * (49) M. H. Faraq, A. Mandal, P.Huo; Phys. Chem. Chem. Phys. 23, 16868, (2021).
  * (50) M. Šindelka, N. Moiseyev, L. S. Cederbaum; J. Phys. B: At. Mol. Opt. Phys. 44, (2011).
  * (51) H.-D. Meyer, U. Manthe, L. S. Cederbaum, Chem. Phys. Lett. 165, 73, (1990).
  * (52) U. Manthe, H.-D. Meyer, L. S. Cederbaum, J. Chem. Phys. 97, 3199, (1992).
  * (53) M. H. Beck, A. Jäckle, G. A. Worth, H.-D. Meyer, Phys. Rep. 324, 1, (2000).
  * (54) H.-D. Meyer, G. A. Worth, Theor. Chem. Acc. 109, 251, (2003).


# Investigating Gated Recurrent Networks for Speech Synthesis

###### Abstract

Recently, recurrent neural networks (RNNs) as powerful sequence models have
re-emerged as a potential acoustic model for statistical parametric speech
synthesis (SPSS). The long short-term memory (LSTM) architecture is
particularly attractive because it addresses the vanishing gradient problem in
standard RNNs, making them easier to train. Although recent studies have
demonstrated that LSTMs can achieve significantly better performance on SPSS
than deep feed-forward neural networks, little is known about why. Here we
attempt to answer two questions: a) why do LSTMs work well as a sequence model
for SPSS; b) which component (e.g., input gate, output gate, forget gate) is
most important. We present a visual analysis alongside a series of
experiments, resulting in a proposal for a simplified architecture. The
simplified architecture has significantly fewer parameters than an LSTM, thus
reducing generation complexity considerably without degrading quality.

Index Terms—  Speech synthesis, acoustic modelling, recurrent network network,
gated recurrent network, long short-term memory

##  1 Introduction

Statistical parametric speech synthesis (SPSS) has quite steadily advanced in
naturalness in the past decade, as witnessed by the series of Blizzard
Challenges [1]. However, the quality of synthetic speech produced by SPSS is
still far below that of the natural human speech, and cannot compete with the
best unit selection systems, which concatenate waveforms [2]. As suggested in
[3], acoustic modelling, which captures the complex relationship between
linguistic and acoustic representations, is a key limiting factor and is the
focus of this work.

###  1.1 Relation to prior work

Neural networks have re-emerged as a potential powerful acoustic model for
SPSS. In [4, 5, 6, 7, 8], feed-forward neural networks are employed to map a
linguistic representation derived from input text directly to acoustic
features. In [9], a deep belief network (DBN) was used to model the
relationship between linguistic and acoustic representations jointly. In [10]
and [11], mixture density networks (MDNs) and real-valued neural
autoregressive density estimators (RNADEs) were proposed, respectively, to
predict acoustic feature distributions given input linguistic features. These
various implementations can be viewed as a replacement of the decision tree in
HMM-based speech synthesis; they map linguistic features to acoustic features
frame by frame through multiple hidden layers. However, the temporal sequence
nature of speech is not explicitly modelled in the network architectures.

To include temporal constraints, we proposed to include contextual information
by stacking low-dimensional bottleneck features from multiple consecutive
frames [12]. Still in the DNN framework, minimum trajectory error training
[13] or sequence error training criterion [14] have been proposed to minimise
the utterance-level trajectory error rather than the frame-by-frame error. On
the other hand, recurrent neural networks (RNNs) directly and elegantly
include temporal information in the network architecture, making them
attractive for modelling speech parameter trajectories. In [15], a standard
RNN was employed to predict prosodic information for speech synthesis. In
[16], two variants on standard RNNs, the Elman RNN and clockwork RNN, were
investigated for speech synthesis.

The most widely used recurrent network in speech processing applications is
the long short-term memory (LSTM) architecture. Because the LSTM addresses the
vanishing gradient problem of the standard RNN, it is easier to train. In
[17], an LSTM was employed to model the F0 contour. In [18], a bidirectional
LSTM was employed to map a sequence of linguistic features to the
corresponding sequence of acoustic features. In [19], an LSTM with a recurrent
output layer was proposed to perform sequence mapping from linguistic to
acoustic representations. These studies all formulate SPSS as sequence-to-
sequence mapping and all demonstrate the effectiveness of LSTMs. However, LSTM
architecture seems rather ad-hoc and it is not obvious what its various
components are actually contributing to performance.

This raises at least two questions that have not been answered in previous
studies: a) _how exactly does the LSTM architecture model a speech parameter
sequence_ ; b) _which components of the LSTM architecture are important, and
which could be discarded_. Answers to these questions may suggest better and
perhaps simpler recurrent network architectures.

###  1.2 The novelty of this work

We attempt to reach a better understanding of the “black-box” LSTM
architecture and our findings lead us to propose a simplified architecture for
speech modelling.

First, we give an analysis of the forget gate and memory cell in the LSTM
architecture. Specifically, we visualise the activation of the forget gate to
understand when the forget gate resets the memory cell state, and how the
forget gate relates to speech structure. We analyse how the cell state
correlates with the trajectory to be predicted. These visualisations enable us
to understand how LSTMs model the temporal structure in speech synthesis. _To
the best of our knowledge, this is the first attempt to visually analyse the
LSTM architecture in predicting a speech parameter sequence_.

Second, we analyse the importance of each LSTM component for speech synthesis
and propose a simplified architecture. The analysis is done empirically with
several variants of the LSTM. Each removes a different component of the
vanilla LSTM. The analysis was inspired by the studies in [20, 21], and we
focus on the speech synthesis application. Based on this analysis, we present
a simplified architecture, which only has the forget gate. The simplified
architecture has significantly fewer parameters than the vanilla LSTM, and so
reduces the computational cost of generation considerably without degrading
the quality of the synthesised speech.

##  2 Long Short-Term Memory

Standard RNNs are hard to train due to the well-known vanishing or exploding
gradient problems [22, 23]. To address the vanishing gradient problem, the
LSTM architecture was proposed, the basic idea of which was presented in [24].
The most commonly used architecture was described in [25], and is formulated
as,

$\displaystyle\mathbf{i}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{i}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{i}}\mathbf{h}_{t-1}+\mathbf{p}^{\mathrm{i}}\odot\mathbf{c}_{t-1}+\mathbf{b}^{\mathrm{i}})$
$\displaystyle\mathbf{f}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{f}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{f}}\mathbf{h}_{t-1}+\mathbf{p}^{\mathrm{f}}\odot\mathbf{c}_{t-1}+\mathbf{b}^{\mathrm{f}})$
$\displaystyle\mathbf{c}_{t}$
$\displaystyle=\mathbf{f}_{t}\odot\mathbf{c}_{t-1}+\mathbf{i}_{t}\odot
g(\mathbf{W}^{\mathrm{c}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{c}}\mathbf{h}_{t-1}+\mathbf{b}^{\mathrm{c}})$
$\displaystyle\mathbf{o}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{o}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{o}}\mathbf{h}_{t-1}+\mathbf{p}^{\mathrm{o}}\odot\mathbf{c}_{t}+\mathbf{b}^{\mathrm{o}})$
$\displaystyle\mathbf{h}_{t}$ $\displaystyle=\mathbf{o}_{t}\odot
g(\mathbf{c}_{t})$

In these formulations, $\mathbf{i}_{t}$, $\mathbf{f}_{t}$, $\mathbf{c}_{t}$,
$\mathbf{o}_{t}$, and $\mathbf{h}_{t}$ are the input gate, forget gate, cell
state, output gate and block output at time instance $t$, respectively;
$\delta(\cdot)$ and $g(\cdot)$ are the sigmoid and tangent activation
functions, respectively; $\mathbf{x}_{t}$ is the input at time $t$;
$\mathbf{W}^{\mathrm{*}}$, and $\mathbf{R}^{\mathrm{*}}$ are the weight
matrices applied on input and recurrent hidden units, respectively;
$\mathbf{p}^{*}$ and $\mathbf{b}^{*}$ are the peep-hole connections and
biases, respectively; and $\odot$ means element-wise product. We will call
this the _vanilla LSTM_.

The central idea of the LSTM is the so-called memory cell $\mathbf{c}$ which
maintains its state over time, and the gating units which are used to regulate
the information flow into and out of the memory cell [20]. More specifically,
the input gate can allow the input signal to adjust the cell state or prevent
that (e.g., setting the input gate to zero); the output gate can allow the
cell state to affect other neurons or block that; and the forget gate enables
the cell to remember or forget its previous state. However, as discussed in
[20, 21], the architecture might not be optimal for all the tasks, and the
relative importance of each component is not at all clear.

##  3 Gated Recurrent Neural Networks

In this section, we present several variants of the LSTM and propose a
simplified version that only has the forget gate; it therefore has
significantly fewer parameters and lower computational cost. As these variants
all share with the LSTM the concept of a memory cell with gates, we will call
them _gated recurrent neural networks_.

###  3.1 Four variants on the LSTM

To assess the importance of each component, we start with four variants of the
LSTM architecture. Each removes one component from the LSTM architecture, so
we can understand how much each component contributes to performance. The
differences with the vanilla LSTM are:

  * •

No Peep-holes (NPH): Set $\mathbf{p}^{\mathrm{i}}$, $\mathbf{p}^{\mathrm{f}}$,
$\mathbf{p}^{\mathrm{o}}$ to zero

  * •

No input gate (NIG): $\mathbf{i}_{t}=1$

  * •

No forget gate (NFG): $\mathbf{f}_{t}=1$

  * •

No output gate (NOG): $\mathbf{o}_{t}=1$

In the NFG variant, the past cell state will still contribute to the current
cell state but without any controlling or scaling by the forget gate. Note
that, when removing the input, forget or output gates, the number of
parameters is reduced.

###  3.2 Gated Recurrent Unit (GRU)

As an alternative to the LSTM, the Gated Recurrent Unit (GRU) architecture was
proposed in [26]. In [27], the GRU was found to achieve better performance
than the LSTM on some tasks. The GRU is formulated as:

$\displaystyle\mathbf{r}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{r}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{r}}\mathbf{h}_{t-1}+\mathbf{b}^{\mathrm{r}})$
$\displaystyle\mathbf{z}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{z}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{z}}\mathbf{h}_{t-1}+\mathbf{b}^{\mathrm{z}})$
$\displaystyle\tilde{\mathbf{h}}_{t}$
$\displaystyle=g(\mathbf{W}^{\mathrm{h}}\mathbf{x}_{t}+\mathbf{r}_{t}\odot(\mathbf{R}^{\mathrm{h}}\mathbf{h}_{t-1})+\mathbf{b}^{\mathrm{h}})$
$\displaystyle\mathbf{h}_{t}$
$\displaystyle=\mathbf{z}_{t}\odot\mathbf{h}_{t-1}+(1-\mathbf{z}_{t})\odot\tilde{\mathbf{h}}_{t}$

From these formulae, we can observed that the GRU architecture is similar to
LSTM but without a separate memory cell. The GRU does not use peep-hole
connections and output activation functions, and combines the input and forget
gates into an update gate $\mathbf{z}_{t}$ to balance between previous
activation $\mathbf{h}_{t-1}$ and the candidate activation
$\tilde{\mathbf{h}}_{t}$. The reset gate $\mathbf{r}_{t}$ allows it to forget
the previous state.

###  3.3 Simplified LSTM (S-LSTM)

As we will see in the experiments reported in the next section, the input
gate, output gate and peep-hole connections can be removed without degrading
speech synthesis performance significantly. Hence, we can propose an even
simpler variant, that removes output gates and peep-hole connections, and
replaces the input gate by the forget gate in the form of $1-\mathbf{f}_{t}$.
In this way, only the forget gate is retained. This simplest variant can be
written as:

$\displaystyle\mathbf{f}_{t}$
$\displaystyle=\delta(\mathbf{W}^{\mathrm{f}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{f}}\mathbf{h}_{t-1}+\mathbf{b}^{\mathrm{f}})$
$\displaystyle\mathbf{c}_{t}$
$\displaystyle=\mathbf{f}_{t}\odot\mathbf{c}_{t-1}+(1-\mathbf{f}_{t})\odot
g(\mathbf{W}^{\mathrm{c}}\mathbf{x}_{t}+\mathbf{R}^{\mathrm{c}}\mathbf{h}_{t-1}+\mathbf{b}^{\mathrm{c}})$
$\displaystyle\mathbf{h}_{t}$ $\displaystyle=g(\mathbf{c}_{t})$

The simplified architecture is similar to the GRU, except that it uses a
memory cell state. The cell state is controlled by the forget gate only, which
trades off between past cell state and current block input. When the
activation of forget gate is small, the cell state will mainly depend on the
block input, otherwise it will mainly copy the past cell state.

Fig. 1: Averaged activations of all 256 forget gates as a function of time (in
frames) are shown by the solid red line, with phoneme boundaries marked as
dashed vertical blue lines.

##  4 Experiments

###  4.1 Experimental setup

A corpus from a British male speaker was employed in our experiments, divided
into three subsets: training, development and testing (2400, 70 and 72
utterances). The sampling rate was 48 kHz, and we used the STRAIGHT vocoder
[28] to extract 60-dimensional Mel-Cepstral Coefficients (MCCs), 25 band
aperiodicities (BAPs), and fundamental frequency ($F_{0}$) on log-scale, all
at 5-ms frame step. All systems used the same acoustic features. $F_{0}$ was
linearly interpolated before modelling and a binary voiced/unvoiced feature
was used to record voicing information. Dynamic features for MCCs, BAPs and
$F_{0}$ were also computed. The acoustic features were mean-variance
normalised before modelling, and the mean and variance was restored at the
generation time. At generation time, maximum likelihood parameter generation
algorithm [29] was applied to smooth parameter trajectories.

All systems used the same input linguistic features comprising 601 features.
592 of these are binary features derived from linguistic context, such as
quin-phone identities, part-of-speech, positional information of phoneme,
syllable, word and phrase, and the number of syllables, words and phrases,
etc. The remaining 9 numerical features capture frame position information,
e.g., frame position in HMM state and phoneme. Linguistic features were
normalised to [0.01 0.99] before modelling.

In all RNNs, we employed a three-layer feed-forward neural network at the
bottom. On top of the feed-forward layers, we used the gated recurrent neural
networks. The bottom feed-forward layers were intended to act as feature
extraction layers, with 512 hidden units using tangent activation function in
each layer. All RNN implementations used 256 units (e.g., LSTM blocks) in the
recurrent layer. Hyperparameters for each system were optimised on the
development set. We fixed the momentum, and only tuned learning rates.

Fig. 2: Comparison between the first MCC trajectory and the $120^{\mathrm{th}}$ cell state (vertically offset for clarity). Table 1: Subjective preference scores (%). $p<0.01$ indicates significant difference between the two systems. | LSTM | NIG | NOG | NFG | NPH | GRU | S-LSTM | Neutral | p-value
---|---|---|---|---|---|---|---|---|---
1 | 46.7 | 42.3 | - | - | - | - | - | 11.0 | 0.1582
2 | 40.7 | - | 47.8 | - | - | - | - | 11.5 | 0.0649
3 | 74.3 | - | - | 20.3 | - | - | - | 5.4 | $<10^{-8}$
4 | 44.0 | - | - | - | 45.3 | - | - | 10.7 | 0.7820
5 | 47.8 | - | - | - | - | 42.7 | - | 9.5 | 0.1086
6 | 46.5 | - | - | - | - | - | 45.3 | 8.2 | 0.7948
7 | - | 46.0 | - | - | - | - | 46.8 | 7.2 | 0.8305
8 | - | - | 48.0 | - | - | - | 44.3 | 7.7 | 0.3828
9 | - | - | - | - | 51.7 | - | 40.5 | 7.8 | 0.0018
10 | - | - | - | - | - | 48.0 | 45.0 | 7.0 | 0.4703

###  4.2 Analysis of LSTM

We first visualised the forget gate and cell state, which are thought to be
the two most important components in modelling long-term temporal structure.
The averaged activations (over the 256 units) of the forget gate as a function
of the frame index is presented in Fig. 1. The red solid line is the forget
gates averaged activations; blue dashed lines show phoneme boundaries. It is
clear that the peaks of the forget gate activation trajectory have a strong
correspondence with the phoneme boundaries; within a phoneme, the contribution
of past cell state decays linearly. The forget gate is capturing some
important temporal structure of speech; this is not surprising, since the
phoneme boundaries are explicitly represented in the input linguistic
features.

Table 2: Objective measures. MCD: Mel-Cepstral Distortion. BAP: distortion of band aperiodicities. F0 RMSE is calculated on a linear scale. V/UV: voiced/unvoiced error. Note that the number of parameters listed is for the recurrent layer only and does not include the bottom three feed-forward layers, which are the same size across all systems. The generation time is to generate all the 142 utterances in both development and testing sets. | MCD | BAP | F0 RMSE | V/UV | # parameters | generation
---|---|---|---|---|---|---
| (dB) | (dB) | (Hz) | (%) |  | time (s)
LSTM | 4.14 | 1.95 | 8.96 | 4.15 | 788, 224 | 214
NIG | 4.18 | 1.95 | 9.10 | 4.15 | 591, 104 | 180
NOG | 4.14 | 1.94 | 8.84 | 4.29 | 591, 104 | 167
NFG | 4.68 | 1.99 | 9.69 | 4.41 | 591, 104 | 174
NPH | 4.14 | 1.95 | 9.02 | 4.20 | 787, 456 | 180
GRU | 4.17 | 1.95 | 9.00 | 4.22 | 590, 592 | 159
S-LSTM | 4.19 | 1.95 | 8.87 | 4.14 | 393, 728 | 154

The memory cell should maintains its state over time [20] and so could store
the trend of the trajectory to be predicted. To analyse the relationship
between the cell states and the MCC trajectories, we computed the correlation
between the cell states and the first MCC trajectory, and found that the
$120^{\mathrm{th}}$ cell state has the highest correlation with the first MCC
trajectory. The correlation is as high as 0.9. A comparison between them is
presented in Fig. 2, which shows that the cell state tracks the shape of the
MCC trajectory.

###  4.3 Objective results

Even though objective measures might not always correlate with human
perception, they offer a way to tune the systems and roughly predict model
performance. The objective results are in Table 2. Compared to LSTM, NIG, NOG
and NPH all achieve similar objective distortion, with considerably fewer
parameters and lower generation time: the input gate, output gate and peep-
hole connections are not necessary. The NFG system increases distortion
considerably: the forget gate is important. This finding is consistent with
[20].

The GRU system achieves similar performance to the LSTM system: even though it
has even fewer parameters, it performs as well as NIG, NOG or NPH. This is
also consistent with studies on other tasks [27, 21]. Although S-LSTM slightly
increases MCD distortion from 4.14 dB to 4.19 dB compared to LSTM, it achieves
similar performance on the other measures. The S-LSTM has about half the
number of parameters in its recurrent layer compared to the LSTM, and reduces
generation time from 214 seconds to 154 seconds. The generation time is the
total time to generate all the 142 utterances in both development and testing
sets.

In summary, the S-LSTM has the smallest number of parameters and achieves the
fastest generation, whilst achieving similar objective results to the LSTM and
GRU architectures.

###  4.4 Subjective results

Subjective preference tests were conducted using 30 paid native English
speakers. Each listener was asked to listen 20 pairs of synthesised
utterances. The sentence was the same in both items within a pair, and was
randomly selected from the 72 test sentences111Samples are available at:
http://homepages.inf.ed.ac.uk/zwu2/demo/icassp16/lstm.html. For each pair, the
listener was asked to decide which one sounded more natural; a “neutral”
option was allowed if the listener had no preference.

Preference results are in Table 1. Comparing against the LSTM system, all the
systems except NFG show no significant difference in preference.

The NFG system achieves only a 20.3% preference score when paired against the
LSTM which is preferred 74.3% of the time. As with the objective results in
Table 2, we conclude that the forget gate is the only critical component in
the LSTM architecture; the input gate, output gate and peep-hole connections
can be omitted.

We also compares the proposed S-LSTM system against with all other systems
(except NFG, since it is worse than LSTM). Consistent with the objective
results, the subjective results also demonstrate that S-LSTM is as good as any
other systems.

##  5 Conclusions

We have analysed the forget gate and cell state of the LSTM architecture, and
examined the performance of several variants of LSTM. We conclude that:

  * •

The forget gate can learn the temporal structure of speech; its activations
have a high correspondence with phone boundaries.

  * •

The memory cell maintains a state over time, which matched the shape of the
trajectory to be predicted.

  * •

For this task, the forget gate is the only critical component of the LSTM;
other components can be omitted with no reduction in naturalness.

From these results, we propose a simplified LSTM architecture that only uses
the critical forget gate. The simplified LSTM has significantly fewer
parameters than the vanilla LSTM, but achieves similar performance in both
objective and subjective evaluations.

Acknowledgements: This research was supported by EPSRC Programme Grant
EP/I031022/1, Natural Speech Technology (NST). The NST research data
collection may be accessed at
http://datashare.is.ed.ac.uk/handle/10283/786.

## References

  * [1] Simon King,  “Measuring a decade of progress in text-to-speech,”  Loquens, vol. 1, no. 1, 2014.
  * [2] Andrew J Hunt and Alan W Black,  “Unit selection in a concatenative speech synthesis system using a large speech database,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 1996.
  * [3] Heiga Zen, Keiichi Tokuda, and Alan W Black,  “Statistical parametric speech synthesis,”  Speech Communication, vol. 51, no. 11, pp. 1039–1064, 2009.
  * [4] Heiga Zen, Andrew Senior, and Mike Schuster,  “Statistical parametric speech synthesis using deep neural networks,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2013.
  * [5] Heng Lu, Simon King, and Oliver Watts,  “Combining a vector space representation of linguistic context with a deep neural network for text-to-speech synthesis,”  Proc. the 8th ISCA Speech Synthesis Workshop (SSW), 2013.
  * [6] Yao Qian, Yuchen Fan, Wenping Hu, and Frank K Soong,  “On the training aspects of deep neural network (DNN) for parametric TTS synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2014.
  * [7] Qiong Hu, Zhizheng Wu, Korin Richmond, Junichi Yamagishi, Yannis Stylianou, and Ranniery Maia,  “Fusion of multiple parameterisations for DNN-based sinusoidal speech synthesis with multi-task learning,”  in Proc. Interspeech, 2015.
  * [8] Cassia Valentini-Botinhao, Zhizheng Wu, and Simon King,  “Towards minimum perceptual error training for DNN-based speech synthesis,”  in Proc. Interspeech, 2015.
  * [9] Shiyin Kang, Xiaojun Qian, and Helen Meng,  “Multi-distribution deep belief network for speech synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2013.
  * [10] Heiga Zen and Andrew Senior,  “Deep mixture density networks for acoustic modeling in statistical parametric speech synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2014.
  * [11] Benigno Uria, Iain Murray, Steve Renals, and Cassia Valentini,  “Modelling acoustic feature dependencies with artificial neural networks: Trajectory-rnade,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2015.
  * [12] Zhizheng Wu, Cassia Valentini-Botinhao, Oliver Watts, and Simon King,  “Deep neural networks employing multi-task learning and stacked bottleneck features for speech synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2015.
  * [13] Zhizheng Wu and Simon King,  “Minimum trajectory error training for deep neural networks, combined with stacked bottleneck features,”  in Proc. Interspeech, 2015.
  * [14] Yuchen Fan, Yao Qian, Frank K. Soong, and Lei He,  “Sequence generation error (SGE) minimization based deep neural networks training for text-to-speech synthesis,”  in Proc. Interspeech, 2015.
  * [15] Sin-Horng Chen, Shaw-Hwa Hwang, and Yih-Ru Wang,  “An RNN-based prosodic information synthesizer for Mandarin text-to-speech,”  IEEE Transactions on Speech and Audio Processing, vol. 6, no. 3, pp. 226–239, 1998.
  * [16] Sivanand Achanta, Tejas Godambe, and Suryakanth V. Gangashetty,  “An investigation of recurrent neural network architectures for statistical parametric speech synthesis,”  in Proc. Interspeech, 2015.
  * [17] Raul Fernandez, Asaf Rendel, Bhuvana Ramabhadran, and Ron Hoory,  “Prosody contour prediction with long short-term memory, bi-directional, deep recurrent neural networks,”  in Proc. Interspeech, 2014.
  * [18] Yuchen Fan, Yao Qian, Fenglong Xie, and Frank K. Soong,  “TTS synthesis with bidirectional LSTM based recurrent neural networks,”  in Proc. Interspeech, 2014.
  * [19] Heiga Zen and Hasim Sak,  “Unidirectional long short-term memory recurrent neural network with recurrent output layer for low-latency speech synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2015.
  * [20] Klaus Greff, Rupesh Kumar Srivastava, Jan Koutník, Bas R Steunebrink, and Jürgen Schmidhuber,  “Lstm: A search space odyssey,”  arXiv preprint arXiv:1503.04069, 2015.
  * [21] Rafal Jozefowicz, Wojciech Zaremba, and Ilya Sutskever,  “An empirical exploration of recurrent network architectures,”  in Proc. IEEE Int. Conf. on Machine Learning (ICML), 2015.
  * [22] Yoshua Bengio, Patrice Simard, and Paolo Frasconi,  “Learning long-term dependencies with gradient descent is difficult,”  IEEE Transactions on Neural Networks, vol. 5, no. 2, pp. 157–166, 1994.
  * [23] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio,  “On the difficulty of training recurrent neural networks,”  arXiv preprint arXiv:1211.5063, 2012.
  * [24] Sepp Hochreiter and Jürgen Schmidhuber,  “Long short-term memory,”  Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [25] Alex Graves and Jürgen Schmidhuber,  “Framewise phoneme classification with bidirectional lstm and other neural network architectures,”  Neural Networks, vol. 18, no. 5, pp. 602–610, 2005.
  * [26] Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio,  “Learning phrase representations using rnn encoder-decoder for statistical machine translation,”  arXiv preprint arXiv:1406.1078, 2014.
  * [27] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio,  “Empirical evaluation of gated recurrent neural networks on sequence modeling,”  arXiv preprint arXiv:1412.3555, 2014.
  * [28] Hideki Kawahara, Ikuyo Masuda-Katsuse, and Alain de Cheveigné,  “Restructuring speech representations using a pitch-adaptive time–frequency smoothing and an instantaneous-frequency-based F0 extraction: Possible role of a repetitive structure in sounds,”  Speech communication, vol. 27, no. 3, pp. 187–207, 1999.
  * [29] Keiichi Tokuda, Takayoshi Yoshimura, Takashi Masuko, Takao Kobayashi, and Tadashi Kitamura,  “Speech parameter generation algorithms for HMM-based speech synthesis,”  in Proc. IEEE Int. Conf. on Acoustics, Speech, and Signal Processing (ICASSP), 2000.