Subscriber access provided by McMaster University Library
Article
Photoinduced Charge Transfer versus Fragmentation Pathways in Lanthanum Cyclopentadienyl Complexes Yulun Han, Qingguo Meng, Bakhtiyor Rasulev, P. Stanley May, Mary T. Berry, and Dmitri S. Kilin J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00050 • Publication Date (Web): 26 Jun 2017 Downloaded from http://pubs.acs.org on June 29, 2017
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
1
Photoinduced
2
Lanthanum Cyclopentadienyl Complexes
3
Yulun Han,†,‡ Qingguo Meng,§ Bakhtiyor Rasulev,∥ P. Stanley May,† Mary T. Berry,† and Dmitri S. Kilin*,†,‡
4
†
5 6 7 8 9
Transfer
versus
Fragmentation
Pathways
in
Department of Chemistry, University of South Dakota, Vermillion, South Dakota 57069, United States
‡ §
Charge
Department of Chemistry and Biochemistry, North Dakota State University, Fargo, North Dakota 58108, United States
Shenyang Institute of Automation, Guangzhou, Chinese Academy of Sciences, Guangzhou 511458, China
∥Center
for Computationally Assisted Science and Technology, North Dakota State University, Fargo, North Dakota 58102,
United States
10 11
Abstract:
12
This study compares two competing pathways of photoexcitations in gas−phase metal−organic
13
complexes: first, a sequence of phonon−assisted electronic transitions leading to dissipation of the energy
14
of photoexcitations and second, a sequence of light-driven electronic transitions leading to photolysis.
15
Phonon−assisted charge carrier dynamics is investigated by combination of the density matrix formalism
16
and on−the−fly nonadiabatic couplings. Light-driven fragmentation is modeled by a time−dependent
17
excited−state molecular dynamics (TDESMD) algorithm based on Rabi theory and principles similar to
18
the trajectory surface hopping approximation. Numerical results indicate that, under high intensity of the
19
laser field, light-driven electronic transitions are more probable than phonon−assisted ones. The
20
formation of multiple products is observed in TDESMD trajectories. Simulated mass spectra are extracted
21
from TDESMD simulations and compared to experimental photoionization time−of−flight (PI−TOF)
22
mass spectra. It is found that several features in the experimental mass spectra are reproduced by the
23
simulations.
24
1 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
25
1. INTRODUCTION
26
Lanthanide oxide thin films have attracted much attention due to their wide range of
27
applications. The high dielectric constant and good interfacial properties make them potential
28
alternatives to SiO2 in metal−oxide−semiconductor field−effect transistors (MOSFETs).1-3
29
Lanthanide oxide thin films have applications as waveguides in optics, due to their high
30
refractive index, large band gap, and low phonon energy.4, 5 Other applications include protective
31
and corrosion−resistive coatings,6, 7 and phosphors.8 Metal−organic chemical vapor deposition
32
(MOCVD) and atomic layer deposition (ALD) techniques have been employed to grow
33
lanthanide oxide thin films using a number of precursors including lanthanide cyclopentadienyl
34
complexes.9-13 The high volatility of lanthanide cyclopentadienyl complexes makes them
35
potential candidates for fabrication of lanthanide oxide thin films by laser−assisted chemical
36
vapor deposition (LCVD), where low film deposition temperatures are allowed. The
37
understanding of the photofragmentation mechanisms of gas−phase lanthanide cyclopentadienyl
38
complexes is essential for LCVD applications.
39
Luminescence spectroscopy and photoionization time−of−flight (PI−TOF) mass
40
spectrometry have been used to study the photofragmentation of gas−phase metal−containing
41
molecules.14-25 The corresponding computational analysis of these processes, however, is
42
complicated and challenging. The laser field can induce both excitation and de−excitation
43
processes between the electronic states of the target molecule. Nonadiabatic transitions couple
44
nuclear dynamics with the electronic states, resulting in the breakdown of the
45
Born−Oppenheimer approximation.26 The nonradiative de−excitation processes are occurring in
46
broad range of molecules, solids, and nanostructures, and often referred by different terminology,
47
such as e.g. hot charge carrier (electron or hole) relaxation in quantum dot spectroscopy. Various 2 ACS Paragon Plus Environment
Page 2 of 52
Page 3 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
48
first−principle approaches beyond Born−Oppenheimer approximations have been developed to
49
describe such hot carrier relaxation dynamics in different molecules,27-29 solids,30,
50
nanostructures.32-34 To address the electron−nuclear correlations, trajectory surface hopping
51
methods have been widely used, where the dynamics of a molecular system is governed by
52
potential energy surfaces.35, 36 The dynamics of a molecular system can be interrupted by hops
53
between different potential energy surfaces36 or nonadiabatic events,37 as reflected in recent
54
implementations of nonadiabatic molecular−dynamics approaches targeting simultaneous
55
propagation of both electronic and nuclear degrees of freedom.38-41 Liouville−von Neumann
56
equations of motion or multilevel Redfield theory42-50 can be used to propagate the density
57
matrix of the electronic degrees of freedom.51-53
58
In
a
recent
study,
the
photofragmentaion
mechanisms
of
31
and
gas−phase
59
tris(η5−cyclopentadienyl)lanthanum, La(Cp)3 were studied by photoionization time−of−flight
60
(PI−TOF) mass spectrometry.23,
61
(TDESMD) based on Rabi theory and principles similar to the trajectory surface hopping
62
approximation were used to describe the photofragmentation reaction.23,
63
eikonal approach by Micha and coworkers finds eikonal wave function and treats quantum phase
64
of vibrations.55 Nonadiabatic excited-state molecular dynamics (NA-ESMD) by Tretiak et al.
65
provides a careful treatment to excited state potential surface and parametrized configuration
66
interaction for the electronic structure. However, they limit the consideration to dynamics created
67
by a single photon excitation event.56 Levine and coworkers uses multi-reference methods but
68
also for a single excitation only.57-59 Bubin and Varga considers continuous excitations by light
69
of high intensity. It requires recalculation of electronic structure with very short time step.60
54
The time−dependent excited−state molecular dynamics
3 ACS Paragon Plus Environment
54
A self-consistent
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 52
70
The TDESMD method is an extension of MD. The difference is that TDESMD involves
71
a non-ground state electronic configuration. The electrons are hopping between ground state and
72
excited state, as dictated by Rabi oscillations. The uniqueness of the TDESMD approach is that it
73
goes beyond single excitation event on one hand and offers approximations reducing the
74
numerical cost to the acceptable level. The implementation of TDESMD uses two major
75
principles: combination of Rabi theory and surface hopping approximation. The sequence of
76
recent papers from our group provides views of the photochemistry modeling from different
77
aspects, tests different regimes, and applies the TDESMD methodology to model molecules of
78
drastically different nature, as explained in what follows. The original paper provides general
79
idea in a very simplified implementation.54 A detailed comparison of reaction mechanisms
80
developed from analysis of TDESMD calculations and analysis of experimental data is provided
81
in reference
82
comparison of alternative limits of excited state MD is analyzed in reference
83
TDESMD to catch qualitatively different pathways of the fragmentation reactions is justified in
84
reference
85
explosive organics at the CNT substrate is established in reference
86
neglected the phonon−assisted electronic transitions and assumed that light-driven electronic
87
transitions play a dominant role in determining the photofragmentation reaction. Both pathways
88
are triggered by a photoexcitation. In a phonon−induced pathway energy of photoexcitation is
89
converted into low amplitude kinetic energy of nuclear motion via a cascade of electronic
90
transitions with small amount of energy transferred from electronic to nuclear degrees of
91
freedom, during each elementary event. In a light-driven pathway energy of photoexcitation is
92
converted into kinetic energy of nuclear motion via walking along excited and ground state
25
23
. A conversion of the TDESMD trajectory into theoretical mass-spectrum and 24
. The ability of
. A role of charge transfer events in the TDESMD trajectory in application to
4 ACS Paragon Plus Environment
61
. In earlier studies, we
Page 5 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
93
potential energy surfaces leading to overcoming of potential energy surface barriers and
94
fragmentation events.
95
Here, we compare two competing pathways of photoexcitations in gas−phase La(Cp)3: a
96
sequence of phonon−induced electronic transitions leading to dissipation of energy of
97
photoexcitations and a sequence of light-driven electronic transitions leading to photolysis. The
98
former pathway is investigated by a density matrix formalism with the combination of Redfield
99
theory and on−the−fly coupling of electron−to−nuclear motion by the molecular dynamics
100
trajectory. The latter pathway is explored again by applying TDESMD algorithm. Numerical
101
results indicate that, under high intensity of the laser field, light-driven electronic transitions are
102
more probable than phonon−induced electronic transitions. The dynamical formation of multiple
103
products is observed in TDESMD trajectories. Simulated mass spectra are extracted from
104
TDESMD simulations by using both La(Cp)3 and the intermediate La(Cp) as a starting point. It
105
is found that several features in the experimental photoionization time−of−flight (PI−TOF) mass
106
spectra are reproduced by the simulations.
107
The paper is organized as follows. The Experimental and Theoretical methods section
108
consists of seven subsections. Experimental details are summarized in subsection 2.1.
109
Computational approaches for phonon-induced pathways, dynamics of electrons, and light-
110
driven fragmentation are shown in subsection 2.2, subsection 2.3, and subsection 2.4,
111
respectively.
112
fragmentation is overviewed in subsection 2.5. Computed observables characterizing light-driven
113
fragmentation are presented in subsection 2.6. Subsection 2.7 describes computational details.
114
The Results section is divided into subsections on electronic structure (subsection 3.1), charge
115
carrier dynamics (subsection 3.2), and fragmentation (subsection 3.3).
Practical
implementation
of
computational
5 ACS Paragon Plus Environment
approaches
for
light-driven
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
116
117
2. EXPERIMENTAL AND THEORETICAL METHODS
118
2.1. PI−TOF Mass Spectrometry
119
La(Cp)3 precursors were purchased from Sigma−Aldrich and used as received. The setup
120
of PI−TOF experiments is described in detail in previous publications.16, 17, 23 La(Cp)3 precursors
121
were sublimed, without carrier gas, into the PI−TOF chamber at ~ 450 K. The pressure in
122
sample holder and PT−TOF chamber was approximately 10 mTorr and 10−4 mTorr, respectively.
123
A photofragmentation wavelength of 266 nm was provided by the fourth harmonic of a Nd:YAG
124
laser (Continuum Surelite II) at 70 mJ/pulse, 10 Hz repetition rate, and a pulse width of
125
approximately 6 ns FWHM. The excitation pulse is focused in the PI volume to produce an
126
approximate average power density of ~ 109 W/cm2, assuming a Gaussian pulse distribution.
127
128
2.2 Computational Approaches for Phonon-induced Pathways
129
The excited−state dynamics of electronic and nuclear degrees of freedom are generally
130
coupled and, therefore, inseparable. In order to save numerical resources, we explore our models
131
employing two methods. The first method focuses on dynamics of electrons while the second
132
method focuses on nuclear motion. The optically driven transitions are included by an
133
approximation to Rabi oscillations. Redfield dissipation terms originate from careful averaging
134
of nonadiabatic coupling computed by ab initio. The electronic transitions induced by nuclear
135
motion are treated by Redfield approach.
136 6 ACS Paragon Plus Environment
Page 6 of 52
Page 7 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
137
Journal of Chemical Theory and Computation
2.3 Computational Approaches for Dynamics of Electrons For practical calculations, the density operator must be represented in a specific basis.62,
138 139 140
63
Here we choose the basis of Kohn−Sham (KS) orbitals with a one−particle density matrix
playing the role of expansion coefficients, defined as
, ∑ ∙ , ∙ ∗ , .
141 142 143 144
145
146 147 148 149 150 151 152 153 154 155
156
(1)
The electronic degrees of freedom obey the equation of motion for the density operator ρ, which can be cast in terms of the Liouville−von Neumann super−operator and Redfield super−operator ,
$ , ρ% &'()*, ρ ρ "F ! '+
'-.
.
(2)
Liouville superoperator acting on arbitrary operator /0 is defined as follows and is equivalent to
commutator of argument operator /0 with Fock operator (Hamiltonian) 1$ , /0 ! "1$ , /0%
! 21$ /0 /01$ 3. In application of Liouville operator to density operator (matrix), and by casting
operators in form of matrix elements in a specific basis this definition reads 9 ∙ : , where ! ∑4214 4 4 14 3 . The Fock operator reads 1$ 1$ 5$ 67 8
9 stands for stands for ground state Fock matrix, 5$ 67 represents nonadiabatic couplings, 8
transition dipole operator, : : ∙ cos>
[email protected] is the electric field of the incident light.64, 65 The second term on the right side of equation (2) describes the influence of nuclear motion on
the electronic degrees of freedom. Elements of the Redfield tensor, 4?C , provide the probability of dissipative electronic transition,
&
'()* '+
,
'..
∑?C 4?C ?C 7
ACS Paragon Plus Environment
(3)
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
157
where &
158
computed along molecular dynamics trajectory. Nonadiabatic couplings between electronic and
159
nuclear degrees of freedom originate from a matrix element of kinetic energy of nuclei in the
160 161
'()* '+
,
Page 8 of 52
'..
basis of electronic wavefunctions. Numerical values of the nonadiabatic couplings, 567 , can be found according to the on−the−fly procedure, in finite differences,66-68
' 567 E! < , G'+G , >.
162
163 164
167
averaging along the duration of the trajectory, J,
84? I L MO 5 I54? N.
P ±,
172 173
(5)
A Fourier transform of coupling autocorrelation function, 84? I, provides components R P4? L MO 84? IS TUVWX NI K
L
T P4? L MO 84? IS TU)*X NI K
169
171
L
K
168
170
(4)
The nonadiabatic couplings are then processed with the autocorrelation function, 84? I,
165 166
are electronic dissipative transitions facilitated by thermal fluctuations of ions
L
(6) (7)
These components provide Redfield tensor, 4? , which controls dissipative dynamics of
density matrix,
R R T T 4? P4? P4? Y? ∑C P4? Y4 ∑C P4?
where Y? and Y4 are Kronecker Y symbols.
8 ACS Paragon Plus Environment
(8)
Page 9 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
174
Solution of equation (2) provides several observables including the charge density
175
distribution, rate of energy distribution, and rate of charge transfer. The nonequilibrium
176
distribution of charge as a function of time and energy reads, Z@,[ \, ∑
@,[
177
Y\ \
(9)
178
where (a,b) corresponding to the initial photoexcitation from state a to b. The change of
179
population from the equilibrium distribution is, ∆Z@,[ \, Z@,[ \, ZA^ \,
180 181 182
where ∆Z > 0 the equation describes the population gain and ∆Z < 0 the population loss. The expectation value of a hot electron can be expressed as follows, 〈∆\A 〉 ∑abc \
183
〈:A 〉
184
185 186
187
(10)
〈∆de 〉+T〈∆de 〉f
〈∆de 〉OT〈∆de 〉f
(11)
.
(12)
Assuming a single exponential fit of the energy dissipation, one can calculate the rate of relaxation for electrons as, gA hI A iTK MO 〈:A 〉N . f
TK
(13)
188
189
2.4 Computational Approaches for Light-driven Fragmentation
190
TDESMD calculations can be performed to simulate the photo-dissociation in the mid
191
laser intensity, where photons sequentially excite the molecule and induce sequence of
192
absorption and stimulated emission events. At lower intensities, these events are rare and the
193
modeling is not needed. At higher intensities, the adequate modeling should include broader 9 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 52
194
range of pathways, highly nonlinear in terms of the incident light intensity. The procedure of
195
TDESMD is described in detail in previous work.23-25, 54, 61
196
In TDESMD simulation, the photo-fragmentation process is explored by periodic
197
excitations and de−excitations of the model due to the light−matter interaction term in total
198
energy and in the Fock matrix. Numerical propagation of equation of motion using
199
time−dependent Hamiltonian Fock operators is challenging, since one has to use a time step
200 201
much smaller than the laser field period, i.e. ∆ ≪ 2l/>
[email protected] , of the order of 10−18 s to track
response of electrons to rapid changes in electric field. Use of small time steps limits the duration
202
of the explored trajectory. This challenge is bypassed by rotating wave approximation (RWA) in
203
the interaction picture. Application of this approximation allows to hope for acceptable precision
204
with use of substantially longer time steps, thus exploring trajectory of longer duration. RWA is
205 206 207 208
∙ : in the basis of based on representing matrix elements of light−matter interaction 8 adiabatic states |Eo, |po, with energies \ , \ , in the interaction picture.69 In this picture, transition
0S U)*+ , 8 dipole matrix follows time evolution with phase accumulation factors 8 where !r \ \ . Such a matrix element is written as, ∙ \ \ O ∙ s 8
abc
s
U)* + w+ S U)*+ S Tw+ 8 ySyzy t|Eoup| v 8 Sy{ xyyyzyyy{ xy y |}~+ABTB}+@+
|}TB}+@+
TU)* + w+ TU)* + Tw+ yy yy |pouE| v8 Syzy Sy{ 8 Syzy xy yy xy ySyy{
. ∗
|}~+ABTB}+@+
∗
|}TB}+@+
209
(14)
210
The rapidly oscillating terms S TU)*+ ∙ S Tw+ and S U)*+ ∙ S w+ are dropped due to the rotating
211
wave approximation.70 The slowly rotating terms S U)*+ ∙ S Tw+ and S TU)*+ ∙ S w+ can be 10
ACS Paragon Plus Environment
Page 11 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
212 213
Journal of Chemical Theory and Computation
approximated as constant, and, therefore, time−independent, in the resonant limit r ≈ >
[email protected] . In this case, the light−matter interaction becomes time−independent, according to ∙ : ≈ 8 ∙ : O . 8
214 215 216 217
In presence of the high−intensity laser field, the electrons transitions E → p and E ← p are induced
between ground state, i, and excited state, j, with the transition−specific Rabi frequency, @[ > /2l ,
@[ ∙ : !> 8
218
219 220
(15)
(16)
is the electric−dipole−moment matrix element for the transition between the initial where 8
state i and the final state j. Solving multi−level Bloch equation provides oscillating solutions.
221
Occupation of orbitals experiences change in time, with the major contribution provided by
222
oscillating with Rabi frequency. : : O cos >
223
(17)
|o "K O | 0o
224
@[ @[ t cos2> 3 , t E sin2> 3.
225
(18) (19)
226
The evolution of density matrix under such driving field is illustrated in Figure 1(b). As an
227
approximation, a stepwise change in occupations of two interacting orbitals is adopted. This
228
approximation can be rationalized by comparing continuous irradiation Eq. (17) with sequence
229
230
231
of l-pulses, maintaining same average flux of energy. :
: O limU→O ∑O exp
[+T& R, U
]
cos >
|o . . . K K O O O | 0o. 11
ACS Paragon Plus Environment
(20)
(21)
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
232
The evolution of density matrix under such driving field is illustrated in Figure 1(c). This
233
approximation to Rabi theory helps to avoid the mean−field paradox.
234 235
Page 12 of 52
Figure 1 shows the schematic diagram for a sequence of stimulated absorption and
emission events (a) and Rabi flops induced by intense CW irradiation (b) and sequence of l-
236
pulses (c). For a sequence of finite duration π pulses, the occupation of ground and excited states
237
swaps from 0 to 1 and from 1 to 0, respectively. The off-diagonal elements of density matrix are
238
nonzero during the pulse activity and zero during the population period. For a setup of short π
239
pulses, the TDESMD approach will be accurate. For sequence of broad pulses or continuous
240
irradiation, the TDESMD is an approximation. However, we hypothesize that continuous
241
transition from sequence of short pulses to CW excitation will keep TDESMD solution as a
242
leading contribution to the dynamics with small corrections.
243
244
245
2.5 Practical Implementations of Computational Approaches for Light-driven Fragmentation The orbitals ¡ 2 , 3 and orbital energies \ 2 3 can be obtained from the
246
geometry optimized model at initial time t = 0, and from any updated geometry at later times,
247
and used to calculate the electronic dynamics for the subsequent increment of time, ∆t. At each
248 249 250 251
step, the updated density is obtained with a known increment: ∆ ∆ .
Here ∆ is a density increment gained during finite time step, which accounts for the transitions
between electronic degrees of freedom. Then, with the new density matrix ∆ three key variables characterizing the model are updated: (i) the new total charge density, (ii) the new
252
orbitals and their energies, and (iii) the new interatomic forces and ion positions. The density
253
matrix can be explicitly propagated in time by diagonalizing the Liouville−von Neumann and 12 ACS Paragon Plus Environment
Page 13 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
254 255
Journal of Chemical Theory and Computation
Redfield super−operator, , at each time step.42, 71-75 Here, propagation of electronic degrees
of freedom is modeled as an approximation of optically driven Rabi oscillations with
256
instantaneous and discrete hops in the elements of density matrix, similar to surface hopping
257
methods but not identical.36, 68, 76-78 It is hypothesized that at initial stage of the photoreaction the
258
nuclear configuration is relatively close to equilibrium geometry and used approximations are
259 260 261 262
263 264
valid. As time goes by system may depart from equilibrium and from validity region. ¢
£ could be a criterion of departure from equilibrium. A stepwise change in ∑ A^
occupations of two interacting orbitals i and j with time introduced in subsection 2.4, due to oscillations with the Rabi frequency, is introduced as, ∆+
1 → 0, ∆ 1 ∆+
0 → 1, ∆ 1 .
(22) (23)
265
The Fock matrix is updated at each time step, following changes in forces, geometry, and
266
density. The Fock potential is updated when the nuclei move. Changes in the electronic
267
configuration and in the potential energy surface determine the changes in the interatomic forces.
268 269 270 271 272 273
274
Thus, the nuclear system represented by the position of the nuclei and the electronic
system represented by density matrix, , evolve in time in a coupled way. As a result, the
time−dependent density matrix, , experiences a rapid change in time due to light
perturbation for transitions that match the laser frequency !>
[email protected] ≈ \ \ ). The total density is then used in the Kohn−Sham self−consistent procedure, and determines the energy gradient and forces imposed on each nuclei of the model, '
'+
1 [ ]/8 . 13
ACS Paragon Plus Environment
(24)
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
275
Page 14 of 52
Light−induced changes in the electronic configuration and forces 1 [ ] , which
276
depend on the phase of incident light, allow the nuclear positions to overcome the dissociation
277
barriers and result in a rich fragmentation dynamics trajectory. An optically−driven change of
278
electronic state has a smaller random component than trajectories driven by electron−to−nuclei
279
interactions, because the probability of electronic transitions driven by photons is higher than
280
that driven by phonons. All trajectories in the swirl would be the same. Trajectories computed
281
for a vanishing change in initial geometry will produce similar results. Therefore, we propagate
282
only one trajectory, where the electronic transitions are induced by the laser field. In case
283
simulation starts from the equilibrium geometry, the multiple trajectories will coincide since the
284
initial conditions for positions and momenta of each ion are the same at 0 K. We can prove that
285
for 0 K multiple trajectories are not needed, since the external perturbation is the same. In
286
addition, we have tried trajectories and stochastically identical trajectories and found that at
287
nonzero temperature trajectories give less accurate mass spectra. See Figure S5-S7 in
288
Supporting Information.
289
290
2.6 Computed Observables Characterizing Light-driven Fragmentation
291
The TDESMD simulation allows computation of a nuclear trajectory for a given initial
292
set of conditions. A trajectory includes several hundred time steps with a different nuclear
293
configuration at each step, representing the intermediates and products of the photofragmentation
294
reaction. At each step, the interatomic bond distance for each pair of nuclei is used to evaluate
295
whether a chemical bond is breaking or forming, and thus to determine the number of fragments,
296
¥C@¦ , number of nuclei in each fragment, §C@¦ , and the chemical composition of each fragment, ¨
14 ACS Paragon Plus Environment
Page 15 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
297
298 299
Journal of Chemical Theory and Computation
with all above being time dependent. A single list of nuclei in the model is replaced by a separate list of nuclei for each fragment, K ..
©ª«
©ª« → ¬ K .. , K .. , … , ®©ª« ®©ª« K ..®©ª« ¯ . xyyyyyzyyyyy{ xyyyyyzyyyyy{ xyyyyyyyyzyyyyyyyy{ ©ª« ©ª« ©ª«
¨K
¨B.+ ¨B@CA+
¨£
¨¨
.A|}' ¨B@CA+
?@.+ ¨B@CA+
300
(25)
301
Here, ¥ is the index for a fragment, such that 1 ≤ ¥ ≤ ¥C@¦ . §C@¦ is the number of nuclei in
302 303 304 305 306 307 308 309
310 311 312
¨
fragment number ¥, and § ¨ is an index labeling nuclei in a fragment number ¥, such that
¨ 1 ≤ § ¨ ≤ §C@¦ . ® is the position of nucleus number § ¨ from the fragment number ¥. Note
that the composition and size of each fragment is reconfigured with time. Although this treatment is general, we apply it for monitoring the development of the heavy fragments
containing the metal ion. §C@¦ 1 ± ², corresponding to the chemical composition ¨K
LaCx(t)Hy(t).
Thus, simulated mass spectra can be extracted from TDESMD trajectories. 8¨ is the
molecular mass of fragment, ¥, and can be calculated according to, 8¨ ∑©ª« ® K 8 ®
®
+
where 8® is the mass of ion number § ¨ belonging to fragment ¥.
The mass spectrum, 8³´µ , , at instantaneous time, t, is expressed as, 15 ACS Paragon Plus Environment
(26)
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
® 8³´µ , ∑¨K Y [8¨ ´µ ]
6 +
313 314 315 316
Page 16 of 52
(27)
where ¶¨ is number of fragments, and ¥ is an index labeling fragments 1 ≤ ¥ ≤ ¶¨ . 8¨ is the mass of a given fragment, f. ´µ is an argument, which scans all possible values of molecular
weight. The Dirac delta function, Y[8¨ ´µ ], is approximated by a finite−width Lorentzian
317
function, Y±
318
the Delta function which gives the width of the distribution. The value σ = 0.1 was used to
319
simulate spectral line broadening in experimental mass spectra. The broadening parameter is set
320
as 0.1 to make simulated mass spectra comparable to experimental photoionization
321
time−of−flight mass spectra. This is a regular practice in studies where computational and
322
experimental spectra are compared. The value of 0.1 was a minimal value to get a close−to−real
323
representation of computed spectra. If no broadening parameter is used, then the simulated mass
324
spectra are shown as individual lines. If larger value of broadening is used, then several peaks
325
merge into unresolved feature. Chosen value of σ is dictated by resolution of mass spectrometer.
326 327
328
K
·
· R¦
, where σ is a parameter with the same dimensions as the argument of
The simulated mass spectra, < 8³ > ´µ , are calculated by performing an average along the duration of the trajectory in equation (27):
< 8³ > ´µ L MO 8³´µ , N. K
L
(28)
329
Later, at each step of reaction, the geometry optimization of the intermediate is
330
performed to remove the kinetic energy. We optimize intermediates and report their energy so
331
that there is no kinetic energy of bond contraction and bond elongation. The TDESMD algorithm
332 333 334
generates “hot” trajectories with substantial momentum for each ion ¸ . During translation of the fragments from the PI volume to the detector, some of the kinetic energy can be
adiabatically dissipated, so that fragments may approach their equilibrium configurations with 16 ACS Paragon Plus Environment
Page 17 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
335
zero momentum and zero kinetic energy. Detailed mechanism of such dissipation may include IR
336
radiation and collisions. The flight time of the heaviest fragment La(Cp)3+ is about 27 µs. A
337
quantitative estimate of the rate of IR radiation would depend on Einstein coefficients for
338
spontaneous emission. Even very long−lived states have radiative lifetimes shorter than 27 µs.
339
So by arrival to detector the fragment is expected to cool down. From the technical point of view,
340
“hot” fragments have a huge uncertainty of the total energy. That is fluctuations of kinetic energy
341
and total energy along hot trajectory exceed the difference of total energy between fragments.
342
The instantaneous “cooling” (optimization) is an artificial procedure needed to analyze
343
intermediates and to minimize structural ambiguity present at finite temperatures. Note that there
344
is not “cooling” along computation of the TDESMD trajectory. The redistribution of kinetic
345
energy is made away from computation of the TDESMD trajectory and does not serve as an
346
input for next step of trajectory. “Cooling” is applied to snapshots of the “hot” trajectory that are
347
“extracted” from the dynamics and never “returned” back to the reaction. This “cooling” models
348
the dissipation that occurs long time after the reaction has been complete. The reaction occurs at
349
the laser focus on ps to ns time scale. “Cooling” occurs at the channel of flight towards detector
350
on µs time scale. The photoreaction occurs in the gas phase, so the energy is not dissipated from
351
the molecule. The dissipative process plays a secondary role in the photo-fragmentation. But,
352
after the photo-fragmentation, there is a long-time period before a fragment reaches the detector.
353
The “cooling” of the fragment relates to fragments arriving to the detector. A schematic diagram
354
about the computational procedure of TDESMD algorithm and processing of the output is given
355
in Figure 2. The TDESMD algorithm, on its own, does not include artificial velocity rescaling
356
and a thermostat, so it is not a subject to “flying ice cube” artifact.79 Any changes in nuclear
357
momentum/velocity originate from moving along potential energy surfaces and are subject of
17 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 18 of 52
358
precision of computation of total energy, potentials, and gradients. Each fragment, introduced by
359
equation (25) can be characterized by their total energy and a record of intact bonds. We
360
hypothesize that there is a correlation between number of bonds and total energy of an
361
intermediate as introduced by linear equation (29),
362
: :0 /[¹0 ¹] º[»0 »] ¼[½0 ½] ¢[N0 N] (29)
363
where the regression coefficients A, B, C, and D stand for proportionality coefficients between
364
number of bonds and the energy for different types of bonds (energy per bond). The variables
365
a(t), b(t), c(t), and d(t) stand for the number of bonds of different type at each step of the reaction,
366
such as a = La−C, b = C−C, c = C−H, and d = La−H. For each reaction, the values of
367
coefficients A, B, C, and D are found by solving a system of linear algebraic equations composed
368
of data points at each step of the reaction.
369
370
2.7 Computations Details
371
The atomic model of La(Cp)3 consists of 31 atoms and 86 valence electrons, with the Cp
372
rings η5 bonded to a lanthanum ion. In the simulation cell, a 7 Å vacuum in each direction, x, y, z,
373
is added to minimize interactions between the fragments and penetrations into the simulation cell.
374
All calculations are done in a basis set of Kohn−Sham (KS) orbitals computed in density
375
functional theory (DFT) using the Perdew−Burke−Ernzerhof (PBE) form of generalized gradient
376
approximation (GGA) with the projected augmented wave (PAW) potentials as implemented in
377
the Vienna ab initio simulation package (VASP).80-83 The spin−restricted approach is adopted,
378
based on the consideration that La3+ has the 4f0 configuration. The geometry optimization of the
379
atomic model is performed to get the starting point for adiabatic MD and TDESMD simulations. 18 ACS Paragon Plus Environment
Page 19 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
380
In TDESMD simulations, the system is under continuous laser perturbation, and the electrons are
381
hopping between the ground state and selected excited state with an inverse Rabi frequency,
382
@[ 2l/> = 10 fs. The time increment for simulation is 1 fs.
383
384
3. RESULTS
385
This section is divided into three subsections. Subsection 3.1 focuses on the basic
386
electronic structure of La(Cp)3, which is necessary for the study of pathways of photoexcitations.
387
Subsection 3.2 covers the phonon−assisted processes of charge density migration in energy and
388
in space. By comparing the transition rate and Rabi frequency, we emphasize the importance of
389
light-driven electronic transitions leading to photolysis. Subsection 3.3 focuses on TDESMD
390
simulations on the photofragmentation processes.
391
3.1 Electronic structure
392
The basic electronic structure of La(Cp)3 is shown in Figure 3. The calculated HO−LU
393
gap is 2.86 eV, where the highest occupied Kohn−Sham orbital is abbreviated as HO and the
394
lowest unoccupied Kohn−Sham orbital is abbreviated as LU. In panel (a), features a − c
395
correspond to occupied orbitals, while features d − j correspond to unoccupied orbitals. Each
396
orbital is labeled based on HO−LU notation and illustrated by isosurfaces of partial charge
397
density. LU+n means the nth orbital above LU and HO−n means the nth orbital below HO. The
398
decomposition of Kohn−Sham orbitals on the atomic basis can be used to analyze the nature of
399
transitions between two orbitals of interest. One can notice that the charge density is
400
concentrated on carbons of the Cp ring for features a (HO−3), b (HO−1), c (HO), g (LU+8), i
19 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
401
(LU+13), and j (LU+16). Features d (LU), e (LU+3), f (LU+7), and h (LU+10) show the highest
402
concentration of the charge density is located on the La ion. Specifically, features d and h show
403
hybridization of La and Cp ring. Features e and f show all seven orbitals of La 4f, with the
404
largest contribution from orbitals with magnetic quantum number values mz = 0 and mz = 3,
405
respectively.
Page 20 of 52
406
The calculation of absorption spectra is based on the non−interacting orbital
407
approximation, where the optical property of the model is described by matrix elements of
408
electron−photon interaction in the basis of auxiliary KS orbitals for the ground state.8, 24, 84-88 The
409
calculated absorption spectrum of La(Cp)3 is shown in panel (b). Features a', b', c', and d'
410
correspond to transitions HO−1 → LU, HO−3 → LU+3, HO−1 → LU+10, and HO−3 → LU+10,
411
respectively. These transitions are typical ligand−to−metal charge−transfer (LMCT) transitions,
412
as illustrated by the partial charge density in panel (a). The charge carrier dynamics induced by
413
these transitions are explored later. Studies on several lanthanide metal−organic complexes
414
indicate LMCT transitions facilitate the photofragmentation.14-17, 23, 24, 54 These transitions are
415
thus explored in TDESMD calculations to simulate the laser perturbation of the model.
416
The geometry−optimized model was heated to 450 K, the experimental sublimation
417
temperature, within 300 fs, prior to MD calculations. This is done by reassigning momentum to
418
each nucleus at each time step, which is a standard procedure of velocity rescaling. In MD
419
calculations, the time step is set to 1 fs and the average temperature is held at approximately 450
420
K. Figure 4 shows the fluctuations of temperature and KS orbital energies as a function of time
421
during MD simulations. The average temperature in MD simulation is 448.7 K with a standard
422
deviation of 57.5 K. In what follows we borrow solid state term bandgap for the energy offset of
423
the HO and LU orbitals, also related to minimal excitation energy. The band gap throughout the 20 ACS Paragon Plus Environment
Page 21 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
424
MD simulation is clearly observed. The energies of HO and LU fluctuate by 0.6 eV and 0.4 eV,
425
respectively. Fluctuations of orbitals energies are indirect signatures of electronic response to
426 427 428
nuclear motion. The amplitude of ∆\ represents electron−phonon coupling. It is interesting that larger amplitudes of fluctuation of orbital energies are observed in KS orbitals near the band gap than in HO or LU.
429
430
3.2 Charge carrier dynamics
431
Figure 5 shows representative examples for the autocorrelation function of nonadiabatic
432
coupling and elements of Redfield tensor based on the microcanonical MD trajectory at 450 K.
433
We have inspected phonon/nuclear assisted transitions at low temperatures corresponding to
434
initial stage of evolution. The system studied here falls within the regime of applicability of
435
Redfield theory, because for initial states temperature is low. Autocorrelation drops abruptly. It is
436
interesting to see that the autocorrelations decay rapidly, within a few femtoseconds. The
437
maximal absolute values of the Redfield tensor elements appear between adjacent orbitals with
438
close energies. In addition, one observes that the transitions in the unoccupied orbitals are
439
quicker than those in the occupied orbitals, because the density of unoccupied orbitals is greater
440
than that of occupied orbitals and there are more than seven lowest unoccupied orbitals residing
441
on metal and do not involve change in space. The high orbital density corresponding to small
442
subgap energies provides more relaxation channels. It should be noted that the elements of
443
Redfield tensor connecting the bandgap states are vanishing, indicating a nonradiative lifetime of
444
the lowest excitation ~ 20 ps. It is interesting and important that the average transition rate
445
between any pair of orbitals is less than 0.05 fs–1, corresponding to a relaxation time of 20 fs,
21 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
446
indicating that the quickest pathway of the radiative transition is slower than the adopted Rabi
447
frequency in the TDESMD simulations.
Page 22 of 52
448
Figure 6 shows nonadiabatic relaxation dynamics of La(Cp)3 for the initial condition,
449
corresponding to feature d' (HO−3 → LU+10) in Figure 3. As shown in Table 1, feature d' has
450
the largest oscillator strength. In Figure 6, red, green, and blue represent the distribution for gain,
451
no change, and loss, respectively, relative to the equilibrium distribution. Panel (a) shows the
452
distribution of charge as a function of energy and time. It is observed that there are intermediate
453
states during charge carrier relaxation. Specifically, the electronic part of the excitation relaxes
454
from orbital LU+10 to a lower electronic state LU+3, and then approaches LU. The terms
455
electron and hole, are used below to refer to either additional charge at the originally empty
456
orbital or a lack of charge at the originally filled orbital. The hole part of the excitation relaxes
457
from HO–3 to HO–1, and then approaches HO. The transitions to lower energy excitations
458
follow from dissipating the energy of electronic excitation into the energy of nuclear motion. In
459
the solvent, or at a substrate, such energy dissipates to the environment. Here, in gas phase or
460
vacuum, energy accumulated by nuclei from electron relaxation may contribute to the increase in
461
the overall kinetic energy and temperature. Panel (b) shows spatial distributions of charge as a
462
function of time. Throughout the relaxation, the positive charge (blue) resides on the organic
463
ligands and the negative charge (red) resides on the metal. Note that electron relaxation is much
464
faster than hole relaxation. Electrons relax within 60 fs, while holes relax within 1500 fs. One
465
possible reason is that the density of states is higher in unoccupied orbitals than occupied orbitals,
466
which provides additional channels for energy dissipation. Optical transitions are facilitated by
467
high energy photon, in response with transition energy. For phonon−assisted transitions there are
468
no phonons in resonance to transitions from HO−to−LU, and similar ones. The detail of the 22 ACS Paragon Plus Environment
Page 23 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
469
charge carrier relaxation dynamics for the rest of initial conditions in Table 1 can be found in the
470
Supporting Information (see Figure S1-S3).
471
472
3.3 Fragmentation
473
Figure 7 shows the comparison between experimental PI–TOF mass spectrum of
474
La(Cp)3 and simulated mass spectra extracted from TDESMD simulations using La(Cp)3 as the
475
starting point in the high mass region. The spectra were extracted from TDESMD simulations
476
using a combination of cheminformatics methods and code written in python. In the high mass
477
region, one clearly observes La(Cp)3+, La(Cp)2+, and La(Cp)+. The experimental PI–TOF mass
478
spectrum suggests the reaction pathway of La(Cp)3 → La(Cp)2 + Cp → La(Cp) + 2 Cp → La0 +
479
3 Cp, with neutral Cp ligand ejection through laser excitation into the LMCT states. In Figure 7
480
all simulated mass spectra clearly show features of La(Cp)3+, La(Cp)2+, and La(Cp)+. It is
481
observed that relative peak intensities in the calculated mass spectra depend on the excitation
482
energy. Specifically, it is found that for spectra 5(d) and 5(e) intensities of predicted peaks match
483
well with that of observed peaks in the experimental mass spectrum 5(a). Used methods allow to
484
find correlation between choice of populated PES, reaction pathway, and distribution of products
485
observed in the mass spectra. More detailed analysis of the PES shape is needed to rationalize
486
direction of the reaction.
487
Figure 8 shows the energy diagram for the relaxed intermediates from the TDESMD
488
simulation with electrons hopping between the orbital pair of (HO–3, LU+10) by using La(Cp)3
489
as the starting point. The atomic model at the end of each Rabi cycle is geometry optimized to
490
get the total energy without the kinetic energy of bond contraction and bond elongation. In the 23 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
491 492 493
Page 24 of 52
first 90 Rabi cycles, one observes the ejection and re−adsorption of the intact Cp ligands. Thus, the total amplitude of energy change ∆:, ∆, is small. It should be noted that the ejected
ligands are not removed from the system due to a “trapping artifact” in the simulation model.
494
Note in this study, the dissociated fragment in the TDESMD trajectory is kept in the simulation
495
cell rather than removed to avoid recombination with other fragments. The unphysical
496
recombination events are referred to as “trapping artifact”. In the following simulation, one
497
observes the cracking of Cp ligands, and thus, a larger amplitude of total energy change. The
498
data represented in Figure 8 are summarized in Table 2 and processed according to equation
499
(29), resulting in the following regression coefficients A = −0.1, B = 3.5, C = 2.0, and D = −4.4.
500
TDESMD simulations are also applied to La(Cp), which is an important intermediate
501
during the photofragmentation process, as illustrated previously. Figure 9 shows a comparison
502
of the experimental PI–TOF mass spectrum and the simulated mass spectra extracted from
503
TDESMD simulations by using La(Cp) as the starting point in the low mass region. In the
504
experimental PI–TOF mass spectrum, the most intense peak is bare metal (m/z = 139). There are
505
many peaks which can be assigned to metal carbides or hydrocarbides (e.g. LaC2+ at m/z = 163
506
and LnC3H+ at m/z = 176). The simulated mass spectra reproduce several features of the
507
experimental mass spectra, such as LaC2H+ and LaC3H2+. There are, however, other observed
508
peaks (e.g., LaC+) which are not reproduced in the simulation. Also, there are additional peaks in
509
the simulated mass spectra which are not observed in the experimental mass spectrum.
510
Figure 10 shows the energy diagram for relaxed intermediates from TDESMD
511
simulation with electrons hopping between the orbital pair of (HO–1, LU+5) by using La(Cp) as
512
the starting point. In the first 20 Rabi cycles, one observes the ring opening of the Cp ligand
513
leading to an increase of the total energy. Local minima of total energy are reached when the 24 ACS Paragon Plus Environment
Page 25 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
514
metal is inserted into the ring, forming a six–membered metallo–heterocyclic ring. Subsequently,
515
one observes the cracking of the six–membered metallo–heterocyclic ring, and thus, an increase
516
of the total energy. After 30 Rabi cycles, a relatively stable structure is formed with a four–
517
membered metallo–heterocyclic ring and a La–C2H3 group. In the following simulations, the
518
metal−containing species undergoes cracking. We report the value of the total energy change
519 520
∆:, ∆ associated with each reaction step. It is found that four–membered metallo– heterocyclic ring structure is energetically more favorable than metal alkyls. The data
521
represented in Figure 10 are summarized in Table 3 and processed according to equation (29),
522
resulting in the following regression coefficients A = 0.6, B = 1.0, C = 0.2, and D = −1.0.
523
524
4. DISCUSSION
525
There are several approximations adopted in TDESMD algorithm. (i) DFT is used as a
526
single determinant method for modelling of dissociation reaction. In principle, one would need to
527
utilize multi-configurational/multi-reference approaches to account for non-dynamic electron
528
correlation. It is possible to convert the TDESMD procedure in the basis of configuration
529
interaction (CI) or coupled cluster (CC) methodology with the requirement of substantial
530
numerical resources. (ii) Spin restricted approach is used. (iii) Optical transition has a fixed value
531
of transition dipole and oscillator strength. At zeroth order approximation and at the initial stages
532
of the trajectory, the oscillator strength and transition dipoles are approximately equal to the
533
equilibrium values. This fact is used for maximal simplification of the numerical procedure and
534
is scheduled to be corrected in the next release. (iv) A sequence of one−photon processes is
535
considered. Simulations don’t involve the instantaneous absorption or emission of multiple
25 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 26 of 52
536
photons at given instant of time. (v) Simulations treat electrons quantum mechanically but nuclei
537
classically as point charges thus adopting classical path approximation (CPA). All these
538
approximations may need to be waived in the future.
539 540
541
542 543
544 545
546
547 548
Singly-excited state formation is often described as superposition of occupied and unoccupied molecular orbital pairs e.g. as in solutions Bethe-Salpeter À À ∑ ¾ ¾ [\ \ Y ¾ Y ¾ ¿ ¾ ¾ ] / ¾ ¾ : / xyyyyyyyzyyyyyyy{
(30)
)*)¾ *¾
with many-electron Hamiltonian Á ¾ ¾ in basis of one electron orbitals |Eo and |po . The eigenstate of solutions read |Âo ∑ /À |Eo|po, and obey normalization ∑Ã/À à 1. £
The description of fragmentation dynamics requires re-calculation of such excited states
at each time along trajectory |Âo, /À . As an approximation, superposition is represented as /À ≈ Y, ¾ À ∙ Y,À¾ 0,
(31)
where excited state  is approximated as pair of orbitals, occupied E  to unoccupied p Â,
both depend on Â, referred as independent orbital approximation (IOA).
549
To solve equation (29), it is important to determine the number of bonds at each step of
550
the reaction. Eggers et al. reported the crystal structure of La(Cp)3 with the average La−C bond
551
length of 2.805 Å and the maximum of 2.999(6) Å.89 Ram and Bernath reported that the
552
equilibrium bond length of LaH is 2.032 Å, compared to the theoretical value of 2.08 Å
553
predicted by Das and Balasubramanian.90, 91 The normal C−H bond length is 1.09 Å. However,
554
the elongation of the C−H bond is also reported by others.92 Bau et al. found an elongated C−H
555
bond of 1.153(6) Å for the nitrosyl complex Cp*W(NO)(CH2CMe3)2.92 The threshold values of 26 ACS Paragon Plus Environment
Page 27 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
556
the bond length in equation (29) for La−C, La−H, C−H, and C−C are thus set to be 2.945, 2.132,
557
1.145, and 1.540 Å, respectively, where all except for C−C have 5% longer distance than the
558
reported or traditional value (2.805 Å for La−C, 2.032 Å for La−H, and 1.09 Å for C−H). At this
559
study the empirically−defined bond length threshold is used to determine the bond formation or
560
bond breaking. We do not expect a significant change in composition of fragments with
561
alternative methods use, since many computational quantum−chemically based approaches
562
correlate well with empirically defined bond length ranges for most of the atoms. Bond lengths
563
depend on bond orders significantly. There is a very high correlation between these two variables
564
and using only bond length is often enough to make a right estimation for fragmentation.
565
The features in the PI–TOF mass spectra are charged species. However, it is assumed,
566
based on PI–TOF experiments on similar complexes, that the photofragmentation reaction most
567
likely starts with the neutral species giving rise to neutral fragments.25 These neutral fragments
568
can then be captured by the mass spectrometer following post−fragmentation photoionization.
569
The minor inconsistency between experiment and simulation is probably due to the limits
570
in the implementation of TDESMD methods and due to starting trajectory with incomplete
571
intermediates. In this work, one uses the spin-restricted DFT as an approximation. However, the
572
photofragmentation of lanthanide cyclopentadienyl complexes is indeed radical chemistry. A
573
more accurate description of the photoreaction of this model needs the use of spin-unrestricted
574
DFT, which requires more computational resources. TDESMD calculations including spin would
575
be a subject for future studies. In addition, there are multiple channel (dissociation / ionization)
576
and multiphoton processes in the experiments. These processes are assumed small by the present
577
computational modeling as an approximation. The trajectories and corresponding mass spectra
578
do have specific ranges of most reliable applicability: (i) shorter trajectories starting with the 27 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
579
complete precursor i.e. La(Cp)3 are adequate for predicting the mass spectra in the high mass
580
region e.g. Figure 7; (ii) longer trajectories of trajectories starting with a stable intermediate are
581
suitable for predicting the mass spectra in the low mass region e.g. Figure 9. We do not expect
582
absolute accuracy of the used method. Instead, we report correlations between (i) choice of
583
approximations to molecular dynamics methodology and (ii) computed distribution of products
584
obtained along the trajectory. Such correlation allows the continuous development and
585
systematic improvement of adequate methods.
586
We recognize the challenge of exploring bond breaking with a single-determinant method
587
and discuss ways to mitigate it. Possible improvement is to use different electronic structure tool.
588
Similar research efforts in the community are based on singe determinant methods.93, 94 There is
589
a state-of-art agreement that closed shell single-determinant methods for photochemical
590
dynamics can describe intermediates but not the final products.56, 95 In future we plan to recode
591
the TDESMD dynamics based on several implementations of multi-reference codes, CI and CC.
592
However, such code would work on infinite long. Due to this reason, there are very few multi-
593
reference MD codes available. In fact, multi-reference molecular dynamics is practically applied
594
to small non-metal molecules for very short trajectories.57, 58, 96-98 The most practical goal is to
595
evaluate energies of intermediates and products with multi-reference methods. See Figure 11
596
and Figure S10 in Supporting Information.
Page 28 of 52
597
Currently, we use number of bonds as structural descriptors to correlate with the total
598
energy of intermediates during the TDESMD trajectory. In the future, we plan to add two more
599
structural descriptors: x(t) and y(t), which are the number of carbons and the number of
600
hydrogens, respectively, directly bound to the metal forming a metal−containing fragment
601
LaCx(t)Hy(t). Observed correlations between the number of bonds, the chemical composition, and 28 ACS Paragon Plus Environment
Page 29 of 52
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Journal of Chemical Theory and Computation
602
the energy at each step of the photoreaction may prompt an ambitious future development. That
603
is (i) each step is encoded with the descriptors x(t), y(t), a(t), b(t), c(t), d(t), and E(t); (ii) each
604 605 606 607
step is encoded with a discrete change of these descriptors ∆x, ∆y, ∆a, ∆b, ∆c, ∆d, and ∆E; and
(iii) a map of all possible reaction steps x(t), y(t), a(t), b(t), c(t), d(t), E(t) and ∆x, ∆y, ∆a, ∆b, ∆ c, ∆ d, ∆ E would allow to predict each reaction step without ab initio modeling; (iv)
accumulation of massive statistics of possible reaction maps will result in forming a complete
608
picture of fragmentation photochemistry.
609
5. CONCLUSIONS
610 611
In summary, we have studied the pathways of photoexcitations in gas−phase La(Cp)3, which involve several competing photophysical and photochemical processes.
612
Photophysical processes: the study of the nonradiative, nonadiabatic, dissipative
613
dynamics of charge density illustrates formation of the LMCT excitations and establishes
614
timescales at which the reaction dynamics occur. Since the dissipative processes compete with
615
the fragmentation dynamics, they establish the limits at which reactions are possible. The
616
timescale of nonradiative dissipation of excitation energy is a limit for photoexcitation to occur.
617
It is found that electron relaxation, ke, is faster than hole relaxation, kh.
618
Photochemical
processes:
TDESMD
algorithm
is
applied
to
simulate
the
619
photofragmentation reaction dynamics induced by laser irradiation. The specific application used
620
here as an example illustrates regime of a non-small (medium) perturbation, when the nuclear
621
motion induced by light is more (larger/noticeable) than the one associated with equilibrium
622
thermal motion of nuclei prior to excitations. The values of laser intensity and Rabi frequency
623
are expected to determine the type of light-driven process. (i) At the limit of low laser intensity 29 ACS Paragon Plus Environment
Journal of Chemical Theory and Computation
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
624 625 626 627 628 629
Page 30 of 52
@[ and slow Rabi oscillations with respect to nonadiabatic dissipation, 2l/> > ´¹±gATK , gÄTK ,
the energy of photoexcitation is expected to dissipate, bringing the model to the ground
electronic state without change of chemical composition. (ii) At the limit of medium laser @[ intensity and fast Rabi oscillations with respect to nonadiabatic dissipation, 2l/>