# On the unusual Stokes shift in the smallest PPE dendrimer building block:
Role of the vibronic symmetry on the band origin?

Joachim Galiana Benjamin Lasorne ICGM, Univ Montpellier, CNRS, ENSCM,
Montpellier, France<EMAIL_ADDRESS>
###### Abstract

1,3-bis(phenylethynyl)benzene is the primary chromophore of light-harvesting
polyphenylene ethynylene (PPE) dendrimers. It is experimentally known to share
the same absorption spectrum as its pair of diphenylacetylene (aka. tolane)
meta-substituted branches, yet exhibits an unusual Stokes shift of about 2000
cm-1 with respect to its band origin (corresponding to the loss of one
vibrational quantum within the antisymmetric acetylenic stretching) in its
emission spectrum. We suggest in the present work the unusual but plausible
involvement of molecular symmetry selection rules in a situation where the
Born-Oppenheimer approximation is far to be valid. Our hypothesis is comforted
with quantum dynamics (MCTDH) simulations of absorption and emission UV-
visible spectra based on quantum chemistry (TD-DFT) data and a diabatic
vibronic coupling Hamiltonian model.

Keywords: light-harvesting molecules, PPE dendrimers, UV-vis spectroscopy,
Stokes shift, nonadiabatic quantum dynamics, vibronic coupling, conical
intersections.

The following article has been submitted to The Journal of Chemical Physics.

††preprint: AIP/123-QED

##  I Introduction

Light-harvesting dendrimers have drawn much interest in many fields of
applications for the past few decades, such as in the design of new
photovoltaic technologies or photosynthetic processes.balzani_harvesting_1995
; gilat_light_1999 ; adronov_light-harvesting_2000 ; balzani_light-
harvesting_2003 Phenylacetylene dendrimers (covalent tree-like macromolecules)
are among the most promising $\pi$-conjugated organic systems for harvesting
and amplifying the energy of sunlight, based on simple polyphenylene
ethynylenes (PPEs) building blocks connected via meta-substitution at
phenylene nodes. Such building blocks are also of interest on their own, as
some of their spectroscopic properties are still challenging to understand and
interpret.

One of the best known phenylacetylene dendrimer is the nano-starxu_design_1994
, for which the building blocks are PPEs with different lengths and
substitutions schemes (para- or meta-substitution). The connectivity of the
building blocks within the nano-star is such that the structure is two-
dimensional, with a sequential threefold branched graph (two donors on each
acceptor). Experimental and theoretical studies so far have suggested two main
results. On the one hand, the absorption spectra of the building blocks of PPE
dendrimers are similar and additive, with mainly local electronic excitations
(LE).kopelman_spectroscopic_1997 ; swallen_dendrimer_1999 ;
swallen_correlated_2000 ; rana_steady-state_2001 ; wu_exciton_2006 ;
palma_electronic_2010 On the other hand, the building blocks are arranged such
that a local excitation at the periphery of the nano-star is transferred to
its center thanks to an intramolecular electronic excitation energy gradient
from the external chromophores to the energy-releasing
core.devadoss_energy_1996 ; shortreed_directed_1997 ; kleiman_ultrafast_2001 ;
melinger_optical_2002 ; fernandez-alberti_nonadiabatic_2009 ; fernandez-
alberti_shishiodoshi_2012 ; fernandez-alberti_non-adiabatic_2016 ;
nelson_electronic_2017

Among these building blocks, the 1,3-bis(phenylethynyl)benzene (fig. 1,
bottom, called m22 in the following) shows interesting steady-state
spectroscopic properties. Alike other building blocks based on
diphenylacetylene (fig. 1, top, called p2 in the following; also known as
tolane), it has an absorption spectrum similar in its vibronic structure to
the absorption spectrum of diphenylacetylene, only with a different intensity.

Figure 1: Diphenylacetylene p2 and 1,3-bis(phenylethynyl)benzene m22
molecules. A scheme of the displacement associated to normal modes of
vibration 81, 87 and 88 is given for m22. The Mulliken axis convention used
here is such that in the C2v symmetry point group, the normal mode of
vibration 87 is of in-plane $y$-type symmetry and belongs to B2.

This peculiarity was assigned to the additivity of local single excitations on
p2 branches in the building blocks, and further rationalized with a pseudo-
fragmentation scheme for m22.thompson_variable_2004 ; ho_diabatic_2019 In
addition, a minimum energy conical intersection (MECI) was characterized
between the first two singlet electronic excited states.ho_diabatic_2019 The
nature of the coupled electronic singlet excited states was further studied,
with two quasi-diabatic representations: one with a locally excited state on
each of the p2 branch (right or left), the other with delocalized excited
states with a well-defined symmetry: B2 (right minus left) or A1 (right plus
left).

However, whereas the absorption spectra are almost identical for m22 and p2,
the emission spectrum of m22 is completely different from the emission
spectrum of p2 experimentally. An unusual Stokes shift
$\Delta\bar{\nu}\simeq$2000\text{\,}{\mathrm{cm}}^{-1}$$ has been reported for
the vibronic spectra of m22 separately by different
authors.thompson_variable_2004 ; gaab_meta-conjugation_2003 ;
chu_vibronic_2004 Let us stress out here that we are not considering the
typical textbook Stokes shift between the maxima of bell-shaped absorption and
emission spectra, essentially due to the influence of a significant
displacement between ground- and excited-state equilibrium geometries on
vertical absorption and emission transition energies (hence on maximal Franck-
Condon factors $0^{\prime}\rightarrow v^{\prime\prime}>0$ _vs._
$0^{\prime\prime}\rightarrow v^{\prime}>0$), but rather focusing on the
relative spectral shift of the two band origins ($0^{\prime}\leftrightarrow
0^{\prime\prime}$), which occur to coincide indeed and be the most intense for
p2 but are strangely shifted in the triangular-shaped emission of m22 and not
in its triangular-shaped absorption spectrum at low temperature.

A first hypothesis was suggested by some authorschu_vibronic_2004 of a missing
band origin associated to the $0-0$ transition in both the absorption and
emission spectra. Such a proposition does not seem satisfactory because the
$0-0$ transition is not missing in the absorption spectrum of m22, which is
similar to the one of p2 and for well-justified reasons. It thus seems more
likely to investigate some particular effect influencing the emission spectrum
of m22 compared to that of p2. This is precisely the objective of the present
work, where we propose some potential theoretical explanations based on
symmetry, yet cannot provide any definitive answer until further time-resolved
experiments are carried out.

Gaab, Thompson and co-workers produced an extensive experimental study for
m22-like molecules, and proposed a variable excitonic coupling model for the
reproduction and interpretation of the Stokes shift in H-terminated meta-
substituted phenylacetylene instead of phenyl-terminated
onesthompson_variable_2004 ; gaab_meta-conjugation_2003 . Their model
qualitatively agrees with _ab initio_ calculations for the ordering of
electronic states of interest in the absorption and emission phenomena for
meta-substituted PPEs. In their work, the theoretical Stokes shift is
reproduced by only looking at vertical transition energies for absorption and
emission. Once again, such a picture implies some type of usual Stokes shift
due to some geometrical displacement between the absorbing and emissive
species and/or incomplete vibrational relaxation. Further along this line, we
could invoke some additional effect due to the relaxed or unrelaxed state of
the solvent. Yet, there seems to some more fundamental effect due to the
nature of the absorbing and emissive species while considering steady-state
spectra, expected to exhibit the _a priori_ common $0^{\prime}\leftrightarrow
0^{\prime\prime}$ band origin connecting via light-matter resonance the long-
lived (stationary) eigenstates of the system.

As such, we propose hereby a minimal (three-dimensional; 3D) vibronic coupling
model, in the fashion of Cederbaum, Köppel, Cattarius, and co-
workerscederbaum_strong_1977 ; koppel_multimode_1984 ; cattarius_all_2001 for
the construction of a model vibronic Hamiltonian, able to reproduce the
vibrational structure of the electronic absorption and emission spectra for
m22 via nonadiabatic quantum dynamics techniques beyond the Born-Oppenheimer
approximation. This has never been tried on the present system to the extent
of our knowledge; we shall use it to discriminate Stokes and non-Stokes shift
contributions to the emission spectrum according to light polarization. Let us
also stress that such a strategy based on the spectral response of quantum
dynamics has been largely tested on photoabsorption or
photoelectron/photoionization spectra, but rarely on photoemission spectra,
which somewhat remain uncharted territory in terms of spontaneous _vs._
stimulated emission contributions.

With the vibronic coupling models used in this work, the coupling explicitly
depends on the geometry of the molecule, through variations of normal mode
coordinates, following the primary idea of a variable excitonic coupling
proposed by Gaab and Thompson. If one may say, note that our interpretation is
more the one of a chemist point of view (vibronic coupling) than of a
physicist point of view (excitonic coupling). Indeed, our vibronic coupling
models allow for the interpretation of the absorption and emission spectra in
terms of vibronic transitions and vibrational overlaps. Because the electronic
excited states of interest are strongly coupled, the simulation of absorption
and emission spectra is a non-trivial procedure for which quantum dynamics was
used to overcome the breakdown of the Born-Oppenheimer approximation. The low-
dimensionality of the model allows us to use the MCTDH (multiconfiguration
time-dependent Hartree) method for the propagation of vibrational wave packets
on the coupled electronic excited states, thus paving the way to further high-
dimensional simulations and/or other types of calculations including the role
of the solvent and/or of soft modes in a more effective way, which are works
in progress.

