Photoinduced Charge Transfer versus Fragmentation Pathways in

Jun 26, 2017 - Simulated mass spectra are extracted from TDESMD simulations and compared to experimental photoionization time-of-flight (PI-TOF) mass ...
3 downloads 13 Views 2MB Size
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/>