Franck–Condon simulation for unraveling vibronic
origin in solvent enhanced absorption and
fluorescence spectra of rubrene†
Ying Hu,‡
ab
Chen-Wen Wang,‡
b
Chaoyuan Zhu,
*
ab
Fenglong Gu
*
a
and Sheng-Hsien Lin
b
Quantum chemistry calculations at the level of (TD)-DFT plus PCM solvent models are employed for
analyzing potential energy surfaces and as a result two local minima with D
2
, two local minima with C
2H
,
and one second-order transition state with D
2H
group symmetry are found in both ground S
0
and
excited-state S
1
potential energy surfaces. Simulated vibronic coupling distributions indicate that only
second-order transition states with D
2H
group symmetry are responsible for observed absorption and
fluorescence spectra of rubrene and vibrational normal-motions related with atoms on the aromatic
backbone are active for vibronic spectra. The Stokes shift 1120 cm
1
(820 cm
1
) and vibronic-band peak
positions in both absorption and fluorescence spectra in non-polar benzene (polar cyclohexane) solvent
are well reproduced within the conventional Franck–Condon simulation. By adding damped oscillator
correction to Franck–Condon simulation, solvent enhanced vibronic-band intensities and shapes are well
reproduced. Four (three) normal modes with vibration frequency around 1550 cm
1
(1350 cm
1
) related
to ring wagging plus CC stretching and CH bend motions on the backbone are actually interpreted for
solvent enhanced absorption (fluorescence) spectra of rubrene in benzene and cyclohexane solutions.
1. Introduction
Franck–Condon (FC) factors play a very important role for
theoretical modeling and interpreting experimental observa-
tions for irradiative processes like electronic spectra of VUV
absorption and uorescence as well as the nonradiative
processes like electron and energy transfer.
1–10
For large poly-
atomic molecules, multidimensional Franck–Condon over-
lapping integrals are conventionally decomposed into a product
of one-dimensional Franck–Condon overlap integrals in terms
of normal mode coordinates under harmonic oscillator
approximation. Harmonic oscillator approximation provides
the leading-term contribution to FC factors, and further higher
order corrections of distorted oscillator approximation and
mode-mixing including Duschinsky effect and anharmonic
oscillators have been introduced to improve FC factors and its
various applications.
11–17
Most of the electronic molecular
spectra and vibrational relaxation processes are experimentally
carried out in solution phases. Therefore, theoretical modeling
can deal with solutes and solvent molecules separately and FC
factors including electronic and nuclear parts are treated
separately as well. The electronic parts involve calculating
electronic transition dipole momentum and nonadiabatic
coupling matrix under solvent environment: the most popular
theoretical approach is represented by continuum solvation
models, especially the polarizable continuum model (PCM)
combined with time-dependent density functional theory
(TDDFT). The (TD)-DFT plus PCM approach provides an accu-
rate simulation of static electronic potential energy surfaces of
both ground and excited states for molecules in solution. The
nuclear parts represent dynamic motion and its FC factors have
been investigated based on various corrections to the leading-
term FC factors.
18–36
By introducing damped harmonic oscil-
lator approximation, Zhu and co-authors
36
formulated the
empirical correction to each normal mode differently in solvent
environment by converting the mass-weighted unperturbed gas-
phase Hessian matrix to solution-phase perturbed Hessian
matrix. This method is based on idea that each normal-mode
motion can be affected differently in solution, however each
local-mode motion can be affected in quite similar way for
highly symmetric molecules like polycyclic aromatic hydro-
carbon molecules in which CH-stretch local modes, for
example, are considered as in the same solvent environment.
This method was previously demonstrated for interpreting
solvent enhanced absorption and uorescence spectra of
a
Key Laboratory of Theoretical Chemistry of Environment, Ministry of Education,
School of Chemistry & Environment, South China Normal University, Guangzhou
510006, P. R. China. E-mail: gu@scnu.edu.cn
b
Department of Applied Chemistry, Institute of Molecular Science and Center for
Interdisciplinary Molecular Science, National Chiao-Tung University, Hsinchu
30050, Taiwan. E-mail: cyzhu@mail.nctu.edu.tw
† Electronic supplementary information (ESI) available. See DOI:
10.1039/c7ra00417f
‡ Y. Hu and C.-W. Chen contributed equally to this work.
Cite this: RSC Adv.,2017,7,12407
Received 11th January 2017
Accepted 14th February 2017
DOI: 10.1039/c7ra00417f
rsc.li/rsc-advances
This journal is © The Royal Society of Chemistry 2017 RSC Adv.,2017,7,12407–12418 | 12407
RSC Advances
PAPER
Open Access Article. Published on 21 February 2017. Downloaded on 8/26/2022 12:18:09 PM.
This article is licensed under a
Creative Commons Attribution 3.0 Unported Licence.
View Article Online
View Journal
| View Issue
perylene in benzene solution by applying a single empirical
parameter to CH-stretch local modes.
36
One of the purposes in the present study is to continue
damped harmonic oscillator FC simulation for much large size
of polycyclic aromatic hydrocarbon like rubrene. The other is to
interpret experimentally observed vibronic spectra in terms of
potential energy surfaces of ground and excited states, transi-
tion patterns, molecular structure, solvent dependence and so
on. For instance, it is traditionally considered that vibronic
spectra simulated by FC factors are based on local minima
between ground- and excited-state potential energy surfaces,
and this is usually true for rigid molecules. However, for exible
molecules, potential energy surfaces are very complicated with
many local minima closed to one another, in this case, (high-
order) transition states that connect (more than) two local
minima can play crucial role in FC simulation in order to
interpret observed vibronic spectra.
Rubrene molecules have been attracted more attention in the
eld of organic semiconductor because of its high mobility at room
temperature and optical properties, and rubrene can be used as
alaserdyeandauorescentdopantinorganiclightemitting
devices, eld effect transistors and solar cell.
37–46
In the present
study we focus on the optical properties of monomer of rubrene in
non-polar benzene and polar cyclohexane solvents. Rubrene
molecule is constructed with four six-member rings as the aromatic
backbone connected by four phenyl rings perpendicular to the
backbone. The backbone is the uorescentcenterthatcanhave
planar or twisted geometry, and its optical spectra, however, are not
well understood in terms of structure of the backbone geometry as
well as phenyl groups. Moreover, it is well-know that rubrene has
the large Stokes shi, but what is dependence in terms of solvent
environment. We investigate how solvent enhanced absorption
and uorescence spectra vary with respect to solvents by damped
FC simulation and identify physical origin of vibronic spectra in
terms of ground- and excited-sate potential energy surfaces.
The rest of this paper is organized as follows: Section 2 briey
introduces damped harmonic oscillator FC factors and discuss
(TD)-DFT+ PCM method for calculating potential energy
surfaces of ground and excited states for rubrene molecule.
Section 3 discusses how to reproduce experimentally observed
absorption and uorescence spectra in benzene and cyclohexane
solution, and to interpret vibronic spectra in terms of molecular
structure, molecular orbitals, vibrational normal modes and
solvent environment. Section 4 presents concluding remarks.
2. Damped Franck–Condon factors
Within displaced harmonic oscillator approximation, absorp-
tion and uorescence coefficients can be analytically derived as
follows:
1,42
haðuÞifu
ð
N
N
dte
itðu
ba
uÞ
D
ba
2
t
2
4
g
ba
jtj
exp
"
X
j
S
j
2
n
j
þ 1
n
j
þ 1
e
itu
j
n
j
e
itu
j
#
(1)
and
h
IðuÞifu
3
ð
N
N
dte
itðju
ba
j
uÞ
D
ba
2
t
2
4
g
ba
j
t
j
exp
"
X
j
S
j
2
n
j
þ 1
n
j
þ 1
e
itu
j
n
j
e
itu
j
#
(2)
where u
ba
represents electronically adiabatic excitation energy
between electronic states b and a.
n
j
¼ðe
h
-
u
j
=k
B
T
1Þ
1
is the
average phonon distribution, g
ba
is the homogeneous broad-
ening parameter and D
ba
is inhomogeneous broadening
parameter that reects static interaction between electronic
states of molecule and solvent, and u
j
is the j-th normal mode
vibrational frequency. The dimensionless Huang–Rhys factor S
j
corresponding to the j-th normal mode is dened as
S
j
¼
1
2ħ
u
j
d
j
2
(3)
where displacement d
j
is given by
d
j
¼ Q
0
j
Q
j
¼
X
n
L
jn
q
0
n
q
n
(4)
The q
0
n
and q
n
in eqn (4) are the mass-weighted Cartesian
coordinates at equilibrium geometries (or transition states) of
the electronic excited and ground states, respectively. Trans-
formation matrix L in eqn (4) can be computed with frequency
analysis using G09 programs, for instance. By applying damped
diatomic harmonic oscillator to local-mode force constant of
gas-phase Hessian matrix H
ij
, we formulate damped FC factors
as transferring the mass-weighted gas-phase unperturbed
Hessian matrix to solvent perturbed Hessian matrix,
36
H
0
ia;jb
¼
k
0
ia;jb
ffiffiffiffiffiffiffiffiffiffi
m
i
m
j
p
/
k
0
ia;jb
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
ðm
i
z
i
Þ
m
j
z
j
q
¼
k
ia;jb
ffiffiffiffiffiffiffiffiffiffi
m
i
m
j
p
¼ H
ia;jb
; (5)
where m
i
and m
j
are atomic mass, and a and b are Cartesian
component (x, y, z) of the position of nucleus i and j. Unper-
turbed local force constant k
0
ia;jb
is converted to perturbed force
constant k
ia;jb
¼ k
0
ia;jb
=
ffiffiffiffiffiffiffiffi
z
i
z
j
p
. The dimensionless scaling
parameter z
j
in eqn (5) is treated as an empirical parameter that
can be considered as a local dynamic interaction of atom j of
solute molecule with solvent molecules. By diagonalizing per-
turbed Hessian matrix in eqn (5), we obtain normal mode
frequencies and transformation matrix (that transfers the mass-
weighted Cartesian coordinates to normal mode coordinates)
corresponding to solute molecule under the solvent environ-
ment. It should be emphasized that solvent effect by damped
oscillator can directly modify the leading-term of FC factors, so
that it is usually larger than the other effects such as distorted
oscillator, Duschinsky effect and anharmonic effect which are
in higher-order perturbation terms.
We have employed the (time-dependent) density functional
theory with functional (TD)-B3LYP
47
(i.e., B3LYP20), (TD)
B3LYP35, (TD)-B3LYP50 (i.e., BHandHLYP),
48
and HF–CIS (i.e.,
B3LYP100) plus PCM condition in benzene and cyclohexane
solvents for calculating potential energy surfaces of ground and
the rst excited states. The numbers following aer B3LYP
12408 | RSC Adv.,2017,7,12407–12418 This journal is © The Royal Society of Chemistry 2017
RSC Advances Paper
Open Access Article. Published on 21 February 2017. Downloaded on 8/26/2022 12:18:09 PM.
This article is licensed under a
Creative Commons Attribution 3.0 Unported Licence.
View Article Online
stand for different hybrid exchange–correlation functionals
containing 20% (B3LYP), 35%, 50%, and 100% of exact Hartree–
Fock exchange in the density functional theory. The basis set is
chosen as 6-31G throughout all calculations and all calculations
are carried out using GAUSSIAN 09 program package.
49
Calcu-
lations from four functionals given above lead to the same
conclusion in which there are two local minima with D
2
, two
local minima with C
2H
, and one second-order transition state
with D
2H
group symmetry in both ground- and excited-state
potential energy surfaces. This conclusion does hold for PCM
condition in both non-polar benzene and polar cyclohexane
solvents although electronic structures can vary slightly. From
preliminary vibronic spectral calculation, we found that (TD)-
B3LYP50 presents the best agreement with experimental
results among four functionals, and this is the same as results
obtained in ref. 42. However, there is important difference in
comparison with ref. 42 in which there are only two local
minima, one is in D
2
and the other is in D
2H
symmetry (note for
this point, we show it is second-order transition state). The
Huang–Rhys factors calculated from our D
2H
structure are
similar to those from their D
2H
symmetry, but the Huang–Rhys
factors calculated from our D
2
structure are very different to
those from their D
2
symmetry. Actually, the second-order tran-
sition point at D
2H
connects two local minima at D
2
as well as
two local minima at C
2H
, and we will discuss this in detail later.
In the following discussions, we only report calculations from
(TD)-BHandHLYP (B3LYP50) as it performs best for vibronic
spectral simulations.
3. Results and analysis
Rubrene molecule plays as a prototype of large-size polycyclic
aromatic hydrocarbon molecules for a better understanding of
molecular structure, spectra, dynamics, solvent effect, and so
on. It's highly symmetrical electronic structures (see atomic
numbering in Fig. 1) in which all hydrogen atoms are approxi-
mately subjected to the same environment interacting with
solvent molecules and this greatly simplies empirical search of
the scaling parameters in damped FC simulation. Geometry of
rubrene molecule is composed with four linearly connected six-
member rings as aromatic backbone connected by four phenyl
rings perpendicular to the backbone. The present (TD)-DFT
calculations reveal that there are two local minima at C
2H
group symmetry where C
2H
(1) minimum is mirror image of
C
2H
(2) minimum as shown in Fig. 2(a) and (b), there is one
second-order transition state at D
2H
group symmetry as shown
in Fig. 2(c), and there are two local minima at D
2
group
symmetry where D
2
(1) minimum is mirror image of D
2
(2)
minimum as shown in Fig. 2(d) and (e). These kind of electronic
conformers displayed in Fig. 2 show similar structures for
ground and excited states in both benzene and cyclohexane
solvents. At D
2H
transition state, there are actually two normal
modes with imaginary frequencies 58.62 cm
1
and 54.49
cm
1
(58.06 cm
1
and 53.21 cm
1
) in Au and B3G symme-
tries, respectively for the ground (excited) state in benzene
solvent. Along Au and B3G normal-mode coordinates, we draw
two-dimensional potential energy surface diagrams in which
second-order D
2H
transition state connects two D
2
minima via
Au and two C
2H
minima via B3G as shown in Fig. 3(a) for ground
state and Fig. 3(b) for excited state in benzene solvent envi-
ronment. Similar potential energy surface diagrams in cyclo-
hexane solvent environment are also shown in Fig. 4(a) and (b),
where two imaginary frequencies are 56.54 cm
1
and 52.85
cm
1
(55.70 cm
1
and 51.29 cm
1
) with Au and B3G mode
symmetries, respectively for the ground (excited) state. Differ-
ence between ground and excited state potential energy surfaces
is not large, and this conrms that the present displaced
harmonic FC simulation should be good approximation to
reproduce experimental vibronic spectra. The adiabatic poten-
tial energy gap between D
2H
-transition state and D
2
-local
minima (C
2H
-local minima) is 0.31 eV (0.12 eV) as shown in
Fig. 3(a) and 4(a) for S
0
state and is 0.51 eV (0.20 eV) as shown in
Fig. 3(b) and 4(b) for S
1
state, respectively in both benzene and
cyclohexane solvents. Potential energy surfaces around D
2H
-
transition state in terms of two imaginary normal-mode coor-
dinates show little dependence on solvent.
3.1. Electronic structures, potential energy surfaces and the
Huang–Rhys factors
The present (TD)-DFT+ PCM calculations reveal that the back-
bone formed by the four six-member rings is planar structure at
local minima C
2H
(1), C
2H
(2) and the second-order transition
state D
2H
, but the backbone is twisted at local minima D
2
(1) and
D
2
(2) for both S
0
and S
1
states. Since there are existing D
2
, C
2H
and D
2H
group symmetries at ve optimized geometries for both
Fig. 1 Atomic numbering (optimized geometries under D
2H
group
symmetry by BHandHLYP/6-31G plus PCM model. The yellow and
blue indicate carbon and hydrogen atoms respectively).
This journal is © The Royal Society of Chemistry 2017 RSC Adv.,2017,7,12407–12418 | 12409
Paper RSC Advances
Open Access Article. Published on 21 February 2017. Downloaded on 8/26/2022 12:18:09 PM.
This article is licensed under a
Creative Commons Attribution 3.0 Unported Licence.
View Article Online
S
0
and S
1
states, the internal coordinates of selected bond
lengths, bond angles and dihedral angles which well represent
geometries of electronic structures contributing to vibronic
spectra are given in Tables 1 and 2. Tables 1 and 2 actually show
difference of bond lengths, bond angles and dihedral angles
between S
0
and S
1
states with the same group symmetry in
Fig. 2 Optimized geometries by BHandHLYP/6-31G plus PCM model in benzene solvent for ground-state S
0
. Local minima (a) C
2H
(1) and (b)
C
2H
(2), (c) second-order transition state D
2H
, and local minima (d) D
2
(1) and (e) D
2
(2).
12410
| RSC Adv.,2017,7,12407–12418 This journal is © The Royal Society of Chemistry 2017
RSC Advances Paper
Open Access Article. Published on 21 February 2017. Downloaded on 8/26/2022 12:18:09 PM.
This article is licensed under a
Creative Commons Attribution 3.0 Unported Licence.
View Article Online
benzene and cyclohexane solvents, respectively, and the corre-
sponding absolute values of those internal coordinates are
shown in Tables S1 and S2† (Full Cartesian coordinates for
optimized geometries shown in Table S3) (ESI†). As is well-
known, the magnitude of Huang–Rhys factor dened in eqn
(3) depends on the displacement of the normal mode between
local minimum on ground state and corresponding local
minimum on excited state, while the displacement of normal
mode is converted from differences of internal coordinates
between local minima of the ground and excited states. Bond-
length differences vary with respect change of group
symmetry from D
2H
to C
2H
and from D
2H
to D
2
within small
magnitude (the largest value is about 0.05
˚
A for C9C10 bond in
benzene solvent) as shown in Tables 1 and 2, and this can lead
the change of Huang–Rhys factor about 0.4. However, bond-
angle differences can vary within as large as 6
, and this can
lead to the change of Huang–Rhys factor about 4. The largest
change happens for dihedral angles as shown in Tables 1 and 2
and it can go as large as ten degree, but these changes affect
change of Huang–Rhys factor in complicated way and it is hard
to predict the magnitude change before carrying out real
calculations. This fact was actually conrmed in which the
Huang–Rhys factors of the some vibrational modes can go as
large as from 5 to 50 due to big change in dihedral angles for the
two isomeric compounds.
50
Dihedral angle differences between
S
0
and S
1
states at D
2H
group symmetry are all equal to zero as
Fig. 3 Potential energy surface along two imaginary normal-mode coordinates (Au and B3G symmetries) starting from D
2H
transition state point
which connects two C
2H
minima via B3G-mode direction, and two D
2
minima via Au-mode direction. (a) Ground state and (b) excited state in
benzene solvent.
Fig. 4 Potential energy surface along two imaginary normal-mode coordinates (Au and B3G symmetries) starting from D
2H
transition state point
which connects two C
2H
minima via B3G-mode direction, and two D
2
minima via Au-mode direction. (a) Ground state and (b) excited state in
cyclohexane solvent.
This journal is © The Royal Society of Chemistry 2017 RSC Adv.,2017,7,12407–12418 | 12411
Paper RSC Advances
Open Access Article. Published on 21 February 2017. Downloaded on 8/26/2022 12:18:09 PM.
This article is licensed under a
Creative Commons Attribution 3.0 Unported Licence.
View Article Online