As usual in this context, the vibronic spectra are then obtained using the
Fourier transform of the autocorrelation functions obtained via the dynamics.
For the absorption spectra, the results are in very good agreement in spite of
the low-dimensionality of the models. However, the emission spectra show at
this stage two types of contributions: one that complies very well with the
observed Stokes shift of the band origin,
$\Delta\bar{\nu}\simeq$2000\text{\,}{\mathrm{cm}}^{-1}$$, and another that
does not and merely provides coincidence of both band origins. The Stokes-
shifted band is the one that corresponds with a coherent flip of light
polarization consistent with the odd character of the interstate coupling and
could be attributed to an emissive Herzberg-Teller contribution or even
perhaps a geometric phase; to the extent of our knowledge, this is not an
effect that has been reported in large molecules. However, this may reflect
some symmetry-driven selection rule concerning the conservation of orbital
momentum that is rarely to be mentioned in this context.

In the first part, we recall our main insights into the electronic structure
calculations on the most important points of the energy landscape of m22, and
we give relevant computational details associated to electronic structure and
quantum dynamics calculations for this work. We also describe the vibronic
coupling Hamiltonian models that were used in the quantum dynamics
calculations hereby. Next, we discuss the results of these models fitted on
accessible ab initio data and describe the vibronic eigenstates produced by
the models and their molecular-symmetry/light-polarization peculiarities.
Finally, computed absorption and emission spectra are provided and discussed,
and further compared to their experimental counterparts, before conclusions
and outlooks are given so as to open to further investigations.

##  II Theoretical background and methods

###  II.1 Conical Intersections and branching-space vectors

The lowest two excited states of the energy landscape of m22 have been studied
previously. ho_diabatic_2019 Their potential energy surfaces (PESs) share a
conical intersection (CoIn) seam, which lies locally within an
$(N-8)$-dimensional subspace (where $N=108$ is the number of internal
coordinates for the molecule). The local complementary subspace to the seam is
denoted the branching space, where the energy difference increases linearly
from zero along the two branching-space vectors. The minimum-energy conical
intersection (MECI) of the seam was characterizedho_diabatic_2019 and will be
used as a central point in this work. For high-symmetry points, the branching-
space vectors can be associated to simply defined quantities: the (halved)
gradient difference vector (GD) and the derivative coupling vector (DC). In
our case, the point group of the m22 molecule is C2v, and the GD and DC expand
along A1 normal modes (totally symmetric) and B2 normal modes (in-plane, non-
totally symmetric), respectively.

The usual Born-Oppenheimer approximation fully breaks down over regions such
as CoIn seams, which requires considering both electronic excited states
simultaneously together with their non-adiabatic couplings. The adiabatic
(eigenstate) representation of the corresponding electronic Hamiltonian is not
practical because of the singularity with respect to the energy difference in
the non-adiabatic couplings. Thus, a quasi-diabatic representation should be
preferred, with smoothly varying functions (with respect to the internal
coordinates) for the diabatic potentials and the inter-state couplings. In
this work, such quantities are, indeed, smoothly varying functions of the
internal coordinates associated to acetylenic and quinoidal normal modes of
vibration of m22 (fig. 1).

###  II.2 Computational details

#### Quantum chemistry

All quantum chemistry calculations were performed with the Gaussian16 package
(revision A.03)g16 using DFT (for the ground state) and TD-DFT (for the
excited states) at the CAM-B3LYP/6-31+G* level of theory. Such computational
choices were already validated for para-conjugated PPEadamo_calculations_2013
; ho_vibronic_2017 and this same meta-substituted molecule for the study of
its absorption spectrum, with good agreement to the
experiment.ho_diabatic_2019 The defect of TD-DFT for the evaluation of charge
transfer (CT) contributions to the excited states is not troublesome in the
excited state of interest in this work for m22 (the first two excited states
are essentially locally excited (LE) states).

The MECI of interest here was taken from the work of E. K. L. Ho _et al._
ho_diabatic_2019 , and has been optimized in the C2v subspace of the molecular
geometries. In the current TD-DFT implementation of the Gaussian package,
analytic non-adiabatic couplings (NAC) between excited states and thus the
derivative-coupling vector (DC) are not routinely available. For this reason,
the branching-space vectors are evaluated at the MECI using the energy-based
and wave-function free method of Gonon _et al._ , along with their lengths,
through the diagonalization of the Hessian of the squared energy
difference.gonon_applicability_2017 The main contribution in GD comes from the
quinoidal mode 81 (43% of the GD), the main contribution in DC comes from the
acetylenic mode 87 (73% of the DC). Because mode 87 mainly is the
antisymmetric elongation of left and right acetylene bonds (B2), its symmetric
counterpart, the acetylenic mode 88 (A1), is added to the description; it also
participates in the GD (10% of the GD). Both A1 modes 81 and 88 are important
for the description of the region between the Franck-Condon point (FC) and the
MECI. Indeed, both modes have non-negligible contributions to the gradients at
the Franck-Condon (FC) point and are associated to important shifts between
the FC point and the MECI ($d_{81}=-7.36\sqrt{m_{e}}a_{0}$ and
$d_{88}=6.78\sqrt{m_{e}}a_{0}$).

Rigid scans of the PESs of S0, S1, and S2 were obtained starting from the MECI
and following the normal modes directions computed at the minimum of the
electronic ground state. The grid points are described in table 1, and the
displacement associated to the chosen normal modes of vibration are
schematized in fig. 1.

Table 1: Frequency calculation results, symmetry, and generated points along modes 81, 87, and 88 originated from the MECI. Mode S0 (S1) |  Frequency (cm-1) | Reduced mass (AMU) | Symmetry |  Grid ($\sqrt{m_{e}}a_{0}$)
---|---|---|---|---
81 | 1656 | 5.465 |  A1 | 15 points [-18,18]
87 | 2367 (2184) | 11.998 (11.979) |  B2 | 15 points [-18,18]
88 | 2367 (2353) | 11.998 (11.997) |  A1 | 15 points [-18,18]

B2 displacements break C2v geometries to Cs geometries (exploring the minima
of the S1 PES) whereas A1 displacements conserve the C2v symmetry of the
molecule allowing for an assignation of the diabatic PES using unambiguous
symmetry arguments.

#### Quantum dynamics

The Quantics package was used for quantum dynamics
calculationsbeck_multiconfiguration_2000 ; worth_quantics_2015 , mainly for
the multiconfiguration time-dependent Hartree (MCTDH) wave packet propagation
method. Details concerning the single-particule functions (SPFs) and the
primitive basis set are given in table 2. The multi-set formulation of the
MCTDH ansatz was used, with different sets of SPFs (8 per state). The
primitive basis set consists in harmonic oscillator wave functions for the
three degrees of freedom. Two types of calculations were done: relaxations and
propagations. Relaxations are used for computing vibronic excited eigenstates
within our three-mode two-state model. As a result, vibronic states expand in
a basis of two electronic wave functions, which in most of this work are
delocalized diabatic electronic wave functions. Propagations are used for
computing autocorrelation functions and associated spectra (for absorption and
emission). The autocorrelation function is computed as:

$C(t)=\Braket{\Psi(0)}{\Psi(t)}=\Braket{\Psi(t/2)^{*}}{\Psi(t/2)}$ (1)

and the intensity of the spectra are:

$I(\omega)\propto\int_{-\infty}^{+\infty}C(t)\exp{(i\omega t)}\mathrm{d}t.$
(2)

Because the propagation time is finite, the autocorrelation function was
multiplied as usual with a decaying function

$g(t)=\cos^{n}(\pi t/2T),$ (3)

where $T$ is the propagation time and $n=1$ in our calculations. Additionally,
the autocorrelation function was multiplied by a Gaussian function
characterized by the damping time $\tau=(2\sqrt{2\log{2}})/(\Delta\omega)$ to
simulate the experimental full width at half-maximum (FWHF), $\Delta\omega$.

Table 2: Single-particle functions and primitive basis parameters for the relaxation and propagation of nuclear wave packets. Reduced masses for the HO primitive basis are 1 because the coordinates are mass-weighted. Coordinates | SPFs per state |  |  |  |
---|---|---|---|---|---
Q87 | multi-set 8,8,8 |  |  |  |
Q81 | multi-set 8,8,8 |  |  |  |
Q88 | multi-set 8,8,8 |  |  |  |
Coordinates | Basis Type | Size of the Basis | Eq. Position | Frequency | Reduced Mass
Q87 | HO | 15 | 0.0 |  2365 cm-1 | 1.0
Q81 | HO | 15 | 0.0 |  1655 cm-1 | 1.0
Q88 | HO | 15 | 0.0 |  2365 cm-1 | 1.0

For absorption spectra calculations (ground state to excited states) and for
vibronic excited states relaxation, the initial wave packet simply consists in
the zero-th order harmonic oscillator wave function with corresponding
frequencies for the three coordinates 81, 87, and 88. For emission spectra
calculations (excited states to ground state), the initial wave functions are
the vibronic excited states obtained from the relaxations in the electronic
excited states. The different contributions to the emission spectra are
discussed in the results section.

