ULTRAFAST TRANSIENT STATES IN NONEQUILIBRIUM QUANTUM SYSTEMS By Bin Hwang A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics – Doctor of Philosophy Computational Mathematics, Science and Engineering – Dual Major 2017 ABSTRACT ULTRAFAST TRANSIENT STATES IN NONEQUILIBRIUM QUANTUM SYSTEMS By Bin Hwang Photo-induced phase transitions (PIPT) in quantum systems are the epitome of challenging non-equilibrium many-body phenomena, that also have a wide range of potential applications. Recently interest in light-matter coupled states with an energy gap have yielded evidence for Floquet topological states. In this study we demonstrate nonequilibrium Floquet band formation under ultrafast optical excitation using a one-dimensional topological insulator. As an example, the effects are illustrated using a new Zig-Zag Su-SchriefferHeeger model of polyacetylene, which is a paradigmatic Hamiltonian exhibiting nontrivial edge states. Our results indicate short optical pulses feasible in experiments can induce novel topological states, local spectral selection and novel pseudospin textures in polyacetylene. Pump-probe photoemission spectroscopy is able to study these states by measuring Floquet band formation and sizeable energy gaps on femtosecond time scales. We find that optically activated nontrivial variations of sublattice mixing could lead to novel topological phenomenon. The rich variety of states induced by lasers have a wide range of potential applications so that control of these states has become a key objective. We present a computational approach to finding optimal ultrafast laser pulse shapes to induce target states and population inversion in pump-probe PIPT experiments. The Krotov approach for Quantum optimal control theory (QOCT) is combined with a Keldysh Green’s function calculation to describe experimental outcomes such as photoemission, transient single particle density of states and optical responses. Results for a simple model charge density wave (CDW) system are presented, including generation of almost complete population inversion and negative temperature states. Copyright by BIN HWANG 2017 This thesis is dedicated to my family. iv ACKNOWLEDGEMENTS It has been a long and exciting journey of pursuing my Ph.D. I would like to thank many people who have guided me, helped me and supported me. First of all, I must thank my advisor and mentor, Phillip Duxbury. I sincerely appreciate him for his help in all aspects of my research career at Michigan State University. He is not only one of the smartest teachers I know but also the most thoughtful friend as well. Because of him, I have investigated deeply the physical world and also learned wisdom in research and life. His wide knowledge and passion in physics kept me wanting to dig out hidden truths behind nature, and his patience and support were great help when I faced challenges in research and life. Second, I appreciate greatly the help from the rest of the PA professors and members of our research group. I would like to thank Jennie Portman, Xukung Xiang, Faran Zhou, Professor Chong-Yu Ruan, and Professor Carlo Piermarocchi for very useful discussions in all kinds of theoretical and computational physics research topics. I also would like to thank Professor Jianliang Qian and Professor Sean Couch for their support during my defense process. This work was supported by the College of Natural Science at Michigan State University (MSU). Computational resources were provided by the High Performance Computer Center at MSU. Finally, I would like to thank my lovely wife Chi-Chun Hsieh and my family members for their continued support and encouragement through all of this. v TABLE OF CONTENTS LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . 1.1 Time-resolved ARPES . . . . . . . . . . . . . . . . . . 1.1.1 Nonequilibrium Regime . . . . . . . . . . . . . 1.2 CDW Materials . . . . . . . . . . . . . . . . . . . . . . 1.2.1 CDW in 1T − T aS2 . . . . . . . . . . . . . . . 1.2.2 Optimal Laser Pulse Shaping . . . . . . . . . . 1.3 Floquet-Bloch States . . . . . . . . . . . . . . . . . . . 1.3.1 Laser Driven Topological Phases . . . . . . . . . 1.4 Chapters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 3 5 6 8 10 13 13 CHAPTER 2 MODELS AND METHODS . . . . . . . . . . 2.1 Introduction to Nonequlibrium Green’s Functions . . . 2.2 Introduction to Topology in Band Structure . . . . . . 2.2.1 Topology . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Berry Phase and Chern Number . . . . . . . . . 2.2.3 Bulk Edge Correspondence . . . . . . . . . . . . 2.3 Tight-binding Models at Equilibrium . . . . . . . . . . 2.3.1 SSH, Rice Mele, Platero Models . . . . . . . . . 2.3.2 Freerick’s Model . . . . . . . . . . . . . . . . . . 2.3.3 Green’s Function Method . . . . . . . . . . . . 2.3.4 Graphene and Dirac Cone . . . . . . . . . . . . 2.4 Tight Binding Models in a Time Varying Field . . . . . 2.4.1 Peierl’s Substitution . . . . . . . . . . . . . . . 2.5 Introduction to Floquet Theory . . . . . . . . . . . . . 2.5.1 Periodic Hamiltonian . . . . . . . . . . . . . . . 2.5.2 Floquet Formalism . . . . . . . . . . . . . . . . 2.6 Calculation of Time Resolved Photoemission . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 16 19 20 21 23 26 26 34 38 40 43 43 48 48 49 52 CHAPTER 3 THE ZIG-ZAG SSH MODEL: FROM FLOQUET TO TRPES 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Model and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Time-dependent Zig-Zag SSH Model . . . . . . . . . . . . . . 3.2.2 Floquet Hamiltonian for Zig-Zag SSH Model . . . . . . . . . . 3.2.3 Calculations from Non-equilibrium Green’s Functions . . . . . 3.3 Analysis of the SSH Floquet Hamiltonian . . . . . . . . . . . . . . . . 3.3.1 Limit of High Frequency . . . . . . . . . . . . . . . . . . . . . 3.3.2 Simple Model for the Case of One Overlap . . . . . . . . . . . 3.3.3 Computational Results and Comparison to Analytic Results . 3.3.4 Effects of Varying Amplitude and Frequency . . . . . . . . . . 3.3.5 Topological Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 56 56 57 59 60 61 61 62 63 66 69 vi 3.4 3.5 Pump-probe Results . . . . . . . . . . . . . . . . . . . . . 3.4.1 Pseudospin Content . . . . . . . . . . . . . . . . . . 3.4.2 Evolution of States and Their Occupancy . . . . . . 3.4.3 Dynamics of Phase Evolution . . . . . . . . . . . . 3.4.4 Effects of Varying Pump and Probe Pulsewidth and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 4 KROTOV OPTIMAL CONTROL THEORY 4.1 Introduction to Optimal Control . . . . . . . . . . . 4.2 Preliminary Preparation of the Krotov Method . . . 4.3 The Tricks of Krotov’s Method . . . . . . . . . . . 4.3.1 Decomposition of Goal Function . . . . . . . 4.3.2 Iterative Algorithm . . . . . . . . . . . . . . 4.3.3 Monotonic Convergence of Krotov Method . 4.4 Construction of ϕ . . . . . . . . . . . . . . . . . . . 4.4.1 First Order in x . . . . . . . . . . . . . . . . 4.4.2 Second Order in x . . . . . . . . . . . . . . 4.4.3 Algorithm . . . . . . . . . . . . . . . . . . . 4.5 Examples . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Example 1: First Order in x . . . . . . . . . 4.5.2 Example 2: Second Order in x . . . . . . . . 4.5.3 Example 3: A Two Level Quantum System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Amplitude . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 71 74 76 78 80 82 82 83 84 84 86 87 88 88 90 91 91 91 92 94 CHAPTER 5 QUANTUM OPTIMAL CONTROL OF A TRANSIENT CDW STATE 98 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.2 Charge Density Wave Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 5.2.1 Equations for the Nonequilibrium Solution . . . . . . . . . . . . . . . 102 5.3 Quantum Optimal Control Method . . . . . . . . . . . . . . . . . . . . . . . 104 5.3.1 Proof of Optimality . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.4 Application to a CDW Model . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 CHAPTER 6 SUMMARY AND OUTLOOK . . . . . . . . . . . . . . . . . . . . . 116 CHAPTER 7 EXPERIMENTS AND APPLICATIONS . . . . . . . . . . . . . . . 118 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 vii LIST OF FIGURES Figure 1.1: The experimental design of time-resolved ARPES, and corresponding ultrafast momentum-dependent dynamics of the melting of the charge density wave state in T bT e3 monitored by trARPES. The figure is taken from [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Figure 1.2: Schematic picture for the comparison between uniform electron density and the Peierl’s transition of a charge density wave in a 1D system. . . . 6 Figure 1.3: (a) The resistivity versus temperature in different charge-density-wave (CDW) phases of 1T − T aS2 . (b) Schematic of the T a atom layer in the CCDW phase. 13-atom David-star clusters with the three T a sites a, b, and c. The black arrows indicate the displacement of the T a atoms from their original positions. Figure is take from [2]. . . . . . . . . . . . 7 Figure 1.4: Pump-induced gap magnitude oscillations in the CDW material 1T − T aS2 near 30K. Using a 1.5eV pump pulse and 6eV probe pulse to induce instantaneous collapse of the Mott gap and to observe its recovery. The oscillation in the band is interpreted as the CDW amplitude mode related to electron-phonon coupling. Figure is taken from [3]. . . . . . . . 8 Figure 1.5: A closed-loop for an optimal control process in quantum systems. First, input initial conditions to the quantum system, such as an initial random field. Second, a current laser control field design is created with a pulse shaper and then applied to the sample. Third, The system evolves and the results are fed to a learning algorithm to suggest an improved laser field for repeated excursions around the loop until the objective is achieved. Figure is taken from [4]. . . . . . . . . . . . . . . . . . . . . . 10 Figure 1.6: An idealized band structure for a topological insulator. The bulk band gap is traversed by topologically-protected surface states [5]. . . . . . . 11 Figure 1.7: Researchers recently showed experimental results of the Floquet-Bloch states [6] through pump-probe trARPES measurements on the topological insulator Bi2 Se3 . Figure is taken from [6]. . . . . . . . . . . . . . . . 12 Figure 2.1: Kadanoff-Baym-Keldysh contour in the complex time plane. . . . . . . . 18 Figure 2.2: A coffee mug can smoothly deform to a donut since they both have one hole or genus g = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 ˆ for a Figure 2.3: The Berry phase is equal to half the sold angle swept out by d(k) two band theory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Figure 2.4: Edge states as skipping cyclotron orbits. . . . . . . . . . . . . . . . . . . 24 viii Figure 2.5: The chiral edge states connect the valence band near K and K 0 with linear or nonlinear dispersion for linear and nonlinear Hamiltonians respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Figure 2.6: Two distinct topological states of the SSH model. (a) is the case δτ > 0 and (b) is the case δτ < 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Figure 2.7: Two distinct topological paths of the SSH model. (a) is the case δτ > 0 and (b) is the case δτ < 0. . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 2.8: Dispersion relation of the SSH model based for R = τ + δτ and r = τ − δτ . 31 Figure 2.9: Topological path of the SSH model depends on R, r and ka. . . . . . . . 31 Figure 2.10: A dimers chain with two sublattices A and B, and τ and τ 0 are hopping terms between two sublattices A and B. a0 is the cell size and b0 is the distance between A-B atoms. . . . . . . . . . . . . . . . . . . . . . . . . 33 Figure 2.11: CDW model with two different sublattices A and B. . . . . . . . . . . . . 35 Figure 2.12: Dispersion relation εk± of the CDW model when µ = 1/2 and U = 1. . . 38 Figure 2.13: Graphene is an atomic-scale hexagonal lattice made of carbon atoms. . 40 Figure 2.14: Dirac cone–the electrons in graphene behave like massless Dirac particles that appear in the electronic band structure as gapless excitations with a linear dispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Figure 3.1: The Zig-Zag SSH model: (a) the trivial topological state δτ > 0 and (b) the nontrivial topological state δτ < 0. Note that the polarization is linear so the electric field is along one bond but perpendicular to the other. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Figure 3.2: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 5. Dashed lines: the exact solution for the 4×4 Floquet Hamiltonian. Solid lines: Numerical results for the 42×42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Figure 3.3: Scheme of band interactions for g0 , g1 and g2 . n is the Floquet band index. g0 is the intraband process; while g1 and g2 are interband process. 65 ix Figure 3.4: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 3. Dashed lines: the exact solution for the 4×4 Floquet Hamiltonian. Solid lines: Numerical results for the 42×42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Figure 3.5: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 1.5. Dashed lines: the exact solution for the 4×4 Floquet Hamiltonian. Solid lines: Numerical results for the 42×42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Figure 3.6: Dispersion relations of quasi-static Floquet Hamiltonian with drive frequency Ω = 2.8. Left: trivial condition δτ = 0.5 and Right: the nontrivial condition δτ = −0.5. (a)(b) is for amplitude A0 = 2, (c)(d) is for amplitude A0 = 4 and (e)(f) is for amplitude A0 = 6. . . . . . . . . 67 Figure 3.7: Quasienergy spectrum properties calculated from the Floquet matrix with m = 10 for the original bands at k = π vs laser amplitude A0 for δτ = ±0.5. Left: are for the trivial condition δτ = 0.5 and Right: are for the nontrivial condition δτ = −0.5 for different driving frequencies Ω. 68 Figure 3.8: Bandwidth calculated from quasi-static Floquet band structures vs laser amplitude A0 for δτ = ±0.5. Left: are for the trivial condition δτ = 0.5 and Right: are for the nontrivial condition δτ = −0.5 for different driving frequency Ω. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 Figure 3.9: Topological phase diagram for different frequencies and amplitudes of the drive pulses from effective Floquet Hamiltonian. (a)(b) shows the cases for δτ = 0.25 on the left and δτ = −0.25 on the right. (c)(d) shows the cases for δτ = 0.5 on the left and δτ = −0.5 on the right. (e)(f) shows the cases for δτ = 0.75 on the left and δτ = −0.75 on the right. The four red points indicate the frequencies and amplitudes we discuss in Fig. (3.16) of the next subsection. . . . . . . . . . . . . . . . 70 Figure 3.10: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = 0 under the pulse frequency Ω = 1, width σ = 50 and amplitude A0 = 1. Dashed lines show the corresponding quasi-static Floquet band structures at the same frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 x Figure 3.11: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = 0.5 under the pulse frequency Ω = 3, width σ = 50 and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. . . . 75 Figure 3.12: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = −0.5 under the pulse frequency Ω = 3 and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. . . . . . . . . . 75 Figure 3.13: Berry phase mapping using the radius r0 = 1 + (1 + k)/2 with δτ = 0.5 on the left and δτ = −0.5 on the right under the pulse frequency Ω = 3, amplitude A0 = 2 and pulse width σ = 50. (a)(b) is at time t = −200 with integral frequency window ωl = 0.0 and ωu = 3.0, (c)(d) is at time t = 0 with integral frequency window ωl = −1.5 and ωu = 0.0 and (e)(f) is at time t = 0 with integral frequency window ωl = 1.5 and ωu = 3.0. . 77 Figure 3.14: Time evolution of pseudospin phase ψ deriving from pseudospin context P eiψ = Px + iPy are presented for time t = −200, −100, −50, −25, 0 with δτ = 0.5 on the top (a)-(e) and δτ = −0.5 on the bottom (f)-(j) under the pulse frequency Ω = 3, amplitude A0 = 2 and pulse width σ = 50. Dashed lines show the corresponding quasi-static Floquet band structures. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Figure 3.15: Time-resolved PES of the electronic SSH model with δτ = 0 at time t = 0 with a pulse of frequency Ω = 3, amplitude E = 1 and width (a) σ = 1, (b) σ = 2, (c) σ = 5, (d) σ = 10, (e) σ = 20, (f) σ = 40. Dashed lines show the corresponding quasi-static Floquet band structures. . . . . 79 Figure 3.16: Time-resolved PES with δτ = 0.5 on the left and δτ = −0.5 on the right at time t = 0. (a)(b) with pulse of frequency Ω = 1.6, amplitude A0 = 0.5 and width is σ = 50, (c)(d) with pulse of frequency Ω = 1.6, amplitude A0 = 2.0 and width is σ = 50, (e)(f) with pulse of frequency Ω = 3, amplitude A0 = 0.5 and width σ = 50, (g)(h) with pulse of frequency Ω = 3, amplitude A0 = 2 and width σ = 50. Dashed lines show the corresponding quasi-static Floquet band structures. . . . . . . 81 Figure 4.1: Left: fidelity versus iteration times. Right: optimal control sequence with respect to time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Figure 4.2: Left: Cost function versus iteration times. Right: Optimal evolution of x with respect to time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Figure 4.3: Left: Cost function versus iteration times. Right: Optimal control pulse of ε respect to time t. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 Figure 4.4: Optimal probability evolution of ψ(t) respect to time t. . . . . . . . . . . 97 xi Figure 5.1: Time-resolved PES of the electronic CDW system before, during and after application of the optimal pulse (Ac (t), As (t), A(t)) from time t = 0 to t = 20 h/t∗ with the initial state a filled lower band at temperature T = 0. The target is excitation of a single state with the goal energy 1t∗ . In this calculation λ(t) = 2πh/∆ is a constant, with ∆ = 0.02 and h = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 Figure 5.2: Time-resolved PES, with the optimal pulse interacting with the system from time t = 0 to t = 20 h/t∗ . The initial state is a filled valence band at temperature T = 0, and the goal is to excite all of the states in the conduction band. The electron density at a given energy is plotted in false color and is high near the band edges as this model has a square root singularity in the density of states at both the valence band and conduction band edges. In this calculation λ(t) = 2πh/∆ is a constant, with ∆ = 0.02 and h = 10. . . . . . . . . . . . . . . . . . . . . . . . . . 114 Figure 5.3: Properties of the system with the optimal pulse acting on the system from time t = 0 to t = 20 h/t∗ (a) The average occupancy of the conduction and valence bands n± (t), (b) The real space order parameter Ω(t), which measures the difference in occupancy of the two sublattices, (c) The energy hH(t)i of the CDW system. . . . . . . . . . . . . . . . . 115 Figure 7.1: Angle-resolved photoemission spectra (APRES) of Bi2 Se3 . (A) A sketch of the experimental geometry for the p-polarized case. kx is defined to be the in-plane electron momentum parallel to the pump scattering plane. (B to F) ARPES data for several pump-probe time delays t (values indicated in the figure) under strong linearly polarized mid-infrared (MIR) excitation of wavelength l = 10 mm. Figure and caption are taken from [7].118 Figure 7.2: (a) Waveguide array fabricated out of PMMA on top of a Cr- and Aucoated glass substrate. Alternating center-to-center separations, 600 and 1000 nm, implement the bulk SSH model. (b) Plasmonic waveguide array incorporating a topological defect where the long separation is repeated twice. Three different excitation sites, I, II, and III, are highlighted. Figure and caption are taken from [8]. . . . . . . . . . . . . 120 Figure 7.3: Absorption images taken after 760µs(≈ 0.78h/t) of evolution following the initialization and quench, corresponding to phases of π (top) and 0 (bottom), respectively, for ∆/t = 0.38(1). Figure and caption are taken from [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Figure 7.4: (a) Sketch of the passive waveguide array acting like a PT-symmetric structure with a topological interface. (b) End facet of the experimentally realized structure in fused silica glass. Image courtesy of Mark Kremer, FSU Jena. Figure and caption are taken from [10]. . . . . . . . 121 xii CHAPTER 1 INTRODUCTION In physics, using spectroscopy to investigate interactions between light and matter has been a long standing and broad area of study. Since the development of the laser as a source for coherent light, many new methods have been utilized in a variety of fields, such as laser pulse shaping and nonlinear optics. With advanced techniques, ultrafast pulses of durations from picoseconds to a few femtoseconds are available; and more recently pulses in the attosecond range have been generated. Source with frequencies from the terahertz to the X-ray are now available. More intensive optical laser pulses have been employed recently to generate tabletop X-ray and extreme-ultraviolet sources enabling the shortest pulses with a duration in the range of attoseconds [11, 12, 13, 14, 15]. These tremendous improvements have enabled spectroscopy of matter with time resolution at the femtosecond scale. Some of the most interesting properties only happen at such time-scales, particularly the dynamics of quantum states. In the last decade, many researchers have used ultrafast spectroscopy to study correlated and nanoscale quantum materials; including superconductors, charge density wave materials and topological insulators, and novel effects and new phases that only occur far from equilibrium are being discovered. Angle-resolved photoemission spectroscopy (ARPES), is of special interest, as it is able to resolve the occupied electronic structure in momentum space; yielding unique insight into many-body correlations and the dispersion of quasiparticles in bulk and low-dimensional superconductors, strongly correlated materials and topological insulators. In the time domain, applying an ultrafast laser pulse to a material and studying the related time-resolved ARPES spectra (trARPES), one can study the dynamics of the electronic structure and quasiparticle occupations; and discover the hidden structure of quantum states in physics. 1 Figure 1.1: The experimental design of time-resolved ARPES, and corresponding ultrafast momentum-dependent dynamics of the melting of the charge density wave state in T bT e3 monitored by trARPES. The figure is taken from [1] 1.1 Time-resolved ARPES Since the early 1980s, time-resolved experimental techniques have been under development [16, 17, 18], and related technology keeps progressing in resolution of momentum, energy and in the time domain. Some of the important developments include generating probe photons that can access large values of electronic crystal momentum by using high-harmonic sources [19, 20, 21, 22], use of 2D hemispherical analyzers to replace time-of-flight electron spectrometers [23, 24], and aiding in selective nonequilibrium photoexcitation through the incorporation of tunable pump frequencies [6]. Time and angular-resolved photoelectron spectroscopy (trARPES) combines ARPES with the ultrahigh time resolution of femtosecond lasers. Related advanced researches in capturing the fastest processes relevant to behaviors at surfaces in materials have been explored. For the experiments, two important processes are included in time-resolved ARPES, 2 namely photoemission and ultrafast pump-probe spectroscopy. First, in a photoemission process, a crystalline material is pumped by ultraviolet (UV) photons in vacuum. When the photon energy is larger than the work function of the system, electrons will be excited and ejected into free space where their exit angles and energies can be observed. Because of the energy and momentum conservation laws, one can ensure that the energy and momentum of the photoelectrons in vacuum are directly related to the crystal momentum and energy of the same electrons when they were still inside the material. Therefore the result can be used to study electronic structure, and has been successfully applied to characterize many different materials [25, 26, 27, 28, 29]. On the other hand, even though ARPES has many advantages and gives a great contribution in the field, there are a few drawbacks. First, some microscopic dynamic processes cannot be decoupled from the basic analysis. Second, it can only be used to measure electronic bands that are occupied initially. Third, specific physical phases only occur after driving the system far beyond the equilibrium state. One of the useful techniques scientists already implemented is to incorporate pump-probe methods. By shining an ultraviolet pump laser pulse on a material, the material will be irradiated into excited states at a nonequilibrium condition and, after that, the second laser pulse, as a probe laser pulse, provides a photon energy within the range (6 to 60 eV) to detect photoemission of electrons from the nonequilibrium state. By studying variations in the photoemission spectrum as the time delay between pump and probe pulses changes, one is able to investigate nonequilibrium dynamics in quantum systems (see Fig. (1.1)). 1.1.1 Nonequilibrium Regime Using near-infrared or terahertz pump pulses to generate nonequilibrium states is one of the most important topics for time-resolved spectroscopy nowadays. In an early finding a pumpprobe scheme was used to acquire the time resolved ARPES data on the high temperature 3 superconductor Bi2212 [30] yielding its quasiparticle dynamics. The results were related to a quasi-thermal system of metallic quasiparticles [31], which was used to study the normalstate electron-phonon coupling. From Fig. (1.1), the experimental design of time-resolved ARPES is shown, and the ultrafast momentum-dependent dynamics of the melting of the charge density wave state in T bT e3 monitored by trARPES is illustrated [1]. Moreover, time-resolved ARPES can help scientists to study the phase transition induced by electron-phonon interactions that have no analogue in equilibrium ARPES. Gaps in the electronic excitation spectrum usually open up simultaneously with the formation of superconducting or charge density wave (CDW) order. With time-resolved ARPES in the ultrafast regime, many novel experiments have revealed new nonequilibrium behaviors. For instance, researchers have discovered a long-range CDW order persists in the nonequilibrium state but the local electronic spectrum becomes gapless for a transient period of time [32, 33, 34, 3, 1, 35, 36, 37, 38, 39]. Similar results for other materials indicate a novel quasiuniversal behavior in the nonequilibrium state. Understanding of the momentum-dependent structure of relaxation rates provides a unique probe of coupling effects in a variety systems; including charge density waves, superconductors, magnetic states and Floquet-Bloch states. Recently, extreme-UV (10-50 eV) pulses have been induced by high-harmonic generation in gases driven by intense laser pulses. Time-resolved ARPES experiments making use of such techniques have been applied to discover nonequilibrium behaviors in multiple materials [40, 41, 42, 43, 44]. For example, the chalcogenide based charge density wave (CDW) materials such as 1T − T aS2 , 1T − T iSe2 and RT e3 [45, 37, 46]. Scientists now can also access the timescales of charge-order collapse, and control the CDW gap through changing frequency and amplitude of the laser pulse. The delay time is of order 600 fs for collapse of the superconducting gap in Bi2212 while the experimental results on 1T − T iSe2 show that the quenching procedure can rapidly appear within 20 fs [47, 45]. Time scales occurring for various other charge density wave systems help to understand the difference between Peierls and Mott effects driving CDW formation [46]. 4 An intriguing result in experimental observations is that of oscillating spectral signatures that indicate the coupling of electrons to collective modes. This pump-induced phenomena lasts longer than what was originally expected and is a very interesting discovery; perhaps due to the fact that the electron-phonon coupling and charge density wave and lattice response is far from the near equilibrium expectation. Coherently driven oscillations appear in CDW materials as well. Early studies of these oscillations in KMO by Demsar et al. applying optical techniques [48] were replicated in a wide range of materials, such as 1T − T aS2 to T bT e3 [37, 34, 1]. 1.2 CDW Materials A CDW is an ordered quantum fluid of electrons normally in a linear chain compound or layered crystal. In 1955, the CDW state was first mentioned by Peierls, and more recently time-resolved ARPES (trARPES) has been widely applied to study several CDW materials [49]. Peierls discovered a one-dimensional metal coupled to lattice vibrations is not stable at low temperature. This instability was later called the Peierls’ CDW state. The argument first considers a one-dimensional metal with one conduction electron per lattice site and with temperature ideally at T = 0 and lattice constant a. The ground state is a non-interacting metal without the electron-electron and electron-phonon interaction. A periodic lattice distortion occurs when electron-phonon coupling is present, and will influence the Fermi wave vector and the overall energy of the material. The reason is that lattice energy is increased less than electronic energy is reduced. And, the size of the gap opening at the Fermi level is related to the amplitude of the periodic lattice distortion. In Fig. (1.2), we show the uniform electron density and lattice distortion of a one-dimensional charge density wave. Note that, in a half filling situation, a dimerization transition occurs and the distortion will double the periodicity of the model as illustrated in Fig. (1.2). More discussion can be found in Ref. ([50]). 5 Figure 1.2: Schematic picture for the comparison between uniform electron density and the Peierl’s transition of a charge density wave in a 1D system. 1.2.1 CDW in 1T − T aS2 One of the original 2D CDW materials discovered is the material 1T − T aS2 [3, 51, 32], which has a sandwich structure consisting of S-Ta-S sheets and the sheets are bound together by Van der Waals forces. Because of the weak coupling between layers, the system has a quasi two-dimensional property, and has a variety of stacking structures. Fig. (1.3)(a) shows the measurements of conductance of 1T − T aS2 with temperature. At several critical temperatures, one can easily find multiple phase transitions between phases such as the commensurate CDW phase; the incommensurate CDW phase and the metallic phase. A commensurate CDW indicates a CDW state that is periodic as an integer times the periodicity of the underlying lattice system. Note that the material has a quasi commensurate CDW phase at temperatures 180K-350K, and a commensurate CDW phase (CCDW) at temperature below 180K. Also an incommensurate CDW phase in the range 350K-550K, and a metallic phase at temperature larger than 550K. In Fig. (1.3), we also show the structure of 6 the material 1T − T aS2 that is commonly studied with pump-probe experiments at room temperature. The schematic shows the T a atom layer in the CCDW phase; consisting of 13-atom David-star clusters with the unique three T a sites a, b, and c. The black arrows indicate the displacement of the T a atoms from their original undistorted positions. Coherently driven oscillations appear in this CDW material as well. Perfetti et al. demonstrated the total photoemission signal in these phases by time resolved photo-emission measurement [52, 48]. They applied a 1.5 eV laser pulse with FWHM of 125 fs to pump the material and a 6 eV near UV pulse to probe the effect. They discovered solid evidence that the commensurate CDW 1T − T aS2 is a Mott insulator and the gap on the Fermi surface has a completely electronic origin by comparing two spectra in the metallic phase and Mott phase. In the Mott insulator phase, they found the quasi-instantaneous collapse of the Mott gap and its recovery at 20 ps time scale. The result shown in Fig. (1.4) indicates coherent oscillations in the photoemission signal with a fixed period. The corresponding frequency is the phonon mode that forms the CDW when the outer lattice sites couple with the inner lattice sites. The mode is normally called a CDW amplitude mode or a breathing mode. Figure 1.3: (a) The resistivity versus temperature in different charge-density-wave (CDW) phases of 1T − T aS2 . (b) Schematic of the T a atom layer in the CCDW phase. 13-atom David-star clusters with the three T a sites a, b, and c. The black arrows indicate the displacement of the T a atoms from their original positions. Figure is take from [2]. 7 Figure 1.4: Pump-induced gap magnitude oscillations in the CDW material 1T − T aS2 near 30K. Using a 1.5eV pump pulse and 6eV probe pulse to induce instantaneous collapse of the Mott gap and to observe its recovery. The oscillation in the band is interpreted as the CDW amplitude mode related to electron-phonon coupling. Figure is taken from [3]. 1.2.2 Optimal Laser Pulse Shaping As discovered in the previous sections, high fidelity characterization of transient excited many-body electron distributions in the ultrafast time domain is now available through a variety of pump-probe experiments. Amongst the rich variety of non-equilibrium responses observed, photo-induced phase transitions (PIPT) are particularly interesting for fundamental and applied reasons [39, 38, 53, 54, 55, 56, 57, 58, 59, 35, 60, 32, 33, 51, 3, 35, 36, 61, 60]. Important advances in optimal control theory were proposed in the twentieth century, and our objective is to extend this approach to treat trARPES as a target. Fig. (1.5) shows the scheme of pulse design, which is through a learning loop that comprises a computer, a pulse shaper, the system and feedback spectroscopy [4]. The computer offers a particular pulse shape, we shine it on the system and afterwards send probe light that analyses the effect of the shaped pulse. An evolutionary algorithm iteratively refines the shape until an optimal shape produces the desired effect. Therefore, by adjusting the laser pulse properties it is possible to tune the non-equilibrium PIPT response from the adiabatic to the non-adiabatic limits. Optimal laser pulse shaping 8 methods have not been applied to PIPT yet and have the potential to control non-equilibrium response in order to isolate selected physical phenomena, and in order to tune response for selected applications such as high speed electronics or optics. In Chapter 5 we introduce a method, based on Krotov optimal control theory [62], to direct photo-induced phase transitions (PIPT) by combining non-equilibrium models with quantum optimal control theory (QOCT). QCOT is a powerful tool based on calculating the optimal pulse shape by minimizing a physical cost function or maximizing a desired physical objective, and it has been developed within a variety of variational frameworks to obtain control sequences [63, 64, 65, 66, 67, 68, 69, 70, 71, 62, 72, 73]. Quantum optimal control methods based on the classical gradient optimization methods provide an alternative to iterative methods based on the Krotov approach. Krotov methods have been applied to the fields of Quantum computing and control of charge transfer processes, [72, 74, 75, 62, 76, 73]. The Krotov approach has several appealing advantages over the gradient methods in the following ways: First, a monotonic increase toward the objective with iteration number. Second, no requirement for a line search and faster convergence to a given target. Third, at each iteration it guarantees macrosteps of the time interval and it can reach the global maximum. To illustrate the approach, we consider PIPT in a simple model [56] for a transient metal-insulator state in a charge density wave system. In experiments, a long-range chargedensity-wave (CDW) is formed in a variety of layered chalcogenide materials, in oxides, in two dimensional materials, and in many other systems. The investigated layered CDW materials have disclosed a new nonequilibrium pattern where the long-range CDW is preserved while the local electronic excitation spectrum becomes gapless (by having subgap states) for a transient period of time, as has been elucidated using the simple model considered here [56]. In the experimental systems the mechanism for gap-closing and population inversion remains an open question, though in some cases there has been significant theoretical 9 progress. Figure 1.5: A closed-loop for an optimal control process in quantum systems. First, input initial conditions to the quantum system, such as an initial random field. Second, a current laser control field design is created with a pulse shaper and then applied to the sample. Third, The system evolves and the results are fed to a learning algorithm to suggest an improved laser field for repeated excursions around the loop until the objective is achieved. Figure is taken from [4]. 1.3 Floquet-Bloch States The light-matter interaction between the coherent electric field of an intense laser pulse and the electronic wave functions of a solid can induce intriguing phenomena such as modifications of the electronic structure and topology. In particular the finding of Floquet-Bloch states in topological insulators has become an important area recently. A topological insulator is a material that has non-trivial topological order. It behaves as an insulator in the bulk but its surface contains conducting states, indicating that electrons can only move along the surface of the material. Fig. (1.6) shows an idealized band structure for a topological insulator. The bulk band gap is traversed by topologically-protected surface states. On the other hand, Floquet-Bloch states can be considered as a time domain version of Bloch waves, where Bloch waves are defined as a periodic wave with fixed crystal momentum and 10 the corresponding Hamiltonian has a discrete translation invariance in space. Similarly, Floquet-Bloch states in a material are periodic in both energy and momentum, and the Hamiltonian has discrete invariance in both space and time. Figure 1.6: An idealized band structure for a topological insulator. The bulk band gap is traversed by topologically-protected surface states [5]. Researchers recently provided interesting experimental characterization of the effect [6] through pump-probe trARPES studies of the topological insulator Bi2 Se3 as show in Fig. (1.7). They found a ladder of replica bands appear in the time resolved photoemission spectrum with equally spaced energies after a stable frequency of the pump laser pulse was input into the system, and measured by a probe pulse within the duration of the pump pulse. Also, the surface state bands are dispersive and lead to different overlaps for some bands. With a strong enough polarized laser pulse, band gaps will open at the crossing points with induced gap proportional to the amplitude of the applied electric field; and new topological properties may be induced [77, 78]. For pulse energies increasing to very high values, the gap will shrink again and connect to the possibility of band gap opening depending on pump fluence. 11 The Floquet-Bloch phenomena is related to photon-dressed electron states in vacuum and has the potential to disclose new physics in a variety of fields. It is also possible to control these properties to design the next generation of quantum devices for new applications. Figure 1.7: Researchers recently showed experimental results of the Floquet-Bloch states [6] through pump-probe trARPES measurements on the topological insulator Bi2 Se3 . Figure is taken from [6]. Moreover, nowadays ultrafast ARPES not only can study novel nonequilibrium states but also can be used to map out the unoccupied band structure in solids by applying two photon photoemission techniques. With these important techniques, one can study dispersion of correlated quasiparticles in energy-momentum space for quantum materials and the conduction bands of nanoscale quantum systems. In topological insulators, the two photon photoemission ARPES yields a mapping of the unoccupied state band structure on the surface, as discussed in Ref. [79]. 12 1.3.1 Laser Driven Topological Phases Optically induced phase transitions of topological states of matter by ultrafast lasers is becoming an important topic in condensed matter physics. By combining both of these lively areas, recent work showed the coupling of short laser pulses to Dirac fermions in the topological insulator Bi2 Se3 as shown in Fig. (1.7) [80]; and in other recent studies [81, 82, 83]. A simple model of conducting polymers, the Su, Schrieffer, and Heeger (SSH) model, describes the dimerization occurring in polyacetylene, and is a classic example of a onedimensional topological insulator [5]. While the topologically trivial or nontrivial character of the dimerized chain is controlled by the relative strength of the nearest neighbor couplings, it was only recently proposed that high-frequency laser light can turn trivial equilibrium bands into topological nonequilibrium Floquet bands [6]. In Chapter 3, we introduce a new Zig-Zag SSH model where a laser pulse can induce novel topological behavior. 1.4 Chapters An outline of the thesis is as follows: In chapter 2, the goal is to introduce key methods and ideas used, particularly nonequilibrium Green’s functions for tight-binding models. The two main models used later in the thesis; for charge density waves and topological insulators, will be introduced. Many techniques for studying physical properties from nonequilibrium Green’s function will be mentioned, such as density of states, number density and order parameter. In chapter 3, we simulate the real-time processing of single-particle energy gaps in polyacetylene coupled to short laser pulses through a Zig-Zag SSH model, applying trivial and nontrivial topological parameters of the chain for time-resolved photoemission spectroscopy (TRPES) analysis. We find that the TRPES band structure presents well-defined Floquet 13 bands and predict a novel band formation of the non-trivial topological chain. This chapter is new work, and will lead to one first author paper and one or two co-auther papers. In chapter 4, we first explain the complete foundation and schemes to solve the control equation of motions and exhibit how to include additional constraints into the goal/cost functions. Second, we introduce the tricks of Krotov optimal control method and the corresponding backpropagator. Third, we present the monotonic convergence property and related proof. Fourth, we discuss the first order and second order cases for the backpropagator format, and the algorithm for implementation on a computing resource. Last, we will show three examples to demonstrate the power of the optimal control method. In chapter 5, we consider the hot-electron model following the simulation study of Ref. [56, 55]. They found ultrafast laser pulses can quickly heat electrons forming a hot quasithermal gas that equilibrates with phonons on much longer time scales compared to the electron relaxation time. We introduce a new QOCT method to find the possible shaped-laser pulses to study the short-time transient phenomenon and to discover new photo-induced phase transitions without melting the system. This method is applicable to general tight binding models, and we illustrate it by controlling negative temperature states in the model of Ref. [56]. In chapter 6, we summarize our progress in the field and give an outlook on promising future directions. 14 CHAPTER 2 MODELS AND METHODS As outlined in chapter 1, the use of lasers to direct photo-induced phase transitions in solids has become a very important research subject in the field of condensed matter physics. Techniques for the ultrafast studies of ferroelectricity, electrical conductivity, magnetism and superconductivity have been widely utilized for different materials. Ultrafast study of strongly correlated electron systems is one of the most interesting topic nowadays; especially for the dynamics of an order parameter’s amplitude and phase, which determines the collective behaviour of novel states emerging in complex materials. Time-resolved photoemisssion spectroscopy can decouple entangled degrees of freedom by visualizing their different responses in the time domain. Topological insulator (TI), a new quantum state of matter, is another new field that can be studied with pump-probe spectroscopy. Their surfaces and interfaces can act as a topological boundary to generate massless Dirac fermions with spinhelical textures and topological order. Investigation of the dynamics of the topology of these materials is crucial and further development of theories and computational approaches is a high priority. In last decade, nonequilibrium physics has become an important area which is widely studied; due to the observation of novel behaviours and phenomena that are not well explained by equilibrium models. Nonequlibrium problems often indicate time reverse symmetry breaking by a rapid disturbance in time, such as an interaction quench or a pulsed field. Tight binding models of electronic structure are expected to capture many of the essential features of the effects even in the nonequilibrium and ultrafast domains. Moreover, high fidelity characterization of transient excited many-body electron distributions in the ultrafast time domain is now available through a variety of pump-probe experiments. Using nonequilibrium Green’s functions, it is possible to fully calculate and understand the nonequilibrium behaviours of interest. 15 In this chapter, the goal is to introduce key methods used, and in particular nonequilibrium Green’s function for tight-binding models. The two main models studied in the thesis for charge density waves and topological insulators are also included. Many techniques for studying physical properties using nonequilibrium Green’s function will be outlined, such as density of states, number density and order parameters. The purpose for this chapter is to cover the key tools and theories we use in the later chapters. 2.1 Introduction to Nonequlibrium Green’s Functions In general, interacting quantum and classical many-body systems in nonequilibrium states can be completely represented by a set of N -particle wave functions entering the timedependent Schrödinger equation (TDSE) as i~ ∂ |ψ(t)i = H(t)|ψ(t)i ∂t (2.1) where H(t) is a time-dependent Hamiltonian and |ψ(t)i is the corresponded time-dependent wave function. Next, it is important to know how expectation values can be calculated at finite temperatures, and how they change in time when the many-body system is affected by a timedependent disturbance leading to deviations of the system properties from equilibrium. To this end, assume the system to be initially, at time t < t0 , in thermodynamic equilibrium corresponding to the time-independent Hamiltonian H0 . The expectation value of an observable O at a given time point t > t0 is then accessible through the trace with the grand canonical density operator as hO(t)i = Tr{ρOH (t)}, 16 (2.2) with e−βH0 ρ = Tr{e−βH0 } . (2.3) Note that the operator OH (t) is in the Heisenberg picture, OH (t) = U (t0 , t)OU (t, t0 ), and the time-evolution operator U is define by U (t, t0 ) = ´ − ~i tt dt̄H(t̄) 0 . Te (2.4) From the above, one can also define the imaginary time operator for temperature as U (t0 − iβ, t0 ) = e−βH0 , (2.5) and, by using the Trotter formula, the time evolution in numerical calculations becomes U (t, t0 ) = U (t, t − ∆t)U (t, t − 2∆t)...U (t0 + ∆t, t0 ), (2.6) where ∆t is a small time step. For each moment time t, we have   i∆t U (t, t − ∆t) = exp − H(t − ∆t/2) . ~ (2.7) This operator propagates the system from the initial time t0 parallel to the imaginary axis to t0 − iβ, Further, this allows us to refine expression Eq. (2.3) for the time-dependent expectation values. Inserting the definition Eq. (2.5) leads to the formula hO(t)i = Tr{U (t0 − iβ, t0 )U (t0 , t)OU (t, t0 )} . Tr{U (t0 − iβ, t0 )} (2.8) While the denominator is just the grand canonical partition function, the numerator of Eq. (2.8) gives rise to the following interpretation: Reading the time-arguments from right to left, one may say that the system first evolves along the real chronological time axis from 17 t0 to time t at which the operator O acts. Then the system anti-chronologically evolves back along this axis from time t to t0 and finally gets propagated parallel to the imaginary axis from t0 to t0 − iβ. Graphically, this leads to a time contour in the imaginary time plane which originally has been introduced by Keldysh [84] and is referred to as the KadanoffBaym-Keldysh contour as Fig (2.1) shows. The imaginary branch of this contour accounts for the ensemble averaging at the given temperature β. Figure 2.1: Kadanoff-Baym-Keldysh contour in the complex time plane. A more general expression for all time-dependent expectation values is consequently   ´ Tr{TC exp(−i C dt̄H(t̄))O(t) } hO(t)i = , (2.9) Tr{U (t0 − iβ, t0 )} where the exponential function is to be understood similarly to Eq. (2.4) as a Dyson series, and TC is now the contour time-ordering operator satisfying TC {a(t)b(t0 )} = θ(t − t0 )a(t)b(t0 ) − θ(t0 − t)b(t0 )a(t) (2.10) with θ(τ ) the step function that equals to one when τ > 0 and is equal to zero when τ < 0. Note that, operator O(t) is of the form O(t) = U (t0 , t)OU (t, t0 ). (2.11) With the definitions and considerations above one is now ready to move on to the definition of the nonequilibrium Green’s functions that we use. In the Heisenberg representation, 18 the one-particle time-ordered Green’s function is defined by G(x, x0 , t, t0 ) = −ihTC = −i h i c(x, t)c† (x0 , t0 ) i h i Tr{U (t0 − iβ, t0 ) c(x, t)c† (x0 , t0 ) } Tr{U (t0 − iβ, t0 )} (2.12) where c(x, t) and c† (x0 , t0 ) are fermion creation and annihilation operators. Through choosing t and t0 at different times on the Kadanoff-Baym-Keldysh contour, one can obtain different kinds of Green’s functions. In fact, ten different kinds of single particle Green’s function are mentioned in the literature [85, 56, 55] but, in the thesis, we will only use two main independent Green’s functions from which we are able to derive most other Green’s functions. One is the retarded Green’s function defined as GR (x, x0 , t, t0 ) = −iθ(t − t0 )h{c(x, t), c† (x0 , t0 )}+ i, (2.13) and the other is the lesser Green’s function defined as G< (x, x0 , t, t0 ) = −ihc† (x0 , t0 )c(x, t)i. (2.14) Note that here {, }+ represents the anti-commutation relation. 2.2 Introduction to Topology in Band Structure In condensed matter physics, a central topic is to identify phases of matter. Symmetry breaking is one of the main reasons for certain phase transitions, for example, in magnets and superconductors. These kind of phenomenon can be understood using ideas related to order parameter symmetry breaking. Topological insulators have a different symmetry breaking where the idea of topology or topological order is essential. In condensed matter physics, many of these ideas developed through the quantum Hall effect which has fundamental topological properties such as the quantized Hall conductivity, and the number of conducting 19 edge modes. These properties are insensitive to adiabatic changes of states and only change when the system passes through a quantum phase transition point. Even though the topological order of the quantum Hall effect has been studied for a long time, the properties of topological insulators have been discovered recently. A topological insulator has a bulk energy gap separating the highest occupied electronic band from the lowest empty band. The system may have the same energy bands but have totally distinct topological phases. One of the main objectives of topological band theory is to classify these different electronic phases. New phenomena occur when there is a spatial interface between two topologically different phases. Somewhere along the way between the two phases, the energy gap has to go to zero. To characterize these gapless interfaces is a key goal. The study of this goal brings us to the relation between boundary topological invariants and the bulk topological invariants, which we will refer to as the bulk-boundary correspondence. In this section, we give a short introduction to the relevant topics in topology, Berry phase and Chern number. 2.2.1 Topology In mathematics, topology is used to describe geometrical properties of objects that can change shape smoothly. A famous example is that a coffee mug can smoothly deform to a donut (as Fig. (2.2) shows) since they both have one hole or genus g = 1. For the condition that one can only change the shape smoothly and cannot add a hole, a donut (g = 1) cannot be the same topology as a soccer ball (g = 0) so they are topologically distinct. Manifolds that can be deformed into one another are topologically equivalent. For manifolds, a mathematical theory that defines an integer topological invariant is the Gauss-Bonnet Theorem: ˆ 1 CdA, g= 2π S 20 (2.15) Figure 2.2: A coffee mug can smoothly deform to a donut since they both have one hole or genus g = 1. with Gaussian curvature C and S is a surface. Similar concepts are applied to band theory wehre related quantities are the Berry phase and the Chern number, as discussed in the following subsection. 2.2.2 Berry Phase and Chern Number Berry phase is an important property in topological band theory [86, 87, 5] since it can be used to understand the intrinsic phase of a quantum wavefunction. We know the Bloch states are invariant under the shift |ψ(k)i → eiφ(k) |ψ(k)i. (2.16) This shift is an electromagnetic gauge transformation. We introduce the idea of the Berry connection ~ = −ihψ(k)|∇k |ψ(k)i. A (2.17) Here A can be considered as analogous the electromagnetic vector gauge transformation. Even though A is not gauge invariant, the corresponding flux in a closed loop is. One can define Berry phase for a given closed loop C in momentum space as ˛ ˆ ~ αC = A · dk = F d2 k, C S 21 (2.18) where ~ F =∇×A (2.19) defines the Berry curvature. Note the surface and contour here can be in any dimension. To be more specific, let us consider the simplest two level Hamiltonian that can be represented in terms of Pauli matrices ~σ as  ~ · ~σ =  H(k) = d(k)   dx − idy  . dx + idy −dz dz (2.20) This Hamiltonian would give us eigenvalues ±|d|. One could have an additional additive term ~ d| ~ so but that would not affect the eigenvectors that depend only on the unit vector dˆ = d/| we can skip it. dˆ here can be considered as a vector point on a sphere S 2 . If one considers a loop C, let dˆ go over a 2π rotation in a plane, then the Berry phase is found to be π [88]. ˆ We therefore can say αC is half the solid angle swept out by d(k). Now we can define the ˆ as Berry curvature as half the solid angle element for the mapping d(k) 1 ˆ F = ij dˆ · (∂i dˆ × ∂j d). 2 (2.21) From above, we have a relation that the Berry curvature integrated over a closed 2D space ˆ circles around the equator as k is a multiple of 2π that is equal to the number of times d(k) goes over a path in the Brillouin zone as Fig. (2.3) shows. We can then define a topological invariant named the Chern number [86, 87, 5] ˆ 1 N= F d2 k. 2π S (2.22) The result of quantization of the Chern number is not restricted to two band systems but also generally applies to most cases. Note that the surface integral should be equal with each other for the inside loop and outside loop up to a multiple of 2π, and quantization of 22 ˆ for a two band Figure 2.3: The Berry phase is equal to half the sold angle swept out by d(k) theory. the Dirac monopole are clearly related to the quantization, such that we can consider F as a curvature number analogous to the Gaussian curvature K mentioned earlier. 2.2.3 Bulk Edge Correspondence A topological insulator has gapped band structures in the bulk and gapless conducting states at interfaces where the topological invariant changes. The states are classified into different topological orders based on their fundamental properties. To understand the edge states, it is easy to consider the interface between the integer quantum Hall state and vacuum [14]. In a semiclassical framework, these edge states can be understood in term of the skipping orbits where electrons have cyclotron orbits that bounce off the edge (Fig. (2.4)). Note that, the electronic states responsible for this behavior can propagate in two different directions along the edge so they are chiral. There are no states available for backscattering which underlies the perfectly quantized electronic transport in the quantum Hall effect, the 23 Figure 2.4: Edge states as skipping cyclotron orbits. states are insensitive to disorder effects. The topology of the bulk quantum Hall state then leads to the chiral surface states and quantized conductance. A simple explanation of the chiral edge states can be shown by applying the two band Dirac model, where at a Dirac point, HD = mσz , (2.23) with σz is a Pauli matrix. Now consider a mass m is a function of y as m(y) at one of the Dirac points of an interface and will change sign with y. For y > 0, m(y) > 0 gives the insulator, and, for y < 0, m(y) < 0 gives the quantum Hall state. Assume m0 > 0 is fixed. Since the states satisfy translation symmetry in the x direction, one can write ψqx = eiqx x ψ(y). The zero energy mode ψ0 (y) is of the form ´ y 0 ψ0 (y) = e− 0 m(y )/vF | ↑i 24 (2.24) where | ↑i is the upper eigenvector of σz . The corresponding eigenenergy is E(qx ) = ~vF qx . (2.25) The band of states has a positive group velocity dE/dqx = ~vF defining a right moving chiral edge mode that intersects with the Fermi energy EF . Figure 2.5: The chiral edge states connect the valence band near K and K 0 with linear or nonlinear dispersion for linear and nonlinear Hamiltonians respectively. In Fig. (2.5), the energy bands versus momentum kx along the edge is shown. The blue areas indicate the bulk conduction and valence bands, which have continuum states and have the energy gap near K and K 0 . Interestingly, a line connecting the valence band at K 0 to the conduction band at K with a positive group velocity describes edge states. The edge states may be different if the Hamiltonian near the surface is described by a non-linear function. For example, the red line in Fig. (2.5) shows the edge states can pass the Fermi level three times if the group velocity is non-linear. However the difference NR − NL between the number of right and left chiral modes is the same, and is an integer topological invariant characterizing the interface. The edge related value of NR − NL is therefore decided by the topological structure of the bulk states. This is summarized by the 25 bulk edge correspondence: ∆n = NR − NL . (2.26) Note that ∆n is the difference in the Chern number across the interface. 2.3 Tight-binding Models at Equilibrium In solid-state physics, tight-binding models (or TB models) are used to compute electronic band structure, using an approximate approach based upon superposition of localized wave functions at each atomic site. The model is closely connected to the LCAO (Linear Combination of Atomic Orbitals) method used in chemistry. Since the electrons in the model are quite tightly bound to the atoms and have limited overlap with surrounding atoms of the solid, the model is particularly well-suited for calculations of solids for which the electronic states remain relatively close to the unperturbed atomic orbitals. The energy of the electrons will also be quite close to the ionization energy of the electron in the free atom or ion because the interaction with potentials and states on neighboring atom is limited. Though complete physical principles are involved, tight-binding models yield a relatively simple mathematical formulation. 2.3.1 SSH, Rice Mele, Platero Models In this subsection, we introduce an important one dimensional model that illustrates some of the fundamental topological properties in solid state physics. The simple but elegant Su Schrieffer Heegar (SSH) model was used as a basic model of the conducting polymer polyacetylene [5]. At half filling it has a Peierls instability to a dimerized state. The most intriguing part is that in a finite chain there are two different dimerized states, as shown in Fig. (2.6). These two distinct topological states lead to different topological properties and responses in the nonequilibrium situation. The SSH model here presents the simplest two 26 band model for introducing these special topological phenomena. Figure 2.6: Two distinct topological states of the SSH model. (a) is the case δτ > 0 and (b) is the case δτ < 0. SSH used a 1D tight-binding Hamiltonian to model polyacetylene as H= X † † (τ + δτ )cAi cBi + (τ − δτ )cAi+1 cBi + h.c. (2.27) i with labels A and B referencing to the two different atoms/sites in a unit cell. The term δτ represents the dimerization and induces a gap opening in the energy band. The operators satisfy the anti-commutation relations † {ci , cj }+ = δij , {ci , cj }+ = 0, † † {ci , cj }+ = 0, (2.28) Note that, δτ can be positive or negative and we would have two distinct patterns or topologies determined by the sign of δτ . Even though the real SSH model should consider spin, we treat spinless electrons here for simplicity. To analyze the Hamiltonian in momentum space through Fourier transform we impose periodic boundary conditions for a chain with an even number of sites. In that case we get a two-by-two matrix for each momentum k, and we 27 define N as the number of sites: H= X † Hab (k)cak cbk , (2.29) k with ~ · ~σ H(k) = d(k) (2.30) and dx (k) = (τ + δτ ) + (τ − δτ ) cos(ka), dy (k) = (τ − δτ ) sin(ka), dz (k) = 0. (2.31) Note that ~σ = {σx , σy , σz } and  0 σx =  1  0 σy =  i  1 σz =  0  1  0  −i  0  0  −1 (2.32) (2.33) (2.34) P ~ so the inner product with d(k) is H(k) = i={x,y,z} di (k)σi . From Eq. (2.31), one writes ~ the two band H(k) in terms of d(k) enabling use of topological ideas introduced in section 2.2.1. An important point is that if dz = 0; the two bands have a chiral symmetry defined by the operator Π = σz , which anticommutes with the Hamiltonian: {H(k), Π} = 0. Then each eigenstate with energy E has a corresponding state |ψE i = Π|ψ−E i with energy −E, so that this chiral symmetry actually indicates a particle-hole symmetric spectrum. While this symmetry is not preserved in real polyacetylene because there are second order effects such as second neighbor hopping; it is still very useful to investigate the effects and topologies for understanding some outstanding physical phenomena. 28 Consider the different signs of δτ , which can lead to different topological phases or geo~ metric results. For δτ > 0, then dx > 0 for all k so d(k) has zero phase if d~0 = ~0 is the center (degeneracy point). One can see this path from Fig. (2.7). On the other hand, for δτ < 0, ~ dx (k ≈ π/a) < 0, so that d(k) circles around the center and has phase 2π. The phase results can be understood as polarization in the strong coupling region, |δτ | ≈ τ . In that limit, one can think that the electrons stay in localized states so, when we shift δτ from positive to negative, the states will move electrons over by half a unit cell, inducing polarization and topological phase. Figure 2.7: Two distinct topological paths of the SSH model. (a) is the case δτ > 0 and (b) is the case δτ < 0. 29 The polarization gains phase of integer n times phase 2π when the chiral symmetry is kept. One can break the symmetry by considering a dz (k) > 0 so the polarization now can vary continuously and topology is lost in this 1d model. Otherwise, keeping chiral symmetry leads to topologically distinct states that are characteristic by their quantized polarization. Note that, in order to get from the positive δτ > 0 state to the negative δτ < 0 state we need to pass a point with d~ = ~0 to obey the chiral symmetry, which is a quantum phase transition point. For the SSH model the eigenenergies are derived from the Hamiltonian   0 R + re−ika  ~ · ~σ =  H(k) = d(k) ,  ika R + re 0 where R = τ + δτ and r = τ − δτ ; so that, q Ek = ± R2 + r2 + 2Rr cos(ka). (2.35) (2.36) Note that, Ek is invariant under r ↔ R, but the topology changes. Here we use R and r to show a relation between the topological phase and the dispersion relation in Eq. (2.36). We also know that the length of R and r are decided by the dimerization term δτ . The geometric paths for the dispersion relation and the topological phase can be found in Fig. (2.8) and Fig. (2.9). Note that when R = 0 or r = 0, we will have a flat band for dispersion relation, and when |δτ | is the same, the dispersion relation will be the same no matter for R > r or R < r. Another simple but significant model is the Rice-Mele model [89], which is the SSH model with an extra staggered onsite potential. The Hamiltonian of the Rice-Mele model on a 1D chain is written as H= N X † † † † u(cAi cAi − cBi cBi ) + vcAi cBi + wcAi+1 cBi + h.c. i 30 (2.37) Energy / τ 2 1 4δτ = 2|R − r| 0 4τ = 2(R + r) −1 −2 -1.0 -0.5 0.0 k(a/π) 0.5 1.0 Figure 2.8: Dispersion relation of the SSH model based for R = τ + δτ and r = τ − δτ . Figure 2.9: Topological path of the SSH model depends on R, r and ka. with labels A and B referencing to the two different atoms/sites in a unit cell, and the onsite potential u, intercell hopping amplitude w, and the intracell hopping amplitude v. Note that, for the SSH model, u = 0, v = τ + δτ and w = τ − δτ . For the Rice-Mele model, the matrix of the Hamiltonian on a 1D chain with N = 4 sites is of the form 31   u v 0 0 0 0 0 0     v −u w 0 0 0 0 0      0 w u v 0 0 0 0         0 0 v −u w 0 0 0  . H =    0 0 0 w u v 0 0       0 0 0 0 v −u v 0        0 0 0 0 0 w u v    0 0 0 0 0 0 v −u (2.38) To analyze the Hamiltonian in momentum space through Fourier transform we impose periodic boundary conditions for a chain with an even number of sites. In that case we get a two-by-two matrix for each momentum k, and we define N as the number of sites: H= X † Hab (k)cak cbk , (2.39) k with ~ · ~σ H(k) = d(k) (2.40) and dx (k) = v + w cos(ka), dy (k) = w sin(ka), dz (k) = u. Note that ~σ = {σx , σy , σz } are Pauli matrices. The eigenvalues are, q Ek = u ± v 2 + w2 + 2vw cos(ka). (2.41) (2.42) The bulk energy eigenstates of a band insulator are delocalized over the whole system. We use as an example the bulk Hamiltonian of the Rice-Mele model. The energy eigenstates are the plane wave Bloch states |Ψ(k)i = |ki ⊗ |u(k)i, 32 (2.43) with N 1 X imk √ |ki = e |mi, N m=1 (2.44) for k = {δk , 2δk , ..., N δk } with δk = 2π/N . The |u(k)i are eigenstates of the bulk momentumspace Hamiltonian from Eq. (2.40). Note that, the Bloch states |Ψ(k)i are spread over the whole chain. They span the occupied subspace, defined by the projector P = X |Ψ(k)ihΨ(k)|. (2.45) k∈BZ The phase of each Bloch eigenstate |Ψ(k)i can be set at will. A change of these phases, a gauge transformation, |u(k)i → eiα(k) |u(k)i gives an equally good set of Bloch states, with an arbitrary set of phases α(k) ∈ R for k = δk , 2δk , .... Because of this freedom, we are able to ensure that in the thermodynamic limit of N → ∞, the components of |Ψ(k)i are smooth and continuous functions of k. While this gauge may not be easy to derive; it is topologically trivial when u 6= 0 since a circle loop does not enclose degeneracy point. Another interesting extension of the SSH model is a dimers chain coupled to an ac electric field discussed by Platero [81] and illustrated in Fig. (2.10). τ and τ 0 are hopping terms between two sublattices A and B, and with a periodic boundary conditions. a0 is the cell size and b0 is the distance between A-B atoms. Figure 2.10: A dimers chain with two sublattices A and B, and τ and τ 0 are hopping terms between two sublattices A and B. a0 is the cell size and b0 is the distance between A-B atoms. 33 The corresponding bulk Hamiltonian in momentum space is therefore of the form   ρF (k)  0 τ̃n,m = τ  (2.46) , ρ̃F (k) 0 with ρF (k) = λe−ikb0 + eik(a0 −b0 ) , (2.47) ρ̃F (k) = ρ̃∗F (k) (2.48) and where λ = τ 0 /τ . 2.3.2 Freerick’s Model A simple model of an electronic CDW in a 2-D system is shown in Fig. (2.11) [56, 55]. Note that, this is a special case of the Rice-Mele model described in the previous subsection on a square lattice. In the model, the charge density wave system is constituted by two different sublattices A and B of equal size. The nearest neighbors for lattice sites in the A sublattice is from the B sublattice and the nearest neighbors for lattice sites in the B sublattice is from the A sublattice. To describe different sublattices, one can use the fixed site potential energy, which is described by an onsite potential U on the sublattice A and onsite potential zero on the sublattice B. Hopping term τ is allowed between nearest sublattice A and sublattice B and is the same for all directions. Note that here we use spinless electrons for simplicity. To start with, let consider the Hamiltonian without an external field so the model now is time-independent and in a second quantization form is H=− X † τij ci cj + (U − µ) ij X † X † ci ci − µ ci ci , i∈A 34 i∈B (2.49) Figure 2.11: CDW model with two different sublattices A and B. † Note that here ci and ci are the creation, annihilation operators at site i for spinless electrons, and µ is the chemical potential. To study the electrons in a periodic lattice, we work in Fourier space and analyze it within the first Brillouin zone. We define the lattice vector P R = i ni ai êi where ni is an integer, ai is the lattice constant, and êi is the basis vector direction in the lattice. We know any periodic function on the lattice R has nonzero Fourier elements only for certain vectors {K}. In momentum space, we know these vectors {K} form a reciprocal lattice. The sites of a reciprocal lattice represent wave vectors that satisfy the condition eiK·r = eiK·(R+r) for any real space vector r. Similar to the way one chooses the primitive unit cell in real space. The first Brillouin zone is decided by the WignerSeitz primitive unit cell centered at the origin in reciprocal space. Note that we apply the reduced zone scheme by reducing all wave vectors into first Brillouin zone so energy levels originating from k + nK are now regarded as belonging to the nth band. We therefore have the transformation from real space to momentum space in the reduce Brillouin zone scheme 35 satisfying † ci = X † e−ik·Ri ck , (2.50) k with k going over the first Brillouin zone, and Ri is a lattice vector for site i. Since the translational symmetry is partially broken in the model so the momentum points k and k + Q are coupled, with Q = (π, π, π, ...). Note that this holds true for hypercubic lattices in general dimension, In the reduced Brillouin zone, the transformation equation becomes † ci = X † † (e−ik·Ri ck + e−i(k+Q)·Ri ck+Q ). (2.51) k Because e−iQ·R is one for lattice sites on the A sublattices and minus one on the B sublattices, we can find; † † † eik·Ri (ck + ck+Q ), (2.52) X ik·R † j (c − c† = e k+Q ). k (2.53) ci∈A = X k † ci∈B k By taking the Hermitian conjugate, the annihilation operator identities are derived. One can show the electronic band structure at U = 0 is then; εk = − X hiji d d X X   t∗ l √ cos(k l a), τij exp −ik · (Ri − Rj ) = −2τij cos(k a) = − d l=1 l=1 (2.54) then restricting k to the reduced Brillouin zone, is equivalent to having εk ≤ 0, and l is the index for spatial component, d is the number of spatial dimensions, a is the lattice constant, t∗ is the re-defined hopping energy between nearest neighbors in each dimension √ and τij → t∗ /2 d. The rescaled hopping term t∗ is the rescaled hopping energy between nearest neighbors in each dimension to preventing divergence in large dimensions. To be √ more specific, tij = t∗ /2 d in the last equation allows us to approach the d → ∞ limit. Note that the reduced Brillouin zone is half the size of the original one (without the sublattice modulation) since, due to the onsite chemical potential, the periodic size in real 36 space is doubled. Now one can write down the Hamiltonian in a two by two [ck , ck+Q ] basis, H= X k    U/2 k + U/2 − µ   ck  † †  . ck ck+Q  U/2 −k + U/2 − µ ck+Q  (2.55) If the system is in equilibrium, one is able to diagonalize the Hamiltonian matrix with the eigenfunction basis, ck+ = αk ck + βk ck+Q , (2.56) ck− = βk ck − αk ck+Q . (2.57) For the upper band and lower band eigenstates at the momentum point k in the reduced Brillouin zone, we know that ck+ and ck− are annihilation operators with αk and βk given by q αk = (U/2)/ 2(2k + U 2 /4 − 2k + U 2 /4) q q βk = ( 2k + U 2 /4 − k )/ 2(2k + U 2 /4 − 2k + U 2 /4) (2.58) (2.59) Note that the Hamiltonian matrix is diagonalized so we are able to rewrite the Hamiltonian in the two by two [ck+ , ck− ] basis as    k+ 0  ck+  † Hk = c†  .  c k+ k− 0 k− ck−   (2.60) Where εk+ and εk− are the eigenstates; εk± = U/2 − µ ± q ε2k + U 2 /4. (2.61) Fig. (2.12) shows the corresponding dispersion relation εk± for the CDW model. Note that, in the Pauli matrix representations d~ · ~σ , it is clear that the Chern number is always zero and hence without light the model is not topological. The corresponding dispersion relation for the CDW model is given in Fig. (2.11). 37 15 10 εk ± 5 0 −5 −10 −15 −10 −5 0 5 εk 10 Figure 2.12: Dispersion relation εk± of the CDW model when µ = 1/2 and U = 1. 2.3.3 Green’s Function Method We calculate the two-time retarded Green’s function and two-time lesser Green’s function defined in subsection 2.1. Because the time-translation symmetry is not broken in the equilibrium system, one only need to consider the relative time t − t0 in the Green’s functions. The local density of states is found from the local retarded Green’s function and defined through ˆ 0 1 1 R Ai (ω) = − Im[Gi (ω)] = − Im( eiω(t−t ) )GR (t, t0 )d(t − t0 )). i π π (2.62) The two-time retarded Green’s function is represented using the creation and annihilation operators for sublattices A and B as a sum over all k points; 0 0 GR i (t, t ) = −iθ(t − t ) X † † h{ck (t) ± ck+Q (t), ck (t0 ) ± ck+Q (t0 )}+ i, (2.63) k when in equilibrium, Ô(t) = eiHt Oe−iHt for any operator. Note that the + sign represents a local retarded Green’s function on the A sublattice and the − sign represents that on the B sublattice. One can re-write the Green’s function in the diagonalized basis so the local 38 retarded Green’s function on the A sublattice becomes 0 0 GR A (t, t ) = −iθ(t − t ) X h{(αk + βk )ck+ (t) + (βk − αk )ck− (t), k (αk + βk )ck+ (t0 ) + (βk − αk )ck− (t0 )}i (2.64) By solving equation of motion for the time evolution of the ck+ (t) and ck− (t) operators, one finds; † † † † ck+ (t) = exp[−ik+ t]ck+ (2.65) ck− (t) = exp[−ik− t]ck− (2.66) The local density of states (LDOS) is then of the following form # q "s ω ± U/2 ρ( ω 2 − U 2 /4). AA,B (ω) = Re ω ∓ U/2 (2.67) Through summing over the local DOS for each sublattice with weight 1/2, the total DOS is calculated. Note that the two sublattice model at equilibrium gives a two band structure with a band gap equal to the onsite interaction U . We also know the LDOS on the A sublattice has a divergence at ω = U/2 and the LDOS on the B sublattice has a divergence at ω = −U/2. Just like the retarded Green’s function, the imaginary part of the Fourier transform of the lesser Green’s function gives 1 1 (ω) = Ai (ω)f (ω) = ImG< Im i 2π 2π ˆ  0 0 0 eiω(t−t ) G< i (t, t )d(t − t ) . (2.68) In our calculation, we set chemical potential µ = 0 which is the half-filling case. That is the reason why the imaginary part of the lesser Green’s function indicates that the electrons fill the lower band. By definition, one can find that imaginary part of local lesser Green’s function corresponds to the number of local electrons at time t0 so the local electron number density on each sublattice can be written as nA,B (t) = Im(G< A,B )(t, t)). 39 (2.69) One can also define the order parameter at equilibrium, which is the difference between the electron number density on the A and B sublattices as Ω(t) = nB (t) − nA (t) . nB (t) + nA (t) (2.70) In equilibrium there are always more electrons on the B sublattice than on the A sublattice as there is a repulsive potential on the A sublattice. In the half filling case, the upper band is completely empty and the lower band is filled at T = 0. So we can conclude that the larger the interaction U , the stronger the electronic CDW order is. 2.3.4 Graphene and Dirac Cone The term graphene first appeared in 1987 to describe single sheets of graphite as a constituent of graphite intercalation compounds. Graphene also is the most typical 2D material with Dirac cones [12]. It has two C atoms per unit cell arranged in a hexagonal lattice (Fig. (2.13)). Figure 2.13: Graphene is an atomic-scale hexagonal lattice made of carbon atoms. 40 Because the C atoms bind together with both σ and π bonds, graphene is rather stable and can be studied through a relatively simple model. A tight-binding model describes the electronic structure of graphene quite well. The Hamiltonian is simplified into a 2 × 2 matrix for each wave vector k using a nearest-neighbor hopping approximation as     X  0 ti ei(k·di )  ε0 0  H(k) =  +  , −i(k·d ) i 0 ε0 0 i=1,2,3 ti e (2.71) where di with (i = 1, 2, 3) are vectors that connect a C atom to its three nearest neighbors, ε0 is the on-site energy, and ti are the corresponding hopping energies. The energy bands are thus found to be E± (k) = ε0 ± | X ti ei(k·di ) |. (2.72) i=1,2,3 For convenience, one can set the Fermi level to be ε0 = 0. For the equilibrium structure, di = r0 , ti = t0 , and the valence and conduction bands contact at K and K 0 points of the hexagonal Brillouin zone. One can expand the energy bands around K or K 0 as E± (q) = ±~vF |q| (2.73) where K = k − q, and vF = 3t0 r0 /2~ is the Fermi velocity. Near K and K 0 points, Eq. (2.73) shows that graphene has a cone-like band structure with linear dispersion, similar to a relativistic particle. The density of states (DOS) per unit cell with a degeneracy of 4 (2 × 2 spin and valley) near Fermi level is √ ρ(E) = 4|E|/ 3πt20 (2.74) Therefore graphene is a gapless semiconductor with zero DOS per site at Fermi level. The Hamiltonian near the K and K 0 points can be transformed into H(q) = vF p · σ (2.75) where σ are the Pauli matrices, and p = −i~∇ is the momentum operator. Note that, Eq. (2.75) is identical to, by replacing c (speed of light) with vF , the massless Dirac equation or Dirac-Weyl equation with spin S = 1/2. For this reason, the K and K 0 points are 41 also called Dirac points, and the linear band structure is called a Dirac cone as Fig. (2.14) shows. These graphene possesses massless Dirac fermions with pseudospins of ±1/2. Here pseudospins indicate which sublattice has a great occupancy. Figure 2.14: Dirac cone–the electrons in graphene behave like massless Dirac particles that appear in the electronic band structure as gapless excitations with a linear dispersion. Due to the Dirac-cone structure, many new electronic properties and physical phenomena have been found in graphene. For instance, when a uniform perpendicular magnetic field B is applied, a particular Landau level form in graphene is q 2 |n|B E(n) = sgn(n) 2e~vF (2.76) where n ∈ Z (the Landau level index) and E(n) ∝ n in normal semiconductors. The Landau levels in 2D systems can be verified by quantum Hall effect (QHE) if the magnetic field is strong enough. Different from the conventional 2D systems, graphene shows a novel 42 half-integer QHE with σxy = (4e2 /h)(N + 1/2) (2.77) which can be described by the pseuospins and the degeneracy of Dirac fermions. 2.4 Tight Binding Models in a Time Varying Field In the past decade, since precise nonequilibrium study is now possible in experiments, extensions of theory to treat pump-probe experiments has become very important. Normally, for the case of pump-probe experiments, the Hamiltonian will be time-dependent due to the effect of a time-varying electric field. We thus need to include electromagnetic fields into the lattice models described earlier, and the way that we do this is described in the next section. 2.4.1 Peierl’s Substitution To include the time-dependent electric field, we use the Peierls’ substitution [90, 56, 55], which is a widely employed approximation for describing tight-binding models in the presence of a slowly changing magnetic vector potential. It is a simplified semi-classical treatment of the electromagnetic field that is nonperturbative. The hopping matrix gains a phase factor with the Peierls’ substitution as # ˆ ie Rj τij → τij exp − A(r, t)dr . ~c Ri " (2.78) The Peierls phase originates from the propagator of an electron in a magnetic field due to the dynamical term qv · A appearing in the Lagrangian. In the path integral formalism, which generalizes the action principle of classical mechanics from the Hamiltonian H0 = 2 1 p − eA + ˜0 , 2m 43 (2.79) and one can derive the Peierls’ phase from the path integral. A scalar potential term P † −e i φ(ri , t)ci ci is also added into the Hamiltonian. φ(ri , t) is the scalar potential for the external field and A(r, t) is the vector potential. In the following, we only work with the gauge that has zero scalar potential and we also assume a spatially uniform time-dependent vector potential so one can neglect the magnetic field effects. The is called the Hamiltonian gauge. In Maxwell’s equations, one can find the corresponding electric field E(t) from the derivative of the vector potential A(r, t) as E(t) = − 1 ∂A(r, t) . c ∂t (2.80) Note that for a spatially uniform field, we have A(t) = A(t)(1, 1, ..., 1) for a vector potential directed along the (1, 1, ..., 1) axis of a hypercubic lattice. In momentum space, the timedependent band structure for the U = 0 case is then i h X e AB εk (t) = − τij exp −i(k − A(t)) · (RiA − RjB ) . ~c (2.81) hiji Therefore, the influence of the Peierls’ substitution is to induce a time-dependent shift to the momentum in the non-interacting electron band structure. At U = 0 we have,    d X t∗ eA(t) l √ cos a k − . εk (t) = − ~c d (2.82) l=1 The time-dependent Hamiltonian for the CDW case becomes (U 6= 0),   X ie † H(t) = − τij exp A(t) · (Ri − Rj ) ci (t)cj (t) ~c hiji X † X † +(U − µ) ci (t)ci (t) − µ ci (t)ci (t). (2.83) i∈B i∈A The Fourier transformation to momentum space can be applied to the above, † ci (t) = X † † [e−ik·Ri ck (t) + e−i(k+Q)·Ri ck+Q (t)], (2.84) k and the corresponding hermitian conjugate. In momentum space and in the Schrödinger representation, the time-dependent Hamiltonian can be represented as      X U/2 k (t) + U/2 − µ   ck  † † Hs (t) =   ck ck+Q  u/2 − (t) + U/2 − µ c k k k+Q 44 (2.85) The band structure εk (t) at U = 0 in time-dependent form can be written with the difference formula of the cosine, εk (t) = cos(eaA(t)/~c)ε(k) + sin(eaA(t)/~c)ε̄(k) (2.86) which depends on the band structure at U = 0 d X t∗ √ cos(ak l ) ε(k) = − d l (2.87) d X t∗ √ sin(ak l ). ε̄(k) = − d l (2.88) and One can consider the equation of motion for the operators ck (t) and ck+Q (t) in the Heisenberg picture, as ∂ck (t) = [HH (t), ck (t)], ∂t (2.89) ∂ck+Q (t) = [HH (t), ck+Q (t)], ∂t (2.90) i and i Note that HH (t) indicates the Heisenberg representation for the Hamiltonian. If one brings in the time-dependent Hamiltonian and evaluates the commutators, we find ∂ck (t) = (εk (t) + U/2 − µ)ck (t) + U/2ck+Q (t), ∂t (2.91) ∂ck+Q (t) = U/2ck (t) + (−εk (t) + U/2 − µ)ck+Q (t), ∂t (2.92) −i and −i Then the time evolution for the creation and annihilation operators satisfy      ck (t)   ck (t0 )    = Uk (t, t0 )   ck+Q (t) ck+Q (t0 ) 45 (2.93) Note that for each momentum the time-evolution operator Uk (t, t0 ) is a time ordered product      ˆ   t U/2 k (t) + U/2 − µ  0 0 dt  (2.94) Uk (t, t ) = T exp i  .     t0 U/2 − (t) + U/2 − µ k Numerically, one needs to slice the continuous time-dependent terms into small pieces for very short time steps ∆t to calculate the time evolution for Uk (t, t0 ) with the time-dependent terms inside the exponential by changing t̄ −→ t − ∆/2. (2.95) The corresponding result is more concise by employing an identity of the exponential of the Pauli matrices {σx , σy , σz }. Since any two by two matrix has four degrees of freedom, we know any two by two matrix can be represented as a linear combination of the Pauli matrices σx , σy , σz and the unit matrix 1. So we expand    a0 + az ax − iay  A = a0 1 + ~a · ~σ =  . ax + iay a0 − az (2.96) If we consider the case that the matrix A has a0 = 0, and define A = λ~a · ~σ , we can prove the generalized Euler identity relation eiλ~a·~σ = cos(λ)1 + i(~a · ~σ ) sin(λ). (2.97) Here ~a is chosen to be a unit vector and λ is the magnitude. A matrix exponential has its Taylor expansion in general as eA ∞ X An = . n! (2.98) n=0 By separating the odd and even orders, we have   ! n n X X λ λ eiλ~a·~σ = (i)n (a · σ)n + i  (i)n−1 (a · σ)n−1  (~a · ~σ ). n! n! n∈even (2.99) n∈odd One also has the useful identity (~a · ~σ )(~b · ~σ ) = ~a · ~b1 + i~σ · (~a × ~b), 46 (2.100) which includes the special case (~a · ~σ )2 = 1 (2.101) with a unit vector ~a = ~b. We see that the expansion of the cosine is equal to the even terms in the expansion, and the expansion of sine is equal to the odd terms in the expansion, hence eiλ~a·~σ = cos(λ)1 + i sin(λ)(~a · ~σ ). (2.102) In this thesis, we only work at half filling where µ = U/2 in momentum space k so the above equations become Uk (t, t − ∆t) = eiλk a~k ·~σ = cos(λk )1 + i sin(λk )(~ak · ~σ ) (2.103) with   U/2 ∆t k (t − ∆t/2)  a~k · ~σ =  , λk U/2 −k (t − ∆t/2) (2.104) and ~ak = (0, U ∆t/2λ, ∆tk (t − ∆t/2)/λ) q λk = ∆t 2k (t − ∆/2) + U 2 /4. (2.105) (2.106) In the real calculations, instead of the time at −∞, one must start from a minimum time t0 so one can calculate the time-evolution operator in the form Uk (t, t0 ) = Uk (t, t − ∆t)Uk (t − ∆t, t − 2∆t)...Uk (t + ∆t, t0 ) (2.107) For each k, one can derive the two-time evolution operator from the relation Uk (t, t0 ) = Uk (t, t0 )Uk (t0 , t0 ). (2.108) Once the time evolution at each time step is derived, we are able to calculate the nonequilibrium Green’s functions and obtain the physical properties of the model. 47 2.5 Introduction to Floquet Theory Floquet states can be considered as a time domain version of Bloch waves, where Bloch waves are defined crystal momentum and the corresponding Hamiltonian has discretized translation invariance in space. Similarly, Floquet states in a material are periodic in energy, and the Hamiltonian has discretized invariance in time [91]. 2.5.1 Periodic Hamiltonian Assume that time-dependent Hamiltonian Ĥ(t) can be written as Ĥ(t) = f (t) + H(t), (2.109) where the scalar function f (t) depends only on time t so that it commutes with H(t). We use |ψ, ti as the time dependent state vectors. One thus can reduce the time dependence from the problem i~ ∂ |Ψ, ti = Ĥ(t)|Ψ, ti ∂t (2.110) i~ ∂ |ψ, ti = H(t)|ψ, ti. ∂t (2.111) to Where the remaining non-trivial Hamiltonian H is a periodic function of time, i.e. H(t + T0 ) = H(t) (2.112) where T0 is the period time. The related angular frequency Ω is then Ω = 2π/T0 . (2.113) This periodicity allows us to rewrite H in the Fourier series as H(t) = inf X H n einΩt , n=− inf 48 (2.114) where the Fourier components H n are of the form Hn 2.5.2 ˆ T 0 1 = dte−inΩt H(t). T0 0 (2.115) Floquet Formalism Floquet theory claims that the solutions of the time dependent Schrödinger equation i~ ∂ |ψ, ti = H(t)|ψ, ti ∂t (2.116) can be written as i − t |ψ, ti = e ~ XX n Fnα einΩt |αi. (2.117) α Note that the set {|αi} is a complete orthonormal basis for the Hilbert space of the periodic Hamiltonian H. To verify the claim, we substitute the expansion into the Schrödinger equation and derive " # X ∂ ∂ − ~i t i~ |ψ, ti = i~ e Fnα einΩt |αi ∂t ∂t nα X X i − t − i t = e ~ Fnα einΩt |αi + i~e ~ Fnα inΩeinΩt |αi = nα X i − t e ~ Fnα einΩt ( − n~Ω)|αi nα nα " = X H m eimΩt − ~i t e m # X Fnα einΩt nα i − t = e ~ X Fnα ei(n+m)Ωt H m |αi. (2.118) nmα One can make a summation index substitution as i − t e ~ X i − t Fnα ei(n+m)Ωt H m |αi −→ e ~ nmα X nmα 49 Fmα einΩt H (n−m) |αi (2.119) which leads to, i − t e ~ X i − t Fnα einΩt ( − n~Ω)|αi = e ~ nα X Fmα einΩt H (n−m) |αi. (2.120) nmα Note that the exponentials {einΩt } are an orthogonal set on the time interval [0, T0 ], and that the basis {|αi} was considered to be an orthonormal one so one can have the equality X Fnβ ( − n~Ω)|βi = β X Fmβ H (n−m) |βi. (2.121) mβ By taking the scalar product from the left with hα| one will have ( − n~Ω)Fnα = X hα|H (n−m) |βiFmβ . (2.122) mβ n as We can define the matrix elements Hαβ n = hα|H n |βi, Hαβ (2.123) and the matrix elements Γnα,mβ become n−m + n~Ωδ Γnα,mβ = Hαβ nm δαβ . (2.124) One can then rewrite this equation as the eigenvalue problem X ΓFmβ = Fnα , (2.125) mβ and the matrix form is ΓF = F. (2.126) By using the property −n = hβ|H −n |αi Hβα ˆ T 0 inΩt 1 e hβ|H(t)|αidt = T0 0 ˆ T 0 inΩt 1 ∗ = e hα|H † (t)|βi dt T #∗ " 0 0ˆ T0 1 = e−inΩt hα|H(t)|βidt T0 0 n ∗ = Hαβ (2.127) 50 it is easy to derive that Γnα,mβ = Γ∗mβ,nα and therefore the matrix Γ is Hermitian. (2.128) Because of the Hermiticity of Γ, the eigen- value problem ΓF = F has a solution. One also knows the eigenvalues  are real and the eigenvectors F form a complete orthonormal basis. Therefore the assumption |ψ, ti = − i t P inΩt |αi solves the periodic time dependent problem with Floquet theory. In e ~ nα Fnα e Chapter 3, Floquet theory is applied to the SSH model. An example we show here is the effect of adding a time-dependent field into the Platero Model mentioned in section 2.3.1. The time-dependent electric field is taken to be E(t) = −∂t A(t) with A(t) = A0 sin(ωt). At steady state the AC field can be interpreted as adding an extra dimension to the model. The corresponding bulk Hamiltonian in momentum space is therefore of the form   ρF (k)  0 τ̃n,m = τ  , ρ̃F (k) 0 (2.129) ρF (k) = λe−ikb0 Jn−m (A0 b0 ) + eik(a0 −b0 ) Jm−n (A0 (a0 − b0 )), (2.130) ρ̃F (k) = λeikb0 Jm−n (A0 b0 ) + e−ik(a0 −b0 ) Jn−m (A0 (a0 − b0 )), (2.131) with and where λ = τ 0 /τ and (n, m) are integers. Note that, the spectrum depends on the intra-dimer distance b0 , and the hopping terms are normalized by the field amplitude. In the high frequency limit ω > τ, τ 0 with the chosen band n = m = 0, the Hamiltonian is block diagonal and can be described by a time independent 2 × 2 matrix Hk0 = τ~g (k) · ~σ 51 (2.132) with ~g (k) = (Re[ρ̃F ], Im[ρ̃F ], 0) for n = m = 0, and ~σ = (σx , σy , σz ) are the Pauli matrices. Since Hk0 here is similar to the SSH model and has chiral symmetry, one can calculate the corresponding Chern number (winding number) directly ˛ huα,k |i∂k |uα,k idk ν1 = = π (1 + sign(J02 (y) − λ2 J02 (x))), 2 (2.133) where y = A0 (a0 − b0 ) and x = A0 b0 , and |uα,ki are the eigenfuctions of Hk0 . This result indicates, in contrast with the undriven model, that one can induce non-trivial topological phases even for λ > 1 (the trivial phase). This is a great example to demonstrate topological phases created by the a laser drive. 2.6 Calculation of Time Resolved Photoemission Photoemission spectroscopy (PES) refers to energy measurement of electrons emitted from materials, such as solids, by the photoelectric effect, in order to find the binding energies of electrons in a substance. It is an important time-resolved method, now in ultrafast region to study the effects of a pump laser in materials properties at the femtosecond time-scale. It is straightforward to understand the pump-prob process in time resolved pump-probe photoemission spectroscopy (TRPES). First the laser pulse pumps the surface of a material to induce an initial nonequilibrium response of the system and, second, the probe pulse would be used to study the excited system to generate the emitted photoelectrons. Here we follow the derivation of the time-resolved photoemission signal from Freericks et al [92, 56]. Time resolved angle resolved photoemission spectroscopy (trARPES) experiments measure the photoelectrons with a momentum ke collected at a solid angle dΩk̂ within an energy e interval dE. The system evolves from an initial condition {|φn i} in equilibrium, then the pump pulse will excite the system to an ensemble {|Ψn i}. After that, the probe pulse is turned on around time t = t0 . One can say {|Ψn i} = {U (t, −∞)|φn i} with U (t, −∞) is the time evolution of the system including the pump. The Hamiltonian has to be modified when 52 the probe pulse is coming as the following H = Hsolid (t) + Hf ree + Hc (t). (2.134) Note that the first term Hsolid (t) only contains the creation c† (t) and annihilation c(t) operators of electrons in the model without the effect of the probe pulse. Hf ree is for the free electron Hamiltonian, and can be written as Hf ree = X † [E(ke ) + W ]ak ake . e (2.135) k † Note that W is the work function, E(ke ) is the free electron kinetic energy and ak is the e creation operator and ake is the annihilation operator for a free electron with momentum ke . The third term Hc shows the coupling between the electrons in the solid and the free electrons out of the solid via an input photon with a wave number q. This term is of the form Hc = s(t − tp ) X † Mq (k, kq , q, t)ei~ωq t ak (t)ck (t)Aq . e (2.136) k We use a matrix element Mq (k, ke , q, t) to describe the absorption of a photon with energy ~ωq and ejection of an electron with momentum k on the surface of the material and momentum ke outside, and Aq is the annihilation operator of the photon. Since we have the time-varying pump laser pulse, the matrix element Mq (k, ke , q, t) is also time-dependent. The interaction between the free electron that leaves the solid and the electrons in the solid is assumed to be weak and can be ignored. The time evolution of the system is derived in terms p of the time-evolution operator Ũ (t, t0 ) and {|Ψn i} = {Ũ (t, t0 )Ψn (t0 )}. In the interaction picture, we can write the time-evolution operator with a probe pulse as the form Ũ (t, t0 ) = ´ −i/~ tt dt0 U (t,t0 )ˆ(H)c U (t0 ,t0 ) 0 U (t, t0 )Tt e . (2.137) Note that Hf ree commutes with the propagator U (t, t0 ) so we do not need to include it. Therefore, one can calculate the probability P (t) that an electron transfers from the ensemble 53 p {|Ψn i} to a free electron state for each time step t when the probe pulse is turned on. From Fermi’s golden rule, the probability is P (ke , t) = X p |hke |Ĥc |Ψn (t)i|2 . (2.138) n |ke i represents the free electron state with momentum ke . The system absorbs a photon of wave vector q and ejects an electron with a wave vector k = (kxy , kz ) inside the system and ke = (kexy , kez ) outside. In the process, since the momentum is conserved parallel to the surface, one can define kexy = kxy . With this conservation rule, one can derive the matrix element of the perturbation from Ĥc between the final state with momentum ke and the initial state |Ψn (t0 )i 1 |hke |Ĥc |Ψn (t)i| = | h ˆ ˆ t dkz Mq (k, ke , q, t) × 0 dt0 s(t0 )e−iωt × t0 † † 0 hke |a |Φm ihΦm |Ũ (t , t0 )ck Ũ (t0 , t0 )|Ψn (t0 )i|. (2.139) Here we made the assumption that once the electron is pumped to the excited state; inside the material, we have hke |a† |Ψm i = 1 |Ψm i with momentum k and a free electron with momentum ke will be generated through the matrix Mq . Here we define ω to be the energy of the excitation left in the system and it satisfies ~ω = ~ωq − (~ke )2 /(2me ) − W . We can expand the time-evolution matrix Ũ as ˆ i t 0 Ũ (t, t0 ) = U (t, t0 ) − dt U (t, t0 )Ĥc U (t0 , t0 ) ~ t0 (2.140) Since in the experiments, the amplitude of the probe pulse is small and the pump pulse is stronger than the probe pulse, one can apply the zeroth order perturbation in the time evolution operator Ũ (t, t0 ) = U (t, t0 ). Then one gets the probability ˆ ˆ 1 0 P (t, ke ) ≈ 2 dkz dkz I(t, w, êk ; kz , kz0 ) ~ (2.141) with I(t, w, êk ; kz , kz0 ) ˆ t = −i t0 dt00 ˆ t t0 dt0 Mq (kz , kez , kxy , t0 )Mq (kz0 , kez , kxy , t00 ) × 00 0 s(t00 )s(t0 )eiω(t −t ) G< (t0 , t00 ). k,k 0 54 (2.142) Note that G< 0 (t0 , t00 ) represents the two time lesser Green’s function in momentum space k,k with time-dependent creation and annihilation operators at k 0 and k. To be able to compare the theoretical results with experimental data, one needs further assumptions. One of the important assumptions is that the matrix element Mq will not change rapidly so we can take it as time-independent and it conserves the momentum parallel to the surface of the system. For these reasons, we take Mq to be constant in the calculations. With this assumption, we can directly calculate the time-resolved photoemission at site i and sum over all k. The above function of G< 0 (t0 , t00 ) works both for the momentum-diagonal Green’s 0 00 function G< k (t , t ) = k,k 0 00 G< k,k (t , t ) and for the local lesser Green’s function since the local Green’s function is calculated by summing over different momentum points as ˆ t Pii (ω, t) = −i 0 dt0 ˆ t 0 0 00 0 00 dt00 s(t0 )s(t00 )eiω(t −t ) G< ii (t , t ). (2.143) Note that, in the following chapters, we take s(t) as a Gaussian shape 1 s(t) = √ exp[−(t − tc )/σ 2 ] σ π Where tc is the central time of the probe pulse. 55 (2.144) CHAPTER 3 THE ZIG-ZAG SSH MODEL: FROM FLOQUET TO TRPES 3.1 Introduction In this chapter a detailed study of light driven PIPT in the Su, Schrieffer, and Heeger (SSH) model is presented [93]. This model was introduced in subsection 2.3.1. It was originally developed to describe transport in conducting polymers and more recently it has become a fundamental model for topological insulators and Majorana bound states [5, 81]. While the topologically trivial or nontrivial character of the SSH dimerized chain is controlled by the relative strength of the nearest neighbor couplings, it was only recently proposed that highfrequency laser light could turn trivial (non-topological) equilibrium bands into topological nonequilibrium Floquet bands [7]. We investigate these problems by simulating the real-time evolution of single-particle spectra in polyacetylene coupled to short laser pulses. We introduce a modified SSH model which allows us to tune topological states with polarized light. We compare the responses of initially trivial and nontrivial topological states, and we also calculate the steady state Floquet bands that are also used for comparisons. 3.2 Model and Methods To understand the nonequilibrium Floquet bands of the SSH model mentioned in Chapter 2, we consider a one dimensional tight-binding model that can be divided into two sites A and B with time-dependent hoping term τij only between the two closest sites as described in section 2.3.1. The time independent Hamiltonian is H= X † † (τ + δτ )cAi cBi + (τ − δτ )cAi+1 cBi + h.c. i In the following, a time-dependent effect will be included into the model. 56 (3.1) Figure 3.1: The Zig-Zag SSH model: (a) the trivial topological state δτ > 0 and (b) the nontrivial topological state δτ < 0. Note that the polarization is linear so the electric field is along one bond but perpendicular to the other. 3.2.1 Time-dependent Zig-Zag SSH Model We consider equal numbers of A and B sites with one electron per site, so that the electrons are at half filling. In addition we include a polarized time-dependent pump pulse incorporated into our model through the Peierls substitution, resulting in a time-dependent modulation of the hopping term. Note that, since a space uniform laser field is not able to change the original topological states [81], we consider a laser field parallel to (τ − δτ ) bond and perpendicular to (τ + δτ ) bond or vise versa and we show that this leads to topological phase transitions. This can be calculated using a Zig-Zag chain with light polarized along one of the band directions as Fig. (3.1) shows. We call this model the Zig-Zag SSH model. The Hamiltonian, using standard notation for creation and annihilation operators, is H(t) = X † † (τ + δτ )cAi cBi + (τ − δτ )(t)cAi+1 cBi + h.c. (3.2) i where the pump laser pulse is treated using the Peierls substitution as (see Section 2.4.1) # " ˆ ie Rj (τ − δτ )(t) = (τ − δτ ) exp A(r, t) · dr . (3.3) ~c Ri 57 Here A(r, t) is the time dependent vector potential in the Hamiltonian gauge, Ri is the position vector of the ith lattice site, and r is a real space vector. The δτ is the dimerization term and τ is taken as the unit of energy in this chapter (τ = 1). In the presence of a pump field, through the vector potential A(r, t), the Floquet bands form, and then disappear when the pump pulse amplitude goes to zero. In equilibrium, the energy bands have a gap of width equal to 4δτ that is symmetric about zero energy as shown in Fig. (2.6). For the time dependent case, the eigenfunctions change as a function of time and the quasi energy bands are modified by the presence of the external field. As described in ~ t) · ~σ Subsection 2.3.1, the Hamiltonian in momentum space can be written as H k (t) = d(k, with H k (t) is 2 × 2 matrix written in the basis of A and B sites, ~σ is the vector of Pauli matrices (σx , σy , σz ), and d~k (t) is a vector as dx (k, t) = (τ + δτ ) + (τ − δτ ) cos(ka − A(t)), dy (k, t) = (τ − δτ ) sin(ka − A(t)), dz (k, t) = 0, (3.4) Where a is unit cell length and, for simplicity, e = ~ = c = 1. Eigenenergies + and − exist as a pair, and this symmetry ensures the existence of two distinct topological phases (δτ > 0 and δτ < 0). We are able to calculate the time evolution of different momentum k states, from which we derive the nonequilibrium Green’s functions to obtain the physical properties of the system at time t (see Chapter 2). In the following calculation, we use ck (t) = (cAk (t), cBk (t)) and the equation of motion can be generally written as i~∂t ck (t) = H k (t)ck (t). (3.5) Since we have evolution equations, we can write down the time-evolution operator Uk (t, t0 )  ˆ i t dt̄H k (t̄) , = exp − ~ t0  58 (3.6) and, by using the Trotter formula, the time evolution in numerical calculations becomes Uk (t, t0 ) = Uk (t, t − ∆t)Uk (t, t − 2∆t)...Uk (t0 + ∆t, t0 ), (3.7) where ∆t is a small time step. For each moment time t, we have   i∆t H k (t − ∆t/2) . Uk (t, t − ∆t) = exp − ~ 3.2.2 (3.8) Floquet Hamiltonian for Zig-Zag SSH Model An effective Floquet Hamiltonian (see Section 2.5) can be used to derive Floquet spectra where for the Zig-Zag SSH model we find, HF = − X mΩ|m, αihm, α| mα + X [gm−n (k)|m, Aihn, B| + h.c.], (3.9) mn where gm−n (k) are the Fourier series expansion coefficients of the time-dependent hopping term g(k, t) = dx (k, t) − idy (k, t) from Eq. (3.4) that can be written as (see Section 2.5) ˆ 1 T gm−n (k) = dtei(m−n)Ωt g(k, t) T 0 π = (τ + δτ )δmn + (τ − δτ )ei[ka− 2 (m−n)] Jm−n (A0 ) (3.10) with a periodic vector potential A(t) = A0 cos(Ωt) and T = 2π/Ω. Note that Jm−n (A0 ) is the Bessel function of the first kind and we achieve convergence of the corresponding spectrum numerically when |m| ≥ 10 for all the cases discussed here. As an example, the 59 6 × 6 Floquet Hamiltonian matrix is   g0 0 g1 0 g2  Ω − λ    g∗ 0  0 g2∗ Ω − λ g1∗   0      0 g1 0−λ g0 0 g1  ,  M −λ=  ∗ ∗ ∗   g 0 g 0 − λ g 0 1 0 1       0 g 0 g −Ω − λ g 2 1 0     g2∗ 0 g1∗ 0 g0∗ −Ω − λ (3.11) where g0 = R + reika J0 (A0 ), g1 = −ireika J1 (A0 ) and g2 = −reika J2 (A0 ). Note that R = τ + δτ , r = τ − δτ and J0 , J1 are Bessel functions. The Bessel functions have the integral form, ˆ π 1 Jn (x) = ei(nτ −x sin(τ )) dτ. 2π −π 3.2.3 (3.12) Calculations from Non-equilibrium Green’s Functions By following the theory presented in Chapter 2 and Ref. [56, 55], the lesser Green’s function is defined as † 0 0 G< ij (t, t ) = ihcj (t )ci (t)i, (3.13) where we take a quantum statistical average of the time dependent creation and annihilation † operators. In the Heisenberg representation, ci (t) and ci (t) are creation and annihilation operators for a spinless fermion at site i. The angle brackets represent a trace over all quantum states in real or momentum space weighted by the equilibrium density matrix ~ t)|/[γ(τ 2 −δτ 2 )] initialized in the far past; and the density of states is 2(∂E(k)/∂k)−1 = 2|d(k, where γ is the normalization term to make sure that the integral of density of states is one. We can employ the Green’s function method to calculate the time-resolved photoemission spectroscopy signal (see Section 2.6), as a probe pulse weighted time Fourier transform of 60 the lesser Green’s function centered at time tp , ˆ PI (k, ω, tp ) = Im ˆ dt 0 0 dt0 s(t)s(t0 )e−iω(t−t ) G< I (k, t + tp , t + tp ). (3.14) with lesser Green function 1 < < G< I (k, t1 , t2 ) ≡ 2 [GAA (k, t1 , t2 ) + GBB (k, t1 , t2 )]. (3.15) The response should be calculated for each sublattice and is averaged over both sites A and B to compare with the experimental response. No extrapolation to large times is needed for this calculation since the probe pulse provides a natural cutoff. The probe pulse is assumed to be, s(t) = 2 2 1 √ e−t /α α π (3.16) with width α. The narrower the probe width, the better the time resolution and the worse the energy resolution. 3.3 3.3.1 Analysis of the SSH Floquet Hamiltonian Limit of High Frequency In high frequency regime E0 /Ω → 0 with E0 = A0 Ω, we only need to consider the g0 (m = n) term and all other terms can be set to zero. For any Floquet sideband centered at mΩ if we assume all gi6=0 (k) have a smaller effect than g0 (k), the energy band is approximately ± (m, k) ' mΩ ± g02 (k) where g02 (k) = R2 + r2 J02 (A0 ) + 2RrJ0 (A0 ) cos(ka) (3.17) with R = (τ + δτ ) and r = (τ − δτ ). We find the gap ∆0 (k = π) = 2g02 (k = π) is dependent on the electric field amplitude of the drive laser. An interesting feature is that A0 can be chosen so that flat sidebands occur, when J0 (A0 ) = 0, and reverse sidebands when 61 J0 (A0 ) < 0. We also note that gap opening depends on the amplitude of the light and also on whether the system is in the trivial or non-trivial phase; i.e. δτ > 0 or δτ < 0. When δτ > 0 we have R > rJ0 (A0 ) so the behavior of the gap is dominated by R. On the other hand, when rJ0 (A0 ) > R the effect of the light through rJ0 (A0 ) would strongly influence the behavior of the gap opening. We show below that the different topological phases can lead to completely different spectral responses when higher order terms gi6=0 (k) are included. Since we know HF from above, we can calculate the Berry phase or winding number from the corresponding eigenvectors uα,k with α the band index, ˛ ν = huα,k |i∂k |uα,k i. (3.18) Ignoring terms gi6=0 , it is possible to calculate the exact winding number as [81] ν0 = 3.3.2 π {1 + sgn[J02 (A0 ) − (R/r)2 ]} 2 (3.19) Simple Model for the Case of One Overlap To understand physics of the light driven SSH model for one overlap, deriving an exact solution is a long standing challenge. Some high frequency limit results have been calculated [81] recently but not an exact solution. While to solve a large Hamiltonian from Floquet theory is not a reasonable approach, we found a small size Floquet Hamiltonian which only considers two interbands is quite enough to describe the dispersion relation of light driven SSH model for one overlap between bands. We would start with a 4 × 4 Hamiltonian for two interband interactions as   g0 0 g1  Ω − λ   ∗  g∗  Ω − λ g 0   0 1 H =  ,  0 g1 −λ g0      ∗ ∗ g1 0 g0 −λ 62 (3.20) For simplicity, one can first make the diagonals symmetric,   Ω −λ g 0 g 0 1  2   Ω ∗ ∗  g 0  g1   0 2 −λ H =  .   0 Ω −λ g − g   1 0 2   −Ω 0 g0∗ g1∗ − λ 2 (3.21) The corresponding eigenvalues can then be derived directly as λ = ±(β ± α1/2 )1/2 (3.22) Ω β = |g0 |2 + |g1 |2 + | |2 2 (3.23) α = (Ω|g0 |)2 + (g0 g1∗ + g0∗ g1 )2 . (3.24) with and Note that, |g0 |2 = R2 + (rJ0 (E))2 + 2RrJ0 (E) cos(ka), (3.25) |g1 |2 = (rJ1 (E))2 , (3.26) (g0 g1∗ + g0∗ g1 )2 = −4(rRJ1 (E))2 (cos(ka))2 . (3.27) |g0 |2 is a term that varies with cos(ka) and |g1 |2 is a constant with respect to ka. Interestingly, (g0 g1∗ + g0∗ g1 )2 changes with (cos(ka))2 and induces a second order modulation of the band (a local minimum/maximum that will curve the band). 3.3.3 Computational Results and Comparison to Analytic Results One can directly compute eigenenergies and eigenvectors by solving the Floquet Hamiltonian. The Floquet Hamiltonian used in the calculations has size 2(2m + 1). m indicates a Floquet frequency nΩ in the range of [−mΩ, mΩ] and m, n ∈ Z. We choose m = 10 to ensure that 63 7 6 6 5 5 4 4 Energy / τ Energy / τ 7 3 2 1 3 2 1 0 0 −1 −1 −2 -1.0 -0.5 0.0 k(a/π) 0.5 −2 -1.0 1.0 -0.5 0.0 k(a/π) 0.5 1.0 Figure 3.2: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 5. Dashed lines: the exact solution for the 4 × 4 Floquet Hamiltonian. Solid lines: Numerical results for the 42 × 42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). the bands near m = 0 converge with respect to variations in m. In the following, comparison between computational results for m = 10 and exact result from the 4 × 4 Hamiltonian will be presented. In Fig. (3.2), the dispersion relation for a no overlap case is presented; for laser amplitude A0 = 1 and frequency Ω = 5. The exact solution of the 4 × 4 Floquet Hamiltonian (dashed lines) and results for the 42×42 Floquet Hamiltonian with m = 10 (solid lines) are compared. One finds that the exact 4 × 4 solution is quite accurate in the high frequency regime. Only for the nontrivial case the outer bands of the 4 × 4 solution are inaccurate indicating that interband affects are important in this case. Fig. (3.3) shows the scheme of band interactions for g0 , g1 and g2 , where n is the Floquet band index. g0 is the intraband process; while g1 and g2 are interband process. g1 is related to the one photon interaction and g2 is related to the two photons interaction. In Fig. (3.4), a one overlap case is presented with laser amplitude A0 = 1 and frequency Ω = 3. The exact solution for the 4 × 4 Floquet Hamiltonian (dashed lines) and results for the 42 × 42 Floquet Hamiltonian with m = 10 (solid lines) are compared for the trivial condition δτ = 0.5 (R > r) and the nontrivial condition δτ = −0.5 (R < r). The inner two bands for the exact 4 × 4 solution are very close to the computational results for the 64 Figure 3.3: Scheme of band interactions for g0 , g1 and g2 . n is the Floquet band index. g0 is the intraband process; while g1 and g2 are interband process. one overlap situation. For the trivial case, however in the left figure in Fig. (3.4) the gap at crossing points is not captured by the 4 × 4 model. By making gn≥2 = 0 in the m = 10 calculation we found that the gap closing effect is caused by higher order mixing through the terms g0 and g1 terms which cannot be derived from the 4 × 4 Hamiltonian as Fig. (3.3) shows. For the nontrivial case the outer bands of the 4 × 4 solution have some shifting indicating effects due to bands above and below. In Fig. (3.5), a two overlap case is presented with laser amplitude A0 = 1 and frequency Ω = 1.5. Neither the trivial nor nontrivial cases are fully described by the 4 × 4 Hamiltonian in this case though the nontrivial case deviates the most, the need for higher order terms to describe the gap of the nontrivial case is again evident. 65 5 4 4 3 3 Energy / τ Energy / τ 5 2 1 0 −1 −2 -1.0 2 1 0 −1 -0.5 0.0 k(a/π) 0.5 −2 -1.0 1.0 -0.5 0.0 k(a/π) 0.5 1.0 4 4 3 3 2 2 Energy / τ Energy / τ Figure 3.4: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 3. Dashed lines: the exact solution for the 4 × 4 Floquet Hamiltonian. Solid lines: Numerical results for the 42 × 42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). 1 0 −1 −2 -1.0 1 0 −1 -0.5 0.0 k(a/π) 0.5 −2 -1.0 1.0 -0.5 0.0 k(a/π) 0.5 1.0 Figure 3.5: Dispersion relation for the laser driven SSH model with laser amplitude A0 = 1 and frequency Ω = 1.5. Dashed lines: the exact solution for the 4 × 4 Floquet Hamiltonian. Solid lines: Numerical results for the 42 × 42 Floquet Hamiltonian with m = 10. Left is for the trivial condition δτ = 0.5 (R > r) and Right is for the nontrivial condition δτ = −0.5 (R < r). 3.3.4 Effects of Varying Amplitude and Frequency In Fig. (3.6), we show the dispersion relations of the quasi-static Floquet Hamiltonian with m = 10 and drive frequency Ω = 2.8. Left: trivial condition δτ = 0.5 and Right: the nontrivial condition δτ = −0.5. (a)(b) is for amplitude A0 = 2, (c)(d) is for amplitude A0 = 4 and (e)(f) is for amplitude A0 = 6. Clearly a variety of different dispersion relations can be achieved through use of polarized light and a well designed lattice, Gap opening and closing can also be changed by different frequencies, amplitudes, and polarization of light, 66 as well as lattice structure. Energy / τ Energy / τ 3 2 1 0 −1 −2 −3 Energy / τ δτ = 0. 5 3 2 1 0 −1 −2 −3 3 2 (e) 1 0 −1 −2 −3 -1.0 -0.5 δτ = (a) (b) (c) (d) − 0. 5 (f) 0.0 k(a/π) 0.5 1.0 -1.0 -0.5 0.0 k(a/π) 0.5 1.0 Figure 3.6: Dispersion relations of quasi-static Floquet Hamiltonian with drive frequency Ω = 2.8. Left: trivial condition δτ = 0.5 and Right: the nontrivial condition δτ = −0.5. (a)(b) is for amplitude A0 = 2, (c)(d) is for amplitude A0 = 4 and (e)(f) is for amplitude A0 = 6. Fig. (3.7) shows the quasienergy properties calculated from the quasi-static Floquet matrix of the original bands at k = π vs laser amplitude A0 for δτ = ±0.5. The left figure is for the trivial condition δτ = 0.5 and the right figure is the nontrivial condition δτ = −0.5 for different drive frequency Ω. At high frequency, they oscillate as the Bessel function of the first kind J0 (A0 ) but the bands have different trends in the trivial and nontrivial cases. For the trivial case, the energy bands remain open and fluctuate with the laser amplitude while the gap increases with increasing frequency. For large frequency the bands become similar to 67 2.0 1.5 1.5 1.0 1.0 0.5 0.5 0.0 Ω = 2. 0 Ω = 2. 4 Ω = 2. 8 Ω = 3. 2 Ω = 3. 6 Ω = 4. 0 −0.5 −1.0 −1.5 −2.0 Energy / τ Energy / τ 2.0 0 2 4 6 A0 8 10 12 0.0 Ω = 2. 0 Ω = 2. 4 Ω = 2. 8 Ω = 3. 2 Ω = 3. 6 Ω = 4. 0 −0.5 −1.0 −1.5 −2.0 14 0 2 4 6 A0 8 10 12 14 Figure 3.7: Quasienergy spectrum properties calculated from the Floquet matrix with m = 10 for the original bands at k = π vs laser amplitude A0 for δτ = ±0.5. Left: are for the trivial condition δτ = 0.5 and Right: are for the nontrivial condition δτ = −0.5 for different driving frequencies Ω. each other. For the nontrivial case, the gap closes at laser amplitude A0 ≈ 1.8, and is always smaller than the trivial case when frequency Ω > 2.8. When amplitude increases all bands become similar (A0 > 6.2). Note that the chosen frequencies are in the range 2.0 ≤ Ω < 4.0 where only one band overlap occurs. 1.0 0.6 0.4 0.2 0.0 Ω = 1. 0 Ω = 2. 0 Ω = 3. 0 Ω = 4. 0 Ω = 5. 0 0.8 Energy / τ 0.8 Energy / τ 1.0 Ω = 1. 0 Ω = 2. 0 Ω = 3. 0 Ω = 4. 0 Ω = 5. 0 0.6 0.4 0.2 0 2 4 A0 6 8 10 0.0 0 2 4 A0 6 8 10 Figure 3.8: Bandwidth calculated from quasi-static Floquet band structures vs laser amplitude A0 for δτ = ±0.5. Left: are for the trivial condition δτ = 0.5 and Right: are for the nontrivial condition δτ = −0.5 for different driving frequency Ω. In Fig. (3.8), we show the bandwidth calculated from quasi-static Floquet band structures vs laser amplitude A0 for δτ = ±0.5. The left are for the trivial condition δτ = 0.5 and the right are for the nontrivial condition δτ = −0.5 for different driving frequency Ω. Note that, the N overlap regime can be decided by ( N4τ+1 , 4τ N ), e.g. one overlap regime is in the frequency range (2τ, 4τ ). In the high frequency domain Ω > 4 (no overlap) the laser 68 amplitude dominates; so the bandwidths oscillate with the Bessel function of the first kind |J0 (A0 )|. However the nontrivial bandwidths are larger than the trivial bandwidths since the variation is decided by [τ − (−δτ )]J0 (A0 ) > [τ − (δτ )]J0 (A0 ) with δτ = 0.5. In the low frequency domain, the influencies from other second order coupling through Ji (A0 ) with i = 1, 2, 3, ... also plays a role in bandwidth variations so the behaviors change. Especially for the case Ω = 3 we have monotonic decreases in bandwidth as amplitude increases. It is then clear that both the amplitude and frequency of the drive laser can be used to control the electronic structure. 3.3.5 Topological Phase Diagram We found the topological phase diagram by evaluating Eq. (3.18) using the Floquet matrix eigenvectors uα,k . Here we discuss the topological phases for the trivial case δτ > 0 and the nontrivial case δτ < 0. In Fig. (3.9), the topological phase diagrams for different frequencies and amplitudes are presented. Note that different topological regions can have different winding numbers ν (0, π, 2π, ...) but the same Berry phase (0 or π). It is clear that the topological phase diagrams are totally different between the trivial and the nontrivial cases (See also [81]). The original phase for the case δτ = 0.25 is trivial and changes to a nontrivial phase at low frequency when Ω < 2.0. On the other hand, The original phase for the case δτ = −0.25 is nontrivial and changes to a trivial phase in the high frequency domain when Ω > 2.0. This is also true for large values of δτ as shown in Fig. (3.9) (c)-(f). (c)(d) shows the cases for δτ = 0.5 on the left and δτ = −0.5 on the right. The high frequency results from Eq. (3.19) describes this regime well. The four red points indicate the frequncies and amplitudes we choose to discuss in Fig. (3.16) of the next subsection. 69 Ω Ω Ω (a) δτ = 0. 25 (b) δτ = − 0. 25 (c) δτ = 0. 50 (d) δτ = − 0. 50 (e) δτ = 0. 75 (f) δτ = − 0. 75 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 A0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 A0 Figure 3.9: Topological phase diagram for different frequencies and amplitudes of the drive pulses from effective Floquet Hamiltonian. (a)(b) shows the cases for δτ = 0.25 on the left and δτ = −0.25 on the right. (c)(d) shows the cases for δτ = 0.5 on the left and δτ = −0.5 on the right. (e)(f) shows the cases for δτ = 0.75 on the left and δτ = −0.75 on the right. The four red points indicate the frequencies and amplitudes we discuss in Fig. (3.16) of the next subsection. 3.4 Pump-probe Results In this section we describe fully time dependent pump-probe calculations using nonequilibrium Green’s functions (NEGF). The vector potential of the pump pulse is of the form A(t) = A0 cos(Ωt) exp(−t2 /2σ 2 ) where A0 is amplitude, Ω is frequency and σ is pulse width. 70 When σ → ∞, this should lead to behavior like that observed in Floquet Theory. However, there is a key difference. The NEGF calculation gives the occupancy of the nonequilibrium states; which is not the same as the states themselves. 3.4.1 Pseudospin Content We would like to use NEGE to make contact with the pseudospin representation of 2 × 2 Hamiltonians that reflect their orbital content and fully determine the Berry phase in and out of equilibrium and in particular for low driving frequencies. When a description in terms of a simple effective Floquet Hamiltonian is not available, one has to find an analogue of a pseudospin analysis in terms of nonequilibrium Green functions. The pseudospin content of Green functions from our numerical simulations is extracted by expanding the Green function matrices in orbital representation using Pauli matrices [83], 1 < [GAB (k, t1 , t2 ) + G< BA (k, t1 , t2 )], 2 i < < G< y (k, t1 , t2 ) ≡ − [GAB (k, t1 , t2 ) − GBA (k, t1 , t2 )], 2 1 < < [GAA (k, t1 , t2 ) − G< Gz (k, t1 , t2 ) ≡ BB (k, t1 , t2 )]. 2 G< x (k, t1 , t2 ) ≡ (3.28) The respective pseudospin content Px,y,z (k, ω, tp ) is obtained by computing the analogue < of the PES response for G< x,y,z (k, t1 , t2 ) instead of GI (k, t1 , t2 ) using Eq. (3.14). In or- der to obtain the Berry phase with pseudospin vector S(k, tp ) = (Sx , Sy , Sz ) shown in Fig.(3.13), we integrate Gx,y,z in a frequency window to cover a single band as Sx,y,z (k, tp ) = ´ ωu ω dωPx,y,z (k, ω, tp ). We note that the pseudospin behavior is dependent on the choice of l Floquet sideband for which it is analyzed. Floquet bands may have different pseudospin content, and the pseudospins in the upper and lower bands within a manifold always point in opposite directions. The proof of pseudospin content can be found in the following. 71 First, consider a time-independent steady-state 2 × 2 Hamiltonian as   0 dx − idy  ~ = d~ · ~σ =  H(d)   dx + idy 0 (3.29) ~ = d = (d2 + d2 )1/2 and corresponding eigenenergies, with |d| x y ~ H(d)|±i = ±d|±i. (3.30) and time-dependent eigenvectors,   e∓idt  1  |±, ti = e∓idt |±i = √  , 2 ±eiφ (3.31) with de±iφ = dx ± idy . Now we show how to derive di from the lesser Green’s functions. First put σj with j = {x, y} and e∓idt on both sides of Eq. (3.30), we will have ~ σj H(d)|±, ti = σj (±d)|±, ti. (3.32) Using Eq. (3.29) this equation can be written as (dj 1 + idl σz )|±, ti = σj (±d)|±, ti, (3.33) where l = {x, y} is not equal to j. We take the inner product of the above equation with another eigenvectors at a different time t0 to find, h±, t0 |(dj + idl σz )|±, ti = ±dh±, t0 |σj |±, ti. (3.34) h±, t0 |σz |±, ti = 0, (3.35) Since we know the above equation can be rewritten as dj h±, t0 |±, ti = ±dh±, t0 |σj |±, ti. 72 (3.36) In this formula the eigenvectors | ↑i and | ↓i of σz can be considered as electron occupancy components in sublattice A and B as shown, in the main text Eq. (3.5) so we can represent the right-hand side in Eq. (3.36) as the pseudospin content of the lesser Green’s 0 functions G< k,j (t, t ; φ(0)) and the left-hand side is an identity for the lesser Green’s function 0 G< k,I (t, t ; φ(0)) with a given initial state φ(0) for momentum k. Without loss of generality Eq. (3.36) can then be written as 0 G< k,j (t, t ; ±) = ± dj (k) < Gk,I (t, t0 ; ±) d(k) (3.37) where j = x, y, z. Note that any state can be a linear combination of eigenvectors so we can apply Eq. (3.37) to any given state. One can also show that for the direction of pseudospin content of the lesser Green’s function is opposite for the two different eigenvectors |±i. To have an exact calculation of pseudospin content from PES response we can put Eq. (3.37) into Eq. (3.28) and assume that the laser probe pulse is of rectangular form in the range of [− σ2 , σ2 ] so we have ˆ σ Pj (k, ω, tp ) ' ±Im × 2 −σ 2 ˆ σ dt 2 −σ 2 0 dt0 e−iω(t−t ) dj (k) < G (t + tp , t0 + tp ). d(k) k,I (3.38) For the initial state |+i we then derive the pesudospin content from PES response as Pj (k, ω, tp ) ' dj (k) sin2 [(d(k) − ω)σ/2] . d(k) (d(k) − ω)2 (3.39) Since the integral term is positive, the only term that would affect the sign of the pseudospin content from PES is the sign of dj ; and the initial state |±i. The intensity of pseudospin content is then proportional to the ratio dj (k)/d(k) and the density of states. Hence, when the Hamiltonian is time-dependent, we can still calculate psedudospin content from PES and the results can be directly connected with an effective time-independent Hamiltonian. Below it is shown that the average outcome shows consistency with steady state calculations from the effective Floquet Hamiltonian. 73 3.4.2 Evolution of States and Their Occupancy In this subsection we present trPES results for the dynamic generation of Floquet states. The initial state in these calculation is a filled lower band and an empty upper band. Energy / τ 3 t= − 200. 0 (a) t= − 100. 0 (b) t = 0. 0 t = 100. 0 (c) (d) t = 200. 0 (e) 2 1.000 1 0.100 0.010 0 0.001 -1 -2 0.000 -3 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 k(a/π) k(a/π) k(a/π) k(a/π) k(a/π) Figure 3.10: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = 0 under the pulse frequency Ω = 1, width σ = 50 and amplitude A0 = 1. Dashed lines show the corresponding quasi-static Floquet band structures at the same frequency. In Fig. (3.10), we study time evolution of time-resolved PES of the electronic SSH model for times t = −200, −100, 0, 100, 200 with δτ = 0 for a pulse frequency Ω = 3, width σ = 50 and amplitude A0 = 1. Dashed lines show the corresponding quasi-static Floquet band structures. At time t = −200 the pulse starts with a small amplitude so the Floquet states are not generated with significant occupancy. At time t = −100 the occupation of the nearest Floquet bands begins to be observed. At time t = 0 occupation of many Floquet bands can be seen and, interestingly, their is strong intensity of occupation at the intersection points of different bands. At times t = 100 and t = 200, we discover that the occupation at the intersection points between the Floquet bands and the original upper band remains while the Floquet states disappear. The anomalous occupation of the upper band is momentum selected and can be considered as a resonant excitation enhanced by Floquet effects. In Fig. (3.11), we study time evolution of time-resolved PES of the electronic SSH model for time t = −200, −100, 0, 100, 200 with δτ = 0.5 for pulse frequency Ω = 3, width σ = 50 74 Energy / τ 3 t= − 200. 0 (a) t= − 100. 0 (b) t = 0. 0 t = 100. 0 (c) (d) t = 200. 0 (e) 2 1.000 1 0.100 0.010 0 0.001 -1 -2 0.000 -3 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 k(a/π) k(a/π) k(a/π) k(a/π) k(a/π) Figure 3.11: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = 0.5 under the pulse frequency Ω = 3, width σ = 50 and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. As for the case δτ = 0, the generation of states is consistent with Floqeut theory with the additional feature of resonant excitation at band crossing points. Energy / τ 3 t= − 200. 0 (a) t= − 100. 0 (b) t = 0. 0 t = 100. 0 (c) (d) t = 200. 0 (e) 2 1.000 1 0.100 0.010 0 0.001 -1 -2 0.000 -3 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 -1.0-0.5 0.0 0.5 1.0 k(a/π) k(a/π) k(a/π) k(a/π) k(a/π) Figure 3.12: Time evolution of time-resolved PES of the electronic SSH model is presented for time t = −200, −100, 0, 100, 200 with δτ = −0.5 under the pulse frequency Ω = 3 and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. In Fig. (3.12), we study time evolution of time-resolved PES of the electronic SSH model 75 for δτ = −0.5 with the pulse frequency Ω = 3, width σ = 50 and amplitude A0 = 2. Dashed lines show the corresponding quasi-static Floquet band structures. The time variation for this nontrivial phase is totally different to the two previous cases. One can also see second order mixing between bands near k = 1 and k = −1 when t = 0 since the electron occupancy near k = 1 and k = −1 moves up with the second order mixing. Interestingly, at time t = 200, the occupancy for the second order mixing stays in the upper band. There are two reasons for these differences: (i) The Floquet bands change with amplitude A0 and the occupation should be with the pulse amplitude at that time. (ii) Excitation only occurs between states of the same phase. As in the cases δτ = 0; and δτ = 0.5; resonant excitation at intersection points occurs. 3.4.3 Dynamics of Phase Evolution Fig. (3.13) shows a Berry phase mapping using the renormalized radius r0 = 1 + (1 + k)/2 for δτ = 0.5 on the left and δτ = −0.5 on the right with pulse frequency Ω = 3, amplitude A0 = 2 and pulse width σ = 50. Note that the renormalized radius r0 is for illustration purposes in order to separate the paths that go around the origin more than once. (a)(b) shows the trivial and nontrivial Berry phases before the laser pulse comes in at time t = −200 with integral frequency window ωl = 0.0 and ωu = 3.0. (c)(d) is at time t = 0 with integral frequency window ωl = −1.5 and ωu = 0.0. For the trivial case, since the bands still mix, the Berry connection moves up and down twice about the original line Sy = 0; leading to a Berry phase of zero. (e) shows the trivial band with integral frequency window ωl = 1.5 and ωu = 3.0 at time t = 0; in this case there is no phase flip. Note that it has reverse direction from (c) since the band has the opposite phase ψ(k 0 ) + π. Some bands can have the same Berry phase but different winding numbers ν, as we seen by comparing (b) and (f). In the latter case the Berry connection goes around the origin twice so the winding number ν = 4π. To further illustrate the novel topological behaviours due to a laser pulse, in Fig. (3.14) 76 2 δτ =0.5 (a) (b) δτ = −0.5 Sy 1 0 −1 −2 2 (c) (d) (e) (f) Sy 1 0 −1 −2 2 Sy 1 0 −1 −2 −2 −1 0 Sx 1 2 −2 −1 0 Sx 1 2 Figure 3.13: Berry phase mapping using the radius r0 = 1 + (1 + k)/2 with δτ = 0.5 on the left and δτ = −0.5 on the right under the pulse frequency Ω = 3, amplitude A0 = 2 and pulse width σ = 50. (a)(b) is at time t = −200 with integral frequency window ωl = 0.0 and ωu = 3.0, (c)(d) is at time t = 0 with integral frequency window ωl = −1.5 and ωu = 0.0 and (e)(f) is at time t = 0 with integral frequency window ωl = 1.5 and ωu = 3.0. we show the time evolution of pseudospin phase ψ deriving from the pseudospin context P eiψ = Px + iPy that are presented for times t = −200, −100, 0, 100, 200 with δτ = 0.5 on the top (a)-(e) and δτ = −0.5 on the bottom (f)-(j) with pulse frequency Ω = 3, amplitude A0 = 1 and pulse width σ = 50. We can see in figures (a)(f) the upper band points in reverse direction of P, which means the phase of the upper band is ψ(k 0 ) + π if lower band is at phase ψ(k 0 ) with momentum k 0 . In (b), however, the occupied Floquet states cross with the original band leading to a change of phase. In contrast, Floquet bands in (g) have the same 77 t= − 100. 0 t= − 50. 0 t= − 25. 0 t = 0. 0 Energy / τ − 200. 0 Energy / τ t= 3 (a) 2 1 0 -1 -2 -3 3 (f) (g) (h) (i) (j) ν=π 2 ν = 2π 1 0 -1 ν=0 -2 -3 -1.0-0.50.0 0.5 1.0 -1.0-0.50.0 0.5 1.0 -1.0-0.50.0 0.5 1.0 -1.0-0.50.0 0.5 1.0 -1.0-0.50.0 0.5 1.0 ν=0 (b) (c) (d) (e) ν=0 3.0 ν=0 1.5 k(a/π) k(a/π) k(a/π) k(a/π) 0.0 −1.5 −3.0 k(a/π) Figure 3.14: Time evolution of pseudospin phase ψ deriving from pseudospin context P eiψ = Px + iPy are presented for time t = −200, −100, −50, −25, 0 with δτ = 0.5 on the top (a)-(e) and δτ = −0.5 on the bottom (f)-(j) under the pulse frequency Ω = 3, amplitude A0 = 2 and pulse width σ = 50. Dashed lines show the corresponding quasi-static Floquet band structures. phase ψ(k 0 )+π as the phase of the upper band. We find the phase of these Floquet sidebands are decided by the topological phase of the original band where δτ > 0 (trivial case) and δτ < 0 (non-trivial case). From (b) we can find the phases are different at the crosspoints so there is no gap opening at time t = 0, while we will have degeneracy at intersection and would induce a state mixing or pauli exclusion between two chiral fermions so to open a gap at time t = 0 as shows in (h). As time t increases the gap in the nontrivial case δτ < 0 becomes larger and causes the bands to move apart. When the bands touch other bands they mix again with a reverse phase band and form a new band with different Berry phase as one can find in (i)(j). On the other hand, the trivial case will not open a gap so we can only see the flattening effects from the pulse amplitude. 3.4.4 Effects of Varying Pump and Probe Pulsewidth and Amplitude Fig. (3.15) shows several time-resolved PES results with δτ = 0 for a pulse with frequency Ω = 3, amplitude A0 = 1, and for a range of different pump pulse widths; for full occupancy of the lower band as the initial state. When the pulse width is small (a) σ = 1 and (b) 78 (b) σ =2.0 1.000 σ =5.0 σ =10.0 Energy / τ Energy / τ σ =1.0 3 (c) 2 1 0 -1 -2 -3 Energy / τ 3 (a) 2 1 0 -1 -2 -3 3 (e) (f) 2 1 0 -1 -2 -3 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 (d) 0.100 0.010 0.001 σ =20.0 σ =40.0 k(a/π) 0.000 k(a/π) Figure 3.15: Time-resolved PES of the electronic SSH model with δτ = 0 at time t = 0 with a pulse of frequency Ω = 3, amplitude E = 1 and width (a) σ = 1, (b) σ = 2, (c) σ = 5, (d) σ = 10, (e) σ = 20, (f) σ = 40. Dashed lines show the corresponding quasi-static Floquet band structures. 79 σ = 2, we find that most of the electrons are only excited to the upper band and spread out to the nearby energy states. As the pulse width becomes larger (c) σ = 5 and (d) σ = 10, we can see a much clear generation of the Floquet bands with a weaker dispersion to the nearby energy states. Note that for the case (c) the pulse width is about 10 times the period of the probe. When the pulse width is larger (e) σ = 20 and (f) σ = 40, we find a clear occupancy at the quasi-static Floquet band structures. Floquet physics is then accessible for pulse widths of order 10 times the period of the probe. Fig. (3.16) shows several TRPES results for full occupancy of the lower band with δτ = ±0.5. (a)(b) gives the energy spectrum with laser pulse of frequency Ω = 1.6, amplitude A0 = 0.5 and width is σ = 50. One can find similar responses both for δτ = ±0.5 but with weaker excitation for δτ = 0.5 of low frequency and small amplitude. (c)(d) shows Floquet band formulation for low frequency but high amplitude region with laser frequency Ω = 1.6, amplitude A0 = 2.0 and width σ = 50. We can see that the occupancy of the Floquet bands increases with the intensity of the laser pulse, and the effect of a high amplitude pulse on the Floquet spectra is clearly observed. (e)(f) shows a high frequency region, where the different topological behaviours for δτ = 0.5 and δτ = −0.5 are evident. A gap opening only occurs for the case δτ = −0.5 at the intersection of the Floquet band with origin band. The occupancy of the electrons in nontrivial case is also split at the crossing points. Note that these different responses to the laser can also be found in (c)(d) but it is not as obvious. 3.5 Conclusion We find that short optical pulses can lead to local spectral and novel pseudospin textures in a one-dimensional topological insulator given by the Zig-Zag SSH model. Pump-probe photoemission spectroscopy can probe these states by measuring sizeable energy gaps and Floquet band formation on femtosecond time scales. Analysing band structures and pseudospin textures, we identify new states with optically induced nontrivial changes of sublattice mixing that leads to novel topological phenomenon. This study reveals the possibility to dis- 80 cover new topological phases driven by optical pulses by turning the lattice structure and polarization of the light. Energy / τ (d) Energy / τ 3 (c) 2 1 0 -1 -2 -3 3 (e) 2 1 0 -1 -2 -3 (f) Energy / τ δτ = −0.5 Energy / τ δτ =0.5 3 (a) 2 1 0 -1 -2 -3 3 (g) (h) 2 1 0 -1 -2 -3 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 (b) 1.000 0.100 0.010 0.001 k(a/π) 0.000 k(a/π) Figure 3.16: Time-resolved PES with δτ = 0.5 on the left and δτ = −0.5 on the right at time t = 0. (a)(b) with pulse of frequency Ω = 1.6, amplitude A0 = 0.5 and width is σ = 50, (c)(d) with pulse of frequency Ω = 1.6, amplitude A0 = 2.0 and width is σ = 50, (e)(f) with pulse of frequency Ω = 3, amplitude A0 = 0.5 and width σ = 50, (g)(h) with pulse of frequency Ω = 3, amplitude A0 = 2 and width σ = 50. Dashed lines show the corresponding quasi-static Floquet band structures. 81 CHAPTER 4 KROTOV OPTIMAL CONTROL THEORY 4.1 Introduction to Optimal Control In many physical areas, quantum control via specially tailored laser pulses is a long-standing goal. In the most recent decade, this goal has been achieved, as sophisticated pulse shaping experiments can now coherently control some quantum states in quantum computing, quantum dots and phase transitions in materials. For example, laser pulses may be applied to create and break a particular bond in a molecule, to control charge transfer within molecules, or to optimize high harmonic generation and phase transitions in topology [39, 38, 53, 54, 55, 56, 57, 58, 59, 35, 60, 32, 33, 51, 3, 35, 36, 61, 60]. Optimal control theory (OCT) can help to theoretically design the laser pulse to transfer an initial state to a given final state. Especially Krotov’s optimal control theory or other similar optimal control theories are widely used in recent years. This chapter provides an introduction to Krotov’s optimal control theory. It indicates how the control parameters and equation of motion define such an optimal control pulse [62]. Krotov’s method has several significant advantages over the traditional gradient methods. First of all, monotonic increase in fidelity or decrease of cost/goal function with iteration number. Secondly, no requirement for an exhaustive line search, and thirdly update in macrosteps in time and at each iteration. Multiple versions of Krotov’s optimization method have been used recently in many areas to deal with the Markovian and the non-Markovian optimal control problem of a quantum Brownian motion model with an exact stochastic equation of motion (master equation) in quantum computing. It has also been applied to condensed matter physics in recent years enabled by progress in spectroscopy and laser control [63, 64, 65, 66, 67, 68, 69, 70, 71, 62, 72, 73]. In this Chapter, We first explain the foundation and schemes to solve the control equation of motion and exhibit how to include additional constraints into the goal/cost functions. 82 Second, We introduce the tricks of the Krotov optimal control method and the corresponding backpropagator. Third, We present the monotonic convergence property and related proof. Fourth, We discuss the first order and second order cases for the backpropagator format, and algorithms for implementation. Last, We will show a few examples to demonstrate the power of the optimal control method. This is a review Chapter providing background for the developments described in Chapter 5. 4.2 Preliminary Preparation of the Krotov Method In this section the fundamentals of Krotov optimal control theory are sketched. For Krotov optimal control theory, one needs to know the equation of motion of a given system, and decide the cost/goal function to indicate the goal to be achieved. The formula of a goal function depends on a target function of the system and a set of control parameters c. Normally one can use square error as the target: (f (T, x, c) − y)2 or correlation metrics such as: y T f (T, x, c) with y is a goal function, x are other parameters and T is terminal time. To be more precise, one can consider an equation of motion of the form ∂x = f [t, x(t), c(t)], ∂t (4.1) and assume we want to minimize the goal/cost function ˆ T I[v] = f 0 (t, x(t), c(t))dt + F [x(T )] −→ min. (4.2) 0 Here x(t) is the time-evolution function or the trajectory of the system, c(t) is a timedependent control parameter. Note that v = (x(t), c(t)) ∈ D with D is the set of permissible process v = (x(t), c(t)) satisfying Eq. (4.2). The functionals f 0 (t, x(t), c(t)) and F [x(T )] are defined for all t, x(t), c(t). Both need to be twice differentiable with respect to c and x. The initial value x(0) = x0 is a given value and x(T ) is the value of x(t) at final time T. Note that c(t) can be chosen within a closed domain U . The general functional F [x(T )] depends only 83 on the final time step of x(t) and f 0 (t, x(t), c(t)) depends on the intermediate time steps of x(t) and c(t) so one can say the goal function I is a general functional based on the terminal and intermediate time steps of x(t). For a quantum system with a multi-dimensional vector space or multi-argument processes and control parameters, we will have more than one equation of motion, ∂xi = f i [t, ~x(t), ~c(t)], ∂t i = 1, 2, ..., n. (4.3) And the minimization problem becomes ˆ T I[v] = f 0 [t, ~x(t), ~c(t)] + F [~x(T )] (4.4) 0 Note that, ~c(t) and ~x(t) now are vectors, e.g. ~x = (x1 , x2 , ..., xn ). 4.3 4.3.1 The Tricks of Krotov’s Method Decomposition of Goal Function To implement the Krotov method, a real and twice differentiable function ϕ[t, x(t)] is considered. The idea is to use the function ϕ[t, x(t)] to help update the control parameters. Note that this function would derive from a reverse time evolution start with the final time t = T to the initial time t = 0. The control parameters are updated by comparing the difference between the original forward evolution function and backpropagator ϕ[t, x(t)]. We define the constructions: R[t, x, c] = ∂ϕ ∂ϕ f [t, x, c] − f 0 [t, x, c] + , ∂x ∂t G[T, x] = F [T, x] + ϕ[T, x], (4.5) (4.6) ˆ T L[v, ϕ] = G[T, x(T )] − R[t, x(t), c(t)]dt − ϕ[0, x(0)]. 0 84 (4.7) To be more specific, L[v, ϕ] = I[v] for any function ϕ[t, x(t)]. The following is the proof: ˆ T R[t, x(t), c(t)] − ϕ[0, x(0)] L[v, ϕ] = G[T, x(T )] − 0 ˆ T ∂ϕ ∂ϕ = G[T, x(T )] − [ f [t, x(t), c(t)] − f 0 [t, x(t), c(t)] + ]dt ∂x ∂t 0 −ϕ[0, x(0)] ˆ T ∂ϕ dx ∂ϕ = G[T, x(T )] − [ + − f 0 [t, x(t), c(t)]]dt ∂x dt ∂t 0 −ϕ[0, x(0)] ˆ T dϕ dt − ϕ[0, x(0)] = F [T, x(T )] + ϕ[T, x(T )] − 0 dt ˆ T + f 0 [t, x(t), c(t)]dt 0 ˆ T f 0 [t, x(t), c(t)]dt = F [T, x(T )] + 0 = I[v]. (4.8) Therefore minimizing I[v] can be achieved by minimizing L[v, ϕ], and this also minimizes G[x(T )] and maximizes R[t, x(t), c(t)]. For a multi-dimensional quantum system or for multi-argument processes, the equations for R and G will be written as R[t, ~x(t), ~c(t)] = ∂ϕ ∂ϕ ~ f [t, ~x(t), ~c(t)] − f 0 [t, ~x(t), ~c(t)] + ∂~x ∂t (4.9) and G[T, ~x(T )] = F [T, ~x(T )] + ϕ[T, ~x(T )]. (4.10) For later use, it is convenient to define the function Φ = ∂ϕ ∂x , and the functional R[t, x(t), c(t)] = H[t, x(t), c(t), Φ(t)] + ∂ϕ , ∂t (4.11) where H[t, x(t), c(t), Φ(t)] = Φ(t)f [t, x(t), c(t)] − f 0 [t, x(t), c(t)]. 85 (4.12) Note that the parameters in H denoted by Φ(t) emphasize that x and Φ should be treated as independent variables with respect to H. 4.3.2 Iterative Algorithm The main purpose of the Krotov method is to find optimal control sequence ck+1 (t) to improve the goal/cost function. In other words, the Krotov method hopes that I[v] is monotonically decreasing with respect to ck (t) as k increases. That is I[vk ] ≥ I[vk+1 ] (4.13) at every iteration. Since ϕ[t, x(t)] is not restricted, we can freely choose the form of ϕ[t, x(t)]. However, if we can construct the function ϕ[t, x(t)] to maximize L[vk , ϕ] at each k then we can randomly choose the next control sequence ck+1 (t) and it will increase the value of L[v, ϕ]. We therefore derive a smaller value of the goal function by the chosen ϕ. Specifically, we suppose that we already found the function ϕ for a problem, then the complete processes will be as follows: (i.) Take an initial control sequence c0 (t), and calculate the trajectory x0 (t) from Eq. (4.1). (ii.) Choose the functional ϕ[t, x(t)] to make L[v0 , ϕ] a maximum with the control c0 (t) and trajectory x0 (t). This requirement is equivalent to the following two conditions: R[t, x0 (t), c0 (t)] = min R[t, x, c0 (t)], (4.14) G[T, x0 (T )] = max G[T, x]. (4.15) x x The conditions imply that the functionals R and G are calculated using the new ϕ[t, x]. As a result the current control sequences c0 (t) will be the worst of all possible c(t) in minimizing the goal functional L[v, Φ] = I[v]. Any new c(t) changing from c0 (t) with a new trajectory x will now improve the minimization of the goal function I[v]. 86 (iii.) Finding a new control sequences c1 that maximizes the functional R. The corresponding conditions are c1 [t, x] = max R[t, x, c] c = max H[t, x, c, Φ], c (4.16) where H is mentioned in Eq. (4.12). Note that the control sequence c1 [t, x(t)] depends on the trajectory function x(t). (iv.) With the new control sequence c1 [t, x] the new trajectory x1 (t) can be derived using the equation of motion Eq. (4.1). (v.) It is now guaranteed that the goal function in Eq. (4.2) has been minimized monotonically, which can be written as I[v1 ] ≤ I[v0 ]. The new control sequences and the trajectory become a starting point of the next iteration and (i.)-(iv.) can be repeated to further decrease the goal function. 4.3.3 Monotonic Convergence of Krotov Method First we outline the proof that the new function I[v1 ] is indeed smaller than the previous I[v0 ]. It is straightforward to show that ∆I = I[t, x0 (t), c0 (t)] − I[t, x1 (t), c1 (t)] = L[t, x0 (t), c0 (t), Φ] − L[t, x1 (t), c1 (t), Φ] ˆ T = R[t, x1 (t), c1 (t)] − R[t, x0 (t), c0 (t)]dt + G[T, x0 (T )] − G[T, x1 (T )] 0 = ∆ 1 + ∆2 + ∆ 3 (4.17) where ∆1 = G[T, x0 (T )] − G[T, x1 (T )] ˆ T ∆2 = R[t, x1 (t), c1 (t)] − R[t, x1 (t), c0 (t)]dt 0 ˆ T ∆3 = R[t, x1 (t), c0 (t)] − R[t, x0 (t), c0 (t)]dt. 0 87 (4.18) (4.19) (4.20) Using the conditions in Eq. (4.14) and Eq. (4.15) one can prove that ∆1 ≥ 0 and ∆3 ≥ 0, and Eq. (4.16) also guarantees ∆2 ≥ 0. Therefore the new goal/cost functional I will be smaller than the previous one and the monotonically convergence is proved. 4.4 Construction of ϕ To carry out the above iteration method, the most important and hardest task is finding a ϕ that satisfies the conditions in Eq. (4.14) and Eq. (4.15) which require the absolute minimum of the functional R and maximum of the functional G with the old control sequences c0 and the old trajectory x0 . In this section, we will show how to construct ϕ to first order in x and to second order in x to treat linear and non-linear problems respectively. 4.4.1 First Order in x Consider that the equations of motion of the system are linear and can be written as ∂x = f [t, x, c] = a[t, c]x + b ∂t (4.21) and the functions f 0 [t, x, c] and F [T, x] are concave with respect to x, which means ∂ 2 f 0 [t, x, c] ≤ 0, ∂x∂x ∂ 2 F [T, x] ≤ 0. ∂x∂x (4.22) In this case, we just need to consider ϕ to first order in x since the second derivative is guaranteed. To be more specific, The first order that implies the functional ϕ needs to satisfy Eq. (4.14) and Eq. (4.15). The function ϕ needs to ensure that the first derivative of the functions R and G with respect to x are equal to zero. Therefore, we can choose the 88 function ϕ[t, x] = Φ(t)x which satisfies the following conditions: ∂R[t, x0 , c0 ] ∂H[t, x0 , c0 , Φ] ∂Φ[t, x0 ] = + ∂x ∂x ∂t = 0, (4.23) ∂F [T, x0 ] ∂ϕ[T, x0 ] ∂G[T, x0 ] = + ∂x(T ) ∂x ∂x ∂F [T, x0 ] = + Φ[T, x0 ] ∂x = 0. (4.24) Therefore, Eq. (4.23) is the equation of motion for the function Φ : ∂H[t, x0 , c0 , Φ] ∂Φ[t, x0 ] =− ∂t ∂x (4.25) with boundary conditions Eq. (4.24) ∂F [T, x0 ] ∂x (4.26) ∂x ∂H[t, x0 , c0 , Φ] = ∂t ∂Φ (4.27) Φ[T, x0 ] = − From Eq. (4.1) and Eq. (4.12) To satisfy the above requirements, the possible choice of ϕ is ϕ = Φ(t)x. In the multiargument process, the similar choice of the functional ϕ would be ϕ[t, ~x(t)] = Φi (t)xi (t). Using the formula of Eq. (4.12), the conditions can be rewritten into the form ∂Φ ∂f 0 [t, x0 , c0 ] = −J T (t)Φ(t) + , ∂t ∂x (4.28) where J = ∂f [t, x0 , c0 ] . ∂x Note that, in the multi-argument process, J T (t) is a transpose matrix. 89 (4.29) 4.4.2 Second Order in x If the equations of motion of the system are not linear, one needs to consider a different form of ϕ. Since functional ϕ needs to satisfy Eq. (4.14) and Eq. (4.15), the simplest choice of functional ϕ is of the form 1 ϕ[t, x] = Φ(t)x + ∆xΣ(t)∆x, 2 (4.30) where the ∆(x) ≡ x − x0 and both the function ϕ[t, x] and the matrix Σ(t) should be found. Here Σ(t) is the Hessian or matrix of the second derivatives of the function ϕ(t, x). The first necessary conditions for inequalities of Eq. (4.14) and Eq. (4.15) are equivalent to Eq. (4.25) and Eq. (4.26), and the second necessary conditions for inequalities of Eq. (4.14) and Eq. (4.15) yield the following differential inequalities: d2 R ≥ 0, d2 G ≤ 0, ∂ 2 R[t, x0 , c0 ] ∆x, ∂x∂x ∂ 2 G[T, x0 ] d2 G = ∆x ∆x. ∂x∂x d2 R = ∆x (4.31) (4.32) Because the functional ϕ can be chosen arbitrarily, in the multi-argument process, one can require that the matrix Σ(t) is a diagonal matrix and satisfies the conditions ∂ 2 R[t, x0 (t), c0 (t)] = 0, i 6= j, i, j = 1, 2, ..., n, ∂xi ∂xj ∂ 2 R[t, x0 (t), c0 (t)] = σii (t), σii (t) ≥ 0, i = 1, 2, ..., n, ∂xi ∂xi (4.33) and ∂ 2 G[T, x0 (T )] = 0, i 6= j, i, j = 1, 2, ..., n, ∂xi ∂xj ∂ 2 G[T, x0 (T )] = σii (T ), σii (T ) ≤ 0, i = 1, 2, ..., n. ∂xi ∂xi (4.34) One therefore determines the equation of motion of Σ(t) with chosen boundary conditions σii (t) and σii (T ) using one of the above linear differential equations. 90 4.4.3 Algorithm In the previous section, Krotov’s optimal control method was introduced. Here I will summarize the algorithm: (1) Choose an initial form of the control function c0 (t). (2) Use Eq. (4.1) and initial condition x(0) = x0 to find the trajectory x0 (t). (3) Find the functional Φ(t) using Eq. (4.23) and Eq. (4.24) or Eq. (4.25) and Eq. (4.26). (4) In the non-linear case, use Eq. (4.33) and Eq. (4.34) to find the matrix Σ(t). (5) With the functional ϕ, the control c1 (t) is found according to Eq. (4.16). (6) Derive the new trajectory x1 (t) using the control sequence c1 (t) by Eq. (4.1). (7) Repeat processes (2) to (6) until the desired optimal value is achieved. 4.5 4.5.1 Examples Example 1: First Order in x Consider a linear problem with φ chosen in the form given in subsection (4.4.1) for the following optimal control problem: ẋ(t) = i(1 + c(t))x(t), x(0) = 1; c(t) is real and one wants to minimize the cost function ˆ T 1 (c(t0 ) − c0 )2 dt0 −→ min. I = Re[(1 + x(T ))] + b 2 0 (4.35) (4.36) where b > 0. We choose the parameters b = 5, T = 2 and substitute the linear form of ϕ = Φ[t, x(t)]x to derive R and G : R = Re[Φ(t)[i(1 + c(t))x(t)] + ∂ 1 (Φ(t)x(t))] − b(c − c0 )2 (t), ∂t 2 G = Re[(1 − x(T )(−1)) + Φ(T )x(T )]. 91 (4.37) (4.38) Using ∂ ∂R = Re[Φ(t)[i(1 + c(t))] + Φ(t)] = 0, ∂x ∂t ∂G = Re[1 + Φ(T )] = 0, ∂x (4.39) (4.40) one can derive the equation of motion of Φ in reverse time t → −t: Φ(T ) = −1. Φ̇(t) = i(1 + c(t))Φ(t), (4.41) Performing the algorithm described in subsection (4.4.3), we require ∂ 2R ≤ 0. (∂c)2 ∂R = 0, ∂c (4.42) We therefore obtain the result c1 (t) = c0 (t) + Re[Φ(t)ix(t)]/b. (4.43) Results are shown in Fig. 4.1, where we have used the Runge Kutta method with the segment of integration partitioned into 200 intervals, and the fidelity is defined as F = Re[x(T )(−1)] and error = 1 − F as shown in cost function. 4.5.2 Example 2: Second Order in x Let us consider the approach from subsection (4.4.2) for the following optimal control problem. The functions x(t) and c(t) are constructed by ẋ = c, |c| ≤ 1, x(0) = 0; (4.44) and one wants to minimize the cost function ˆ T I = 1 (c2 − x2 )dt + bx2 (T ) −→ min, 2 0 where b > 0. 92 (4.45) 1.00 0.66 0.95 0.64 0.62 0.60 0.85 c(t) Fidelity 0.90 0.80 0.56 0.75 0.54 0.70 0.65 0 0.58 0.52 5 10 15 Iteration 20 0.50 0.0 0.5 1.0 t 1.5 2.0 Figure 4.1: Left: fidelity versus iteration times. Right: optimal control sequence with respect to time t. Now we choose the parameters b = 20, T = 4 substitute Eq. (4.30) into Eq. (4.5) and Eq. (4.6) to derive R and G of the form 1 R = Φ̇(t)x(t) + Σ̇(t)(∆x)2 + 2Φ(t)c(t) 2 +2Σ∆x(t)c(t) − c2 (t) + x2 (t), 1 1 G = Φ(T )x(T ) + Σ(T )(∆x(T ))2 + bx2 (T ). 2 2 (4.46) (4.47) Since ∂ 2R = Σ̇(t) + 1 ≥ 0, (∂x)2 ∂ 2G = Σ(T ) + b ≤ 0, (∂x)2 (4.48) we first choose that Σ̇(t) = 0, Σ(T ) = −b − 4. 93 (4.49) We therefore bring it into ∂R = Φ̇(t) + Σ̇(t)∆x + 2Σc(t) + 2x(t) = 0, ∂x ∂G = Φ(T ) + Σ(T )∆x(T ) + bx(T ) = 0. ∂x (4.50) (4.51) Performing the algorithm described in subsection (4.4.3), we require ∂R = 0, ∂c ∂ 2R ≤ 0. (∂c)2 (4.52) We obtain the result shown in Fig. 4.2, where we have used the Runge Kutta method with the integration partitioned into 200 intervals. For comparison, the known solution of the problem is shown below    ±t, t ≤ τ1 ,    x(t) = ±kcos(t − T /2), τ1 ≤ t ≤ τ2 ,      ±T ∓ t, τ2 ≤ t, (4.53) where T is the final time and τ1 , τ2 and k are chosen according to smoothness conditions: ẋ = ±1 for t = τ1 , ẋ = ∓1 for t = τ2 , ±t = ±kcos(t − T /2) at t = τ1 and ±kcos(t − T /2) = ±T ∓ t at t = τ2 . Note that the result of the Krotov optimal method in Fig. 4.2 is equal to the known solution. We will extend the Krotov method to investigate a quantum gate in the next example. 4.5.3 Example 3: A Two Level Quantum System The time dependent Schrödinger equation for the evolution equation (forward-propagator) |ψ(t)i of a quantum system including a time-dependent control term µε(t) where ε(t) is the control function i~ ∂ |ψ(t)i = (H + µε(t))|ψ(t)i. ∂t Suppose that one want to minimize the cost function ˆ T I = 1 − Re[hψG |ψ(T )i] + λ (ε(t0 ) − ε0 )2 dt0 −→ min. 0 94 (4.54) (4.55) 160 2.0 140 120 1.5 80 x(t) Fidelity 100 60 1.0 40 0.5 20 0 −20 0.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 1 2 3 4 5 6 7 8 9 10 Iteration t Figure 4.2: Left: Cost function versus iteration times. Right: Optimal evolution of x with respect to time t. where λ > 0 and ε0 is an initial energy that can be considered as the restriction for the optimal control sequences; furthermore, ε0 can also be time-dependent. |ψG i is a target goal for the forward-propagator |ψ(T )i. Now we choose the parameters λ = 1, T = 1 and substitute ϕ = hβ(t)|ψ(t)i into Eq. (4.5) and Eq. (4.6) to derive R and G :  ∂hβ(t) R = Re hβ(t)|(H + µε(t))|ψ(t)i + |ψ(t)i − λ(ε(t) − ε0 )2 , ∂t i h G = Re hψG |ψ(T )i − hβ(t)|ψ(T )i T0 .  (4.56) (4.57) Using Eq. (4.25) and Eq. (4.26), one can derive the equation of motion of β(t) (similar process to Example 1): i~ ∂ hβ(t)| = hβ(t)|(H + µε(t)), ∂t hβ(T )| = hψG |. (4.58) To optimize the control sequence, we require ∂R = 0, ∂ε ∂ 2R ≤ 0. (∂ε)2 95 (4.59) Therefore the optimal control sequence is of the form ε(t) = ε0 + 1 Re[hβ(t)|µ|ψ(t)i]. 2λ (4.60) Using the algorithm in subsection (4.4.3) with the above conditions where hβ(t)| depends on the old ε(t), and |ψ(t)i is built from the new ε(t). Note that for a better performance, one can substitute ε0 with the old ε(t) to derive new ε(t). The system under consideration is a quantum qubit or a two-level quantum dot. Under appropriate conditions, the Hamiltonian of the qubit reads H(t) = −ε(t)σz /2 − Ωσx /2 where H = Ωσx /2 and µ = σz /2 and Ω is a bias voltage. If we consider a target    1  |ψG i = | ↑i =   , 0 (4.61) (4.62) we obtain the result shown in Fig. 4.3 and Fig. 4.4, where we use the Euler method with the segment of integration partitioned into 320 intervals (dt = 0.01), Ω = 1, and fidelity is defined as Re[hψG |ψ(T )i] with initial state    0  |ψ(t = 0)i = | ↓i =   . 1 (4.63) We can achieve fidelity 10−4 or lower with learning rate λ = 10. Note that this is a quantum optimal control problem and is the fundamental in quantum computation. Chapter 5 will use similar optimal techniques to a charge density wave (CDW) tight-binding model. 96 1.0 1.0 0.5 0.6 ε(t) Fidelity 0.8 0.0 0.4 −0.5 0.2 0.0 −1.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0 20 40 60 80 100 120 140 Iteration t Figure 4.3: Left: Cost function versus iteration times. Right: Optimal control pulse of ε respect to time t. 1.0 Probability 0.8 0.6 up down 0.4 0.2 0.0 0.0 0.5 1.0 1.5 Time 2.0 2.5 Figure 4.4: Optimal probability evolution of ψ(t) respect to time t. 97 3.0 CHAPTER 5 QUANTUM OPTIMAL CONTROL OF A TRANSIENT CDW STATE 5.1 Introduction As outlined in Chapter 1, 2, 3, High fidelity characterization of transient excited many-body electron distributions in the ultrafast time domain is now available through a variety of pump-probe experiments. Amongst the rich variety of non-equilibrium responses observed, photo-induced phase transitions (PIPT) are particularly interesting for fundamental and applied reasons [39, 38, 53, 54, 55, 56, 57, 58, 59, 35, 60, 32, 33, 51, 3, 35, 36, 61, 60]. For example, by adjusting the laser pulse properties it is possible to tune the non-equilibrium PIPT response from the adiabatic to the non-adiabatic limits. Optimal laser pulse shaping methods have not been applied to PIPT yet and have the potential to control non-equilibrium response in order to isolate selected physical phenomena, and in order to tune response for selected applications such as high speed electronics or optics. We introduce a method, based on Krotov optimal control theory, to direct photo-induced phase transitions (PIPT) by combining non-equilibrium models with Quantum optimal control theory (QOCT). To illustrate the approach, we consider PIPT in a simple model [56] for a transient metal-insulator state in a charge density wave system. In experiments, a long-range chargedensity-wave (CDW) is formed in a variety of layered chalcogenide materials, in oxides, in two dimensional materials, and in many other systems. The investigated layered CDW materials have disclosed a new nonequilibrium pattern where the long-range CDW is preserved while the local electronic excitation spectrum becomes gapless for a transient period of time, as has been eluciated using the simple model considered here [56]. In the experimental systems the mechanism for gap-closing and population inversion remains an open question, though in some cases there has been significant theoretical progress. QCOT is a powerful tool based on calculating the optimal pulse shape by minimizing a physical cost function or maximizing a desired physical objective, and it has been developed 98 within a variety of variational frameworks to obtain control sequences [63, 64, 65, 66, 67, 68, 69, 70, 71, 62, 72, 73]. Quantum optimal control methods based on the classical gradient optimization methods provide an alternative to iterative methods based on the Krotov approach. Krotov methods have been applied to the fields of Quantum computing and control of charge transfer processes [72, 74, 75, 62, 76, 73]. The Krotov approach has several appealing advantages over the gradient methods in following ways: First, a monotonic increase toward the objective with iteration number. Second, no requirement for a line search and faster convergence to a given target. Third, at each iteration it guarantees macrosteps of the time interval and it can reach the global maximum. As shown here, the Krotov approach can easily be combined with CDW models for high efficiency calculations. To illustrate the methods, we consider the hot-electron model following the simulation study of Ref. [56, 55]. They found ultrafast laser pulses can quickly heat electrons forming a hot quasithermal gas that equilibrates with phonons on much longer time scales compared to the electron relaxation time. We introduce QOCT in this scenario to optimize shaped-laser pulses to study the short-time transient phenomenon and to discover new photo-induced phase transitions, and to control the outcome of pump-probe PIPT experiments to achieve strong population inversion. 5.2 Charge Density Wave Model To model the CDW system, we follow the work by Freericks et al. [56, 55]. We consider a lattice that can be divided into two sublattices A and B with a hoping term τij only between the two sublattices. For on site energy, sublattice A is chosen equal to U and sublattice B is choose equal to zero. We consider equal numbers of A and B sites with one electron per site, so that the electrons are at half filling and form an insulator. As outlined in Chapter 2, we include a spatially uniform time-dependent pump pulse incorporated into our model through the Peierls substitution, resulting in a time-dependent modulation of the hopping 99 term. The Hamiltonian, using standard notation for creation and annihilation operators, is H(t) = − X † τij (t)ci cj + (U − µ) ij X † X † ci ci − µ ci ci , i∈A (5.1) i∈B where the chemical potential is µ = U/2 for half filling and the time-dependent hopping in the presence of the laser pulse using Peierls substitution is # " ˆ ie Rj A(r, t) · d~r . τij (t) = τij exp ~c Ri (5.2) Here A(r, t) is the time dependent vector potential in the Hamiltonian gauge, Ri is the position vector of the ith lattice site, and r is a real space vector. ~ t), the electronic band In the presence of a pump field, through the vector potential A(r, gap can be tuned [56], even though τij (t) only affects the phase but not the amplitude of the hopping term of the Hamiltonian. In equilibrium, the density of states presents a gap of width equal to U that is symmetric about zero energy. For the time dependent case, the eigenfunctions change as a function of time and the nonequilibrium density of states is modified by the presence of the external field. For this reason the non equilibrium energy states are always influenced by the interaction of the system with the pulse field, and hence it is possible to control the local charge density. Translation symmetry is broken when U is nonzero, as the conversion from real space to momentum space is more complicated than in a system with one atom per unit cell. The momentum points k and k + Q are coupled, where Q = (π, π, π, ...) due to the presence of the CDW order. The transformation from reciprocal space to real space becomes † ci = X † † (e−ik·Ri ck + e−i(k+Q)·Ri ck+Q ), (5.3) k The electronic band in momentum space for the U = 0 case is time dependent and can be written as εk (t) = − X τij exp[−i(~k − ij 100 e A(t)) · (RiA − RjB )] ~c (5.4) where uniform vector potential is A(t) = A(t)(1, 1, ..., 1) and is directed along the diagonal of the hypercubic lattice in d-dimensions, and we can rewrite as εk (t) = εk Ac (t) + ε̄k As (t), (5.5) Ac (t) = cos(eaA(t)/~c), (5.6) As (t) = sin(eaA(t)/~c), (5.7) d X τ∗ √ cos(akl ), d l=1 (5.8) d X τ∗ √ sin(akl ). ε̄k = − d l=1 (5.9) with and εk = − In this work we take the limit that dimension d goes to infinity as in the model of Ref. [56]. In our work, the two functions Ac (t) and As (t) are used as control parameters for optimal control. To simplify the notation, we set ~ = e = c = 1 and choose U = t∗ , where t∗ is the renormalized hopping for a hypercubic lattice in infinite dimensions. We also set the initial temperature of the system equal to zero before the field is turned on. We can now derive equations of motion for the creation operators in momentum space as dck (t) = [εk (t) + U2 − µ]ck (t) + U2 ck+Q (t) dt dck+Q (t) i~ = U2 ck (t) + [−εk (t) + U2 − µ]ck+Q (t). dt i~ (5.10) (5.11) We are able to calculate the time evolution of different momentum k states, from which we derive the nonequilibrium Green’s functions to obtain the physical properties of the system at time t. In the following calculations, we use ck (t) = (ck (t), ck+Q (t)) and define a P superposition operator a(t) = k ck (t) in extended Hilbert space. Then the time evolution following Eq. (5.10) and (5.11) can be generally written as i~∂t ck (t) = H k (t)ck (t) or 101 P i~∂t a(t) = H(t)~a(t) with H(t) = k H k (t). and   U U εk (t) + 2 − µ  2 H k (t) =  . U U −µ −ε (t) + k 2 2 (5.12) Since we have evolution equations, we can write down the time-evolution operator as Uk (t, t0 )  ˆ i t¯ dtH k (t̄) , = exp − ~ t0  (5.13) and, by using the Trotter formula, the time evolution in numerical calculation becomes Uk (t, t0 ) = Uk (t, t − ∆t)Uk (t, t − 2∆t)...Uk (t0 + ∆t, t0 ), (5.14) where ∆t is a small time step. For each moment time t, we have   i∆t Uk (t, t − ∆t) = exp − H k (t − ∆t/2) . ~ 5.2.1 (5.15) Equations for the Nonequilibrium Solution In the paper, by following the theory of [56, 55], the retard and lesser Green’s function are defined as † 0 0 0 GR ij (t, t ) = −iθ(t − t )h{ci (t), cj (t )}+ i, † 0 0 G< ij (t, t ) = ihcj (t )ci (t)i, (5.16) (5.17) where we take a quantum statistical average of time dependent creation and annihilation operators. θ(t − t0 ) is the unit step function. In the Heisenberg representation, ci (t) and † ci (t) are creation and annihilation operators for a spinless fermion at lattice site i. The angle brakets represent a trace over all quantum states in real or momentum space weighted by the equilibrium density matrix as is initialized in the far past. We can employ the Green’s function to calculate the election concentration on the A and B sublattices, the DOS, the current, the total energy and the time-resolved PES signal. 102 We calculate the time-resolved PES response function as a probe pulse weighted time Fourier transform of the lesser Green’s function centered at time tp (see Section 2.6), ˆ ˆ 0 0 P (ω, tp ) = −i dt dt0 s(t)s(t0 )e−iω(t−t ) G< ii (t + tp , t + tp ). (5.18) The response should be calculated for each sublattice and is then averaged over both sublattices to compare with the experimental response. No extrapolation to large times is needed for this calculation since the probe pulse provides a natural cutoff. The probe pulse is assumed to be a gaussian function 2 2 1 s(t) = √ e−t /σ σ π (5.19) with width σ. The narrower the probe width, the better the time resolution and the worse energy resolution. We also calculate the total energy, the current and the occupancy in the upper/lower bands of the instantaneous band structure varies with time. Only equal time expectation values are needed considered for these quantities, and They are found in a straightforward manner from the equal time lesser Green’s function. We write the calculations here in terms of the equal time expectation values. Time-dependent number density for different † momentum k and k + Q are defined as nk1 ,k2 (t) = hck (t)ck2 (t)i with (k1 , k2 ) is (k, k + Q). 1 The current in the case satisfies X j(t) =  ∇ε(k; t) nk,k (t) − nk+Q,k+Q (t) , (5.20) nk+Q,k (t) − nk,k+Q (t) . nk,k (t) + nk+Q,k+Q (t) (5.21) k:ε(k)<0 The order parameter is of the form X Ω(t) = k:ε(k)<0 and the total energy becomes E(t) = X   ε(k; t) nk,k (t) − nk+Q,k+Q (t) k:ε(k)<0   U + n (t) + nk,k+Q (t) . 2 k+Q,k 103 (5.22) The occupancy of the upper and lower instantaneous bands become n+ (t) = X h α2 (k; t)nk,k (t) + β 2 (k; t)nk+Q,k+Q (t)) k:ε(k)<0  + α(k; t)β(k; t) nk+Q,k (t) + nk,k+Q (t) (5.23) and n− (t) = h X β 2 (k; t)nk,k (t) + α2 (k; t)nk+Q,k+Q (t) k:ε(k)<0  − α(k; t)β(k; t) nk+Q,k (t) + nk,k+Q (t) (5.24) with the two time-dependent coefficients are given by U/2 α(k; t) = q p 2[ε2 (k; t) + U 2 /4 − ε(k; t) ε2 (k; t) + U 2 /4] (5.25) and p ε2 (k; t) + U 2 /4 β(k; t) = q . p 2 2 2 2 2[ε (k; t) + U /4 − ε(k; t) ε (k; t) + U /4] −ε(k; t) + (5.26) 5.3 Quantum Optimal Control Method In this section, we develop an optimal control method that maximizes the overlap between the state of the system and a target (goal) state at time tf as implemented by QOCT. If we denote an arbitrary target operator (state) as g, then the goal is to control the system in such a way that the final operator is as close as possible to g. Fidelity is one way to quantify performance of the optimal pulse. Fidelity is defined as the real part of the inner product between the desired target operator g and the actual operator a(tf ) at time tf (see e.g. [69, 72, 73, 71, 62]), F = Re[hg † a(tf )i]N , 104 (5.27) with N denotes normalization, so that the maximum fidelity is F = 1. In many optimal theories and machine learning, we restrict control parameters to a finite range is necessary to prevent divergence or discontinuity, and, in physical problems, the power of the pump pulse should be considered in a finite range as well. We are therefore looking for a solution of the optimal control method that can achieve high fidelity while keeping the total power transferred to the system from the control parameters finite. To achieve the goal, we require that the power close to the power of the pulse that we use as previous iteration. This is expressed by [72] K= ˆ t f 0 dt0 λ(t0 )[ X 0 l 0 2 (Al+1 j (t ) − Aj (t )) ], (5.28) j=c,s with l, l + 1 represent the l, l + 1 iteration times of the quantum optimal control algorithm. We set A0c (t) → 1 and A0s (t) → 0 for the 0th iteration because of the initial condition of A(t) = 0. λ(t0 ) is a positive function chosen empirically that can be adjusted during the optimal control procedure. It influences the rate of change of the control parameters at each iteration step. We choose λ(t0 ) as a constant proportional to the number of energy points used in the simulation. Since we know (Ac (t))2 + (As (t))2 = 1, we add an additional constraint through a Lagrangian multiplier R to ensure that this equation is satisfied for all time steps t. Defining β(t) as another adjustable and positive function, we have R= ˆ t f dt0 β(t0 )((Ac (t0 ))2 + (As (t0 ))2 − 1). (5.29) 0 When everything is put together, an objective function is introduced of the form J = F − K − R, (5.30) where F , K and R are defined in Eqs. (5.27)-(5.29). Our goal is to maximize the objective function given the equation of motion for the operator a(t). In practice the algorithm consists of the following steps: First, choose an initial control sequence Ac (t) and As (t), second, given the initial condition the operator a at very beginning, 105 apply the equations of motion to find the forward propagated operator al (t). Third, calculate an auxiliary backward evolved operator bl (t) analogous to a(t) with the condition b† (tf ) = g † . Fourth, propagate al+1 (t) forward in time and update the control sequence iteratively with the relation l Al+1 c (t) = (Ac (t) + δAc (t))/α(t), (5.31)   ∂H(t) l+1 1 l † Re hb (t) a (t)i , δAc (t) = 2λ(t) ∂Ac (t) (5.32) and 2 where α(t) = 1+β(t)/λ(t) is an positive and adjustable term, which ensures that (Al+1 c (t)) + 2 (Al+1 s (t)) = 1 is satisfied for every time t. Specifically α(t) is an arbitrary function greater or equal to one and it can be arbitrarily assigned at every iteration and every time t since the contribution of the Lagrangian multiplier in the objective function J is 0. The equation for Al+1 s (t) has the same form, in Eq. (5.32). Steps (3) and (4) with Als , Al+1 and can be obtained by replacing Alc , Al+1 s c are repeated until a desired fidelity F is reached or until a given number of iterations has been accomplished. Note that here we derive δAc,s (t) base on the initial conditions, while in other cases, one can derive state-independent results by using a trace instead of an inner product. We terminate our program when 1-F is smaller than 10−4 or when the fidelity |F l+1 − F l |/|1 − F l | is smaller than 10−3 . 5.3.1 Proof of Optimality Using the cost function mentioned in the text ˆ t X f † J =Re[h~g ~a(tf )i] − dtλ(t)[ (Aj (t0 ) − A0j (t0 ))2 ] 0 − ˆ t f 0 j=c,s dtβ(t)(A2c (t) + A2s (t) − 1), (5.33) we will now prove that the iteration algorithm for the cost function J exhibits monotonic convergence. The simple proof for two-level quantum system can be found in Ref. [69]. First, 106 following Krotov’s method (see Chapter 4), we group the terms in J in the following way J¯ = G(tf ) + ˆ t f dtR(t). (5.34) 0 Here G(tf ) depends only on the terminal time tf and is defined as  G(tf ) ≡ Re  tf h~g †~a(tf )i − h~b† (t)~a(t)i 0 , (5.35) where ~b(t) is an arbitrary continuously differentiable propagator which can be considered as a Lagrange multiplier function constraining the system to obey the equation of motion. The second term R(t) is related to the time integral part and is of the form " # ~b† (t) ∂ R(t) ≡ Re h~b† (t)H(t)~a(t) + ~a(t)i ∂t X (Aj (t0 ) − A0j (t0 ))2 ] −λ(t)[ j=c,s −β(t)(A2c (t) + A2s (t) − 1). (5.36) In our case we have two control parameters in the Hamiltonian H(t) = H 0 + H c Ac (t) + H s As (t). To maximize J¯ one can simply maximize G and R independently. Note that if R is maximized at each time t the integral of R will be maximized as well. Second, to prove that J¯ converges at every iteration, it is straightforward to show that J¯(l+1) − J¯(l) = ∆1 + ∆2c + ∆2s + ∆3 ≥ 0 (5.37) where ∆1 ≡ G(~a(l+1) (tf )) − G(~a(l) (tf )) = Re[h(~g † − ~b† (tf ))∆~a(tf )i], 107 (5.38) ∆2c ≡ ˆ t f 0 (l+1) dt[R(t, ~a(l+1) (t), Ac (l) (t), As (t)) (l) (l) − R(t, ~a(l+1) (t), Ac (t), As (t))] ˆ t  f † (l+1) = Re dth~b (t)(Ĥ c ∆Ac (t))~a (t)i 0 ˆ t f (l+1) dt2λ(t)(Ac (t) − A0c (t))∆Ac (t) − 0 ˆ t f (l+1) dt2β(t)Ac (t)∆Ac (t) − 0 +λ(t)∆A2c (t) + β(t)∆A2c (t), (5.39) and ∆2s also follows the same form by substituting Ac (t) with As (t). ∆3 ≡ ˆ t f 0 (l) (l) dt[R(t, ~a(l+1) (t), Ac (t), As (t)) (l) (l) − R(t, ~a(l) (t), Ac (t), As (t))] ˆ t f = Re dth~b† (t)(H 0 + H c Aic (t) + H s Ais (t))∆~a(t)i 0 # ˆ t f ∂~b† (t) − dth )∆~a(t)i . ∂t 0 (l+1) Here ∆Ac = Ac (5.40) (l) (t) − Ac (t) and ∆~a(t) = ~a(l+1) (t) − ~a(l) (t) , and in deriving these expressions we also have chosen ~b(t) ≡ ~b(l) (t) in the expression for J¯(l+1) as well as in that for J¯(l) . We can write the equation of motion and initial state of ~b(t) by making the choice ~b† (tf ) = ~g † , ∂~b† (t) = ~b† (t)(H 0 + H c Ac (t) + H s As (t)), ∂t (5.41) (5.42) Therefore we obtain the result ∆1 = 0, ∆3 = 0. (l+1) Finally, the control parameter Ac (5.43) l (t) can be chosen as Al+1 c (t) = (Ac (t)+δAc (t))/α(t) and   l 1 ∂H(t) l+1 ~ † δAc (t) = Re hb (t) ~a (t)i , 2λ(t) ∂Ac (t) 108 (5.44) which implies that ∆2 = (λ(t) + β(t))[(∆Ac )2 + (∆As )2 ] ≥ 0. (5.45) Here λ(t) is a positive function that can be decided empirically and α(t) = 1 + β(t)/λ(t) is 2 l+1 2 a normalization term, which ensures that (Al+1 c (t)) + (As (t)) = 1 is satisfied at every time t. This proves that the iteration algorithm for cost function J exhibits monotonic convergence, given the proper choice of the equation of motion and initial condition for ~b† (t), and using Eq. (5.44) for the optimal control parameter at the next time step. Note that A0c (t) in (l) the l + 1 iteration can be substituted with the optimal control parameter Ac (t) acquired in l the l iteration. Similarly, the result of the As (t) is of the form Al+1 s (t) = (As (t)+δAs (t))/α(t) and   l 1 ∂H(t) l+1 ~ † δAs (t) = Re hb (t)| |~a (t)i , 2λ(t) ∂As (t) 5.4 (5.46) Application to a CDW Model We illustrate the Krotov method above using a solvable spinless fermion model on hypercubic lattices described by a hopping term τij and a periodic potential Ui , which is equal to U on the A sublattice and equal to zero on the B sublattice (See chapter 2). H(t) = − X † τij (t)ci cj + (U − µ) ij X † X † ci ci , ci ci − µ i∈A (5.47) i∈B where the chemical potential is µ = U/2 for half filling, and the time dependence of the hopping parameter is as given in Eq. (5.2). This model provides a useful starting point in understanding of the electronic response of CDW systems to an ultrafast laser pulse [56, 55]. For U = 0, the time dependence of the energy levels in the problem is given by ~ ε(~k − eA(t)/~c) = εk Ac (t) + ε̄k As (t) 109 (5.48) where Ac , As are given in Eq. (5.6). In order to compare with previous work [56], we P t∗ cos(r k ), √ take the limit of a hypercubic lattice in infinite dimensions εk = − d=∞ 0 l l=1 d Pd=∞ t∗ ε̄k = − l=1 √ sin(r0 kl ). d When U is finite the Hamiltonian breaks up into 2 × 2 blocks for each value of k, as the unit cell is now of size 2. When the laser is off, the energy eigenvalues of the system are q given by εk± = U/2 − µ ± ε2k + U 2 /4, where εk are the eigenvalues for U = 0, and the destruction operators after diagonalization may be written in terms of the U = 0 destruction operators ck and ck+Q where Q = (π, π, π, ...). The energy spectrum breaks up into two bands with band gap U , and the destruction operators are ck+ = αk ck + βk ck+Q for states in the conduction band and ck− = βk ck − αk ck+Q for states in the valence band. To impose the condition that an arbitrary target state with quantum number q+ is occupied at the † † final time tf we need hcq+ cq+ i = 1 and hcq− cq− i = 0. This is achieved by defining the goal P † operator ~g = q cq+ cq− , where the sum is over the target excited states. In the remainder of our discussion, the units are taken as r0 = t∗ = U = 1, the speed of light and Planck’s constant are in natural units, the chemical potential µ = U/2, the temperature of the initial state is T = 0, while the time and energy steps are ∆t = 0.01 and ∆ε ≤ 0.02, as used in previous work [56] on this model. The equations of motion for the operators in momentum space are given by, dck (t) = [εk (t) + U2 − µ]ck (t) + U2 ck+Q (t) dt dck+Q (t) i~ = U2 ck (t) + [−εk (t) + U2 − µ]ck+Q (t) dt i~ (5.49) (5.50) ~ where k is a wave vector and and εk (t) ≡ ε(~k − eA(t)/~c). It is useful to introduce the more compact notation ~ck (t) = (ck (t), ck+Q (t)) and a superposition operator ~a(t) = ⊕k~ck (t). The time evolution in Eq. (5.49) and (5.50) can then be written as i~∂t~ck (t) = H k (t)~ck (t) or L i~∂t~a(t) = H(t)~a(t) with H(t) = k H k (t), where H k (t) is deduced by comparing with Eq. (5.49) and Eq. (5.50) [55, 56]. The calculations are directed at optimizing the fidelity from Eq. (5.27), and we terminate the calculations when 1 − F is smaller than 10−3 or when the difference ratio |F l+1 − F l |/|1 − F l | is smaller than 10−3 . The initial condition before 110 arrival of the laser pulse is taken to be the half filled state so that the Fermi energy is at 1.0 (a) 0.5 0.0 0.5 1.0 1.0 (b) 0.5 0.0 0.5 1.0 0.4 0.2 (c) 0.0 0.2 0.4 0.6 0.8 40 ∗ ∗ Energy (unit/ tof t ) Energy A(t) A s (t) A c (t) zero energy and in the middle of the band gap, which is U in this model. 20 0 20 40 60 80 3 0.368 2 0.135 1 0.050 0 0.018 -1 0.007 0.003 -2 -3 -40 -20 0 20 40 Tim e (h/ tof∗)h/ t ∗ ) Tim e (unit 60 80 Figure 5.1: Time-resolved PES of the electronic CDW system before, during and after application of the optimal pulse (Ac (t), As (t), A(t)) from time t = 0 to t = 20 h/t∗ with the initial state a filled lower band at temperature T = 0. The target is excitation of a single state with the goal energy 1t∗ . In this calculation λ(t) = 2πh/∆ is a constant, with ∆ = 0.02 and h = 100. We calculate several physical properties in the final state that is produced by the optimal pulse, and we first focus on the time resolved photoemission spectra (tr-PES), P (ω, t), where i ´t P h´ t iω(t −t ) < 2 1 P (ω, t) = Im α −∞ dt2 −∞ dt1 s(t1 )s(t2 )e Gα,α (t2 , t1 ) , where s(t) is the probe 111 † laser pulse shape and α = A, B. G< α,β (t1 , t2 ) = ihcα (t2 )cβ (t1 )i is the lesser Green’s function; α and β can be either of the two sublattices A, B of a hypercubic lattice; and t1 , t2 are two times on the Keldysh contour (see [55] for details). The ultrafast time resolved PES presented in Fig. (5.1) is found when we choose an excitation band centered on Energy t∗ . In this calculation U = t∗ which is the band gap. This solvable system consists of a set of two-level systems (TLS) labeled by k with energy splitting (for q = k) δq = q+ − q− ~ ~ q · E/~ is the with associated generalized Rabi frequency Ωq = (χ2q + ∆2q )1/2 where χq = µ ~ the electric field and the detuning is Rabi frequency, µ ~ q the transition dipole moment, E ∆q = ω − δq /~ with ω the laser frequency. The optimal pulse found in this case (see Fig. (5.1) is in the weak field limit so that the dominant frequency is 2 which corresponds to the excitation energy δq for the target TLS. However some TLS near the target are also excited and the optimal pulse attempts to minimize excitation of TLS other than the target. The resulting optimal pulse is clearly not single mode with a Gaussian envelope, reflecting the fact that the target state is embedded in a continuum. Fig. (5.2) presents the tr-PES and optimal pulse found when the target is a fully excited state, i.e. excitation of all electrons from the valence band to the conduction band. This is an extreme excitation condition that can be interpreted as a state with a very large negative temperature. The QOCT method yields a laser pulse where the primary mode no longer corresponds to a single photon excitation of the TLS at the midband condition. Instead the optimal pulse is in the strong field limit where higher order processes are dominant. Moreover the optimal pulse is complex and asymmetric. An intriguing feature of this system is that the optimal control procedure acts on the parameters Ac (t), As (t) which are harmonic functions, and the form of the vector potential A(t) is deduced from them. This leads to a multitude of solutions for A(t) → A(t) + 2nπ (n is an integer) corresponding to the same functions As (t) and Ac (t). The result for A(t) presented in Fig. (5.2) is chosen to be the most continuous amongst this family of solutions. In Fig.5.3 we present several quantities to measure the time evolution of the real space and 112 momentum space occupancies of the CDW model. The first panel shows that the occupancy of the valence band (n− (t)) and conduction band (n+ (t)) are almost completely inverted by the optimal pulse. This is also evident in the real space occupancies, as in the initial state the electron occupancy on sublattice A (nA (t)) is much lower than that on sublattice B (nB (t)) due to the higher potential (U > 0) on the A sublattice. The order parameter Ω(t) = (nB (t) − nA (t))/(nB (t) + nA (t)) is presented in the second panel showing that the optimal laser pulse fully “melts" the electronic CDW state producing an almost complete inversion of the sublattice occupancy. The third panel gives the total energy of the electronic system as a function of time (hH(t)i), illustrating the fact that the energy in the electronic system is not monotonically increasing as the ensemble of Rabi oscillators in the model both absorb and emit photons into the radiation field. 5.5 Conclusion To conclude, we developed a new QOCT for tight binding models. Effective laser pulses have been found for the Photo-induced phase transitions (PIPT) in non-equilibruim timedependent charge density wave (CDW) systems based on a promising Quantum optimalcontrol theory (QOCT). The optical and structural properties as well as the temporal evolution of such states provide insight into the mutual dependence of electronic and atomic structure. Through simulations based on non-equilibrium Keldysh Green’s functions we find that the optimal laser pulse can achieve population inversion PIPT. This novel study can apply to many other complex materials and may-body dynamic problems. 113 A c (t) A s (t) A(t) ∗ Energy / t of Energy (unit t∗) 1.0 (a) 0.5 0.0 0.5 1.0 1.0 (b) 0.5 0.0 0.5 1.0 10 (c) 5 0 5 10 40 20 0 20 40 60 80 3 0.368 2 0.135 1 0.050 0 0.018 -1 0.007 0.003 -2 -3 -40 -20 0 20 40 e (h/ of t ∗)h/ t ∗ ) TimTim e (unit 60 80 Figure 5.2: Time-resolved PES, with the optimal pulse interacting with the system from time t = 0 to t = 20 h/t∗ . The initial state is a filled valence band at temperature T = 0, and the goal is to excite all of the states in the conduction band. The electron density at a given energy is plotted in false color and is high near the band edges as this model has a square root singularity in the density of states at both the valence band and conduction band edges. In this calculation λ(t) = 2πh/∆ is a constant, with ∆ = 0.02 and h = 10. 114 n ± (t) 0.50 n+ n- (a) 0.25 (t) 0.00  0.5 (b) 0.0  0.5 H (t) 0.4 (c) 0.0  0.4  40  20 0 20 40 60 80 ∗ TimTim e (unit e (h/oft ∗h/) t ) Figure 5.3: Properties of the system with the optimal pulse acting on the system from time t = 0 to t = 20 h/t∗ (a) The average occupancy of the conduction and valence bands n± (t), (b) The real space order parameter Ω(t), which measures the difference in occupancy of the two sublattices, (c) The energy hH(t)i of the CDW system. 115 CHAPTER 6 SUMMARY AND OUTLOOK This thesis presents two main contributions to the field of photo-induced topological states of matter. (i) A new variant of the SSH model is utilized to show how the interplay of polarized light and lattice structure can lead to tunable topological properties. (ii) A new quantum optimal control method was developed and demonstrated for a simple model of charge density wave response at femtosecond time scales. Moreover, we find that short optical pulses can lead to local spectral and novel pseudospin textures in one-dimensional topological insulators such as the SSH model. Pump-probe photoemission spectroscopy can probe these states and can demonstrate tunable energy gaps and Floquet band formation on femtosecond time scales. Analyzing band structures and pseudospin textures, we identify new states with optically induced nontrivial changes of sublattice mixing that leads to novel topological phenomenon. This study reveals the possibility to discover new topological phases driven by optical pulses in many insulators. As a first step toward optimal control of topological states, effective laser pulses have been found for the photo-induced phase transitions (PIPT) in non-equilibruim time-dependent charge density wave (CDW) systems based on a new and promising Quantum Optimal Control Theory (QOCT). Through simulations based on non-equilibrium Keldysh Green’s function we find the optimal laser pulse can achieve a high degree of population inversion. This new QOCT approach is applicable to many other complex materials and many-body dynamic problems. Future work has a variety of possible directions: One can include disorder terms into the tight-binding models [94, 95]. The effect of disorder may localize electrons and induce 116 Anderson localisation. Topological order may be influenced by the disorder effects if it is large enough to break the topologically protected states. Moreover, the interaction between light and disorder is also an interesting problem. Using disorder to describe noise in the system can also lead to real applications in quantum devices. Adding phonon effects to the tight binding models is already a very vivid field in physics. The effect of electron-phonon interactions on optical absorption spectra requires a special treatment. By incorporating key phonon modes into the charge density wave response or in topological materials, many novel responses and new physics have been discovered [96, 97, 98, 99]. Tuning these response with optimized light pulses is a promising future direction. Interaction effects between electrons is important for many charge density wave systems and for some topological orders so this is an important and challenging topic to study [100, 101, 102]. By introducing the Hubbard model, which is an approximate model used in solid state physics to describe the transition between conducting and insulating systems, one can discuss electron behaviors with Hubbard interaction and discover new phenomena in charge density wave and topological phases; and the effect of light on these phenomena. Finally quantum optimal control is important for operating quantum dots in quantum computing and advanced quantum devices. One of the possible applications is to control Majorana edge states in the time domain, and in general to understand and control FloquetBloch topological states of matter [103, 104]. 117 CHAPTER 7 EXPERIMENTS AND APPLICATIONS Recently many novel experiments have been designed for studying the field of photo-induced phase transitions and topological states of matter. In this Chapter a brief overview of four very different experimental systems are described. In Ref. [7], they used trARPES to study photo-induced band gaps. When the photon energy of the laser is less than 300 meV, the coherent interaction between light and the TI surface states is the main effect. Therefore they used polarized photons at midinfrared (MIR) wavelengths to investigate the photon-dressed surface states in Bi2 Se3 . The laser pulses are focused to a 300 mm diameter spot on the single-crystal Bi2 Se3 sample at an angle of 45◦ , and the pulses are tunable in wavelength from 4 mm to 17 mm with 1 mJ peak energy; with the amplitude of the electric field to be E0 = 2.5 (T1) × 107 V/m and the pulsewidth estimated to be 250 fs (FWHM) on the surface of the sample: Fig. 7.1 shows energy-momentum spectra of Bi2 Se3 obtained at several time delays after arrival of the intense linearly polarized MIR excitation. Note that the probe pulse is a linearly polarized ultraviolet (UV) pulse. Figure 7.1: Angle-resolved photoemission spectra (APRES) of Bi2 Se3 . (A) A sketch of the experimental geometry for the p-polarized case. kx is defined to be the in-plane electron momentum parallel to the pump scattering plane. (B to F) ARPES data for several pumpprobe time delays t (values indicated in the figure) under strong linearly polarized midinfrared (MIR) excitation of wavelength l = 10 mm. Figure and caption are taken from [7]. 118 In Ref. [8], arrays of evanescently coupled dielectric-loaded surface-plasmon polariton (SPP) waveguides was used to study the SSH model. The arrays are built by negativetone gray-scale electron beam lithography on top of a chromium (10 nm) and gold (60 nm) coated glass substrate. The waveguides consist of polymethylmethacrylate (PMMA) ridges with a height of 140 nm and a width of 250 nm. They used weak and strong bonds (1000 and 600 nm) as in the SSH model by alternating different separations between neighboring waveguides, as depicted in Fig. 7.2, giving in a = 1600 nm for the SSH unit cell separation distance. These geometrical parameters ensure, for the vacuum wavelength λ = 980 nm, single-mode operation of the waveguides and sufficient coupling among them. In Ref. [9], they studied on the imaging and probing of topological bound states in the SSH model through adiabatic preparation, quench dynamics and phase-sensitive injection by using an atom-optics realization of lattice tight-binding models. With an optical lattice potential formed by lasers of wave number k = 2π/λ and wavelength λ = 1064 nm, momentum-space dynamics of Rb87 condensate atoms was started through controlled, time-dependent driving. The lasers couple 21 discrete atomic momentum states coherently, inducing a momentum-space lattice of states such that atomic population may change. The momentum states are defined by momenta pn = 2n~k and site indices n(relative to the lowest momentum value). Through 20 distinct two-photon Bragg diffraction processes, the coupling between these states is fully controlled, letting them to simulate tight-binding models with local, arbitrary and time-dependent control of all tunnelling terms and site energies. Fig. 7.3 is the absorption images taken after 760µs(≈ 0.78h/t) of evolution following the initialization and quench, corresponding to phases of π (top) and 0 (bottom), respectively, for ∆/t = 0.38(1). In Ref. [10], photonic waveguide arrays were used (see Fig. 7.4), which are an exceptional system for implementing PT-symmetric non-hermitian Hamiltonians. By wiggling 119 Figure 7.2: (a) Waveguide array fabricated out of PMMA on top of a Cr- and Au-coated glass substrate. Alternating center-to-center separations, 600 and 1000 nm, implement the bulk SSH model. (b) Plasmonic waveguide array incorporating a topological defect where the long separation is repeated twice. Three different excitation sites, I, II, and III, are highlighted. Figure and caption are taken from [8]. the waveguides, which causes radiative loss (coupling to the continuum states), the loss in photonic lattices can be precisely designed. The frequency and amplitude of the wiggling can be tuned to generate a particular loss. The four examples above are a small subset of experimental systems where photo-induced phase transitions and photo-induced topological properties are being explored; as outlined in Chapter 1. 120 Figure 7.3: Absorption images taken after 760µs(≈ 0.78h/t) of evolution following the initialization and quench, corresponding to phases of π (top) and 0 (bottom), respectively, for ∆/t = 0.38(1). Figure and caption are taken from [9]. Figure 7.4: (a) Sketch of the passive waveguide array acting like a PT-symmetric structure with a topological interface. (b) End facet of the experimentally realized structure in fused silica glass. Image courtesy of Mark Kremer, FSU Jena. Figure and caption are taken from [10]. 121 BIBLIOGRAPHY 122 BIBLIOGRAPHY [1] F. Schmitt, P. S. Kirchmann, U. Bovensiepen, R. G. Moore, L. Rettig, M. Krenz, J.-H. Chu, N. Ru, L. Perfetti, D. H. Lu, M. Wolf, I. R. Fisher, and Z.-X. Shen. Transient electronic structure and melting of a charge density wave in tbte3. Science, 321(5896):1649–1652, 2008. ISSN 0036-8075. doi: 10.1126/science.1160778. URL http://science.sciencemag.org/content/321/5896/1649. [2] J. D. Rameau, S. Freutel, L. Rettig, I. Avigo, M. Ligges, Y. Yoshida, H. Eisaki, J. Schneeloch, R. D. Zhong, Z. J. Xu, G. D. Gu, P. D. Johnson, and U. Bovensiepen. Photoinduced changes in the cuprate electronic structure revealed by femtosecond time- and angle-resolved photoemission. Phys. Rev. B, 89:115115, Mar 2014. doi: 10.1103/PhysRevB.89.115115. URL https://link.aps.org/doi/10.1103/PhysRevB. 89.115115. [3] L Perfetti, P A Loukakos, M Lisowski, U Bovensiepen, M Wolf, H Berger, S Biermann, and A Georges. Femtosecond dynamics of electronic states in the mott insulator 1t-tas 2 by time resolved photoelectron spectroscopy. New Journal of Physics, 10(5):053019, 2008. URL http://stacks.iop.org/1367-2630/10/i=5/a=053019. [4] Herschel Rabitz, Regina de Vivie-Riedle, Marcus Motzkus, and Karl Kompa. Whither the future of controlling quantum phenomena? Science, 288(5467):824–828, 2000. ISSN 0036-8075. doi: 10.1126/science.288.5467.824. URL http://science.sciencemag. org/content/288/5467/824. [5] M. Z. Hasan and C. L. Kane. Colloquium : Topological insulators. Rev. Mod. Phys., 82:3045–3067, Nov 2010. doi: 10.1103/RevModPhys.82.3045. URL http://link.aps. org/doi/10.1103/RevModPhys.82.3045. [6] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik. Observation of floquetbloch states on the surface of a topological insulator. Science, 342(6157):453–457, 2013. ISSN 0036-8075. doi: 10.1126/science.1239834. URL http://science.sciencemag.org/ content/342/6157/453. [7] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik. Observation of floquetbloch states on the surface of a topological insulator. Science, 342(6157):453–457, 2013. ISSN 0036-8075. doi: 10.1126/science.1239834. URL http://science.sciencemag.org/ content/342/6157/453. [8] Felix Bleckmann, Zlata Cherpakova, Stefan Linden, and Andrea Alberti. Spectral imaging of topological edge states in plasmonic waveguide arrays. Phys. Rev. B, 96: 045417, Jul 2017. doi: 10.1103/PhysRevB.96.045417. URL https://link.aps.org/doi/ 10.1103/PhysRevB.96.045417. [9] Eric J. Meier, Fangzhao Alex An, and Bryce Gadway. Observation of the topological soliton state in the su-schrieffer-heeger model. 7:13986 EP –, Dec 2016. URL http: //dx.doi.org/10.1038/ncomms13986. Article. 123 [10] S. Weimann, M. Kremer, Y. Plotnik, Y. Lumer, S. Nolte, K. G. Makris, M. Segev, M. C. Rechtsman, and A. Szameit. Topologically protected bound states in photonic parity-time-symmetric crystals. Nat Mater, 16(4):433–438, Apr 2017. ISSN 1476-1122. URL http://dx.doi.org/10.1038/nmat4811. Article. [11] A. L. Schawlow and C. H. Townes. Infrared and optical masers. Phys. Rev., 112: 1940–1949, Dec 1958. doi: 10.1103/PhysRev.112.1940. URL https://link.aps.org/doi/ 10.1103/PhysRev.112.1940. [12] Maiman T. H. Nature, 187(4736):493, 1960. URL http://dx.doi.org/10.1038/187493a0. [13] P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich. Generation of optical harmonics. Phys. Rev. Lett., 7:118–119, Aug 1961. doi: 10.1103/PhysRevLett.7.118. URL https://link.aps.org/doi/10.1103/PhysRevLett.7.118. [14] Kun Zhao, Qi Zhang, Michael Chini, Yi Wu, Xiaowei Wang, and Zenghu Chang. Tailoring a 67 attosecond pulse through advantageous phase-mismatch. Opt. Lett., 37(18): 3891–3893, Sep 2012. doi: 10.1364/OL.37.003891. URL http://ol.osa.org/abstract. cfm?URI=ol-37-18-3891. [15] Christopher L. Smallwood, Robert A. Kaindl, and Alessandra Lanzara. Ultrafast angle-resolved photoemission spectroscopy of quantum materials. EPL (Europhysics Letters), 115(2):27001, 2016. URL http://stacks.iop.org/0295-5075/115/i=2/a=27001. [16] J. BOKOR. Ultrafast dynamics at semiconductor and metal surfaces. Science, 246 (4934):1130–1134, 1989. ISSN 0036-8075. doi: 10.1126/science.246.4934.1130. URL http://science.sciencemag.org/content/246/4934/1130. [17] Richard Haight. Electron dynamics at surfaces. Surface Science Reports, 21(8):275 – 325, 1995. ISSN 0167-5729. doi: http://dx.doi.org/10.1016/0167-5729(95)00002-X. URL http://www.sciencedirect.com/science/article/pii/016757299500002X. [18] H. Petek and S. Ogawa. Femtosecond time-resolved two-photon photoemission studies of electron dynamics in metals. Progress in Surface Science, 56(4):239 – 310, 1997. ISSN 0079-6816. doi: http://dx.doi.org/10.1016/S0079-6816(98)00002-1. URL http: //www.sciencedirect.com/science/article/pii/S0079681698000021. [19] M. Bauer, C. Lei, K. Read, R. Tobey, J. Gland, M. M. Murnane, and H. C. Kapteyn. Direct observation of surface chemistry using ultrafast soft-x-ray pulses. Phys. Rev. Lett., 87:025501, Jun 2001. doi: 10.1103/PhysRevLett.87.025501. URL https://link. aps.org/doi/10.1103/PhysRevLett.87.025501. [20] P. Siffalovic, M. Drescher, M. Spieweck, T. Wiesenthal, Y. C. Lim, R. Weidner, A. Elizarov, and U. Heinzmann. Laser-based apparatus for extended ultraviolet femtosecond time-resolved photoemission spectroscopy. Review of Scientific Instruments, 72(1):30–35, 2001. doi: 10.1063/1.1329904. URL http://dx.doi.org/10.1063/1.1329904. 124 [21] S. Mathias, L. Miaja-Avila, M. M. Murnane, H. Kapteyn, M. Aeschlimann, and M. Bauer. Angle-resolved photoemission spectroscopy with a femtosecond high harmonic light source using a two-dimensional imaging electron analyzer. Review of Scientific Instruments, 78(8):083105, 2007. doi: 10.1063/1.2773783. URL http: //dx.doi.org/10.1063/1.2773783. [22] G. L. Dakovski, Y. Li, T. Durakiewicz, and G. Rodriguez. Tunable ultrafast extreme ultraviolet source for time- and angle-resolved photoemission spectroscopy. Review of Scientific Instruments, 81(7):073108, 2010. doi: 10.1063/1.3460267. URL http: //dx.doi.org/10.1063/1.3460267. [23] Christopher L. Smallwood, Christopher Jozwiak, and Alessandra Lanzara. An ultrafast angle-resolved photoemission apparatus for measuring complex materials. Review of Scientific Instruments, 83(12):123904, 2012. doi: 10.1063/1.4772070. URL http://dx. doi.org/10.1063/1.4772070. [24] J. A. Sobota, S. Yang, J. G. Analytis, Y. L. Chen, I. R. Fisher, P. S. Kirchmann, and Z.-X. Shen. Ultrafast optical excitation of a persistent surface-state population in the topological insulator bi2 se3 . Phys. Rev. Lett., 108:117403, Mar 2012. doi: 10. 1103/PhysRevLett.108.117403. URL https://link.aps.org/doi/10.1103/PhysRevLett. 108.117403. [25] Hüfner S. Photoelectron Spectroscopy: Principles and Applications, 2003. URL http: //dx.doi.org/10.1007/978-3-662-09280-4. [26] Andrea Damascelli, Zahid Hussain, and Zhi-Xun Shen. Angle-resolved photoemission studies of the cuprate superconductors. Rev. Mod. Phys., 75:473–541, Apr 2003. doi: 10.1103/RevModPhys.75.473. URL https://link.aps.org/doi/10.1103/RevModPhys. 75.473. [27] M. Z. Hasan and C. L. Kane. Colloquium. Rev. Mod. Phys., 82:3045–3067, Nov 2010. doi: 10.1103/RevModPhys.82.3045. URL https://link.aps.org/doi/10.1103/ RevModPhys.82.3045. [28] A. A. Kordyuk. Arpes experiment in fermiology of quasi-2d metals (review article). Low Temperature Physics, 40(4):286–296, 2014. doi: 10.1063/1.4871745. URL http: //dx.doi.org/10.1063/1.4871745. [29] D. N. Basov, M. M. Fogler, A. Lanzara, Feng Wang, and Yuanbo Zhang. Colloquium. Rev. Mod. Phys., 86:959–994, Jul 2014. doi: 10.1103/RevModPhys.86.959. URL https: //link.aps.org/doi/10.1103/RevModPhys.86.959. [30] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Eisaki, and M. Wolf. Ultrafast electron relaxation in superconducting bi2 sr2 cacu2 o8+δ by time-resolved photoelectron spectroscopy. Phys. Rev. Lett., 99:197001, Nov 2007. doi: 10. 1103/PhysRevLett.99.197001. URL https://link.aps.org/doi/10.1103/PhysRevLett. 99.197001. 125 [31] Philip B. Allen. Theory of thermal relaxation of electrons in metals. Phys. Rev. Lett., 59:1460–1463, Sep 1987. doi: 10.1103/PhysRevLett.59.1460. URL https://link.aps. org/doi/10.1103/PhysRevLett.59.1460. [32] J. Demsar, L. Forró, H. Berger, and D. Mihailovic. Femtosecond snapshots of gapforming charge-density-wave correlations in quasi-two-dimensional dichalcogenides 1t− tas2 and 2h − tase2 . Phys. Rev. B, 66:041101, Jun 2002. doi: 10.1103/PhysRevB.66. 041101. URL http://link.aps.org/doi/10.1103/PhysRevB.66.041101. [33] V. Brouet, W. L. Yang, X. J. Zhou, Z. Hussain, N. Ru, K. Y. Shin, I. R. Fisher, and Z. X. Shen. Fermi surface reconstruction in the cdw state of CeTe3 observed by photoemission. Phys. Rev. Lett., 93:126405, Sep 2004. doi: 10.1103/PhysRevLett.93. 126405. URL http://link.aps.org/doi/10.1103/PhysRevLett.93.126405. [34] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf. Time evolution of the electronic structure of 1t−tas2 through the insulator-metal transition. Phys. Rev. Lett., 97:067402, Aug 2006. doi: 10.1103/PhysRevLett.97.067402. URL https://link.aps.org/doi/10.1103/ PhysRevLett.97.067402. [35] S. Hellmann, M. Beye, C. Sohrt, T. Rohwer, F. Sorgenfrei, H. Redlin, M. Kalläne, M. Marczynski-Bühlow, F. Hennies, M. Bauer, A. Föhlisch, L. Kipp, W. Wurth, and K. Rossnagel. Ultrafast melting of a charge-density wave in the mott insulator 1t-tas2 . Phys. Rev. Lett., 105:187401, Oct 2010. doi: 10.1103/PhysRevLett.105.187401. URL http://link.aps.org/doi/10.1103/PhysRevLett.105.187401. [36] F Schmitt, P S Kirchmann, U Bovensiepen, R G Moore, J-H Chu, D H Lu, L Rettig, M Wolf, I R Fisher, and Z-X Shen. Ultrafast electron dynamics in the charge density wave material tbte 3. New Journal of Physics, 13(6):063022, 2011. URL http://stacks. iop.org/1367-2630/13/i=6/a=063022. [37] J. C. Petersen, S. Kaiser, N. Dean, A. Simoncig, H. Y. Liu, A. L. Cavalieri, C. Cacho, I. C. E. Turcu, E. Springate, F. Frassetto, L. Poletto, S. S. Dhesi, H. Berger, and A. Cavalleri. Clocking the melting transition of charge and lattice order in 1t−tas2 with ultrafast extreme-ultraviolet angle-resolved photoemission spectroscopy. Phys. Rev. Lett., 107:177402, Oct 2011. doi: 10.1103/PhysRevLett.107.177402. URL https: //link.aps.org/doi/10.1103/PhysRevLett.107.177402. [38] Timm Rohwer, Stefan Hellmann, Martin Wiesenmayer, Christian Sohrt, Ankatrin Stange, Bartosz Slomski, Adra Carr, Yanwei Liu, Luis Miaja Avila, Matthias Kallane, Stefan Mathias, Lutz Kipp, Kai Rossnagel, and Michael Bauer. Collapse of longrange charge order tracked by time-resolved photoemission at high momenta. Nature, 471(7339):490–493, Mar 2011. ISSN 0028-0836. doi: 10.1038/nature09829. URL http://dx.doi.org/10.1038/nature09829. [39] Maximilian Eichberger, Hanjo Schaefer, Marina Krumova, Markus Beyer, Jure Demsar, Helmuth Berger, Gustavo Moriena, German Sciaini, and R. J. Dwayne Miller. 126 Snapshots of cooperative atomic motions in the optical suppression of charge density waves. NATURE, 468(7325):799–802, DEC 9 2010. ISSN 0028-0836. doi: 10.1038/nature09539. [40] Georgi L. Dakovski, Yinwan Li, Steve M. Gilbertson, George Rodriguez, Alexander V. Balatsky, Jian-Xin Zhu, Krzysztof Gofryk, Eric D. Bauer, Paul H. Tobash, Antoinette Taylor, John L. Sarrao, Peter M. Oppeneer, Peter S. Riseborough, John A. Mydosh, and Tomasz Durakiewicz. Anomalous femtosecond quasiparticle dynamics of hidden order state in uru2 si2 . Phys. Rev. B, 84:161103, Oct 2011. doi: 10.1103/PhysRevB. 84.161103. URL https://link.aps.org/doi/10.1103/PhysRevB.84.161103. [41] Isabella Gierz, Jesse C. Petersen, Matteo Mitrano, Cephise Cacho, I. C. Edmond Turcu, Emma Springate, Alexander Stöhr, Axel Köhler, Ulrich Starke, and Andrea Cavalleri. Snapshots of non-equilibrium dirac carrier distributions in graphene. Nat Mater, 12(12):1119–1124, Dec 2013. ISSN 1476-1122. URL http://dx.doi.org/10.1038/ nmat3757. Article. [42] Jens Christian Johannsen, Søren Ulstrup, Federico Cilento, Alberto Crepaldi, Michele Zacchigna, Cephise Cacho, I. C. Edmond Turcu, Emma Springate, Felix Fromm, Christian Raidel, Thomas Seyller, Fulvio Parmigiani, Marco Grioni, and Philip Hofmann. Direct view of hot carrier dynamics in graphene. Phys. Rev. Lett., 111:027403, Jul 2013. doi: 10.1103/PhysRevLett.111.027403. URL https://link.aps.org/doi/10.1103/ PhysRevLett.111.027403. [43] Steve M. Gilbertson, Tomasz Durakiewicz, Georgi L. Dakovski, Yinwan Li, Jian-Xin Zhu, Steven D. Conradson, Stuart A. Trugman, and George Rodriguez. Ultrafast photoemission spectroscopy of the uranium dioxide uo2 mott insulator: Evidence for a robust energy gap structure. Phys. Rev. Lett., 112:087402, Feb 2014. doi: 10. 1103/PhysRevLett.112.087402. URL https://link.aps.org/doi/10.1103/PhysRevLett. 112.087402. [44] Søren Ulstrup, Jens Christian Johannsen, Federico Cilento, Jill A. Miwa, Alberto Crepaldi, Michele Zacchigna, Cephise Cacho, Richard Chapman, Emma Springate, Samir Mammadov, Felix Fromm, Christian Raidel, Thomas Seyller, Fulvio Parmigiani, Marco Grioni, Phil D. C. King, and Philip Hofmann. Ultrafast dynamics of massive dirac fermions in bilayer graphene. Phys. Rev. Lett., 112:257401, Jun 2014. doi: 10. 1103/PhysRevLett.112.257401. URL https://link.aps.org/doi/10.1103/PhysRevLett. 112.257401. [45] Timm Rohwer, Stefan Hellmann, Martin Wiesenmayer, Christian Sohrt, Ankatrin Stange, Bartosz Slomski, Adra Carr, Yanwei Liu, Luis Miaja Avila, Matthias Kallane, Stefan Mathias, Lutz Kipp, Kai Rossnagel, and Michael Bauer. Collapse of longrange charge order tracked by time-resolved photoemission at high momenta. Nature, 471(7339):490–493, Mar 2011. ISSN 0028-0836. doi: 10.1038/nature09829. URL http://dx.doi.org/10.1038/nature09829. [46] S. Hellmann, T. Rohwer, M. Kalläne, K. Hanff, C. Sohrt, A. Stange, A. Carr, M. M. Murnane, H. C. Kapteyn, L. Kipp, M. Bauer, and K. Rossnagel. Time127 domain classification of charge-density-wave insulators. 3:1069 EP –, Sep 2012. URL http://dx.doi.org/10.1038/ncomms2078. Article. [47] Christopher L. Smallwood, Wentao Zhang, Tristan L. Miller, Chris Jozwiak, Hiroshi Eisaki, Dung-Hai Lee, and Alessandra Lanzara. Time- and momentum-resolved gap dynamics in bi2 sr2 cacu2 o8+δ . Phys. Rev. B, 89:115126, Mar 2014. doi: 10.1103/ PhysRevB.89.115126. URL https://link.aps.org/doi/10.1103/PhysRevB.89.115126. [48] J. Demsar, B. Podobnik, V. V. Kabanov, Th. Wolf, and D. Mihailovic. Superconducting gap ∆c , the pseudogap ∆p , and pair fluctuations above Tc in overdoped Y 1−x Ca x Ba 2 Cu 3 O 7−δ from femtosecond time-domain spectroscopy. Phys. Rev. Lett., 82:4918–4921, Jun 1999. doi: 10.1103/PhysRevLett.82.4918. URL https: //link.aps.org/doi/10.1103/PhysRevLett.82.4918. [49] R. Peierls. Zur theorie der elektrischen und thermischen leitfähigkeit von metallen. Annalen der Physik, 396(2):121–148, 1930. ISSN 1521-3889. doi: 10.1002/andp. 19303960202. URL http://dx.doi.org/10.1002/andp.19303960202. [50] G. Grüner. The dynamics of charge-density waves. Rev. Mod. Phys., 60:1129–1181, Oct 1988. doi: 10.1103/RevModPhys.60.1129. URL https://link.aps.org/doi/10.1103/ RevModPhys.60.1129. [51] L. Perfetti, P. A. Loukakos, M. Lisowski, U. Bovensiepen, H. Berger, S. Biermann, P. S. Cornaglia, A. Georges, and M. Wolf. Time evolution of the electronic structure of 1t-tas2 through the insulator-metal transition. Phys. Rev. Lett., 97:067402, Aug 2006. doi: 10.1103/PhysRevLett.97.067402. URL http://link.aps.org/doi/10.1103/ PhysRevLett.97.067402. [52] I. M. Vishik, W. S. Lee, F. Schmitt, B. Moritz, T. Sasagawa, S. Uchida, K. Fujita, S. Ishida, C. Zhang, T. P. Devereaux, and Z. X. Shen. Doping-dependent nodal fermi velocity of the high-temperature superconductor bi2 sr2 cacu2 o8+δ revealed using highresolution angle-resolved photoemission spectroscopy. Phys. Rev. Lett., 104:207002, May 2010. doi: 10.1103/PhysRevLett.104.207002. URL https://link.aps.org/doi/10. 1103/PhysRevLett.104.207002. [53] A. Castro, J. Werschnik, and E. K. U. Gross. Controlling the dynamics of manyelectron systems from first principles: A combination of optimal control and timedependent density-functional theory. Phys. Rev. Lett., 109:153603, Oct 2012. doi: 10. 1103/PhysRevLett.109.153603. URL http://link.aps.org/doi/10.1103/PhysRevLett. 109.153603. [54] J. K. Freericks, V. M. Turkowski, and V. Zlatić. Nonequilibrium dynamical mean-field theory. Phys. Rev. Lett., 97:266408, Dec 2006. doi: 10.1103/PhysRevLett.97.266408. URL http://link.aps.org/doi/10.1103/PhysRevLett.97.266408. [55] Wen Shen, A. F. Kemper, T. P. Devereaux, and J. K. Freericks. Exact solution for high harmonic generation and the response to an ac driving field for a charge density wave insulator. Phys. Rev. B, 90:115113, Sep 2014. doi: 10.1103/PhysRevB.90.115113. URL http://link.aps.org/doi/10.1103/PhysRevB.90.115113. 128 [56] Wen Shen, Yizhi Ge, A. Y. Liu, H. R. Krishnamurthy, T. P. Devereaux, and J. K. Freericks. Nonequilibrium “melting” of a charge density wave insulator via an ultrafast laser pulse. Phys. Rev. Lett., 112:176404, May 2014. doi: 10.1103/PhysRevLett.112. 176404. URL http://link.aps.org/doi/10.1103/PhysRevLett.112.176404. [57] Patrick Doria, Tommaso Calarco, and Simone Montangero. Optimal control technique for many-body quantum dynamics. Phys. Rev. Lett., 106:190501, May 2011. doi: 10. 1103/PhysRevLett.106.190501. URL http://link.aps.org/doi/10.1103/PhysRevLett. 106.190501. [58] T. Caneva, A. Silva, R. Fazio, S. Lloyd, T. Calarco, and S. Montangero. Complexity of controlling quantum many-body dynamics. Phys. Rev. A, 89:042322, Apr 2014. doi: 10.1103/PhysRevA.89.042322. URL http://link.aps.org/doi/10.1103/PhysRevA. 89.042322. [59] Hideo Aoki, Naoto Tsuji, Martin Eckstein, Marcus Kollar, Takashi Oka, and Philipp Werner. Nonequilibrium dynamical mean-field theory and its applications. Rev. Mod. Phys., 86:779–837, Jun 2014. doi: 10.1103/RevModPhys.86.779. URL http://link.aps. org/doi/10.1103/RevModPhys.86.779. [60] Tzong-Ru T. Han, Zhensheng Tao, Subhendra D. Mahanti, Kiseok Chang, ChongYu Ruan, Christos D. Malliakas, and Mercouri G. Kanatzidis. Structural dynamics of two-dimensional charge-density waves in cete3 investigated by ultrafast electron crystallography. Phys. Rev. B, 86:075145, Aug 2012. doi: 10.1103/PhysRevB.86. 075145. URL http://link.aps.org/doi/10.1103/PhysRevB.86.075145. [61] N. Dean, J. C. Petersen, D. Fausti, R. I. Tobey, S. Kaiser, L. V. Gasparov, H. Berger, and A. Cavalleri. Polaronic conductivity in the photoinduced phase of 1t-tas2 . Phys. Rev. Lett., 106:016401, Jan 2011. doi: 10.1103/PhysRevLett.106.016401. URL http: //link.aps.org/doi/10.1103/PhysRevLett.106.016401. [62] Vadim F. Krotov. Global methods in optimal control theory. Marcel Dekker, Inc, 1996. [63] M. Wenin and W. Pötz. Minimization of environment-induced decoherence in quantum subsystems and application to solid-state-based quantum gates. Phys. Rev. B, 78(16): 165118, Oct 2008. doi: 10.1103/PhysRevB.78.165118. [64] Robert Roloff and Walter Pötz. Time-optimal performance of josephson charge qubits: A process tomography approach. Phys. Rev. B, 79(22):224516, Jun 2009. doi: 10.1103/ PhysRevB.79.224516. [65] Shabnam Safaei, Simone Montangero, Fabio Taddei, and Rosario Fazio. Optimized single-qubit gates for josephson phase qubits. Phys. Rev. B, 79(6):064524, Feb 2009. doi: 10.1103/PhysRevB.79.064524. [66] M. Wenin, R. Roloff, and W. Pötz. Robust control of josephson charge qubits. Journal of Applied Physics, 105(8):084504, 2009. doi: 10.1063/1.3080242. URL http://link. aip.org/link/?JAP/105/084504/1. 129 [67] Feng Shuang and Herschel Rabitz. Cooperating or fighting with decoherence in the optimal control of quantum dynamics. The Journal of Chemical Physics, 124(15):154105, 2006. doi: 10.1063/1.2186644. URL http://link.aip.org/link/?JCP/124/154105/1. [68] Jean Christophe Tremblay and Peter Saalfrank. Guided locally optimal control of quantum dynamics in dissipative environments. Phys. Rev. A, 78(6):063408, Dec 2008. doi: 10.1103/PhysRevA.78.063408. [69] Bin Hwang and Hsi-Sheng Goan. Optimal control for non-markovian open quantum systems. Phys. Rev. A, 85:032321, Mar 2012. doi: 10.1103/PhysRevA.85.032321. URL http://link.aps.org/doi/10.1103/PhysRevA.85.032321. [70] Jung-Shen Tai, Kuan-Ting Lin, and Hsi-Sheng Goan. Optimal control of quantum gates in an exactly solvable non-markovian open quantum bit system. Phys. Rev. A, 89:062310, Jun 2014. doi: 10.1103/PhysRevA.89.062310. URL http://link.aps.org/ doi/10.1103/PhysRevA.89.062310. [71] Reuven Eitan, Michael Mundt, and David J. Tannor. Optimal control with accelerated convergence: Combining the krotov and quasi-newton methods. Phys. Rev. A, 83: 053426, May 2011. doi: 10.1103/PhysRevA.83.053426. URL http://link.aps.org/doi/ 10.1103/PhysRevA.83.053426. [72] D. J. Tannor, V. Kazakov, and V. Orlov. Control of Photochemical Branching: Novel Procedures for Finding Optimal Pulses and Global Upper Bounds. In J. Broeckhove & L. Lathouwers, editor, NATO ASIB Proc. 299: Time-Dependent Quantum Molecular Dynamics, pages 347–360, 1992. [73] Ivan I. Maximov, Zdenĕk Tošner, and Niels Chr. Nielsen. Optimal control design of nmr and dynamic nuclear polarization experiments using monotonically convergent algorithms. The Journal of Chemical Physics, 128(18):184505, 2008. doi: 10.1063/1. 2903458. URL http://link.aip.org/link/?JCP/128/184505/1. [74] José P. Palao and Ronnie Kosloff. Quantum computing by an optimal control algorithm for unitary transformations. Phys. Rev. Lett., 89(18):188301, Oct 2002. doi: 10.1103/ PhysRevLett.89.188301. [75] José P. Palao and Ronnie Kosloff. Optimal control theory for unitary transformations. Phys. Rev. A, 68(6):062308, Dec 2003. doi: 10.1103/PhysRevA.68.062308. [76] Ruixue Xu, YiJing Yan, Yukiyoshi Ohtsuki, Yuichi Fujimura, and Herschel Rabitz. Optimal control of quantum non-markovian dissipation: Reduced liouville-space theory. The Journal of Chemical Physics, 120(14):6600–6608, 2004. doi: 10.1063/1.1665486. URL http://link.aip.org/link/?JCP/120/6600/1. [77] Fahad Mahmood, Ching-Kit Chan, Zhanybek Alpichshev, Dillon Gardner, Young Lee, Patrick A. Lee, and Nuh Gedik. Selective scattering between floquet-bloch and volkov states in a topological insulator. Nat Phys, 12, Apr 2016. 130 [78] G. Saathoff, L. Miaja-Avila, M. Aeschlimann, M. M. Murnane, and H. C. Kapteyn. Laser-assisted photoemission from surfaces. Phys. Rev. A, 77:022903, Feb 2008. doi: 10.1103/PhysRevA.77.022903. URL https://link.aps.org/doi/10.1103/PhysRevA.77. 022903. [79] J. A. Sobota, S.-L. Yang, A. F. Kemper, J. J. Lee, F. T. Schmitt, W. Li, R. G. Moore, J. G. Analytis, I. R. Fisher, P. S. Kirchmann, T. P. Devereaux, and Z.-X. Shen. Direct optical coupling to an unoccupied dirac surface state in the topological insulator bi2 se3 . Phys. Rev. Lett., 111:136802, Sep 2013. doi: 10.1103/PhysRevLett.111.136802. URL https://link.aps.org/doi/10.1103/PhysRevLett.111.136802. [80] J. A. Sobota, S. Yang, J. G. Analytis, Y. L. Chen, I. R. Fisher, P. S. Kirchmann, and Z.-X. Shen. Ultrafast optical excitation of a persistent surface-state population in the topological insulator bi2 se3 . Phys. Rev. Lett., 108:117403, Mar 2012. doi: 10.1103/PhysRevLett.108.117403. URL http://link.aps.org/doi/10.1103/ PhysRevLett.108.117403. [81] A. Gómez-León and G. Platero. Floquet-bloch theory and topology in periodically driven lattices. Phys. Rev. Lett., 110:200403, May 2013. doi: 10.1103/PhysRevLett. 110.200403. URL http://link.aps.org/doi/10.1103/PhysRevLett.110.200403. [82] V. Dal Lago, M. Atala, and L. E. F. Foa Torres. Floquet topological transitions in a driven one-dimensional topological insulator. Phys. Rev. A, 92:023624, Aug 2015. doi: 10.1103/PhysRevA.92.023624. [83] M. A. Sentef, M. Claassen, A. F. Kemper, B. Moritz, T. Oka, J. K. Freericks, and T. P. Devereaux. Theory of floquet band formation and local pseudospin textures in pump-probe photoemission of graphene. Nat Commun, 6, May 2015. Article. [84] R. van Leeuwen, N.E. Dahlen, G. Stefanucci, C.-O. Almbladh, and U. von Barth. Introduction to the Keldysh Formalism, pages 33–59. Springer Berlin Heidelberg, Berlin, Heidelberg, 2006. [85] Van-Nam Do. Non-equilibrium green function method: theory and application in simulation of nanometer electronic devices. Advances in Natural Sciences: Nanoscience and Nanotechnology, 5(3):033001, 2014. URL http://stacks.iop.org/2043-6262/5/i=3/ a=033001. [86] M. V. Berry. Quantal phase factors accompanying adiabatic changes. 392(1802):45–57, 1984. ISSN 0080-4630. doi: 10.1098/rspa.1984.0023. [87] J. Zak. Berry’s phase for energy bands in solids. Phys. Rev. Lett., 62:2747–2750, Jun 1989. doi: 10.1103/PhysRevLett.62.2747. URL http://link.aps.org/doi/10.1103/ PhysRevLett.62.2747. [88] Quantal phase factors accompanying adiabatic changes. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 392(1802):45–57, 1984. ISSN 0080-4630. doi: 10.1098/rspa.1984.0023. 131 [89] Anna Przysiężna, Omjyoti Dutta, and Jakub Zakrzewski. Rice–mele model with topological solitons in an optical lattice. New Journal of Physics, 17(1):013018, 2015. URL http://stacks.iop.org/1367-2630/17/i=1/a=013018. [90] N. Goldman, J. C. Budich, and P. Zoller. Topological quantum matter with ultracold gases in optical lattices. Nat Phys, 12(7):639–645, Jul 2016. ISSN 1745-2473. URL http://dx.doi.org/10.1038/nphys3803. Progress Article. [91] G. Floquet. Sur les équations différentielles linéaires à coefficients périodiques. Annales scientifiques de l’École Normale Supérieure, 12:47–88, 1883. URL http://eudml.org/ doc/80895. [92] J. K. Freericks, H. R. Krishnamurthy, and Th. Pruschke. Theoretical description of time-resolved photoemission spectroscopy: Application to pump-probe experiments. Phys. Rev. Lett., 102:136401, Mar 2009. doi: 10.1103/PhysRevLett.102.136401. URL http://link.aps.org/doi/10.1103/PhysRevLett.102.136401. [93] A. J. Heeger, S. Kivelson, J. R. Schrieffer, and W. P. Su. Solitons in conducting polymers. Rev. Mod. Phys., 60:781–850, Jul 1988. doi: 10.1103/RevModPhys.60.781. URL https://link.aps.org/doi/10.1103/RevModPhys.60.781. [94] A. Gold. Disorder effects on the transport properties in the fractional quantized hall regime. EPL (Europhysics Letters), 1(5):241, 1986. URL http://stacks.iop.org/ 0295-5075/1/i=5/a=007. [95] J. T. Chalker, T. S. Pickles, and Pragya Shukla. Anderson localization in tight-binding models with flat bands. Phys. Rev. B, 82:104209, Sep 2010. doi: 10.1103/PhysRevB. 82.104209. URL https://link.aps.org/doi/10.1103/PhysRevB.82.104209. [96] Jasper van Wezel, Paul Nahai-Williamson, and Siddarth S. Saxena. Exciton-phonondriven charge density wave in tise2 . Phys. Rev. B, 81:165109, Apr 2010. doi: 10.1103/ PhysRevB.81.165109. URL https://link.aps.org/doi/10.1103/PhysRevB.81.165109. [97] Jasper van Wezel, Paul Nahai-Williamson, and Siddarth S. Saxena. Exciton-phonon interactions and superconductivity bordering charge order in tise2 . Phys. Rev. B, 83: 024502, Jan 2011. doi: 10.1103/PhysRevB.83.024502. URL https://link.aps.org/doi/ 10.1103/PhysRevB.83.024502. [98] G. Q. Huang. Surface lattice vibration and electron-phonon interaction in topological insulator bi 2 te 3 (111) films from first principles. EPL (Europhysics Letters), 100(1): 17001, 2012. URL http://stacks.iop.org/0295-5075/100/i=1/a=17001. [99] Sébastien Giraud, Arijit Kundu, and Reinhold Egger. Electron-phonon scattering in topological insulator thin films. Phys. Rev. B, 85:035441, Jan 2012. doi: 10.1103/ PhysRevB.85.035441. URL https://link.aps.org/doi/10.1103/PhysRevB.85.035441. [100] Jize Zhao, Kazuo Ueda, and Xiaoqun Wang. Insulating charge density wave for a halffilled su(n) hubbard model with an attractive on-site interaction in one dimension. Journal of the Physical Society of Japan, 76(11):114711, 2007. doi: 10.1143/JPSJ.76. 114711. URL http://dx.doi.org/10.1143/JPSJ.76.114711. 132 [101] R. Roldán, M. P. López-Sancho, and F. Guinea. Effect of electron-electron interaction on the fermi surface topology of doped graphene. Phys. Rev. B, 77:115410, Mar 2008. doi: 10.1103/PhysRevB.77.115410. URL https://link.aps.org/doi/10.1103/PhysRevB. 77.115410. [102] Hiroshi Watanabe, Kazuhiro Seki, and Seiji Yunoki. Charge-density wave induced by combined electron-electron and electron-phonon interactions in 1t−tise2 : A variational monte carlo study. Phys. Rev. B, 91:205135, May 2015. doi: 10.1103/PhysRevB.91. 205135. URL https://link.aps.org/doi/10.1103/PhysRevB.91.205135. [103] Amrit De and Alexey A. Kovalev. Control of majorana edge modes by a g-factor engineered nanowire spin transistor. Solid State Communications, 198:66 – 70, 2014. ISSN 0038-1098. doi: http://dx.doi.org/10.1016/j.ssc.2013.08.004. URL http://www. sciencedirect.com/science/article/pii/S0038109813003669. SI: Spin Mechanics. [104] M. T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus. Majorana bound state in a coupled quantum-dot hybrid-nanowire system. Science, 354(6319):1557–1562, 2016. ISSN 00368075. doi: 10.1126/science.aaf3961. URL http://science.sciencemag.org/content/354/ 6319/1557. 133