[7] Brent Donovan, Larry A Walker, Charles F Yocum, and Roseanne J Sension. Transient absorption ... [24] Jeffrey A Myers, Kristin LM Lewis, Franklin D Fuller, Patrick F Tekavec, Charles F. Yocum, and ...... [158] Ian P Mercer. Angle-resolved ...
OBSERVATIONS OF COHERENCE IN OXYGENIC PHOTOSYNTHESIS
by Franklin D. Fuller
A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Biophysics) in The University of Michigan 2014
Doctoral Committee: Associate Professor Jennifer P. Ogilvie, Chairperson Professor Eitan Geva, Chairperson Associate Professor Kevin J. Kubarych Professor Roseanne J. Sension
©
Franklin D. Fuller
2014
To my parents, Frank and Cindy Fuller
ii
ACKNOWLEDGMENTS I have had the great pleasure and fortune to legitimately enjoy my PhD work. Of course, it has not all been a bed of roses, no pleasure cruise, but I have had a lot of fun figuring out the problems laid before me and, of course, finding a thousand new problems in the process of solving the first ones. The work here is not done. It will never be done. But I have no regrets. I have many people who helped me along in this journey that I would like to thank. I would like to thank my adviser, Jennifer Ogilvie, who has been a font of positive energy and ideas. She has a great intuition for finding new directions and making suggestions of things to try next when the despair of an unsolved problem sets in. Most importantly she genuinely cares both about the work the lab is doing and who is doing it. Jennifer has gone above and beyond what I expected to help me. Professor Charlie Yocum and his post-doc at the time Hana Popelka are, in many ways, responsible for this work even being possible. They knew how to spin Aunt Mid’s spinach into gold and were kind enough to teach me the ways of the force. I benefited enormously from Daniel Wilcox’s knowledge, thoughtful approach to dissecting problems, and scientific moral compass. The lab as a whole has grown significantly and in the “right” direction thanks to a number of his efforts. I had the good fortune to inherit a 2DES setup and knowledge of that setup from Jeff Meyers and Kristin Lewis. While we ultimately have changed a lot of things and those two would likely roll over in their proverbial graves to see what is now missing or changed, there is no doubt that we are where we are today thanks to their hard foundational work. I wish to thank Seckin Senlik and Jie Pan both for their excellent work ethic and no-questions-asked enthusiasm in helping me collect the data in this work. The “next generation” of students: Anton, Andy, and Veronica have all helped me discover where I lacked in explaining things, and for this I am thankful. Writing this dissertation was made significantly easier by having those discussions. Anton in particular iii
has also opened my eyes to new coding and analysis tricks, which I believe have not only benefited this work, but will also come in handy in the future. And I would be remiss not to mention that Anton has also been a great friend, willing to listen to whatever hair brained idea I currently have while getting coffee or soda. The other half of Jennifer’s group, while less intimately involved with 2D spectroscopy, have been helpful. Amar Bhagwat in particular is very knowledgeable about a great many things that continually catch me off guard and I am grateful he is willing to share that knowledge. Finally, I would like to thank the non-scientific people in my life. Mackenzie Martin has been a constantly supportive and wonderful friend over the last two years while I put in many long hours to generate the data presented here. I attribute a large portion of my remaining sanity to her ebullient personality. My parents are two people which I feel I can trust completely and have always been both supportive and quick to give excellent guidance.
iv
CONTENTS
DEDICATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ii
ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iii
LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
vii
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
ix
LIST OF ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . .
x
ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . I. NONLINEAR SPECTROSCOPY OF PHOTOSYNTHETIC PROTEINS 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 PSII Structure, Function, and Nonlinear Spectroscopic Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Non-linear Spectroscopy of PSII . . . . . . . . . . . . . 1.3 Overview of 2D Measurements . . . . . . . . . . . . . . . . . . .
xi 1 1 2 8 13
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . II. EXPERIMENTAL IMPLEMENTATION OF 2DES . . . . . . . . . . . 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Review of 2DES setups and Signal Isolation . . . . . . . . . . . . 2.3 Experimental Implementation of a 5 pulse setup . . . . . . . . . . 2.3.1 Experimental Results . . . . . . . . . . . . . . . . . . 2.3.2 Analytic Phasing of 2DES data . . . . . . . . . . . . . . 2.3.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Appendix: Improvements to the NOPA-based Laser Source . . . . 2.6 Appendix: Comprehensive Analysis of Signals Generated in Hybrid Setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19 24 24 26 31 35 37 40 43 43
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . III. COHERENT DYNAMICS OF PSII REACTION CENTER . . . . . . . 3.1 Two basic spectroscopic models . . . . . . . . . . . . . . . . . . . 3.1.1 Electronic Dimer Model . . . . . . . . . . . . . . . . . 3.1.2 Displaced Harmonic Oscillator Model . . . . . . . . . .
50 58 58 61 65
v
49
. . . . . . .
68 71 72 76 78 82 88
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . IV. POST ANALYSIS OF 2DES DATA . . . . . . . . . . . . . . . . . . . . 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 The Variable Projection Functional . . . . . . . . . . . . . . . . . 4.3 Using 2D model functions to reduce the problem size . . . . . . . 4.4 Using Orthogonal Decomposition to Reduce Problem Size . . . . . 4.5 Goodness of Fit: Taylor Expansion of Objective Function . . . . . 4.6 Sum of Lorenzians Model . . . . . . . . . . . . . . . . . . . . . . 4.7 Optimal Model Order . . . . . . . . . . . . . . . . . . . . . . . .
89 96 96 97 104 107 110 112 113
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . V. EFFICIENT SIGNAL EXTRACTION FOR 2DES . . . . . . . . . . . 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Adaptive Filtering Algorithm . . . . . . . . . . . . . . . 5.1.2 Extending Adaptive Filtering to Two and Three Dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Applying Adaptive Filtering to Non-Lorenzian Lineshapes . . . . 5.3 Applying Adaptive Filtering to Data . . . . . . . . . . . . . . . . 5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
115 116 116 118
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . VI. CONCLUSIONS AND FUTURE WORK . . 6.1 Summary and Innovations of this Work 6.2 Future Work . . . . . . . . . . . . . .
131 133 133 134
3.2
3.3
3.1.3 Phase of Waiting Time Coherences . . . . . . . 3.1.4 Predicting the effect of static disorder . . . . . . 3.1.5 Literature Review on coherence origin protocols Coherence Data of d1d2 . . . . . . . . . . . . . . . . . . . 3.2.1 High Frequency Modes of d1d2 . . . . . . . . . 3.2.2 Low Frequency Modes of d1d2 . . . . . . . . . Conclusions and Outlook . . . . . . . . . . . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
125 127 129 129
BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
vi
LIST OF FIGURES
1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
Anatomy of a Chloroplast . . . . . . . . . . . . . . . . . . . . . . Selective Survey of Thylakoid Proteins involved in Photosynthesis Crystal structure of pigments in the d1d2 reaction center . . . . . . Visible Absorption Spectrum of PSII RC . . . . . . . . . . . . . . Initial Studies of PS2 RC with 2DES . . . . . . . . . . . . . . . . Feynman Diagram Tutorial . . . . . . . . . . . . . . . . . . . . . Cartoon 2DES Spectrum . . . . . . . . . . . . . . . . . . . . . . . Cartoon 2DES Dynamics . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
3 3 6 7 12 16 17 19
2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8
box-CARS experimental setup . . . . . Pump-probe experimental setup . . . . . Experimental Setup and Signal Diagrams Phase Stability Measurement . . . . . . Demonstration Experiment on Chl a . . Phasing Protocols . . . . . . . . . . . . Optimal NOPA external crossing angles Double Pass NOPA optical layout . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
29 29 32 35 36 41 46 48
3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 3.10 3.11 3.12 3.13 3.14 3.15
Energy level diagram of Electronic Dimer . . . . . . . . . . ED Model: Peak Locations . . . . . . . . . . . . . . . . . . Displaced Harmonic Oscillator Model . . . . . . . . . . . . Vibrational coherence signal pathways . . . . . . . . . . . . Inhomogeneous broadening and vibrational coherences . . . A model for electronic and vibrational coherences . . . . . . Example 2DES of d1d2 . . . . . . . . . . . . . . . . . . . . High frequency rephasing amplitude maps of d1d2 . . . . . . High frequency non-rephasing amplitude maps of d1d2 . . . High frequency rephasing phase maps of d1d2 . . . . . . . . Literature Exciton Models of d1d2 . . . . . . . . . . . . . . Low frequency rephasing coherence amplitude maps . . . . . Low frequency rephasing coherence phase maps . . . . . . . Simulated Coherence Maps and Effect on Charge Separation Initial study on 800 nm anion band . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
63 64 67 68 71 75 77 80 81 83 85 86 87 90 91
5.1 5.2
Simulated inhomogeneously broadened lineshape . . . . . . . . . . . . . 128 Reconstruction of simulated lineshapes with various algorithms . . . . . . 128 vii
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
5.3
Reconstruction of Chl a Rephasing using Adaptive Filtering . . . . . . . . 129
viii
LIST OF TABLES
2.1 2.2 2.3
Possible Pulse Interactions Table 1 . . . . . . . . . . . . . . . . . . . . . 51 Possible Pulse Interactions Table 2 . . . . . . . . . . . . . . . . . . . . . 52 Possible Pulse Interactions Table 3 . . . . . . . . . . . . . . . . . . . . . 53
3.1 3.2 3.3 3.4
Coherence Map Literature Review 1 Coherence Map Literature Review 2 High frequency coherences of d1d2 . Low frequency coherences of d1d2 .
ix
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
73 74 78 85
LIST OF ABBREVIATIONS 2DES Two Dimensional Electronic Spectroscopy AFD Adaptive Fourier Decomposition BBO β -Barium Borate BRC Bacterial Reaction Center CCD Charge Coupled Device Chl Chlorophyll CMOS Complementary metal–oxide–semiconductor CP Cross Peak CS Charge Separation CT Charge Transfer DO Diffractive Optic ESA Excited State Absorption ET Energy Transfer FFT Fast Fourier Transform GSB Ground State Bleach LO Local Oscillator NOPA Non-collinear Optical Parametric Amplifier DOPA Degenerate Optical Parametric Amplifier Pheo Pheophytin PSII Photosystem II PSI Photosystem I PSII RC Photosystem II Reaction Center SVD Singular Value Decomposition TA Transient Absorption TRF Time Resolved Fluorescence
x
ABSTRACT
Observations of Coherence in Oxygenic Photosynthesis
Chair: Jennifer P. Ogilvie The field of two dimensional electronic spectroscopy (2DES) is rapidly advancing, both in theory and implementation to tackle increasingly complex and delicate problems. In the past seven years, observations of coherent or wave-like dynamics in 2D spectra of photosynthetic antenna has captured the imagination of many practitioners in the field, from theorists to experimentalists. Two questions are being raised: what is the origin of coherent dynamics in photosynthesis and, more importantly, do they matter for the function of biological systems? For certain photosynthetic antenna systems there is now considerable evidence and theoretical backing to suggest that coherent dynamics have a positive functional impact on energy transfer. Less explored is how such dynamics may influence charge separation, the primary purpose of photosynthetic reaction centers. Coherent signals are typically weak and difficult to resolve from population dynamics. To address this issue, we developed a method to collect 2DES which has dramatically improved the signal to noise over previous implementations. The new method has been applied to the photosystem II reaction center (PSII RC). The PSII RC is the photosynthetic enzyme uniquely capable of using solar energy to split water. As such it is an important system both for basic plant science and renewable energy generation. With this technique, we find eight coherent modes in PSII RC in the first such report of coherent dynamics on this system. Most of the wave-like motions are assigned to be of vibrational character while four are assigned to a mixture of vibrational and electronic character. Based on supporting simulations it is shown that charge separation is enhanced by the inclusion of such mixed character modes.
xi
CHAPTER I
NONLINEAR SPECTROSCOPY OF PHOTOSYNTHETIC PROTEINS
1.1
Introduction
Photosynthesis is the process by which plants and other photosynthetic organisms convert solar energy into chemical energy. Although the photosynthetic machinery varies among organisms, the basic architecture consists of light-harvesting antennae arrays that gather solar energy and funnel it to reaction centers[3]. Within reaction centers the energy is converted to a charge separated state that drives the later stages of photosynthesis. In oxygenic photosynthesis in plants, algae, and cyanobacteria the Photosystem II Reaction Center (PSII RC) is responsible for transforming solar energy into a charge separation event which in turn powers a catalytic reaction to split water into hydrogen and electrons for use in producing ATP. The oxygen we breathe and depend upon is a by-product of this catalytic reaction. Although the ability to split water and respire oxygen is shared among many completely different organisms, accomplishing this task requires a delicately tuned device which evolution has perfected and conserved in structure for billions of years[39]. Given the difficulty of making a molecular machine that efficiently catalyzes water splitting using only solar energy and its obvious facility in the production of solar fuels, understanding the function and structural design concepts of PSII is of interest for designing artificial photosynthetic devices[1]. The specific structure of PSII conserved across species. Pigments in the RC are organized specifically to produce a shallow energy landscape which facilitates transfers 1
of energy and charge in a millionth of a millionth of a second (or a picosecond, which is called the ultrafast time scale). These initial energy transfer and charge separation events are not only fast, but efficient, with quantum efficiencies exceeding 95%[3]. Once a hole is created in the reaction center, it is transferred via a tyrosine residue to a cluster of four Manganese atoms in the Oxygen Evolving Complex (OEC) where, after a cycle of four such transfers, water is finally split. An intense amount of research spanning decades has been spent on understanding all aspects of PSII, and while the broad strokes of what happens in PSII is understood, many finer questions remain that are critical for implementing a device of comparable performance ourselves. For instance, even the nuclear structure of the catalytic site in the OEC remains an open question[32] (although the structure is highly constrained and with new techniques may soon be solved[18]). The excitonic states of the RC, which describe the coupling between electron orbitals of the pigments, are also not known with certainty. A number of empirical models exist to describe the kinetics of excitonic interactions[19, 120, 27, 107, 33], but they have not yet been able to fully account for features in recent two-dimensional spectroscopic measurements on the RC[23]. More specifically, the sequence of charge transfer events leading up to the transmission of electrons to the quinone pool still disputed[19, 14, 120]. While the nuclear structure of the entire system is characterized[36], along with the principle nuclear vibrational modes[111], whether those vibrations interact strongly with excitons and how that interaction might affect the function PSII is an open question.
1.2
PSII Structure, Function, and Nonlinear Spectroscopic Measurements
PSII is just one small piece in a much larger organization. Beginning at the chloroplast level (the organelle responsible for energy production in higher plants), the photosynthetic apparatus is essentially a system of nested membranes. The nest of membranes inside chloroplasts are called the thylakoid membranes for their vaguely sack-like appearance. 2
Stroma Thylakoid
Grana Thylakoid
Chloroplast 1.1 Figure 1.1: A cartoon of a chloroplast is shown, with a cutaway to expose the stroma and grana thylakoid membranes inside.
Figure 1.2: A selected survey of Thylakoid proteins and a graphical depiction of the processes they control in the thylakoid membrane. Adapted from [5]
3
Two main types of thylakoid membranes exist in the chloroplast: grana, which are columns of tightly packed disks, and stroma thylakoids, which are longer and more loosely packed. See figure 1.1 for reference. Cyanobacteria lack organelles, but even so their external membranes are also contorted into many folds. The reason for these membrane structures, in both cyanobacteria and higher plants, is to increase the surface to volume ratio. Membranes are a means to isolate one region of the thylakoid stroma from an inner region (the lumen) in order to create a chemical potential (across the membrane). In photosynthesis, as with many other biological systems, chemical potential is stored in the form of a pH difference, i.e. the quantity of dissolved hydrogen ions. As the reactions occur at the interface between the two regions, in order to increase the number of energy producing units it is desirable to have a high surface to volume ratio. Ultimately, the proteins living in the thylakoid membranes exist to produce ATP, the cellular energy currency. At the protein level, there are still many mechanisms, length and time scales, and functions of photosynthesis to understand. In contrast to solar cells, it would be inaccurate to say that the study of energy transfer or charge transfer captures the majority of what happens at these membrane interfaces. Converting sunlight to charge carriers and then chemical potential requires a number of supporting systems, that enable to adaptation to changing light conditions, means to repair proteins that are damaged, and machinery to further process the created chemical energy into a form usable by the rest of the organism. Many of these support systems are either embedded in or extrinsic to the thylakoid membrane. A cartoon summary of a portion of the proteins involved in the generation and conversion of chemical potential to ATP is given in figure 1.2, adapted from [5]. This study will focus on PSII, a protein complex whose main purpose is to collect energy, transfer it to the reaction center, and produce a stable charge separation that drives the later stages of the photosynthesis cycle. PSII is found primarily in the grana thylakoid membranes of higher plants, and
4
particularly in spinach, this can be exploited to separate PSII from the other main photosystem: PSI[2]. PSII consists of two bound antenna proteins CP43 and CP47, two tightly bound subunits D1 and D2, which house the reaction center pigments, and PsbO, which stabilizes the Manganese cluster mentioned in section 1.1. Bound within these proteins are optically active chromophores, which include 35 chlorophyll, 2 pheophytin, and 12 beta carotene. The reaction center itself houses 6 chlorophyll, the two pheophytin, and 2 beta-carotene which are tightly packed together with distances ranging from as little as 8 angstroms between PD1 and PD2 to about 70 angstroms from ChlZ D1 to ChlZ D2 [36]. It is this tight packing which causes the valence orbitals of the chromophores, whose transition energies are in the visible, to be coupled together. This means that an incident photon will not simply excite one chlorophyll or pheophytin, but a collection of them simultaneously. Such a simultaneous excitation is called an exciton. The 8 pigments in the reaction center are depicted in figure 1.3. As these pigments will be often referred to by name, we describe them here: PD1 and PD2 are two chlorophyll that are the closest analog to the special pair in the bacterial reaction center (BRC). ChlD1 , ChlD2 , ChlZ D1 and ChlZ D2 are also chlorophyll, where D1 or D2 denotes which protein they are bound to. PheoD1 and PheoD2 are pheophytin (similar to chlorophyll, but lacking a coordinated Magnesium). It is actually possible to selectively purify the PSII core complex including CP43 and CP47 and the reaction center, where these antenna proteins are removed from spinach in a preparation called d1d2-cyt.b559[37], which we colloquially shorten to d1d2. In this preparation, the only optically active pigments are only 6 chlorophyll, 2 pheophytin, 2 carotenoids, and the heme in cytochrome. Thus, the structure in figure 1.3 is representative of the d1d2 cyt. b559 complex, which the preparation studied in this work. Other preparations, from cyanobacteria Synechocystis[30] or T. Elongatus[21] also exist and are popular in mutagensis and crystallization respectively, but are typically more difficult to produce. The reaction center has a number of unique features in the visible spectral region,
5
Pheo D1
Pheo D2
Chl D2
Chl D1 PD2
Chlz D1
Chlz D2
PD1 Yz
Figure 1.3: A crystal structure of the main pigments in d1d2 taken from [36]. Also shown is Yz , the tyrosine residue which conducts holes to the Manganese cluster. Not shown are two carotenoids on the periphery of the protein and the phytol tails of the chlorophyll/pheophytin. Color coding is arbitrary, but helps the viewer see which pigments belong to the D1 protein or the D2 protein at a glance.
which yield different information. Chlorophyll absorbs primarily in the deep blue (Soret band) and the red (Qy band), giving rise to the green appearance of chlorophyll and most plant life. Most studies focus on the Qy band, as all the processes associated with charge separation are energetically associated with this band (see figure 1.4 panel B). An absorption spectrum of the purified reaction center is shown in panel A of figure 1.4, where a number low amplitude bands are clearly visible. The Soret band is primarily used for pigment stoichiometry measurements, where in the case of purified PSII RC, an absorption ratio of A416 /A435 = 1.2 indicates a sample with 6 Chlorophyll and 2 pheophytin[9]. At 545 nm, there is an isolated absorption feature associated with the Qx transition of pheophytin and is an important marker for determining both pheophytin excitation and the formation of a pheophytin anion. Studies of chemically reduced chlorophyll and PSII reaction centers have found two other anion bands, one broad absorption peaking at 455 nm and another from 790 to 820 nm[38, 15]. It is noted in [13], however that all of these anion bands are spectrally overlapped with excited state absorption features and thus are not perfect markers for observing charge transfer.
6
293 K Linear Absorption of D1D2 RC
A
Absorption (O.D)
0.25
0.2
0.15
0.1
0.05
0 400
1.75
450
B
500
550 600 Wavelength (nm)
650
700
P* PD1+PheoD1-
1.50
Free Energy (eV)
3
P
PD1+QA-
Yz+QA-
1.25
1.00
0.75
0.50
0.25
0
P
Figure 1.4: A: Visible absorption spectrum of the purified PSII reaction center, taken at room temperature. The optical density (O.D) here is representative of that used in studies to be presented later. B: An energy diagram adapted from [4], depicting the first 3 main steps of the reaction center up to the transfer of a hole to + Pheo− Yz en route to the Manganese cluster. Note that a PD1 D1 can readily decay to a triplet state in reaction center preparations where transfer to QA is blocked and this triplet state is nearly isoenergetic with the product state of the reaction center.
7
1.2.1
Non-linear Spectroscopy of PSII
In order to compete with the known nanosecond excited state lifetimes of chlorophyll, any of the initial photoreactions constructed from chlorophyll aggregates must proceed on timescales that are at least 10-100 times faster[34]. Understanding the processes which immediately follow photoexcitation in photosynthetic systems therefore requires an instrument capable of resolving sub nanosecond dynamics. A number of methods exist to probe such these dynamics, but the most commonly employed methods in the study of the PSII core and reaction center complexes have been time resolved fluorescence and transient absorption, with the last being the most popular and enduring method as it has been used to resolve dynamics down to femtosecond timescale[8]. Both of these methods give slightly different, but complementary information. As mentioned in the beginning of this section, anion bands offer a chance to directly monitor the formation of charge separated products and so these bands have been studied extensively. In the following I will review the works on each band and the conclusions reached by those studies. Transient absorption studies of ground state bleach in the pheophytin Qx band at 545 nm has been monitored under a number of excitation conditions over the years. In the early 1990s a 20 ps charge separation time was found for a room temperature study using selective and non-selective excitation of the Qy band[17] and confirmed in the mid 1990s[7, 22] along with a weaker 3 ps component. The same band probed with numerous excitation conditions at 77 K in [120], is slightly faster with the main component at about 17-18 ps, and faster components down to 1.5 ps. It is interesting to note that in this study, 680 nm excitation results in a dramatically slower 40 ps lifetime of the Pheo Qx band. In a comprehensive study[19] relying heavily on kinetics from this band, 1.5 ps components were identified as energy trapping in the reaction center and primary charge separation was assigned to a 5.5 ps component. Fluorescence of the reaction center monitored by time resolved fluorescence (TRF) measurements or stimulated emission of a the Qy(0,1) band indirectly give information 8
about charge separation through the depopulation of singlet states which are able to fluoresce to charge separated states which do not. Streak cameras offer the fastest measurement rates for TRF, with response times as small as 4.5 ps[6]. Faster response times can be achieve through fluorescence upconversion, though the number of these studies on d1d2 or core complexes has been limited[40]. The conclusions from these studies, particularly for early ones, was that the charge separation was occurring on a time scale of 30 ps or less. By using deconvolution of the instrument response, researchers in [6] concluded that both a 1.5 ps component and a longer, multi-exponential component best described the time scale of charge separation, with the faster component being significantly larger than reported for the Qx band (50% vs 30%). Stimulated emission at low temperature (15 K) corroborates with a shorter time component (2 ps)[38]. Transient absorption studies have examined the Qy band, and as one might expect the kinetics are multi-exponential. 20-27 ps components were identified in the Qy band in [28], along with 3 ps components[38]. More modern studies at 77K have identified a myriad of components in the Qy and other bands, ranging from 1.2 to 2 ps, circa 20 ps, 250 ps and nanosecond components[120]. Early magic angle polarized studies in the mid 1990s identified a 30 ps energy transfer component in the Qy band[12] at room temperature. Markers of charge separation in the infrared were studied primarily by Groot and co-workers[14]. In this study it was concluded, based on vibrational modes from 1585 to 1775 cm−1 , that photoreduction of pheophytin occurs faster than other studies with a + is estimated from the data rate of 0.6-0.8 ps. A circa 6 ps lifetime for the formation of PD1
and the conclusion is that ChlD1 and PheoD1 complex to form the initial charge separation event. It is clear that a long history of measurements on reaction centers and core complexes have been made, and the brief review here is hardly doing justice to them. An excellent review of early work is found in [13], and another of work leading up to the early 2000s in
9
[40] which also includes a review of mutant studies not mentioned here, and more recent reviews in [29, 20]. Despite all this work, there are still opposing views on what model best describes the initial charge separation pathway. The 20-50 ps lifetimes found in nearly all studies are thought to result from energy transfer to the RC and the multi-exponential time components seen in anion bands reflect the cascade of energy transfers from antenna and peripheral chlorophyll. The fast components observed intuitively result from direct excitation of the reaction center. The dispersion of reported fast rates, ranging from sub picosecond through 3 picoseconds, however open the door for a great deal of speculation. Two main models are present in the current transient absorption literature. One group strongly advocates both pathways to charge separation[120] shown in (1.1) and (1.2), evoking a multi-pathway picture that is understood to be true in the bacterial reaction center[16]. Others advocate only the ChlD1 pathway in equation (1.1) as the single dominant pathway[19, 20].
+ + − Pheo− [ChlD1 , PheoD1 ]∗ → ChlD1 D1 → PD1 PheoD1
(1.1)
+ + − [PD1 PD2 ,ChlD1 ]∗ → PD1 Chla− D1 → PD1ChlaD1
(1.2)
In 2010, Myer and co-workers[24] reported a two dimensional spectroscopic study of the PSII RC examining the Qy band that tentatively supports a multi-pathway model by observing a spectrally dependent dispersion of rates. Subsequently I made measurements of d1d2, expanding the probe bandwidth to cover from 700 to 455 nm with the aim of observing dynamics in anion bands and the pheophytin Qx band[10]. The Qy band of this data was analyzed via simulations and again, tacit support for multiple pathways was found[100]. While signal was visible for the anion bands, it was too weak to offer interpretable results along the excitation axis. In figure (1.5), the 2D spectrum of d1d2 is shown for the regions from 450 to 640 nm in panel A and from 650 to 690 nm for comparison in panel B. As in the original study[24], the Qy band exhibited excellent signal
10
to noise, but the signals in the blue, particularly at the Pheo Qx bleach band, are nearly 100 fold weaker, which places the signal to noise of our measurement for these regions on the order of unity. In the next chapter we develop a new method to collect 2D spectra that dramatically improves the signal to noise achievable, with the aim to study weak, yet information rich transitions.
11
660 640
A
Detection frequency (THz)
620 600 580 560 540 520 500 480 425
430
435
440
445
450
455
460
450
455
460
Excitation frequency (THz)
460
Detection frequency (THz)
455
B
450 445 440 435 430 425 425
430
435
440
445
Excitation frequency (THz)
Figure 1.5: A 2DES measurement employing a Calcium Fluoride generated continuum probe examining the PS2 RC at a waiting time delay of 3.5 picoseconds. A: the probe region from 455 nm to 640 nm. B: the probe region from 650 nm to 690 nm.
12
1.3
Overview of 2D Measurements
The introduction of two dimensional electronic spectroscopy (2DES) in section 1.2.1 deserves some exposition, as for the remainder of this work I will be using 2DES to examine the PSII RC. Two-dimensional electronic spectroscopy (2DES) has become a popular extension of visible transient absorption (TA) spectroscopy. By resolving the excitation frequency in addition to the detection frequency, 2DES aids in understanding condensed matter processes such as solvation dynamics[56] and excitonic energy transfer in photosynthetic complexes[46, 62], J-aggregates[48] and semiconductor materials[75, 76]. In concept, 2DES expands TA by using two pump pulses with variable time delay t1 between them, rather than a single pulse, and then probes the system after a waiting time t2 . The result of splitting the pump pulse into two is that the measured transient signal oscillates as a function of t1 , and the Fourier transform with respect to this delay produces the ω1 axis (or excitation axis). As in frequency-resolved TA, the signal is typically detected in the frequency domain, yielding the ω3 axis (detection axis). The reason the detected signal oscillates with t1 is because the first pump pulse puts the system into a coherence: a super-position between the initial (usually called ground) state and some excited state. Left alone, this coherence will oscillate between a mixture of mostly ground or mostly excited state or somewhere in between. When the system is “observed” by the dipole operator via interaction with the probe field a coherence will emit a photon with energy equal to the difference between the two states that it was in a coherence with, which is to say the photon emitted will have a well defined frequency. This is equivalent to saying that, in the time domain, the output polarization will oscillate with a frequency equal to ωeg , where h¯ ωeg is the aforementioned photon energy. The concept that a field interaction produces a coherence, which can then evolve and potentially re-emit a photon comes from a perturbative picture of quantum mechanics and is conveniently described for open or closed quantum systems by a density matrix formalism. The density matrix is essentially an outer product of all possible quantum states in the system in a basis for 13
those states that is convenient. If a closed quantum system comprised of only two states A and B is considered, then a convenient basis for such a system would be {|A , |B}, and so the density matrix would take the form: ⎡
⎤
⎢ |A A| |B A| ⎥ ρ = ⎣ ⎦ |A B| |B B| The off diagonal elements of ρ describe a coherent superposition of the states A and B. In so far as the basis completely describes the system, then so too does the density matrix describe the system. 2DES, like transient absorption, is a third order process, which means that three field interactions are required to generate the desired signal. The term “field interaction” refers to a time dependent perturbation theory of the system response with respect to the exciting field (field-matter interaction) mentioned earlier in this section. Each field interaction can can change only states of the density matrix which correspond to energy differences within the bandwidth of the laser. In the two state example above, a field interaction would only be able to produce a transition from |A A| to |A B| or |B A|. Transitioning from |A A| to |B B| would require two field interactions, one on the ket side and one on the bra side. Field interactions can come in any allowed order with any time difference between them. If the field interactions are impulsive, the system propagates between field interactions along the zeroth-order Hamiltonian, which for spectroscopic measurements includes system and bath dynamics[171]. During the infinitesimally short field interaction the system propagates along the field-matter Hamiltonian alone and the system propagation during this time can be neglected. Using a delta-function field interaction allows for the dynamics of a system to be described by a Green’s operator[105]. A “diagrammatic” approach can make perturbative density matrix treatments of spectroscopic signals simple to understand. The diagrammatic approach uses so called double sided Feynman diagram to visually represent the density matrix element that 14
describes the system after each excitation of the system via field or vacuum interactions. A brief tutorial on diagrams is depicted in figure 1.6, adapted from [35]. These diagrams can be used to anticipate a number of features in the 2D spectrum, as shown in figure 1.7. From such a diagrammatic approach, it is clear that numerous combinations of signals exist. For 2DES the two combinations of interest are called the “rephasing” and “non-rephasing” signals, which are depicted in figure 1.7. These two signals are two of three pathways which produce signals that have a frequency within the laser bandwidth. The third such pathway is called the “double quantum” signal and arises when two positive frequency interactions happen followed by a negative frequency interaction. As will be demonstrated in the next chapter, the double quantum signal can be isolated from the rephasing or non-rephasing signal via a phenomenon called phase matching.
15
Signal Emission
Ket-side Interaction contributes: +AG, +k
Direction of Time
Field Interaction
G
G
A
G
A
B
A
G
G
G
Bra-side Interaction contributes: -AG, -k
A
G
G
A
G
G
A
A
Ket-side Interaction contributes: -AG, -k
Bra-side Interaction contributes: +AG, +k
G
A
A
G
A
A
A
A
Figure 1.6: A tutorial on the components of double sided Feynman diagrams and how they encode the energy (ω) and momentum (k) contributions to the signal.
16
SE
ESA
G
G
A
G
B
B
A
G
A
G
F
B
A
A
G
G
B
B
G
A
G
B
G
B
G
G
G
G
G
G
G
G
G
G
A
A
A
G
B
G
F
A
A
A
G
G
A
A
A
G
B
G
A
G
G
G
G
G
G
G
Rephasing
Non-Rephasing
GSB
Rephasing
Non-Rephasing
B
B
A
A
A
A
B
B
Figure 1.7: The three main signals observed with 2DES: Excited State Absorption (ESA), Ground State Bleach (GSB), and Stimulated Emission (SE) are depicted with double sided Feynman diagrams color coded to their peak locations in a cartoon 2DES spectrum of the rephasing and non-rephasing signals. Not all signal pathways or peak positions are represented here. For ESA, the anharmonicity of the excited state F is such that ωBG > ωFA > ωAG .
17
The 2D spectrum offers a wealth of information through the peak shape in addition to peak position. At early waiting times, the peak elongation along the diagonal is informative. Inhomogeneous broadening manifests itself as the diagonal linewidth at short waiting times, while homogeneous linewidth at short waiting times is represented by the anti-diagonal linewidth of the peak[105]. As the waiting time increases, one can measure the loss of correlation of the detected frequency from the excitation frequency by observing the anti-diagonal linewidth broaden, producing a round peak shape. Just like transient absorption, population kinetics can be resolved by examining the spectra as a function of waiting time. The advantage of a 2D measurement for determining population kinetics is that transfers between different detected frequencies are straight-forwardly resolved by watching the evolution of the cross-peak associated with the initial and final state. Second order processes, involving two transfer steps between initial and final states, however remain difficult to resolve. Depictions of energy transfer via Feynman diagrams and cartoon spectra, along with how inhomogeneous broadening manifests itself in 2DES is illustrated in figure 1.8.
18
Rephasing B
B
F
B
B
G
G
A
G
B
A
A
A
A
B
B
G
A
G
B
G
G
G
G
A
B B
A
A
B
T=0
T >> 0
D
C B
B
A
A
A
B
A
B
Figure 1.8: A. Rephasing and Feynman diagrams for depicting an energy transfer event over the population time. B. Cartoon rephasing spectrum of an energy transfer event. C. Inhomogeneously broadened lineshape at early waiting times. D. Inhomogeneously broadened lineshape at later waiting times, illustrating the loss of coherence between initial excitation and final detected state.
[1] James Barber. Photosystem ii: the engine of life. Quarterly reviews of biophysics, 36(01):71–89, 2003.
19
[2] Deborah A Berthold, Gerald T Babcock, and Charles F Yocum. A highly resolved, oxygen-evolving photosystem ii preparation from spinach thylakoid membranes: Epr and electron-transport properties. Febs Letters, 134(2):231–234, 1981. [3] Robert E Blankenship. Molecular mechanisms of photosynthesis. John Wiley & Sons, 2008. [4] Tanai Cardona, Arezki Sedoud, Nicholas Cox, and A William Rutherford. Charge separation in photosystem ii: a comparative and evolutionary overview. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1817(1):26–43, 2012. [5] Public Domain. Light-dependent reactions of photosynthesis at the thylakoid membrane. 2012. [6] Brent Donovan, Larry A Walker, Daniel Kaplan, Marcel Bouvier, Charles F Yocum, and Roseanne J Sension. Structure and function in the isolated reaction center complex of photosystem ii. 1. ultrafast fluorescence measurements of psii. The Journal of Physical Chemistry B, 101(26):5232–5238, 1997. [7] Brent Donovan, Larry A Walker, Charles F Yocum, and Roseanne J Sension. Transient absorption studies of the primary charge separation in photosystem ii. The Journal of Physical Chemistry, 100(5):1945–1949, 1996. [8] Juan Du, Takahiro Teramoto, Kazuaki Nakata, Eiji Tokunaga, and Takayoshi Kobayashi. Real-time vibrational dynamics in chlorophyll< i> a studied with a few-cycle pulse laser. Biophysical journal, 101(4):995–1003, 2011. [9] Camiel Eijckelhoff, Henny van Roon, Marie-Louise Groot, Rienk van Grondelle, and Jan P Dekker. Purification and spectroscopic characterization of photosystem ii reaction center complexes isolated with or without triton x-100. Biochemistry, 35(39):12864–12872, 1996. [10] FD Fuller and JP Ogilvie. Continuum probe two-dimensional electronic spectroscopy of the photosystem ii reaction center. In EPJ Web of Conferences, volume 41, page 08018. EDP Sciences, 2013. [11] Andrius Gelzinis, Leonas Valkunas, Franklin D Fuller, Jennifer P Ogilvie, Shaul Mukamel, and Darius Abramavicius. Tight-binding model of the photosystem ii reaction center: application to two-dimensional electronic spectroscopy. New Journal of Physics, 15(7):075013, 2013. [12] Scott R Greenfield, Michael Seibert, Michael R Wasielewski, et al. Wavelength and intensity dependent primary photochemistry of isolated photosystem ii reaction centers at 5 c. Chemical physics, 210(1):279–295, 1996. [13] Scott R Greenfield and Michael R Wasielewski. Excitation energy transfer and charge separation in the isolated photosystem ii reaction center. Photosynthesis research, 48(1-2):83–97, 1996. 20
[14] Marie Louise Groot, Natalia P Pawlowicz, Luuk JGW van Wilderen, Jacques Breton, Ivo HM van Stokkum, and Rienk van Grondelle. Initial electron donor and acceptor in isolated photosystem ii reaction centers identified with femtosecond mid-ir spectroscopy. Proceedings of the National Academy of Sciences of the United States of America, 102(37):13087–13092, 2005. [15] Marie-Louise Groot, EJ Peterman, PJ van Kan, IH van Stokkum, Jan P Dekker, and Rienk van Grondelle. Temperature-dependent triplet and fluorescence quantum yields of the photosystem ii reaction center described in a thermodynamic model. Biophysical journal, 67(1):318, 1994. [16] ALM Haffa, S Lin, JC Williams, BP Bowen, AKW Taguchi, JP Allen, and NW Woodbury. Controlling the pathway of photosynthetic charge separation in bacterial reaction centers. The Journal of Physical Chemistry B, 108(1):4–7, 2004. [17] Gary Hastings, James R Durrant, James Barber, George Porter, and David R Klug. Observation of pheophytin reduction in photosystem two reaction centers using femtosecond transient absorption spectroscopy. Biochemistry, 31(33):7638–7647, 1992. [18] Johan Hattne, Nathaniel Echols, Rosalie Tran, Jan Kern, Richard J Gildea, Aaron S Brewster, Roberto Alonso-Mori, Carina Glöckner, Julia Hellmich, Hartawan Laksmono, et al. Accurate macromolecular structures using minimal measurements from x-ray free-electron lasers. Nature methods, 2014. [19] AR Holzwarth, MG Müller, M Reus, M Nowaczyk, J Sander, and M Rögner. Kinetics and mechanism of electron transfer in intact photosystem ii and in the isolated reaction center: pheophytin is the primary electron acceptor. Proceedings of the National Academy of Sciences, 103(18):6895–6900, 2006. [20] Ryszard Jankowiak. Probing electron-transfer times in photosynthetic reaction centers by hole-burning spectroscopy. The Journal of Physical Chemistry Letters, 3(12):1684–1694, 2012. [21] J Kern, B Loll, C Luneberg, D DiFiore, J Biesiadka, K-D Irrgang, and A Zouni. Purification, characterisation and crystallisation of photosystem ii from thermosynechococcus elongatus cultivated in a new type of photobioreactor. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1706(1):147–157, 2005. [22] David R Klug, Thomas Rech, D Melissa Joseph, James Barber, James R Durrant, and George Porter. Primary processes in isolated photosystem ii reaction centres probed by magic angle transient absorption spectroscopy. Chemical physics, 194(2):433–442, 1995. [23] Kristin LM Lewis and Jennifer P Ogilvie. Probing photosynthetic energy and charge transfer with two-dimensional electronic spectroscopy. The Journal of Physical Chemistry Letters, 3(4):503–510, 2012.
21
[24] Jeffrey A Myers, Kristin LM Lewis, Franklin D Fuller, Patrick F Tekavec, Charles F Yocum, and Jennifer P Ogilvie. Two-dimensional electronic spectroscopy of the d1-d2-cyt b559 photosystem ii reaction center complex. The Journal of Physical Chemistry Letters, 1(19):2774–2780, 2010. [25] Vladimir I Novoderezhkin, Jan P Dekker, and Rienk Van Grondelle. Mixing of exciton and charge-transfer states in photosystem ii reaction centers: modeling of stark spectra with modified redfield theory. Biophysical journal, 93(4):1293–1311, 2007. [26] Erwin JG Peterman, Herbert Van Amerongen, Rienk Van Grondelle, and Jan P Dekker. The nature of the excited state of the reaction center of photosystem ii of green plants: A high-resolution fluorescence spectroscopy study. Proceedings of the National Academy of Sciences, 95(11):6128–6133, 1998. [27] Grzegorz Raszewski and Thomas Renger. Light harvesting in photosystem ii core complexes is limited by the transfer to the trap: can the core complex turn into a photoprotective mode? Journal of the American Chemical Society, 130(13):4431–4446, 2008. [28] Thomas Rech, James R Durrant, D Melissa Joseph, James Barber, George Porter, and David R Klug. Does slow energy transfer limit the observed time constant for radical pair formation in photosystem ii reaction centers? Biochemistry, 33(49):14768–14774, 1994. [29] Thomas Renger and Eberhard Schlodder. Primary photophysical processes in photosystem ii: bridging the gap between crystal structure and optical spectra. ChemPhysChem, 11(6):1141–1153, 2010. [30] M Rogner, Peter J Nixon, and Bruce A Diner. Purification and characterization of photosystem i and photosystem ii core complexes from wild-type and phycocyanindeficient strains of the cyanobacterium synechocystis pcc 6803. Journal of Biological Chemistry, 265(11):6189–6196, 1990. [31] Elisabet Romero, Ivo HM van Stokkum, Vladimir I Novoderezhkin, Jan P Dekker, and Rienk van Grondelle. Two different charge separation pathways in photosystem ii. Biochemistry, 49(20):4300–4307, 2010. [32] Kenneth Sauer, Junko Yano, and Vittal K Yachandra. X-ray spectroscopy of the photosynthetic oxygen-evolving complex. Coordination chemistry reviews, 252(3):318–335, 2008. [33] Joseph J Shiang, Laurie M Yoder, and Roseanne J Sension. Structure and function in the isolated reaction-center complex of photosystem ii. 2. models for energy relaxation and charge separation in a protein matrix. The Journal of Physical Chemistry B, 107(9):2162–2169, 2003. [34] Villy Sundström. Femtobiology. Annu. Rev. Phys. Chem., 59:53–77, 2008. 22
[35] Andrei Tokmakoff. Diagrammatic perturbation theory. MIT Open Course Work. [36] Yasufumi Umena, Keisuke Kawakami, Jian-Ren Shen, and Nobuo Kamiya. Crystal structure of oxygen-evolving photosystem ii at a resolution of 1.9 angstrom. Nature, 473(7345):55–60, 2011. [37] Peter J van Leeuwen, Maaike C Nieveen, Erik Jan van de Meent, Jan P Dekker, and Hans J van Gorkom. Rapid and simple isolation of pure photosystem ii core and reaction center particles from spinach. Photosynthesis research, 28(3):149–153, 1991. [38] Michael R Wasielewski, Douglas G Johnson, Michael Seibert, et al. Determination of the primary charge separation rate in isolated photosystem ii reaction centers with 500-fs time resolution. Proceedings of the National Academy of Sciences, 86(2):524–528, 1989. [39] Thomas John Wydrzynski and Kimiyuki Satoh. Photosystem II: the light-driven water: plastoquinone oxidoreductase, volume 22. Springer, 2006. [40] Laurie M Yoder, Allwyn G Cole, and Roseanne J Sension. Structure and function in the isolated reaction center complex of photosystem ii: energy and charge transfer dynamics and mechanism. Photosynthesis research, 72(2):147–158, 2002.
23
CHAPTER II
EXPERIMENTAL IMPLEMENTATION OF 2DES
In this chapter, we discuss a number of important details surrounding the realization of a 2DES experiment. While a 2DES apparatus was already present when I began my work, I describe a new method for collecting 2DES data that offers a dramatic improvement over our previous pump probe setup and opened new doors for the sort of questions that could be addressed in d1d2 and other systems with our setup. Specifically the signal to noise increased by nearly 20 fold and we increased the bandwidth of our sources so that we could approach 10 fs pulses, allowing the clear resolution of dynamics on the 20-30 fs time scale. In the appendices of this chapter, some information on improvements to our NOPA setup are presented.
2.1
Introduction
2DES was first experimentally demonstrated[57, 55] in 1998, and since then there have been a number of groups demonstrating different methods for collecting 2DES spectra[49, 45, 81, 77, 79, 52, 65, 71, 53, 66, 64, 68, 80, 42]. Each of the various 2DES methods have advantages and disadvantages, as they attempt to solve different problems associated with the realization of 2DES. Some of the common practical problems addressed in the literature are: achieving phase stability, reducing instrument and measurement complexity, lowering acquisition time and minimizing component costs. Most setups use a fully non-collinear box-CARS arrangement in some form or another, since it gives a background free signal, and use motorized delay lines to scan 24
t1 . Compared to fully non-collinear methods, 2DES in the pulse-shaped pump-probe geometry[52, 65, 73] provides a reduction in measurement and instrument complexity. A significant benefit of the pump-probe geometry itself is that the collected signal is automatically real and absorptive, as are transient absorption measurements made in the same geometry[50]. In typical non-collinear geometries, rephasing and non-rephasing signals are collected separately and must be combined together and “phased” to unmix dispersive and absorptive components of the spectra. As will be discussed later, the unmixing process can be challenging and collecting the pathways at separate times leaves one vulnerable to noise from long-term laser fluctuations. In addition to the simplicity of the geometry itself, the use of a pulse-shaper to create the two excitation pulses eliminates uncertainties in coherence time zero that can be present in motorized setups[41] and enables phase-cycling which can aid in scatter removal and the isolation of signals of interest[79, 65, 72]. The use of a pulse-shaper also offers the flexibility to explore the effect of pulse-shape (amplitude and phase) on the measurement[78, 68]. Despite these advantages, pulse-shaped 2DES in the pump-probe geometry suffers from high component costs and offers less freedom to optimize the signal-to-noise ratio (S/N) compared to background-free geometries. In the pump-probe geometry, the signal emerges collinear with the probe and both are sent directly into the camera. As a result, the probe must be made to be weak in order to avoid saturating the detector. As the signal strength is proportional to the probe field, its reduction lowers the achievable S/N. Background free geometries are popular because they do not send the probe directly into the detector, so the S/N may be improved by increasing the probe strength and optimizing the relative strengths of signal and the local oscillator used for heterodyne detection[61]. To exploit this fact, polarization schemes have been employed in the pump probe geometry to control the local oscillator power[84, 65]. Polarization controlled pump-probe geometry 2DES, however, is limited to measuring certain tensor elements of the third order response, and cannot measure the all parallel response. Traditional box-CARS background free 2DES
25
setups have no restriction on what tensor element is measured. We note that the achievable S/N of background-free detection is reduced if phase instability between the signal and local oscillator is significant[70], a problem that can be overcome with the use of diffractive optics[49]. Here we present a hybrid diffractive optic and pulse-shaping based approach to 2DES. It combines the advantages of background-free box-CARS detection with the precise time-delays and phase-cycling capabilities of pulse-shaping. The setup can be readily inter-converted between the background free and pump-probe geometries to fit the demands of the system being studied.
2.2
Review of 2DES setups and Signal Isolation
In what follows we will consider a “hybrid” setup which employs five input pulses, while typical 2DES setups employ three input pulses. To understand how a five pulse setup could give the same information as a three pulse setup, we begin with a review of third order signals. All 2DES setups seek to measure a specific set of signal pathways in the third order response of a material to an applied optical field, as explained in the first chapter. The phrase “third order response” refers to a Volterra expansion of the output signal with respect to its input fields, where the desired signal is then the third term in such an expansion. It is important to note that the Green’s function formalism which underpins the diagrammatic description of 3rd order signals presented in the first chapter is essentially a special case of a Volterra series where the input fields are all delta functions. The linear, or first term in the Volterra series describes the absorption and dispersion of the sample, both of which can be straight-forwardly combined into a complex index of refraction. It can be shown[43] that the second order term in such an expansion is zero for isotropic media, such as the typical liquid (or glassed) sample. Typically, in the measurement of third order signals, the linear response term is ignored in the analysis by presuming that the third order signal that is generated is only weakly affected by the linear response both in its generation and propagation out of the sample. These assumptions are valid for low 26
optical densities (low absorption) and thin samples (low dispersion). In reality, one often does not actually have a low enough optical density or thin enough sample to legitimately ignore these effects, and so the signals one measures are not a pure representation of the third order response, but a distorted one[47, 63, 85]. Nevertheless, exposition is greatly simplified by ignoring the linear propagation effects and the relevant intuition gained by this approximation unaltered. Assuming all parallel polarized light, the third order polarization we are seeking in a 2DES measurement can be written in the frequency domain as:
E1 (ω) = |A (ω)| exp k (ω) ·r + iφ (ω) ˆ∞ ˆ∞ ˆ∞ ˆ∞ Esignal (ωs ) =
χ (ωs , ω1 , ω2 , ω3 )
(2.1)
−∞ −∞ −∞ −∞ 3
× E (ω1 ) E (ω2 ) E (ω3 ) δ (−ω1 + ω2 + ω3 ) ∏ dω j dr j=1
Here δ (−ω1 + ω2 + ω3 ) implies that we are looking for a signal frequency ωs = ω3 + ω2 − ω1 . Other combinations exist, such as the third harmonic signal: ωs = ω1 + ω2 + ω3 or the double quantum signal ωs = ω1 + ω2 − ω3 . The wave-vector k (ω) contribution of the field is complicated if treated rigorously. If we again assume that none of the excitation fields experience any linear propagation effects, approximate each input beam as a single ray with no divergence, and then ignore the frequency dependence of the wave vector, equation 2.1 simplifies to something easily used in practice. In particular, we are interested in describing the signal direction in terms of input beam directions b j and the phase of the output field:
ksignal = −b1 +b2 +b3 φsignal (ω) = −φ1 (ω1 ) + φ2 (ω2 ) + φ3 (ω3 )
27
(2.2) (2.3)
A number of setups employ what is called a “box-CARS” arrangement of the input beams, where each input beam is placed on the corner of a square and then all beams are focused into the sample forming a tetrahedron with the point striking the sample. In this case, it’s convenient to write the input beam vectors b j by unit displacements along each coordinate axis. For the box-CARS arrangement three corners of the square serve as the inputs so that the signal emerges in the fourth corner of the box: (1, 1, 1) − (−1, 1, 1) + (−1, −1, 1) = (1, −1, 1) = ksignal . The pump-probe geometry involves just two input beam directions, IE two corners of the box instead of three, so the signal emerges in the same direction as the beam which only interacted once: (1, 1, 1) − (1, 1, 1) + (−1, 1, 1) = (−1, 1, 1) = ksignal . If more pulses are added to the system and we want to know the phase and beam direction of the resultant third order signals, we still use the above expressions, except now we need to consider all possible combinations of three inputs fields and beam directions. For five input fields it may be tempting to think that fifth order signals would be significant. This is actually a concern regardless of how many pulses are sent in and what matters is only the excitation probability per pulse. A single pulse might produce fifth order signals if it had a high probability to excite the sample. Consequently, we aim for an excitation probability per pulse of
E1+E2
|E2>
E1
|G1>
B
eg+EED
A
eg-EED
|G2>
G
G
0
q
Figure 3.1: The coupling of two two-level systems depicted in panel b results in the excitonic energy level diagram in panel c. In panel a, the electronic dimer model is depicted as an unshifted set of harmonic oscillators (see text).
63
Figure 3.2: Oscillating signal pathways from rephasing and non-rephasing pathways in the Electronic Dimer Model and their corresponding peak locations in a 2D map. In the upper left corner of each peak location, the phase of the oscillations (0 or π) is encoded as a + or - with multiple signs indicating contribution from multiple pathways. Phase signs were taken from [110]. See section 3.1.3 for more discussion on the phase of the waiting time coherences. A label for each peak color is included in the lower right of the peak for easy reference in the text.
64
3.1.2
Displaced Harmonic Oscillator Model
In this section we will discuss a model called the displaced harmonic oscillator (DHO), following the work in [94, 128], to explain what the 2D coherence map would look like if a single vibrational mode is included in the system. The DHO model is also a general model used to describe fluorescence and absorption properties of single chromophores. Adaptations of the model are also used to derive Marcus theory for electron transfer. In this model, the potential energy surface of the ground electronic state is presumed to be harmonic, which means that it is parabolic with eigen states evenly spaced at some specific frequency ω0 . The excited state is also assumed to be harmonic with mode spacing ω0 , but it is shifted along some nuclear coordinate q by a distance d and raised in energy. In the case of zero displacement, we have no coupling of the electronic state to the nuclear states and the transition energy between the ground state is given by the vertical displacement of the excited state. Zero displacement is then equivalent to the TLS used to model each chromophore in the EDM. When the excited state is displaced from the ground state along q the inner product (overlap) between the ground state vibrational wavefunction is no longer unit valued (complete overlap). As a consequence, higher vibrational states of the excited state now have some non-zero transition dipole amplitude, and so a series of peaks shifted by the integer multiples of ω0 from the main peak with decaying amplitude forms in the absorption spectrum. This is called a vibronic progression, and is observed in chlorophyll a with one main side peak shifted from the Qy band, which is usually called the Qy(1,0) band. The main peak of the absorption, relative to a zero displacement oscillator, is also shifted by an amount called the “reorganization energy” denoted in figure 3.3 by the symbol λ . The name “reorganization energy” comes from the idea that the solvent is re-organizing around the excited state (which has shifted along the nuclear coordinate q), and this re-organization of molecules carries an additional energy cost. The reorganization energy is given by λ = 12 mω02 d 2 , which can also be written in terms of the dimensionless Huang-Rhys factor S as λ = S¯hω0 . The displaced harmonic oscillator model can be 65
generalized to multiple dimensions, where each vibrational mode ω j acts along a different coordinate and then the electronic states are coupled to a sum of oscillators. This model is frequently used in fluorescence line narrowing spectroscopy to estimate Huang-Rhys factors for a given vibrational frequency. In those experiments, the integrated intensity of the narrowed lines (measured at very low temperatures) is approximately proportional to an overall S, the Huang-Rhys factor along all reaction coordinates [111]. This data can be used, along with linear absorption and CD measurements to obtain S j the Huang-Rhys factor for the oscillator frequency ω j [106]. To determine where oscillations will occur in the 2D spectrum under this model, we again employ a diagrammatic approach (see Figure). In this case, we see oscillations showing up in 6 different peaks and now both the rephasing and non-rephasing signals exhibit oscillations off and on the diagonal. Since all the data presented in the chapter was collected at 77K and kT is about 0.7 cm−1 K −1 , we can ignore the possibility of the system beginning in a hot vibrational state for any mode greater than ~55 cm−1 . Several interesting features are predicted by this model: • Assuming the number of pathways contributing to a given peak is proportional to the amplitude of the peak, then the rephasing map will be asymmetric about the diagonal, unlike the what is seen in the EDM. • The rephasing spectrum will have a peak located on a diagonal line that displaced twice the coherence energy from the main diagonal (magenta peak in the 2D map). This is a unique feature to vibrational coherence and requires the involvement of a ground-state vibrational mode (though the coherence may be on the excited state manifold). • If the amplitude contributions from all pathways are equal, then the lower diagonal peak in both rephasing and non-rephasing may cancel due to oppositely signed phase contributions. This may produce rephasing maps which look like electronic 66
coherences in that they have no diagonal amplitude, and only cross peak amplitude. In this case, one should look either for the magenta peak for clarification. It may not be possible to distinguish between electronic and vibrational in the case that the transition dipole to the vibrationally excited ground state is weak, however.
a
b
E'
E+
E
E
G'
G
0
eg+½m02d2
{
{
q
{
0
d
Figure 3.3: The displaced harmonic oscillator model is depicted in panel a and a 4-level energy level approximation of the same system is shown in panel b.
67
Figure 3.4: Vibrational coherence signal pathways. On the left, all signal pathways that give rise to coherent oscillations during the waiting time for rephasing and non-rephasing signals. On the right, the corresponding peak locations (color coded) for rephasing and non-rephasing signals. In the upper left corner of each peak location, the phase of the oscillations (0 or π) is encoded as a + or - with multiple signs indicating contribution from multiple pathways. Phase signs were taken from [110]. See section 3.1.3 for more discussion on the phase of the waiting time coherences. A label for each peak color is included in the lower right of the peak for easy reference in the text.
3.1.3
Phase of Waiting Time Coherences
In section 3.1.1 a simple model for the time evolution of the density matrix during a coherence was introduced, and it was shown that at least phenomenologically one could describe the time evolution by a damped sinusoid. This idea can be extended for all three coherences in the Liouville pathways considered in this chapter and then the evolution of a 68
2D lineshape S (ω3 ,t2 , ω1 ) can be calculated as follows:
S (ω3 ,t2 , ω1 ) = An ∑ k
ˆ∞ ˆ∞ dt1 dt3 exp [iωt3 ] exp [iωt1 ]
−∞ −∞
× [±Gk (t3 ) Gk (t2 ) Gk (t1 )] ˆ∞ ˆ∞ Sk (ω3 ,t2 , ω1 ) ≈ ± dt1 dt3 exp [iωt3 ] exp [iωt1 ] −∞ −∞
× exp [−iε3t3 − γ3t3 ] exp [−iε2t2 − γ2t3 ] exp [iε1t1 − γ1t1 ] = ±θ (t2 ) exp [−γ2t2 − iε2t2 ] 1 2π (γ1 + i (ε1 − ω1 )) (γ3 + i (ε3 − ω3 )) 1 = ± θ (t2 ) exp [−γ2t2 − iε2t2 ] 2π 1 × (γ1 + iΔ1 ) (γ3 + iΔ3 ) ×
Here An is a complex number incorporating the excitation fields and the transition dipole strength. The ± is included for the signal type: + for ground state bleaching and stimulated emission, − for excited state absorption. For the real part of the signal we get the following expression for the signal:
Re [Sk (ω3 ,t2 , ω1 )] = 2πAn exp [−γ2t2 ] L cos [|ε2 |t + φ ] −1 2 2 L = (γ1 γ3 − Δ1 Δ3 ) + (Δ1 γ3 + Δ3 γ1 ) φ = −sign (ε2 ) Arg [(γ1 + iΔ1 ) (γ3 + iΔ3 )]
(3.2)
From this we can see that phase of the coherence, interestingly, depends on the initial and final coherences prepared by the first and third excitation pulses. That is to say, the phase will vary continuously for different points in (ω1 , ω3 ) space. To associate a phase with the peak locations we found via diagrammatic modeling we note that at the peak Δ1 = 0 and Δ3 = 0. Thus the phase contribution from φ is zero, and all peak centers should
69
be in phase except for the contribution from An . For transform limited pulses, the only phase contribution we get from An comes from the sign of the transition dipole. In the case of a displaced harmonic oscillator, we can compute the transition dipole analytically by calculating the inner product between a shifted and unshifted harmonic oscillator wavefunction (see equation 3.3), which can be negative for some values of displacement. 2 /2 − (q − Δ)2 /2 H [ j − 1, q] H [k − 1, q] dq exp −q −∞ π 2 j−1 2k−1 ( j − 1)! (k − 1)!
´∞ An =
(3.3)
Here H [ j, q] refers to the Hermite polynomials, and j, k are vibrational quantum numbers. Without knowing the Huang-Rhys factor (strength of the displacement) or, in the case of electronic coupling, the electronic wavefunctions, a general statement about the phase of the peaks cannot be made. It is known, however, that the diagrammatic modeling approach will only predict a 0 or π phase shift between peaks[94]. In [110], some phase guidelines are given for Huang Rhys values < 1 and for simple models of electronic dimer measured with all parallel polarization. Those phases are shown in figures 3.2 and 3.4 for the ED and DHO models respectively. Researchers have reported observing peaks with phases that are not 0 or π, however[109]. The simplest way to explain this with the framework above is spectral overlap. Due to the fact that each peak produces a continuously varying phase throughout the 2D spectrum (see equation 3.2), if two peaks come close together their phases will combine (via sum of angles identity) and produce a value which is neither 0 or π. Quantum transport, where coherences relax to populations or other coherences during their free evolution is another way that the phase may change from 0 or π. In [109], quantum transport was invoked to explain an observed π/2 phase difference between a diagonal peak and cross peak. Finally, it has been demonstrated via simulation [110, 93] that mixing between vibrational and electronic degrees of freedom can also cause non 0 or π phase shifts. 70
3.1.4
Predicting the effect of static disorder
An approximate way to treat static disorder is to say that the electronic degrees of freedom are affected, but the nuclear modes are not. This essentially amounts to a statement that vibrations are oblivious to small changes to electrostatic environment. Of course, this is not completely accurate – one can neatly observe small (no more than 10s of wavenumbers) vibrational shifts in response to dipole shifts of neighboring chromophores[91]. Nevertheless, if this approximation is made we can say that vibrational coherence maps will simply be broadened as though the diagonal position of the electronic transition had been moved. This gives us the picture for vibrational coherence maps shown in figure 3.5. For electronic coherences, static disorder of exciton energies would create maps at different coherence frequencies and if many such frequencies were closely spaced with, for example, a Gaussian amplitude distribution then one would expect significant amplitude reduction due to cancellation unless there is a reason for the phase of the different static states to be correlated to one another in such a way as to produce a constructive interference.
Figure 3.5: A toy model to describe the effect of inhomogeneous broadening on the vibrational coherence map for both rephasing and non-rephasing spectra. As a consequence of this effect, we label the diagonal lines over which the peaks smear for future reference: D+, D, D-, and D–.
71
3.1.5
Literature Review on coherence origin protocols
The models presented in sections 3.1.1 and 3.1.2 are not state of the art anymore, although they were just a couple of years ago[127, 92]. Now the focus has pushed to understand how coherence maps will behave when both excitonic and vibrational coupling effects are mixed in the same system. Naively, one might expect just a linear super-position of the two models, but of course reality is slightly more complicated. In particular, the simple signed phase relationships are no longer present [110, 94]. Adding disorder is not as simple as the model presented in section 3.1.4, but the intuition is reasonably close. Disorder has extremely deleterious effects on the amplitude of electronic coherences and complicates the phase relationships of vibrational coherence peaks beyond all recognition[94]. The peak shapes for coherences in the presence of electronic disorder do elongate as anticipated. A tabulated list of protocols published in the literature is presented in tables 3.1.5 and 3.1.5. The tables runs approximately chronologically and it can be seen that as time progresses, the interpretation of coherence maps has become more and more involved. Figure 3.6 partially reproduces the results presented in [94], which demonstrate peak positions for an electronic dimer coupled to a harmonic ladder of vibrational states for the rephasing signal. Only a limited subset of possible signal pathways are considered (see figure caption). It is clear that the number of signal pathways for a vibronic system explodes, which explains in part why interpretation of even toy reductions of realistic multi-chromophoric systems is so complicated.
72
Citation [113]
[95]
[127]
Discrimination Rules Electronic coherence is characterized by an anti-correlation of peak width in the cross peak to diagonal amplitude. It is argued that this particular behavior would not result from vibrational wave packet motion. Electronic coherence can be identified by a lack of oscillation amplitude on the diagonal in the rephasing concomitant with diagonal amplitude in the non-rephasing. This is in contrast with vibrational coherence which have significant diagonal amplitude oscillations. Presence of amplitude detected at e − ν is an indication of a ground-state vibrational coherence. However, the lack of amplitude here is not necessarily and indication of electronic coherence, as the transition diple strength for these e − ν transitions may be weak. Additionally, the rephasing spectrum shows a peak at (e − ν, e + ν), while this does not appear in the rephasing.
Table 3.1: Review of various protocols found in the literature to discriminate between vibrational and electronic coherences
73
Citation [94]
[110]
Discrimination Rules Continued Simulations reveal that system disorder has a strong effect on the amplitude of the coherence map and electronic coherences are more strongly affected than vibrational. In addition, the phase of the oscillation in the presence of disorder and spectral overlap is so complex that any protocol depending upon observing 180 degree phase shifts for identification is unrealistic. It is shown that when a vibrational frequency is in resonance with an excitonic transition, two main signatures are predicted. 1. Stronger cross-peak amplitude than on the diagonal in rephasing maps, and furthermore the amplitude of the lower cross peaks is greater (asymmetric amplitude across the diagonal). 2. Cross peaks at (e1 , e2 + ν) and (e2 + ν, e1 ). Peak locations similar to [127], but now relative phases are examined with sophisticated modeling to explain phase shifts that are less than 180 degrees. Temperature is suggested as a discriminating agent. Vibrational phase is more sensitive to temperature than electronic degrees of freedom. Of particular note, the authors point out that arbitrary phase relationships between cross-peaks and diagonal peaks are more likely an indication of vibronic coupling than non-secular “quantum transport”,
Table 3.2: Review of various protocols found in the literature to discriminate between vibrational and electronic coherences
74
Figure 3.6: A list of Feynman pathways and peak locations for the rephasing signal of an electronic dimer with a single vibrational state on each exciton. Not shown are the excited state absorptions, phases, or pathways from a second vibrational state.
75
3.2
Coherence Data of d1d2
To study the potential coherent dynamics of d1d2, we collected 185 2DES spectra from 80 fs to 1920 fs in 10 fs steps, with about 90 seconds of integration per spectrum for a total acquisition time of about 5 hours. Excitation probability per reaction center per pulse for each pump and probe pulse was 4%. The sample was not moved during this time, since we found changing scatter environments would ruin our ability to resolve the small coherent effects riding on top of the much larger population dynamics. Photo-damage is a concern when the sample is not moved, which is one of the reasons for employing such a low bleach rate. Furthermore, photodamage is thought to primarily result from destruction of the sample by the presence of radical oxygen, which is created by absorbing the triplet states created by d1d2. Thus, to reduce photo-damage the sample was vacuum degassed immediately prior to use in a bid to reduce dissolved oxygen content. For each 2D spectrum, t1 was scanned to a maximum delay of 300 fs in increments of 1.85 fs, phase locked at 592.4 nm, which meant that we were oversampling the t1 axis by at least a factor of 3. Oversampling is not exactly a good thing, which will be demonstrated in the final chapter, but for those unfamiliar with phase locking, it demonstrates that no information was lost by stepping the t1 delay more than a half optical period. An example 2D spectrum is illustrated in figure 3.7. If the real part of the 2D spectra are collected into a 3D solid with the third dimension being the waiting time, then we can find oscillating peaks by explicitly fitting the entire solid with a functional form. I have chosen to fit the linewidths and frequencies of the main oscillations with a sum of real Lorenzian functions, which have the functional form (in waiting time domain t2 ):
fm (t2 ) =
∑ exp [−βn,1t2] (αm,n,1 cos [βn,2t2] + αm,n,2 sin [βn,2t2])
(3.4)
n
The model also incorporates zero-frequency “oscillators”, which are then essentially 76
165
Detection frequency cm-1x102
100
80
160
60 155
40
20 150 0 145
-20 145
150
155
160
165
Excitation frequency cm-1x102
Figure 3.7: A real phased absorptive 2DES spectrum of d1d2 at 170 fs waiting time delay. This is shown primarily to illustrate the new features we observe with the expanded pump and probe bandwidth (and improved signal to noise). Of particular note is the positive cross-peak at 16000 cm−1 .
fitting the population kinetics with sums of exponential decays. The index m here indicates a linear “pixel” index, where m can refer to any frequency-frequency point in the 2D map. In principle one could fit all M pixels in the spectrum independently, but this dramatically reduces the degrees of freedom for the fitting procedure and does not make sense if we want to see the response of all frequency-frequency points at a given population frequency. Instead we fit all pixels with the same exact waiting time function, allowing only the linear coefficients αm,n to vary from pixel to pixel. The linear coefficients must vary because they encode the phase, which was demonstrated in section 3.1 to vary across the 2D spectrum. The non-linear variables, which encode lifetimes, frequencies, and linewidths are “global” to the entire spectral solid. This sort of fitting is similar in spirit to “multi-way” or “global” kinetic analysis performed on transient absorption[129].
77
Mode cm−1 494 730 854 971
Rel Err % 0.12 0.02 0.04 0.02
Dephasing Time f s 219 723 400 1054
Rel Err % 2.4 2.0 2.5 3.6
Table 3.3: This table shows coherent modes from a ~2 ps 2DES scan of d1d2 that were found to have frequencies > 450 cm−1 . The error shown here is a relative error in units of percent, given by √ ( σν /ν) × 100, where σ is a diagonal element of the estimated covariance matrix and ν is the mean value of the estimated parameter.
3.2.1
High Frequency Modes of d1d2
A few exciton models exist for d1d2 (see section 3.2.2) and most exciton difference frequencies in them are lower than 450 cm−1 , so we will define “high frequency” modes as those which are above 450 cm−1 . These frequencies are very unlikely to be electronic in character or be in resonance with an excitonic difference. Thus we expect these frequencies to be good candidates for having vibrational character. Indeed we find that most coherence maps are in agreement with the vibrational model in section 3.1.2, though two may fit better with the vibronic model in figure 3.6. In the fit of the d1d2 data, the following frequencies, dephasing times, and their estimated errors are shown in table 3.3. To identify the character (vibrational, electronic, or vibronic) of the high frequency modes, we will first look at the rephasing amplitude maps to see if the peak positions match the models presented in 3.1. From figure 3.8, the clearest map is given in panel b for the 730 cm−1 mode. Peaks show up in this map at every single predicted point, although with varying amplitude. The difference in amplitude for peaks with equal pathway contributions is due to differences in transition dipole strength. It should be noted, however, that no correction in these spectra was made for optical density or laser profile, so some additional work would need to be done in order to really measure relative transition dipole strengths from this map. The lower peaks (LP10 and LP20 in figure 3.4) are not visible in panels c and d because the detection bandwidth did not span that region, but otherwise the 854 and 971 cm−1 modes appear to match vibrational amplitude profiles,
78
with a weak peak above the diagonal, some lower diagonal amplitude and high amplitude on the cross peak below the diagonal. The 494 cm−1 mode does not appear to fit the vibrational model very well. The dephasing time fit for this mode was only 219 fs (which is unusually short for a vibrational dephasing time) and there appears to be significant amplitude for this mode detected in the blue around 16500 to 17000 cm−1 (as is the case also for the 854 mode). The long diagonal peak in this map is not well explained by the vibrational model presented in section 3.1, because there should not be any significant oscillation at the upper diagonal peak position DP22. In [110], it is pointed out that some waiting time amplitude may appear at DP22 if a second vibrational quantum level is allowed. To assist in seeing the relation between the models presented in section 3.1, stick “chairs” are super imposed on all panels of figure 3.8, where the corners or end points of the chair correspond to a peak position in the vibrational model. Vibronic character is identified by an appreciable amplitude on the lowest diagonal line, detected at the same energy as CP21, which in relation to the “chair” we call the “ottoman”. It can be seen that both the 494 and 730 cm−1 modes have some small amplitude at the ottoman position, indicating that they may be vibronic in character. Non-rephasing signals for d1d2 are a factor of 2-3 times weaker, so the signal to noise ratio for these maps is very poor compared to the rephasing amplitude maps. Consequently, assignment of coherence origin from these maps alone is not feasible. Non-rephasing amplitude maps for the high frequency modes can be seen in figure 3.9 for reference. Further discussions of coherence maps will focus on rephasing only as a consequence of the low signal to noise found in the non-rephasing maps. To look further support a vibrational or vibronic origin for the high frequency modes, we should look to see if the phase of the oscillations match the appropriate model. Here we are expecting cross peaks below the diagonal (CP21 and LP10) in the rephasing to have the opposite phase of diagonal peak DP11. Due to the presence of a strong excited state absorption, the phase of the upper cross peak CP12 should be opposite that of CP21.
79
3
b 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
16 15.5 15 14.5 14 14
17
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
d 17
17
3
16.5
16.5
-1
-1
Detection frequency in cm x10
3
c Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
3
a 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
16 15.5 15 14.5 14 14
Figure 3.8: Rephasing maps of the 4 high frequency modes of d1d2: a: 494 cm−1 , b: 730 cm−1 , c: 854 cm−1 , d: 971 cm−1 . Black lines have been drawn on the maps to guide the eye (in the shape of a chair, see text). In the case of the 854 and 971 cm−1 modes, the bottom rectangle goes outside the detection axis range. All maps are thresholded so that only regions of the map with at least 33% of the peak amplitude within the viewing window are non-white. In black contours the mean rephasing spectrum is shown to illustrate where in the usual 2D spectrum the oscillations are occuring.
80
3
b 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
15.5 15 14.5
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
d 17 3
16.5
16.5
-1
-1
Detection frequency in cm x10
3
16
14 14
17
c 17 Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
3
a 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
16 15.5 15 14.5 14 14
Figure 3.9: Non-Rephasing maps of the 4 high frequency modes we observe in d1d2: a: 494 cm−1 , b: 730 cm−1 , c: 854 cm−1 , d: 971 cm−1 . The signal to noise of the non-rephasing is significantly lower than the rephasing signal, rendering the interpretation of these maps difficult.
Finally, LP20 is predicted to have a phase opposite of CP21. In figure 3.10, the phase plots for the rephasing only are shown, since the poor signal to noise of the non-rephasing renders the phase difficult to interpret. The phase maps require some careful examination if one is not accustomed to looking at them. As the 730 cm−1 mode has the highest signal to noise and clearest amplitude map, we begin there. The lower cross peak CP21 of the 730 cm−1 map shows a phase transition from green to red, centering about on orange, while the lower diagonal (DP11) is more red, but about the same phase within half a radian. This phase difference is not anticipated by a vibrational model, but is possible in a vibronic model[110]. The phase transitions to deep blue and cool blue for the upper cross peak, which is about π radians out of phase from the lower cross-peak – exactly what is expected for given that the upper cross peak signal is from an excited state absorption. The lower cross peak LP10 is green to yellow, which corresponds to a half radian distance from the lower cross peak at CP21. It’s not clear whether this phase deviation from the expected model should be considered significant (the phase at CP21 and LP10 are expected to be 81
the same). The cross peak at LP20 is clearly pi shifted from CP21, which matches the expectations from the vibrational model. It is interesting to note that while the center of gravity for peak amplitudes do not always line up exactly with the parallel lines, the phase transitions do quite well, which is anticipated in [94]. The 971 cm−1 mode appears to follow the vibrational model well, at least for the peaks we have the detection bandwidth to see. There the upper cross peak is about π out of phase with the lower cross peak. The diagonal appears to be nearly a full radian out of phase from the lower cross peak at the diagonal, again not quite matching the expected π phase shift expected between CP21 and DP11. An interesting feature present in all phase maps is a small peak slightly below the main diagonal peak which has oppositely signed phase from the peak on the diagonal. I assign this feature to an excited state absorption pathway where the spectral shift seen is a reflection of the anharmonicity of the electronic state. It’s possible that the presence of this peak is the reason why the phase difference between DP11 and CP21 are not matching expectations, while the other phase relationships are. Even the 494 cm−1 mode has reasonable phase agreement to the vibrational model (CP12 and LP20 are π shifted from CP21), although DP11 is more phase shifted from CP21 than the other modes. It is noted that 494 cm−1 mode has amplitude in the ottoman position and that the phase here matches expectations from a vibronic model (π phase shifted from CP21).
3.2.2
Low Frequency Modes of d1d2
For modes below 450 cm−1 there are a myriad of excitonic differences in d1d2 which could either resonantly enhance vibrational modes here or even exhibit electronic coherence. A brief review of two excitonic models from literature and vibrational modes from fluorescence line narrowing and resonance Raman are shown in figure 3.11. We find four modes ranging from 88 cm−1 to 339 cm−1 , the frequencies and dephasing times of which are recorded in table 3.4. As for the high frequency modes, we first examine the rephasing amplitude maps,
82
2
16
1
15.5
0
15
-1
14.5
-2
3
16.5
2
16
1
15.5
0
2
16
1
15.5
0
14.5
-2 -3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
d 17
3
16.5
2
16
1
15.5
0
3
16.5
-1
3
-1
14 14
17
17
15
-1
3
c
-3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
Detection frequency in cm x10
14 14
Detection frequency in cm x10
3
16.5
17
b -1
3
Detection frequency in cm x10
17
-1
Detection frequency in cm x10
3
a
15
-1
14.5
-2
14 14
-3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
15
-1
14.5
-2
14 14
-3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
Figure 3.10: Rephasing phase maps of the 4 high frequency modes of d1d2: a: 494 cm−1 , b: 730 cm−1 , c: 854 cm−1 , d: 971 cm−1 . All maps are noise-floor thresholded so that only regions of the map with at least 33% of the peak amplitude within the viewing window are non-white. The color scale is such that −π and π are the same color, while all other values between are uniquely encoded. Since no fitting of time zero was performed (it was selected to be approximately the peak of the non-resonant response), the absolute value of the phase here is meaningless, only the relative value phase of peaks within the map is meaningful (see text).
83
which are summarized in figure 3.12. An immediately distinguishable characteristic of the low frequency maps, particularly for the 246 and 339 cm−1 maps is that there is a significant amplitude on D– detected at the same frequency as CP21. In other words, the amplitude at the ottoman of the chair super-imposed on the maps in figure 3.12 is strong. This feature is not anticipated from the vibrational model and the strength of the amplitude relative to CP21 is much higher than what we see in the high frequency maps. While one could say that the inhomogeneous broadening of LP21 would cause its amplitude to smear upwards in detected energy, the extent to which we see it smear upwards (and not downwards) makes this explanation seem unlikely. A more likely explanation is offered in [94], where they examine an electronic dimer coupled with vibrational modes. In this model, they predict a peak detected at the lower energy of the two dimer exciton states and along the D– line which has opposite phase of the CP21 cross peak. Remarkably, both peak position and phase agree from this model agree excellently with the data here (see figure 3.13 for the phase along D–). In the 339 cm−1 mode, the lower cross peak CP21 is shifted by over 1 radian from the diagonal line, which again is not anticipated in a purely vibrational model, but can be explained by a vibronic model[110]. The presence of any significant amplitude on D–, however, outright excludes the possibility that either the 246 or 339 cm−1 are of pure electronic origin. The 88 cm−1 mode suffers from extreme spectral overlap, so it is not possible to distinguish its character. The 88 cm−1 does exhibit a long dephasing time (nearly 2 ps) so it seems probable that it is not a purely electronic coherence. The identity of the 117 cm−1 mode is less clear, as it shows significant amplitude above the D+ line. The vibronic model of [94] does anticipate a peak on what would be a D++ line (that is a line displaced by 2 vibrational quanta above the diagonal), but the phase they predict is opposite to what we observe. It is proposed in [98] that this mode is a consequence of mixing between a low frequency electronic exciton difference and the vibrational mode at 339 cm−1 . This conclusion was reached by, again, simulating a vibronic dimer but with parameters matching the PD1-PD2
84
Mode cm−1 87.8 117 246 339
Rel Err % 0.13 0.50 0.10 0.06
Dephasing time f s 1919 184 439 441
Rel Err % 4.1 2.0 2.1 1.84
Table 3.4: This table shows coherent modes from a ~2 ps 2DES scan of d1d2 that were found to have frequencies < 450 cm−1 . The error shown here is a relative error in units of percent, given by √ ( σν /ν) × 100, where σ is a diagonal element of the estimated covariance matrix and ν is the mean value of the estimated parameter.
b
Amplitude a.u.
Novoderezhkin Model
c
Gelzinis Model
a
Vibrational Frequencies
100
200
300
400
500 600 2 cm-1
700
800
900 1,000
Figure 3.11: Exciton models of d1d2 from the literature.a: Exciton difference frequencies for the Novoderezhkin[107] model of the PSII RC. b: Exciton difference frequencies for the Gelzinis[100] model of the PSII RC. c: vibrational frequencies for the PSII RC from fluorescence line-narrowing[111] (black) and surface-enhanced resonance Raman[112] (green) experiments. The coherences observed in the 2DES data are super imposed by red-dotted lines over all models and vibrational lines.
Chlorophyll pair. The new observation is that upon performing such a simulation, they saw that vibronic modes produce a number of discrete frequencies in the waiting time response with diminished amplitude, not just a single frequency as would be anticipated by a diagrammatic approach, and that the maps of these secondary peaks do not match expectations for the main resonance.
85
3
b 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
16 15.5 15 14.5 14 14
17
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
d 17 3
17
16.5
16.5
-1
-1
Detection frequency in cm x10
3
c Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
16.5
-1
Detection frequency in cm x10
3
a 17
16 15.5 15 14.5 14 14
14.5
15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
16 15.5 15 14.5 14 14
Figure 3.12: Low frequency rephasing coherence amplitude maps for the following waiting time frequencies: a: 87.8 cm−1 , b: 117 cm−1 , c: 246 cm−1 , and d: 339 cm−1 . Black lines have been drawn on the maps to guide the eye (in the shape of a chair, see text). The amplitude maps are thresholded so that only amplitude responses greater than 33% of the peak value of the viewing window are non-white. In black contours the mean rephasing spectrum is shown to illustrate where in the usual 2D spectrum the oscillations are occuring.
86
0 15.5
-3
16.5
2
16
1
15.5
0
16.5
2
16
1
15.5
0
15
-1
14.5
-2 -3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
d 17
3
16.5
2
16
1
15.5
0
-1
3
-1
17
3
14 14
17
3
14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
Detection frequency in cm x10
3
-2
14.5 14 14
Detection frequency in cm x10
-1
15
17
-1
1 16
Detection frequency in cm x10
3
2
16.5
c
b
17
-1
Detection frequency in cm x10
3
a
15
-1
14.5
-2
14 14
-3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
15
-1
14.5
-2
14 14
-3 14.5 15 15.5 16 16.5 -1 3 Excitation frequency in cm x10
17
Figure 3.13: Rephasing coherence phase maps. Low frequency rephasing coherence phase maps for the following waiting time frequencies: a: 87.8 cm−1 , b: 117 cm−1 , c: 246 cm−1 , and d: 339 cm−1 . The phase maps are thresholded so that only amplitude responses greater than 33% of the peak value of the viewing window are non-white. In black contours the mean rephasing spectrum is shown to illustrate where in the usual 2D spectrum the oscillations are occuring.
87
3.3
Conclusions and Outlook
In this chapter we examined the coherent dynamics of the photosystem II reaction center and came to the conclusion that both vibrational and vibronic features exist. Of note we find that the 246, 339, 494, and 730 cm−1 vibrational modes, which can be seen in fluorescence line narrowing spectra, are coupled to the exciton states of d1d2. The existence of coupling between excitonic states and the vibrational modes of Chlorophyll a is perhaps not surprising, but needed experimental verification. 2DES is an incivie method for identifying such coupling, as the character of each mode can in many cases be qualitatively assigned by matching the coherent amplitude and phase “map” to maps created from simple models. From these results we conclude that non-adiabatic effects (i.e. the unseparable mixing of vibrational and electronic wavefunctions) are important in explaining the dynamics of the PSII RC. The existence of non-adiabatic effects can have a profound impact on the rate of electron transfer in d1d2[90], thus these coherent dynamics do actually speak to the functional purpose of the system. Furthermore, the role of vibronic coupling (or non-adiabatic effects) in energy transfer, a second purpose of the PSII RC and primary purpose of antenna proteins, is gaining support[126, 96]. The question which remains, of course, is how important these vibronic modes are to the function of d1d2. Does the probability of electron transfer significantly increase in the presence of vibronic coupling? Thus far simulations of the PD1 PD2 chlorophyll pair indicate that this is true (see figure3.14, taken from [98]). Panel d of figure 3.14, shows the results of a Hierarchical Equation of Motion (HEOM) simulation of a coupled electronic dimer allowing for charge transfer via a third excited “charge transfer” state, which is coupled to either monomeric excitation. To study the impact of a single vibrational mode, an Ohmic spectral density spectral density containing either a heavily damped or weakly damped 339 cm−1 mode was added, and it was found that an enhancement in the charge separation efficiency is seen when the mode is weakly damped. In further simulations (not shown) it was found that the proximity of the vibrational mode energy to an excitonic 88
difference determines the strength of the charge separation enhancement. If the vibrational mode is far from an excitonic difference, the effect is weak, but near resonance (as is the case for the 339 cm−1 mode shown) it is strong. The relevance of vibronic coupling to charge transfer remains to be demonstrated experimentally. To begin to address this question experimentally we have examined an anion band from 790-820 nanometers, identified first by Wasielewski and co-workers[134]. Initial attempts to resolve coherent oscillation in this band was inconclusive, due to low signal to noise (see figure 3.15). Transient grating measurements are not ideal, however, since our TG measurements of Qy band have lower contrast coherent dynamics when compared to those measured with 2DES, likely because of phase cancellation on integration over excitation frequency. An attractive approach to demonstrating the effect of vibrational or vibronic coherence is to preferentially prepare the system in such a state via coherent control techniques[122] and then examine a product state. A coherent control experiment remains to be attempted, but it may offer improved signal to noise and will be more incisive to uncover whether or not a vibronic coherence has a measurable effect on the anion yield. Finally, additional work should be done to confirm our assignment of the 246 cm−1 , 339 cm−1 , 500 cm−1 , and 730 cm−1 as vibronic modes. As mentioned in section 3.1.5, temperature control can be used to strongly affect the phase relationship of vibronic states and so the vibronic character of these modes could be identified by examining coherence phase maps as a function of temperature.
[87] R Almeida and RA Marcus. Dynamics of electron transfer for a nonsuperexchange coherent mechanism. 2. numerical calculations. Journal of Physical Chemistry, 94(7):2978–2985, 1990. [88] Koji Ando and Hitoshi Sumi. Nonequilibrium oscillatory electron transfer in bacterial photosynthesis. The Journal of Physical Chemistry B, 102(52):10991– 11000, 1998.
89
16
Ȧ2: 250 cm
b
-1
0.9 0.8
-1
Detection frequency cm x10
3
Ȧ2: 130 cm
16
aa
-1
15.5
15.5
0.6 0.5 0.4
15
15
0.3 0.2 0.1
14.5
14.5 14.5
15 15.5-1 3 Excitation frequency cm x10
16 3
Ȧ2: 340 cm
-1
16
14.5
15 15.5 -1 3 Excitation frequency cm x10
0.5
c
Γ = 10 cm−1 Γ = 100 cm−1
16
d
0.4
-1
Detection frequency cm x10
0.7
15
0.3
C’’(Z) cm-1 100 10,000
P+D2PíD 1 population
15.5
0.2
1
0.1 500 1000-1 1500 ω cm
14.5 0.0 14.5
15 15.5 -1 3 Excitation frequency cm x10
16
0
200
400
600 Time fs
800
1000
Figure 3.14: Simulated coherence amplitude maps (filled contours) derived from the simulated real rephasing 2D spectra for the special pair dimer model with a 339 cm-1 vibrational mode, shown at frequencies a: ω2 =130 cm-1, b: ω2 =250 cm-1 and c: ω2 =340 cm-1. To better capture the mixed origin of the 250 cm-1 mode, the map shown here in Fig. 4b includes an equal mixture of contributions from the 339 cm-1 and 251 cm-1 simulations, while the other maps are derived from the 339 cm-1 simulations only. The dashed black lines indicate the diagonal and parallel lines offset from the diagonal by ω2 and -2ω2 . Overlaid open contours show the real rephasing 2D spectrum, averaged over waiting time t2 . d: Simulated population of the charge transfer state for coherent and incoherent cases of the 339 cm-1 vibrational mode. The bath spectral densities for the coherent (black line) and incoherent (red line) 339 cm-1 mode are shown in the inset. Blue lines mark optical excitonic splittings.
90
í
Frequency Integrated TG of d1d2 Anion Band
í í í í í í í í 0
1000
2000
3000
4000
5000
6000
Waiting Time (fs) Figure 3.15: a: The frequency integrated real absorptive transient grating signal for d1d2 in the 790-820 nanometer band from 150 to 900 fs with global kinetic fit in solid line. This feature is an excited state absorption (hence the negative sign) and it grows in with a 1.2 ps time component, corresponding to the time scale of charge separation in the PD1PD2 pathway[120].
91
[89] Koji Ando and Hitoshi Sumi. Nonequilibrium oscillatory electron transfer in bacterial photosynthesis. The Journal of Physical Chemistry B, 102(52):10991– 11000, 1998. [90] S Aubry. A nonadiabatic theory for electron transfer and application to ultrafast catalytic reactions. Journal of Physics: Condensed Matter, 19(25):255204, 2007. [91] Carlos R Baiz and Kevin J Kubarych. Ultrafast vibrational stark-effect spectroscopy: Exploring charge-transfer reactions by directly monitoring the solvation shell response. Journal of the American Chemical Society, 132(37):12784–12785, 2010. [92] Vytautas Butkus, Leonas Valkunas, and Darius Abramavicius. Molecular vibrations-induced quantum beats in two-dimensional electronic spectroscopy. The Journal of chemical physics, 137(4):044513, 2012. [93] Vytautas Butkus, Leonas Valkunas, and Darius Abramavicius. Vibronic phenomena and exciton–vibrational interference in two-dimensional spectra of molecular aggregates. The Journal of Chemical Physics, 140(3):034306, 2014. [94] Vytautas Butkus, Donatas Zigmantas, Darius Abramavicius, and Leonas Valkunas. Distinctive character of electronic and vibrational coherences in disordered molecular aggregates. Chemical Physics Letters, 587:93–98, 2013. [95] Yuan-Chung Cheng and Graham R Fleming. Coherence quantum beats in two-dimensional electronic spectroscopy. The Journal of Physical Chemistry A, 112(18):4254–4260, 2008. [96] Niklas Christensson, Harald F Kauffmann, Tonu Pullerits, and Tomas Mancal. Origin of long-lived coherences in light-harvesting complexes. The Journal of Physical Chemistry B, 116(25):7449–7454, 2012. [97] Gregory S Engel, Tessa R Calhoun, Elizabeth L Read, Tae-Kyu Ahn, Tomáš Mancal, Yuan-Chung Cheng, Robert E Blankenship, and Graham R Fleming. Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems. Nature, 446(7137):782–786, 2007. [98] Franklin D Fuller, Jie Pan, Andrius Gelzinis, Vytautas Butkus, S. Seckin Senlik, Daniel E. Wilcox, Charles F. Yocum, Leonas Valkunas, Darius Abramavicius, and Jennifer P. Ogilvie. Vibronic coherence in oxygen photosynthes. Nature, under review, 2014. [99] John N Gehlen, Massimo Marchi, and David Chandler. Dynamics affecting the primary charge transfer in photosynthesis. Science-AAAS-Weekly Paper Edition-including Guide to Scientific Information, 263(5146):499–501, 1994. [100] Andrius Gelzinis, Leonas Valkunas, Franklin D Fuller, Jennifer P Ogilvie, Shaul Mukamel, and Darius Abramavicius. Tight-binding model of the photosystem ii reaction center: application to two-dimensional electronic spectroscopy. New Journal of Physics, 15(7):075013, 2013. 92
[101] Akihito Ishizaki and Graham R Fleming. Theoretical examination of quantum coherence in a photosynthetic system at physiological temperature. Proceedings of the National Academy of Sciences, 106(41):17255–17260, 2009. [102] John M Jean, Graham R Fleming, and Richard A Friesner. Classical and quantum models of activationless reaction dynamics. Berichte der Bunsengesellschaft für physikalische Chemie, 95(3):253–258, 1991. [103] Venugopal Karunakaran, Ilia Denisov, Stephen G Sligar, and Paul M Champion. Investigation of the low frequency dynamics of heme proteins: Native and mutant cytochrome p450cam and redox partner complexes. The Journal of Physical Chemistry B, 115(18):5665–5677, 2011. [104] RA Marcus and R Almeida. Dynamics of electron transfer for a nonsuperexchange coherent mechanism. 1. Journal of Physical Chemistry, 94(7):2973–2977, 1990. [105] Shaul Mukamel. Principles of nonlinear optical spectroscopy, volume 29. Oxford University Press New York, 1995. [106] Vladimir I Novoderezhkin, Elena G Andrizhiyevskaya, Jan P Dekker, and Rienk van Grondelle. Pathways and timescales of primary charge separation in the photosystem ii reaction center as revealed by a simultaneous fit of time-resolved fluorescence and transient absorption. Biophysical journal, 89(3):1464–1481, 2005. [107] Vladimir I Novoderezhkin, Jan P Dekker, and Rienk Van Grondelle. Mixing of exciton and charge-transfer states in photosystem ii reaction centers: modeling of stark spectra with modified redfield theory. Biophysical journal, 93(4):1293–1311, 2007. [108] Vladimir I Novoderezhkin, Andrey G Yakovlev, Rienk van Grondelle, and Vladimir A Shuvalov. Coherent nuclear and electronic dynamics in primary charge separation in photosynthetic reaction centers: a redfield theory approach. The Journal of Physical Chemistry B, 108(22):7445–7457, 2004. [109] Gitt Panitchayangkoon, Dmitri V Voronine, Darius Abramavicius, Justin R Caram, Nicholas HC Lewis, Shaul Mukamel, and Gregory S Engel. Direct evidence of quantum transport in photosynthetic light-harvesting complexes. Proceedings of the National Academy of Sciences, 108(52):20908–20912, 2011. [110] Vaclav Perlik, Craig N Lincoln, Frantisek Sanda, and Jurgen Hauer. Distinguishing electronic and vibronic coherence in 2d spectra by their temperature dependence. The Journal of Physical Chemistry Letters, 2013. [111] Erwin JG Peterman, Herbert Van Amerongen, Rienk Van Grondelle, and Jan P Dekker. The nature of the excited state of the reaction center of photosystem ii of green plants: A high-resolution fluorescence spectroscopy study. Proceedings of the National Academy of Sciences, 95(11):6128–6133, 1998.
93
[112] Rafael Picorel, George Chumanov, Elena Torrado, Therese M Cotton, and Michael Seibert. Surface-enhanced resonance raman scattering spectroscopy of plant photosystem ii reaction centers excited on the red-edge of the q y band. The Journal of Physical Chemistry B, 102(15):2609–2613, 1998. [113] Andrei V Pisliakov, Tomáš Manˇcal, and Graham R Fleming. Two-dimensional optical three-pulse photon echo spectroscopy. ii. signatures of coherent electronic motion and exciton population transfer in dimer two-dimensional spectra. The Journal of chemical physics, 124(23):234505, 2006. [114] Martin B Plenio and Susana F Huelga. Dephasing-assisted transport: quantum networks and biomolecules. New Journal of Physics, 10(11):113019, 2008. [115] Patrick Rebentrost, Masoud Mohseni, and Alan Aspuru-Guzik. Role of quantum coherence and environmental fluctuations in chromophoric energy transport. The Journal of Physical Chemistry B, 113(29):9942–9947, 2009. [116] Patrick Rebentrost, Masoud Mohseni, Ivan Kassal, Seth Lloyd, and Alan AspuruGuzik. Environment-assisted quantum transport. New Journal of Physics, 11(3):033003, 2009. [117] Nicolas Renaud, Daniel Powell, Mahdi Zarea, Bijan Movaghar, Michael Wasielewski, and Ratner Mark. Quantum interferences and electron transfer in photosystem i. 117(29):5899–5908, 2013. [118] Gernot Renger. Mechanism of light induced water splitting in photosystem ii of oxygen evolving photosynthetic organisms. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1817(8):1164–1176, 2012. [119] Christian Rischel, Diane Spiedel, Justin P Ridge, Michael R Jones, Jacques Breton, Jean-Christophe Lambry, Jean-Louis Martin, and Marten H Vos. Low frequency vibrational modes in proteins: Changes induced by point-mutations in the protein-cofactor matrix of bacterial reaction centers. Proceedings of the National Academy of Sciences, 95(21):12306–12311, 1998. [120] Elisabet Romero, Ivo HM van Stokkum, Vladimir I Novoderezhkin, Jan P Dekker, and Rienk van Grondelle. Two different charge separation pathways in photosystem ii. Biochemistry, 49(20):4300–4307, 2010. [121] Gregory D Scholes. Biophysics: Green quantum computers. Nature Physics, 6(6):402–403, 2010. [122] Yaron Silberberg. Quantum coherent control for nonlinear spectroscopy and microscopy. Annual review of physical chemistry, 60:277–292, 2009. [123] S Spörlein, W Zinth, and J Wachtveitl. Vibrational coherence in photosynthetic reaction centers observed in the bacteriochlorophyll anion band. The Journal of Physical Chemistry B, 102(38):7492–7496, 1998. 94
[124] Robert J Stanley and Steven G Boxer. Oscillations in the spontaneous fluorescence from photosynthetic reaction centers. The Journal of Physical Chemistry, 99(3):859–863, 1995. [125] AM Streltsov, TJ Aartsma, AJ Hoff, and VA Shuvalov. Oscillations within the b< sub> l absorption band of< i> rhodobacter sphaeroides reaction centers upon 30 femtosecond excitation at 865 nm. Chemical physics letters, 266(3):347–352, 1997. [126] Vivek Tiwari, William K Peters, and David M Jonas. Electronic resonance with anticorrelated pigment vibrations drives photosynthetic energy transfer outside the adiabatic framework. Proceedings of the National Academy of Sciences, 110(4):1203–1208, 2013. [127] Daniel B Turner, Rayomond Dinshaw, Kyung-Koo Lee, Michael Belsley, Krystyna E Wilk, Paul MG Curmi, and Gregory Scholes. Quantitative investigations of quantum coherence for a light-harvesting protein at conditions simulating photosynthesis. Phys. Chem. Chem. Phys., 2012. [128] Daniel B Turner, Krystyna E Wilk, Paul MG Curmi, and Gregory D Scholes. Comparison of electronic and vibrational coherence measured by two-dimensional electronic spectroscopy. The Journal of Physical Chemistry Letters, 2(15):1904– 1911, 2011. [129] Ivo HM van Stokkum, Delmar S Larsen, and Rienk van Grondelle. Global and target analysis of time-resolved spectra. Biochimica et Biophysica Acta (BBA)-Bioenergetics, 1657(2):82–104, 2004. [130] Marten H Vos, Michael R Jones, C Neil Hunter, Jacques Breton, and Jean-Louis Martin. Coherent nuclear dynamics at room temperature in bacterial reaction centers. Proceedings of the National Academy of Sciences, 91(26):12701–12705, 1994. [131] Marten H Vos, Jean-Christophe Lambry, Steven J Robles, Douglas C Youvan, Jacques Breton, and Jean-Louis Martin. Direct observation of vibrational coherence in bacterial reaction centers using femtosecond absorption spectroscopy. Proceedings of the National Academy of Sciences, 88(20):8885–8889, 1991. [132] Marten H Vos, Fabrice Rappaport, Jean-Christophe Lambry, Jacques Breton, and Jean-Louis Martin. Visualization of coherent nuclear motion in a membrane protein by femtosecond spectroscopy. Nature, 363:320–325, 1993. [133] Haiyu Wang, Su Lin, James P Allen, JoAnn C Williams, Sean Blankert, Christa Laser, and Neal W Woodbury. Protein dynamics control the kinetics of initial electron transfer in photosynthesis. Science, 316(5825):747–750, 2007. [134] Michael R Wasielewski, Douglas G Johnson, Christopher Preston, Michael Seibert, et al. Determination of the primary charge separation rate in photosystem ii reaction centers at 15 k. Photosynthesis research, 22(1):89–99, 1989. 95
CHAPTER IV
POST ANALYSIS OF 2DES DATA
4.1
Introduction
To extract parameters characterizing the dependence of a 2D spectrum with population time, a straight-ahead approach is to use non-linear least-squares (NLLS) fitting. This chapter details how the coherence maps in chapter 3 were generated using a NLLS approach. At first glance it may seem strange to dedicate an entire chapter to the methodology of non-linear fitting, when numerous packages for such task can be found in Matlab, Python, C++ or other computer languages. What those packages actually implement are methods to minimize or maximize objective functions, and if your objective function is simple and well behaved, then there is little need to discuss the issue. However, the problem of fitting a sum of Lorenzians to a 3D stack of data turns out to be a rather ill-posed problem. The objective function is rugged even for infinite precision much less double precision, due to the strong coupling between of phase and frequency. The problem size if naively implemented is enormous, as the entire 2DES data tensor consumes nearly 2 gigabytes. Performing operations which require hundreds or thousands of manipulations on data of this size is not only prohibitive on desktop computers, but also unnecessary as we will show. What this chapter details are: • Methods to improve the numerical condition of the objective function through a method called variable projection. • A means to analytically compute the gradient, Jacobian, and Hessian of the 96
variable projection functional, which are necessary for analyzing goodness of fit and obtaining rapid convergence to a solution when using standard minimization packages. • Ways to dramatically reduce the problem size so that that hundreds of data fits can be generated in seconds, rather than a single fit in hours. • A strategy for fitting data that gives a rigorous means to selecting the statistically significant number of parameters that should be used in the model.
4.2
The Variable Projection Functional
Given any arbitrary functional form of the model f θ with a N vector of parameters θ = {θ1 , · · · θi , · · · θN }, an objective function of the form (4.1) may be used to estimate θ
O =
2 f θ − y l l ∑ L
(4.1)
l=1
θ = argmin [O] θ Where the data at a time index l ∈ [1, L] is given by yl , and the model at that time is given by fl θ . For three-dimensional data, such as 2DES, we can write:
O =
J,K,L
∑
2 f j,k,l θ j,k − y j,k,l
j,k,l=1
=
M,L
∑
2 fm,l θm − ym,l
m,l
m = 1, · · · , JK Where m is used to represent the linearization of the matrix of j, k pairs. This linearization of multi-dimensional indices is extensible beyond three dimensions, and would correspond to a specific unfolding of a tensor of d1 , d2 , · · · , dn tuples. Once the 97
objective function is specified, one is ostensibly done and the estimate θˆ can be found using one of the myriad of non-linear optimizers that exist. For non-linear optimization, I use NLopt [140], which gives one access to a large a library of algorithms that may be used to employ in the above search for a global optimum of the objective function. In particular, I use a multi-start algorithm called Multi-Level Single Leakage employing a Low Discrepancy Sequence in lieu of random starting points [142], which controls a local optimizer. The local optimizer is either a derivative-free algorithm BOBYQA[146] (requiring only the objective function) or a low-storage BFGS algorithm [143] which requires that an estimate of the gradient be supplied. One can quickly find, however, that simply using a non-linear fitting package on a sum of Lorenzians model with a data set of the size found in fitting 2D spectra will be disastrous. As posed, the fitting procedure can be highly unstable due to strong coupling between linear and non-linear terms in the model. A method called “variable projection”, seminally demonstrated in [138], alleviates the numerical instability of the naive fitting problem by separating the parameter vector θ into parameters which enter linearly in the model functional and those which enter non-linearly, such that we have θm = αm , βm . Under this separation, the model function fl,m can be written as fl,m = Φm βm αm . The fitting process under variable projection is broken into two steps: a set of non-linear parameters βm are guessed, and then αm are estimated via linear regression. An important simplification can be made if we can assume that the nonlinear parameters βm are independent of m, i.e. βm = β . In a sum of Lorenzians model, this simplification implies that the resonant frequencies observed along the waiting time are independent of excitation and detection frequency. In our experiment, we cannot feasibly sample the waiting time with a sampling period much less than 10 fs, so the maximum resolvable waiting time resonances will be 50 terahertz (1600 cm−1 ) or less. In other words, we will only be resolving IR resonances along the waiting time coordinate axis. Under an adiabatic approximation, we can assume that the vibrational resonances are independent of the electronic excitation or detection frequency. Thus, in
98
what follows we will make the simplification βm = β . The amplitudes and phases of the oscillators, which are encoded in α are allowed to vary from pixel to pixel. Assuming the linear regression estimation of αm remains stable, the question that remains is how we will vary β such that O is minimized. The precise implementation details of how to traverse a multi-dimensional objective function are beyond the scope of this thesis, but any method apart from brute-force grid-searching or genetic-algorithms is going to employ either an analytic or numerically estimated gradient or Jacobian of the objective function. Convergence towards an optimal solution can proceed more rapidly if an accurate analytic gradient is used. In addition, understanding the validity of fit parameters can be efficiently estimated in terms of the Hessian, which as will be shown can be approximated from the Jacobian. By separating the parameter space into linear and non-linear parameters, we are making the statement that the linear parameters depend upon the non-linear parameters: θm = αm β , β . So, in order to determine an analytic gradient of the objective function with respect to β , the parametric dependence of αm β needs to be known. While the normal equations immediately give an equation for α β −1 Φ β ym,l = α (β ) , the inverse of Φ Φ can be numerically Φ β Φ β unstable[137]. Thus many authors who have developed the variable projection method pursue more stable solutions, which typically involve either QR[141] or SVD[144] factorization of Φ, with QR being the faster if somewhat less reliable[139]. Both the SVD and QR decomposition methods can be used to construct an estimate of the pseudo-inverse. The pseudo-inverse Φ+ of Φ, can be constructed for any full rank rectangular matrix, which upon left multiplication with the observation vector y results in the minimum norm solution for α β . More generally, Φ+ falls into a class of matrices called general matrix inverses which satisfy (at minimum) the properties ΦΦ− Φ = Φ and Φ− ΦΦ− = Φ− . The pseudo-inverse also satisfies (Φ− Φ) = Φ− Φ and (ΦΦ− ) = ΦΦ− . †
†
Equipped with a suitable general inverse, the objective function may be written as: 99
2 O = ∑m |Φm Φ− mym −ym |2 , where |·|2 denotes the usual L2 norm. The notation may be 2
further compactified by introducing the projector operator Pm = Φm Φ− m and its orthogonal ⊥ 2 , so that our problem may be represented as O = complement Pm⊥ = I − Φm Φ− ∑ m Pm ym 2 . m Our task to compute the gradient and Hessian of the objective function must therefore involve differentiating the projector.
O =
L,M
∑
2 − fl,m αm , β + yl,m
l,m=1 L,M ∂f ∂O l,m = − ∑ 2 − fl,m αm , β + yl,m ∂ βn ∂ βn l,m=1
= −
L,M
∑
2rl,m
l,m=1
= − ∂ 2O ∂ βn ∂ βn
∂ fl,m ∂ βn
∂
M
∑ 2rm ∂ βn (Pmym)
(4.2)
m=1 L,M
∂f ∂ l,m − fl,m αm , β + yl,m = −2 ∑ ∂ β ∂ β n n l,m=1 =
L,M
∑
l,m=1
= 2
2
L,M ∂ 2 fl,m ∂ fl,m ∂ fl,m − ∑ 2rl,m ∂ βn ∂ βn l,m=1 ∂ βn,x ∂ βn
M ∂ ∂2 ∂ (P y ) (P y ) − 2 r m m m m ∑ ∑ m ∂ βn∂ βn (Pmym) ∂ βn m=1 ∂ βn m=1 M
(4.3)
Equation (4.2) gives elements of the gradient in terms of the derivatives of the projector times the observation vector ym , and equation (4.3) gives elements of the Hessian. Golub and Pereyra in [138] derived the derivative of the projector Pm using assuming the two aforementioned properties of general matrix inverses, (ΦΦ− ) = ΦΦ− , and the product †
rule: DPm =
Pm⊥ DΦΦ− +
⊥ − Pm DΦΦ
(4.4)
The use of generalized inverse matrices (Φ− ) rather than the specific case of the
100
pseudo-inverse is convenient, as it allows us to employ low-rank approximations available via SVD or QR decomposition, in case the columns of the model matrix Φ become nearly collinear. For column pivoted QR decomposition, a rank k general left inverse is given by: ⎡ ⎢ Φ− = Z ⎣
R−1 (1 : k, 1 : k) Q (:, 1 : k) |
⎤ k
⎥ ⎦
0|N−k Where Z permutes the rows of Φ− according to the pivoting scheme employed in the QR algorithm, and the augmentation of N − k zero rows occurs when k < N. This particular general inverse shares the same properties as the Moore-Penrose pseudo inverse when the rank of Φ is full (they are the same to machine precision), but no longer satisfies the condition that Φ− Φ is hermitian when the rank of Φ is not full. The condition (ΦΦ− ) = ΦΦ− for this general inverse is satisfied for any rank k, thus the †
projector derivative of Golub and Pereyra for the QR-based general inverse is valid. The main reasons for selecting this particular general inverse are that it is fast to compute and coincides with the solution returned by Matlab’s backslash operator. Using the pseudo-inverse (which is often most easily constructed from an SVD decomposition of Φ) is equally valid, though slower. The Jacobian of the model function Pmym = fm is given by D (Pmym ) = (DPm )ym , under the assumption that ym has no dependence on β . Thus, from equation (4.4) we have:
101
Jm =
∂ fm ∂ fm ∂ fm ,··· , ,··· ∂ β1 ∂ βi ∂ βN
∂ fm ∂ βi
Jm,i =
∂ Φm ⊥ Φ Pm ym ∂ βi − ∂ Φm − ∂ Φm αm + Φ I − ΦΦ rm ∂ βi ∂ βi ∂ Φm ∂ Φm − αm − Φ Φ αm ∂ βi ∂ βi − ∂ Φm Φ rm ∂ βi ∂ Φm − Pm⊥ Φ ym + ∂ βi
= = = +
−
(4.5)
The full Jacobian over all m is given by the concatenation of J1···m . The ordering of the parenthesis is important here to ensure that only matrix-vector multiplications are made. =y P⊥ with the When computing the gradient, which is the dot product of the residualrm m
Jacobian, the second term of (4.5) can be dropped, due to the fact that Pm⊥ [Φ− ] = 0[138], giving:
∇βi O = −
M
∑
rm
m=1
∂ Φm ∂ Φm − αm − rm Φ Φ αm ∂ βi ∂ βi
(4.6)
The formulation above for both Jacobian and gradient allows for a potentially full rank 3 tensor ∂ Φm /∂ β = DΦm . To be explicit, DΦm is a three dimensional array, which when cut along the third dimension gives “faces” that are matrices whose columns are the derivative of each column of Φ with respect to the jth nonlinear variable. The jth face is denoted [DΦm ] j : [DΦm ] j =
∂ Φm (1, ·) ∂ Φm (i, ·) ∂ Φm (A, ·) ,··· , ,··· , ∂βj ∂βj ∂βj
In total, there will be A = length [α ] columns in the jth tensor face, forming a L × A matrix. 102
Often, DΦm is very sparse. In the case of a sum of Lorenzians model, each face contains only two non-zero columns. To avoid calculations of the Jacobian or gradient scaling poorly, the sparsity of DΦm should be exploited. A practical means of using the potential sparsity of DΦm is to store a pattern matrix S of size N × Q, where Q is the total number of non-zero columns in DΦm , which encodes the structure of the sparsity. An example structure for a sum of Lorenzians model containing 3 oscillators (6 nonlinear parameters) is shown below: ⎡
⎤
⎢ 1 0 1 ⎢ ⎢ 1 0 1 ⎢ ⎢ ⎢ ⎢ 1 0 1 S = ⎢ ⎢ ⎢ 1 0 1 ⎢ ⎢ ⎢ 1 0 1 ⎢ ⎣ 1 0 1
⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
Once the sparsity structure is known, only the non-zero derivatives of Φ need be stored in columns of a L × Q matrix Φ . The tensor product DΦmαm can then be efficiently implemented using S and Φ :
DΦmαm = Φ (S α )
(4.7)
Where indicates a column-wise scalar multiplication, implemented in Matlab by the built-in function bsxfun. When S is scaled column-wise by α , then the columns of the L × N matrix Φ (S α ) are the sum of the scaled derivatives of all columns in Φ with respect to a given non-linear parameter. This is equivalent to what the tensor product would return: a series of L × 1 column vectors, one for each face of the tensor (and there are N faces in the tensor). Equation (4.5) for the Jacobian contains a different tensor product in its second term, 103
which needs to be treated separately. In (DΦm )rm , each face of the tensor is returning a A × 1 column vector, so the result of the tensor-vector product is a A × N matrix. In this case, the resulting rows of the product retain the sparsity of each tensor face, so it is only necessary to compute those elements.
4.3
Using 2D model functions to reduce the problem size
So far, we have discussed a simplified picture where we’ve assumed that the data for each frequency-frequency “pixel” in a 2D spectrum is given by a sum of Lorenzians with linear scaling factors that are completely independent from pixel to pixel. In other-words, we have only considered a 1-dimensional non-linear problem, spread over multiple “ways” or pixels. It can be advantageous for model parsimony to fit a model function which contains nonlinear parameters in both data dimensions. This sort of model also applies to standard transient absorption data, which when time resolved is also a two-dimensional function. In general a function of the form f (x, y, βx , βy ) could be used to describe two-dimensional data with some list of parameters βx to describe the x-dimension and βy to describe the y-dimension. We consider here a special case where the function f (x, y, βx , βy ) = Φx (x, βx ) Φy (y, βy ), i.e. it is separable. This assumption, while obviously restrictive, will make the math considerably simpler. The separable case can be considered to be two separate non-linear problems, where one seeks to find an optimal non-linear basis for each dimension. Discretizing these functions to fit discrete data results in the following matrix equation:
104
⎡
hi (x, βx ) αi, j g j (y, βy )
i, j
⎢ ⎢ ⎢ ⎢ ⎢ = ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ × ⎢ ⎢ ⎣
Φx,1 (x0 , βx ) .. . Φx,1 x j , βx .. .
···
···
Φx,MX (x0 , βx ) .. . Φx,MX x j , βx .. .
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
Φx,1 (xNX , βx ) · · · Φx,MX (xNX , βx ) ⎤ α1,1 · · · α1,My ⎥ ⎥ .. .. .. ⎥ . . . ⎥ ⎦ αMX ,1 · · · αMx ,My
⎡
Φ (y , β ) · · · Φy,My (y0 , βy ) ⎢ y,1 0 y ⎢ .. .. ⎢ . . ⎢ ⎢ × ⎢ ⎢ Φy,1 (ys , βy ) · · · Φy,My (ys , βy ) ⎢ .. .. ⎢ . . ⎢ ⎣ Φy,1 yNy , βy · · · Φy,My yNy , βy
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
Yˆ = Φx AΦ y The data fitting procedure now involves two projectors:
Px = Φx Φ− x xP
= Φ− x Φx
Py = Φy Φ− y yP
= Φ− y Φy
ˆ 2 . |·| here denotes the Frobenius norm, which Our problem is O = min Y − Y f f is the sum of the square root of the absolute value for each matrix element. For fixed nonlinear variables, we know the solution to the linear problem is given by:
105
y − Φx AΦy = 0 y = Φx AΦy − Φ− x yΦy =
x PAPy
= Aˆ So then we get this solution that no longer contains the linear terms A:
2 − Y Φ Φ O = min Y − Φx Φ− y x y f 2 = min Y − PxYy P f R,C 2 = min ∑ [Y − PxYy P]r,c r,c
Where ∑R,C r,c [·]r,c is meant to indicate that the matrix argument is to be summed over all rows and columns. The derivative of the objective function with respect to elements of βx or βy are:
∂O = ∂ βx,i
R,C
∑ r,c
∂ Px 2 (Y − PxYy P) [Yy P] ∂ βx,i
r,c
∂ Px [Yy P] ∂ βx,i ∂ yP = −2R [PxY ] ∂ βy,i = −2R
∂O ∂ βy,i
Where R is the matrix residual. We need to supply the derivatives of the projectors Px and y P. From [138], lemma 4.1:
⊥ − + P DΦ Φ DPx = Px⊥ DΦx Φ− x x x x
106
The derivative of the row-space projector y P is also found in [138] as a corollary to lemma 4.1 found there:
Dy P =
⊥ Φ+ y DΦy y P
+ ⊥ + Φy DΦy y P
The requirement that one use the pseudo-inverse rather than a more general inverse is + + = A A = A+ A, requires because the beginning of the derivation: A A that the transpose operator · commute with pseudo-inversion, which it does. This may not be true in the general inverse case. In any event, we can easily use the pseudo-inverse here for the derivative of both projectors. It carries a factor of 3-4 greater computational cost to compute the pseudo-inverse rather than the QR-based general inverse discussed earlier, but the scaling is the same. So then we have: ∂O + ⊥ + ⊥ = −2R [PxY ] Φy DΦy y P + Φy DΦy y P ∂ βy,i ∂O ⊥ + ⊥ + [Yy P] = −2R Px DΦx Φx + Px DΦx Φx ∂ βx,i
4.4
Using Orthogonal Decomposition to Reduce Problem Size
In section 4.3, we illustrated a strategy to reduce the number of linear coefficients to be estimated by encoding the “pixel” dimension with some continuous functional form that makes sense, Gaussian for instance. While this strategy can give more physically relevant parameters, the physicality of the parameters one estimates really depends strongly on how valid the functional form being fit to the data was in the first place. If, for instance, excitons in the model are not well described by a single Gaussian function, then one still has multiple “pixels”, which are now Gaussian functions rather than delta functions in frequency-frequency space, describing the same excitonic dynamics. This fact is what 107
gives rise to the need for fitting multiple exponentials in the first place. If we knew the true excitonic basis that would diagonalize the matrix exponential describing the population dynamics, then we could fit the population time data with single exponentials. Assuming, of course, that the process one is looking at is well described by first order differential equations. Given the potential difficulty of finding the true excitonic basis, we can avoid non-linear fitting all together by simply seeking some basis which offers more parsimonious representation than the delta function basis. Discrete Gram polynomials, also known as Chebyshev polynomials, are one such example basis and were chosen because the smoothness of the data makes polynomial interpolation attractive and more importantly because existing software packages are available which can compute these polynomials to double precision for large numbers of basis functions[145]. We find that a high resolution data set containing ~1000x500 frequency-frequency “pixels” can be represented with around 0.1% error using only 30x30 Gram polynomial “pixels”. Thus the problem size can be reduced from fitting 500,000 independent pixels to just 900 pixels. As the Gram basis is orthonormal, if the number of Gram polynomials is equal to the number of data points, we are evaluating essentially finding the Fourier spectrum of the data, in the Gram basis, rather than the Fourier basis. So reducing the number of gram pixels to something less than the number of data points is equivalent to filtering the data in the Fourier domain with a rectangular apodization filter. The legitimacy of such a filter is related to how quickly the basis functions converge to represent the data. We find that for frequency windowed signals, which one has if they are fitting a region of interest in the 2D spectrum, the Gram polynomials converge more rapidly than the Fourier basis. The fact that Gram polynomials have no physical interpretation is irrelevant, we know that the dynamics of these “Gram” pixels will contain a super-position of excitonic states and their dynamics, so we must fit a multi-exponential to the decay of each Gram pixel. The orthogonality of the Gram polynomial basis is important, since it ensures that (unlike a Gaussian decomposition of pixel space), the linear coefficients of neighboring Gram pixels
108
will not be coupled. The problem presented in section 4.2 can be re-written in the Gram G basis as follows:
Y˘ = GY O = ∑ f β α˘ m − Y˘m f β α˘ m − Y˘m m
βˆ = argmin [O] β ,α˘
αˆ = Gα˘ ˆ Ym = f β αˆ Where now the index m refers to the “Gram” pixels, rather than frequency-frequency pixels. The inverse of G = G due to the orthogonality of the Gram polynomials. In English, this means one decomposes the data Y into the Gram basis and then minimizes the difference between the model function f that is independent of the pixel basis. After the data is fit in the Gram basis, the fit is then re-expanded into the original pixel basis. Since the parameters of the Gram basis are fixed, they never enter the gradient of the variable projection problem, which allows us to avoid the dual projector issue in section 4.3. Thus, the only reason to use the two-dimensional variable projection problem is when one believes that the functional form one can fit to the frequency-frequency pixel axis will result in significant decoupling of the population time dynamics, i.e. that each Gaussian “compartment” will exhibit a different rate or frequency than the others. If one is not sure of this, then the approach describe in this section is significantly easier to implement while still offering a great degree of model parsimony over the naive pixel-basis approach.
109
4.5
Goodness of Fit: Taylor Expansion of Objective Function
After running an optimization, it is necessary to know to what degree the parameters found are “valid”. Validity can be measured in a number of ways. One can ask for the sensitivity of the objective function to small changes in the fit parameters. Another metric is to test how well the model predicts values withheld from the data set, something called validation. While the second metric is perhaps the most intuitive and robust, it is computationally intensive to compute. We can estimate the sensitivity of the objective function by expanding the objective function in the vicinity of the optimal solution S using the gradient and Hessian we derived in section 4.2. N ∂O O β ≈ O S + ∑ S (βn − Sn ) n=1 ∂ βn 1 N N ∂ 2O + (βn − Sn ) · (βn − Sn ) ∑∑ 2 n=1 ∂ βn ∂ βn S n =1
(4.8)
An intuitive notion of the sensitivity of the objective function to changes in parameter is the difference: O β − O S = ΔO, where S is the solution vector found via minimization. Under the second order approximation (4.8), and the assumption that the optimal solution lies at a stationary point of the objective function (such that ∂ O/∂ βn = 0), we have: ∂ 2O 1 N N ΔO ≈ S (βn − Sn ) · (βn − Sn ) ∑∑ 2 n=1 ∂ β ∂ β n n n =1 1 ≈ HO S Δβ Δβ 2 1 = Δβ HO S Δβ 2
110
−1 1 u Δβ = Δβ HO S 2 −1 1 −1 1 Δβ HO S Δβ = 2 u −1 1 1 HO S Δβ Δβ = 2 u If the objective function has found a deep minimum, we expect the difference ΔO to require a small displacement β − S = Δβ to achieve a unit of range, while for a shallow displacement it will take a large displacement to achieve a unit of range. Thus, solving ΔO2 = u for Δβ will yield an estimate of the error, where u is the unit of range. At the moment u is a place-holder for the confidence interval level we would like to report. Stated another way, if we assume β to be unit normally distributed and we look one standard deviation away from the solution, we have:
−1
1 HO S 2
=
1 σn σn , u
Which has the form of a covariance matrix. Plugging our expression for the Hessian into the above, we see: L,M ∂ 2 fl,m 1 L,M ∂ fl,m ∂ fl,m 2 + 2r ∑ ∂ βn ∂ βn ∑ l,m ∂ βn2 2 l,m=1 l,m=1
−1 S
=
1 σn σn u
(4.9)
Dropping the second term in the Hessian, containing the second derivative, we arrive at the form known as the outer product gradient estimate for the covariance:
−1
L,M
∂ fl,m ∂ fl,m 1 σn σn = ∑ ∂ βn ∂ βn S u l,m=1 −1 1 cov β = J J S u 111
Where J = (L × M) × N, or the Jacobian of the model over a linearized index z of length L × M. For small residuals this approximation is justifiable, and in practice, near a local minimum the correction that the second derivative gives to the Hessian can be less than 10% or less even for moderate signal to noise of 4:1. This observation can be verified by computing a finite difference estimate of the Hessian. Equipped with an estimate of the Hessian, whether it be numerically calculated or via analytic approximation, if we want to know the standard error, we should scale the assumed unit normally distributed variance by the explained variance: σˆ = RSS/ (# degrees of freedom).
cov β =
−1 RSS J J S L×M−N −M×N
From this expression, the standard variance of β may be extracted as the εstd = ! diag cov β . If the off-diagonal elements of the covariance matrix are small, then the standard variance gives a reasonable idea of how “valid” the extracted parameters β are, i.e. we can put an error bar on the numbers extracted in the fit. A related and even easier to understand metric is the relative standard error, which is the ratio of the standard ! error to the estimated parameter: εrel = βi / diag cov β .
4.6
Sum of Lorenzians Model
Throughout this chapter we have alluded to a sum of Lorenzians model, but have kept the exposition in terms of a general model. For typical minimization packages, the objective function must be supplied in terms of real numbers only. Thus, we break the model down as follows:
112
= fl,m α , β
N
∑ exp [−βn,1tl ] (αm,n,1 cos [βn,2tl ] + αm,n,2 sin [βn,2tl ])
n=1
∂ fl.m = tl exp [−βn,1t] (αm,n,2 cos [βn,2tl ] − αm,n,1 sin [βn,2tl ]) ∂ βn,2 ∂ 2 fl.m = −tl2 exp [−βn,1t] (αm,n,1 cos [βn,2tl ] + αm,n,2 sin [βn,2tl ]) 2 ∂ βn,2 ∂ 2 fl.m = −tl2 exp [−βn,1tl ] (−αm,n,1 sin [βn,2tl ] + αm,n,2 cos [βn,2tl ]) ∂ βn,1 ∂ βn,2 ∂ fl.m = −tl exp [−βn,1tl ] (αm,n,1 cos [βn,2tl ] + αm,n,2 sin [βn,2tl ]) ∂ βn,1 ∂ 2 fl.m = tl2 exp [−βn,1tl ] (αm,n,1 cos [βn,2tl ] + αm,n,2 sin [βn,2tl ]) 2 ∂ βn,1 Here I’ve explicitly encoded the fit parameters α and β with sub-indices to clarify their dependency on the time index l, the frequency-frequency index m, and their real or imaginary part (denote 1 or 2 respectively). The n index on β and α refer to the index of the oscillator, where each oscillator has an intrinsic frequency and linewidth. In total, each oscillator requires 2 non-linear parameters and 2 linear parameters.
4.7
Optimal Model Order
Even if a given fit returns what seem to be reasonable parameter estimates, it is not necessarily true that the model employed is the best. A particularly insidious problem that arises when fitting a sum of Lorenzian models is that if too few oscillators are specified in the model, several true frequencies can be averaged into a single fit frequency, rather than the fit only selecting a subset of the true frequencies. Such an error could result in a dramatic change in the interpretation of the fitted data. To insure against these potential mis-interpretations, one should run many possible models and select which is the best. Information criteria can be used to perform this selection task. Although there are many different information criteria in the literature, the corrected Akaike Information Criteria (AICc) offers a simple and theoretically well grounded means for model selection[135]. 113
In information theory, one is seeking to minimize the “information” that is lost when a model is used to parametrize a data set. AIC is a specific implementation of this idea, which forms an estimate of the Kullback-Leiber information “distance” from an unknown true model. An exposition of the metric and its underpinnings is beyond the scope of this thesis, but the main assumptions and limitations of the particular implementation of AICc used in the analysis are discussed here. First, the log-likelihood and the maximum likelihood estimate are given by the least-squares estimate. Second, and more restrictive, the residual of a given fit is assumed to be a normally distributed random variable with constant variance. While bootstrap estimates could give much more accurate estimates of the residual distribution and variance, for the data size used in 2DES, such estimates are computational expensive. With these assumptions, the AIC takes on a simple form: 2K (K + 1) RSS + 2K + AICc = n log n n−K −1 K = numel β + numel (α ) + 1
RSS = r r Where RSS is the residual sum of squares, K describes the number of parameters in the model (2 real numbers per complex oscillator, and 2 per complex amplitude), and n is effective the number of data points (2 real numbers per record if fitting complex data). With respect to the adaptive filtering method, presented in chapter 5, a signal of length N can be perfectly reconstructed as soon as N poles are selected, since this gives a full rank basis for the signal. However, it should be obvious that as soon as one uses N/2 poles, the number of parameters will be equal to the number of data points! The third term in AICc has a singularity at
N 2
poles as well. Thus fitting model orders beyond this point is useless
and so whether one is using variable projection or adaptive filtering the maximum model order one should select is
N 2
− 1. In terms of implementation, one simply runs the fitting
software many times, once for each model order. This is, at least, parallelizable so its a 114
good idea to run the fitting on a multi-core processor with code that can take advantage of this.
[135] Kenneth P Burnham and David R Anderson. Multimodel inference understanding aic and bic in model selection. Sociological methods & research, 33(2):261–304, 2004. [136] Jianhan Chen, AJ Shaka, and Vladimir A Mandelshtam. Rrt: The regularized resolvent transform for high-resolution spectral estimation. Journal of Magnetic Resonance, 147(1):129–137, 2000. [137] G Golub. Numerical methods for solving linear least squares problems. Numerische Mathematik, 7(3):206–216, 1965. [138] Gene H Golub and Victor Pereyra. The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. SIAM Journal on numerical analysis, 10(2):413–432, 1973. [139] Gene H Golub and Charles F Van Loan. Matrix computations, volume 3. JHU Press, 2012. [140] Steven G Johnson. The nlopt nonlinear-optimization package. http://abinitio.mit.edu/nlopt. [141] Linda Kaufman. A variable projection method for solving separable nonlinear least squares problems. BIT Numerical Mathematics, 15(1):49–57, 1975. [142] Sergei Kucherenko and Yury Sytsko. Application of deterministic low-discrepancy sequences in global optimization. Computational Optimization and Applications, 30(3):297–318, 2005. [143] Jorge Nocedal. Updating quasi-newton matrices with limited storage. Mathematics of computation, 35(151):773–782, 1980. [144] Dianne P OLeary and Bert W Rust. Variable projection for nonlinear least squares problems. Computational Optimization and Applications, pages 1–15. [145] Paul O’Leary and Matthew Harker. An algebraic framework for discrete basis functions in computer vision. In Computer Vision, Graphics & Image Processing, 2008. ICVGIP’08. Sixth Indian Conference on, pages 150–157. IEEE, 2008. [146] Michael JD Powell. The bobyqa algorithm for bound constrained optimization without derivatives. Cambridge NA Report NA2009/06, University of Cambridge, Cambridge, 2009. [147] Tao Qian, Wolfgang Sprößig, and Jinxun Wang. Adaptive fourier decomposition of functions in quaternionic hardy spaces. Mathematical Methods in the Applied Sciences, 35(1):43–64, 2012. 115
CHAPTER V
EFFICIENT SIGNAL EXTRACTION FOR 2DES
5.1
Introduction
The need to scan over the coherence time delay dramatically increases the acquisition time of the 2DES experiment relative to more traditional transient absorption measurements, which only scan the waiting time delay. This simple problem, among others, has motivated the development of single shot 2DES setups [152, 158]. Single shot setups, which encode the coherence delay spatially, unfortunately require high speed large area CCD or CMOS detectors (which are expensive) and significantly more pulse energy (since the energy is spatially spread out) than traditional scanning setups. Furthermore the spatial mode quality of the laser source for single shot setups must meet more stringent requirements to avoid aberrations from spatial temporal coupling. So when single shot setups are not feasible to implement (for a variety of reasons) we should instead seek ways to take fewer scan measurements in order to recover the same information in less time. The simplest, demonstrated[73], way to reduce the number of data points required to reconstruct a full 2D spectrum is to use what is called phase locking. Phase locking is a pulse shaping method which allows the measurement to be conducted in a rotating frame. To accomplish phase locking, a phase −ωlockt1 is added to the pulse as it is scanned through t1 , and the result is that the detected signal will now oscillate at a reduced frequency ω − ωlock in t1 . When the lock frequency is brought into resonance with an optical transition then that optical transition will not oscillate at all as a function of t1 . For frequencies away
116
from resonance the oscillation frequency is not zero, but still low. The main experimental implication of phase locking is that the sampling rate necessary to satisfy the Nyquist sampling criterion is significantly reduced. For 120 nm bandwidths in the visible, this means we can safely step through the coherence time in 10 femtosecond steps. For typical electronic linewidths, only 300 femtoseconds need be measured before the signal decays beneath our signal to noise. Thus we require about 30 measurements along the excitation axis to extract a signal. As the pulse shaper we use is capable of stepping the time delay and phase with each laser shot at a 1 KHz repetition rate, this means that an entire 2D spectrum can be sweeped out in as little as 30 ms. But if a further reduction of the acquisition time is desired, we need to look for ways to either break the Nyquist criterion altogether or find ways to reliably extrapolate the data based on a truncated measurement. In this chapter we will explore the latter idea. If the signal we are trying to observe can be sparsely represented in some basis, then that sparsity can be exploited to reduce the number of coefficients we need to estimate – and thus the number of measurements we need to make in order to estimate them. Discovering the basis in which a given signal is sparse is a general problem found in many fields. The simplest way to state the sparse solution problem for a given basis Φ is that we would like to minimize the quadratic error |y − Φx|22 subject to the constraint that the l0 norm of x,|x|0 , is also minimized. Which is to say: find the coefficient vector with the least number of non-zero elements. Unfortunately, this direct problem statement is NP hard to solve[160]. In practice, sparse solutions can be found by solving alternative, more well behaved problems like L1 -minimization, which is called basis pursuit, or greedy algorithms like matching pursuit. In this chapter we pursue an algorithm in the latter class, but with modifications that can make it similar to a basis pursuit approach. A number of papers have recently developed accelerated acquisition methods for multi-dimensional spectroscopies, ranging from NMR[154], to 2DIR[149], and even 2DES[164]. All of these methods seek to exploit a sparsity in the Fourier basis, which is applicable mainly to systems
117
with narrow line-shapes relative to the excitation bandwidth, while the NMR papers are using linear predictive methods. For the broad line-shapes commonly found in electronic spectroscopy of the condensed phase, the Fourier basis obviously does not admit a sparse representation. We demonstrate that these broad 2DES signals can, however, be efficiently represented by an orthonormalized basis of rational polynomials in the frequency domain. To find this representation, we reformulate an existing method called Adaptive Fourier Decomposition[147] as a filtering method for physical signals and develop its discrete time-domain analog. The result is a simple, easily implemented and stable algorithm with few necessary tunable parameters that yields considerable compression relative to the discrete Fourier basis.
5.1.1
Adaptive Filtering Algorithm
2DES signals of condensed phase samples in the frequency domain are typically broad, with smooth overlapping irregularly shaped features. In other words, the signals are inefficiently represented by the delta-function (in frequency) Fourier basis. The ideal basis would be smooth and causal, be capable of simultaneously describing the real and imaginary parts of the nonlinear signal, and also have finite support that smoothly approaches zero as it approaches infinity in either direction. Complex rational polynomials (given by ratios of complex polynomials), satisfy all these requirements and have been extensively developed in both the System Identification and function approximation communities for close to a hundred years[167]. By partial fraction decomposition, any rational polynomial can be represented as a sum of complex Lorenzian functions in the frequency domain or sum of complex exponentials in the time domain. Because Lorenzian functions encode both a bandwidth and central frequency, it should be clear that the problem of broad line-shapes can be handled by such a basis. What is less clear at the outset is how sparsely experimental data will be represented, since Lorenzian lineshapes have significantly broader tails than what is observed. As will be shown in section 5.2,
118
the bulk of the amplitude for simple inhomogeneously broadened lineshapes can be well represented with just a few basis functions. As most 2DES setups involve scanning a time delay, we focus on formulating the problem of extracting the frequency domain from a time series directly. Of course, in 2DES one collects a collection of time series for each 2D spectrum, but here the focus is on a single time series. In section 5.2, the adaptation of this method to two dimensions is made. The problem of fitting time series to sums of complex exponentials is in general non-linear and therefore potentially troublesome. Linear predictive algorithms, which date back to the late eighteenth century beginning with Prony’s method, can reduce the non-linear problem to a linear problem[165]. Modern examples of such “linearizing” algorithms are Linear Predictive Singular Value Decomposition[150], Multiple Signal Classification[166], Filter Diagonalization Method[157], and Regularized Resolvent Transform[136]. Despite the considerable success of linear methods, we consider an adaptive filtering algorithm here primarily because it is straight-forwardly understood and implemented and exhibits excellent numerical stability in the presence of noise without recourse to regularization. We do not make any general claim about the superiority of an adaptive filtering approach over linear methods, however. To begin we introduce a parametrized orthogonal rational polynomial basis Φ known as the Takenaka-Malmquist (TM) basis in System Identification literature[153]. The nth column of Φ is given by: Φn (z) =
1 − |an |2 n−1 1 − a¯ j z z − an
∏
j=1
z−aj
, k∈N
This basis accepts as parameters a collection of complex values a j , called poles inside the unit disk (D) excluding the boundary of the disk (T). The TM basis is orthogonal on the Lebesgue measure, which can be expressed in the complex plane as a contour integral, though it is more commonly known on the unit circle. 119
˛ f , g =
T
f (z) g (1/¯z)
1 ∀z ∈ T = f , g = 2π
ˆπ
dz z f eit g (eit )dt
−π
Because all the poles of the TM basis lie inside the unit circle, the TM basis is orthogonal from [0, ∞) in the discrete time domain, which is understood by examining the radius of convergence of the z-transform for a complex function. Furthermore, the numerator polynomial degree is strictly less than that of the denominator, making the TM basis appropriate for strictly proper transfer functions, which implies that it is capable of describing physical signals. Note that a “pole”, which mathematically just means a singularity of a certain kind, in this context also implies a certain lifetime and frequency of the Lorenzian oscillator. Orthogonality is important for a number of theoretical reasons[153], particularly in estimating convergence and variance properties, but also a key practical one: numerical conditioning. While the orthogonality in theory holds only for infinite record lengths, in double precision the vectors become effectively orthogonal within a few lifetimes of the slowest pole (the lifetime of the least damped oscillator). Even for record lengths on the order of one lifetime of the slowest pole, the basis functions are nearly orthogonal, such that the condition number of Φ Φ is near unity. The basis of individual complex exponentials described by the same poles and record length, on the other hand, can very quickly become ill-conditioned. Thus using this basis acts, at the least, as a preconditioner for the problem of fitting complex exponentials. Recall that our objective is to find the basis in which the signal may be sparsely represented. This now amounts to finding the optimal poles for the TM system such that as few as possible basis functions are needed to achieve a requested accuracy. A number of approaches may be taken, the most straight-forward of which is to use simplex 2 minimization of the variable projection[151] functional ε = y − Φ (a) Φ (a)+ y 2 . Where
120
Φ (a)+ denotes the pseudo-inverse of the basis give the pole vector a and y is the data. Unfortunately, it was found that the best fit done in this manner was often less than satisfactory due to the ruggedness of this objective function. The objective function can be made to be more well behaved by adopting an iterative algorithm called Adaptive Fourier Decomposition (AFD)[162], which is outlined below in brief. Details on the theoretical development and optimality of the algorithm is discussed in[162] and references therein. AFD is an iteration over three steps. At the nth step, 1. Select the pole a j that maximizes the squared absolute value of the inner product of √1−|an |2 the signal Gn (z) with a parametric dictionary element e z, a j = z−an . This is called “Maximal Selection” by the original authors of the AFD algorithm. # " 2. Computing the remainder: Rn = Gn (z) − Gn (z) , e z, a j e z, a j 3. Orthogonalize the remainder with respect to the pole a j , to create the signal for the z−a
next iteration: Gn+1 (z) = Rn 1−a¯ jjz . This process can be repeated N times until the remainder is as small as desired. The idea here can be connected to Non-Linear Least Squares by saying that steps 2 and 3 constitute a different “functional” to minimize. Unlike the variable projection functional, we are approximating the linear coefficients with dot products, rather than general inverses. Furthermore, the non-linear portion of AFD is different in the sense that one does not fit a matrix of functions, but rather a single function and the residual is then filtered at each function fitting step. It’s not a priori clear which method is superior, though we find in practice that the AFD “functional” seems to give better extrapolative performance than variable projection offers. In any event, The final signal estimate Gˆ can be synthesized by:
n−1 1 − a¯
Tn (z) = Gn (z) , e (z, an ) e (z, an ) ∏
j=1
Gˆ (z) =
N
∑ Tn
n=1
121
jz
z−aj
In the frequency domain a closed form for the inner product can be found by complex contour integration and using the property that f , g = g, f . Gn (z) , e (z, an ) = 1 − |an |2 (1/a¯n ) Gn (1/a¯n ). While the frequency domain algorithm here may be used by computing the z-transform of the time series data, there is one problem with this approach (see equation 5.1):
F (z) =
∞
∑ f (k) z−k
(5.1)
k=0
Z-transforming the input data requires that we measure the time series until the response has converged to zero ( f (k) → 0) so that the z-transform will converge to its true value. Doing this would defeat the purpose of estimating the signal with fewer data points than the Fourier transform requires. In fact, the discrete Fourier transform IS the z-transform, for z ∈ T. This problem can be avoided by transforming the AFD iteration to operate directly in time domain, at which point it begins to more closely resemble filtering approaches used for Laguerre modeling[168]. To do this we will create a recursive solution to the time domain:
Gn (z) = Rn−1 (z)
z − an−1 1 − a¯n−1 z
Gn (z) − a¯n−1 zGn (z) = zRn−1 (z) − an−1 Rn−1 (z) gn (k) − a¯n−1 gn (k + 1) = rn−1 (k + 1) − an−1 rn−1 (k)
122
gn (k) =
1
gn (k − 1) a¯n−1 1 − [rn−1 (k) − an−1 rn−1 (k − 1)] , k ∈ [1, N] a¯n−1 i+1 k−1 1 [−rn−1 (k − i)] , k ∈ [1, N] gn (k) = ∑ i=0 a¯n−1 i+1 k−1 1 an−1 rn−1 (k − i − 1) + ∑ i=0 a¯n−1 g1 (k) = y (k)
(5.2) (5.3) (5.4) (5.5) (5.6)
In the closed form expression (5.4) we have used the fact that gn (0) = 0, because the signal is assumed to have a strictly proper transfer function. Even if the time series does not begin with zero value at the first index, a zero can be concatenated at the beginning without loss of generality to satisfy this condition. g1 (k) is given by the input time series. In practical computation, the recursive formula (5.2) is much faster to compute, but degrades in accuracy more rapidly than the closed form version. In calculations presented here we employ the recursive form. The inner product between the signal and the dictionary can be approximated by a discrete dot product gn (k) , e (k, an ) ≈ [gn (k) , e (k, an )], where e (k, an ) is given by (in closed and recursive form): k−1 1 − |an |2 ∑ ain u (k − i) e (k, an ) =
(5.7)
i=0
e (k, an ) = an e (k − 1, an ) + 1 − |an |2 u (k)
(5.8)
Here u (k) is the time domain of the pump pulse. For an impulsive pump, the dictionary element simplifies to 1 − |an |2 akn S (k − 1), where S (k) is the unit step function. To summarize the algorithm in the time domain, we again iterate over three steps: 1. Find an such that the squared absolute value of the discrete dot product between the dictionary (5.7) and gn (k) is maximized. Call this the Discrete Maximal Selection 123
step. 2. Compute the nth remainder rn (k) = gn (k) − [gn (k) , e (k, an )] e (k, an ) 3. Factor the remainder to find signal for the next iteration using equation (5.4). In either the frequency domain or time domain, it is necessary to find which value of an maximizes the overlap of the dictionary element with the signal. This can be accomplished a couple of ways, the simplest of which is a brute-force grid search. We found that a multi-start subplex algorithm implemented in NLopt[163, 155] finds the global maximum of each iteration much more quickly and with higher accuracy than a grid search. One should observe, however, that the algorithm proposed here is finding a greedy solution and therefore may not be the globally optimal solution if all possible pole orderings are considered. One simply might not bother ensuring one has found a global optimum and to live with the representation one gets from the greedy solution. It is significantly faster to compute the greedy solution and it can be useful for exploratory work. A safer, albeit much slower, approach is to find the global solution by selecting the entire vector of poles all at once and then minimize residual of the AFD “functional”. A subplex optimizer again performs well here, by now minimizing the norm of the nth remainder. We note again that this is different than minimizing the variable projection functional. Global optimization is done by optimizing a population of starting points, which is conveniently done through the Multi-Level Single Linkage (MLSL) algorithm in NLopt[156]. A third option is to employ the recently published cyclic AFD variant[161], described below, which acts to refine the output of the greedy algorithm iteratively. (0) (0) (0) 1. Select a vector of poles P(0) = a1 , a2 , · · · , an using greedy adaptive filtering. (0) Remove the first pole from the set to create P = P(0) \ a1 . 2. Using the fixed set of poles P run the AFD iteration to obtain the (n − 1)th (1)
remainder rn−1 starting with the original function to be approximated g (k) by 124
looping over steps 2-3 of the usual AFD algorithm. (I.E., skip Discrete Maximal Selection and use the poles of P ). (1)
3. Select a new pole using Discrete Maximal Selection, operating on rn−1 to obtain a (1)
new pole a1 . 4. Construct the new pole vector
P(1)
(1) = P ∪ a1
(1) (1) (1) ≡ a1 , a2 , · · · , an .
(1) (1) (1) (1) 5. Remove the second pole from P(1) to obtain P = P(1) \ a2 = a1 , a3 , · · · , an . Go to step 2 and repeat through step 5, replacing the kth pole from the set at the kth iteration. Once the entire pole vector is traversed a “cycle” is complete. One can then go back to the beginning of the pole vector and repeat the entire process for N cycles. A good stopping criterion is to specify the relative value of the pole vector norm from cycle to cycle. i.e. Prel = norm P(n) − P(n−1) /norm P(n−1) . Of these three options, we find that the greedy algorithm typically gives an excellent starting point and that both cyclic AFD and the MLSL algorithm eventually refine that starting point to the same answer. The cyclic algorithm converges more slowly than MLSL, but its stopping criteria is well defined, which makes it less “fiddly” to work with than MLSL. Thus for a user who is patient to wait for a good answer, but not patient enough to try optimization parameters, I recommend using cyclic AFD with the relative pole vector norm stopping criteria mentioned in step 5 of the cyclic AFD algorithm.
5.1.2
Extending Adaptive Filtering to Two and Three Dimensions
In order to estimate the entire two dimensional spectrum, or even a sequence of 2D spectra, we use singular value decomposition (SVD) to break down the estimation problem into a sequence of 1D problems, which are solved by the algorithms developed in section 5.1.1. If the 2D data is given by a matrix with the first dimension describing the detection axis and the second axis describing the t1 axis, then the left singular matrix returned by SVD gives an orthonormal basis for the detection axis and the right singular matrix gives 125
an orthonormal basis for the time axis. The bases that SVD returns are optimal in the sense that they capture the maximal variation possible and are typically ordered from most important to least important. It is well known that truncating the SVD basis at some number of columns M < N results in a low-rank approximation of the N column input matrix[169]. Unfortunately, the basis functions returned by SVD are non-parametric, so they cannot be used to extrapolate the time axis further. The strategy is then to independently filter each column of the right singular matrix in order to extrapolate it to longer time. As a result, this strategy is parallelizable and, by filtering only the first M columns of the right singular matrix, one can reduce the problem size considerably. As a consequence, even sequences of hundreds of 2D spectra can be quickly filtered and extrapolated in less than a minute on modern desktop computers. In brief, when given a data matrix D, with the first dimension containing the detection axis and the second dimension containing the t1 axis, we have:
D = W3 ΣT1
(5.9)
T1 = [T1,1 , T1,2 , · · · , T1,N ] Tˆ1, j = adaptive filter T1, j Dˆ = W3 ΣTˆ1 Equation (5.9) is the SVD decomposition of D and thus W3 and T1 are the left and right singular matrices respectively. In order to fit a sequence of 2D spectra, described by a 3-dimensional array of measurements, one can unfold the 3D array into a 2D array, where the first dimension becomes a concatenation of the detection axes for each population time point, while the second dimension remains the t1 axis. SVD is performed on the unfolded array and the filtering procedure proceeds as in the 2D dimensional case. Note, however, that this entire strategy results not in a new basis for the whole data, but in a set of bases for each column of the right singular vector. As a result, it may be that a more parsimonious 126
representation could be constructed if one were to generate a multi-dimensional filtering method. Such methods are possible, but not explored here. In particular for AFD, it is required that one describe the data in a higher dimensional algebra like quaternions[162].
5.2
Applying Adaptive Filtering to Non-Lorenzian Lineshapes
A simple model for 2D lineshapes is given in Chapter 6 [159], which allows one to incorporate homogeneous linewidths and inhomogeneous frequency distributions as phenomenological parameters. In this model, the inhomogeneous broadening is created by convolving the frequency response of the purely homogeneously broadened spectrum with some distribution, which equates to a product of the homogeneously broadened spectrum with that distribution in the time domain. This sort of model is valid when the bath dynamics are much slower than the system dynamics. In general, static inhomogeneity of, for instance, site energies or couplings, will result in more complicated lineshapes as the time scale separation of bath and system dynamics is not generally valid. For the sake of demonstrating how this algorithm treats lineshapes which are not Lorenzian, however, it will suffice. For Gaussian distributed inhomogeneous broadening we have:
Ieg (t) = exp [iωegt − Γegt] S (t1 ,t2 ,t3 ) ∝ Ieg (t1 ) Ieg (t3 ) (1 + exp [−Γeet2 ]) Δ2eg Sr (t1 ,t2 ,t3 ) ∝ S (t1 ,t2 ,t3 ) exp − (t3 − t1 )2 2 Δ2eg Snr (t1 ,t2 ,t3 ) ∝ S (t1 ,t2 ,t3 ) exp − (t3 + t1 )2 2 Selecting a larger inhomogeneous broadening parameter relative to the homogeneous linewidth, we generate a example rephasing spectrum shown in Figure 5.1. In this case it should be clear that fitting this lineshape with Lorenzians is not optimal, since the
127
Figure 5.1: Simulated rephasing line shape with homogeneous linewidth of 1/ (2ΔT ) and inhomogeneous Gaussian distributed linewidth of 1/ (5ΔT ). a
c
b
Figure 5.2: Reconstruction error relative to the maximum intensity of the full spectrum of the three methods mentioned in section 5.1.1. All methods run using between 1-7 poles using the first 20 data points (of 100). a. Greedy algorithm, b. Cyclic AFD with 1E-6 relative pole norm stopping criteria, c. MLSL algorithm. Equation for the error is: ε = (yˆ − y) /max (y).
tails of the spectrum in the frequency domain will damp more quickly than a purely homogeneously broadened lineshape. Nevertheless, it is still possible to reconstruct the spectrum using few data points (and few basis functions). In essence the adaptive filtering is approximating the continuous distribution of homogeneous linewidths present in the spectrum by a weighted discrete sum. In figure 5.2, 20 t1 time points are used to reconstruct the spectrum using between 1-7 basis functions, depending on the optimal order determined via an information criteria (see section 4.7). It can be seen that the greedy algorithm has higher reconstruction error, while the refinements of the greedy solution, either via cyclic filtering or MLSL, are nearly the same.
128
a
b
Figure 5.3: Real rephasing spectrum of chlorophyll a in a 50/50 mixture of isopropanol and glycerol at 77K taken at 170 fs population time delay. a. FFT of a truncated data set containing 20 data points sampled every 5 fs. b. Cyclic adaptive filter reconstruction of the same truncated data set in a. Note that the excitation axis modulation is eliminated in the adaptive filtered version without broadening the peakshape, as a fixed Fourier filter would.
5.3
Applying Adaptive Filtering to Data
To demonstrate this method on real data, some broad-band 2DES data of chlorophyll a in isopropanol and glycerol at 77K is examined. We show that for truncated data sets, adaptive filtering allows one to extract 2D spectra that eliminate Gibbs phenomenon that results from abruptly truncating and zero-padding an FFT. The ability to collect truncated data (which decreases acquisition time) is useful, since it means more spectra can be collected in the same period of time, which allows one to average out the differing laser noise environments present between spectra. Furthermore, the fact that we are fitting a function which does not completely span the discrete space of the data means that noise which does not fit within the model is removed to some degree. In other words, poles which represent pure noise terms are rejected by selecting a model order that is less than the number of data points.
5.4
Conclusions
We have demonstrated that complex Lorenzians form a reasonably sparse basis to expand 2D electronic signals. The difficulty, however, is to find the appropriate parameters for the Lorenzian basis so that the data can be well represented. We have proposed the use 129
of an adaptive filtering scheme, based on existing frequency domain system identification techniques, and have adapted it to the time domain, which is relevant for typical 2DES data collection setups. Once the signal is expanded in some sparse basis, we can reduce the acquisition time required to represent the full frequency domain data. As mentioned in the introduction of this chapter, significant reductions in acquisition times can be accomplished through the use of phase locking. We find that such measurements can be reduced further by nearly about a factor of 2 by using adaptive filtering to extrapolate the data. The reason for pursuing a reduction in acquisition time so doggedly is because improved acquisition speed results in higher signal to noise along the population time axis by integrating the measurement over many scan instances. Integrating many waiting time scans is a popular strategy in transient absorption, but one which was often denied to us because the acquisition time of a single scan was many hours. The amount of acquisition time reduction one can achieve depends on how well the algorithm can find such a sparse representation and more fundamentally how well the basis approximates the true form of the data. The need for a flexible and adaptive basis was the original reason for pursuing AFD. Unfortunately, even though the idea of using Lorenzians to fit the data appears to be good, the adaptive filtering algorithm, like linear predictive algorithms, requires that the Nyquist sampling criteria still be met and that the sample spacing be linear in order to find a sparse representation. True “compressive sensing” techniques[148] avoid both of these issues, but there new and more esoteric constraints arise and the minimization process becomes much more expensive. In order to reduce the acquisition time further than is possible with phase locking alone, it was necessary to demonstrate that 2DES signals admit a sparse representation in some basis. The adaptive filtering scheme I have presented here suggests is that a subset of the space of all possible Takenaka-Malmquist basis functions admits a sparse representation of the data, which is something not a priori obvious in the presence of inhomogeneous broadening. Exploitation of this observed sparsity in more modern compressive sensing approaches remains to be
130
explored in this context.
[148] Emmanuel J Candes, Yonina C Eldar, Deanna Needell, and Paige Randall. Compressed sensing with coherent and redundant dictionaries. Applied and Computational Harmonic Analysis, 31(1):59–73, 2011. [149] Josef A Dunbar, Derek G Osborne, Jessica M Anna, and Kevin J Kubarych. Accelerated 2d-ir using compressed sensing. The Journal of Physical Chemistry Letters, 4(15):2489–2492, 2013. [150] Jorge L Galazzo and James E Bailey. Application of linear prediction singular value decomposition for processingin vivo nmr data with low signal-to-noise ratio. Biotechnology techniques, 3(1):13–18, 1989. [151] Gene H Golub and Victor Pereyra. The differentiation of pseudo-inverses and nonlinear least squares problems whose variables separate. SIAM Journal on numerical analysis, 10(2):413–432, 1973. [152] Elad Harel, Andrew F Fidler, and Gregory S Engel. Single-shot gradient-assisted photon echo electronic spectroscopy. The Journal of Physical Chemistry A, 115(16):3787–3796, 2010. [153] Peter SC Heuberger, Paul MJ Van den Hof, and Bo Wahlberg. Modelling and identification with rational orthogonal basis functions. Springer, 2005. [154] Daniel J Holland, Mark J Bostock, Lynn F Gladden, and Daniel Nietlispach. Fast multidimensional nmr spectroscopy using compressed sensing. Angewandte Chemie, 123(29):6678–6681, 2011. [155] Steven G Johnson. The nlopt nonlinear-optimization package. http://abinitio.mit.edu/nlopt. [156] Sergei Kucherenko and Yury Sytsko. Application of deterministic low-discrepancy sequences in global optimization. Computational Optimization and Applications, 30(3):297–318, 2005. [157] Vladimir A Mandelshtam. Fdm: the filter diagonalization method for data processing in nmr experiments. Progress in Nuclear Magnetic Resonance Spectroscopy, 38(2):159–196, 2001. [158] Ian P Mercer. Angle-resolved coherent optical wave mixing. Physical Review A, 82(4):043406, 2010. [159] Shaul Mukamel. Principles of nonlinear optical spectroscopy, volume 29. Oxford University Press New York, 1995. [160] Balas Kausik Natarajan. Sparse approximate solutions to linear systems. SIAM journal on computing, 24(2):227–234, 1995. 131
[161] Tao Qian. Cyclic afd algorithm for the best rational approximation. Mathematical methods in the applied sciences, 2013. [162] Tao Qian, Liming Zhang, and Zhixiong Li. Algorithm of adaptive fourier decomposition. Signal Processing, IEEE Transactions on, 59(12):5899–5906, 2011. [163] Thomas Harvey Rowan. Functional stability analysis of numerical algorithms. 1990. [164] Jacob N Sanders, Semion K Saikin, Sarah Mostame, Xavier Andrade, Julia R Widom, Andrew H Marcus, and Alan Aspuru Guzik. Compressed sensing for multidimensional spectroscopy experiments. The Journal of Physical Chemistry Letters, 3(18):2697–2702, 2012. [165] Louis L Scharf. Statistical signal processing, volume 98. Addison-Wesley Reading, MA, 1991. [166] Ralph O Schmidt. Multiple emitter location and signal parameter estimation. Antennas and Propagation, IEEE Transactions on, 34(3):276–280, 1986. [167] Satoru Takenaka. On the orthogonal functions and a new formula of interpolation. Japanese Journal of Mathematics, 2:129–145, 1925. [168] Bo Wahlberg. System identification using laguerre models. Automatic Control, IEEE Transactions on, 36(5):551–562, 1991. [169] Jieping Ye. Generalized low rank approximations of matrices. Machine Learning, 61(1-3):167–191, 2005.
132
CHAPTER VI
CONCLUSIONS AND FUTURE WORK
6.1
Summary and Innovations of this Work
This thesis presents my work to improve both the detected signal bandwidth for 2DES measurements of d1d2[172] and to improve the signal to noise of the experiment in general[173]. The former study showed the possibility of resolving excitation dependent kinetics in weak blue transitions of d1d2, but the signal to noise needed improvement. With the new setup that was subsequently developed it is now possible to return to that study, as was shown in preliminary data on an anion band at 790-820 nm. Combined with work to improve pump source bandwidth and other instrumental changes, we are now able to resolve nearly twice the excitation bandwidth of previous 2DES studies on d1d2 done in this group[174]. The extended bandwidth along with improved pulse compression and signal to noise has made the high quality observations of coherent dynamics of d1d2 possible. Improvements to data fitting procedures compared to those presented in [174] has enabled us resolve multiple kinetic rates and eight coherent oscillations from the data. Additionally, the move away from a pump-probe geometry to a non-collinear box-CARS geometry has greatly facilitated the separation of rephasing and non-rephasing pathway signals, which is important for analysis of coherent effects. While the box-CARS geometry requires phasing, we have clarified methods to implement the phasing (be it through minimization procedures or through the new analytic method presented) and in so doing have achieve better results than published elsewhere[170]. Finally, we have
133
developed a method of reducing the number of data points needed to generate 2DES spectra in a way that is general enough to describe both broad and narrow lineshapes. Of the coherent oscillations in d1d2 that we observe, the 246 cm−1 , 339 cm−1 , 494 cm−1 , and 730 cm−1 modes were identified as having mixed vibrational/excitonic character. The higher frequence modes observed (854 and 971cm−1 ) were identified as purely vibrational. Two of the lower frequency modes (89 and 117 cm−1 ) have less clear origin, though an interesting speculation for the 117 cm−1 mode from simulations of a dimer system suggest that coherent mixing between excitons and vibrational states can give rise to numerous peaks in the waiting frequency response. Simulations further demonstrate that mixing between vibrational and excitonic states leads to an increase in charge transfer rate, which is an indication that the observations of coherence here are relevant for the function of d1d2 and potentially other charge transfer systems.
6.2
Future Work
An excellent ground-work for future experimentation has been established here, and the observation of coherent dynamics in the PSII RC raises a number of questions. An initial focus of this work was to develop the instrumentation to extend the bandwidth to detect many transitions. Some initial work in probing anion bands at 890-820 has been presented, but more work needs to be done here. Given the design principle behind the improved signal to noise of this setup, though, returning to measurements which span from the Qy to pheophytin anion bands in the blue will require the development of an amplified probe over this bandwidth. Some exciting work to create a strong 200 nm bandwidth probe over the 700-900 nm region using a degenerate optical parametric amplifier[176] has already well under way in the lab. Such a source would be able to resolve the weak stimulated emission bands, excited state absorptions, and anion bands of d1d2 in that region. While simulations suggest that coherent effects can enhance charge separation, it remains to be experimentally demonstrated. Given the pulse shaping capabilities this lab 134
has, a coherent control experiment is appropriate to attempt. Such an experiment would involve pumping d1d2 with a train of pulses, rather than just a single pulse or pair of pulses. When the time spacing between pulses in a pulse train is equal to the period of a vibrational coherence, the mode can be “driven” or preferentially excited over other modes. This sort of experiment would be ideal in the sense that it allows one to directly answer whether exciting a particular coherence has an impact on the generation of a photoproduct relative to exciting off-resonance. Investigating how general the coherent dynamics are is another important direction. In other aggregates of chlorophyll, what do the coherence maps look like? How closely do they match the data here? Do protein mutations affect coherences? Though it will be difficult, a temperature dependent study is straight-forward to accomplish and may yield significant insight into the nature of the coherences observed in d1d2 or other systems[175]. And, as mentioned in the introduction, a data-based model of long-time kinetics both for the PSII RC and Core complexes to rival those created from transient absorption data needs to be developed. The main question is whether or not a more complex model of the kinetics is justified. It is hoped that the global data fitting techniques demonstrated here will be of utility in that study, since target models generally involve a large number of non-linear rates components. Finally, some work to understand the physical nature of vibrations and how they impact the electronic state would be helpful. From molecular dynamics simulations and Raman studies the 261 cm−1 and 350 cm−1 modes have been identified and characterized as originating from both in and out of plane deformations delocalized over the entire chlorophyll molecule[177]. Ultimately electronic structure calculations capable of describing the interplay of vibrational motion and excited state electron density could be useful for grounding the concepts of vibrational electronic coupling to a more physical picture.
135
[170] J. M. Anna, E. E. Ostroumov, K. Maghlaoui, J. Barber, and G. D. Scholes. Two-dimensional electronic spectroscopy reveals ultrafast downhill energy transfer in photosystem i trimers of the cyanobacterium thermosynechococcus elongatus. Journal of Physical Chemistry Letters, 3(24):3677–3684, 2012. [171] Minhaeng Cho. Two-dimensional optical spectroscopy. CRC press, 2010. [172] FD Fuller and JP Ogilvie. Continuum probe two-dimensional electronic spectroscopy of the photosystem ii reaction center. In EPJ Web of Conferences, volume 41, page 08018. EDP Sciences, 2013. [173] Franklin D Fuller, Daniel E Wilcox, and Jennifer P Ogilvie. Pulse shaping based two-dimensional electronic spectroscopy in a background free geometry. Optics Express, 22(1):1018–1027, 2014. [174] Jeffrey A Myers, Kristin LM Lewis, Franklin D Fuller, Patrick F Tekavec, Charles F Yocum, and Jennifer P Ogilvie. Two-dimensional electronic spectroscopy of the d1-d2-cyt b559 photosystem ii reaction center complex. The Journal of Physical Chemistry Letters, 1(19):2774–2780, 2010. [175] Vaclav Perlik, Craig N Lincoln, Frantisek Sanda, and Jurgen Hauer. Distinguishing electronic and vibronic coherence in 2d spectra by their temperature dependence. The Journal of Physical Chemistry Letters, 2013. [176] AM Siddiqui, G Cirmi, D Brida, FX Kärtner, and G Cerullo. Generation of< 7 fs pulses at 800 nm from a blue-pumped optical parametric amplifier at degeneracy. Optics letters, 34(22):3592–3594, 2009. [177] Chengli Zhou, James R Diers, and David F Bocian. Q y-excitation resonance raman spectra of chlorophyll a and related complexes. normal mode characteristics of the low-frequency vibrations. The Journal of Physical Chemistry B, 101(46):9635– 9644, 1997.
136