###  II.3 The Quadratic and Linear vibronic coupling Hamiltonian models (QVC,
LVC)

As stated above, three modes are chosen for modelling the potential energy
surfaces. They have a specific symmetry, which makes more natural the
definition of a diabatic vibronic coupling Hamiltonian
model.cederbaum_strong_1977 ; koppel_multimode_1984 ; cattarius_all_2001 Two
levels of description are considered here: a quadratic vibronic coupling (QVC)
model, relying on a second-order expansion of the diabatic potential energy
surfaces, or a linear vibronic coupling (LVC) model, tuning and coupling the
two diabatic potential energy surfaces through linear terms only. In the
following, electronic state 1 (notation superscript (1)) corresponds to the
singlet excited A1 state and electronic state 2 (notation superscript (2))
corresponds to the singlet excited B2 state. The diabatic representation of
the Hamiltonian in the QVC model for the three modes of the m22 molecule
reads:

$\displaystyle\hat{H}^{\text{dia}}$
$\displaystyle=\left(\hat{T}_{\text{nu}}(\bm{Q})\right)\mathbb{1}_{2}+\begin{bmatrix}E^{(1)}(\bm{Q}=0)&0\\\
0&E^{(2)}(\bm{Q}=0)\end{bmatrix}$ (4)
$\displaystyle+\frac{1}{2}\begin{bmatrix}k_{81}^{(1)}Q_{81}^{2}+k_{88}^{(1)}Q_{88}^{2}+k_{87}^{(1)}Q_{87}^{2}&0\\\
0&k_{81}^{(2)}Q_{81}^{2}+k_{88}^{(2)}Q_{88}^{2}+k_{87}^{(2)}Q_{87}^{2}&\end{bmatrix}$
$\displaystyle+\begin{bmatrix}\kappa_{81}^{(1)}Q_{81}+\kappa_{88}^{(1)}Q_{88}&0\\\
0&\kappa_{81}^{(2)}Q_{81}+\kappa_{88}^{(2)}Q_{88}\end{bmatrix}$
$\displaystyle+\begin{bmatrix}0&\lambda_{87}Q_{87}\\\
\lambda_{87}Q_{87}&0\end{bmatrix}+\begin{bmatrix}\gamma_{81,88}^{(1)}Q_{81}Q_{88}&0\\\
0&\gamma_{81,88}^{(2)}Q_{81}Q_{88}\end{bmatrix}$
$\displaystyle+\begin{bmatrix}0&\mu_{87,81}Q_{87}Q_{81}+\mu_{87,88}Q_{87}Q_{88}\\\
\mu_{87,81}Q_{87}Q_{81}+\mu_{87,88}Q_{87}Q_{88}&0\\\ \end{bmatrix}.$

In this expression:

  1. 1.

the first matrix represents the Kinetic Energy Operator (KEO), identical for
each electronic state;

  2. 2.

the second matrix represents the adiabatic energies (vertical excitation
energies) at the Franck-Condon point (which is the reference point for the
expansion);

  3. 3.

the third matrix represents the distortion matrix in the framework of the
Franck-Condon principle and of the harmonic approximation for vibronic
transitions;

  4. 4.

the fourth matrix represents the first-order description of the MECI with
respect to the tuning modes;

  5. 5.

the fifth matrix represents the first-order description of the MECI with
respect to the coupling mode;

  6. 6.

the sixth matrix represents the second-order description of the MECI with
respect to the tuning modes, corresponding to the _primary Duschinsky_
transformation or mode mixing between modes of same symmetry;

  7. 7.

the seventh matrix represents the second-order description of the MECI with
respect to the coupling mode, corresponding to the _secondary Duschinsky_
transformation or mode mixing between modes of different symmetry.

In this context, we call _primary Duschinsky_ parameters the second-order
coupling parameters between normal modes of identical symmetry (81 with 88)
and _secondary Duschinsky_ parameters the second order coupling parameters
between normal modes of different symmetry (87 with 81 or 88). The ab initio
Duschinsky coefficients obtained for the correlation between the first
electronic excited state and the electronic ground state are given in the
section V.1, along with the shift vector.

Such a full-QVC 3D model should grasp most of the physics of the problem, with
a total of 17 parameters. However, as will be pointed out in the results
section, further simplifications can be made without significant loss of
information. Hence, in the parent LVC model, all bilinear terms are neglected,
such that the Hamiltonian reads:

$\displaystyle\hat{H}^{\text{dia}}$
$\displaystyle=\left(\hat{T}_{\text{nu}}(\bm{Q})\right)\mathbb{1}_{2}+\begin{bmatrix}E^{(1)}(\bm{Q}=0)&0\\\
0&E^{(2)}(\bm{Q}=0)\end{bmatrix}$ (5)
$\displaystyle+\frac{1}{2}\begin{bmatrix}k_{81}^{(1)}Q_{81}^{2}+k_{88}^{(1)}Q_{88}^{2}+k_{87}^{(1)}Q_{87}^{2}&0\\\
0&k_{81}^{(2)}Q_{81}^{2}+k_{88}^{(2)}Q_{88}^{2}+k_{87}^{(2)}Q_{87}^{2}&\end{bmatrix}$
$\displaystyle+\begin{bmatrix}\kappa_{81}^{(1)}Q_{81}+\kappa_{88}^{(1)}Q_{88}&0\\\
0&\kappa_{81}^{(2)}Q_{81}+\kappa_{88}^{(2)}Q_{88}\end{bmatrix}$
$\displaystyle+\begin{bmatrix}0&\lambda_{87}Q_{87}\\\
\lambda_{87}Q_{87}&0\end{bmatrix}$

In both expansions of the model, diagonal terms are the diabatic potential
energies associated to delocalized electronic states, with specified symmetry.
Thus, the electronic states 1A1 and 1B2 are coupled through the
$\text{B\textsubscript{2}}=\text{A\textsubscript{1}}\bigotimes\text{B\textsubscript{2}}$
normal mode of vibration, 87. The eigenvalues of the model Hamiltonian are to
be compared with the adiabatic energies of the first two electronic excited
states computed with _ab initio_ quantum chemistry calculations.

###  II.4 Evaluating the parameters of the LVC and QVC models

In the present L/QVC models, the diabatic states 1 and 2 have distinct
symmetries when $Q_{87}=0$, that is in the C2v subspace of the molecular
geometries. This subspace is partly described by the profiles of 81 and 88,
and energies from quantum chemistry calculation can be assigned to either the
A1 or the B2 diabatic states (which fully identify to either adiabatic states
S1 or S2 depending on the side of the MECI the calculation was made). Thus,
all parameters except for the one related to the mode 87 can be optimized by
fitting directly the diabatic potential energies onto the symmetry-assigned
potential energy surfaces from _ab initio_ calculations, such that
$\\{H_{11}(Q_{81},Q_{88}),H_{22}(Q_{81},Q_{88})\\}$ are simultaneously fitted
to the data
$\\{E_{\text{A\textsubscript{1}}},E_{\text{B\textsubscript{2}}}\\}$ when
$Q_{87}=0$. The fitting procedure uses the minimization of the squared-
residual functions defined as:

$\displaystyle
L(E^{(1)},E^{(2)},\bm{\kappa},\bm{\gamma},\bm{k}_{81},\bm{k}_{88})=$ (6)
$\displaystyle\sum_{n,Q_{87}(n)=0}\left(H_{11}[E^{(1)},\bm{\kappa}^{(1)},\bm{\gamma}^{(1)}](\bm{Q}(n))-E_{\text{A\textsubscript{1}}}(n)\right)^{2}$
$\displaystyle+$
$\displaystyle\sum_{n,Q_{87}(n)=0}\left(H_{22}[E^{(2)},\bm{\kappa}^{(2)},\bm{\gamma}^{(2)}](\bm{Q}(n))-E_{\text{B\textsubscript{2}}}(n)\right)^{2}.$

Then, taking the resulting parameters as granted, the eigenvalues of the full
LVC and QVC models are fitted to the data from maps (87,81) and (87,88) so as
to optimize the parameters related to mode 87. Again, the parameters are
obtained through a least-square fitting procedure minimizing the function:

$\displaystyle L(\lambda_{87},\bm{\mu},\bm{k}_{87})=$ (7)
$\displaystyle\sum_{n}\left(V_{1}[\lambda_{87},\bm{\mu},\bm{k}_{87}](\bm{Q}(n))-E_{\text{S\textsubscript{1}}}(n)\right)^{2}$
$\displaystyle+$
$\displaystyle\sum_{n}\left(V_{2}[\lambda_{87},\bm{\mu},\bm{k}_{87}](\bm{Q}(n))-E_{\text{S\textsubscript{2}}}(n)\right)^{2},$

where $V_{1}(\bm{Q})$ and $V_{2}(\bm{Q})$ are the eigenvalues of the model
Hamiltonians.

For the LVC model, only separate profiles along 81, 87, and 88 are required,
whereas the QVC model requires 2D-surfaces (87,81), (87,88), and (81,88) to
optimize the second-order parameters. Both yield parameters that allow the
reconstruction of the PESs in the 3D-space but the LVC model does not account
for full correlation between the coordinates, which will be discussed in the
following section. The fitting procedure with the determination of the L/QVC
parameters was done using the functions of the SciPy library in Python3.4.

##  III Results and discussion

###  III.1 _Ab initio_ calculations and LVC parameters

In this section we focus on the use of the LVC model only, for which the
parameters are given in table 3. Note that although the fitting procedure was
carried out from data in atomic units and corresponding to mass-weighted
coordinates, the parameters are given in understandable equivalent quantities,
such as frequencies (actually, wavenumbers) and geometrical shifts instead of
curvatures and gradients, respectively.

Table 3: LVC parameters obtained upon fitting ab initio calculations. Parameters associated to the first-order expansion ($\lambda_{i}$ or $\kappa_{i}$) are given by the associated characteristic shift $d_{i}=\pm\frac{\lambda_{i}}{k_{i}}\text{ or }-\frac{\kappa_{i}}{k_{i}}^{(k)}$ in the mass-weighted framework for displacement along the normal modes directions. Parameter | Equivalent | Value | Unit
---|---|---|---
$E^{(1)}$ | – | 4.405 | $\mathrm{eV}$
$E^{(2)}$ | – | 4.380 | $\mathrm{eV}$
$\kappa_{81}^{(1)}$ | ${d_{81}^{(1)}}$ | 7.900 | $a_{0}\sqrt{m_{e}}$
$\kappa_{81}^{(2)}$ | ${d_{81}^{(2)}}$ | 3.094 | $a_{0}\sqrt{m_{e}}$
$\kappa_{88}^{(1)}$ | ${d_{88}^{(1)}}$ | -7.205 | $a_{0}\sqrt{m_{e}}$
$\kappa_{88}^{(2)}$ | ${d_{88}^{(2)}}$ | -8.274 | $a_{0}\sqrt{m_{e}}$
$\lambda_{87}$ | $d_{87}$ | 7.417 | $a_{0}\sqrt{m_{e}}$
$k_{81}^{(1)}$ | $\omega_{81}^{(1)}$ | 1515 | ${\mathrm{cm}}^{-1}$
$k_{81}^{(2)}$ | $\omega_{81}^{(2)}$ | 1552 | ${\mathrm{cm}}^{-1}$
$k_{88}^{(1)}$ | $\omega_{88}^{(1)}$ | 2274 | ${\mathrm{cm}}^{-1}$
$k_{88}^{(2)}$ | $\omega_{88}^{(2)}$ | 2281 | ${\mathrm{cm}}^{-1}$
$k_{87}^{(1)}$ | $\omega_{87}^{(1)}$ | 2200 | ${\mathrm{cm}}^{-1}$
$k_{87}^{(2)}$ | $\omega_{87}^{(2)}$ | 2201 | ${\mathrm{cm}}^{-1}$

Figure 2 shows the adiabatic energies of states S1 and S2 from ab initio
calculations and the potential energy profiles computed with the latter
optimized parameters. The profiles for modes 81 and 88 are analogous.
Displacements along these A1 normal modes lead to neither even nor odd
potential energy surfaces with respect to the associated A1 coordinates. These
displacements preserve the C2v geometry and thus the delocalized quasi-
diabatic electronic states (based on symmetry labelling) maintain their
nature. Thus, we can assign ab initio adiabatic energies
($E(\text{S\textsubscript{1,2}})$ symbol + in fig. 2, a and b) directly to the
diabatic potentials ($E^{(1)}=E(\text{A\textsubscript{1}})$ and
$E^{(2)}=E(\text{B\textsubscript{2}})$) of the LVC model.

Figure 2: Adiabatic energies from quantum chemistry calculations (symbol
$\bm{+}$) and diabatic potential energies from the delocalized model (blue and
red plain lines) along the 81th (a), 88th (b) and 87th (c) normal modes of
vibration. Diabatic potential energies from the localized model are also given
for 87 (c, purple dashed lines). For unspecified coordinates, the values are
those at the MECI. All coordinates are mass-weighted and given in atomic
units.

Note that displacements along the normal mode 81 lift more efficiently the
degeneracy than along 88, which is consistent with the fact that the gradient
difference vector GD overlaps more with the direction associated to the normal
mode 81 than 88. This assignment is not so simple for the profile along the
mode 87, where there are two relevant quasi-diabatic representations. For mode
87, the eigenvalues of the LVC model Hamiltonian are directly fitted to the
adiabatic energies from ab initio calculations. Thus, the symmetry-labelled
quasi-diabatic representation is not necessarily the most relevant for
comparison with the ab initio calculations. A more relevant representation is
the localized representation fitting left and right wells to the S1 potential
energy surface, which is obtained through half-sum and difference of the
symmetry-labelled quasi-diabatic states.ho_diabatic_2019

Finally, these parameters give a three-dimensional picture in the vicinity of
the MECI, with exact description of the MECI (starting point of our
calculations) and accurate description of the minima in the first electronic
excited state. In fig. 3, we give the LVC model PESs of the first two
electronic excited states in the planes (87,81), (87,88), and (81,88). The
minima of S1 are correctly reproduced, with $E=$4.21\text{\,}\mathrm{eV}$$
which is in a good agreement with the ab initio calculations
($E=$4.12\text{\,}\mathrm{eV}$$). The deviation from the ab initio optimized
minima of S1 is reasonable considering the following: (i) the starting point
for the rigid scans is a MECI and not one of the optimized minima, and (ii)
the explored space is only a restricted three-dimensional subspace within all
possible molecular geometries. The positions of the minima are consistent with
the shift vector, given in section V.1, with two equivalent minima with
respect to the non-totally symmetric coordinate 87. The positions are also in
good agreement for the totally symmetric coordinates 81 and 88, with the
minima being in the positive region of 81 (quinoidal elongation) and negative
region of 88 (acetylenic elongation). The comparison of the critical points
energy in the ab initio PES and in the parameterized PES is summarized in
table 4. The most important deviation is found for the vertical transition
energies and for the minimum of the S1 state.

Figure 3: LVC potential energy surfaces (in $\mathrm{e}\mathrm{V}$) of electronic states S1 (a) and S2 (b) in the planes (87,81), (87,88), and (81,88) from left to right. For unspecified coordinates, the values are those at the MECI. All coordinates are mass-weighted and given in atomic units, $\sqrt{m_{e}}a_{0}$. Table 4: Energies of the critical points in the ab initio PESs (CAM-B3LYP/6-31+G* level of theory from Ho _et al._ ho_diabatic_2019 ) and in the parameterized LVC model PESs.  Critical Point | MinS1 |  TSB2 |  TSA1 | $E_{S\textsubscript{1}}(\bm{0})$ | $E_{S\textsubscript{2}}(\bm{0})$ | MECI
---|---|---|---|---|---|---
CAM-B3LYP/6-31+G* | 4.12 | 4.25 | 4.29 | 4.43 | 4.47 | 4.29
LVC Model | 4.209 | 4.288 | 4.272 | 4.380 | 4.405 | 4.29

As mentioned earlier, the LVC model is limited for reproducing the correlation
between the coordinates, in particular for taking into account the Duschinsky
mixing of the modes between the ground state and the excited states. The
importance of reproducing this correlation with the QVC model is discussed in
the next section.

###  III.2 Comparison of linear and quadratic VC models

The LVC and QVC models show almost identical values for the parameters that
they share and no significant effects of the extra second-order parameters of
the QVC model. Similar figures and tables as the one provided for the LVC
model earlier are given in the Supplementary Information (SI) for the QVC
model, see tab. SI2 and fig. SI1.

The main difference lies in the ellipsoidal shapes of the isoenergy slices
through the potential energy surfaces with respect to the direction of the A1
normal mode 88. Indeed, with the linear description only (LVC), the ellipsoids
are aligned with the directions of the ground state normal modes of vibration
(chosen for the generation of the grids), which is in contradiction with the
Duschinsky transformation of the normal modes from the electronic ground state
to excited states. This transformation is taken into account in the quadratic
description (through bilinear terms), for which the ellipsoids have different
main axes, as shown in fig. 4.

Figure 4: Potential surfaces of the first excited state in the LVC (left) and
QVC (right) models with respect to the acetylenic coordinates 87 and 88.
Central arrows correspond to the directions of the main axes of the
ellipsoids. All coordinates are mass-weighted and given in atomic units,
$\sqrt{m_{e}}a_{0}$.

However, there is no significant effect of this transformation on the physics
of the (dimensionally-reduced) system and in particular for the absorption and
emission spectra. This can be explained by the almost equivalent frequencies
for modes 87 and 88, which make the rotation of the potential wells chemically
motivated (pseudofragmentation, consistent with left or right local modes on
each tolane moiety) but physically (or, say, numerically) not significant.

###  III.3 Vibronic states in the excited-state landscape

To calculate and interpret the details of vibronic transitions within UV-
visible absorption and emission spectra, we must access the vibronic
eigenstates of the molecule and their eigenenergies. In the electronic ground
state, the vibronic states are simple because the Born-Oppenheimer
approximation is essentially valid there. We describe those states as simple
products of the electronic ground state wave function and the 3D harmonic
oscillator wave functions, for which the parameters are the frequencies and
reduced masses of the normal modes of vibration 81, 87, and 88 in the
electronic ground state. For the excited states, the description of the
vibronic eigenstates is more involved. We are in fact facing a situation
whereby the Born-Oppenheimer approximation is completely off. Due to the
quasi-degeneracy of the two states of interest, already in the Franck-Condon
region, and the presence of an accessible conical intersection in its close
vicinity and inducing a strong diabatic mixing, the first two excited
adiabatic electronic states are strongly coupled nondiabatically and cannot at
all be treated separately for the description of the associated vibronic
states. In other words, and using the language of spectroscopy, we are in the
unusual situation where the first-order Herzberg-Teller effect along mode 87
is clearly not a perturbation.

Within a description beyond Born-Oppenheimer the vibronic states expand as
linear combinations of products of electronic and nuclear wave functions. We
choose the delocalized, symmetry-labelled, diabatic electronic basis set, as
described earlier for the LVC Hamiltonian model. Viewed as a good
approximation of a crude adiabatic basis for the Herzberg-Teller effect, it
presents the great advantage that transition dipoles will not vary much with
displacements from C2v molecular geometries (such as the Franck-Condon ground-
state equilibrium geometry); in particular, A1 is polarized along $z$ and B2
along $y$. Within this representation, we thus get:

$\displaystyle\Psi^{\text{(exc)}}_{k}(\bm{q},\bm{Q})$
$\displaystyle=\chi_{A_{1},k}(\bm{q})\phi_{A_{1}}(\bm{q};\bm{Q})$ (8)
$\displaystyle+\chi_{B_{2},k}(\bm{q})\phi_{B_{2}}(\bm{q};\bm{Q})$

where $\bm{Q}=(Q_{81},Q_{87},Q_{88})$ are the internal coordinates and
$\bm{q}$ are the electronic coordinates. The electronic states
$\phi_{A_{1}}(\bm{q};\bm{Q})$ and $\phi_{B_{2}}(\bm{q};\bm{Q})$ are quasi-
diabatic states, which vary moderately and smoothly (ideally not at all to
first order) with the internal coordinates $\bm{Q}$. The associated diabatic
state populations are simply the integral of the vibrational contribution over
the internal coordinates $\bm{Q}$. In this work, such vibronic eigenstates are
obtained through quantum dynamics calculations on the “coupled” PESs from the
LVC Hamiltonian model. The calculation consists in the propagation in
imaginary time of an initial nuclear wave packet until variational
convergence. The so-called relaxation method is analogous to variational
energy minimisation (parametric optimisation). We consider two choices of
initial states (they will be the same as those used later on for the
absorption spectra simulations). They are the vibrational ground state for the
electronic ground singlet state (fundamental of the 3D harmonic oscillator for
81, 87, and 88) multiplied by either the diabatic state A1 or the diabatic
state B2 (coupled electronic excited singlet states).

After an imaginary relaxation time of about $20\text{\,}\mathrm{fs}$, which is
short and indicates little internal vibrational redistribution, two quasi-
degenerate ($\Delta E=$0.0076\text{\,}\mathrm{eV}$$) but distinct vibronic
states are obtained (they could be viewed as a tunnelling pair, symmetric and
antisymmetric, in an adiabatic context). Their characteristics are collected
in table 5. We notice that both vibronic states are complementary. Vibrational
contributions of the vibronic states in the planes (87,81), (87,88), and
(81,88) are given (see fig. 5 for the lowest, and fig. SI3 for the second
vibronic state). The first vibronic state is mainly B2-populated (80%) while
the second one is mainly A1-populated (78%) which is consistent with the ab
initio energy and gradient information regarding the two delocalized diabatic
states. Indeed the minimum of B2 is lower in energy than the minimum of A1 in
our 3D model, which is in good agreement with the energy ordering of the
optimized transition states A1 and B2 (see table 4). These results also
indicate that with this model, in each vibronic state, about 80% of the
initial wave packet is preserved in the stationary vibronic states, which
represents about 80% of totally symmetrical vibrational contributions with
respect to the three coordinates (see fig. 5, (a)). As a result, we notice
about 20% of population transfer into the second diabatic state, not-initially
targeted, allowed by the coupling coordinate 87. This 20% transferred
population corresponds to an odd shape of the vibrational contributions with
respect to the coupling coordinate 87 (see fig. 5, (b)). We stress here that
the part of the nuclear wave packet that transfers from one diabatic state to
another must undergo a change in its parity or symmetry according to the odd
behavior of the diabatic interstate coupling function with respect to mode 87.
This would be termed a geometric (Berry or Longuet-Higgins) phase effect if
the adiabatic representation were used instead. This change is easily
identified in the delocalized and symmetry-labelled diabatic states basis set:
the nuclear wave packet turns odd when leaving its initial diabatic state, due
to the off-diagonal coupling being odd. Let us notice that the first two
vibronic states are distinct, but similar with only inversion of the roles of
the diabatic states A1 and B2.

Table 5: Vibronic states energy and diabatic populations. The transition energy is computed by setting the energy reference to the zero-point energy of the 3D model, $E=$0.396\text{\,}\mathrm{eV}$$ Vibronic state | Total Energy (eV) | Transition Energy (eV) | $\mathcal{P}_{A_{1}}$ | $\mathcal{P}_{B_{2}}$
---|---|---|---|---
$\Psi^{\text{(exc)}}_{0}$ | 4.576 | 4.180 | 0.20 | 0.80
$\Psi^{\text{(exc)}}_{1}$ | 4.584 | 4.188 | 0.78 | 0.22
Figure 5: Nuclear wave packet contributions in diabatic electronic states B2
(a) and A1 (b) in the first vibronic state, in planes (87,81), (87,88), and
(81,88) from left to right. For unspecified coordinates, the values are those
at the MECI. All coordinates are mass-weighted and given in atomic units,
$\sqrt{m_{e}}a_{0}$. The blank figure corresponds to the wave functions having
the B2 coordinate 87 in the nodal plane. The second vibronic state is
described in fig. SI3

We propose to interpret again the vibronic states looking mostly at their
character with respect to the non-totally symmetrical coordinate 87. Figure 6
gives a schematic representation of the vibronic eigenstates in the ground and
excited electronic states. The pair of vibronic eigenstates in the electronic
excited states manifold can be seen as vibronic eigenstates in diabatic states
B2 and A1 interacting through the non-totally symmetrical coordinate 87. This
representation is analogous to an interaction diagram for electronic orbitals,
only with coupled vibrational orbitals here. The first two pairs of vibronic
eigenstates are the results of four vibrational orbitals interacting and
mixing their character. As a result, the lower vibronic eigenstate has two
vibrational contributions of different parity, odd and even with respect to
the coordinate 87.

Figure 6: Representation of the vibronic eigenstates as interacting diabatic
vibronic orbitals. The vibrational number $v$ refers to the excitation number
of non-totally symmetrical normal mode of vibration 87. The nuclear wave
packets are given in the plane (87,88).

For absorption, the initial wave packet in the ground electronic state
overlaps with the even contribution of each of the vibronic eigenstates. For
emission, two cases arise. The initial wave packets have even and odd
vibrational contributions, that overlap with the first and second vibronic
ground states, respectively. Two contributions to the emission spectrum are
thus obtained. One is sharing the $0-0$ band origin with absorption and one
shifted by the vibrational energy of the second vibrational state of the
electronic ground state. Analogous transitions occur for the second vibronic
eigenstate. The transition energies to the vibronic eigenstates are ordered as
the transition energies of the non-interacting diabatic vibronic eigenstates.
The latter ordering is consistent with the ab initio potential energy minima
of the electronic states in the C2v sub-space which are transition states of
the first excited state in the CS space of molecular geometries.

###  III.4 Absorption and emission spectra

In this section, we describe the procedure for producing absorption and
emission spectra (vibronic spectra) using quantum dynamics for our model of
the m22 molecule. The different contributions to the spectra will be discussed
and compared to the low-temperature experimental spectra from Chu and Pang
experiments. Vibronic spectra are computed through Fourier transformation (FT)
of the autocorrelation function associated to the propagation of the nuclear
wave packet in the electronic excited states (for absorption) or in the
electronic ground state (for emission). The fact that only the autocorrelation
function of the nuclear wave packet is required as a good approximation is
helped by the use of the delocalized and symmetry-labelled diabatic basis set
for the electronic states. Indeed, in this basis set, the transition dipole
moments have very little dependence on the nuclear coordinates (Franck-Condon-
type approximation):

$\displaystyle\bm{\mu}_{\text{A\textsubscript{1}}}(\bm{Q})$
$\displaystyle\simeq\bm{\mu}_{\text{A\textsubscript{1}}}(\bm{0})=$-1.8303$\hat{\bm{z}}(a_{0}\sqrt{m_{e}})$
(9) $\displaystyle\bm{\mu}_{\text{B\textsubscript{2}}}(\bm{Q})$
$\displaystyle\simeq\bm{\mu}_{\text{B\textsubscript{2}}}(\bm{0})=$3.9656$\hat{\bm{y}}(a_{0}\sqrt{m_{e}}).$

In other words, we can safely neglect Herzberg-Teller terms in the symmetry-
adapted delocalized diabatic representation, which, in contrast, would be
unusually large with respect to variations in the coordinate 87 in the
strongly mixed adiabatic representation, as already pointed out. As a
consequence, the intensity spectra are here simply proportional to the FT of
the autocorrelation function of the propagated wave function. Here, the wave
functions must be propagated in real-time to access the vibronic spectra, as
opposed to the previous section where propagation was done in imaginary time
to access relaxed vibronic states.

#### Modelling absorption

For absorption and the initial wave packet, we choose the vibrational ground
state (fundamental of the 3D harmonic oscillator for 81, 87, and 88)
multiplied by either the diabatic state A1 or the diabatic state B2 as for the
vibronic states relaxations. The choices of A1 or B2 has little or no
influence on the shape of the absorption spectrum obtained in the end, even
though they correspond to distinct polarization processes.

Before discussing the absorption spectrum, let us comment the early dynamics
(first $25\text{\,}\mathrm{fs}$ of the $200\text{\,}\mathrm{fs}$ of simulation
time) of the real-time propagation of the initial wave packets in the
electronic excited states. We discuss first the propagation of the
A1-populated initial state. The trajectories of the centers (coordinate
expectation values) of the wave packets in the (81,88) plane for the two
diabatic states are given in fig. 7 (a), where the PESs of the A1 and B2
states are shown, along with the positions of the conical intersection seam
and the MECI. We also report diabatic and adiabatic populations (computed
using the Quantics package), fig. 7 (b). First, let us notice that the initial
state has a population of about 0.6 and 0.4 for the adiabatic states S1 and
S2. This almost 50/50 population of the first two adiabatic states is not
troublesome because, in the FC region, their energy difference is small
($\Delta E=$0.04\text{\,}\mathrm{eV}$$ at our level of theory) and both are
optically bright (oscillator strengths $f_{1}=1.71$, $f_{2}=0.37$ in atomic
units). Besides, in our 3D model, the ordering of the diabatic states B2 and
A1 corresponding to adiabatic states S1 and S2 respectively is consistent with
the ab initio data. It indicates the importance of taking into account both A1
modes 81 and 88 for modelling m22 as the conical intersection seam and the FC
gradients are oblique within the (81,88) plane. Using a minimal model on the
acetylenic vibrations with only modes 87 and 88 would initiate the propagation
on the wrong side of the conical intersection seam, with an ordering of the
diabatic states that is wrong compared to the ab initio data. The multiple
crossings of the conical intersection seam in the early dynamics stress the
importance of having a consistent ordering of the diabatic states at the FC
point and thus the importance of taking into account the coordinate 81. On a
final note, the numerical kink of the trajectory in the B2 state is
unconsequential because of the almost-zero diabatic population of the state B2
at this time of the simulation, and seems to be of no relevance for resuming a
satisfactory trajectory for times with non-zero diabatic population. Analogous
interpretations of the early dynamics can be done for the other choice of
initial diabatic state, with a B2-initially populated state, for which figures
are given in fig. SI2.

Figure 7: Early dynamics: trajectories of the center of the wave packets in
the (81,88)-plan. Diabatic potential energy surfaces A1 and B2 in blue and red
respectively, adiabatic PESs S1 and S2 in dashed lines and plain lines
respectively (a). Diabatic (plain lines) and adiabatic (dashed lines)
populations associated to these trajectories (b). Propagation for an initial
state on A1 is given here, propagation for an initial state on B2 is given in
fig. SI2

The two propagations performed here yield two absorption spectra according to
the populated excited diabatic state and corresponding light polarization. Let
us notice that two representations are always given for each of our
theoretical spectra. A representation obtained with FT on the autocorrelation
damped with $\tau=$100\text{\,}\mathrm{fs}$$ and another one one with FT on an
autocorrelation function damped with $\tau=$19\text{\,}\mathrm{fs}$$. The
first representation resembles to a stick-spectrum while the second one shows
realistically broadened bands. The damping time $\tau$ was chosen so as to
fulfil the experimental spectra, as described in the last sub-section. With
our model, we access two distinct but almost identical absorption spectra
(fig. 8). This is consistent with the results obtained for the vibronic states
relaxation. If analogous vibronic states are found from the two different
initial states, then absorption spectra must also be analogous because the
same eigenstates are accessible for the system. Here we discard for now the
differences in oscillator strength between the two diabatic states, but as
both are bright, both simulations and both absorption spectra are relevant.

Figure 8: Simulated absorption spectra from electronic ground state to
electronic excited state B2 (plain lines) and A1 (dashed lines).

Notice that within the 3D LVC model, the $0-0$ band corresponds to a
transition at $296\text{\,}\mathrm{nm}$. The second most intense band is at
$278\text{\,}\mathrm{nm}$, which corresponds to a difference of
$\Delta\bar{\nu}=$2190\text{\,}{\mathrm{cm}}^{-1}$$. The latter energy
difference is expected for the acetylenic mode 87, optically active in the
absorption, for which the frequency is about $2200\text{\,}{\mathrm{cm}}^{-1}$
in the 3D LVC model of the electronic excited states. The other transition at
$284\text{\,}\mathrm{nm}$, corresponding to
$\Delta\bar{\nu}=$1430\text{\,}{\mathrm{cm}}^{-1}$$ with respect to the $0-0$
band, is the expected transition involving the quinoidal mode 81. The
absorption spectra are also reported in fig. 10 for completeness of the
theoretical results and in fig. 11 for comparison with the experiments.

#### Modelling the emission spectrum

The emission spectrum is simulated by propagating the initial nuclear wave
packet in the electronic ground state, for which the PES is described as a 3D
harmonic potential with frequency parameters for mode 81, 87, and 88 extracted
from ab initio calculations at the minimum of the electronic ground state (see
table 1). For the initial state, the choice is more complicated than for
simulating absorption. Indeed, as mentioned earlier, the first two vibronic
states in the twofold excited electronic structure are quasi-degenerate and
share similar characteristics in vibrational shape and populations. The
assumption that the emission would occur only from the lowest vibronic state
should not be completely valid, but we focus first on this case as the second
one is analogous. Emitting from the lowest excited vibronic state (initially
B2-populated), we identify two simple initial states for the quantum dynamics
which will result in the end to two contributions to the emission spectrum.
The two initial states are: $\chi_{A_{1},0}\phi_{0}$ and
$\chi_{B_{2},0}\phi_{0}$. As described earlier the first vibrational part
$\chi_{A_{1},0}\phi_{0}$ is odd with respect to coordinate 87 while the second
one $\chi_{B_{2},0}\phi_{0}$ is even with respect to coordinate 87 (see fig.
5).

For the first vibronic eigenstate, two contributions to the emission spectrum
are obtained (see fig. 9). The $0-0$ band for emission is recovered from the
contribution of the vibrational part on B2, at $297\text{\,}\mathrm{nm}$. A
vibronic progression ($319\text{\,}\mathrm{nm}$ and $345\text{\,}\mathrm{nm}$)
analogous to absorption is obtained as expected for frequencies in the
acetylenic regions, with
$\Delta\bar{\nu}=$2320\text{\,}{\mathrm{cm}}^{-1}2360\text{\,}{\mathrm{cm}}^{-1}$$.
Such a contribution to the emission spectrum will be called a non-Stokes
contribution in the following. The contribution of the vibrational part on A1,
on the contrary, is shifted in the red wavelengths with an origin at
$319\text{\,}\mathrm{nm}$, corresponding to a Stokes shift
$\Delta\bar{\nu}=$2320\text{\,}{\mathrm{cm}}^{-1}$$, and maintains the
expected vibronic progression from this origin band to larger wavelengths.
Such a contribution to the emission spectrum will be called a Stokes
contribution in the following. The positions of the first three transitions
and the associated shifts with respect to the most intense band are summarized
in table 6.

Spectrum contribution |  $\lambda_{\text{max}}$ ($\mathrm{nm}$) |  First vibronic progression (${\mathrm{cm}}^{-1}$)
---|---|---
Absorption $\rightarrow$ |  296, (284), 278 |  –, (1430), 2190
Emission $\Rsh$ |  |
$\Psi_{0}$, non-Stokes |  297, (312), 319 |  –, (1620), 2320
$\Psi_{0}$, Stokes |  319, (337), 345 |  –, (1670), 2360
$\Psi_{1}$, non-Stokes |  296, (311), 318 |  –, (1630), 2340
$\Psi_{1}$, Stokes |  319, (336), 344 |  –, (1590), 2280
Table 6: Spectroscopic data for theoretical absorption and emission spectra,
before any comparison with the experiment. One row is given for absorption as
both computed spectra are identical for reading the transition energies. The
four contributions to the emission spectrum are given next. The most intense
transition for each contribution is given in bold. The first next transition
with low intensity is given in parenthesis. Associated vibronic progression
are given with respect to the most intense transition of each contribution.

Let us notice that all these results are expected given the form of the
vibrational contributions in the vibronic states. Indeed, with a time-
independent reasoning, the vibrational part $\chi_{B_{2},0}$ of the vibronic
state is expected to overlap with the vibrational ground state $\chi_{0}$
resulting in an band origin in the same energy region than for the band origin
of the absorption. However, the vibrational part $\chi_{A_{1},0}$ of the
vibronic state is expected to overlap with the first vibrational excited state
$\chi_{1}$ resulting in the shifted origin band at $319\text{\,}\mathrm{nm}$.
The same interpretation holds for the second vibronic state, only
interchanging the roles of A1 and B2 diabatic states.

Figure 9: Simulated emission spectra from vibrational contributions of B2
(plain lines) and A1 (dashed lines) of the first vibronic eigenstate to the
electronic ground state. Simulated emission spectra from the second vibronic
eigenstate are given in fig. SI4.

All the transitions observed between the ground vibronic states and the
excited vibronic states are schematized in the interaction diagram fig. 6, all
the resulting contributions to the absorption and emission spectra (two for
each vibronic state) are collected altogether in fig. 10.

Figure 10: Absorption spectra from the electronic ground state to both
diabatic states (blue lines) and emission spectra from the different
contributions of the first vibronic eigenstate (red lines). Simulated emission
spectra from the second vibronic eigenstate are given in fig. SI5.
Realistically broadened bands are obtained using a damping time
$\tau=$19\text{\,}\mathrm{fs}$$.

#### Discussion and comparison to the experiment

We chose to compare our theoretical results to the experiments of Chu and Pang
because of the exact same substitution of the molecule in both studies and
because of the temperature studychu_vibronic_2004 . Indeed, the measurement of
the spectra was done at
$25\text{\,}\mathrm{\SIUnitSymbolCelsius}-108\text{\,}\mathrm{\SIUnitSymbolCelsius}-198\text{\,}\mathrm{\SIUnitSymbolCelsius}$
and we choose the latter measurement for comparison. The experimental spectra
at $-198\text{\,}\mathrm{\SIUnitSymbolCelsius}$ are reproduced in fig. 11,
with courtesy of Chu and Pang, 2004chu_vibronic_2004 . The theoretical spectra
(considering the absorption to and emission from the first vibronic state) are
superimposed to the experiment, shifting the $x$-axis of the theoretical
spectra so that the band origins from theoretical absorption
($297\text{\,}\mathrm{nm}$) and experimental absorption
($304\text{\,}\mathrm{nm}$) match which represents a shift of
$7.3\text{\,}\mathrm{nm}$ at the $0-0$ band of absorption. The shift applied
to our theoretical spectra for the $0-0$ band of absorption to match with the
experiment is thus $\Delta E=$0.10\text{\,}\mathrm{eV}$$ and
$\Delta\bar{\nu}=$818\text{\,}{\mathrm{cm}}^{-1}$$. This difference between
the experiment and our model can be understood as an underestimated zero-point
energy difference for the minimal 3D model, because other vibrations should be
included to the complete description of the system. Besides, our model is
built on data at a TD-DFT level of theory, which may imperfectly reproduce
vertical transition energies. We stress that we shift the $x$-axis of all our
theoretical spectra of the same values of $818\text{\,}{\mathrm{cm}}^{-1}$ and
thus do not introduce new or artificial Stokes shifts, between emission and
absorption.

Figure 11: Experimental absorption and emission spectra (dashed lines),
courtesy of Chu and Pang, 2004chu_vibronic_2004 , reproduced in (a) and (b).
Theoretical absorption spectrum to the first vibronic eigenstate reproduced in
(a) and (b), theoretical emission spectrum from the first vibronic eigenstate,
contributions of B2 (a) and of A1 (b). Realistically broadened bands are
obtained using a damping time $\tau=$19\text{\,}\mathrm{fs}$$. Simulated
absorption and emission spectra from the second vibronic eigenstate are given
in fig. SI6.

The broadening of the theoretical spectra was chosen so that the FWHM
$\Delta\omega$ of the origin band for the absorption is the same as in the
experiment. To ensure this, we damp the autocorrelation function with the
damping time $\tau=$19\text{\,}\mathrm{fs}$$.

For absorption, the vibronic progression is well reproduced with our model for
transitions involving acetylenic vibrations. The missing transitions in-
between ($294\text{\,}\mathrm{nm}277\text{\,}\mathrm{nm}$) should correspond
to transitions involving other vibrations, not taken into account in our
minimal 3D model, such as quinoidal and triangular vibrations, expected and
identified in the absorption spectra of p2 and m22 in previous
workho_vibronic_2017 . In the emission, the Stokes shift contribution is
always identified as the contribution from the vibrational part of the
vibronic states that underwent a transfer from one diabatic state to another.
The complementary vibrational part shows no contribution to a spectrum with a
significant Stokes shift. In the case of the Stokes contribution, the
theoretical spectrum also reproduce the expected vibronic progression, with
the missing bands expected to result from transitions involving other
vibrations.

Strictly speaking, the spectrum to be compared to the experiment should be a
combination of the four theoretical contributions described, involving the
transition dipole moments for states A1 and B2 along with factors to take into
account the population of these two states in the vibronic states. However,
searching such a combination within this model and reasoning would never
result in a predominant Stokes contribution to the emission spectrum. Indeed,
the vibrational part of the vibronic states that overlaps with the $\chi_{1}$
ground vibrational state never correspond to more than 20% of the diabatic
population. In addition, the oscillator strength for the B2 state, in the
direction of $y$, is much greater than the oscillator strength of the other
state (along $z$). As it stands, our model is still limited by several
approximations on the system and preconceptions of how to simulate the
phenomenon under study. For the reasons mentioned above, it cannot explain why
the Stokes contribution dominates over the non-Stokes one, but it shows that
there is a rationale for the Stokes contribution that could be based on
symmetry considerations and thus may provide some tentative avenues for
further explorations.

##  IV Conclusions and Outlooks

In this work, we have devised diabatic linear and quadratic vibronic coupling
models of Hamiltonians for the study of unusually strongly nondiabatically
coupled excited states in a symmetric meta-substituted PPE (m22) and used them
for spectral simulations based on quantum dynamics. We showed that molecular
symmetry may play a major role in the unusual Stokes shift observed for this
system regarding its band origin; other effects are still under study, such as
the role of soft modes, solvent relaxation dynamics, and perhaps the effect of
possible competing photo-isomerizations.

On practical terms, our study is expected to serve two main purposes: laying
the basic foundations for a more general model of nonadiabatically coupled
PESs for symmetrical or asymmetrical meta- or tri-substituted PPEs, and using
it for providing a possible geometric explanation to the unusual spectroscopic
properties of m22 and its derivative species.

As pointed out, our modelling of m22 suggests some plausible rationale (based
on molecular vibronic symmetry and light polarization transfer or
conservation) for distinct different Stokes and non-Stokes contributions to
the emission spectrum, with no fully conclusive elucidation yet on why the
Stokes contributions should be predominant, as indeed observed in the
experiments.

It must be understood that modelling emission (as opposed to absorption)
spectra is not a routine task and potentially raises other fundamental
questions about how we should consider the very role and nature of the
stationary states involved in the experimental acquisition of absorption and
emission spectra against the duration of the excitation light source and the
time for the environment (external: solvent; internal: soft modes) to fully
relax compared to the time it takes for an excited molecule to emit light.

Within our description, for the Stokes contribution to overcome the non-Stokes
contribution, the effect of the coupling would have to be stronger, hence
fully forcing the electronic population to be affected quantitatively by the
non-totally symmetrical coordinate 87 and its odd character (according to the
linear behavior of the coupling with respect to it), along the way of a
geometric (Berry or Longuet-Higgins) phase effect such that the nuclear
density goes to the other side of the conical intersection with a nodal line.

This is a possibility, but it obviously is very sensitive to the level of
theory used for producing _ab initio_ data. However, in the extreme case of a
completely transferred population to a non-totally symmetrical wave packet,
the $0-0$ band in the absorption spectrum should also disappear, which thus
seems to raise some contradiction.

Such hypotheses are thus to be evaluated carefully in future works,
considering for example additional but relevant degrees of freedom for the
vibronic transitions observed experimentally in the steady-state spectra, or
dissipation effects (for example in the context of the hierarchical equations
of motion that we recently explored together with our colleagues on the same
system but within another context)jaouadi_laser-controlled_2022 . More
extensive simulations are currently under study and will be addressed soon in
future work.

Finally, it must be noted as a last hypothesis that the non-Stokes
contributions to the emission identified in this work seem to be associated to
some unconventionally ultrafast coherently polarized resonant Raman process,
while usual fluorescence spectra are typically viewed as incoherent, slow, and
spontaneous, and they should exhibit both Stokes and non-Stokes contributions.
Such a mystery certainly now awaits for further detailed comparisons with
time-resolved spectroscopic experimental studies.

## Acknowledgments

The authors warmly thank Y. Pang and Q. Chu for allowing the reproduction of
the experimental absorption and emission spectra of
1,3-bis(phenylethynyl)benzene (m22). The authors would also like to personally
thank M. Desouter-Lecomte and E. Mangaud for very fruitful discussions on the
system and the emission phenomenon, and their incentive for plausible
extensions to further studies with complementary viewpoints. J. G.
acknowledges the French MESR (Ministère de l’Enseignement supérieur, de la
Recherche) and the ENS (Ecole Normale Supérieure) of Lyon for funding his PhD
grant, hosted at the University of Montpellier.

##  V Appendices

###  V.1 Duschinsky matrix and shift vectors

Normal modes of vibration are computed at the minima of the ground electronic
state and the first excited electronic state. The transformation from the
normal modes of vibration at the minimum of the first electronic excited
state, $\bm{Q}^{\prime}$ to the normal modes of vibration at the minimum of
the electronic ground state $\bm{Q}^{\prime\prime}$ is as follows:

$\bm{Q}^{\prime}=\bm{J}\bm{Q}^{\prime\prime}+\bm{K}$ (10)

where $\bm{J}$ and $\bm{K}$ are the so-called Duschinsky matrix and shift
vector, respectively. For modes 81, 87, and 88, the Duschinsky matrix is:

$\begin{bmatrix}-0.118694&-0.000010&0.000187\\\ -0.050807&0.733369&0.662382\\\
-0.005732&-0.671423&0.740842\end{bmatrix}$ (11)

and the shift vector is

$\begin{bmatrix}-4.01414\\\ 7.25456\\\ 7.59322\end{bmatrix}.$ (12)

Note that the correlation between modes 81 and 87 or 88 is negligible (less
than 0.3% in any case). These small correlations suggest the validity for
neglecting the primary Duschinsky paramaters, $\gamma_{81,88}^{(k)}$ in the
QVC model. The correlation between 87 and 88 correspond to a
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ angle of the ground state normal
modes to obtain the excited state normal modes, which is consistent with the
results found for the orientation of the ellipsoids in the plane (87,88) of
the QVC model PESs.

## References

## References

  * (1) V. Balzani, S. Campagna, G. Denti, A. Juris, S. Serroni, and M. Venturi, Harvesting sunlight by artificial supramolecular antennae, Solar Energy Materials and Solar Cells, 1995, 38 (1) 159–173.
  * (2) S. L. Gilat, A. Adronov, and J. M. J. Fréchet, Light Harvesting and Energy Transfer in Novel Convergently Constructed Dendrimers, Angewandte Chemie International Edition, 1999, 38 (10) 1422–1427.
  * (3) A. Adronov and J. M. J. Fréchet, Light-harvesting dendrimers, Chemical Communications, 2000, 0 (18) 1701–1710.
  * (4) V. Balzani, P. Ceroni, M. Maestri, and V. Vicinelli, Light-harvesting dendrimers, Current Opinion in Chemical Biology, 2003, 7 (6) 657–665.
  * (5) Z. Xu and J. S. Moore, Design and synthesis of a convergent and directional molecular antenna, Acta Polymerica, 1994, 45 (2) 83–87.
  * (6) R. Kopelman, M. Shortreed, Z.-Y. Shi, W. Tan, Z. Xu _et al._ , Spectroscopic Evidence for Excitonic Localization in Fractal Antenna Supermolecules, Physical Review Letters, 1997, 78 (7) 1239–1242.
  * (7) S. F. Swallen, R. Kopelman, J. S. Moore, and C. Devadoss, Dendrimer photoantenna supermolecules: energetic funnels, exciton hopping and correlated excimer formation, Journal of Molecular Structure, 1999, 485-486 585–597.
  * (8) S. F. Swallen, Z. Zhu, J. S. Moore, and R. Kopelman, Correlated Excimer Formation and Molecular Rotational Dynamics in Phenylacetylene Dendrimers, The Journal of Physical Chemistry B, 2000, 104 (16) 3988–3995.
  * (9) D. Rana and G. Gangopadhyay, Steady-state spectral properties of dendrimer supermolecule as a light harvesting system, Chemical Physics Letters, 2001, 334 (4) 314–324.
  * (10) C. Wu, S. V. Malinin, S. Tretiak, and V. Y. Chernyak, Exciton scattering and localization in branched dendrimeric structures, Nature Physics, 2006, 2 (9) 631–635.
  * (11) J. L. Palma, E. Atas, L. Hardison, T. B. Marder, J. C. Collings _et al._ , Electronic Spectra of the Nanostar Dendrimer: Theory and Experiment, The Journal of Physical Chemistry C, 2010, 114 (48) 20702–20712.
  * (12) C. Devadoss, P. Bharathi, and J. S. Moore, Energy Transfer in Dendritic Macromolecules: Molecular Size Effects and the Role of an Energy Gradient, Journal of the American Chemical Society, 1996, 118 (40) 9635–9644.
  * (13) M. R. Shortreed, S. F. Swallen, Z.-Y. Shi, W. Tan, Z. Xu _et al._ , Directed Energy Transfer Funnels in Dendrimeric Antenna Supermolecules, The Journal of Physical Chemistry B, 1997, 101 (33) 6318–6322.
  * (14) V. D. Kleiman, J. S. Melinger, and D. McMorrow, Ultrafast Dynamics of Electronic Excitations in a Light-Harvesting Phenylacetylene Dendrimer, The Journal of Physical Chemistry B, 2001, 105 (24) 5595–5598.
  * (15) J. S. Melinger, Y. Pan, V. D. Kleiman, Z. Peng, B. L. Davis, D. McMorrow, and M. Lu, Optical and Photophysical Properties of Light-Harvesting Phenylacetylene Monodendrons Based on Unsymmetrical Branching, Journal of the American Chemical Society, 2002, 124 (40) 12002–12012.
  * (16) S. Fernandez-Alberti, V. D. Kleiman, S. Tretiak, and A. E. Roitberg, Nonadiabatic Molecular Dynamics Simulations of the Energy Transfer between Building Blocks in a Phenylene Ethynylene Dendrimer, The Journal of Physical Chemistry A, 2009, 113 (26) 7535–7542.
  * (17) S. Fernandez-Alberti, A. E. Roitberg, V. D. Kleiman, T. Nelson, and S. Tretiak, Shishiodoshi unidirectional energy transfer mechanism in phenylene ethynylene dendrimers, The Journal of Chemical Physics, 2012, 137 (22) 22A526.
  * (18) S. Fernandez-Alberti, D. V. Makhov, S. Tretiak, and D. V. Shalashilin, Non-adiabatic excited state molecular dynamics of phenylene ethynylene dendrimer using a multiconfigurational Ehrenfest approach, Physical Chemistry Chemical Physics, 2016, 18 (15) 10028–10040.
  * (19) T. Nelson, S. Fernandez-Alberti, A. E. Roitberg, and S. Tretiak, Electronic Delocalization, Vibrational Dynamics, and Energy Transfer in Organic Chromophores, The Journal of Physical Chemistry Letters, 2017, 8 (13) 3020–3031.
  * (20) A. L. Thompson, K. M. Gaab, J. Xu, C. J. Bardeen, and T. J. Martínez, Variable Electronic Coupling in Phenylacetylene Dendrimers: The Role of Förster, Dexter, and Charge-Transfer Interactions, The Journal of Physical Chemistry A, 2004, 108 (4) 671–682.
  * (21) E. K. L. Ho and B. Lasorne, Diabatic pseudofragmentation and nonadiabatic excitation-energy transfer in meta-substituted dendrimer building blocks, Computational and Theoretical Chemistry, 2019, 1156 25–36.
  * (22) K. M. Gaab, A. L. Thompson, J. Xu, T. J. Martínez, and C. J. Bardeen, Meta-Conjugation and Excited-State Coupling in Phenylacetylene Dendrimers, Journal of the American Chemical Society, 2003, 125 (31) 9288–9289.
  * (23) Q. Chu and Y. Pang, Vibronic structures in the electronic spectra of oligo(phenylene ethynylene): effect of m-phenylene to the optical properties of poly(m-phenylene ethynylene), Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy, 2004, 60 (7) 1459–1467.
  * (24) L. Cederbaum, W. Domcke, H. Köppel, and W. Von Niessen, Strong vibronic coupling effects in ionization spectra: The “mystery band” of butatriene, Chemical Physics, 1977, 26 (2) 169–177.
  * (25) H. Köppel, W. Domcke, and L. S. Cederbaum, Multimode Molecular Dynamics Beyond the Born-Oppenheimer Approximation, in Advances in Chemical Physics, 1984, 59–246.
  * (26) C. Cattarius, G. A. Worth, H.-D. Meyer, and L. S. Cederbaum, All mode dynamics at the conical intersection of an octa-atomic molecule: Multi-configuration time-dependent Hartree (MCTDH) investigation on the butatriene cation, The Journal of Chemical Physics, 2001, 115 (5) 2088–2100.
  * (27) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb _et al._ , Gaussian 16 Revision A.03, 2016.
  * (28) C. Adamo and D. Jacquemin, The calculations of excited-state properties with Time-Dependent Density Functional Theory, Chemical Society Reviews, 2013, 42 (3) 845–856.
  * (29) E. K.-L. Ho, T. Etienne, and B. Lasorne, Vibronic properties of para-polyphenylene ethynylenes: TD-DFT insights, The Journal of Chemical Physics, 2017, 146 (16) 164303.
  * (30) B. Gonon, A. Perveaux, F. Gatti, D. Lauvergnat, and B. Lasorne, On the applicability of a wavefunction-free, energy-based procedure for generating first-order non-adiabatic couplings around conical intersections, The Journal of Chemical Physics, 2017, 147 (11) 114114.
  * (31) M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets, Physics Reports, 2000, 324 (1) 1–105.
  * (32) G. Worth, K. Giri, G. Richings, I. Burghardt, M. Beck, A. Jäckle, and H. Meyer, The quantics package, University of Birmingham: Birmingham, UK, 2015.
  * (33) A. Jaouadi, J. Galiana, E. Mangaud, B. Lasorne, O. Atabek, and M. Desouter-Lecomte, Laser-controlled electronic symmetry breaking in a phenylene ethynylene dimer: Simulation by the hierarchical equations of motion and optimal control, Phys. Rev. A, 2022, 106 043121.