Interference in Complex Quantum SystemsFrom localization in high-dimensional lattices to surface spin echo with moleculesbyJoshua Tyler CantinB. ASc., University of Waterloo, 2014A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Chemistry)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)September 2019c© Joshua Tyler Cantin, 2019The following individuals certify that they have read, and recommend to the Faculty of Graduateand Postdoctoral Studies for acceptance, the dissertation entitled:Interference in Complex Quantum Systems: From localization in high-dimensional latticesto surface spin echo with moleculessubmitted by Joshua Tyler Cantin in partial fulfillment of the requirements for the degree ofDoctor of Philosophy in ChemistryExamining Committee:Supervisor Roman V. Krems, ChemistrySupervisory Committee Member Grenfell N. Patey, ChemistryUniversity Examiner Takamasa Momose, Chemistry and PhysicsUniversity Examiner Marcel Franz, PhysicsAdditional Supervisory Committee Members:Supervisory Committee Member Mona Berciu, PhysicsSupervisory Committee Member Allan K. Bertram, ChemistryiiAbstractThis thesis considers (1) novel manifestations and applications of quantum interference in complex systemsand (2) development of new approaches to study complex quantum systems. First, I examine a generalmodel of particles with long-range hopping amplitudes. For at least 30 years, it has been widely acceptedthat these particles do not undergo Anderson localization in 3D lattices. We show that these particles doundergo Anderson localization in 3D lattices if their hopping amplitudes are isotropic. In contrast, particleswith anisotropic long-range hopping appear to follow the widely held expectations. We show these resultsby demonstrating that cooperative shielding extends to 3D cubic lattices with isotropic long-range hopping,but not with anisotropic long-range hopping and by computing the scaling behaviour of participation ratiosand energy level statistics.Secondly, I develop a fully quantum mechanical model of molecular surface spin-echo experiments,which study surface properties and dynamics by scattering molecules off the sample surface. This model,based on the transfer matrix method, incorporates molecular hyperfine degrees of freedom, allows for theefficient calculation of the experimental signal given a molecule-surface scattering matrix, and permits us tobegin addressing the inverse scattering problem. This fully quantum model is required to properly describethese experiments as the semi-classical methods used to describe experiments using helium-3 atoms donot take magnetic-field induced momentum changes into account. We apply our method to the case ofortho-hydrogen and then apply Bayesian optimization to determine the molecule-surface scattering matrixelements from a calculated signal, for a scattering matrix defined by three parameters. Our work sets thestage for Bayesian optimization to solve the inverse scattering problem for these experiments.Finally, I propose using Bayesian model calibration to improve the convergence of Monte Carlo cal-culations in regions where the sign problem or critical slowing down are an issue. Specifically, Bayesianiiimodel calibration would correct poorly converged Monte Carlo calculations with the information from asmall number of well-converged Monte Carlo calculations. As a simple proof of principle demonstration,we apply Bayesian model calibration to a diagrammatic Monte Carlo calculation of the scattering length ofa spherical potential barrier.ivLay SummaryFor at least 30 years, it has been accepted that particles with long-range hopping are not localized by dis-order in 3D lattices. However, we show that these particles can be localized by disorder if their hoppingamplitudes are isotropic, while we find particles with long-range anisotropic hopping amplitudes followthese prior expectations. We also develop a fully quantum model of molecular hyperfine interferometryexperiments. These experiments study surface properties and dynamics by scattering molecules from sur-faces. Our model allows us to interpret these experiments by giving us an efficient tool needed to addressthe inverse scattering problem. We use Bayesian optimization to begin determining the molecule-surfacescattering matrix elements from a calculated signal. Finally, we propose using Bayesian model calibration toimprove the convergence of Monte Carlo calculations in regions where the sign problem or critical slowingdown cause these calculations to be impractical.vPrefaceChapter 2 is based on the article by Joshua T. Cantin, Tianrui Xu and Roman V. Krems, Effect of theanisotropy of long-range hopping on localization in three-dimensional lattices, Physical Review B, 98,014204 (2018). Specifically, the introduction of the article was expanded with background information tomake sections Sections 2.2 – 2.4, Section 2.5 contains the body of the article, and the conclusion of thearticle was expanded into Sections 2.6 and 2.7. Also, a modified form of the abstract of the article was usedin the thesis abstract and in Chapter 1. The original draft of the article was written by me and edited byRoman V. Krems, Tianrui Xu, and me. I wrote the data analysis code for this work as well as the wrapperfor the diagonalization algorithm. I also wrote a Python version of the original C++ code to check theoriginal code. I performed most of the data analysis. The project was identified and designed by Roman V.Krems, Tianrui Xu, and me. Copyrights to all sections of this chapter taken from the article are owned bythe American Physical Society, 2018.Chapter 3 is based on the article by Joshua T. Cantin, Gil Alexandrowicz, and Roman V. Krems, Transfermatrix theory of surface spin echo experiments with molecules, arXiv:1906.04846 (2019), which is to besubmitted to a peer-reviewed journal. Specifically, the introduction of the article was modified to formSections 3.1 and 3.2, Sections 3.3-3.8 contain the body of the article, and the conclusion of the article wasmodified to form Section 3.9. Also, a modified form of the abstract of the article was used in the thesisabstract and in Chapter 1. The original draft of the article was written by me and edited by Roman V.Krems, Gil Alexandrowicz, and me. I derived the formulae, with some conceptual aid by both Roman andGil, and wrote the corresponding code in Python and Octave. I also performed the numerical calculationsfor the article and for the comparison between the semi-classical and quantum methods. The semi-classicalcode was provided to me by Helen J. Chadwick and was originally written by Oded Godsi. I made somevimodifications to the semi-classical code to allow for a useful comparison. I performed most of the dataanalysis for this article. The project was identified and designed by Roman V. Krems, Gil Alexandrowicz,and me.Chapter 3 also frequently references the article by Oded Godsi, Gefen Corem, Yosef Alkoby, Joshua T.Cantin, Roman V. Krems, Mark F. Somers, Jo¨rg Meyer, Geert-Jan Kroes, Tsofar Maniv and Gil Alexandrow-icz, A general method for controlling and resolving rotational orientation of molecules in molecule-surfacecollisions, Nature Communications 8, 15357 (2017). I developed and performed the transfer matrix analysisand the corresponding numerical calculations for this article. Roman V. Krems contributed to the devel-opment of this transfer matrix analysis. I also wrote the section titled Time independent approach of thisarticle’s Supplementary Information, which Roman edited. My contribution to this article was identifiedand designed by Roman V. Krems, Gil Alexandrowicz, and me.Chapter 4 is original, unpublished work. The project was identified and designed by Roman V. Kremsand me. I wrote the Python code, in consultation of the code provided to me by Rodrigo A. Vargas-Herna´ndez, and performed the numerical calculations and data analysis.Chapter 5 is original, unpublished work. The project was identified and designed by Roman V. Kremsand me. I wrote the Python code for the Bayesian model calibration and performed the numerical calcula-tions and data analysis. The diagrammatic Monte Carlo Fortran code used is freely available on the websiteof Nikolay Prokof’ev [1].viiContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viContents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Introduction and Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.1 Anderson Localization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1.2 Surface Spin Echo Experiments with Molecules . . . . . . . . . . . . . . . . . . . . 21.1.3 Monte Carlo Enhanced with Bayesian Model Calibration . . . . . . . . . . . . . . . 41.2 Main Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Anderson localization for particles with long-range hopping? . . . . . . . . . . . . . . . . . . 62.1 Short- vs. Long-range Hopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6viii2.2 Anderson Localization with Short-range Hopping . . . . . . . . . . . . . . . . . . . . . . . 82.3 Anderson Localization with Long-range Hopping . . . . . . . . . . . . . . . . . . . . . . . 132.4 Cooperative Shielding in One Dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . 172.5 Anderson Localization with Anisotropic Hopping . . . . . . . . . . . . . . . . . . . . . . . 182.5.1 Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 192.5.2 Cooperative Shielding in Three Dimensions . . . . . . . . . . . . . . . . . . . . . . 212.5.3 Participation Ratios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.5.4 Energy Level Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 382.7 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.8 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413 Molecular hyperfine interferometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 423.2 From Helium-3 Spin Echo to Molecular Hyperfine Interferometry . . . . . . . . . . . . . . 433.3 Description of a Molecular Hyperfine Interferometry Experiment . . . . . . . . . . . . . . . 443.3.1 Molecular Hyperfine Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . 473.4 Impact of the Magnetic Lens on the Molecular States . . . . . . . . . . . . . . . . . . . . . 503.5 Wavepacket Propagation and Signal Calculation . . . . . . . . . . . . . . . . . . . . . . . . 523.6 Transfer Matrix Formalism with Internal Degrees of Freedom . . . . . . . . . . . . . . . . . 553.6.1 Propagation and Discontinuity Matrices . . . . . . . . . . . . . . . . . . . . . . . . 563.6.2 Rotation Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.6.3 Scattering Transfer Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 623.6.4 Calculation of Eigenstate Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . 643.7 Application to ortho-Hydrogen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 653.7.1 Ortho-Hydrogen Hyperfine Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . 653.7.2 Experiment and Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 663.8 Comparison with a Semi-Classical Method . . . . . . . . . . . . . . . . . . . . . . . . . . 703.9 Conclusion and Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75ix3.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 774 Bayesian optimization to determine the molecule-surface scattering matrix . . . . . . . . . . 784.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.2 The Inverse Scattering Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 784.3 Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 794.3.1 Gaussian Processes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804.3.2 Gaussian Process Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 834.3.3 Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 854.4 Preliminary Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 884.5 Conclusion and Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 934.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 955 Monte Carlo enhanced with Bayesian model calibration . . . . . . . . . . . . . . . . . . . . . 965.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.2 Monte Carlo Calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 965.3 Bayesian Model Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 985.4 Simplified Bayesian Model Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 995.5 Proof of Principle Demonstration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1015.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1035.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1056 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1066.1 Summary of Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1066.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1076.2.1 Anderson Localization with Long-range Hopping . . . . . . . . . . . . . . . . . . . 1086.2.2 Molecular Hyperfine Interferometry . . . . . . . . . . . . . . . . . . . . . . . . . . 1086.2.3 Monte Carlo Enhanced with Bayesian Model Calibration . . . . . . . . . . . . . . . 109Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110xA Supplementary material for Anderson localization for particles with long-range hopping . . 124A.1 Outline of the Derivation for Equation 2.26 . . . . . . . . . . . . . . . . . . . . . . . . . . 124B Supplementary material for molecular hyperfine interferometry . . . . . . . . . . . . . . . . 129B.1 Schro¨dinger Equation for Eigenstate Coefficients . . . . . . . . . . . . . . . . . . . . . . . 129B.2 Coefficient Relations Across a Discontinuity . . . . . . . . . . . . . . . . . . . . . . . . . . 130B.3 Computational Parameters Used for the Application to ortho-Hydrogen . . . . . . . . . . . 131xiList of TablesTable B.1 Relative Probabilities of the State Selector and the Detector . . . . . . . . . . . . . . . . 131xiiList of FiguresFigure 2.1 Localized Wavefunction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10Figure 2.2 Energy Gap for 1D, 3D Isotropic and 3D Anisotropic Hopping Systems . . . . . . . . . 22Figure 2.3 Energy Level Structures in the Case of Isotropic and Anisotropic Hopping . . . . . . . . 26Figure 2.4 Scaling of Participation Ratios for a 3D System with Isotropic Hopping and α = 0 . . . 28Figure 2.5 Scaling of Participation Ratios for a 3D System with Anisotropic Hopping and α = 0 . . 31Figure 2.6 Scaling of Participation Ratios for a 3D System with α = 1 for Isotropic and AnisotropicHopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 2.7 Scaling of Participation Ratios for a 3D System with α = 3 for Isotropic and AnisotropicHopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35Figure 2.8 Scaling of Mean Energy Level Spacing Ratios for α = 3 with Isotropic and AnisotropicHopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36Figure 2.9 Scaling of Mean Energy Level Spacing Ratios for Isotropic Hopping with α = 1,3 NearZero Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 2.10 Mean Energy Level Spacing Ratios for α = 3 as a Function of Filling Fraction andDisorder Strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38Figure 3.1 Schematic of a Generic Molecular Hyperfine Interferometer . . . . . . . . . . . . . . . 48Figure 3.2 Abstraction of a Generic Field Profile of a Molecular Hyperfine Interferometer . . . . . 56Figure 3.3 Approximate Magnetic Field Profile for an Experiment Using oH2. . . . . . . . . . . . 66Figure 3.4 Theoretical and Experimental Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . 67Figure 3.5 Theoretical and Experimental 2D Fourier Amplitude Plots . . . . . . . . . . . . . . . . 68xiiiFigure 3.6 Calculated Signals with Various Scattering Matrices . . . . . . . . . . . . . . . . . . . 71Figure 3.7 Fourier Transform of Calculated Signals with Various Scattering Matrices . . . . . . . . 72Figure 3.8 Error between Semi-classical and Fully Quantum Methods . . . . . . . . . . . . . . . . 74Figure 3.9 Comparison of Fourier Amplitudes Between Semi-classical and Fully Quantum Methods 75Figure 4.1 Illustration of a Gaussian Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81Figure 4.2 Representation of a Gaussian Process . . . . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.3 Drawn Functions for Various Kernels . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 4.4 Fitting a Gaussian Process Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 4.5 Example of Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87Figure 4.6 Gaussian Process Model of the Mean Square Error for Three Unknown Phases . . . . . 91Figure 4.7 Convergence of the Root Mean Square Error and the Determined Phases . . . . . . . . . 92Figure 5.1 Monte Carlo Enhanced with Bayesian Model Calibration Proof of Principle . . . . . . . 103xivAcknowledgementsChrist! Tianrui! Mama! Papa! Roman! Krems group! Gil! Alexandrowicz group! Takuto! Gord! Ute!Birgitta! Whaley group! Norm! Geoff! Michelle! Claude!Je veux remercie mon pe`re Joseph Re´jean Rene´ Cantin, ma gand-me`re Marie Marguerite Cantin (ne´e Pom-inville), et mon grand-pe`re Joseph Ernest Cantin.Auch, wu¨rde Ich danken meine Mutter Sylvia Anna Cantin (geb. Morell), meine Großmutter ChristineMorell (geb. Weirauch), und mein Großvater Anton Morell.我也想给我妻子徐天睿说谢谢。她支持我支持了几年了。她也鼓励我鼓励了几年了。I would also like to thank the Natural Sciences and Engineering Research Council of Canada, the Horizon2020 Research and Innovation Programme, and the International Research Training Group (IRTG 2079)“Cold Controlled Ensembles in Physics and Chemistry” for their financial support.Thank you also to the members of IRTG 2079, the members of the Canadian Center for Research on Ultra-Cold Systems, and to everyone else with whom I interacted during my PhD, the number of whom is alwaystoo large to list.S.D.G.xvIn memory of Ernest Cantin and Tony MorellxviChapter 1Introduction1.1 Introduction and OverviewOver the past several decades, researchers have achieved an unprecedented level of control of atomic andmolecular single and many-body systems [2]. Researchers now routinely capture ensembles of atoms inoptical dipole traps [3–5], magneto-optical traps [5, 6], and optical lattices [5, 7–9] at milli- to microkelvintemperatures [5, 10]. Ensembles of molecules can also be placed in various traps [11–14], though withgreater difficulty. The large de Broglie wavelengths of these atoms and molecules gives rise to interestingquantum interference phenomena, from the formation of atomic [9, 15, 16] and molecular [17, 18] Bose-Einstein condensates to Anderson localization [19, 20] (the trapping of a quantum particle by disorder).1.1.1 Anderson LocalizationThese new experiments with ensembles of atoms and molecules are stimulating the research field to revisitmodels that were not thoroughly studied. For example, recent experiments have shown the ability to tunethe range of interactions in a one dimensional lattice [21, 22] and there is an experiment that can tune theanisotropy of the interactions of atoms dressed by Rydberg excitations [23]. Thus, one model with renewedinterest describes the motion of particles on a disordered lattice that permits long-range hopping amplitudeswith various forms of anisotropy. Interestingly, it has been widely accepted for at least 30 years that particleswith such long-range hopping do not undergo Anderson localization in three dimensional lattices [24–26],despite a subtle hint otherwise [27]. However, several recent studies demonstrated localization of particles1with long-range hopping in one dimensional lattices [28–30]. In particular, it was recently shown that theeffect of long-range hopping in one dimensional lattices can be mitigated by cooperative shielding, whichmakes the system behave effectively as one with short-range hopping [30].In Chapter 2, we demonstrate that particles with long-range isotropic hopping amplitudes can undergoAnderson localization in three dimensional lattices and that whether the particle localizes depends on theanisotropy of the hopping amplitudes. We show that cooperative shielding extends to three dimensionallattices with isotropic long-range r−α hopping, but not to three dimensional cubic lattices with anisotropiclong-range hopping. We demonstrate the presence of localization in three dimensional lattices with uniform(α = 0) isotropic long-range hopping and the absence of localization with uniform anisotropic long-rangehopping by using the scaling behaviour of eigenstate participation ratios. We compare the scaling behaviourof participation ratios and energy level statistics to show that the existence of delocalized, non-ergodicextended, or localized states in the presence of disorder depends on both the exponents α and the anisotropyof the long-range hopping amplitudes.1.1.2 Surface Spin Echo Experiments with MoleculesGreater control over cold atomic and molecular beams allows researchers to use these beams as tools to studyproblems in other areas of physics. Researchers now use these beams to investigate systems as different asthe electric dipole moment of the electron [31] and the dynamics of material surfaces [32]. For instance, 3Hespin-echo experiments use beams of 3He atoms to study surface morphology, molecular and atomic surfacediffusion [33–39], phonon dispersions [40, 41], phason dispersions [42] and even phase transitions of ionicliquids [43]. However, the interactions between 3He atoms and surfaces or their adsorbates are typicallyisotropic and weak. To overcome these limitations, we propose to use generic closed shell molecules insteadof 3He in surface spin-echo experiments. The molecular degrees of freedom, such as rotation, may beexploited to provide additional insight into surfaces and the behaviour of their adsorbates. Indeed, a recentexperiment has shown that ortho-hydrogen can be used as a probe that is sensitive to the orientation ofa Cu(115) surface [44]. However, the additional degrees of freedom offered by molecules also pose atheoretical challenge: a large manifold of molecular states and magnetic field-induced couplings betweeninternal states.To address this challenge, in Chapter 3, we develop a fully quantum mechanical approach to model2these molecular surface spin-echo experiments, which are also known as molecular hyperfine interferometryexperiments. We thereby connect the experimental signal to the elements of the time-independent molecule-surface scattering matrix. We derive a one-dimensional transfer matrix method that includes the molecularhyperfine degrees of freedom and accounts for the spatial separation of the molecular wavepackets due tothe magnetic control fields. This work gives us an efficient tool that we need to address the inverse scatteringproblem for these experiments and sets the stage for a framework that describes molecular surface spin-echoexperiments to study dynamic surfaces. We apply the fully quantum method to the case of ortho-hydrogen,show that the calculated experimental signal is sensitive to the scattering matrix elements, and perform apreliminary comparison to experiment.With this new framework, we also set the stage for machine learning techniques to determine the scat-tering matrix elements from experimental measurements. These machine learning techniques are gainingincreasing popularity not only in the technology industry, but also in physics and chemistry [45–50]. Neu-ral networks have been used in potential energy surface fitting [51], to identify phase transitions [52, 53],and for quantum state tomography [54]. Support vector machines have been used to identify order param-eters [55] and Gaussian processes have been used in potential energy surface fitting [56–58], prediction ofa phase transition [59], and scattering calculations [60–62]. Bayesian optimization has been used to tunedensity functionals [63] and for the inverse scattering problem in quantum reaction dynamics [64].In Chapter 4, we apply Bayesian optimization to the inverse scattering problem for molecular hyperfineinterferometry experiments. The goal of this inverse scattering problem is to determine the several dozenscattering matrix elements from the experimental signal. Such a large parameter space makes it difficult touse traditional optimization methods. Our approach to this problem is to efficiently minimize the differencebetween the experimental signal and the results of the transfer-matrix computations by using the Bayesianoptimization algorithm to intelligently vary the scattering matrix elements. In this way, we apply Bayesianoptimization to determine the scattering matrix corresponding to a target calculated signal for a diagonalmatrix determined by three or fewer parameters. We also demonstrate that Bayesian optimization is sig-nificantly more efficient than grid searches of comparable resolution and lay the foundation for applyingBayesian optimization to determine scattering matrices defined by several tens of parameters.31.1.3 Monte Carlo Enhanced with Bayesian Model CalibrationGiven the complexity of new cold atomic and molecular experiments, such as those studying the two di-mensional Fermi-Hubbard model [65, 66], it is important to have a solid theoretical understanding of thephysics involved. As such, there are several many-body numerical techniques that are used to study thesesystems. These include the density matrix renormalization group [67–69], dynamical mean field theory [70],path integral Monte Carlo [71, 72], and diagrammatic Monte Carlo [73–75]. However, there are still manysystems, such as the two dimensional Hubbard model [76] just mentioned, that are not easily addressed bythese numerical techniques. Because of these difficulties, researchers are beginning to use machine learningtechniques to improve the efficiency of various numerical methods. For example, restricted Boltzmann dis-tributions have been used to learn effective Hamiltonians in classical Monte Carlo calculations [77, 78] andCui and Krems have used Bayesian model calibration to improve the efficiency of scattering calculations bycombining quantum scattering calculations with classical scattering calculations [60].In Chapter 5, we propose the idea to use Bayesian model calibration to improve the efficiency of MonteCarlo calculations in regions where critical slowing down or the sign problem inhibit the efficient calculationof observables. Particularly, we propose to perform several poorly converged calculations in these regionsand then use Bayesian model calibration to correct these calculations with the information obtained fromjust a few well-converged calculations. Such a procedure would reduce the computational cost required toachieve the desired level of accuracy. As a simple proof of principle demonstration, we apply Bayesianmodel calibration to a simple diagrammatic Monte Carlo calculation of the scattering length of a sphericalpotential barrier.In Chapter 6 of this thesis, I summarize the results of these different projects related to complex quantumsystems and discuss their future directions.1.2 Main ResultsIn this section, I list the main results of this thesis.• Showed that particles with long-range isotropic hopping can localize in three dimensional latticesand that, generally, the presence or absence of delocalized, extended non-ergodic, or localized statesdepends on both the range and the anisotropy of a particle’s hopping amplitudes (Chapter 2)4• Showed that cooperative shielding extends to three-dimensional systems, if the hopping amplitudesare isotropic, while cooperative shielding appears absent if the hopping is anisotropic with the angularform of the interaction between dipoles aligned along one of the principal axes of the cubic lattice(Chapter 2)• Proposed using generic closed shell molecules for molecular hyperfine interferometry experimentsand developed an efficient fully quantum formalism needed to address the inverse scattering problem(Chapter 3)• Applied the fully quantum formalism to the description of molecular hyperfine interferometry ex-periments using oH2 molecules and showed that various features of the scattering matrix affect theexperimental signal (Chapter 3)• Applied Bayesian optimization to the inverse problem for molecular hyperfine interferometry exper-iments and determined a scattering matrix parametrized by three phases, from data generated fromknown target parameters (Chapter 4)• Laid the foundation to efficiently determine, via Bayesian optimization, matrices defined by morethan three parameters (Chapter 4)• Proposed applying Bayesian model calibration to improve the efficiency of Monte Carlo calculationsin regions where the sign problem or critical slowing down are problematic (Chapter 5)• Provided a proof of principle demonstration that Bayesian model calibration can improve the ac-curacy of many poorly converged Monte Carlo calculations with the information from just a fewwell-converged calculations (Chapter 5)5Chapter 2Anderson localization for particles withlong-range hopping?I begin this chapter by precisely defining the terms long- and short-range hopping in Section 2.1. I thendiscuss the original insights into Anderson localization in Section 2.2. Next, I describe one of the primaryinfluential ideas related to Anderson localization of particles with long-range hopping in Section 2.3. Idiscuss the recently discovered phenomenon of cooperative shielding in Section 2.4. In Section 2.5, I discussour results concerning the effect of anisotropy of long-range hopping amplitudes on Anderson localizationand the connection with cooperative shielding. I then conclude the chapter and discuss future directions inSections 2.6 and 2.7, respectively. Finally, I list the key results in Section 2.8.2.1 Short- vs. Long-range HoppingBefore discussing Anderson localization with short-range hopping, I will first define short- and long-rangehopping. A common model of electrons moving in a material is the tight-binding Hamiltonian with nearestneighbour hopping [79]. In this model, each atom is modelled as a lattice site, between two of whichthe electron can “hop”. This hopping is parametrized by a hopping amplitude that describes the probabilityamplitude of the electron (or other particle or quasiparticle) to move from one lattice site to another (i.e. fromatom to atom). Nearest neighbour hopping describes the situation where the hopping amplitude betweentwo sites is a constant for neighbouring lattice sites, and zero otherwise; this is the shortest possible hopping6range (excepting a range of zero, that is).There can also be next-nearest neighbour hopping where the particle can hop two sites simultaneously.While longer-ranged than nearest-neighbour hopping, I also consider next-nearest neighbour hopping short-ranged for the purposes of this thesis. Extending this notion further, a particle could in principle move anynumber of sites in a single hop. Thus, the definition of “short-” or “long-” ranged hopping is, in a sense,arbitrary. Is a hop of 10 sites long? 100 sites?Given these difficulties, we could instead look at the decay profile of the likelihood the particle hops aparticular distance. That is, how does the hopping amplitude t(~r) decay as a function of |~r|, where~r is thevector between the initial and final sites? Clearly, if t(~r) does not decay as a function of distance, such thatthe likelihood to hop to the next site or to a site infinitely far away are equal, the hopping can be consideredlong-ranged. Again, however, we are left with arbitrarily defining a cut-off between short- and long-rangedhopping. A common, but by no means universal, definition is to choose a cut-off in line with whether theintegral over all space of the hopping amplitude diverges or converges [24, 26, 28, 80].We examine the integralI =∫ddr t(~r), (2.1)where d is the dimension of the lattice. We then declare the hopping to be short-ranged if the integralconverges and long-ranged if the integral diverges. Such a definition aligns well with the notion of long- andshort-ranged forces [80, 81].For an exponentially decaying t(~r), I is clearly finite. However, if t(~r) follows the power-law profile|~r|−α , the convergence of I depends on the value of α . In particular, if α > d, the integral converges and thehopping is short-ranged. If α < d, the integral diverges and the hopping is long-ranged. Finally, if α = d,the integral also diverges and the hopping is considered critically long-ranged. Taking this definition oflong- and short-ranged hopping, I will now explain the Anderson localization of particles with short-rangedhopping.72.2 Anderson Localization with Short-range HoppingDisorder is present throughout materials, from dislocations in crystal lattices to amorphous mixes of organicmolecules; it is critical that we understand the impact of this disorder on material physics and energy trans-port through these materials. In 1958, Anderson [24] was the first to show that a quantum particle remainslocalized if there is sufficient disorder present within the system. That is, a quantum particle placed at aparticular location within a material will eventually diffuse away, unless sufficient disorder exists to trap itin or near its original location. In particular, the disorder scatters the wave-like particle to such an extentthat the resultant interference pattern is constructive primarily near the particle’s initial location.It was later realized that this localization phenomenon is not restricted to quantum particles, but to awhole array of classical wave phenomena. Experiments have shown the presence of Anderson localizationin light [82, 83] (even with overhead transparencies! [84]) and sound [85–87], in addition to ultra-cold atoms[19, 20, 88]. Indeed, an optical fibre has been developed that does not use total internal reflection to trap thelight within the fibre, but rather localization in the radial direction caused by engineered disorder in the fibrematerial [89].Let us look in more detail at a particle, such as an electron, a quasi-particle, or an exciton, movingthrough a material. First, we model the system as a generic particle moving on a lattice. Such a modelapplies to electrons moving from atom to atom in a crystal, rotational excitations moving from molecule tomolecule in an optical lattice, or to electronic excitations moving from molecule to molecule in an organic-semiconductor. For now, we focus on a particle that only has nearest neighbour hopping in the tight-bindingmodel:H = ∑〈i, j〉ti jc†jci+∑iεic†i ci, (2.2)where ti j is the hopping amplitude of a particle to move from lattice site i to site j, c†i is the creation operatorfor a particle on site i and εi is the energy of a particle on site i. The sums run over all sites of a d-dimensionallattice, while the angle brackets indicate the inclusion of only nearest neighbours in the sum.If the on-site energies εi and the hopping amplitudes are all identical, H can be diagonalized in themomentum basis and a particle placed on a single lattice site will ballistically transport away from its orig-8inal site. However, we can introduce disorder into the model by randomizing either εi, known as on-site ordiagonal disorder, or by randomizing ti j, known as off-diagonal disorder. Many different probability distri-butions can be used for this, from uniform, to Gaussian, to binary. Typically, the “strength” of the disorderis characterized by some measure of the width of the distribution; that is, the strength is characterized byhow much variation is likely from site to site or hop to hop.Typically, as disorder is added to the system, the transport changes from ballistic to diffusive. Thediffusion coefficient in this regime is a function of the disorder strength. It is initially just renormalized toa smaller value (weak localization), but eventually becomes zero at a particular disorder strength (stronglocalization). Note, however, that not all systems have this diffusive regime. Particularly, while a threedimensional system exhibits the transition from ballistic to diffusive motion and then from diffusive motionto localization, one and two dimensional systems can be localized by any small amount of disorder [90].Once the system is in the strong localization regime, the particle can no longer traverse the lattice but retainsa profile that decays exponentially from the initial site. That is,∣∣ψ (x)∣∣2 ∼ e− x−x0ξ for large x, (2.3)where ψ (x) is the particle wavefunction as a function of position x, x0 is the original location of the particle,and ξ is termed the localization length [91].There are two things to note. First, though the wavefunction decays exponentially at long distances, itdoes not do so smoothly but varies widely within an exponential envelope, see Figure 2.1. The wavefunctionis also not truly exponentially decaying near x0. Secondly, the localization length is rather loosely defined;differences of factors of two are common in the definition of ξ and differences arise in practical calculationas the definition of “large” x is rather arbitrary. Thus, the localization length is better viewed as a conceptand a referential length scale than an actual material property.At this point, I will discuss the importance of quantum coherence in Anderson localization and demon-strate one of the many subtleties of Anderson localization. Following the Les Houches lecture notes byMu¨ller and Delande [91], let us consider a single particle in a continuous one dimensional system that hasno potential except for delta-functions spaced apart at random intervals (i.e. a Dirac comb with randomintervals instead of equal intervals).90 50 100 150Site0.000.050.100.150.200.25||20 50 100 150Site10 610 410 2||2Figure 2.1: Probability distribution |ψ|2 of a particle initially placed on the centre site of a disorderedone dimensional lattice. The Hamiltonian corresponds to Eqn. (2.2) with the hopping amplitudesconstant and the on-site energies drawn from a uniform distribution of width 2.5t. Both panelsshow the probability distribution at time T = 50t−1 and only differ in the scaling of the ordinate.Exponential, or near exponential, decay of the wave function is observed at large distances.If we look at two neighbouring delta-functions, spaced apart by a distance ∆z, the transmission T12 of aquantum particle from the left of the first delta-function to the right of the second is given by [91]:T12 =T1T2∣∣1−√R1R2eik∆z∣∣2 , (2.4)where Ti is the transmission through the ith delta-function, Ri is the reflection from the ith delta-function,and k is the linear momentum of the particle.If we decide to ignore any coherence between pairs of delta-functions, we can treat T12 as the represen-tative transmission of all the pairs and average the phase θ ≡ k∆z over the interval [0,2pi) to get an ensembleaverage [91]:〈T12〉= T1T21−R1R2 , (2.5)which is the classical solution. From this, we can get the resistance to be proportional to the system lengthand obtain a result consistent with Ohm’s law [91]. Clearly, ignoring the coherence between pairs of scat-terers does not result in an answer consistent with Anderson localization (cf. Eqn. 2.3). This highlights the10necessity of coherence for Anderson localization; the phenomenon of localization is weakened or destroyedat high temperatures or by the presence of other degrees of freedom (such as phonons) that couple to theparticle in question [91].Thus, to obtain a result consistent with Eqn. 2.3, we would seemingly need to keep track of all of thephase relationships between scattering events and then perform a suitable averaging over all phases. This isno easy task. Intriguingly, there is a simpler route. Let us take the ensemble average of lnT12 [91]:lnT12 = lnT1+ lnT2+ ln∣∣∣1−√R1R2eik∆z∣∣∣2 (2.6)〈lnT12〉= lnT1+ lnT2+〈ln∣∣∣1−√R1R2eik∆z∣∣∣2〉= lnT1+ lnT2, (2.7)where the third term on the right hand side of the second line is zero as the integrand is analytic inside theintegration contour (i.e. a circle of radius less than 1 and centred at 1).From the above, we can infer that the total log-averaged transmission lnTtot through a set of N delta-functions is given by [91]:〈lnTtot〉=N∑ilnTi= N lnT1= nL lnT1, (2.8)where it is assumed that all Ti = T1 and that the delta-functions appear at a linear density n in a system oflength L, such that the number of delta-functions N = nL. Thus, the typical transmission Ttyp through a11region of length L is given by [91]:Ttyp ≡ e〈lnTtot〉= enL lnT1= e−nL|lnT1|≡ e− Lξ , (2.9)where the third lines follows as T1 < 1 and the localization length is identified as ξ−1 ≡ n |lnT1|.This shows that while the transmission itself is a self-averaging quantity for a classical or incoherentquantum system, it is the log-transmission that is self-averaging in the coherent quantum system (or classicalwave system). Indeed, the transmission through a disordered sample varies from sample to sample and if onemeasures the transmission of a set of disordered systems, the resultant probability distribution is log-normal,rather than normal, and the representative or “typical” transmission of the ensemble is the geometric meane〈lnTtot〉 rather than the arithmetic mean 〈Ttot〉 [91].The above analysis of the transmission through a disordered system not only highlights the importanceof coherence for Anderson localization, but also the subtleties that can appear when trying to find anddetermine the representative quantities of a disordered system. These subtleties can hinder the analysisof these disordered systems, but also provide for a rich field that includes many different approaches andtechniques. These techniques range from locator expansions [24], to scaling theory [90], to replica fieldtheory [92].Finally, it is also important to understand how Anderson localization manifests itself in different dimen-sions. In 1979, Abrahams et al. [90] used scaling theory to show that a single particle always localizes inone and two dimensional systems, assuming zero temperature, an infinite system, and no correlation in thedisorder. With correlation in the disorder, however, delocalized states can still exist [93–96]. In contrastto these lower dimensions, Abrahams et al. also showed that in three dimensions there is a critical disorderstrength that separates a diffusive regime from a localized regime [90]. It is also known that this diffusion-localization transition is a critical transition where the localization length ξ and the diffusion coefficientD have the relations ξ ∼ (w−wc)−ν and D ∼ (wc−w)ν with the critical exponent ν = 1.58± 0.01, the12disorder strength w, and the critical disorder strength wc [90, 97, 98].I would like to emphasize that all of the above results are for systems whose particles have short-rangehopping amplitudes or for systems of classical waves. If the hopping amplitudes are long-range, however,many of the above results can change drastically, while others remain the same. In the next section, I willdiscuss some of the known results and how well they actually apply to systems with long-range hopping.2.3 Anderson Localization with Long-range HoppingIn Section 2.2, the presence or absence of Anderson localization in non-interacting systems with nearestneighbour hopping was examined for various dimensionalities. Here, I discuss situations where the hoppingamplitude ti j of Eqn. (2.2) is beyond nearest-neighbour and decays as r−α , where α ∈ R≥0 is the decayexponent or the hopping “range”. To reiterate, I consider the hopping “short-ranged” when α > d, thesystem dimension, “long-ranged” when α ≤ d, and critically long-ranged when α = d. There is increasinginterest in disordered lattice models with long-range hopping, as excitons with hopping range α = 3 arefound in molecular crystals with topological disorder [99], organic semiconductors with impurities [100],diluted ensembles of atoms and molecules trapped in optical lattices [101–106], J-aggregates [107], photo-synthetic complexes [108, 109], and ensembles of Rydberg atoms [110–112].Systems whose hopping amplitudes exhibited such a power-law decaying profile were also already ex-amined by Anderson in his seminal paper [24]. Through use of a locator expansion, Anderson concluded thatlocalization is present when the hopping is short-ranged and absent otherwise. Later Levitov [25, 26, 113]approached the system with a resonance argument and came to the same conclusion. I will describe thisargument in detail, with some slight generalizations.In essence, a particle with long-range hopping can diffuse away from its initial lattice site because ofresonances between sites that are an arbitrary distance away from each other [25]. Such resonances occurwhenever the hopping amplitude ti j between sites i and j is similar or larger in magnitude to the differencein on-site energy between the two sites. That is, a resonance occurs whenti j '∣∣E j−Ei∣∣ , (2.10)where Ei is the energy of the particle at site i [25].13Let’s assume that the particle is initially localized to some single lattice site i. To prove that this particledoes not remain localized, Levitov [25] showed that site i is in resonance with another lattice site j suchthat the distance between the two sites∣∣~r j−~ri∣∣> R, where R is an arbitrarily large distance. We now followthe main ideas of Levitov’s proof, but generalized to a d-dimensional system [25]. We begin by takingconcentric (d−1)-spheres centred on the site i such that [25]2kR <∣∣~r j−~ri∣∣< 2k+1R, (2.11)where k ∈ Z and~r is some distance vector. Now, the volume Vk,d of the kth region isVk,d = vold(2k+1R)−vold(2kR)∝(2k+1R)d−(2kR)d∝(2kR)d, (2.12)where vold (r) gives the volume of a (d−1)-sphere of diameter r. Taking the lattice site density n, thenumber of sites between the concentric spheres is given by nVk,d [25].Now, assuming that the on-site energies are identically and independently drawn from a uniform distri-bution centred about zero and of width W , the probability distribution of the difference in energy betweentwo sites ∆≡ ∣∣E j−Ei∣∣ is triangular and given byP(∆) =2W(1− ∆W). (2.13)Taking the hopping amplitude ti j = 1rαi j , where ri j ≡∣∣~r j−~ri∣∣, the probability that sites i and j are in resonance14is given byPres = P(∆≤ ti j)=∫ ti j0d∆ P(∆)=2ti jW− t2i jW 2(2.14)=2Wrαi j+O(1r2αi j)≈ 2W(2kR)α , (2.15)where ti j = 1rαi j is substituted in the penultimate line and the final line arises as~r j is taken to be in the regionbetween the kth and (k+1)th concentric (d−1)-spheres, see Eqn. (2.11). Now, the number of resonancesNk in the kth shell is given byNk = nVk,dPres=2n(2kR)dW(2kR)α=2nW(2kR)α−d . (2.16)From the above, we can see that three particular cases arise:(i) If α > d, then Nk → 0 as k→ ∞. That is, there are resonances only a finite distance away and thesystem remains localized.(ii) If α < d, then Nk → ∞ as k→ ∞. That is, there are an infinite number of resonances infinitely faraway that cause the system to delocalize.(iii) If α = d, then Nk = const. as k→ ∞. That is, there are still at least some resonances infinitely faraway and the system delocalizes.We note that the neglected O(1r2αi j)term does not change the above conclusions as, when α ≥ d, theneglected term approaches zero as k→ ∞ and, when α < d, the neglected term cannot change the fact that15Nk diverges as k→ ∞.The above analysis thus predicts that systems with short-range hopping localize, while those with long-range hopping delocalize. A similar analysis has also been recently applied to study many-body localizationof particles with long-range interactions [114]. The analysis by both Anderson [24] and Levitov [25, 26,113] have led to it becoming widely accepted that non-interacting particles with long-range hopping donot undergo Anderson localization. However, this is despite the fact that, in the same year as Levitov’swork [25], Burin and Maksimov [27] use a renormalization group analysis to show that one may expecta system with long-range hopping to localize under certain conditions. Burin and Maksimov indicate thattheir differing prediction from Levitov’s analysis arises as they treat a system whose hopping amplitudes areisotropic with respect to angle, while the hopping amplitudes in Levitov’s system have a dipolar-like angulardependence.Interestingly, the analysis I presented above, as distilled from Ref. [25], does not appear to depend onthe angular dependence of the hopping amplitudes. This calls into question why the paper by Levitov [25]and that of Burin and Maksimov [27] come to different conclusions. Unfortunately, this difference has notbeen appreciated by the community. Indeed, the paper by Burin and Maksimov [27] has had only a singlecitation [115] until 2018, when it received two more [116, 117].Furthermore, an increasing number of papers have identified situations where long-range hopping (orinteractions) does not result in delocalization. Deng et al. [28] have seen algebraic localization of a particlewith long-range hopping in a 1D system with off-diagonal disorder. Nandkishore and Sondhi [118] haveshown using field-theory techniques that many-body systems with long-range interactions in one and twodimensions can localize and hypothesized that the same is true in three dimensions. Nandkishore and Sondhieven go so far as to question the validity of the locator expansion and resonance arguments for systems withlong-range hopping.Furthermore, Santos et al. [119] and Celardo et al. [30] have discovered a phenomenon, cooperativeshielding, that causes effective short-range behaviour in a system with long-range interactions or hopping.Cooperative shielding precludes the resonance arguments used to demonstrate delocalization in the presenceof on-site disorder and allows localization to take place. Celardo et al. [30] have in particular demonstratedlocalization in one dimension for particles with long-range hopping. Ossipov has also demonstrated local-16ization of particles with uniform hopping on a d simplex, via a similar mechanism [29]. In the next section,I discuss the phenomenon of cooperative shielding in one dimension and how it impacts localization be-haviour.2.4 Cooperative Shielding in One DimensionCooperative shielding in one dimensional systems was first discovered by Santos et al. [119] and Celardo etal. [30] in 2016. In the non-interacting case, cooperative shielding describes the situation where the long-range terms of the Hamiltonian cause the Hilbert space to split into a long-range and a short-range subspace.The short-range subspace can be described by an effective short-ranged Hamiltonian. In the thermodynamiclimit, if a wavefunction begins in the short-range subspace, it will be described by the effective short-rangedHamiltonian for all time. Adding disorder to this short-ranged effective Hamiltonian can then cause theparticle to localize.Following Ref. [119], one can take the HamiltonianHˆ = Hˆ0+Vˆ, (2.17)which consists of a short-ranged part Hˆ0 and a long-ranged part Vˆ such that[Hˆ0,Vˆ]= 0. If Vˆ has a degen-erate eigensubspace V such that Vˆ |νk〉= ν |νk〉∀|νk〉 ∈ V , then|ψ (t)〉= e−i Hˆth¯ |ψ (0)〉= e−iνth¯ e−iHˆ0th¯ |ψ (0)〉 , (2.18)if |ψ (0)〉 = ∑k ck |νk〉, where ck are the expansion coefficients of the initial wavefunction |ψ (0)〉. In thissituation, where the wavefunction is initially contained within the eigensubspace V , one can see that theonly effect of the long-range part of the Hamiltonian is to add a global phase; the effective Hamiltonian isHˆ0 and the wavefunction is “shielded” from the long-range terms of Hˆ [119].Now, however, if[Hˆ0,Vˆ]6= 0, but is such that the Hilbert subspace V is still quasi-degenerate, theshielding effect still occurs, albeit for finite time [119]. The shielding time is positively related to the gapbetween the Hilbert subspaces and, if the gap grows with the system size, diverges in the limit of infinite17system size. This divergence with the system size results in cooperativity [30, 119].Thus, we see that systems with long-range hopping can, in principle, localize, at least in one-dimension.Indeed, Celardo et al. [30] demonstrate single particle localization for a one dimensional system with long-range hopping. In the next section, I discuss how we have extended the notion of cooperative shielding tothree dimensions and demonstrate that the presence of cooperative shielding is dependent on the angularprofile of the hopping amplitudes.2.5 Anderson Localization with Anisotropic HoppingAs mentioned in prior sections, the effect of anisotropy of long-range hopping on localization in high-dimensional lattices has not been examined. We examine this here by studying the impact of anisotropy oncooperative shielding in 3D lattices. Specifically, we show that the cooperative shielding demonstrated pre-viously for 1D lattices extends to particles with isotropic long-range hopping in 3D cubic lattices. However,this cooperative shielding does not extend to particles in 3D cubic lattices with the anisotropic hoppingthat corresponds to the dipole-dipole interaction between dipoles aligned along one of the lattice axes. Wealso provide numerical and analytical evidence for the localization of particles with long-range isotropichopping in 3D lattices.The remainder of this chapter is organized as follows: Following the description of the models, we dis-cuss the phenomenon of cooperative shielding in 3D systems. We analytically diagonalize the Hamiltonianwith isotropic hopping t ∝ 1/rα for 3D cubic lattices for arbitrary α and illustrate that it exhibits the sameenergy level structure as that of a 1D system with cooperative shielding. We find this to be the case forboth periodic and open boundary conditions. The presence of cooperative shielding suggests the presenceof localization. Contrastingly, we find no clear evidence for cooperative shielding at most values of α foranisotropic dipolar hopping, for either periodic or open boundary conditions.We use the scaling behaviour of participation ratios to demonstrate the existence of localized statesfor isotropic hopping when α = 0 (i.e. uniform infinite-range hopping), even for weak disorder. We showthat the anisotropic hopping case we consider contains only non-ergodic extended and delocalized states,even for strong disorder. This provides further evidence that cooperative shielding does not exist for theanisotropic hopping considered here, as cooperative shielding is expected to be strongest when α = 0.18We then examine finite α , and show that anisotropic hopping again only supports delocalized and non-ergodic extended states for α = 1. Interestingly, for anisotropic hopping with α = 3 (i.e. dipolar inter-actions), we observe no delocalized states given sufficient disorder. Whether the states are localized ornon-ergodic extended is indeterminate, however. The isotropic hopping case supports delocalized states andeither non-ergodic extended or localized states (or both) for α = 1 and 3.Following this, we use the scaling behaviour of energy level statistics to investigate the physically im-portant case of α = 3. Our results indicate the presence of at least some localized states near zero energy forisotropic hopping, while the anisotropic case is inconclusive. For α = 1 and isotropic hopping, energy levelstatistics also indicate the presence of, at least, some localized states near zero energy. To connect with cur-rently attainable experiments, we also examine the phase diagram for both diagonal and binary off-diagonaldisorder (as per the quantum percolation model) when α = 3. This provides a basic map that can guideexperiments with relevant finite systems, such as polar molecules on an optical lattice [105].2.5.1 ModelsWe consider a single particle in a disordered and diluted cubic lattice with N sites per dimension, with thehopping amplitude∝ 1rα , where α sets the hopping range. The angular dependence of the hopping amplitudeis described below. The Hamiltonian we consider has the following general structure:Hˆ =∑iωicˆ†i cˆi+∑i∑j 6=iti jcˆ†i cˆ j, (2.19)where the operator cˆi removes the particle from site i, ωi is the energy of the particle in site i, and ti j is theamplitude for particle tunnelling from site j to site i. We introduce disorder by randomizing both the valuesof ωi and ti j, which makes Eqn. (2.19) relevant for both disordered lattices and amorphous systems.We randomize the values ωi ∈ [−ω/2,ω/2] by drawing them from a uniform distribution. We definethe disorder strength asW ≡ ωtmax, (2.20)where tmax ≡ max |ti j|, in order to allow a direct comparison between the isotropic and anisotropic models19defined below by normalizing the disorder amplitude to the largest hopping amplitude present.We randomize the values ti j as in the site percolation model. To do this, we define ti j = di jτi j, where di jis the dilution parameter, and divide the lattice sites into two subsets P and Q, with pN3 sites in the P subsetand (1− p)N3 in the Q subset. For a given value p, the lattice sites are assigned to the subsets at random.The dilution parameter is then defined as:di j =1, i ∈ P and j ∈ P0, i ∈ Q and/or j ∈ Q(2.21)With ti j and di j thus defined, the Hamiltonian (2.19) describes a generic particle in a disordered, dilutedlattice with pN3 sites. The angular dependence of the hopping amplitude is determined by the magnitudesof τi j. We consider two types of models for τi j.Isotropic HoppingFor the case of isotropic hopping, we defineτ Ii j =γ|~ri j|α , (2.22)where r i j = r i− r j is the distance between sites i and j in the 3D lattice. With α > 3 the hopping is short-range, with α ≤ 3, long-range. The value of γ is chosen such that τ Ii j ≡ t˜ = 1 for nearest neighbour (NN)sites. Thus, tmax is 1 for isotropic hopping.Anisotropic HoppingWe choose the tensorial form of the anisotropic hopping to be the same as the dipole-dipole interactionbetween polar molecules subjected to an electric field along the z-direction, see, for example, Ref. [120]:τAi j =γ(1−3cos2 θi j)|~ri j|α , (2.23)where θi j is the angle between~ri j and the z-axis. We choose γ as defined above, which makes τAi j = 1 for NNsites in the x-y plane and τAi j =−2 for NN sites along the z-axis. Thus, tmax = 2 for the anisotropic hopping20considered here.Note that, for the remainder of this chapter, the term anisotropic refers specifically to the form inEqn. (2.23), where θi j is the angle between ~ri j and the z-axis, with the z-axis chosen to lie along one ofthe principal axes of the cubic lattice. Changing the alignment of the dipoles to be along another direction(i.e. changing the direction of the electric field) may result in different localization behaviour. Furthermore,many other forms of anisotropy are possible (i.e. quadrupolar-like, octopolar-like, etc.), but are beyond thescope of this chapter.In this section (2.5), the system parameters of interest are p, W , N, α and whether the hopping isiso- or anisotropic. For all of the numerical results, we perform exact diagonalization to obtain all of theeigenvalues and eigenstates of the dense Hamiltonian matrix. These eigenvalues and eigenstates are neededto perform the energy level statistics without any approximations and to examine the participation ratiosacross the entire spectrum. This restricts the size of the system to be smaller than what can be achievedwith state-of-the-art Jacobi-Davidson algorithms, which only obtain a small subset of the total number ofeigenvalues and eigenstates [121, 122]. However, we are still able to observe scaling behaviour for someproperties.2.5.2 Cooperative Shielding in Three DimensionsRecent work has shown the presence of cooperative shielding in 1D single- and many-body systems withlong-range hopping or interactions [30, 119]. Cooperative shielding allows the dynamics of a Hamiltonianwith long-range features to be described by an effective short-range Hamiltonian for a finite time, i.e. thedynamics are effectively “shielded” from the long-range components for a finite time.This phenomenon occurs because of the formation of short- and long-range subspaces in the Hilbertspace, in turn caused by the long-range terms in the Hamiltonian. In many instances, the short-range sub-spaces have many more states than, and are separated by an energy gap from, the long-range subspace.When this gap increases with the system size, it increases the “shielding” time and makes the shieldingcooperative. In the infinite size limit, the dynamics of a wavefunction initially in the short-range subspaceare then completely governed by the effective short-range Hamiltonian. Cooperative shielding allows for thelocalization of particles with long-range hopping [30] and provides an explanation as to why localization isobserved in various systems with long-range hopping or long-range interactions [27, 28, 118].21Here, we show that the energy level structure and gap behaviour conducive to cooperative shielding in1D are also present in 3D cubic lattices, if the long-range hopping is isotropic. If the hopping is of theanisotropic form (2.23), the energy level structure becomes different and does not exhibit well-separatedHilbert subspaces, precluding the cooperative shielding described in [30].10−310−210−1100101102103α10−1010−810−610−410−2100102104106Δ1Dα←α=11D(A)NΔ=Δ113NΔ=Δ153NΔ=Δ213NΔ=Δ253NΔ=Δ313NΔ=Δ101310−310−210−1100101102103α10−210−1100101102103104105106Δα←α=3Isotropic,Δ3D(B)NΔ=Δ11NΔ=Δ15NΔ=Δ21NΔ=Δ25NΔ=Δ31NΔ=Δ10110−1100101Hopping Parameter Exponent, α10−410−310−210−1100101102103104Gap Energy (̃t̃←α=3Anisotropic, 3D(C̃N=10N=15N=20N=25N=301.5 2.0 2.5 3.0 3.5Hopping Parameter Exponent, α10−510−410−310−210−1100101102Gap Energy (̃t̃←α=3Anisotropic, 3D(D̃N=10N=15N=20N=25N=30Figure 2.2: The energy gap at the top of the spectrum as a function of the hopping range exponent α .Panel A is for long-range hopping in a 1D lattice with periodic boundary conditions. ∆1Dα is deter-mined analytically; see Eqn. (2.29) and associated text. Panel B is for a 3D lattice with isotropichopping. ∆α is determined analytically under periodic boundary conditions; see Eqn. (2.27).Panels C and D show the energy gap for the Hamiltonian with anisotropic hopping. The gapenergy is determined numerically and under open boundary conditions. Panel D is an expandedview of panel C. N is the system side-length (total of N3 sites in the 3D case). In all plots, p = 1and W = 0. Note that we examine the top of the spectrum as the effective mass is negative.Isotropic HoppingTo illustrate cooperative shielding in 3D, we first diagonalize analytically the Hamiltonian for an ideal latticewith hopping of arbitrary range. The model with the isotropic hopping (2.22) can be written, in the absence22of disorder and with periodic boundary conditions, asHˆIα =N−12∑i, j,h,l,m,n=−N−12ζrα(1−δilδ jmδhn)cˆ†lmncˆi jh, (2.24)where ζ > 0 is the hopping parameter, the lattice side-length N is assumed odd for simplicity, cˆi jh is theannihilation operator on site (i, j,h), r =√(l− i)2+(m− j)2+(n−h)2 is the unit-less distance betweensites, and the term (1−δilδ jmδhn)rα is understood to be zero when i = l, j = m, and h = n (i.e. the particle canhop anywhere except to its original site). Note that ζ > 0 implies a negative effective mass, such as is seenfor holes in a semi-conductor. Given that we examine only non-interacting particles, whether the effectivemass is positive or negative does not change the final results, except for where one would expect to find therelevant shielding gap (i.e. at the top or bottom of the spectrum).The case of α = 0 produces the uniform, isotropic Hamiltonian HˆI0, which can be diagonalized in themomentum basis to yieldHˆI0 = ζN3cˆ†000k cˆ000k −ζ ∑k1k2k3cˆ†k1k2k3 cˆk1k2k3 , (2.25)where the subscript of cˆ000k means k1 = k2 = k3 = 0. The energy gap between this specific state |000k〉 andthe (N3−1)-fold degenerate states |k1k2k3〉, with one or more of k1,k2,k3 not equal to zero, is thus ∆= ζN3for the uniform, isotropic Hamiltonian. This is the energy gap responsible for cooperative shielding. Inthe 1D case [30], this gap is ∆ = ζN, which allows us to surmise that in the 2D case, ∆ = ζN2 , and thatfor any lattice geometry or dimension, ∆ = ζM, where M is the total number of sites. This occurs becausethe Hamiltonian (2.24) with α = 0 for an arbitrary lattice with M sites is proportional to the adjacencymatrix of a complete graph with M nodes. The spectrum of this adjacency matrix is known [123] to containone (M− 1)-fold degenerate eigenvalue and one non-degenerate eigenvalue, which are separated by a gapproportional to M. This further implies that all lattice dimensions and geometries are equivalent for a particlewith uniform isotropic hopping. In agreement with this, the above energy level structure (2.25) has beenobserved by Ossipov, who studied the Anderson model on the d-simplex [29].The limiting case of α→∞ produces the tight-binding model, also diagonal in momentum space. Giventhat the eigenstates of these two limits are momentum states and that the 1rα factor does not break any23additional symmetries, the eigenstates of HˆIα with finite α must also be momentum states. Evaluating thematrix elements of HˆIα in momentum space, we findE Ik1,k2,k3(α) = 〈k1k2k3|HˆIα |k1k2k3〉=2ζN−1∑l=1(1− lN) cos(2pik1lN )+ cos(2pik2lN )+ cos(2pik3lN )lα+4ζN−1∑l=1N−1∑m=1(1− lN)(1− mN) cos(2pik1lN )cos(2pik2mN )+cos(2pik2lN )cos(2pik3mN )+cos(2pik3lN )cos(2pik1mN )(l2+m2)α/2+8ζN−1∑l=1N−1∑m=1N−1∑n=1(1− lN)(1− mN)(1− nN) cos(2pik1lN )cos(2pik2mN )cos(2pik3nN )(l2+m2+n2)α/2 , (2.26)where ki ∈{k|k ∈ Z∧ k ∈ [−N−12 , N−12 ]}refers to the component of the reciprocal lattice vector in the di-rection i, such that the crystal momentum is 2piki/Na, a is the lattice constant, and |k1k2k3〉 is a momentumeigenstate. Details about the derivation can be found in Appendix A. It can be shown that Eqn. (2.26) re-duces to 〈k1k2k3|HˆI0|k1k2k3〉, see Eqn. (2.25), in the limit where α→ 0 and to the tight-binding model resultin the limit where α → ∞ and N→ ∞.The gap for generic α can be defined as∆α = E I0,0,0(α)−E I0,0,1(α) (2.27)Given Eqn. (2.26) and ζ > 0, it is clear that E I0,0,0(α) is the largest eigenvalue and that the next highest energylevel has states with a single unit of momentum. This energy level is six-fold degenerate (i.e. E I0,0,±1 =E I0,±1,0 = EI±1,0,0,∀α), as can be seen from the symmetries of E Ik1,k2,k3(α) in ki; the choice of EI0,0,1(α) in thedefinition of ∆α is arbitrary. Also, since ζ > 0, EI0,0,0 ≥ E I0,0,1, and ∆α ≥ 0.For comparison, in 1D the Hamiltonian is given by:HˆI,1Dα =N−12∑i, j=−N−12ζ|i− j|α(1−δi j)cˆ†i cˆ j. (2.28)24The corresponding eigenvalues are:EI,1Dk (α) = 〈k|HˆI,1Dα |k〉= 2ζN−1∑l=1(1− lN) cos(2piklN )lα, (2.29)where k ∈{k|k ∈ Z∧ k ∈ [−N−12 , N−12 ]}is the component of the reciprocal lattice vector. A correspondinggap can be defined as ∆1Dα = EI,1D0 (α)−E I,1D1 (α) and is shown in Figure 2.2 (A). This 1D Hamiltonian(2.28) has been shown to exhibit cooperative shielding in Ref. [30].Figure 2.2 shows ∆α as a function of α for various system sizes. Panel (A) of Figure 2.2 plots the resultsfor the 1D lattice. Panel (B) shows that, for a 3D lattice with isotropic hopping, the dependence of the gapseparating the long-range and short-range subspaces on α is qualitatively the same as in 1D, except that thegap becomes independent of the system size at α = 3, instead of α = 1. We have confirmed by numericalcalculations that this behaviour is the same for even and odd N and with open boundary conditions. Basedon the 1D and 3D results of Figure 2.2, it should be expected that a 2D system with isotropic hopping likelyhas a similar gap behaviour, exhibiting the transition at α = 2.The Anderson transition is known to exist in 3D for α > 3 [24], but questions arise for the case withα ≤ 3. Figure 2.2 (B) shows that cooperative shielding can be expected for isotropic hopping with α < 3,suggesting that the Anderson transition could also occur in this region. As a side note, the localization maynot be exponential, as seen in [28].Note that the curves in Figure 2.2 (B) corresponding to different lattice sizes all converge at the transitionpoint. Thus, as in 1D [30], the gap is independent of lattice size when α equals the dimension (in our case,three). If the disorder strength is smaller than the gap, then the gap between the long- and short-rangesubspaces remains open. A finite gap would seemingly indicate a finite shielding time, perhaps precludinglocalization at infinite times for α = 3. However, consider a quantum particle placed at an individual latticesite. The contribution of the delocalized state |000k〉 to the wave packet of this particle is proportional to1N , indicating a zero contribution in the infinite lattice size limit. This leaves the possibility of cooperativeshielding and localization open in the infinite size limit, but does not confirm them. In the following sections,we examine the effect of cooperative shielding on localization in 3D lattices.2510−310−210−1100101102103Hopping Parameter Exponent, α100101102103104Energy ( t̃Triplet →←α=3Isotropic, 3D10−310−210−1100101102103Hopping Parameter Exponent, α100101102103104Energy (̃t̃Doublet←α=3Anisotropic, 3DFigure 2.3: Energy of several top eigenstates as a function of α for isotropic hopping and anisotropichopping with open boundary conditions (left and right panels, respectively). For the left panel,the red dashed line is triply degenerate; the other is non-degenerate. For the right panel, theorange line is doubly degenerate for α . 1; the rest are non-degenerate. In both plots, N3 = 303,p = 1, and W = 0.Anisotropic HoppingThe dependence of the gap separating the short- and long-range subspace is markedly different for themodel with anisotropic hopping (see Figure 2.2 (C) and (D)). It is more difficult to obtain the analyticalresults for the model with anisotropic hopping because of the cos2 θ term, even if the eigenstates of thesystem with periodic boundary conditions are still momentum states (rotational symmetry is broken, butnot translational). Here, we study the behaviour of the gap at the top of the spectrum as a function of αnumerically, using open boundary conditions. The results shown in Figure 2.2 (C) illustrate the qualitativedifference from the isotropic case. In particular, for the anisotropic model, there are several values of α < 3where the gap becomes zero. At these values of α , cooperative shielding may not occur, thus precludinglocalization.From examining the energy level structure, it can be concluded that the zero-gap points arise from levelcrossings. It is thus clear that there is no separation of subspaces at the top of the spectrum, in contrast tothe case of isotropic hopping. While there is a gap that grows with the system size for α . 2, the state at thetop of the spectrum is not an equal superposition of all sites, as it is when cooperative shielding is knownto occur. Furthermore, the top energy level is shown to be a doublet under periodic boundary conditions,26indicating that this particular gap actually disappears in the infinite size limit when translational invarianceis restored.There is, however, a gap between the highest-energy doublet and the next doublet. We examine this gapfor N = 11 and periodic boundary conditions. The curve (not shown) appears similar to those in Figure 2.2(C) and (D), though the number of zero-gap points is reduced. This is likely because the higher symmetrypermits the formation of doublets, which produce fewer level crossings (i.e. one for each doublet insteadof one for each state). Given these level crossings, the same conclusions can be reached for the periodicboundary case as for the open boundary case: there is no clear separation of subspaces at the top of thespectrum and no evidence for cooperative shielding. While there is a gap at small values of α , the states atthe top of the spectrum are not equal superpositions of all sites and not necessarily fundamentally differentfrom the rest of the eigenstates. Cooperative shielding, of the form discussed in [30], thus appears precluded.This is investigated more fully further below by examining the scaling behaviour of the participation ratios.In addition to the gap behaviour, the overall energy level structure appears different between the twocases, as exemplified in Figure 2.3 for ideal lattices of side-length N = 30. Here, the top several eigenstatesare shown for both isotropic (left panel) and anisotropic (right panel) hopping with open boundary condi-tions. The values are obtained numerically. Both the degeneracies of the levels and the overall behaviouras a function of α are very different. The energy level crossings in the anisotropic hopping case are in theregion 2 . α . 3. Their exact locations are best observed as the zero-gap points in Figure 2.2 (D). Notethat the differences between Figure 2.3 (left panel) and Eqn. (2.26), particularly in terms of the degeneracylayout, are due to the differing boundary conditions (open vs. periodic).2.5.3 Participation RatiosWe have established the presence of cooperative shielding in 3D systems with isotropic long-range hopping,which suggests the possibility of localization. We now examine this possibility, for both isotropic (2.22)and anisotropic (2.23) hopping with open boundary conditions and α = 0,1, and 3, from a participationratio point of view. In a later section, we examine the energy level statistics for α = 3.The participation ratio (PR) is an effective count of the number of sites occupied by a wavefunction. It27Isotropic, α = 0Figure 2.4: Participation ratio (PR) as a function of energy for the eigenstate of a 3D system withisotropic hopping, W = 16.5, p = 1, and α = 0. The number of disorders included ranges from148 to 2100, as required to obtain sufficiently small error bars. Left Panel: Plot of the PR ofeach eigenstate vs its eigenvalue for several disorder realizations for N = 11, 21, and 31 (green, orange , and blue , respectively). Note the symmetric logarithmic scaling of the abscissa,where the linear portion is within ±10t˜. Large black circles merely highlight the location of thestates at the top of the band. Right Panel: PR averaged within energy bins and over disorderrealizations. The plot is focused on the band near zero energy. Note the linear scale of theabscissa. Error bars for the average PR are 95% confidence intervals and smaller than the markersize where not seen.is defined asPR =1∑i |ψ(~xi)|4, (2.30)where the sum is over all lattice sites, ψ is a normalized wavefunction and~xi is the position of lattice site i.Typically, PR ∝ Nβ [124–126]. This scaling behaviour gives information on the nature of the eigenstates.In particular, if β = 0, the eigenstates are localized, if β = 3, delocalized and, if 0 < β < 3, the statescan be labelled as extended non-ergodic [124–126]. Formally, extended non-ergodic states are states wherePR−1 6→ 〈PR−1〉 as N → ∞, where 〈· · · 〉 indicates the average over disorder realizations. Informally andin the context of this thesis, extended non-ergodic states span the system, but do not fill the entirety of theaccessible space. Extended non-ergodic states are typically multi-fractal, though this property is definitivelyidentified by examining the generalized participation ratios [125–127]. We focus on the scaling behaviourof the PR (2.30), rather than the generalized participation ratios, as the limited system sizes accessible whenobtaining the entire spectrum via exact diagonalization makes a full quantitative analysis difficult.28We analyze the scaling behaviour of the PR of the eigenstates via two approaches. The first obtains aquantitative estimate of β (E) by fitting the above scaling equation to the behaviour of 〈PR〉, where 〈PR〉is the PR averaged within a particular energy window and over multiple disorder realizations. The sec-ond approach obtains a qualitative understanding of the system by examining the scaling behaviour of thedistributions of both PR and PR/N3 at all energies. Three cases are discerned:(i) If the PR does not change with the system size, but PR/N3 does, then β = 0.(ii) If the PR scales with the system size, but PR/N3 does not, then β = 3.(iii) If both PR and PR/N3 change with the system size, then 0 < β < 3.No other cases arise as β ∈ [0,3].We note that the spectra of the systems we examine have an interesting dual-scale behaviour: there is alarge density of states near zero energy, but also many states at energies very far away from zero. Thus, thespectra typically have one or more “wings”. To capture the behaviour of the PR across all of these energyscales, a symmetric logarithmic scale is used for the abscissa, where the portion of the axis near zero isscaled linearly, while the portions far from zero are scaled logarithmically. The size of the linear portion isdenoted on the abscissa of the graph and in the captions.Scaling of Participation Ratios when α = 0: Isotropic HoppingFigure 2.4 (left panel) shows the PR of every eigenstate, for multiple disorder realizations, plotted againstthe energy of the eigenstate for W = 16.5. The band at low energy is of width 16.5t˜ and arises from the(N3− 1)-fold degenerate ground state of the ideal lattice. The points at the top of the spectrum arise fromthe |000k〉 state of the ideal system. Figure 2.4 (right panel) shows the average PR as a function of energyfor the band at low energies. There is clearly no scaling of 〈PR〉 with the system size. The fit to the 〈PR〉 ofthe high energy |000k〉 state shows the scaling exponent to be 3. Identical behaviour has been observed forW = 1 and 35.Evidently, the band at low energy is composed of localized states, while the |000k〉 state remains delo-calized, even in the presence of strong disorder. Such behaviour is easily understood from the energy levelstructure of the ideal system (described in Section 2.5.2). The extreme degeneracy of the ground state al-lows any small amount of disorder to strongly couple the states, leading to localization [29]. The |000k〉 state29remains delocalized, however, as the shielding gap is significantly larger than the disorder strength, mini-mizing coupling. Given that the shielding gap diverges with the system size, the |000k〉 state must remaindelocalized in the thermodynamic limit for any finite disorder strength, while the low energy states remainlocalized. The same behaviour is observed for particles with isotropic uniform hopping on the d-simplex[29].Given that cooperative shielding allows the system to be described by an effective short-range Hamilto-nian [30], one could expect to see the localization transition near W = 16.5, the critical disorder strength ofthe Anderson transition for nearest-neighbour hopping in 3D lattices. Indeed, the dynamics of the systemstudied by Celardo et al. [30] is governed by an effective nearest-neighbour hopping Hamiltonian. However,they explicitly include a nearest-neighbour hopping term in their Hamiltonian, in addition to the long-rangeterm. As it would be inconsistent with our generalization of the dipolar hopping amplitude to arbitrary α , wedo not include an additional nearest-neighbour hopping term in our Hamiltonian (2.19). Thus, the groundstate of our system when α = 0 is extremely degenerate (essentially flat-band) and localizes given any finitedisorder strength.Scaling of Participation Ratios when α = 0: Anisotropic HoppingThe upper panels of Figure 2.5 show the PR vs. energy for W = 1 and 16.5 (left and right, respectively),while the lower panels show the averaged PR at the same disorder strengths. One can clearly see scalingwith the system size for both disorder strengths. Given that the bandwidth of the system is a function ofthe system size, the scaling behaviour cannot be quantitatively determined for all band energies. Clearlyacceptable choices of energies for comparison include the top and bottom states of the band and the statesnear the centre. The wing regions, however, do not overlap well enough for a quantitative analysis. Fitsto the scaling equation reveal that the top and bottom states at both disorders scale as N3. The states inthe energy region [−1,1] for the system with W = 1 also scale as N3. Interestingly, the states in the centreportion of the system with W = 16.5 have energy-dependent scaling exponents that vary from 2.6 to 2.8.We also examine the system W = 35, but find it qualitatively identical to the W = 16.5 case: the scalingexponents of the top and bottom states are 3, while those of the states near the band centre vary from 2.4 to2.7.We note that the errors on the scaling exponents can be as large as 0.2. We obtain this upper bound on30Anisotropic, α = 0Figure 2.5: Participation ratio (PR) as a function of energy for the eigenstate of a 3D system withanisotropic hopping; W = 1 and 16.5; p = 1; α = 0; and N = 11, 21, and 31 (green , orange, and blue , respectively). The number of disorders included ranges from 50 to 2350, asrequired to obtain sufficiently small error bars. Note the symmetric logarithmic scaling of theabscissa, where the linear portion is within ±1t˜ or within ±16t˜. Top Panels: Plot of the PR ofeach eigenstate vs its eigenvalue for several disorder realizations. Lower Panels: PR averagedwithin energy bins and over disorder realizations for W = 1 and 16.5. Error bars for the averagePR are 95% confidence intervals and smaller than the marker size where not seen.the error from one fit to the scaling equation that resulted in a value of 3.2, while β cannot be larger than 3.Given that we only have three system sizes that span a relatively small range of sizes, our scaling exponentshave low, but not negligible, precision. Qualitatively, an examination of PR/N3 (not shown) shows that thewings, for both disorder strengths, scale as N3. The central states also scale as N3 for W = 1, while theyscale with 0< β < 3 for the central states at W = 16.5. These qualitative results support the accuracy of the31quantitative results discussed above, despite their low precision.The above results show that the states in the wings are delocalized, even for strong disorder. Further-more, the states near the centre of the band become extended non-ergodic for sufficient disorder strength.We also see that there are no localized states in this system, in great contrast to the isotropic case.α = 1Figure 2.6: Top Panels: Plot of the Participation ratio (PR) of each eigenstate vs its eigenvalue forseveral disorder realizations. Lower Panels: Plot of the PR/N3 of each eigenstate vs its eigenvaluefor several disorder realizations. Plots are for a 3D system with isotropic and anisotropic hopping;W = 16.5; p = 1; α = 1; and N = 11, 21, and 31 (green , orange , and blue , respectively).The number of disorders included ranges from 100 to 350, as required to obtain well-formeddistributions. Note the symmetric logarithmic scaling of the abscissa, where the linear portion iswithin ±10t˜ or within ±16t˜.32Scaling of Participation Ratios when α = 1,3We examine the scaling behaviour of the PR for isotropic and anisotropic hopping with α = 1 and α = 3 byexamining the qualitative scaling behaviour of PR and PR/N3. We do not examine the quantitative behaviouras the limited system sizes accessible causes fits to the scaling equation to be even less reliable than whenα = 0. Figure 2.6 demonstrates that, for α = 1 and W = 16.5, there is a single wing for the isotropic case,while there are two for the anisotropic case. These wings contain PR values that increase with the systemsize, while the PR/N3 values are constant: the values scale as N3 and the states are delocalized.In the isotropic case near zero energy, however, the behaviour of the PR is not as clear. A portion of thePR values in this energy range does scale with the system size, but the distribution asymmetrically increasesin width instead of shifting vertically. This obfuscates the interpretation. However, given that the PR/N3values all reduce in magnitude as the system size increases, it can be concluded that none of these centralstates are delocalized. Whether the states are extended non-ergodic or localized is not clear, however.In the anisotropic case with α = 1, the states near zero energy are clearly extended non-ergodic. This isseen as both the PR and the PR/N3 distributions shift vertically as the system size increases, in the appropri-ate directions. Additionally, note that the rate of scaling with the system size near zero energy is differentbetween the isotropic and anisotropic cases. The isotropic case has PR distributions that scale slowly, whilethe PR/N3 distributions scale quickly. The reverse is true for the anisotropic case. We can thus concludethat the isotropic hopping scaling exponent is smaller than the anisotropic hopping scaling exponent for thestates near zero energy.Figure 2.7 demonstrates that, for α = 3 and W = 16.5, there is a single wing for the isotropic case, whilethere are none for the anisotropic case. Also note that the bandwidths are much smaller, allowing for a linearscaling of the x-axis. The behaviour of the isotropic case is similar to that seen for α = 1 and points to thesame conclusion: the wing contains delocalized states, while the energy region near zero contains states thatare not delocalized.At this level of disorder, the spectrum of the anisotropic case is dominated by the disorder, producinga rather uniform distribution without wings. Interestingly, the scaling behaviour is similar to the near-zero-energy portion of the isotropic case: the PR distributions widen asymmetrically but do not shift in locationas the size increases, while the PR/N3 distributions shift downward as the system size increases. There are33thus no delocalized states, though there may be a combination of localized and extended non-ergodic states.This is not necessarily in contrast with the prior work of Levitov [25, 26, 113] as the extended non-ergodicstates still span the system, allowing weak transport.The above analysis of the participation ratios for various values of α and for different isotropies showsthat single-particle systems with long-range hopping display a rich landscape of states, from delocalized toextended non-ergodic to localized. We expect that the extended non-ergodic states we typically find nearzero energy in this work are likely multi-fractal. To find and quantify the multi-fractal dimensions andmulti-fractal spectra [125] of these states, however, will require more sophisticated numerical or analyticaltechniques capable of accessing larger system sizes. There are algorithms that can reach system sizes ofN ≈ 300 for targeted energy ranges [121, 122], but it is not clear whether even these would be enough toproperly study the scaling behaviour of the eigenstates.2.5.4 Energy Level StatisticsGiven the ambiguity observed in the scaling behaviour of the participation ratios for the physically importantcase of α = 3, we now examine the system for both isotropic (2.22) and anisotropic (2.23) hopping andopen boundary conditions from an energy level statistics point of view.We use the mean level spacing ratio, defined by Oganesyan and Huse [128] as〈r〉= 〈min(r′n,1r′n)〉 (2.31)wherer′n =En+1−EnEn−En−1 ,En is the nth eigenvalue, sorted in order of increasing magnitude and 〈·〉 denotes the average over all eigen-values and over many instances of disorder. For a system whose energy level spacing distribution is Pois-sonian, 〈r〉 = 2ln2− 1 ≈ 0.38629. If, however, the system belongs to the Gaussian Orthogonal Ensemble(GOE), 〈r〉 = 4− 2√3 ≈ 0.53590 [129]. As discussed in Refs. [128, 130], Poissonian level statistics arecharacteristic of localization and GOE statistics are characteristic of diffusion.34α = 3Figure 2.7: Top Panels: Plot of the Participation ratio (PR) of each eigenstate vs its eigenvalue forseveral disorder realizations. Lower Panels: Plot of the PR/N3 of each eigenstate vs its eigenvaluefor several disorder realizations. Plots are for a 3D system with isotropic and anisotropic hopping;W = 16.5; p = 1; α = 3; and N = 11, 21, and 31 (green , orange , and blue , respectively).The number of disorders included ranges from 99 to 350, as required to obtain well-formeddistributions. Note the linear scaling of the abscissa.Figure 2.8 presents 〈r〉 vs W for the isotropic and dipolar Hamiltonians with α = 3 and open boundaryconditions. It can be seen in both cases that 〈r〉 drops as W increases, though not completely to the Poisso-nian value. However, in the isotropic case, 〈r〉 drops with the increasing lattice size towards the Poissonianlimit. As illustrated in Figure 2.8 (left panel), we observe a significant change with the system size. Incontrast, in the anisotropic case, the curves display no discernible scaling with the system size, even in thecase of extremely strong disorder. There is thus no indication that the Poissonian limit will be attained.It is noteworthy that the scaling behaviour in the isotropic case begins at W ≈ 15, which is the expected35location, given cooperative shielding. Cooperative shielding causes the Hamiltonian to be effectively short-ranged, such as the tight-binding model. For the tight-binding model, the diffusion-to-localization transitionoccurs at W ≈ 16.5 [131]. At disorder strengths between 5 and 15, the curves in the isotropic case collapseonto one another, preventing the possibility of scaling to the infinite size limit. At W = 1, it is clear that thesystem is approaching the GOE value as the system size increases.Given that the scaling behaviour of the PRs reveals two clear regimes in the isotropic case, one wherethe states are delocalized and one where they are localized or extended non-ergodic (see Figure 2.7), it isuseful to focus on the region near zero energy. Figure 2.9 (left panel) demonstrates that the scaling behaviourobserved previously is enhanced when the energy level spacing ratio is averaged only over the energy region−t˜ ≤ E ≤ t˜. Figure 2.9 (right panel) demonstrates that the same scaling behaviour also occurs when α = 1.These results provide numerical evidence for localization in the infinite size limit for the case with isotropiclong-range hopping and α = 1 and 3. This, combined with the likely presence of cooperative shielding andthe observed scaling behaviour of the participation ratios, strongly suggests that some, to many, of the statesat low energy are localized.Full Spectrum0 10 20 30 40 50W0.380.400.420.440.460.480.500.520.54⟨r⟩NIsotropic⟩⟨3DGOEPoissonianN⟨=⟨11N⟨=⟨21N⟨=⟨310 10 20 30 40 50W0.380.400.420.440.460.480.500.520.54⟨r⟩No ScalingAnisotropic, 3DGOEPoissonianN = 11N = 21N = 31Figure 2.8: Mean energy level spacing ratio 〈r〉 as a function of disorder strength W = ωtmax for theisotropic and anisotropic hopping Hamiltonians with open boundary conditions. The hoppingrange exponent is α = 3 and the filling fraction is p = 1 in both cases. N denotes the lattice side-length (N3 sites). The horizontal dashed lines denote the values of 〈r〉 for the Poisson distribution(∼0.38629) and the GOE (∼0.53590). The error bars are 95% confidence intervals based on 96disorders and are smaller than the marker size where not seen.36Near Zero Energy0 10 20 30 40 50W0.380.400.420.440.460.480.500.520.54⟨r⟩Nα=3⟩ Isotropic, 3D−̃t≤ E≤̃tGOEPoissonianN = 11N = 21N = 310 10 20 30 40 50W0.380.400.420.440.460.480.500.520.54⟨r⟩Nα=1⟩ Isotropic, 3D−̃t≤ E≤̃tGOEPoissonianN = 11N = 21N = 31Figure 2.9: Left panel: Same as Figure 2.8 (left panel), but with the mean energy level spacing ratio〈r〉 including only the eigenvalues in the energy range−t˜ ≤ E ≤ t˜. The results are averaged over1071 disorders for N3 = 113 and 213, and 96 disorders for N3 = 313. Right panel: Same as leftpanel, but for α = 1. The results are averaged over 1100 disorders for N3 = 113 and 213, and 100disorders for N3 = 313.To connect with realistic experimental conditions, we also examine the effect of a diluted lattice. Thisis relevant for experiments with polar molecules in an optical lattice [104, 105] or for amorphous molecularsolids [100]. A diagram of the mean level spacing ratio, 〈r〉, as a function of both p and W for both isotropicand dipolar hopping with α = 3 and open boundary conditions is presented in Figure 2.10.While the two plots are qualitatively similar, 〈r〉 falls off more quickly with the disorder strength andfilling fraction in the isotropic case. This can be seen by comparing the areas of the two diagrams where〈r〉> 0.49. The value of 〈r〉 remains high along the p-axis, which is in agreement with the results of Denget al. [124], who studied the case of dipolar hopping in a diluted lattice.Figure 2.10 (right panel) shows that in experiments with polar molecules in optical lattices one shouldnot expect to see localization of single rotational excitations due to the dilution alone, if it exists at all.Additional on-site disorder, such as from an optical speckle potential, is required. Exploring this phasediagram is within reach of current experiments. For example, Figure 2.10 (right panel) shows that formolecules on an optical lattice with N ≈ 30 and a lattice population of 30 % [105] (a typical filling fractionaimed at in current experiments), the region where 〈r〉 drops significantly can be explored by varying the37optical speckle potential from below to above W = 5. Scaling behaviour of 〈r〉 can also be investigated asoptical lattices can have N up to 60 (a typical size of an optical lattice).0.05 0.2 0.4 0.6 0.8 1.0p05101520WIsotropic, 3D⟨r⟩ ⟩ 0.4⟨0.380.400.420.440.460.480.500.52⟨r⟩0.05 0.2 0.4 0.6 0.8 1.0p05101520WAnisotropic, 3D⟨r⟩ ⟩ 0.4⟨0.380.400.420.440.460.480.500.52⟨r⟩Figure 2.10: Mean energy level spacing ratio 〈r〉 as a function of disorder strength W = ωtmax and fillingfraction p for the isotropic and anisotropic hopping Hamiltonians with open boundary condi-tions. The hopping range exponent is α = 3 and the lattice size is N3 = 313 in both cases.The colour varies from the value of 〈r〉 for the Poisson distribution (∼0.38629) to that for theGOE (∼0.53590). The white circles denote the values obtained via exact diagonalization ofthe Hamiltonian; piecewise cubic interpolation is used for the rest of the plot. The black linesindicate where 〈r〉= 0.49. The results are averaged over 96 realizations of disorder, giving 95%confidence intervals that are at most ±0.002. The green data points mark the location of theAnderson transition for the 3D tight-binding model in the infinite size limit, as determined byRoot and Skinner [131].2.6 ConclusionWe have illustrated that the cooperative shielding discovered for 1D lattices with long-range hopping [30,119] also occurs in 3D lattice models with isotropic long-range hopping. This suggests the possibility ofAnderson localization in 3D systems with long-range hopping (i.e. α ≤ 3). Given the form of the results inthree dimensions, the same is likely to be true of a 2D system with isotropic long-range hopping (α ≤ 2).We have also presented evidence indicating the lack of cooperative shielding in models with anisotropichopping. The specific form of the anisotropy we consider here has the angular dependence of the dipole-dipole interaction between dipoles aligned along one of the principal axes of the lattice.We have demonstrated that there are fundamental differences between disordered lattice systems withisotropic long-range hopping and those with anisotropic long-range hopping, particularly with regards to38the impact of disorder. In addition to the difference in cooperative shielding, we have shown that the energylevel structures of systems with isotropic hopping are qualitatively and quantitatively different from those ofsystems with dipolar anisotropic hopping (whether with periodic or open boundary conditions).We have used the scaling behaviour of the eigenstate participation ratios to demonstrate the presence oflocalized states in the isotropic case with uniform hopping (α = 0) and the absence of localized states inthe anisotropic case with uniform hopping (α = 0). We have also demonstrated that energy level statisticssupport the presence of localized states in 3D cubic lattices with isotropic hopping that varies as r−3, whilethey are inconclusive for systems with the anisotropic dipolar hopping considered here. We have furthershown that energy level statistics indicate the presence of localized states for Coulomb-like (α = 1) isotropichopping, while the scaling behaviour of the eigenstate participation ratios illustrates the absence of localizedstates when α = 1 and the hopping is of the anisotropic form (2.23).We have shown that the localization properties (or lack thereof) found for systems with dipolar hopping[25, 26, 113, 114] cannot be assumed to apply to systems with isotropic hopping, in accordance with thesuggestion of Burin and Maksimov [27]. More generally, we have demonstrated that the presence or absenceof delocalized, extended non-ergodic, or localized states depends on not only the range of the hoppingamplitudes, but also their anisotropy.2.7 Future DirectionsFuture studies that can access significantly larger system sizes are required to fully characterize the shapeof the localized states, to determine whether the extended non-ergodic states are truly multi-fractal, and todetermine how the multi-fractal dimensions and spectra change with hopping range and isotropy. Futuretheoretical and experimental studies can also probe the dependence of localization on the direction of thequantization axis (as defined by the direction of an external field) relative to the principal axes of the under-lying cubic lattice. Indeed, it would be interesting to investigate whether cooperative shielding in 3D latticescan be induced by tilting the quantization axis. Additionally, it would be interesting to explore the effectof the lattice geometry on cooperative shielding. Finally, this work opens up the question of how differenttypes of anisotropy, such as quadrupolar-like, can impact localization behaviour.Furthermore, while some of our results indicate that cooperative shielding likely also exists in two39dimensional systems, it is important to directly confirm this and to determine the impact of anisotropy oncooperative shielding and localization in this case. Examining the impact of different lattice geometries andof stacked two dimensional lattices on cooperative shielding and localization would also be interesting. Thismay give insights into the localization behaviour (or lack thereof) of excitons in stacked monolayers [132].Additionally, it is not currently known why anisotropic hopping seems to preclude cooperative shielding.It may be because the average of the hopping amplitudes over angle is zero, which may indicate disruptionof the constructive interference required to induce localization. That is, the two paths on a loop may nolonger interfere in the same manner as in the isotropic case. This possibility needs to be explored further,potentially by taking a diagrammatic approach or simply examining the impact of a hopping amplitudeanisotropy that does not average to zero on cooperative shielding and localization.Also, a recent paper published in 2019 by Nosov, Khaymovich and Kravtsov [133] discusses the no-tion of correlation-induced localization where increasing correlation in the long-range hopping amplitudesdecreases the disorder strength needed to localize the system, to the point where localized states exist forany finite disorder strength. Indeed, it may be that changing the anisotropy of the hopping amplitudes is ef-fectively changing the correlations present in the hopping amplitudes and thus the constructive interferencepattern required for localization. This is a direction worth pursuing.This now concludes my discussion on Anderson localization and how it relates to many-body molecularsystems. Namely, ensembles of molecules can have excitations that effectively behave as particles withlong-range hopping. The nature of this long-range hopping can greatly impact the localization properties ofthese quasi-particles and thus energy transport within the molecular ensemble.402.8 SummaryIn this section, I list the key results presented in Chapter 2. More detailed results can be found in Section2.6.• The presence or absence of delocalized, extended non-ergodic, or localized states depends on both therange and the anisotropy of the hopping amplitudes.• Cooperative shielding extends to three-dimensional systems, if the hopping amplitudes are isotropic(Section 2.5.2)• Cooperative shielding is also likely to exist in two dimensional systems, if the hopping amplitudes areisotropic (Section 2.5.2)• Cooperative shielding appears absent if the hopping is anisotropic with the angular form of the dipole-dipole interaction between dipoles aligned along one of the principal axes of the cubic lattice (Sec-tion 2.5.2)• Isotropic hopping on a disordered cubic lattice permits the formation of localized states, even if thehopping is long-ranged (α ≤ 3), as evidenced by the scaling of eigenstate participation ratios and ofenergy level statistics (Sections 2.5.3 and 2.5.4)• The anisotropic hopping case we consider appears to prevent the formation of localized states whenthe hopping is long-ranged (α < 3), as evidenced by the scaling of participation ratios and apparentlack of cooperative shielding (Sections Section 2.5.2 and 2.5.3)• The evidence provided by the scaling of participation ratios and of energy level statistics is inconclu-sive regarding whether the states are localized or extended non-ergodic for the anisotropic case whenthe hopping is critically long-ranged (α = 3). Extended non-ergodic states would be consistent withthe literature. (Sections 2.5.3 and 2.5.4)41Chapter 3Molecular hyperfine interferometry3.1 IntroductionThis chapter discusses our proposal to use closed shell molecules as probes in surface spin echo experi-ments and the fully quantum framework we have developed to model these experiments. Particularly, Idiscuss 3He spin echo experiments and their extension to molecular hyperfine interferometry experiments atthe start of Section 3.2. In Section 3.3, I describe a generic molecular hyperfine interferometry experiment.I then discuss, in Section 3.4, the molecular state after the state-selecting magnetic lens. In Section 3.5, wetime-evolve the molecular state and integrate the result over the length of the detection window to obtainthe relationship between the system eigenstates and the detector current. To obtain the system eigenstates,we derive and apply, in Section 3.6, a transfer matrix formalism that includes internal degrees of freedom. Ialso discuss the rotation and scattering transfer matrices used to account for the apparatus geometry and themolecule-surface interaction, respectively. In Section 3.7, we demonstrate the application of this theoreticalframework to the case of oH2, illustrate the need to integrate over the velocity distribution, illustrate the sen-sitivity of the calculated signal to various features of the scattering transfer matrix and perform a preliminarycomparison with experiment. We compare the method of the present chapter to the semi-classical methoddiscussed by Godsi et al. [44] in Section 3.8. I conclude the chapter in Section 3.9 and give a bulletedsummary of the chapter results in Section 3.10.423.2 From Helium-3 Spin Echo to Molecular Hyperfine InterferometryUnderstanding the dynamics of molecule-surface interactions is of central importance for many applicationsin science and technology. One of the most powerful experimental techniques for studying surface dynamicsis 3He spin echo (HeSE) [134, 135]. In a HeSE experiment, a beam of spin-polarized 3He atoms is passedthrough two regions of strong magnetic field before and after the 3He atoms meet a sample surface. Thedependence of the spin-echo signal on one or both of the magnetic fields encodes information about theatom-surface scattering event. HeSE experiments and related helium atom scattering experiments fulfil aunique role amongst other experimental techniques for studying surfaces as they can probe both fast time-scales (pico- to nanosecond) and small distances (angstroms to tens of nanometers) and they are explicitlysurface sensitive [134].HeSE experiments [134, 135] have been used to study surface morphology [136], molecular and atomicsurface diffusion [33–39], inter-adsorbate forces [37, 137], phonon dispersions [40, 41], phason dispersions[42], structures and phase transitions of ionic liquids [43] and friction between adsorbates and surfaces[138, 139]. HeSE experiments have improved our understanding of potential energy surfaces [140–142]and surface-adsorbate interactions [143–146] and are frequently combined with microscopic calculations toboth test theory and gain insight into surface-adsorbate interactions [147, 148]. However, the interactionsbetween 3He and surfaces or their adsorbates are typically isotropic and relatively weak.One can potentially obtain further information about surface dynamics by using molecules instead of3He atoms in a surface spin-echo experiment. Molecules offer rotational degrees of freedom, which couldbe exploited to gain new insights into surface dynamics. Indeed, Godsi et al. [44] have recently shown thatortho-hydrogen (oH2) can be used as a sensitive probe of surface morphology: their experiment was ableto discern how the interaction between an oH2 molecule and a Cu(115) surface depends on the orientationof the rotational plane of the hydrogen molecule relative to the surface. In addition, one could exploitthe transfer of rotational energy from the probe molecules to surface adsorbates (or vice versa) in orderto study the relative effects of the rotational and translational motion on the dynamics of the adsorbates.However, the increased complexity of molecules (compared to 3He atoms) makes the analysis of the spin-echo experiments complicated and requires one to account for the interplay of the translational, nuclearspin and rotational degrees of freedom in strong magnetic fields of differing orientations, in addition to the43molecule-surface scattering event.In this chapter, I present a general and efficient theoretical framework that can be used to analyze surfaceinterferometry experiments that use closed shell molecules to study static surfaces. The large manifoldof molecular states and magnetic field-induced couplings between internal states preclude a simple semi-classical description for these molecular hyperfine interferometry (MHI) experiments, in contrast to HeSEexperiments. We develop a fully quantum mechanical model of surface-sensitive MHI experiments byderiving a one-dimensional transfer matrix method that includes the internal hyperfine degrees of freedomof the probe molecules. This allows us to efficiently solve the time-independent Schro¨dinger equation witha multiplication of matrices. We account for the experimental geometry with rotation matrices and describethe molecule-surface interaction with a scattering transfer matrix (a transformed version of the standardscattering matrix).We apply our method to an oH2 hyperfine interferometry experiment and illustrate the importance ofintegrating over the velocity distribution of the molecules. We further show that information about thescattering matrix elements is encoded in the experimental signal. In particular, we demonstrate that theexperimental signal is sensitive to the magnitude and phase of the diagonal elements of the scattering transfermatrix. We also show that the signal is sensitive to inelastic scattering events that change the projectionquantum numbers of the molecular hyperfine states, described by scattering transfer matrices with non-zerodiagonal and off-diagonal matrix elements. This sets the stage for determining, in part or in whole, thescattering transfer matrix elements of a particular molecule-surface interaction by comparing the computedand experimentally-measured signals. Finally, through comparison with the semi-classical method discussedby Godsi et al. [44], we demonstrate that our method opens a path towards an even more powerful frameworkthat can analyze surface interferometry experiments using closed shell molecules to study dynamic, insteadof static, surfaces.3.3 Description of a Molecular Hyperfine Interferometry ExperimentA surface-sensitive molecular hyperfine interferometer uses a beam of molecules to probe various surfaceproperties. To do this, a set of magnetic fields are used to simultaneously manipulate the internal hyperfinestates of the probe molecules and create a spatial superposition of molecular wavepackets. These wavepack-44ets sequentially impact the sample surface and scatter in all directions. A second set of magnetic fieldscollects the molecules scattered in a narrow solid angle. This second set of magnetic fields further manip-ulates the molecular wavepackets, partially recombining them and allowing for molecular self-interference.Wavepackets with particular hyperfine states are then passed into a detector. A schematic of the experimentis depicted in Figure 3.1. We now discuss the different stages of the experiment in more detail.The beam source must produce a continuous (or pulsed) beam of molecules with a sufficiently narrowvelocity profile, mean velocity suitable for a particular experiment, sufficiently high flux, and a densitylow enough to ensure that the molecules are non-interacting. One current apparatus [44] uses a supersonicexpansion to produce such a beam. One can also envision experiments with slow molecular beams pro-duced by extraction (sometimes with hydrodynamic enhancement) from a buffer-gas cooled cell [149] orwith molecular beams controlled by electric-field [150] or magnetic-field deceleration [151]. Decelerationprovides control over the mean velocity and narrows the velocity spread [152], which could be exploited fornovel interferometry-based applications.The experiment selects molecules in particular hyperfine states by employing a magnetic lens whosemagnetic field has a gradient in the radial direction. A cylindrically-symmetric field gradient is used toensure sufficient molecular flux. The lens focusses molecules with low-field seeking states and defocussesmolecules with high-field seeking states, allowing for purification of the molecular beam. After the lens,a unique quantization axis for the internal states is developed by using an auxiliary field that adiabaticallyrotates all magnetic moments until they lie along a single direction perpendicular to the beam path. Theend of this auxiliary field is a strong dipolar field aligned along the z direction that defines the quantizationaxis. Hexapole magnets can be used as a magnetic lens as their magnetic field gradients are sufficientlycylindrically symmetric [44, 153, 154]. More details about the internal states of the molecules immediatelyafter the magnetic lens can be found in Section 3.4.Solenoids whose magnetic fields are parallel to the beam propagation path are used to manipulate themolecular hyperfine states. These solenoids are labelled as the control fields in Figure 3.1. Arbitrary mag-netic field profiles can be obtained by changing the solenoid winding patterns and/or using multiple succes-sive solenoids.The hyperfine states of a molecule change energy as the molecule enters a magnetic field. These changes45to the hyperfine energy levels cause simultaneous changes in the molecular momenta, as the total energy isconserved. That is, when molecules enter a solenoid, molecules in low-field seeking states slow down andthose in high-field seeking states speed up. Furthermore, because the direction of the magnetic field in acontrol field is not along the z axis, the molecules are in a superposition of hyperfine states, with respect tothe quantization axis defined by the magnetic field. Thus, the differences in momenta cause the differentcomponents of each molecular wavepacket to spatially separate as the wavepacket traverses the solenoid.Upon exiting the solenoid, the components of each wavepacket return to their original momenta, but remainspatially separated. That is, each wavepacket is now in an extended spatial superposition.Each of these spatially separated wavepacket components comprise a superposition of the field-freehyperfine states. The exact superpositions of each wavepacket component, as well as the spatial separationsbetween the components, depend on the magnetic field profile of the first branch. Each of the wavepacketcomponents sequentially impacts the sample surface and scatters in all directions. However, the experimentonly captures those molecules that pass through a particular solid angle. While a current experiment [44]fixes the angle between the two branches, one can in principle explore many different scattering geometriesby varying both the angle between the two branches of the apparatus and the orientation of the sample.After scattering, the collected molecules enter another set of control fields in the second branch of theapparatus. The hyperfine states again change in energy and momenta. In a helium-3 spin echo experiment,if the second magnetic field profile is identical but opposite in direction to the magnetic field profile of thefirst branch, the spatially separated wavepacket components realign (to first order) as they traverse the mag-netic field(s), producing a spin echo. This allows the wavepacket components to interfere with each other.Interestingly, it has recently been shown [155] that echoes are also produced when the device operates withthe fields in the same direction. With an arbitrary hyperfine Hamiltonian, such a realignment is only par-tial, though still useful. Experiments can be performed that explore either this spin-echo region or differentrelationships between the two magnetic field profiles, which may allow for a variety of insights about thesample surface. For example, the two field profiles can be different or the field magnitudes can be variedsimultaneously, keeping B1 = −B2. These different regimes of operation may produce different echoes,which can be collectively analyzed to provide more insight into molecule-surface interactions.Additionally, as the spatially separated wavepacket components hit the surface sequentially, rather than46simultaneously, any temporal changes in the surface that are on the time scale of the impact-time separationcan differentially impact the phases of each wavepacket component. This may result in different interferencepatterns or even loss of coherence. This loss of coherence is the basis for the sensitivity of HeSE measure-ments to surface motion [41]. Here, as in the recent experiment by Godsi et al. [44], we focus on surfaceswhose dynamics are either much faster or much slower then the molecule-surface or wavepacket-surfaceinteraction times. Note, however, that the current framework is suitable for extension to interaction regimeswhere the surface dynamics are comparable to these time scales.After leaving the last solenoid of the second branch, the wavepackets pass through another auxiliaryfield that begins with a strong dipolar field in the z′ direction. The auxiliary field then adiabatically connectsmagnetic moments aligned along the quantization axis to the radial direction of the final hexapole lens.This hexapole lens then focusses wavepackets with low-field seeking hyperfine states into the ionizationdetector and defocusses the rest. Finally, the ionization detector produces a current that is proportional tothe molecular flux into the detector port. We describe how to calculate the molecular flux that enters thedetector port in Section 3.5 and the related transfer matrix formalism in Section 3.6.Analyzing the detector current as a function of the magnetic field profiles, the apparatus geometry,and the sample orientation can provide information about the interaction of the molecules with the samplesurface. We discuss one possible analysis scheme in Section 3.7.3.3.1 Molecular Hyperfine HamiltonianIn principle, the only requirement for a molecular species to be suitable for molecular hyperfine interferome-try is that the molecule have internal degrees of freedom whose energies are magnetic-field dependent. Sucha requirement could be fulfilled by the presence of a nuclear spin, a rotational magnetic moment, or even anelectronic spin. In practice, however, if the energy dependence on the magnetic field is too weak relative tothe kinetic energy, state selection and state manipulation is difficult. On the other hand, if the dependenceis too strong, the molecules may be difficult to control. Given these restrictions, we deem molecules thathave a closed shell and are in an electronic state with zero orbital angular momentum to be most suitable formolecular hyperfine interferometry. In this case, the dominant interactions induced by magnetic fields aredue to the nuclear magnetic spins of the molecules.The hyperfine states of such a closed shell molecule with zero orbital angular momentum arise from47SampleθBeam SourceMagnetic LensDetectorMagnetic LensControlFieldsControlFields𝐳𝐱𝐳′𝐱′Auxiliary Field Auxiliary FieldFigure 3.1: A generic molecular hyperfine interferometer consists of a beam source (green), magneticlenses (dark blue), auxiliary fields (light blue), control fields (purple), the sample (hatched rect-angle) in an ultra-high vacuum chamber, and the detector (red). See Section 3.3 for more detailson each component. The arrows and dashed line indicate the direction and path of the molecularbeam, which is initially along the +x direction and then along the −x′ direction after scattering.The two branches of the apparatus are separated by an angle θ . z and z′ denote the direction ofthe quantization axes before and after scattering, respectively. This definition of the quantizationaxes has been chosen to match the experiment by Godsi et al. [44] and to simplify rotating thequantization axes in the transfer matrix method. The y and y′ axes are identical and point into thepage.coupling between the nuclear spin and the rotational degrees of freedom. Interactions of these hyperfinestates with a magnetic field arise from the response of the nuclear and rotational magnetic moments to theexternal magnetic field. We assume that the hyperfine Hamiltonian, also referred to here as the RamseyHamiltonian [156], is of the following form:HˆR(~B) =U(Iˆ2, Jˆ2, Iˆ · Jˆ, I,J)+V(Iˆ2, Iˆ ·~B, I,~B2)+Q(Jˆ2, Jˆ ·~B,J,~B2), (3.1)where ~B is the vector of the external magnetic field, assumed to be uniform across the molecule; Iˆ and Jˆare the nuclear spin and rotational angular momentum operators, respectively; I and J are the nuclear spinand rotational angular momentum quantum numbers, respectively; U contains all spin-rotational couplings(such as Iˆ · Jˆ or Iˆ2Jˆ2); V contains all interactions of the nuclear spins with the magnetic field (such as Iˆ ·~B);and Q contains all interactions of the rotational angular momentum with the magnetic field (such as Jˆ ·~B).Both V and Q are assumed to be proportional to positive powers of |~B|.At large magnetic fields, V and Q dominate, making the eigenbasis |ImIJmJ〉, where mI/J is the projec-48tion of the angular momentum ~I/~J onto the external magnetic field direction. At zero field, HˆR is diago-nalized by |IJFM〉, where Fˆ = Iˆ+ Jˆ is the total angular momentum operator and M is the projection of ~Fonto a chosen quantization axis. At intermediate fields, the eigenbasis is a function of the magnetic field andcan be represented as a superposition of either |IJFM〉 or |ImIJmJ〉 states. Note that M is a good quantumnumber at all field strengths. We call an eigenstate of HˆR a Ramsey state, which we denote as |R〉 and whichhas the energy ER. The number of eigenstates of HˆR is NR, such that 1≤ R≤ NR.We treat the apparatus as a one dimensional system and account for the actual geometry by rotating thebasis of the hyperfine states at the appropriate locations (see Section 3.6.2). The total Hamiltonian can thusbe written asHˆ(x) =pˆ22m+ HˆR(~B(x))(3.2)where pˆ is the centre of mass momentum operator, m is the molecular mass, and x is the position of themolecule in the apparatus. The magnetic field ~B(x) is now spatially dependent, reflecting the magneticfield profiles of the two branches of the apparatus. In principle, the total Hamiltonian should incorporatemolecule-surface interaction terms. However, we include these molecule-surface interactions through theuse of a scattering transfer matrix (see Section 3.6.3). This allows for the coherent propagation of themolecular wavepacket through the apparatus to be determined analytically and for the molecule-surfaceinteractions to be described phenomenologically.The system eigenstates |ER〉 are defined by the total Hamiltonian (3.2) through Hˆ |ER〉= E |ER〉. Notethat the system eigenstate |ER〉 is NR degenerate and that any linear combination of these states with thesame label E is also an eigenstate of Hˆ. This degeneracy occurs as, while the NR different Ramsey statesmay have different energies, the kinetic energy can always be selected to maintain the same total energy. Forthe sake of convenience, we choose the orthonormal basis to be that defined by HˆR(~B(x))|ER〉= ER |ER〉for x ≤ 0−. The zero of x is defined to be immediately after the magnetic lens, while y± ≡ limδ→y± δ . Weuse these limit definitions as we will deal with discontinuities in the magnetic field when working with thetransfer matrix formalism (Section 3.6). As an example of the use of this notation, the statement that bothone-sided limits are equal at the point x(i.e. lima→x− f (a) = limb→x+ f (b))can be written as f (x−)= f (x+).The above definition of |ER〉 produces, for all x, a unique labelling of the system eigenstate |ER〉 by49the total energy E and the internal state R, where R is a Ramsey state in the high magnetic field locatedimmediately after the magnetic lens (i.e. at x= 0−). Note that, because of this definition, HˆR(~B(x))|ER〉 6=ER |ER〉 for x≥ 0+; that is, the system eigenstates are superpositions of the local Ramsey states for x≥ 0+.This unique labelling of the system eigenstates is valid for all x as the eigenstate wavefunctions have awell-defined phase relationship throughout the entire apparatus. See Section 3.6.1 for more details on thespecifics of this phase relationship.3.4 Impact of the Magnetic Lens on the Molecular StatesThe magnetic lenses are designed to focus molecules with certain hyperfine states either onto the sample orinto the detector. The remaining molecules are either defocussed or insufficiently focussed and contributesignificantly less to the experimental signal. Roughly, high-field seeking states are defocussed, some of thelow-field seeking states are well focussed and the rest of the low-field seeking states are partially focussed.The actual proportions of each hyperfine state in the molecular beam must be measured or calculated fromsimulation. These magnetic lenses typically use large magnetic fields and large magnetic field gradients toperform this focussing [153, 154].In general, magnetic lenses may take different forms, but we will consider lenses that have one keyfeature: the internal degrees of freedom of the outgoing molecular wavepackets are decohered in the high-magnetic field basis (i.e. |ImIJmJ〉). More precisely, we assume that the wavepacket exiting the magneticlens is a mixed state of the form:ρ0 =∑R0PR0 |ΨR0k0〉〈ΨR0k0 | , (3.3)where|ΨR0k0〉=∫dr ψR0k0(r) |rR0〉 ; (3.4)ψR0k0(r) ≡ 〈r|ΨR0k0〉 is the wavefunction of a molecule in state |R0〉; ρ0 is the initial (time t = 0) densitymatrix; |rR0〉 ≡ |r〉 |R0〉; |r〉 is an eigenstate of the position operator; |R0〉 is an eigenstate of HˆR(~Blens); ~Blensis a high magnitude, z-aligned magnetic field; k0 is the experimentally-determined mean wavenumber of the50wavepacket; and PR0 is the probability that the hyperfine statevector of the molecule is |R0〉. Note that ρ0is diagonal in |R0〉 but not in |r〉 (or |k〉, the momentum basis). Also, ~Blens = ~B(x = 0−)corresponds to thefinal portion of the auxiliary field (i.e. a strong, z aligned, dipolar field), not the field inside the hexapolemagnet itself (see Section 3.3).That such a form of the wavepacket is valid follows from the work by Utz et al. [157]. The authorsshow that the two wavepackets arising from a spin– 12 particle passing through a Stern-Gerlach apparatusare quickly decohered with respect to one another, even before they separate spatially. That is, the quantumdynamics themselves cause decoherence between the spin degrees of freedom (but not the spatial); a mea-surement or coupling to an external bath is not required. This decoherence occurs as the large magnetic fieldgradients cause a rapid oscillation in the off-diagonal terms of the extended Wigner distribution. That is, thephase relationship between the spin-up and spin-down components oscillates heavily in both the positionand momentum bases, destroying coherence.Given that the magnetic lenses we consider act like a Stern-Gerlach apparatus for the molecular hyperfinestates, it is reasonable to assume that the internal hyperfine degrees of freedom will also decohere. Thus, weneed only determine the values of PR0 for a specific magnetic lens. These can be found via semi-classicalcalculations [44, 158], may be measured experimentally [158] or may potentially be determined by solvingthe full 3D Schro¨dinger equation within the lens.The mean velocity v0 and velocity spread σv of the molecules in the molecular beam can be measuredexperimentally [44]. Both of these values are determined from the position and profile of scattering peaksobtained from the scattering of the probe molecules by appropriate sample surfaces [44]. We assume that theinitial wavefunction of a molecule ψR0k0(r) is Gaussian and is characterized by k0 ≡mv0/h¯ and σk ≡mσv/h¯,where m is the mass of the molecule. More precisely,ψR0k0(r) =∫dk1(2piσ2k) 14e−(k−kR00)24σ2keikr√2pi=√σk(2pi) 14eikR00 re−r2σ2k (3.5)where kR00 is taken to be k0. Though kR00 may in fact depend slightly (on the order of ppm) on R0, we show51later that the experimental signal is insensitive to small changes in kR00 .3.5 Wavepacket Propagation and Signal CalculationThe primary measured value of the experiment is a current that is proportional to the molecular flux en-tering the detector. This measured current is a function of the magnetic fields, the scattering geometry,and the surface properties. The molecular flux entering the detector can be calculated as the product ofthe molecular flux incident to the apparatus and the probability that a molecule entering the apparatus willsuccessfully pass through the apparatus and be detected. It is this probability of detection Pdetection that issensitive to the experimental parameters and surface properties. Note that the incident molecular flux couldbe either continuous or pulsed, as long as the density is low enough that the molecules can be considerednon-interacting.As the detector has a finite time-response, the probability of detection is given byPdetection =1τ∫ t2t1dt〈Cˆ(t)〉, (3.6)where t1 and t2 are the initial and final times of the detection window τ = t2−t1, and 〈Cˆ(t)〉 is the expectationvalue of the detector measurement operator Cˆ. This expectation value is given by〈Cˆ(t)〉= Tr(ρˆ(t)Cˆ), (3.7)where ρˆ(t) ≡ Uˆρ0Uˆ† = ∑R0 PR0 |ΨR0k0(t)〉〈ΨR0k0(t)| is the time evolved density matrix, Uˆ ≡ e−iHˆh¯ t is thetime evolution operator, ρ0 is the density matrix (3.3) at t = 0, and |ΨR0k0(t)〉 ≡ Uˆ |ΨR0k0〉.Given that the detector consists of a magnetic lens that focusses molecules with particular states into ameasuring apparatus, such as an ionization detector [44], and that the internal degrees of freedom of thesemolecules are decohered by the second magnetic lens (see Section 3.4), we can model the detector with adiagonal operatorCˆ =∑RD∫dx cRD(x) |xRD〉〈xRD| (3.8)The matrix elements of Cˆ are the probabilities cRD(x) of detecting, at position x, a molecule whose internal52state is a high-field eigenstate |RD〉 of HˆR. Note that cRD(x) = 0 for x < xD, the detector position.Using the time evolution operator, we determine the time-dependence of the density matrix ρ(t) to beρ(t) = ∑R0RR′∫dE∫dE ′ PR0e− ih¯ (E−E ′)tαERk0R0α∗E ′R′k0R0 |ER〉〈E ′R′| , (3.9)where αERk0R0 ≡∫dr ψR0k0(r)Φ∗ERR0 (r) is the overlap between the initial wavefunction ψR0k0(r) and the systemeigenstate wavefunction ΦERR0 (r)≡ 〈rR0|ER〉.We can evaluate 〈Cˆ(t)〉 by inserting a resolution of the identity∑RD∫dr |rRD〉〈rRD|, where HˆR(~B(x)) |RD〉=ERD |RD〉 for x≥ x+D and xD is the starting location of the detector (see Figure 3.2). In other words, |RD〉 is aRamsey state in the strong dipolar magnetic field of the detector auxiliary field. The result is〈Cˆ(t)〉= ∑RD,R′D∫dr∫dr′ 〈r′R′D|ρ(t) |rRD〉〈rRD|Cˆ |r′R′D〉 (3.10)where we have evaluated the trace in the |r′R′D〉 basis and〈r′R′D|ρ(t) |rRD〉= ∑R0RR′∫dE∫dE ′ PR0e− ih¯ (E−E ′)tαERk0R0α∗E ′R′k0R0 ΦERR′D(r′)Φ∗E′R′RD (r). (3.11)We emphasize that R and RD are indices of different sets of Ramsey states, i.e. 〈R|RD〉 6= δRRD , unless themagnetic fields at the first magnetic lens (x = 0−) and the detector magnetic lens (x = x+D) happen to beidentical.We also have〈rRD|Cˆ |r′R′D〉=∑R′′D∫dz cR′′D(z)δ (r− z)δRDR′′Dδ (r′− z)δR′DR′′D= cRD(r)δ (r′− r)δR′DRD (3.12)which, when inserted with Eqn. (3.11) into Eqn. (3.10), results in〈Cˆ(t)〉= ∑R0RR′∫dE∫dE ′ PR0e− ih¯ (E−E ′)tαERk0R0α∗E ′R′k0R0(∑RD∫dr ΦERRD (r)Φ∗E ′R′RD (r)cRD(r)). (3.13)53The wavepacket is almost entirely confined to the region r≤ 0−, as ψR0k0(r) has a Gaussian profile (3.5)with spatial width on the order of 10 A˚ (as determined from the measured velocity distribution for oH2 [44]).Thus, we can evaluate αERk0R0 ≡∫dr ψR0k0(r)Φ∗ERR0 (r) if we know ΦERR0 (r) for r ≤ 0−. Given the definitionof the eigenstate |ER〉, discussed in Section 3.3.1, we show in Section 3.6.1 that ΦERR0 (r) = AReirkERδRR0for r ≤ 0− (cf. Eqn. (3.34)), where kER0 is defined in Eqn. (3.23). Combined with the definition (3.4) ofψR0k0(r),αERk0R0 ≈∫dr δRR0A∗R√σk(2pi) 14ei(kR00 −kER)re−r2σ2k= δRR0ΓERk0R0 (3.14)where ΓERk0R0 = A∗R(2pi)14√σke−(kER−kR00)24σ2k . Thus,〈Cˆ(t)〉=∑R0∫dE∫dE ′ PR0e− ih¯ (E−E ′)tΓER0k0R0Γ∗E ′R0k0R0(∑RD∫dr ΦER0RD (r)Φ∗E ′R0RD (r)cRD(r)), (3.15)where we have performed the sums over R and R′.If the detection window t2− t1 is large enough that the entire wavepacket passes through the detectionregion defined by cRD(z) we havePdetection =1τ∫ t2t1dt〈Cˆ(t)〉 ≈ 1τ∫ ∞−∞dt〈Cˆ(t)〉=∑R0∫dE∫dE ′ PR02pi h¯τδ(E−E ′)ΓER0k0R0Γ∗E ′R0k0R0(∑RD∫dr ΦER0RD (r)Φ∗E ′R0RD (r)cRD(r))=∑R0∫dE PR0∣∣∣ΓER0k0R0∣∣∣2(∑RD∫dr2pi h¯τ∣∣∣ΦER0RD (r)∣∣∣2 cRD(r)), (3.16)where 2pi h¯τ δ(E−E ′)= 1τ ∫ ∞−∞ dte− ih¯ (E−E ′)t .Physically, one can see that the probability of detection is proportional to the overlap∣∣∣ΓER0k0R0∣∣∣2 of theinitial wavepacket and a system eigenstate multiplied by the overlap∫dr 2pi h¯τ∣∣∣ΦER0RD (r)∣∣∣2 cRD(r) of the samesystem eigenstate and the detection region, as expected.54Substituting for∣∣∣ΓER0k0R0∣∣∣2 and given thatΦER0RD (r)≡ 〈rR0|ER〉= eikRD r 〈RD|ER0〉≡ eikRD rβER0RDfor r ≥ x+D (cf. Eqn. (3.34)), we havePdetection =∑R0PR0∣∣AR0∣∣2 ∫ dE (2pi) 12σk e−(kER0−kR00)22σ2k ∑RDcRD2pi h¯τ∣∣∣βER0RD ∣∣∣2 (3.17)where cRD ≡∫dr cRD(r) and βER0RD ≡ 〈RD|ER0〉, the projection of the system eigenstate |ER0〉 onto thedetector eigenstate |RD〉 at x+D . For the purposes of comparing to experiment, only the dependence of Pdetectionon the experimental parameters is needed, not its absolute value. Also, the value of AR0 = 1 as AReirkER ≡〈rR0|ER0〉 = eirkER(for r ≤ 0−) because of the specific definition of the system eigenstates (see Section3.3.1). Additionally, one can see that Pdetection is not sensitive to minor (on the order of ppm) changes in kR00as σk ∝ k0 in experiment [44]. Finally, in Eqn. (3.17), only βER0RD is dependent on the magnetic fields, thescattering geometry, and the surface properties. It is thus sufficient to work with the following equation:Pdetection ∝∑R0PR0∫dE e−(kER0−k0)22σ2k ∑RDcRD∣∣∣βER0RD ∣∣∣2 (3.18)To determine the values of βER0RD , we derive and apply the transfer matrix method with internal degreesof freedom (Section 3.6).3.6 Transfer Matrix Formalism with Internal Degrees of FreedomThe transfer matrix method turns the solution of the time-independent Schro¨dinger equation of a 1D systeminto a product of matrices [159]. The present problem has two unique features: (i) the propagating moleculeshave many internal degrees of freedom which may be mixed as the molecule transitions from one local fieldto another and (ii) molecules change their propagation direction after scattering by the surface. Problem55(i) is addressed in Section 3.6.1, while (ii) is addressed in Section 3.6.2. The impact of scattering on theinternal degrees of freedom is accounted for by using a scattering transfer matrix (Section 3.6.3).DetectorFirst Branch Second BranchSampleState SelectorFigure 3.2: Generic field profile of a molecular hyperfine interferometry experiment. The actual mag-netic field profiles of the experiment are approximated by N+2 regions of length Li and constantmagnetic field ~Bi (black line). The true field profile is asymptotically approached as N→ ∞. Weassume large magnetic fields in the regions of the state selector (large arrow) and the detector(eye), which, when combined with the dephasing discussed in Section 3.4, allows us to neglectpropagation in the selector and detector regions. That is, the exact locations of x0 and xD areunimportant as long as x0 is in the high-field region of the state selector, xD is in the high-fieldregion of the detector, and all propagation is treated coherently between the two points. Theinitially Gaussian wavepacket propagates from x−0 along the first branch to the sample surface(cross) at xS then, after scattering, propagates along the second branch to x+D . The two branchesare separated by an angle θ . The vertical axis indicates the magnitude of the magnetic field|~B| (the direction is not depicted for clarity), with |~B| = 0 indicated by the grey solid line. |Ri〉denotes the set of eigenstates of HˆR(~Bi), Eqn. (3.1), in each region.3.6.1 Propagation and Discontinuity MatricesWe first break up the arbitrary magnetic field profiles of the apparatus into rectangular regions of constantfield, as shown in Figure 3.2. We then solve the Schro¨dinger equation for a single eigenstate in a singleregion of constant field. Subsequently, we determine the impact of the boundary conditions that exist atthe discontinuity between two regions of constant field. Using these solutions, we determine matrices thatdescribe the spatial dependence of the eigenstate wavefunction coefficients within a region of constant field(propagation matrices) and matrices that describe how these coefficients change across the discontinuitybetween two regions of constant field (discontinuity matrices). Note that while we derive these matrices for56molecules whose internal degrees of freedom are described by the Ramsey Hamiltonian (3.1), the formalismis not limited to this Hamiltonian.Within a region of uniform magnetic field, the Ramsey Hamiltonian HˆR is constant, which allows usto derive the propagation matrix that includes the internal degrees of freedom. We begin by expanding asystem eigenstate |ER˜〉 as|ER˜〉=∑R∫dx ΦER˜R (x) |xR〉 , (3.19)where ΦER˜R (x)≡ 〈xR|ER˜〉, we define |xR〉 ≡ |x〉 |R〉, and |R〉 is one of the NR Ramsey states of a molecule insome magnetic field ~B. Note that ~B is not necessarily the local magnetic field ~Bloc of the current region andthus |R〉 is not necessarily an eigenstate of HˆR(~Bloc) at this point. Also, the eigenstates |ER˜〉 are labelledby their energy E and a particular Ramsey index R˜, such that HˆR(˜~B) |ER˜〉 = ER˜ |ER˜〉, with ˜~B an arbitrarilychosen magnetic field.Using Eqn. (3.19), the Schro¨dinger equation with the total Hamiltonian (3.2) can be shown to be (Ap-pendix B.1):− h¯22m∂ 2∂x2ΦER˜R0 (x) =ΦER˜R0 (x)E−∑RHRR0RΦER˜R (x), (3.20)where HRR0R = 〈R0| HˆR(~Bloc) |R〉. Eqn. (3.20) is in general difficult to solve because of the coupling of theinternal degrees of freedom by HˆR(~Bloc). However, if we choose the eigenbasis of the internal degreesof freedom to satisfy HˆR(~Bloc) |R〉 = ER |R〉 (that is, |R〉 is now a Ramsey state of a molecule in the localmagnetic field ~Bloc), the equations decouple and we obtain∂ 2∂x2ΦER˜R (x) =−2mh¯2(E−ER)ΦER˜R (x). (3.21)The solution isΦER˜R (x) = AReikRx+BRe−ikRx, (3.22)57where AR and BR are R-dependent coefficients andkR ≡√2m(E−ER)h¯. (3.23)As per the single channel transfer matrix method [159], given that ΦER˜R (x+∆x) = AReikRxeikR∆x +BRe−ikRxe−ikR∆x, we can collect the AR and BR coefficients into a 2NR-dimensional coefficient vector ~φx =(A1,A2, ...,ANR,B1,B2, ...,BNR)T and write~φx2 =Πx2−x1~φx1 , (3.24)where Πx is the 2NR×2NR propagation matrixΠx ≡[⊕ReikRx]⊕[⊕Re−ikRx]=eik1x. . . 0eikNR xe−ik1x0. . .e−ikNR x, (3.25)where ⊕ denotes the direct sum.Following the derivation of Ref. [159], we can determine how the coefficients transform across a stepdiscontinuity in the magnetic field. Using the propagation matrix (3.25) and a relabelling of the coordinatesystem, we can always set the discontinuity to appear at x = 0. Given that Eqn. (3.20) applies everywhere,the coefficients ΦER˜R (x) and their derivatives are continuous across the discontinuity (i.e. ΦER˜R (x) ∈C1(x)),for each value of R. However, the coefficients are only known when |R〉 is an eigenstate of HˆR(~Bloc), whichdiffers on each side of the discontinuity (that is, ~B(0−) 6= ~B(0+)). Note that the wavevector |ER˜〉 is the sameeverywhere in the system. Thus, by writing the wavevector |ER˜〉 in the two different bases corresponding tothe eigenstates of HˆR on each side of the field, we see that the coefficients at a specific value of x are related58by a basis transformation:|ER˜−〉= |ER˜+〉∑R−∫dx ΦER˜R− (x) |xR−〉=∑R+∫dx ΦER˜R+ (x) |xR+〉∑R−R+∫dx ΦER˜R− (x)〈R+|R−〉 |xR+〉=∑R+∫dx ΦER˜R+ (x) |xR+〉=⇒ ΦER˜R+ (x) =∑R−ΦER˜R− (x)〈R+|R−〉 , (3.26)where |ER˜±〉 is the wavevector written in the basis of |R±〉, |R±〉 are the eigenstates of HˆR(~B(0±)) on theleft (−) and right (+) sides of the discontinuity at x= 0, respectively, and ∑R+ |R+〉〈R+| was inserted in thethird line (recall that |xR−〉 ≡ |x〉 |R−〉). The values 〈R−|R+〉 are recognized as the matrix elements SR−R+ ofthe matrix SR+R− whose columns are the eigenstates of HˆR(~B(0+)) written in the |R−〉 basis.Since ΦER˜R (x) ∈C1(x) for each value of R separately, we can equate the two limits limx→0∓ΦER˜R+ (x) andthe two limits of the derivative limx→0∓ ∂∂xΦER˜R+ (x). Solving the resultant equations for the coefficients AR+and BR+ , we obtain (Appendix B.2):AR+ =∑R−S∗R−R+∆+R+R−AR−+∑R−S∗R−R+∆−R+R−BR− (3.27)BR+ =∑R−S∗R−R+∆−R+R−AR−+∑R−S∗R−R+∆+R+R−BR− (3.28)where S∗R−R+ ≡ 〈R+|R−〉, ∆±R+R− ≡ 12(1± kR−kR+), kR± ≡√2m(E−ER±)h¯ , and ER± ≡ 〈R±|HˆR(~B(0±))|R±〉. Thereare NR such sets of equations, one for each value of R+. Working again with ~φx =(A1,A2, ...,ANR ,B1,B2, ...,BNR)T ,one can write the matrix equation~φx+ = K~φx− , (3.29)where x∓ indicates the location just before (−) or just after (+) the discontinuity located at x and K is the592NR×2NR discontinuity matrixK≡SR+R−† ◦∆+ SR+R−† ◦∆−SR+R−† ◦∆− SR+R−† ◦∆+ (3.30)where ◦ denotes the element-wise Hadamard product, such that (SR+R−†◦∆±)R+R− ≡ S∗R−R+∆±R+R− . This matrixallows one to calculate the coefficients of the wavefunction as one moves from one region of constantmagnetic field to another through a discontinuity. Thus, if one breaks up any magnetic field profile into aseries of constant regions separated by discontinuities, one can systematically approach a perfect descriptionof the propagation of a molecule with internal degrees of freedom through a magnetic field of arbitraryprofile through repeated application of K and Πx. Furthermore, this approach is not restricted to moleculesmoving through magnetic fields. Many other types of quantum objects moving in a single dimension withinternal degrees of freedom that couple to an external static potential can also be analyzed in this way.The above analysis indicates that one needs to keep track of 2NR components to build up the eigenstatesof the system exactly. However, for the current application in mind, one only needs NR components as themagnetic fields typically change the linear molecular momentum by such a small amount that the amplitudesBR of the reflected part of the wavefunction are negligible.For example, a typical velocity of the oH2 molecules in the experiment of Ref. [44] is vH2 = 1450m/s.This corresponds to the kinetic energy EH2 =12 mH2v2H2 = 5.31× 109 kHz. Ramsey [156] indicates that themaximum energy change for the hyperfine states of oH2 at 500G is approximately -2550 kHz. The experi-ment of Ref. [44] has magnetic fields up to about 1000G. For such fields, the energy changes are approxi-mately linear, so we expect the maximum change in energy to be ∆E ≈−5100 KHz. In the field-free regionbefore the discontinuity, kR− ≈mH2vH2/h¯ and after the discontinuity in the field, kR+ ≈√2mH2(E−∆E)/h¯,as per Eqn. (3.23). Then, |∆−R+R− | ≈ 2.4× 10−7 and |∆+R+R− | ≈ 1, making K approximately diagonal andillustrating the decoupling of the forward and backward channels under typical experimental conditions.We thus only need to keep track of the AR components and can define a new coefficient vector~ψx ≡(A1,A2, ...,ANR)T. (3.31)60The corresponding NR×NR propagation Px and discontinuity D matrices arePx ≡⊕ReikRx=eik1x 0. . .0 eikNR x (3.32)D≡ SR+R−† ◦∆+≈ SR+R−†, (3.33)where the matrix elements of SR+R−†are S∗R−R+ ≡ 〈R+|R−〉, HˆR(~B(0±)) |R±〉 = ER± |R±〉, 0± indicates theposition just to the left (−) or right (+) of the discontinuity, and kR is defined as in Eqn. (3.23). Specifically,D changes the basis of the vector ~ψx from |R−〉 to |R+〉. That is, ~ψx is always in the eigenbasis of HˆR(~B(x)).Finally, given that BR ≈ 0, the eigenstate coefficients are nowΦER˜R (x) = AReikRx. (3.34)Given that a generic transfer matrix M has the property MσzM† =σz[159], the decoupling of the forwardand backward channels implies that the forward channel matrix MF (composed of a product of Px and Dmatrices) is now unitary.3.6.2 Rotation MatricesScattering off of the sample surface changes both the propagation direction and the internal states of themolecule. To take into account the change in the direction of the propagation path when applying the transfermatrix formalism, we need only change the orientation of the quantization axis. To do so coherently, weapply an NR×NR rotation matrix R(φ ,Θ,χ) to ~ψx, where φ , Θ, and χ are the Euler angles in the ZY Zconvention. In particular, given the definition of the axes shown in Figure 3.1, a passive rotation of the statevector about the y axis by the angle θ is required to account for a change of angle θ in the propagation61direction. This passive rotation modifies the basis of ~ψx, but leaves the physical state unchanged. Weperform this rotation by applying the equivalent active rotation of angle −θ to ~ψx.The matrix elements of R(0,−θ ,0), when written in the |R〉 eigenbasis of HˆR(~Bloc) and where ~Bloc isthe local magnetic field, are:RR(0,−θ ,0)≡[〈R′| Rˆ(0,−θ ,0) |R〉]=[∑FMF ′M′〈R′|F ′M′〉〈F ′M′|Rˆ(0,−θ ,0)|FM〉〈FM|R〉]= SRFM†RFM(0,−θ ,0)SRFM, (3.35)where |FM〉 ≡ |IJFM〉 is an angular momentum state with total angular momentum F , z axis projectionM, total nuclear spin angular momentum I and total rotational angular momentum J; the subscripts of RRand RFM denote the basis of the matrix representation, |R〉 and |FM〉, respectively; SRFM is the matrix whosecolumns are the eigenstates |R〉 written in the |FM〉 basis, Rˆ(α,β ,γ) is the rotation operator (also in the ZY Zconvention) andRFM(0,β ,0) =[δFF ′DFM′M(0,β ,0)]=[δFF ′dFM′M(β )], (3.36)where DFM′M(α,β ,γ) are the Wigner D-matrices and dFM′M(β ) are the Wigner small d-matrices [160]. Notethat RFM(0,−θ ,0) is diagonal in F , because of conservation of angular momentum, but not diagonal in M.Thus, one must be careful to also perform a passive rotation on the local magnetic field vector if rotations areperformed in a region with non-zero field. Typically, however, the sample chamber is magnetically shielded.3.6.3 Scattering Transfer MatricesScattering by the sample surface can involve many complex phenomena that may change the internal state,the momentum, and the total energy of the scattering molecule. For this thesis, we focus on scattering pro-cesses that conserve the total energy of the molecules. The interactions of the molecules with the samplesurface can be phenomenologically described with the total scattering transfer matrix. This matrix is the622NR× 2NR matrix Σ˜ that relates the wavefunctions on the “left” side of the scattering event to those onthe “right” (as opposed to the scattering matrix, which relates the incoming wavefunctions to the outgo-ing). However, because the initial wavepacket (3.5) does not contain any negative momentum states, themagnetic fields of the solenoids do not cause significant backscattering (Section 3.6.1), and the detectoronly detects molecular flux in the forward scattering direction, we need only work with the NR×NR matrixΣ ≡ PfwdΣ˜P†fwd, where Pfwd is an NR×2NR projection matrix onto the forward scattering states. We defineΣ in the |ImIJmJ〉 basis, where the |ImIJmJ〉 states are themselves defined with respect to the quantizationaxis that is normal to the surface sample. We choose this basis to relate to scattering calculations, which arefrequently carried out in the |JmJ〉 basis with a quantization axis normal to the sample surface.In general, the scattering transfer matrix elements ΣImIJmJI′m′IJ′m′J are functions of the incident energy E,the outgoing energy E ′, the incident momentum~k, and the outgoing momentum~k′. As we are restrictingourselves to iso-energetic processes, E = E ′. Also, Eqn. (3.23) defines the magnitudes of the momentumbefore and after the scattering event. This leaves the scattering transfer matrix elements as functions of onlyenergy and the four angles that define the scattering geometry. These angles are: (1) the angle between thetwo branches, (2) the angle between the surface normal and the scattering plane, (3) the angle between thefirst branch and the projection of the surface normal on the scattering plane, and (4) the azimuthal angle ofthe sample. The scattering plane is the plane defined by the two branches of the apparatus.Given that the experiment only probes a single scattering direction at a time, the scattering transfermatrix will not, in general, be unitary. This incorporates state-dependent loss channels into the formalism.Additionally, the scattering transfer matrix is, in general, time-dependent. Here, we assume that the time-scales of the surface dynamics are significantly different from the molecule-surface or wavepacket-surfaceinteraction time-scales and assume Σ to be time-independent.Because Σ is defined with respect to the surface normal, we use rotation matrices to appropriately changethe basis of ~ψ before and after applying the scattering transfer matrix. We ensure that the total rotationcorresponds to the change in propagation direction induced by scattering off of the sample surface and thatthe quantization axis is again coplanar with the two branches of the apparatus.The scattering transfer matrix elements for a specific molecule-surface interaction can be determinedfrom scattering calculations [44, 161]. Alternatively, they can be treated as free parameters and determined63from the experimental measurements by solving the inverse scattering problem. Such a problem can poten-tially be solved efficiently using machine learning based on Bayesian optimization [45, 64].3.6.4 Calculation of Eigenstate CoefficientsTo determine the dependence of the probability of detection (3.18) on the magnetic fields and the surfaceproperties, we must determine the coefficients βER0RD . This can be done by multiplying the initial coefficientvector ~ψER0x0 (3.31) of a system eigenstate |ER0〉 by a succession of transfer matrices to obtain the finalcoefficient vector ~ψER0xD ≡(βER01 ,βER02 , · · · ,βER0NR)T:~ψER0xD = SRDRN†M2MΣM1~ψER0x0 , (3.37)where M1/2 describes the propagation through the first/second branch of the apparatus, MΣ describes thescattering, and SRDRN†changes the basis of the coefficient vector to the eigenbasis |RD〉 of HˆR(~B(x+D)) at thelocation of the detector xD. The M matrices are defined asM1 = PLnSRnRn−1† · · ·PL2SR2R1†PL1SR1Rini†(3.38)MΣ = SRnFM†RFM(0,−θ2,0)ΣFMRFM(0,θ1,0)SFMRn†(3.39)M2 = PLN SRNRN−1† · · ·PLn+1SRn+1Rn†(3.40)where Rini refers to the eigenbasis |Rini〉 of HˆR(~B(0−)) at the initial location of the wavepacket; Ri refers tothe eigenbasis |Ri〉 of HˆR(~Bi) in region i of the apparatus, as depicted in Figure 3.2; FM refers to the |IJFM〉basis where ~F ≡~I+ ~J and M is the projection on the local z axis; θ1 is the angle between the z axis and thesample surface normal (see Figure 3.1); θ2 is the angle between the sample surface normal and the z′ axis;Li is the length of region i, as depicted in Figure 3.2; N is the total number of regions between x0 and xD (seeFigure 3.2); n is the number of regions between the initial position of the wavepacket x0 = 0 and the sampleposition xS; ΣFM ≡ SFMRIJ†ΣSRIJFM†is the scattering transfer matrix written in the |IJFM〉 basis; RIJ ≡ ImIJmJ;and Σ is the scattering transfer matrix in the |ImIJmJ〉 basis. Note that θ1 +(−θ2) ≡ −θ , where θ is theangle between the branches of the apparatus (see Figure 3.1). Note also that Eqn. (3.39) assumes that thesurface normal is in the scattering plane; additional rotation matrices are needed if this is not the case.64By defining a matrix ΨExi ≡(~ψE1x0 , ~ψE2x0 , · · · , ~ψENRx0), all NR×NR coefficients βER0RD can be simultaneouslyobtained fromΨExD = SRDRN†M2MΣM1ΨEx0= SRDRN†M2MΣM11NR , (3.41)where ΨEx0 ≡ 1NR because of the specific definition of the system eigenstates (see Section 3.3.1). UsingEqns. (3.38–3.41), we can obtain βER0RD , and thus Pdetection (3.18), as functions of the magnetic field profile,the scattering matrix elements, and the scattering geometry.3.7 Application to ortho-HydrogenThe theoretical framework described in Sections 3.3–3.6 connects the scattering transfer matrix elementsΣImIJmJI′m′IJ′m′J to the experimentally observed signal, which is proportional to Pdetection (3.18). By changingthe magnetic field profiles in the two arms of the apparatus, one can obtain information about how thescattering affects various hyperfine states. To illustrate our theoretical framework and to demonstrate theimpact of the scattering transfer matrix on the experimentally observed signal, we consider a beam of oH2 inits ground state and a simplified apparatus that contains only a few regions of constant magnetic field, asdepicted in Figure 3.3.3.7.1 Ortho-Hydrogen Hyperfine HamiltonianThe Hamiltonian describing the relevant internal degrees of freedom of oH2 is [156]:HˆRoH2(~B)h=−α Iˆ ·~B−β Jˆ ·~B− cIˆ · Jˆ+ 5d(2J−1)(2J+3) [3(Iˆ · Jˆ)2+32Iˆ · Jˆ− Iˆ2Jˆ2] (3.42)where, for simplicity, we have neglected magnetic shielding of the nuclear and rotational magnetic momentsby the molecule and diamagnetic interactions of the molecule with the magnetic field; ~B is the local magneticfield; Iˆ is the nuclear spin operator; Jˆ is the rotational angular momentum operator; α ≡ µIhI ≈ 4.258kHz;β ≡ µJhJ ≈ 0.6717kHz; c ≈ 113.8kHz; d ≈ 57.68kHz; I = 1 is the total nuclear spin angular momentum inunits of h¯; J = 1 is the total rotational angular momentum in units of h¯; µI is the nuclear magnetic moment65𝜃 =𝜋4𝐁1 𝐁2𝐳𝐱Variable0G 0G 0G 0G𝐳′𝐱′DetectorFirst Branch Second BranchSampleVariable1000GState SelectorField Magnitude: 1000G1m20cm 30cm 30 cm 20cm1mRegion Length:Figure 3.3: A magnetic field profile that approximates the true magnetic field profile of an experimentusing oH2. We combine this approximate field profile with the transfer matrix formalism tocalculate the observed signal. Bi refers to the different magnetic field vectors of the controlfields. z′ and x′ refer to the new coordinate system defined to align with the second branch ofthe apparatus (see Figure 3.1). The sample is located at the cross in the centre of the diagram.The surface normal of the sample is assumed to bisect the angle between the two branches of theapparatus. The propagation direction is x before scattering and −x′ after scattering. The anglebetween x and −x′ (i.e. the angle between the two arms of the apparatus) is θ = 45◦. B1 isdirected along x and B2 is directed along −x′, as per the arrows. The fields just after the stateselector and just before the detector are directed toward the z and z′ directions, respectively and asper the arrows. Additional computational parameters not shown above can be found in AppendixB.3.of a single nucleus; and µJ is the magnetic moment due to molecular rotation. The first two terms describethe interaction of the nuclear and rotational magnetic moments with the external magnetic field, the thirdterm describes the nuclear spin-rotational magnetic interaction [156, 162, 163], and the terms proportionalto d describe the magnetic spin-spin interaction of the two nuclei [156, 162, 163].3.7.2 Experiment and ObservablesWhile there are many possible experimental protocols, we focus on the full interferometer mode used byGodsi et al. [44]. The experiment is performed by initiating a continuous flux of oH2 molecules throughthe apparatus and measuring the current of the ionization detector while varying the first and second controlfields (B1 and B2 in Figure 3.3).In particular, B1 is set to a specific value while B2 is varied around the point−B1 (i.e. about the spin-echocondition). In principle, B2 could also be set to vary around +B1, where spin echoes have also been observed66TheorySingle Velocity430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.) (a)0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) (d)TheoryIntegrated Over Velocity430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.) (b)0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) (e)Experiment430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Measured Signal (a.u.) (c)0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) (f)Figure 3.4: Upper Panels: Calculated and experimental signals close to the spin echo condition versusthe magnetic field of the second coil |B2|. Lower Panels: Fourier amplitudes of the upper panelsversus the generalized gyromagnetic ratio γ . For panels (a), (b), (d), and (e), the field profile isdepicted in Figure 3.3; B1 = 440gauss; the scattering transfer matrix Σ = 19 and is constant forall energies; and the signal is sampled at a rate of 300 points per 20 gauss. Panels (a) and (d) onlyinclude a single velocity (or, equivalently, a single value of energy) in Eqn. (3.18) while panels(b) and (e) include the full integral. For the experimental data shown in panel (c), B1 = 437gauss,the sample was the (111) surface of Cu and the signal was sampled every 0.065 gauss (a samplingrate of ∼308 points per 20 gauss). Panel (f) shows data for oH2 scattering off of Cu(111) (bluefull circles) and Cu(115) (red open circles). All experimental data was obtained from Godsi etal. [44].[155], but we choose to vary B2 about −B1 to match the relevant experiment by Godsi et al. [44]. Thisvariation of the magnetic fields results in oscillatory curves of the detector current versus B2, as shown inFigure 3.4 (a–c). These oscillations reflect the interference pattern that occurs when the various wavepacketsrecombine after passing through the final control field (see Section 3.3). This interference pattern containsinformation about how the individual hyperfine states of the molecule interact with the sample surface.The x directed magnetic fields of a solenoid changes the energies of all of the NR = 9 hyperfine statesand induces all(NR2)= 36 possible transitions. The frequencies of these transitions depend on the magnitudeof the magnetic fields. By changing the magnitude of the second magnetic field, we are able to probethe rates of change of these transition frequencies with the magnetic field: the (generalized) gyromagnetic67ratios γi j(B) =∣∣∣dωi j(B)dB ∣∣∣, where ωi j ≡ 1h∆Ei j = Ei−E jh , and Ei is the energy of Ramsey state i [44]. The Fouriertransforms of the oscillatory curves that give these gyromagnetic ratios are shown in Figure 3.4 (d–f). Toobtain these results, we assumed that the surface normal of the sample bisects the angle defined by the twobranches and that the scattering transfer matrix is the identity matrix and is independent of energy, i.e. weassume for the present calculation that the only impact of scattering is the change of propagation direction,as modelled with rotation matrices (Section 3.6.2).TheorySingle Velocity0 200 400 600B1 (gauss)246810 (kHz/gauss)10 310 210 1100TheoryMany Velocities0 200 400 600B1 (gauss)246810 (kHz/gauss)10 310 210 1100Experiment0 200 400 600B1 (gauss)246810 (kHz/gauss)10 810 610 410 2100Figure 3.5: 2D Fourier amplitude plots formed by the concatenation of spectra plots (such as Figure3.4 (d-f)) for various values of the magnetic field of the first solenoid B1. Colour indicates theFourier amplitude. For the theory plots, the field profile is depicted in Figure 3.3; the scatteringtransfer matrix Σ = 19 and is constant for all energies; B2 was varied from −(B1− 10gauss) to−(B1 +10gauss); the signal was sampled at a rate of 300 points per 20 gauss; and all data witha value less than 10−3.5 has been replaced with 10−3.5 for clarity. For the experimental plot, thesample was the (111) surface of Cu; all data with a value of less then 10−8 has been replaced with10−8 for clarity; the dashed lines indicate transitions identified by Godsi et al. [44]; and the datawas obtained from Godsi et al. [44].The location of each feature in the spectra is reflective of a gyromagnetic ratio and is independent ofthe molecule-surface interactions, being only a function of the hyperfine energy level structure of oH2. Therelative height of each feature, however, is dependent on the molecule-surface interactions, as exemplifiedin the experimental spectrum shown in Figure 3.4 (f). From Figure 3.4, one can see that integrating over thevelocity distribution is important to produce the spin-echo effect and to bring the observed signal closer toexperiment.A different spectrum can be obtained for every possible value of B1 and then combined to form a 2Dmap of the generalized gyromagnetic ratios and their contributing amplitudes as a function of B1, as shownin Figure 3.5. This protocol is equivalent to observing the scattering of molecules with different internal68hyperfine states as different values of the magnetic field in the first branch produce different superpositionsof the hyperfine states. One can clearly see both the magnetic field-dependence of the gyromagnetic ratios,the impact of integrating over the velocity distribution, and the stark similarities and differences between theexperimental and theory plots.These differences (also seen in Figure 3.4) occur as we have not taken the details of the molecule-surfaceinteractions into account in the theoretical calculations. That is, the identity scattering transfer matrix thatwe use to calculate the theoretical plots does not well-represent the molecule-surface interactions presentin the experiment. A scattering transfer matrix that well-represents the molecule-surface interactions stillneeds to be determined; one approach to determining this scattering transfer matrix is discussed in Chapter 4.Furthermore, differences between the theoretical and experimental plots also arise because the field profileused for the theoretical calculations (Figure 3.3) only approximates the real magnetic field profile; we cansystematically reduce the differences arising from this approximation by simply increasing the number offield regions in the field profile used, however.We now examine the sensitivity of the calculated signals to various changes in the scattering transfermatrices. Figures 3.6 and 3.7 demonstrate the impact of random variations of the scattering transfer matrixΣ on the oscillatory plots (for B1 = 440gauss) and their spectra, respectively. For simplicity, we keep thematrix elements of Σ independent of energy.The first row of each figure (labelled RP) reflects the impact of differing phases imparted to each hy-perfine state after scattering. Specifically, Σ =⊕9i=1 eiθi is a diagonal unitary matrix whose nine phases θiare randomly chosen from a uniform distribution of width 2pi . Such a form of scattering would result frompurely elastic scattering where the different hyperfine states probe the surface for different lengths of time(i.e. each state penetrates to a different depth or encounters a resonance with a different lifetime). Significantdifferences in the relative peak amplitudes can already be seen at this point, indicating that the calculatedsignal is sensitive to these phases.The second row of each figure (labelled RDA) reflects the impact of differing state losses due to scat-tering. Specifically, Σ =⊕9i=1 Ai is a diagonal matrix whose diagonal elements are randomly chosen froma uniform distribution on the interval [0,1). This form models the impact of different losses of each hyper-fine state to different scattering directions, reactions with the surface, or adsorbtion to the surface. Again,69significant changes are observed, indicating sensitivity to these features.The third and fourth rows (labelled ROM and RUM, respectively) probe the impact of inelastic (pro-jection mImJ-changing) scattering on the calculated signal. For the third row, Σ is an orthogonal matrixrandomly drawn according to the Haar measure on O(9), while Σ is a unitary matrix randomly drawn ac-cording to the Haar measure on U(9) for the fourth row of each figure. Here, randomly drawing accordingto the Haar measure can be understood as analogous to drawing from the “uniform distribution” over thespace of possible matrices [164]. The orthogonal matrices model inelastic scattering where no relative phasechanges occur, while the unitary matrices model inelastic scattering where relative phase changes do occur.In both cases, there is no loss of total population during scattering. Clearly, the calculated signals are alsosensitive to inelastic scattering events, both with and without relative phase changes.3.8 Comparison with a Semi-Classical MethodThe present approach is fully quantum mechanical, while Godsi et al. [44] have described a semi-classicalmethod for calculating Pdetection that they used to model the propagation of oH2 in their molecular hyperfineinterferometer. This semi-classical method treats the internal degrees of freedom of the molecules quantummechanically and the centre-of-mass motion classically. As a result, the momentum changes induced bythe magnetic field are ignored and every internal state is described as propagating at the initial velocity v0of the molecule. The internal degrees of freedom are treated by applying the time evolution operator forthe time period ti ≡ Liv0 spent in each magnetic field of length Li. That is, the propagation is calculated inthe molecular reference frame with a time-dependent Hamiltonian. Here, we compare the results of the twoapproaches.As per the application of the semi-classical method in Godsi et al. [44], we compare the methods in the“flux detector” mode of the molecular hyperfine interferometer. This flux detector mode has a field profile asshown in Figure 3.3, but with the second arm assumed to be of zero length and B2 = 0. The field B1 is varied.For the sake of comparison, we also set the state selector and detector fields in the transfer matrix methodto 100000 gauss so that the basis changes performed by the transfer matrix method out of and into theseregions match well the Clebsch-Gordon transformation from |mImJ〉 to |Fm〉 and its inverse, as used by thesemi-classical method. Note that the off-diagonal elements of the full discontinuity matrix K (3.30) are still70430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)RP430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)RP430 435 440 445 450Magnetic Field, |B2| (gauss)1.00.50.00.51.0Calculated Signal (a.u.) RP430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)RDA430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)RDA430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)RDA430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.)ROM430 435 440 445 450Magnetic Field, |B2| (gauss)1.00.50.00.51.0Calculated Signal (a.u.) ROM430 435 440 445 450Magnetic Field, |B2| (gauss)1.00.50.00.51.0Calculated Signal (a.u.)ROM430 435 440 445 450Magnetic Field, |B2| (gauss)1.00.50.00.51.0Calculated Signal (a.u.)RUM430 435 440 445 450Magnetic Field, |B2| (gauss)0.50.00.51.0Calculated Signal (a.u.) RUM430 435 440 445 450Magnetic Field, |B2| (gauss)1.00.50.00.5Calculated Signal (a.u.)RUMFigure 3.6: Calculated signals close to the spin echo condition as functions of the magnetic field ofthe second coil |B2|. The field profile is depicted in Figure 3.3; B1 = 440gauss; and the signalwas sampled at a rate of 300 points per 20 gauss. Each of the plots was created with identicalparameters, except for the scattering transfer matrices Σ. The scattering transfer matrices areidentical for all energies and are randomly chosen for each plot as follows. First row – RandomPhases (RP): Σ =⊕9i=1 eiθi is a diagonal unitary matrix whose nine phases θi are randomlychosen from a uniform distribution of width 2pi . Second row – Random Diagonal Amplitudes(RDA): Σ =⊕9i=1 Ai is a diagonal matrix whose diagonal elements are randomly chosen from auniform distribution on the interval [0,1). Third Row – Random Orthogonal Matrices (ROM): Σis an orthogonal matrix randomly drawn according to the Haar measure on O(9). Fourth Row –Random Unitary Matrices (RUM): Σ is a unitary matrix randomly drawn according to the Haarmeasure on U(9). Here, randomly drawing according to the Haar measure can be understood asanalogous to drawing from the “uniform distribution” over the space of possible matrices [164].710 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RP0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RP0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RP0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RDA0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RDA0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RDA0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) ROM0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) ROM0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) ROM0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RUM0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RUM0 2 4 6 8 10 (kHz/gauss)0.00.20.40.60.81.0Fourier Amplitude (a.u.) RUMFigure 3.7: Fourier transforms of the signals shown in Figure 3.6 as functions of the generalized gyro-magnetic ratio γ . The panel labels are described in the figure caption to Figure 3.6 .only ∼10−5 at 100000 gauss, such that their neglect still does not invalidate our fully quantum formalismat these large field strengths. We also retain the rotation from the first branch to the second branch and setthe scattering matrix to 19. To maximize sensitivity of the comparison, we use only a single velocity whencalculating Pdetection in both methods. All other parameters, including the state selector and detector relativestate probabilities ηmImJ and κmImJ , are as per Appendix B.3. This allows for a test that includes all incomingand outgoing states and their relative phases at experimentally relevant conditions.72Before discussing the results of the comparison between the two methods, however, I will comment ontheir relative computational cost. Without performing a detailed analysis, the computational complexities ofthe two methods appear to be similar as the computational costs of both methods should scale proportionallyto the number of field regions included in the magnetic field profile used. As such, it appears that the detailsof how each method is programmatically implemented will determine whether one method out-performs theother computationally.I will now return to the results of the comparison between the two methods. We calculate the signalsPdetection (B1) from B1 to B1 + 10gauss for various values of B1. We include 1500 data points in these ten-gauss intervals. The calculated signals are compared between the two methods by calculating their relativeabsolute difference at identical conditions. This produces a relative absolute difference at each of the 1500magnetic field points. We then calculate the maximum, mean, and median relative absolute difference overthe ten-gauss interval. Figure 3.8 shows how these maximum, mean, and median values vary as a functionof B1. At low fields, there is no significant dependence on the magnetic field and the relative absolutedifference is below the expected experimental error. This lack of dependence on B1 is possibly due to someinherent numerical error present in the implementation of one or both methods, which masks any underlyingfield-dependence. At approximately 460 gauss, however, the error begins to increase with the magnetic fielduntil it saturates at approximately 46000 gauss at a relative absolute difference of approximately one. Thisincrease in error as a function of magnetic field points to a systematic difference between the two methods.To illustrate the difference between the two approaches in more detail, we plot the Fourier transforms ofthe calculated signals at various magnetic field strengths in Figure 3.9. At field strengths below 1000 gauss,little difference is observed. At higher field strengths, the feature locations agree, while the Fourier ampli-tudes differ. The feature locations are determined by the eigenvalues of the Hamiltonian, identical in bothmethods, while the amplitudes are a function of the relative phases and amplitudes of the wavefunctioncomponents. These amplitudes and phases are expected to differ between the two methods at sufficientlyhigh fields because of the approximations made in the semi-classical method.In particular, the semi-classical method accounts for most of the relative phase and amplitude changesinduced by the controlling magnetic fields. It does this by time-evolving the internal state vector for timesthat correspond to the time ti spent in each magnetic field by a molecule moving at its unchanged initial73100 101 102 103 104 105B1 (gauss)10 510 410 310 210 1100Relative Absolute DifferenceMaximumMedianMeanFigure 3.8: Maximum, median, and mean relative absolute difference between the calculated signalsobtained by the semi-classical method discussed in the Supplementary Material of Godsi etal. [44] and the present method, for various values of the controlling magnetic field. See Sec-tion 3.8 for a description of the semi-classical method and the field profiles used. The relativeabsolute difference between the calculated signals is calculated point-by-point as a function ofthe magnetic field. The mean, median, and maximum values are then calculated over the mag-netic field interval spanned by the calculated signal. The calculated signal is sampled at a rate of1500 points per 10 gauss; the magnetic field varies from B1 to B1 +10gauss for each calculatedsignal; a single velocity was included in the calculations; the scattering transfer matrix Σ = 19and is constant for all energies. All other parameters are listed in Appendix B.3.velocity. However, the semi-classical method ignores the small changes in the molecular velocity causedby the magnetic fields. These changes to the velocity modify the time spent in each magnetic field for eachindividual component of the internal state vector. Thus, ti should depend on the internal state |R〉. It is notimmediately clear how to include these state-dependent velocity changes into the semi-classical method,however.At low fields, these velocity changes and the dependence of ti on |R〉 are negligible and the fully quantumcalculations agree with the semi-classical results to at least 0.1% for fields below 1000 gauss. However, thisagreement can only be expected to occur for surfaces that do not change between the surface-impact eventsof the spatially-separated wavepacket components (discussed in Section 3.3). The maximum temporal sep-aration between these impact events, caused by the velocity changes, varies from a few to several hundredsof picoseconds. Many surfaces do change on this timescale, as has been measured in several helium-3 spinecho experiments [33–39]. In other words, the semi-classical method cannot be used to probe the dynam-ics of surfaces, while the method presented in this chapter opens the possibility to account for the surface74dynamics with molecular scattering experiments.0.0 2.5 5.0 7.5 10.0 (kHz/gauss)10 210 1100Fourier Amplitude (a.u.)1 gaussPresent MethodSemi-Classical Method0.0 2.5 5.0 7.5 10.0 (kHz/gauss)10 210 1100Fourier Amplitude (a.u.)1000 gaussPresent MethodSemi-Classical Method0.0 2.5 5.0 7.5 10.0 (kHz/gauss)10 210 1100Fourier Amplitude (a.u.)10000 gaussPresent MethodSemi-Classical Method0.0 2.5 5.0 7.5 10.0 (kHz/gauss)10 210 1100Fourier Amplitude (a.u.)46000 gaussPresent MethodSemi-Classical MethodFigure 3.9: Fourier amplitudes of the calculated signals at various magnetic field values computedwith the semi-classical method discussed in the Supplementary Material of Godsi et al. [44](orange) and the present method (blue) as functions of the generalized gyromagnetic ratio γ . Thecalculation conditions are identical to those of Figure 3.8.3.9 Conclusion and Future DirectionsIn this chapter, we have developed a theoretical framework for simulating a surface-sensitive molecularhyperfine interferometer. The approach treats the interferometer as an effective one-dimensional system,accounting for the real experimental geometry by rotating the quantization axis of the hyperfine states at thescattering point. The time evolution of the molecular states is described fully coherently and accounts for themixing of the hyperfine states and momentum changes induced by the magnetic fields in the experiment. Thepresent approach is fully quantum mechanical and includes a full description of the internal-state-dependent75spatial superpositions imposed on the molecular wavepackets by the controlling magnetic fields. This opensthe door for a description of molecular scattering experiments that aim to probe surface dynamics on thepicosecond to hundreds of picosecond time scale. To build the framework, we have derived and implementeda transfer matrix formalism that accounts for the internal (hyperfine) degrees of freedom of molecules andthat allows for efficient computation of the experimental signal.In the present work, the molecule-surface interaction is accounted for by a scattering transfer ma-trix (a transformed version of the scattering matrix) that is suitable for the description of experimentswhere the surface changes either much more slowly or much more quickly than the molecule-surface orwavepacket-surface interaction times. The extension to arbitrary surface dynamics (currently under investi-gation) requires a time-dependent scattering transfer matrix that reflects the underlying time-dependence ofthe molecule-surface interaction potential. We have demonstrated, using the specific case of oH2, how thedifferent features of the time-independent scattering transfer matrix, such as the phases of the diagonal ele-ments, impact the experimental signal. In addition, we have shown that the experimental signal is sensitiveto inelastic scattering events that change the projection quantum numbers of the molecular hyperfine states.The present approach also sets the stage for solving the inverse scattering problem in molecular hyper-fine interferometry by means of machine learning approaches, such as Bayesian optimization [45, 64]. Forexample, the difference between the experimental observations and the results of the transfer-matrix compu-tations can be efficiently minimized by using the Bayesian optimization algorithm to intelligently vary thescattering matrix elements. We discuss this approach in Chapter 4.The formalism presented here is general to all closed shell molecules and is flexible to describe variousexperimental setups. It can be used to explore various experimental protocols and evaluate their effective-ness at determining various molecule-surface interactions and surface properties. Thus, this chapter providesthe theoretical framework necessary to interpret a wide range of molecular hyperfine interferometry experi-ments. These new experiments are poised to provide new information about molecule-surface interactions,surface morphologies, and surface dynamics, such as the impact of molecular orientation on surface reactionrates and the details of molecule-surface potential energy surfaces.763.10 SummaryThis section contains a summary of the key results presented in Chapter 3:• Proposed using generic closed shell molecules for molecular hyperfine interferometry experiments(Section 3.2)• Developed a fully quantum formalism to describe molecular hyperfine interferometry with closedshell molecules (Sections 3.3.1-3.6)• Derived a transfer matrix formalism with internal degrees of freedom (Section 3.6)• Applied the fully quantum formalism to the description of MHI experiments using oH2 molecules(Section 3.7)• Showed that various features of the scattering matrix affect the experimental signal (Section 3.7)• Compared the quantum formalism with a semi-classical formalism and determined that have the sameinformation at low fields, for static surfaces (Section 3.8)• Showed how the quantum formalism opens the door for a formalism that can account for dynamicsurfaces (Section 3.8)• Set the stage for the determination of the molecule-surface scattering matrix via machine learningtechniques (Sections 3.3.1-3.6)77Chapter 4Bayesian optimization to determine themolecule-surface scattering matrix4.1 IntroductionThis chapter discusses our first applications of Bayesian optimization to the inverse scattering problem formolecular hyperfine interferometry (MHI) experiments. I begin this chapter by discussing how we approachthe inverse scattering problem from an optimization perspective in Section 4.2. Then, I give a non-technicaloverview of Bayesian optimization by first describing Gaussian processes (Section 4.3.1), then Gaussianprocess models (Section 4.3.2), and finally the Bayesian optimization algorithm itself (Section 4.3.3). Idiscuss our initial applications of Bayesian optimization to determine the scattering matrix elements from agiven signal in Section 4.4. I then conclude the chapter in Section 4.5 and give a bulleted summary of thechapter results in Section 4.6.4.2 The Inverse Scattering ProblemIn Sections 3.2–3.6, we developed a framework that describes the MHI experiments and that allows us toconnect the scattering transfer matrix elements to the detector signal. In general, however, we do not knowthe scattering matrix elements. We would like to determine these matrix elements from experimental datain order to either gain information about the surface and/or improve surface scattering calculations. This isthe inverse scattering problem for MHI.78One simple way to determine the scattering matrix elements is to vary the values of the matrix elementsuntil the calculated signal matches the experimental signal. (Note that, in this chapter, I use the terms“scattering matrix” and “scattering transfer matrix” interchangeably, for if we know the matrix elements ofthe scattering matrix, we can calculate the elements of the scattering transfer matrix, and vice versa.) If therewere only a small number of parameters that define the scattering matrix and if the calculation of the signalPdetection were fast, such a naive approach would work rather well. However, a simple grid search of a three-parameter space requires 1000 evaluations for even a low-resolution 10×10×10 grid. If each calculationtakes several minutes, this begins to be computationally expensive; if we need to explore a parameter spaceof possibly more than 81 dimensions, such an approach is completely prohibitive.As such, we need a more sophisticated approach. One such approach is Bayesian optimization, whichis a machine learning technique particularly well-suited to situations where every evaluation of the functionto be optimized is computationally expensive [45]. This approach has been used for the inverse scatteringproblem in quantum reaction dynamics [64]. To determine the scattering matrix parameters in MHI experi-ments via optimization, we minimize the mean square error (MSE) between the calculated and experimentalsignals. Ideally, we want to find the global minimum of this MSE and all local minima that also give a smallMSE. We want these local minima as there may be multiple scattering matrices that can explain the exper-imental signal and because the global minimum may be over-fit to noise in the experimental signal. Thisnoise in the experimental signal could arise from imperfections in the magnetic fields, errors in the mea-surement of the magnetic field profile, errors or variation in the current supplied to the solenoids, variationof the molecular velocity distribution between experimental runs, and variation of the detector sensitivity.Ultimately, we may not be able to determine a unique scattering matrix, but we should, at least, be able todetermine classes of matrices that match the data and/or put bounds on the values of the matrix elements.I will now give a non-technical overview of Bayesian optimization and then present the results of ourinitial applications of Bayesian optimization to the inverse scattering problem for MHI experiments.4.3 Bayesian OptimizationAs previously mentioned, Bayesian optimization is most advantageous when each evaluation of the functionto be optimized is rather expensive. Bayesian optimization relies on using a Gaussian process model to guide79the algorithm and, in turn, a Gaussian process model utilizes a Gaussian process to give an interpolatingmodel of the data obtained from function evaluations. Additionally, the Gaussian process model providesan estimate of the uncertainty in the model, which is vital for the algorithm to function correctly. In thissection, I will give a high-level introduction to Gaussian processes, Gaussian process models, and Bayesianoptimization.4.3.1 Gaussian ProcessesEssentially, a stochastic process (or simply process) can be thought of as a continuum limit of a joint proba-bility distribution, such that there is a probability distribution for every point on the real line (there are moregeneral and more rigorous definitions [165], but this suffices for the present purposes). A Gaussian processis then a process where this distribution is Gaussian. To build up an intuition for what a Gaussian processis, let us first look at a one dimensional Gaussian distribution P(x) = 1√2piσ2e(x−µ)22σ2 [166]. This distributionis characterized by both a mean µ and a standard deviation σ .If we look at a Gaussian distribution over two variables, we have the joint probability distributionP(x1,x2) = (2pi)−1(σ21σ22 −σ21σ22ρ2)− 12e11−ρ2((x1−µ1)22σ21+(x2−µ2)22σ22− ρ(x1−µ1)(x2−µ2)σ1σ2)[166]. Here, the distributionhas five parameters: two means µ1,µ2; two standard deviations σ1,σ2; and the covariance σ212 ≡ ρσ1σ2,where ρ is the correlation. We can also define the symmetric covariance matrix Σcov =σ21 σ212σ212 σ22 andthen write P(x1,x2) = P(~x) = (2pi)−1[det(Σcov)]− 12 e− 12(~x−~µ)TΣ−1cov(~x−~µ), where~x≡ (x1,x2) and~µ ≡ (µ1,µ2)[166].We can then expand this to the d-dimensional joint distribution P(~x)= (2pi)−d2[det(Σcov)]− 12 e− 12(~x−~µ)TΣ−1cov(~x−~µ)where~x≡ (x1,x2, ...,xd) and ~µ ≡ (µ1,µ2, ...µd) are both d-dimensional and there are d standard deviationsσi, and(d2)covariance terms σ2i j. That is Σcov is the symmetric d by d matrix:Σcov =σ21 σ212 · · · σ21dσ212 σ22 · · · σ22d....... . ....σ21d σ22d · · · σ2dd.8020 10 0 10 20z3210123x(z)0.0×1008.0×10 41.6×10 32.4×10 33.2×10 34.0×10 34.8×10 35.6×10 36.4×10 3P(x)3 0 33 0 3Figure 4.1: Illustration of a Gaussian process. The colour indicates the probability density P(x) that agiven function x˜(z) drawn from the Gaussian process would pass through the point (z,x(z)). Eachvertical slice of this distribution has the form of a Gaussian, as shown in the two insets. Each ofthese insets corresponds to the slice indicated by the neighbouring dashed line. The Gaussian dis-tribution within every vertical slice is parametrized by the mean and standard deviation functionsµ(z) and σ(z). The blue line shows µ(z) for this Gaussian process.Now, we could already choose the number of dimensions to be infinite, but this would still be a countablyinfinite number of dimensions. Instead, let us now promote ~x from having components with a discreteindex to those with a continuous index. That is, the number of dimensions of the joint distribution is nowuncountably infinite.This is rather like the difference in the matrix representations of operators that act on bound states andon continuum states in quantum mechanics. The hydrogen atom potential has an infinite number of discretebound states, but it also has a continuum of scattering states [167]. This continuum of scattering states hasa Hamiltonian whose “matrix” representation has an uncountably infinite number of dimensions, while thebound state Hamiltonian has a matrix with a countable, yet infinite, number of dimensions. Ultimately, thecomponents xi of~x change from being labelled by i ∈ Z to z ∈ R, such that we have x(z).81Thus, we now have a situation where there is a Gaussian distribution over x(z) at every point z on the realline, with correlations in the joint probability distributions between points z and z′. This can now be calleda Gaussian process. Figure 4.1 gives an illustration of a Gaussian process with the insets showing the Gaus-sian distributions along a vertical slice of the process, while Figure 4.2 gives an alternative representation ofthe same Gaussian process. Alternatively, we can consider x(z) as no longer describing a randomly chosenvector, but a randomly chosen function. We may then be tempted to write the following probability distribu-tion over this function: P(x) ∝ e12∫dzdz′[(x(z)−µ(z))Σ−1cov(z,z′)(x(z′)−µ(z′))], which would be characterized by theparametrizing functions µ (z), and Σcov(z,z′). Unfortunately, while potentially intuitive, this is ill-defined,even if functional forms of this sort are useful in quantum field theory [92].At the moment, however, the exact definition is less important than the notion that the Gaussian processrepresents an ensemble of functions, whose ensemble mean is µ(z) and whose correlations are governedby Σcov(z,z′). That is, a process can be thought of as a probability distribution over a class of functions;the random variable is a function, rather than a vector or scalar. More precisely and more generally, wecan define the Gaussian process to be “a collection of random variables, any finite number of which havea joint Gaussian distribution” [165]. In the situation we are describing, this collection of random variablesis uncountably infinite and indexed by a real number. While this formal definition may appear unintuitive,it is made clearer in light of the discussion above. So far, we have described a stochastic process thatis a distribution over a singly-valued function; however, processes for multi-valued functions can also bedefined.It happens that the correlation functions Σcov(z,z′)correspond to the so-called “kernels” that appear inthe context of the Gaussian process models discussed in the next section [45]; as such, I will often refer to thecorrelation functions as kernels for the remainder of this thesis. Furthermore, there are many such kernelsthat are used to define Gaussian processes, each of which bestows particular properties on the Gaussianprocess. For example, the squared exponential kernel ensures that the stochastic process itself is infinitelymean square differentiable [165]. This mean square derivative can be thought of as the generalization ofdifferentiation to stochastic processes [165]. Figure 4.3 illustrates functions drawn from Gaussian processesdefined with various kernel functions.8220 10 0 10 20z3210123x(z)Figure 4.2: An alternative representation of the Gaussian process illustrated in Figure 4.1. The blueline is the same as that in Figure 4.1 and corresponds to the mean µ(z). The shaded blue arearepresents the spread of the Gaussian process and encompasses µ(z)±2σ(z), where σ(z) is thestandard deviation at the point z.4.3.2 Gaussian Process ModelOne application of Gaussian processes is to interpolation. While the true underlying approach to this typeof interpolation is that of Bayesian inference, I will give an intuitive description of the results and refer toRasmussen and Williams [165] for a more detailed explanation.Let us assume that we have a set of n data points D = {xi,yi} between which we want to interpolate.As such, we are assuming that there is some functional relationship between x and y. We will hold fewassumptions about this function, but let us assume that it is smooth.Now, let us take an initial Gaussian process, with mean µ(x) = 0 and the square exponential kernelΣcov(x,x′)= e−|x−x′ |22l2 . This kernel describes correlations in the process that decay with a Gaussian profile(and length scale l) [165]. See the upper left panel of Figure 4.3 for examples of functions that can come froma Gaussian process defined with this kernel. Presumably, there exists a drawn function that approximatesthe underlying relationship in the data. However, our goal is not to determine this function. Instead, wewill use the mean of the Gaussian process as the prediction and fit the process as a whole to the data. Thisfitting will narrow the distribution over the functions near the data points, but leave it broad in regions withlittle data. In this way, we will obtain both a reasonable estimate of the underlying functional relationship8320 10 0 10 20x3210123ySquare Exponential Kernel20 10 0 10 20x3210123yMatern Kernel, = 1.520 10 0 10 20x3210123yExponential Sine Squared Kernel20 10 0 10 20x3210123yMatern Kernel, = 0.5Figure 4.3: Demonstration of the different functions that can arise when drawn from Gaussian pro-cesses that are defined with different kernels. Note that while each Gaussian process has µ(z) = 0and σ(z) = 1, the behaviour of each drawn function varies drastically between the different ker-nels. The mathematical form of the square exponential kernel, the Mate´rn kernels, and the expo-nential sine squared kernel can be found in Ref. [165]. Note that the exponential sine squaredkernel is also known as the periodic kernel [168].between x and y and information about the uncertainty of our estimate.Before performing this fitting procedure, we initially have some data D and a presumed initial distribu-tion over the functions (our prior). This data and the initial distribution (i.e. the Gaussian process) are shownin the upper left panel of Figure 4.4. As we are generating the data from a known functional relationshipbetween x and y, I also plot this relationship for reference. Let us first perform a fit with only two data points.Now, assuming no noise in the system, we know that the true function must pass through the data points.As such, our fitting procedure forces µ(x3) = y3 and µ(x4) = y4 and removes any uncertainty in the modelat x3 and x4. This results in the upper right panel of Figure 4.4. I refer to Rasmussen and Williams for the84mathematical details behind this fitting procedure [165], as they are not needed for the present discussion.One can see that the model is beginning to approximate the underlying function, albeit extremely poorly.Furthermore, there is little uncertainty in the process near (x3,y3) and (x4,y4), while the uncertainty growsas one moves away from each of these two points.If we perform the fitting procedure with more data points, we obtain a much better fit; see the bottompanel of Figure 4.4. While there is still uncertainty between the data points, it is small. Interestingly, themodel does include some extrapolation past the region spanned by the data points. However, the uncertaintygrows rather quickly outside of this region; the influence of the last data point is decaying as a functionof distance. Thus, the model is telling us where we can be more confident in its prediction (there is smalluncertainty in the model) and where we should be less confident (there is large uncertainty in the model).We have now made a Gaussian process model for our collection of data. If we were to add more datapoints, the model would become increasingly more accurate. As the kernel parameter l is varied duringthe fitting procedure to improve the fit, we also receive an automatic measure of a relevant length scale inthe data. Furthermore, we are not restricted to modelling functions of a single variable, but can performthis analysis with high dimensional data rather efficiently, even to around 100 dimensions [169, 170]. Thenumber of data points needed for a good fit actually scales linearly in the number of dimensions [171],making the method even more efficient for larger dimensions (for up to around 100 dimensions [169, 170]).4.3.3 Bayesian OptimizationIn this section, I now wish to turn to the method of Bayesian optimization. We start with some function fwith d input dimensions that we wish to optimize. One (simplified) algorithm for Bayesian optimization isas follows [45]:1. Evaluate f at n points~xi, randomly chosen (see panel (i) of Figure 4.5)2. Fit a Gaussian process model to the set of data points {~xi, f (~xi)} (see panel (ii) of Figure 4.5)3. Find the parameters ~xnew that minimize the Gaussian process model, taking the uncertainty into ac-count (i.e., ~xnew minimizes the so-called “acquisition function”; see the red cross in panel (ii) ofFigure 4.5)4. Evaluate f (~xnew) and add the result to the set of data points8520 10 0 10 20x3210123y20 10 0 10 20x3210123y20 10 0 10 20x3210123yFigure 4.4: An example of fitting a Gaussian process model. The upper left panel shows the truefunction (dashed line); the sample points (black circles), whose x values are drawn via Latinhypercube sampling; and the unfit Gaussian process model, whose mean function µ(x) = 0 (blueline) and whose standard deviation function σ(x) = 1. The shaded blue area encompasses µ(x)±2σ(x). The upper right panel shows the Gaussian process model fit to two sampling points. Thebottom panel shows the Gaussian process model fit to all five sampling points. Note that µ(x) istaken as the prediction of the Gaussian process model, while σ(x) is taken as the uncertainty inthe model at x.5. Repeat from step 2 until~xnew converges to a particular point (see panels (iii-xii) of Figure 4.5)Figure 4.5 illustrates this process. There are a few key things to note about the above algorithm. First,instead of evaluating the computationally expensive function f many thousands of times, as is often requiredfor optimization, the algorithm is evaluating the computationally cheaper Gaussian process model thousandsof times. (This occurs when minimizing the acquisition function in step 3.) The trade off is that we mustre-perform this minimization every time we update the Gaussian process model.8620 10 0 10 20x3210123y(i)20 10 0 10 20x3210123y(ii)20 10 0 10 20x3210123y(iii)20 10 0 10 20x3210123y(iv)20 10 0 10 20x3210123y(v)20 10 0 10 20x3210123y(vi)20 10 0 10 20x3210123y(vii)20 10 0 10 20x3210123y(viii)20 10 0 10 20x3210123y(ix)20 10 0 10 20x3210123y(x)20 10 0 10 20x3210123y(xi)20 10 0 10 20x3210123y(xii)Figure 4.5: Demonstration of the Bayesian optimization algorithm in one dimension. Panel (i) showsthe initial preparatory step of evaluating the function f at seven initial points xi, chosen via Latinhypercube sampling. The dashed black line shows the true function, while the black circles are(xi, f (xi)). The blue line shows the mean of the unfit Gaussian process model, while the blueshaded area shows the uncertainty in the model by encompassing µ ± 2σ . Panel (ii) shows theprediction and uncertainty of the Gaussian process model that is fit to the initial seven points.Also shown are the acquisition function a(x) = µ(x)− 2.5σ(x) (dashed orange line) and theminimum (xnew,a(xnew)) of this acquisition function in the interval [−20,20] (red cross). Panel(iii) shows the same features, but with the Gaussian process model now fit to eight sample points.This eighth point is (xnew, f (xnew)), where xnew minimizes the acquisition function of the previouspanel. Panels (iv-xii) show the continuation of the Bayesian optimization algorithm, where a newsample point (xnew, f (xnew)) is added every panel.87Secondly, taking the uncertainty into account when minimizing the model is a key step. The Gaussianprocess model gives us not only an interpolation, but also information about the uncertainty of that interpo-lation. We can incorporate this uncertainty by finding the minimum not of the mean of the Gaussian processmodel, but of the acquisition function a(~x) ≡ µ(~x)−κσ(~x), where µ is the mean of the Gaussian processmodel, σ its standard deviation, and κ a tuning parameter. By including the model uncertainty in the min-imization, we are improving the likelihood that we reach the global minimum; the algorithm is forced toexplore regions of the parameter space where there is insufficient information to conclude the absence of theglobal minimum.Finally, the flip side of incorporating this uncertainty is that we can ignore regions where there is suf-ficient information to conclude the absence of the global minimum. Thus, this algorithm evaluates thefunction at points more likely to provide useful information than would algorithms based on simple randomor grid sampling. This reduces the number of required function evaluations and improves the efficiency ofthe algorithm.To conclude this section, I will state that the algorithm I mentioned above describes only the core ideasbehind Bayesian optimization. While following the above algorithm will work, there are various modifi-cations one can make to improve the efficiency of the algorithm. For example, there are many possibleforms that the acquisition function can take; what I mention above is known as the lower confidence bound(LCB) acquisition function. Other forms include expected improvement and probability of improvementacquisition functions [172]. One can also improve the fitting of the Gaussian process model to the pointwhere it can be proven that the Bayesian optimization algorithm converges to the global minimum, thoughthis makes the algorithm more complex [173].4.4 Preliminary ResultsIn this section, I will describe our preliminary application of Bayesian optimization to the problem of deter-mining the scattering matrix from the MHI experimental data. We are particularly considering experimentsthat utilize ortho-hydrogen as the probe molecule, and thus working with a 9×9 scattering transfer matrix.In the ideal case, we would like to obtain all 81 complex matrix elements from the experimental data. Thisis an ambitious goal, however, and we begin by simply considering diagonal scattering matrices, which88describe elastic scattering. These matrices have up to nine possible parameters, the nine variables θi thatdefine the phases of the nine diagonal elements eiθi . The RP matrices used for the first rows of Figures 3.6and 3.7 are of this form, for example.However, even determining these nine parameters is not straight-forward. As such, we do not work withexperimental data immediately, but rather choose a set of target parameters θ ∗i . We take these target phasesand calculate the expected signal. Then, we attempt to determine the parameters from this signal (i.e., we“pretend” that this signal is experimental data). This approach allows us to tune the algorithm, to diagnosewhether the Gaussian process model fit is well-behaved and stable, and to know whether the algorithm isreliably converging to the correct values. Once we are confident that the algorithm is working under theseknown conditions, we can work with the real experimental data, where there will be additional issues to bedealt with, such as noise in the signal.Our computational setup is the same as that presented in Section 3.7 and has the field profile shown inFigure 3.3. We choose B1 = 400gauss. To determine the scattering matrices, we match the calculated signalto the “experimental” signal. That is, we minimize the MSE between the two curves. Note that we workwith the MSE instead of the root mean square error (RMSE) as the RMSE has a cusp that occurs at thetarget values. That is, near ~θ ∗, the RMSE∼ |~θ −~θ ∗|, which produces difficulties when fitting the Gaussianprocess model. Instead, near the target value ~θ ∗, the MSE∼ (~θ −~θ ∗)2, which is differentiable at least onceand improves the stability of the fitting procedure. One disadvantage is that the MSE is flatter near the globalminimum, which increases the number of algorithmic steps required to converge to a particular precision.Additionally, we calculate Pdetection at only ten values of B2 instead of the usual ∼300 values, in order toexpedite this refinement stage of the algorithm.For the Bayesian optimization algorithm itself, the kernel is of the form Σcorr(~θ , ~θ ′) = C(~θ − ~θ ′)×Mate´rnν(~θ − ~θ ′), with the parameter ν = 1.5. This includes the Mate´rn kernel multiplied by the constantkernel C(~θ − ~θ ′). The Mate´rn kernel with ν = 1.5 causes the stochastic process to be once mean squaredifferentiable [165], while the constant kernel allows for scaling of the function range. We choose the Mate´rnkernel with ν = 1.5 as the MSE is also differentiable once at the global minimum. We tried various otherforms of the kernel, but obtain the most reasonable performance of the algorithm with the above kernel, sofar; it is likely that the above choice of kernel is still not optimal. We select the N initial random points89~θ j via Latin hypercube sampling [174]. Latin hypercube sampling shows a greater uniform sampling ofthe parameter space, in comparison to drawing the points according to a uniform distribution. This greatersampling uniformity reduces the amount of redundant information arising from sample point clustering,which is likely when sampling from a uniform distribution.In order to proceed in a relatively systematic way, we first applied Bayesian optimization to determinea scattering matrix parametrized by a single phase θ1, with all other phases held constant. As this met withsuccess, we increased the number of parameters to two (θ1 and θ2); this also met with success. Given thesesuccesses, we now increase the number of parameters to three and present the results here.As the underlying core of the Bayesian optimization algorithm is the Gaussian process model, we vi-sualize this model to ensure that the obtained fits are reasonable. As the mean of the Gaussian process isa function of three variables, we plot the mean on slices of the full parameter space. We choose each sliceto pass through the known target phases (θ ∗1 ,θ ∗2 ,θ ∗3 ). In the case where the target phases are unknown, wewould choose these slices to pass through the current best estimate of these parameters. In Figure 4.6, weshow the mean of a Gaussian process model for an optimization run where the (randomly chosen) targetphases are (θ ∗1 ,θ ∗2 ,θ ∗3 )≈ (3.468,2.482,4.488); the remaining six phases are set to zero. To obtain this plot,we choose 300 initial sample points and perform 55 steps of the Bayesian optimization algorithm, whichresults in a total of 355 evaluations of the MSE. We set the LCB acquisition function tuning parameterκ = 0.01, in line with the values used in Ref. [64]. One can clearly see that the mean of the Gaussian pro-cess model is rather smooth and that the minimum of the model is near the known target phases (θ ∗1 ,θ ∗2 ,θ ∗3 ),shown by the orange circles.In addition to plotting the mean of the Gaussian process model, we can monitor the progress of thealgorithm by looking at the RMSE and the determined phases (θ1,θ2,θ3) as a function of the algorithm step.We show these observables in Figure 4.7. One can see that the algorithm converges rather quickly, with areasonable RMSE of 0.31% between the target and calculated signals occurring after just 20 optimizationsteps. The determined phases also converge rather quickly, with the phases determined to be approximately(3.472,2.494,4.484) after 30 optimization steps. The relative errors between these values and the targetvalues are (0.11%,0.46%,−0.07%), presumably well within expected experimental errors. These resultsshow that we are able to at least partially determine the scattering matrix from a given signal. Furthermore,90Figure 4.6: Gaussian process model of the mean square error between a calculated and a target Pdetectioncurve as a function of the three phases θ1,θ2,θ3. These angles are the phases of a 9×9 diagonalscattering matrix whose diagonal elements Σ j j are of the form eiθ j , with θ4 to θ9 set to zero. Thegoal is to find the phases θ1,θ2,θ3 that minimize the mean square error and thereby determinethe target scattering matrix. Each subplot of this panel is a 2D slice of the 3D surface, with eachslice passing through the known target phases (θ ∗1 ,θ ∗2 ,θ ∗3 ). The orange circles indicate the knowntarget phases. The Gaussian process model is fit with 300 randomly selected points and 55 pointsselected by the Bayesian optimization algorithm.91Figure 4.7: Left panel: The phases as determined by the Bayesian optimization algorithm as a functionof the number of algorithm steps. The dashed lines indicate the known target phases θ ∗1 ,θ ∗2 ,θ ∗3 .Right panel: The root mean square error between ten points of the target Pdetection curve and thePdetection curve calculated with the determined phases (as indicated in the left panel), for each stepof the Bayesian optimization algorithm.we are able to determine the target phases to an error of less than 0.5% with just 330 function evaluations.This is significantly more efficient than a simple 10× 10× 10 grid search, which could at most determineeach phase to a resolution of 2pi/10≈ 0.63.Interestingly, the algorithm as presented does not always perform this well. Under identical conditions,except for another set of target phases (3.190,5.808,3.002), it takes 760 steps of the algorithm (a total of1060 function evaluations) to reach a RMSE of 1.7% between the calculated and target signals, with therelative errors in the determined phases (0.65%,1.05%,0.10%). We are still determining the reason for thisrelatively poor performance (though it is still significantly better than a grid search) and still determiningwhich of these two examples are reflective of the typical behaviour of the algorithm as presented. Certainly,more tuning of the algorithm is required. This may include modifying κ , selecting a different kernel, orchanging the Gaussian process model fitting procedure. Indeed, one can modify the Gaussian process modelfitting procedure to ensure that the algorithm always converges to the global minimum [173], which may bean avenue worth exploring.We have also tried applying the algorithm to problems with larger numbers of unknown phases, but haveobtained no successful results so far. We will first understand and improve the behaviour of the algorithmfor problems with three unknown phases before again applying the algorithm to larger numbers of phases.92Once we have tuned the algorithm to perform well at determining the known target phases, we will attemptto fit the diagonal scattering matrix to true experimental data.Interestingly, we may not need to determine nine phases from the experimental signal, as symmetriesin the molecule-surface interaction may allow us to work with only three phases. That is, the primaryinteraction of an oH2 molecule with the surface is steric. This would modify the rotational degrees offreedom while leaving the nuclear spin degrees of freedom unchanged, such that Σ−1mJ = Σ0mJ = Σ1mJ ,where the fist index corresponds to mI . Finally, while we do not know if a diagonal matrix is the correctform of the scattering matrix, the results, whether positive or negative, will guide us towards the best nextsteps.4.5 Conclusion and Future DirectionsIn this chapter, we have discussed our initial applications of Bayesian optimization to the inverse scatteringproblem for MHI experiments. The goal of this inverse scattering problem is to determine the scatteringmatrix elements from the experimental signal. Our approach to this problem is to minimize the differencebetween the experimental signal and the results of the transfer-matrix computations by varying the elementsof the scattering matrix. We performed this minimization efficiently by using the Bayesian optimizationalgorithm to intelligently vary the scattering matrix elements. In this way, we have successfully appliedBayesian optimization to determine the scattering matrix corresponding to a target calculated signal, whenthat matrix is diagonal and determined by three or fewer parameters.Building on this foundation, we will tune the algorithm to improve performance and then apply it tosituations where the diagonal scattering matrix is defined by all nine phases. In addition to these diagonalscattering matrices, we will also treat different classes of scattering matrices, such as the RDA, ROM andRUM matrices discussed in Section 3.7.2, as well as fully random matrices. Also, all of the work discussedin this chapter has assumed that the scattering matrices are constant as a function of incident energy. This isknown to not be the case and we may expand each scattering matrix element Σi j(E) in a Taylor series; or, aLaurent series, if we expect the presence of some resonances.Furthermore, there is certain to be some experimental error or noise in the measured signals and it isimportant to take this into account. As such, we will determine how errors in the experimental signal affect93our determination of the scattering matrix parameters. Indeed, given the complexities involved, we mayonly be able to determine the parameters to low precision or only determine the classes of matrices that fitthe experimental data. Even if we can only determine classes of matrices, however, this would still give usinformation about what symmetries in the molecule-surface interaction are preserved or broken.Additionally, we may explore other types of optimization algorithms. If we are able to analyticallydetermine the gradient of the MSE from the transfer matrix formalism discussed in Chapter 3, we could takeadvantage of gradient-based methods and achieve improved performance [175, 176].Finally, we may attempt a full Bayesian approach to this problem and begin with a prior distributionover the scattering matrix parameters, which can be updated via Bayesian inference, given the experimentaldata [177]. This would allow us to place error bars on the scattering matrix parameters.The work presented in this chapter serves as the first step on the path to determine scattering matricesdefined by more than three parameters and to use Bayesian optimization to determine the properties ofthe scattering matrix elements compatible with a given experimental measurement. These scattering matrixproperties can then be used to gain physical insight into molecule-surface interactions and surface propertiesand used to test approximations used in ab initio calculations.944.6 SummaryThis section contains a summary of the key results presented in Chapter 4:• Applied Bayesian optimization to determine a scattering matrix parametrized by three phases, fromdata generated from known target parameters (Section 4.2)• Demonstrated that Bayesian optimization applied to the inverse scattering problem for MHI experi-ments is significantly more efficient than a grid search with comparable resolution, for a three dimen-sional parameter space (Section 4.2)• Demonstrated that tuning of the Bayesian optimization algorithm is required to improve performancefor the inverse scattering problem (Section 4.2)• Laid the foundation to efficiently determine, via Bayesian optimization, matrices defined by morethan three parameters95Chapter 5Monte Carlo enhanced with Bayesian modelcalibration5.1 IntroductionIn this short chapter, I propose the application of Bayesian model calibration to reduce the negative impactof critical slowing down and/or the sign problem in classical and quantum Monte Carlo calculations. InSection 5.2, I discuss Monte Carlo calculations and some of their limitations. In Section 5.3, I give a non-technical description of Bayesian model calibration. Then, in Section 5.4, I discuss a simplified version ofBayesian model calibration that is less powerful, but easier to implement. In Section 5.5, I provide a proof ofprinciple demonstration that Bayesian model calibration can correct the values of poorly converged MonteCarlo calculations with the information from only a few well-converged calculations. I conclude and discussfuture applications of Bayesian model calibration in Section 5.6 and give a chapter summary in Section 5.7.5.2 Monte Carlo CalculationsSince the landmark paper by Metropolis et al. in 1953 [178], Markov chain Monte Carlo calculations (hence-forth referred to as Monte Carlo calculations) have been widely used throughout science and engineering.In terms of complex quantum systems, these Monte Carlo calculations have been used to determine theproperties of superfluids [71, 179, 180], calculate the electronic structure of atoms, molecules, and solids[181–183], study the Hubbard model [73, 184], and investigate Fermi-polarons [185, 186], for example.96Despite their massive success, however, Monte Carlo calculations do have some, rather famous, limita-tions. For example, when calculating observables near a second-order phase transition, basic Monte Carlomethods exhibit what is known as critical slowing down [187]. This occurs as sequential Monte Carlosamples become too strongly correlated near the phase transition. More specifically, for Monte Carlo calcu-lations the variance of an observable often scales as (N)−12 , where N is the number of Monte Carlo samples[188]. However, N properly counts the number of independent Monte Carlo samples. Near a second-orderphase transition, the correlation length of the system increases, which also causes an increase in the numberof sequential Monte Carlo samples that are correlated (that is, the autocorrelation time increases) [189, 190].Hence, many more Monte Carlo samples must be generated to obtain reasonable accuracy; this can quicklybecome prohibitive.Another major limitation of Monte Carlo methods is known as the sign problem. For example, indiagrammatic Monte Carlo, the algorithm calculates observables by generating and summing a series ofFeynman diagrams [191]. Unfortunately, for fermionic systems these diagrams often have an oscillatingsign and nearly cancel each other out in the sum [74]. As such, all, or almost all, diagrams of the same orderneed to be summed to obtain an accurate result. (Roughly, the diagram order is the number of vertices inthe diagram.) However, the number of diagrams in a given order n scales as n! and calculations that do notconverge at low diagram orders are difficult to complete [74, 192].Many different approaches have been made to tackle these limitations of Monte Carlo methods. Ap-proaches to the issue of critical slowing down often involve updating the system state more efficiently, suchas with cluster updates [189, 190] or loop algorithms [184]. An approach to the sign problem in diagram-matic Monte Carlo involves using the Dyson equation to reformulate the summation of the bare Feynmandiagrams into a self-consistent formulation that requires explicit summation over a smaller set of diagrams;this is bold diagrammatic Monte Carlo [74, 185, 186].In summary, as one nears a second order phase transition or works with a fermionic system, it canbecome difficult to reduce the uncertainties in the calculated observables because of the large number ofMonte Carlo samples or diagram orders required; to then map out these observables as a function of theinput parameters is prohibitive. However, these difficulties are only true for well-converged calculations.Calculations with smaller numbers of samples or including only low order diagrams can still be performed97easily, albeit with an inaccurate final answer. At this point, we ask: Instead of trying to converge the MonteCarlo calculations by increasing the number of samples, can we just correct the poorly converged calcula-tions in an efficient manner? Here, we propose to do this correction with machine learning techniques.In recent years, machine learning techniques have been applied to improve the efficiency of various nu-merical methods. For example, restricted Boltzmann distributions have been used to learn effective Hamil-tonians that are then used to propose more efficient updates in classical Monte Carlo calculations [77, 78].Also, Cui and Krems have used Bayesian model calibration to improve the efficiency of scattering calcu-lations by combining quantum scattering calculations with classical scattering calculations [60]. That is,classical calculations are used to calculate scattering cross sections at many points of the parameter space asthey are computationally low-cost. Then, with just a few computationally expensive quantum calculations,Cui and Krems correct the classical calculations and obtain a reasonably accurate prediction of the scatteringcross section as a function of the input parameters.Here, we propose to use Bayesian model calibration [193] to partially alleviate both the sign problemand critical slowing down in Monte Carlo calculations by using Bayesian model calibration to correct easily-obtained, but poorly-converged, Monte Carlo calculations with the information gleaned from a small numberof difficult, but well-converged, calculations. I give a non-technical overview of Bayesian model calibrationin the next section.5.3 Bayesian Model CalibrationBayesian model calibration (BMC) is a machine learning technique based on the Gaussian process modeldiscussed in Section 4.3.2. BMC was proposed in 2001 by Kennedy and O’Hagan [193] and was originallyformulated to combine experimental measurements and computer model predictions in order to improvemodel predictions, even if not all important physical factors were taken into account in the computer model.Note that Kennedy and O’Hagan proposed two similarly-named and closely-related methods in the samepaper [193]: Bayesian calibration and Bayesian model calibration. While both of these methods use the samecombination of Gaussian process models, the goal of Bayesian calibration is to determine the parametersthat maximize the accuracy of a mathematical model (while including error correction) and the goal of BMCis to correct the errors of a given mathematical model.98BMC was originally proposed to handle a scenario like the following [193]: We have measurementdata, such as of the NOx concentrations in the atmosphere over the course of a year, and we would liketo predict these concentrations in future years. Unfortunately, the best mathematical model we have is apoor fit to the measurement data and we do not know of any physically-motivated ways to improve the fit.However, we do know that the model must correlate with the physical situation to some extent, given thatthe model is physically motivated. We then exploit this knowledge by (1) fitting a Gaussian process modelto the correlations between the predictions of the mathematical model and the measurement data and then(2) using the Gaussian process model to correct the predictions of the mathematical model.Kennedy and O’Hagan [193] showed how to perform such a correction and obtain information aboutthe uncertainty in the final prediction. One begins by training two Gaussian process models simultaneously,one to model the mathematical model and another to model the difference between the mathematical modeland the experimental data. We then have the following relationship:Y (·) = ρF (·)+G(·) , (5.1)where Y (·) is the final prediction, F (·) is the Gaussian process model that describes our mathematicalmodel, G(·) is the Gaussian process model that describes the difference between the mathematical modeland the experimental prediction, (·) represents the input parameters, and ρ is a tuning parameter.In principle, we can treat these two Gaussian process models in a fully Bayesian manner, such that wecan rigorously determine a final estimate of the uncertainty in the prediction and use the data to determineall hyperparameters, including ρ [193]. However, because this can be rather complex to implement, we firstimplement a simplified version of this algorithm. I discuss this simplified algorithm in the next section.5.4 Simplified Bayesian Model CalibrationIdeally, we would like to implement the full Bayesian form of BMC. However, for the present work, wehave implemented a simplified version that is both easier and faster to program. As the core idea of BMCis for the algorithm to learn the correlations between a mathematical model and experimental data, we havechosen to implement only this feature and neglect the rest. That is, instead of the full Bayesian treatmentwhere both the mean and standard deviations of the two Gaussian process models are properly taken into99account and the two Gaussian process models are trained simultaneously, we only take the means of theGaussian process models into account and train the two Gaussian process models individually.Thus, in our simplified Bayesian model calibration (sBMC) method we set ρ = 1 in Eqn. (5.1) and firstfit only F (·) to the numerical model. We then fit G(·) to χ(·) = y(·)− µF(·), where y(·) represents theexperimental data and µF is the mean of the Gaussian process F . Finally, we define the prediction of sBMCas Y (·) = µF(·) + µG(·), where µG is the mean of the Gaussian process G. While this is a rather crudeapproach and any final results should rely on the full BMC method, this sBMC method has the advantage ofsimplicity and the ability to still provide proof of principle for our ideas.Instead of working with a mathematical model and experimental data, as did Kennedy and O’Hagan[193], we will work with the data from poorly converged Monte Carlo calculations and well-convergedMonte Carlo calculations, respectively. As both kinds of calculations are physically motivated, there shouldbe a correlation between the inaccurate poorly converged data and the accurate well-converged data. Ourgoal is to learn these correlations at just a few points in parameter space and then use these correlations tocorrect the poorly converged calculations in the remainder of the parameter space of interest.To enact this sBMC correction scheme, we do the following:1. Perform several poorly-converged calculations that span the parameter space of interest2. Fit a Gaussian process model F(·) to these data points3. Perform a sparse set of well-converged calculations in the same parameter space4. Fit a Gaussian process model G(·) to the difference between µF(·) and the well-converged data points5. Calculate the final prediction Y (·) = µF(·)+µG(·)In this way, we obtain a Gaussian process model F(·) of the poorly-converged calculations and a Gaussianprocess model G(·) of the difference between µF(·) and the well-converged data points. Our final predic-tion consists of µF(·), which models the poorly converged results, corrected by µG(·), which models thecorrelations between the poorly- and well-converged calculations.1005.5 Proof of Principle DemonstrationWe now apply this sBMC correction scheme to the bold diagrammatic Monte Carlo method described byProkof’ev and Svistunov in Ref. [185]. One application of bold diagrammatic Monte Carlo that they describein Ref. [185] is to the calculation of the s-wave scattering length of the spherical potential barrier. That is, thepotential whose value is U0 for r< 1 and zero otherwise, with U0≥ 0. As this problem has a known analyticalsolution, it is a useful test bed for the initial proof-of-principle application of sBMC. Also, Prokof’ev hasmade the relevant Fortran code for this problem publicly available [1]. Using this publicly-available code,we calculate the scattering length as a function of the potential strength U0 and of the number of Monte Carlosteps included in the algorithm. We work with the number of Monte Carlo steps instead of the diagram orderas the calculation of the s-wave scattering length of the spherical potential barrier involves only two diagramsin this version of the bold diagrammatic Monte Carlo algorithm.While the details of the diagrammatic Monte Carlo algorithm are not needed for the present discussion,those interested in learning more about diagrammatic Monte Carlo can consult Prokof’ev’s excellent intro-duction to Monte Carlo methods [1]. Furthermore, Prokof’ev, Svistunov, and their co-authors have manyexcellent papers that describe diagrammatic Monte Carlo in significant detail for applications that includethe Hubbard model [73], the unitary Fermi gas [75, 194], Holstein polarons [195], Dirac liquids [196], andFermi-polarons [185, 186]. I find Refs. [73–75] particularly useful for understanding the method. Also,excellent introductions to many-body Feynman diagrams (on which diagrammatic Monte Carlo methodsare frequently based) can be found in Chapters 4 and 5 of Condensed Matter Field Theory by A. Altlandand B. Simons [92] and Chapters 3 and 7 of Quantum Theory of Many-Particle Systems by A. L. Fetter andJ. D. Walecka [197].We first calculate the s-wave scattering length a0 for the spherical potential barrier with 1×109 MonteCarlo steps for 100 values of U0 chosen via Latin hypercube sampling from the interval [10−3,30] [174]. Wethen calculate the scattering length with 2× 1010 Monte Carlo steps for ten values of U0 chosen via Latinhypercube sampling from the interval [10−3,20]. We compare the Monte Carlo calculations with the exactsolution for the scattering length of the attractive spherical potential barrier [198]:a0 = 1− tanh√2U0√2U0, (5.2)101where U0 ≥ 0. The relative absolute difference between the Monte Carlo calculations and the exact answer(5.2) are shown in Figure 5.1 for the low Monte Carlo step count results (blue) and the high step countresults (orange). It can be seen that the well-converged results show a more consistent low error than thepoorly converged results. While there are regions where the poorly-converged calculations appear to bemore accurate than the well-converged calculations, the convergence of these diagrammatic Monte Carlocalculations is not always monotonic.We then run the sBMC algorithm with the many low step count results and the few high step countresults. The relative absolute difference between the final prediction Y (·) and the exact result (5.2) is shownas the green curve in Figure 5.1. We see that the error of the sBMC result mostly resembles or improves uponthe error of the high step count results, at least within the range spanned by the high step count results. Wethus see how the computationally low-cost calculations are improved with comparatively few higher-cost,but more accurate, calculations. Clearly, however, this correction only occurs when interpolating betweenthe sparse, well-converged calculations; in regions outside the span of the well-converged calculations, theerror grows rapidly. Admittedly, there are regions where the error is higher than either the low or highstep count results, which shows the imperfect nature of these results. Overall, however, we can see animprovement that justifies further refinement and exploration of this approach.Importantly, we should improve the results by performing the full BMC algorithm, as opposed to thesBMC scheme used here. Particularly, the full BMC algorithm would provide an estimate of the uncertaintyin the final prediction. However, in the absence of easily usable and publically-accessible codes that im-plement the full BMC algorithm (at least in the Python programming language), the sBMC algorithm canbe used by researchers wishing to test the principles of the full BMC algorithm for their particular situationbefore investing the time to implement the full BMC algorithm. Any final results should rely on the fullBMC algorithm, however.As mentioned, the above results only provide a proof-of-principle result and are not yet a demonstratedalgorithmic advance. However, our proposal can be extended to situations where there are more dependentvariables and to situations where the exact answer is not known. In the absence of the exact result, it isimportant to determine ways to assess whether the BMC algorithm is indeed improving upon the poorlyconverged calculations. One simple way is to perform a few more well-converged calculations to compare10210 3 10 1 101 103Barrier Height, U010 510 410 310 210 1100Relative Absolute DifferenceFigure 5.1: The relative absolute difference between the numerical results and the exact scatteringlength (5.2) of a repulsive spherical barrier versus the height of the barrier U0. The numericalresults are obtained from a Monte Carlo calculation with 1× 109 Monte Carlo steps (blue), aMonte Carlo calculation with 2×1010 Monte Carlo steps (orange), and the simplified Bayesianmodel calibration result (green) discussed in Section 5.5. More details are provided in the text.to the BMC predictions. Finally, as the Gaussian process models can be modified to take into account theuncertainty in the data points (i.e. the models can be fit to the error bars, not just the data points) [45],the BMC algorithm could be used to incorporate the variance of calculated observables. Thus, we neednot throw away the statistical information generated by the Monte Carlo algorithm, but can include thisinformation when applying the full BMC algorithm.5.6 ConclusionWe have proposed applying Bayesian model calibration (BMC) to reduce the negative impact of criticalslowing down and the sign problem in Monte Carlo calculations, whether classical or quantum. We havedemonstrated an initial proof of principle calculation that we can use BMC to improve the accuracy ofpoorly converged Monte Carlo calculations with only a few well-converged calculations. In this demonstra-tion, we applied the simplified Bayesian model calibration (sBMC) algorithm to the problem of calculatingthe scattering length of the repulsive spherical barrier via diagrammatic Monte Carlo. This initial success103prompts us to implement the full BMC algorithm and to apply the algorithm to various other Monte Carlocalculations.For example, we could apply the BMC algorithm to more sophisticated diagrammatic Monte Carloschemes, such as those for the Hubbard model [73], angulons [199], or the Fermi-polaron problem [185,186]. In these cases, we may approach the problem by correcting low Feynman diagram order calculationswith a few high diagram order calculations. As another example, we could correct low bead number calcula-tions with a few high bead number calculations in path integral Monte Carlo calculations [71, 72]. And, wecould combine poorly and well-converged calculations for classical Monte Carlo calculations near a secondorder phase transition.Moreover, it appears that BMC may be useful for many types of calculation, not only Monte Carlocalculations. Particularly, if there is a calculation with a tuning parameter that governs both accuracy andcomputational cost, the BMC algorithm may be able to improve calculation efficiency by using a few compu-tationally high-cost, “well-tuned” calculations to correct several computationally low-cost, “poorly-tuned”calculations. Examples of tuning parameters include the bond dimension (or number of eigenstates retained)in density matrix renormalization group calculations [67–69], the number of baths in dynamical mean fieldtheory [70], and the step or mesh sizes in the integration of differential equations. Finally, we need notonly use data from a single type of calculation, but could combine data from different methods and/or an-alytical results in order to improve overall accuracy or interpolate between well-known limits. All of theseapplications are speculative, but warrant further investigation.1045.7 SummaryThis section contains a summary of the key points presented in Chapter 5:• Proposed applying Bayesian model calibration to improve the efficiency of Monte Carlo calculationsin regions where the sign problem or critical slowing down are problematic (Sections 5.2-5.4)• Provided a proof of principle demonstration that Bayesian model calibration can improve the ac-curacy of many poorly converged Monte Carlo calculations with the information from just a fewwell-converged calculations (Section 5.5)• Discussed other types of calculations where Bayesian model calibration may improve calculationefficiency (Section 5.6)105Chapter 6Conclusion6.1 Summary of ResultsIn this thesis, I have explored novel manifestations and applications of quantum interference in complex sys-tems and developed new approaches to study complex quantum systems. We have uncovered new physics,modelled a new experimental system and proposed a way to improve the convergence of Monte Carlo meth-ods.In Chapter 2, we showed that particles whose long-range hopping amplitudes are isotropic can localize indisordered three dimensional lattices. We also showed that cooperative shielding, which mitigates the effectsof long-range terms of a Hamiltonian [30, 119], occurs in three-dimensional lattices for these particles. Incontrast, if the long-range hopping amplitudes have the same anisotropy as the interaction between twodipoles aligned along a single axis, we see the absence of cooperative shielding and no concrete evidenceof localized states. We have shown that the existence of delocalized, non-ergodic extended, or localizedstates in the presence of disorder depends on both the hopping range exponents α and the anisotropy of thelong-range hopping amplitudes.In Chapter 3, we proposed using molecules instead of 3He atoms in surface spin echo experiments.That is, instead of using the nuclear spin states of 3He to probe sample surfaces and their adsorbates, weproposed using the molecular hyperfine states of generic closed shell molecules, whose additional degreesof freedom may provide us with more insight into surface properties and dynamics. We also developed afully quantum model of the experiment that involved re-deriving the transfer matrix formalism to include106the internal degrees of freedom. This formalism is needed to address the inverse scattering problem and tointerpret these molecular hyperfine interferometry experiments. We also applied our model to the case ofoH2 and showed that the experimental signal does indeed include information about the molecule-surfacescattering matrix elements.In Chapter 4, we applied Bayesian optimization to the inverse scattering problem for these experiments.Bayesian optimization models the function to be minimized with a Gaussian process model and uses themodel, along with its uncertainty, to intelligently sample the function to be minimized. Such a samplingprocedure significantly reduces the number of function evaluations required for global minimization. UsingBayesian optimization, we determined a scattering matrix defined by three parameters from data generatedfrom known target parameters. We also demonstrated that the approach is significantly more efficient thangrid searches with a comparable resolution.In Chapter 5, we proposed using Bayesian model calibration to improve the efficiency of Monte Carloalgorithms in regions where the sign problem or critical slowing down inhibit the calculation of observables.Bayesian model calibration learns the correlations between several poorly converged calculations and a fewwell-converged calculations and uses these correlations to correct the computationally low-cost, but poorlyconverged results. We provided a proof of principle demonstration by applying a simplified version ofBayesian model calibration to a simple diagrammatic Monte Carlo calculation of the scattering length of aspherical potential barrier.In conclusion, this thesis has improved our understanding of the quantum localization of particles withlong-range hopping, provided the theoretical framework needed to interpret molecular hyperfine interfer-ometry experiments, began to address the inverse scattering problem for these experiments by applyingBayesian optimization, and proposed an application of Bayesian model calibration to Monte Carlo calcula-tions to reduce the impact of the sign problem and critical slowing down.6.2 OutlookThe results of this thesis motivate many possible future directions, for both experimental and theoreticalwork. This section contains a list, necessarily incomplete, of these possible directions.1076.2.1 Anderson Localization with Long-range Hopping• Our analytical results appear to imply that cooperative shielding would also exist in two dimensionalsystems, but it is important to actually confirm this.• Our work examined only two possible forms of angular dependence of the hopping amplitudes. It isimportant to examine how different forms of anisotropy impact cooperative shielding and localization.Different lattice geometries may also impact cooperative shielding and localization. Furthermore, theimpact of the lattice geometry may also be a function of the anisotropy and vice versa.• In a recent experiment, researchers could tune the anisotropy of interactions between atoms dressedby Rydberg excitations, for atoms that are in a two dimensional optical or magnetic lattice [23]. Thiswould be a useful test bed for the ideas presented in Chapter 2.• The anisotropic form studied in Chapter 2 was that of dipoles aligned along a principal axis of thelattice. This form arises from a model of polar molecules in an optical lattice that are aligned by anexternal electric field. One may be able to change the form of the anisotropy, or at least its impact oncooperative shielding and localization, by changing the angle of the aligning electric field.6.2.2 Molecular Hyperfine Interferometry• Our plan is to expand the fully quantum formalism to also account for the dynamics of surfaces thatare on the molecule-surface interaction time scales. That is, we plan to incorporate a time-dependentscattering matrix, which may require the numerical evaluation of the energy integrals and/or couplingbetween transport channels of differing energy.• Our formalism can help design experimental protocols that produce a signal which may more clearlyelucidate particular features of the scattering matrix or which may be easier to interpret.• We plan to continue applying Bayesian optimization to the inverse scattering problem of molecularhyperfine interferometry experiments. We will work with different classes of scattering matrices andpotentially add energy-dependence to the scattering matrices. Ideally, we will be able to determine all81+ degrees of freedom of the scattering matrix from the scattering data, but even just finding classes108of matrices compatible with the experimental data will help to elucidate the symmetries at play in themolecule-surface interactions.• Other types of optimization algorithms may help address this inverse scattering problem, such asgradient-based methods [175, 176].• We may be able to approach the inverse scattering problem from a more Bayesian perspective. That is,we could begin with a prior over the scattering matrix elements and propagate this uncertainty throughthe fully quantum model we have developed to get an uncertainty in the output signal. Then, we couldcompare this result to the experimental data and use Bayesian optimization to vary the parameters ofthe prior until the generated signal and its uncertainty match the experimental data and its uncertainty.6.2.3 Monte Carlo Enhanced with Bayesian Model Calibration• We plan to implement the full Bayesian model calibration algorithm, which we anticipate will providebetter results than those obtained from the simplified Bayesian model calibration algorithm discussedin Chapter 5• We plan to apply the full Bayesian model calibration algorithm to many different types of MonteCarlo calculations and thus show the general applicability of the approach• We plan to apply Bayesian model calibration to Monte Carlo calculations in order to push into param-eter regimes of important models where the convergence is currently too slow to calculate observableswith small uncertainties, such as near phase transitions• Bayesian model calibration may also improve the efficiency of other types of calculations that havea tuning parameter that governs both accuracy and computational cost, such as the density matrixrenormalization group [67–69] and dynamical mean field theory [70].109Bibliography[1] N. V. Prokof’ev, Monte Carlo Methods Theory, http://mcwa.csi.cuny.edu/umass/lectures.html,Accessed: 2019-05-10. → pages vii, 101[2] I. Bloch, J. Dalibard, and S. Nascimbene, Quantum simulations with ultracold quantum gases, Nat.Phys. 8, 267 (2012). → page 1[3] R. Grimm, M. Weidemu¨ller, and Y. B. Ovchinnikov, Optical dipole traps for neutral atoms, inAdvances In Atomic, Molecular, and Optical Physics (Elsevier 2000), volume 42, pp. 95–170. →page 1[4] N. Friedman, A. Kaplan, and N. Davidson, Dark optical traps for cold atoms, in B. Bederson andH. Walther, eds., Advances In Atomic, Molecular, and Optical Physics (Academic Press 2002),volume 48, pp. 99–151.[5] V. I. Balykin, V. G. Minogin, and V. S. Letokhov, Electromagnetic trapping of cold atoms, Rep.Prog. Phys. 63, 1429 (2000). → page 1[6] J. A. Rushton, M. Aldous, and M. D. Himsworth, Contributed Review: The feasibility of a fullyminiaturized magneto-optical trap for portable ultracold quantum technology, Rev. Sci. Instrum. 85,121501 (2014). → page 1[7] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357,995 (2017). → page 1[8] I. Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1, 23 (2005).[9] O. Morsch and M. Oberthaler, Dynamics of Bose-Einstein condensates in optical lattices, Rev. Mod.Phys. 78, 179 (2006). → page 1[10] Y. Wang, P. Surendran, S. Jose, T. Tran, I. Herrera, S. Whitlock, R. McLean, A. Sidorov, andP. Hannaford, Magnetic lattices for ultracold atoms and degenerate quantum gases, Sci. Bull. 61,1097 (2016). → page 1[11] M. R. Tarbutt, Laser cooling of molecules, Contemp. Phys. 59, 356 (2018). → page 1[12] B. Gadway and B. Yan, Strongly interacting ultracold polar molecules, J. Phys. B-At. Mol. Opt.Phys. 49, 152002 (2016).[13] D. McCarron, Laser cooling and trapping molecules, J. Phys. B-At. Mol. Opt. Phys. 51, 212001(2018).110[14] M. A. Baranov, Theoretical progress in many-body physics with ultracold dipolar gases, Phys.Rep.-Rev. Sec. Phys. Lett. 464, 71 (2008). → page 1[15] P. W. Courteille, V. S. Bagnato, and V. I. Yukalov, Bose-Einstein condensation of trapped atomicgases, Laser Phys. 11, 659 (2001). → page 1[16] A. L. Fetter, Rotating trapped Bose-Einstein condensates, Rev. Mod. Phys. 81, 647 (2009). → page1[17] M. Greiner, C. A. Regal, and D. S. Jin, Emergence of a molecular Bose-Einstein condensate from aFermi gas, Nature 426, 537 (2003). → page 1[18] T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy, M. Teichmann, L. Tarruell,S. Kokkelmans, and C. Salomon, Experimental study of the BEC-BCS crossover region in lithium-6,Phys. Rev. Lett. 93, 050401 (2004). → page 1[19] G. Roati, C. D’Errico, L. Fallani, M. Fattori, C. Fort, M. Zaccanti, G. Modugno, M. Modugno, andM. Inguscio, Anderson localization of a non-interacting Bose–Einstein condensate, Nature 453, 895(2008). → pages 1, 8[20] J. Billy, V. Josse, Z. Zuo, A. Bernard, B. Hambrecht, P. Lugan, D. Cle´ment, L. Sanchez-Palencia,P. Bouyer, and A. Aspect, Direct observation of Anderson localization of matter waves in acontrolled disorder, Nature 453, 891 (2008). → pages 1, 8[21] R. Islam, C. Senko, W. C. Campbell, S. Korenblit, J. Smith, A. Lee, E. E. Edwards, C.-C. J. Wang,J. K. Freericks, and C. Monroe, Emergence and frustration of magnetism with variable-rangeinteractions in a quantum simulator, Science 340, 583 (2013). → page 1[22] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, Quasiparticleengineering and entanglement propagation in a quantum many-body system, Nature 511, 202(2014). → page 1[23] A. W. Glaetzle, M. Dalmonte, R. Nath, C. Gross, I. Bloch, and P. Zoller, Designing frustratedquantum magnets with laser-dressed Rydberg atoms, Phys. Rev. Lett. 114, 173002 (2015). → pages1, 108[24] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958). →pages 1, 7, 8, 12, 13, 16, 25[25] L. S. Levitov, Absence of localization of vibrational modes due to dipole-dipole interaction,Europhys. Lett. 9, 83 (1989). → pages 13, 14, 16, 34, 39[26] L. S. Levitov, Delocalization of vibrational modes caused by electric dipole interaction, Phys. Rev.Lett. 64, 547 (1990). → pages 1, 7, 13, 16, 34, 39[27] A. L. Burin and L. A. Maksimov, Localization and delocalization of particles in disordered latticewith tunneling amplitude with R−3 decay, Pis’ma Zh. Eksp. Teor. Fiz. 50, 304 (1989), [J. Exp.Theor. Phys. Lett. 50, 338 (1989)]; Online: http://www.jetpletters.ac.ru/ps/1129/article 17116.shtml.→ pages 1, 16, 21, 39111[28] X. Deng, V. Kravtsov, G. Shlyapnikov, and L. Santos, Duality in power-law localization indisordered one-dimensional systems, Phys. Rev. Lett. 120, 110602 (2018). → pages 2, 7, 16, 21, 25[29] A. Ossipov, Anderson localization on a simplex, J. Phys. A 46, 105001 (2013). → pages17, 23, 29, 30[30] G. L. Celardo, R. Kaiser, and F. Borgonovi, Shielding and localization in the presence of long-rangehopping, Phys. Rev. B 94, 144206 (2016). → pages 2, 16, 17, 18, 21, 22, 23, 25, 27, 30, 38, 106[31] ACME Collaboration, Improved limit on the electric dipole moment of the electron, Nature 562, 355(2018). → page 2[32] A. P. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison, and J. Ellis, Helium-3 spin-echo:Principles and application to dynamics at surfaces, Prog. Surf. Sci. 84, 323 (2009). → page 2[33] G. Alexandrowicz, P. R. Kole, E. Y. M. Lee, H. Hedgeland, R. Ferrando, A. P. Jardine, W. Allison,and J. Ellis, Observation of uncorrelated microscopic motion in a strongly interacting adsorbatesystem, J. Am. Chem. Soc. 130, 6789 (2008). → pages 2, 43, 74[34] G. Alexandrowicz, A. P. Jardine, H. Hedgeland, W. Allison, and J. Ellis, Onset of 3D collectivesurface diffusion in the presence of lateral interactions: Na/Cu(001), Phys. Rev. Lett. 97, 156103(2006).[35] G. Alexandrowicz, A. P. Jardine, P. Fouquet, S. Dworski, W. Allison, and J. Ellis, Observation ofmicroscopic CO dynamics on Cu(001) using 3He spin-echo spectroscopy, Phys. Rev. Lett. 93,156103 (2004).[36] H. Hedgeland, M. Sacchi, P. Singh, A. J. McIntosh, A. P. Jardine, G. Alexandrowicz, D. J. Ward,S. J. Jenkins, W. Allison, and J. Ellis, Mass transport in surface diffusion of van der Waals bondedsystems: Boosted by rotations?, J. Phys. Chem. Lett. 7, 4819 (2016).[37] O. Godsi, G. Corem, T. Kravchuk, C. Bertram, K. Morgenstern, H. Hedgeland, A. P. Jardine,W. Allison, J. Ellis, and G. Alexandrowicz, How atomic steps modify diffusion and inter-adsorbateforces: Empirical evidence from hopping dynamics in Na/Cu(115), J. Phys. Chem. Lett. 6, 4165(2015). → page 43[38] P. Fouquet, H. Hedgeland, A. Jardine, G. Alexandrowicz, W. Allison, and J. Ellis, Measurements ofmolecule diffusion on surfaces using neutron and helium spin echo, Physica B 385, 269 (2006).[39] A. P. Jardine, G. Alexandrowicz, H. Hedgeland, R. D. Diehl, W. Allison, and J. Ellis, Vibration anddiffusion of Cs atoms on Cu(001), J. Phys. Condens. Matter 19, 305010 (2007). → pages 2, 43, 74[40] P. R. Kole, A. P. Jardine, H. Hedgeland, and G. Alexandrowicz, Measuring surface phonons with a 3He spin echo spectrometer: a two-dimensional approach, J. Phys. Condens. Matter 22, 304018(2010). → pages 2, 43[41] G. Alexandrowicz and A. P. Jardine, Helium spin-echo spectroscopy: studying surface dynamicswith ultra-high-energy resolution, J. Phys. Condens. Matter 19, 305001 (2007). → pages 2, 43, 47[42] E. M. McIntosh, P. R. Kole, M. El-Batanouny, D. M. Chisnall, J. Ellis, and W. Allison, Measurementof the phason dispersion of misfit dislocations on the Au(111) surface, Phys. Rev. Lett. 110, 086103(2013). → pages 2, 43112[43] E. M. McIntosh, J. Ellis, A. P. Jardine, P. Licence, R. G. Jones, and W. Allison, Probing liquidbehaviour by helium atom scattering: surface structure and phase transitions of an ionic liquid onAu(111), Chem. Sci. 5, 667 (2014). → pages 2, 43[44] O. Godsi, G. Corem, Y. Alkoby, J. T. Cantin, R. V. Krems, M. F. Somers, J. Meyer, G.-J. Kroes,T. Maniv, and G. Alexandrowicz, A general method for controlling and resolving rotationalorientation of molecules in molecule-surface collisions, Nat. Comm. 8, 15357 (2017). → pages2, 42, 43, 44, 45, 46, 47, 48, 51, 52, 54, 55, 60, 63, 66, 67, 68, 70, 74, 75, 131[45] R. V. Krems, Bayesian machine learning for quantum molecular dynamics, Phys. Chem. Chem.Phys. 21, 13392 (2019). → pages 3, 64, 76, 79, 82, 85, 103[46] B. Sanchez-Lengeling and A. Aspuru-Guzik, Inverse molecular design using machine learning:Generative models for matter engineering, Science 361, 360 (2018).[47] A. Radovic, M. Williams, D. Rousseau, M. Kagan, D. Bonacorsi, A. Himmel, A. Aurisano,K. Terao, and T. Wongjirad, Machine learning at the energy and intensity frontiers of particlephysics, Nature 560, 41 (2018).[48] D. M. Dimiduk, E. A. Holm, and S. R. Niezgoda, Perspectives on the impact of machine learning,deep learning, and artificial intelligence on materials, processes, and structures engineering, Integr.Mater. Manuf. Innov. 7, 157 (2018).[49] K. T. Butler, D. W. Davies, H. Cartwright, O. Isayev, and A. Walsh, Machine learning for molecularand materials science, Nature 559, 547 (2018).[50] V. Dunjko and H. J. Briegel, Machine learning & artificial intelligence in the quantum domain: Areview of recent progress, Rep. Prog. Phys. 81, 074001 (2018). → page 3[51] J. Behler, Neural network potential-energy surfaces in chemistry: a tool for large-scale simulations,Phys. Chem. Chem. Phys. 13, 17930 (2011). → page 3[52] J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nat. Phys. 13, 431 (2017). →page 3[53] K. Ch’ng, J. Carrasquilla, R. G. Melko, and E. Khatami, Machine learning phases of stronglycorrelated fermions, Phys. Rev. X 7, 031038 (2017). → page 3[54] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-networkquantum state tomography, Nat. Phys. 14, 447 (2018). → page 3[55] P. Ponte and R. G. Melko, Kernel methods for interpretable machine learning of order parameters,Phys. Rev. B 96, 205146 (2017). → page 3[56] J. Cui and R. V. Krems, Efficient non-parametric fitting of potential energy surfaces for polyatomicmolecules with Gaussian processes, J. Phys. B: At. Mol. Opt. Phys. 49, 224001 (2016). → page 3[57] A. Kamath, R. A. Vargas-Herna´ndez, R. V. Krems, T. Carrington, and S. Manzhos, Neural networksvs Gaussian process regression for representing potential energy surfaces: A comparative study of fitquality and vibrational spectrum accuracy, J. Chem. Phys. 148, 241702 (2018).113[58] A. Christianen, T. Karman, R. A. Vargas-Herna´ndez, G. C. Groenenboom, and R. V. Krems,Six-dimensional potential energy surface for NaK–NaK collisions: Gaussian process representationwith correct asymptotic form, J. Chem. Phys. 150, 064106 (2019). → page 3[59] R. A. Vargas-Herna´ndez, J. Sous, M. Berciu, and R. V. Krems, Extrapolating quantum observableswith machine learning: Inferring multiple phase transitions from properties of a single phase, Phys.Rev. Lett. 121, 255702 (2018). → page 3[60] J. Cui and R. V. Krems, Gaussian process model for collision dynamics of complex molecules, Phys.Rev. Lett. 115, 073202 (2015). → pages 3, 4, 98[61] J. Cui, Z. Li, and R. V. Krems, Gaussian process model for extrapolation of scattering observablesfor complex molecules: From benzene to benzonitrile, J. Chem. Phys. 143, 154101 (2015).[62] D. Vieira and R. V. Krems, Rate constants for fine-structure excitations in O–H collisions with errorbars obtained by machine learning, Astrophys. J. 835, 255 (2017). → page 3[63] R. A. Vargas-Herna´ndez, Bayesian optimization for tuning and selecting hybrid-density functionals,arXiv:1903.10678 (2019). → page 3[64] R. A. Vargas-Herna´ndez, Y. Guan, D. H. Zhang, and R. V. Krems, Bayesian optimization for theinverse scattering problem in quantum reaction dynamics, New J. Phys. 21, 022001 (2019). →pages 3, 64, 76, 79, 90[65] W. Xu, W. R. McGehee, W. N. Morong, and B. DeMarco, Bad-metal relaxation dynamics in a Fermilattice gas, Nat. Commun. 10, 1588 (2019). → page 4[66] P. T. Brown, D. Mitra, E. Guardado-Sanchez, R. Nourafkan, A. Reymbaut, C.-D. He´bert,S. Bergeron, A.-M. S. Tremblay, J. Kokalj, D. A. Huse, P. Schauß, and W. S. Bakr, Bad metallictransport in a cold atom Fermi-Hubbard system, Science 363, 379 (2019). → page 4[67] S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69,2863 (1992). → pages 4, 104, 109[68] S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48,10345 (1993).[69] S. Al-Assam, S. R. Clark, and D. Jaksch, The tensor network theory library, J. Stat. Mech. 2017,093102 (2017). → pages 4, 104, 109[70] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of stronglycorrelated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996). →pages 4, 104, 109[71] R. Rodriguez-Cantano, T. Gonzalez-Lezana, and P. Villarreal, Path integral Monte Carloinvestigations on doped helium clusters, Int. Rev. Phys. Chem. 35, 37 (2016). → pages 4, 96, 104[72] J. Grotendorst, ed., Quantum simulations of complex many-body systems: from theory to algorithms:winter school, 25 February - 1 March 2002, Rolduc Conference Centre, Kerkrade, the Netherlands,number 10 in NIC series (NIC-Secretariat, Ju¨lich 2002). → pages 4, 104114[73] K. Van Houcke, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov, Diagrammatic Monte Carlo, Phys.Procedia 6, 95 (2010). → pages 4, 96, 101, 104[74] K. Van Houcke, E. Kozik, N. V. Prokof’ev, and B. V. Svistunov, Diagrammatic Monte Carlo,arXiv:0802.2923 (2008). → page 97[75] K. Van Houcke, F. Werner, N. V. Prokof’ev, and B. V. Svistunov, Bold diagrammatic Monte Carlofor the resonant Fermi gas, arXiv:1305.3901 (2013). → pages 4, 101[76] J. P. F. LeBlanc, A. E. Antipov, F. Becca, I. W. Bulik, G. K.-L. Chan, C.-M. Chung, Y. Deng,M. Ferrero, T. M. Henderson, C. A. Jime´nez-Hoyos, E. Kozik, X.-W. Liu, A. J. Millis, N. V.Prokof’ev, M. Qin, G. E. Scuseria, H. Shi, B. V. Svistunov, L. F. Tocchio, I. S. Tupitsyn, S. R. White,S. Zhang, B.-X. Zheng, Z. Zhu, and E. Gull, Solutions of the two-dimensional Hubbard model:Benchmarks and results from a wide range of numerical algorithms, Phys. Rev. X 5, 041041 (2015).→ page 4[77] L. Wang, Exploring cluster Monte Carlo updates with Boltzmann machines, Phys. Rev. E 96,051301(R) (2017). → pages 4, 98[78] L. Huang and L. Wang, Accelerated Monte Carlo simulations with restricted Boltzmann machines,Phys. Rev. B 95, 035105 (2017). → pages 4, 98[79] N. W. Ashcroft and N. D. Mermin, Solid State Physics (Holt, Rinehart and Winston, New York1976). → page 6[80] T. Dauxois, ed., Dynamics and Thermodynamics of Systems With Long-Range Interactions, number602 in Lecture notes in physics (Springer, Berlin; New York 2002). → page 7[81] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Oxford University Press, Oxford,United Kingdom 2017), 2nd edition. → page 7[82] A. Z. Genack and A. A. Chabanov, Signatures of photon localization, J. Phys. A: Math. Gen. 38,10465 (2005). → page 8[83] T. Schwartz, G. Bartal, S. Fishman, and M. Segev, Transport and Anderson localization indisordered two-dimensional photonic lattices, Nature 446, 52 (2007). → page 8[84] M. V. Berry and S. Klein, Transparent mirrors: rays, waves and localization, Eur. J. Phys. 18, 222(1997). → page 8[85] R. L. Weaver, Anderson localization of ultrasound, Wave Motion 12, 129 (1990). → page 8[86] H. Hu, A. Strybulevych, J. H. Page, S. E. Skipetrov, and B. A. van Tiggelen, Localization ofultrasound in a three-dimensional elastic network, Nat. Phys. 4, 945 (2008).[87] Y. Ye, M. Ke, J. Feng, M. Wang, C. Qiu, and Z. Liu, Transversal Anderson localization of sound inacoustic waveguide arrays, J. Phys.: Condens. Matter 27, 155402 (2015). → page 8[88] J. Chabe´, G. Lemarie´, B. Gre´maud, D. Delande, P. Szriftgiser, and J. C. Garreau, Experimentalobservation of the Anderson metal-insulator transition with atomic matter waves, Phys. Rev. Lett.101, 255702 (2008). → page 8115[89] S. Karbasi, R. J. Frazier, K. W. Koch, T. Hawkins, J. Ballato, and A. Mafi, Image transport througha disordered optical fibre mediated by transverse Anderson localization, Nat. Commun. 5, 3362(2014). → page 8[90] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Scaling theory oflocalization: Absence of quantum diffusion in two dimensions, Phys. Rev. Lett. 42, 673 (1979). →pages 9, 12, 13[91] C. A. Mu¨ller and D. Delande, Disorder and interference: localization phenomena, arXiv:1005.0915(2010). → pages 9, 10, 11, 12[92] A. Altland and B. Simons, Condensed Matter Field Theory (Cambridge University Press,Cambridge; New York 2010), 2nd edition. → pages 12, 82, 101[93] D. H. Dunlap, K. Kundu, and P. Phillips, Absence of localization in certain statically disorderedlattices in any spatial dimension, Phys. Rev. B 40, 10999 (1989). → page 12[94] D. H. Dunlap, H.-L. Wu, and P. W. Phillips, Absence of localization in a random-dimer model, Phys.Rev. Lett. 65, 88 (1990).[95] D. H. Dunlap and P. Phillips, Superdiffusion in structurally disordered systems in any spatialdimension, J. Chem. Phys. 92, 6093 (1990).[96] J. C. Flores, Transport in models with correlated diagonal and off-diagonal disorder, J. Phys.Condens. Matter 1, 8471 (1989). → page 12[97] K. Slevin and T. Ohtsuki, Corrections to scaling at the Anderson transition, Phys. Rev. Lett. 82, 382(1999). → page 13[98] G. Lemarie´, B. Gre´maud, and D. Delande, Universality of the Anderson transition with thequasiperiodic kicked rotor, EPL 87, 37007 (2009). → page 13[99] D. E. Logan and P. G. Wolynes, Localizability and dephasing of dipolar excitons in topologicallydisordered systems, J. Chem. Phys. 87, 7199 (1987). → page 13[100] M. Mladenovic´ and N. Vukmirovic´, Charge carrier localization and transport in organicsemiconductors: Insights from atomistic multiscale simulations, Adv. Funct. Mater. 25, 1915 (2015).→ pages 13, 37[101] A. Frisch, M. Mark, K. Aikawa, S. Baier, R. Grimm, A. Petrov, S. Kotochigova, G. Que´me´ner,M. Lepers, O. Dulieu, and F. Ferlaino, Ultracold dipolar molecules composed of strongly magneticatoms, Phys. Rev. Lett. 115, 203201 (2015). → page 13[102] K.-K. Ni, S. Ospelkaus, M. H. G. d. Miranda, A. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova,P. S. Julienne, D. S. Jin, and J. Ye, A high phase-space-density gas of polar molecules, Science 322,231 (2008).[103] M. H. G. de Miranda, A. Chotia, B. Neyenhuis, D. Wang, G. Que´me´ner, S. Ospelkaus, J. L. Bohn,J. Ye, and D. S. Jin, Controlling the quantum stereodynamics of ultracold bimolecular reactions,Nat. Phys. 7, 502 (2011).116[104] A. Chotia, B. Neyenhuis, S. A. Moses, B. Yan, J. P. Covey, M. Foss-Feig, A. M. Rey, D. S. Jin, andJ. Ye, Long-lived dipolar molecules and Feshbach molecules in a 3D optical lattice, Phys. Rev. Lett.108, 080405 (2012). → page 37[105] B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye,Observation of dipolar spin-exchange interactions with lattice-confined polar molecules, Nature501, 521 (2013). → pages 19, 37[106] K. R. Hazzard, B. Gadway, M. Foss-Feig, B. Yan, S. A. Moses, J. P. Covey, N. Y. Yao, M. D. Lukin,J. Ye, D. S. Jin, and A. M. Rey, Many-body dynamics of dipolar molecules in an optical lattice,Phys. Rev. Lett. 113, 195302 (2014). → page 13[107] F. Wu¨rthner, T. E. Kaiser, and C. R. Saha-Mo¨ller, J-aggregates: From serendipitous discovery tosupramolecular engineering of functional dye materials, Angew. Chem. Int. Ed. 50, 3376 (2011). →page 13[108] G. D. Scholes, G. R. Fleming, A. Olaya-Castro, and R. van Grondelle, Lessons from nature aboutsolar light harvesting, Nat. Chem. 3, 763 (2011). → page 13[109] A. F. Fidler, V. P. Singh, P. D. Long, P. D. Dahlberg, and G. S. Engel, Dynamic localization ofelectronic excitation in photosynthetic complexes revealed with chiral two-dimensionalspectroscopy, Nat. Commun. 5, 3286 (2014). → page 13[110] F. Robicheaux and N. M. Gill, Effect of random positions for coherent dipole transport, Phys. Rev.A 89, 053429 (2014). → page 13[111] O. Mu¨lken, A. Blumen, T. Amthor, C. Giese, M. Reetz-Lamour, and M. Weidemu¨ller, Survivalprobabilities in coherent exciton transfer with trapping, Phys. Rev. Lett. 99, 090601 (2007).[112] H. Schempp, G. Gu¨nter, S. Wu¨ster, M. Weidemu¨ller, and S. Whitlock, Correlated exciton transportin Rydberg-dressed-atom spin chains, Phys. Rev. Lett. 115, 093002 (2015). → page 13[113] L. S. Levitov, Anderson model for localization with long-range transfer amplitudes: Delocalizationof vibrational modes caused by their electric dipole interaction, Ann. N. Y. Acad. Sci. 581, 95(1990). → pages 13, 16, 34, 39[114] N. Yao, C. Laumann, S. Gopalakrishnan, M. Knap, M. Mu¨ller, E. Demler, and M. Lukin, Many-bodylocalization in dipolar systems, Phys. Rev. Lett. 113, 243002 (2014). → pages 16, 39[115] A. Burin, The luminescence of excitons in the disordered media, Fiz. Nizk. Temp. 20, 454 (1994). →page 16[116] J. T. Cantin, T. Xu, and R. V. Krems, Effect of the anisotropy of long-range hopping on localizationin three-dimensional lattices, Phys. Rev. B 98, 014204 (2018). → page 16[117] S. Bera, G. De Tomasi, I. M. Khaymovich, and A. Scardicchio, Return probability for the Andersonmodel on the random regular graph, Phys. Rev. B 98, 134205 (2018). → page 16[118] R. M. Nandkishore and S. Sondhi, Many-body localization with long-range interactions, Phys. Rev.X 7, 041021 (2017). → pages 16, 21117[119] L. F. Santos, F. Borgonovi, and G. L. Celardo, Cooperative shielding in many-body systems withlong-range interaction, Phys. Rev. Lett. 116, 250402 (2016). → pages 16, 17, 18, 21, 38, 106[120] T. Xu and R. V. Krems, Quantum walk and Anderson localization of rotational excitations indisordered ensembles of polar molecules, New J. Phys. 17, 065014 (2015). → page 20[121] L. J. Vasquez, A. Rodriguez, and R. A. Ro¨mer, Multifractal analysis of the metal-insulatortransition in the three-dimensional Anderson model. I. Symmetry relation under typical averaging,Phys. Rev. B 78, 195106 (2008). → pages 21, 34[122] A. Rodriguez, L. J. Vasquez, and R. A. Ro¨mer, Multifractal analysis of the metal-insulatortransition in the three-dimensional Anderson model. II. Symmetry relation under ensembleaveraging, Phys. Rev. B 78, 195107 (2008). → pages 21, 34[123] P. Van Mieghem, Graph Spectra for Complex Networks (Cambridge University Press, Cambridge2011). → page 23[124] X. Deng, B. Altshuler, G. Shlyapnikov, and L. Santos, Quantum Levy flights and multifractality ofdipolar excitations in a random system, Phys. Rev. Lett. 117, 020401 (2016). → pages 28, 37[125] A. De Luca, B. Altshuler, V. Kravtsov, and A. Scardicchio, Anderson localization on the Bethelattice: Nonergodicity of extended states, Phys. Rev. Lett. 113, 046806 (2014). → pages 28, 34[126] A. D. Mirlin, Statistics of energy levels and eigenfunctions in disordered systems, Phys. Rep. 326,259 (2000). → page 28[127] M. Janssen, Statistics and scaling in disordered mesoscopic electron systems, Phys. Rep. 295, 1(1998). → page 28[128] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev.B 75, 155111 (2007). → page 34[129] Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive levelspacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013). → page 34[130] B. I. Shklovskii, B. Shapiro, B. R. Sears, P. Lambrianides, and H. B. Shore, Statistics of spectra ofdisordered systems near the metal-insulator transition, Phys. Rev. B 47, 11487 (1993). → page 34[131] L. J. Root and J. L. Skinner, Localization phase diagram for the energetically and substitutionallydisordered Anderson/quantum percolation model, J. Chem. Phys. 89, 3279 (1988). → pages 36, 38[132] I. C. Gerber and X. Marie, Dependence of band structure and exciton properties of encapsulatedWSe2 monolayers on the hBN-layer thickness, Phys. Rev. B 98, 245126 (2018). → page 40[133] P. A. Nosov, I. M. Khaymovich, and V. E. Kravtsov, Correlation-induced localization, Phys. Rev. B99, 104203 (2019). → page 40[134] A. P. Jardine, H. Hedgeland, G. Alexandrowicz, W. Allison, and J. Ellis, Helium-3 spin-echo:Principles and application to dynamics at surfaces, Prog. Surf. Sci. 84, 323 (2009). → page 43[135] A. P. Jardine, G. Alexandrowicz, H. Hedgeland, W. Allison, and J. Ellis, Studying the microscopicnature of diffusion with helium-3 spin-echo, Phys. Chem. Chem. Phys. 11, 3355 (2009). → page 43118[136] G. Corem, P. R. Kole, J. Zhu, T. Kravchuk, J. R. Manson, and G. Alexandrowicz, Ordered H2Ostructures on a weakly interacting surface: A helium diffraction study of H2O/Au(111), J. Phys.Chem. C 117, 23657 (2013). → page 43[137] P. R. Kole, H. Hedgeland, A. P. Jardine, W. Allison, J. Ellis, and G. Alexandrowicz, Probing thenon-pairwise interactions between CO molecules moving on a Cu(111) surface, J. Phys. Condens.Matter 24, 104016 (2012). → page 43[138] H. Hedgeland, P. R. Kole, H. R. Davies, A. P. Jardine, G. Alexandrowicz, W. Allison, J. Ellis,G. Fratesi, and G. P. Brivio, Surface dynamics and friction of K/Cu(001) characterized by helium-3spin-echo and density functional theory, Phys. Rev. B 80, 125426 (2009). → page 43[139] H. Hedgeland, P. Fouquet, A. P. Jardine, G. Alexandrowicz, W. Allison, and J. Ellis, Measurement ofsingle-molecule frictional dissipation in a prototypical nanoscale system, Nat. Phys. 5, 561 (2009).→ page 43[140] A. P. Jardine, S. Dworski, P. Fouquet, G. Alexandrowicz, D. J. Riley, G. Y. H. Lee, J. Ellis, andW. Allison, Ultrahigh-resolution spin-echo measurement of surface potential energy landscapes,Science 304, 1790 (2004). → page 43[141] D. J. Riley, A. P. Jardine, S. Dworski, G. Alexandrowicz, P. Fouquet, J. Ellis, and W. Allison, Arefined He-LiF(001) potential from selective adsorption resonances measured with high-resolutionhelium spin-echo spectroscopy, J. Chem. Phys. 126, 104702 (2007).[142] D. J. Riley, A. P. Jardine, G. Alexandrowicz, H. Hedgeland, J. Ellis, and W. Allison, Analysis andrefinement of the Cu(001)c(2X2)CO-He potential using 3He selective adsorption resonances, J.Chem. Phys. 128, 154712 (2008). → page 43[143] G. Fratesi, G. Alexandrowicz, M. I. Trioni, G. P. Brivio, and W. Allison, Crucial electroniccontributions to measures of surface diffusion by He atom scattering, Phys. Rev. B 77, 235444(2008). → page 43[144] A. P. Jardine, H. Hedgeland, D. Ward, Y. Xiaoqing, W. Allison, J. Ellis, and G. Alexandrowicz,Probing molecule-surface interactions through ultra-fast adsorbate dynamics: propane/Pt(111),New J. Phys. 10, 125026 (2008).[145] F. E. Tuddenham, H. Hedgeland, J. Knowling, A. P. Jardine, D. A. MacLaren, G. Alexandrowicz,J. Ellis, and W. Allison, Linewidths in bound state resonances for helium scattering from Si(111)-(1x 1)H, J. Phys. Condens. Matter 21, 264004 (2009).[146] A. P. Jardine, E. Y. M. Lee, D. J. Ward, G. Alexandrowicz, H. Hedgeland, W. Allison, J. Ellis, andE. Pollak, Determination of the quantum contribution to the activated motion of hydrogen on a metalsurface: H/Pt(111), Phys. Rev. Lett. 105, 136101 (2010). → page 43[147] A. Tamtoegl, D. Campi, M. Bremholm, E. M. J. Hedegaard, B. B. Iversen, M. Bianchi, P. Hofmann,N. Marzari, G. Benedek, J. Ellis, and W. Allison, Nanoscale surface dynamics of Bi2Te3(111):observation of a prominent surface acoustic wave and the role of van der Waals interactions,Nanoscale 10, 14627 (2018). → page 43119[148] M. Sacchi, P. Singh, D. M. Chisnall, D. J. Ward, A. P. Jardine, W. Allison, J. Ellis, andH. Hedgeland, The dynamics of benzene on Cu(111): a combined helium spin echo anddispersion-corrected DFT study into the diffusion of physisorbed aromatics on metal surfaces,Faraday Discuss. 204, 471 (2017). → page 43[149] N. R. Hutzler, H.-I. Lu, and J. M. Doyle, The buffer gas beam: An intense, cold, and slow source foratoms and molecules, Chem. Rev. 112, 4803 (2012). → page 45[150] H. Bethlem, A. van Roij, R. Jongma, and G. Meijer, Alternate gradient focusing and deceleration ofa molecular beam, Phys. Rev. Lett. 88, 133003 (2002). → page 45[151] A. W. Wiederkehr, H. Schmutz, M. Motsch, and F. Merkt, Velocity-tunable slow beams of cold O2 ina single spin-rovibronic state with full angular-momentum orientation by multistage Zeemandeceleration, Mol. Phys. 110, 1807 (2012). → page 45[152] R. V. Krems, Molecules in Electromagnetic Fields: From Ultracold Physics to Controlled Chemistry(Wiley, New York 2019). → page 45[153] A. P. Jardine, P. Fouquet, J. Ellis, and W. Allison, Hexapole magnet system for thermal energy 3Heatom manipulation, Rev. Sci. Instrum. 72, 3834 (2001). → pages 45, 50[154] S. Dworski, G. Alexandrowicz, P. Fouquet, A. P. Jardine, W. Allison, and J. Ellis, Low aberrationpermanent hexapole magnet for atom and molecular beam research, Rev. Sci. Instrum. 75, 1963(2004). → pages 45, 50[155] I. Litvin, Y. Alkoby, O. Godsi, G. Alexandrowicz, and T. Maniv, Parallel and anti-parallel echoes inbeam spin echo experiments, Results in Physics 12, 381 (2019). → pages 46, 67[156] N. F. Ramsey, Theory of molecular hydrogen and deuterium in magnetic fields, Phys. Rev. 85, 60(1952). → pages 48, 60, 65, 66[157] M. Utz, M. H. Levitt, N. Cooper, and H. Ulbricht, Visualisation of quantum evolution in theStern–Gerlach and Rabi experiments, Phys. Chem. Chem. Phys. 17, 3867 (2015). → page 51[158] C. Kru¨ger, E. Lisitsin-Baranovsky, O. Ofer, P.-A. Turgeon, J. Vermette, P. Ayotte, andG. Alexandrowicz, A magnetically focused molecular beam source for deposition of spin-polarisedmolecular surface layers, J. Chem. Phys. 149, 164201 (2018). → page 51[159] J. S. Walker and J. Gathright, A transfer-matrix approach to one-dimensional quantum mechanicsusing Mathematica, Computers in Physics 6, 393 (1992). → pages 55, 58, 61[160] R. N. Zare, Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics, TheGeorge Fisher Baker non-resident lectureship in chemistry at Cornell University (Wiley, New York1988). → page 62[161] C. Diaz, E. Pijper, R. A. Olsen, H. F. Busnengo, D. J. Auerbach, and G. J. Kroes, Chemicallyaccurate simulation of a prototypical surface reaction: H2 dissociation on Cu(111), Science 326,832 (2009). → page 63[162] J. M. B. Kellogg, I. I. Rabi, N. F. Ramsey, and J. R. Zacharias, The magnetic moments of the protonand the deuteron. The radiofrequency spectrum of H2 in various magnetic fields, Phys. Rev. 56, 728(1939). → page 66120[163] J. M. B. Kellogg, I. I. Rabi, N. F. Ramsey, and J. R. Zacharias, An electrical quadrupole moment ofthe deuteron the radiofrequency spectra of HD and D2 molecules in a magnetic field, Phys. Rev. 57,677 (1940). → page 66[164] F. Mezzadri, How to generate random matrices from the classical compact groups, Notices Am.Math. Soc. 54, 13 (2007). → pages 70, 71[165] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning, AdaptiveComputation and Machine Learning (MIT Press, Cambridge, Mass 2006). → pages80, 82, 83, 84, 85, 89[166] Y. L. Tong, The Multivariate Normal Distribution, Springer series in statistics (Springer, New York1990). → page 80[167] D. J. Griffiths, Introduction to Quantum Mechanics (Pearson Prentice Hall, Upper Saddle River, NJ2005), 2nd ed edition. → page 81[168] D. K. Duvenaud, Automatic model construction with Gaussian processes, Ph.D. Thesis, Universityof Cambridge (2014). → page 84[169] T. W. Simpson, J. D. Poplinski, P. N. Koch, and J. K. Allen, Metamodels for computer-basedengineering design: Survey and recommendations, Engineering with Computers 17, 129 (2001). →page 85[170] V. Dubourg, B. Sudret, and F. Deheeger, Metamodel-based importance sampling for structuralreliability analysis, Probabilistic Engineering Mechanics 33, 47 (2013). → page 85[171] J. L. Loeppky, J. Sacks, and W. J. Welch, Choosing the sample size of a computer experiment: Apractical guide, Technometrics 51, 366 (2009). → page 85[172] E. Brochu, V. M. Cora, and N. de Freitas, A tutorial on Bayesian optimization of expensive costfunctions, with application to active user modeling and hierarchical reinforcement learning (2010),arXiv:1012.2599. → page 88[173] F. Berkenkamp, A. P. Schoellig, and A. Krause, No-regret Bayesian optimization with unknownhyperparameters (2019), arXiv: 1901.03357. → pages 88, 92[174] M. D. McKay, R. J. Beckman, and W. J. Conover, A comparison of three methods for selectingvalues of input variables in the analysis of output from a computer code, Technometrics 21, 239(1979). → pages 90, 101[175] R. A. Vargas-Herna´ndez, (private communication, 2019). → pages 94, 109[176] T. Joyce, (private communication, 2019). → pages 94, 109[177] C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York 2006). → page 94[178] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of statecalculations by fast computing machines, J. Chem. Phys. 21, 1087 (1953). → page 96[179] T. Zeng, G. Guillon, J. T. Cantin, and P.-N. Roy, Probing the superfluid response of para-hydrogenwith a sulfur dioxide dopant, J. Phys. Chem. Lett. 4, 2391 (2013). → page 96121[180] T. Zeng and P.-N. Roy, Microscopic molecular superfluid response: theory and simulations, Rep.Prog. Phys. 77, 046601 (2014). → page 96[181] P. Huang and E. A. Carter, Advances in correlated electronic structure methods for solids, surfaces,and nanostructures, Annu. Rev. Phys. Chem. 59, 261 (2008). → page 96[182] J. Kolorenc and L. Mitas, Applications of quantum Monte Carlo methods in condensed systems,Rep. Prog. Phys. 74, 026502 (2011).[183] B. M. Austin, D. Y. Zubarev, and W. A. Lester, Quantum Monte Carlo and related approaches,Chem. Rev. 112, 263 (2012). → page 96[184] H. G. Evertz, The loop algorithm, Adv. Phys. 52, 1 (2003). → pages 96, 97[185] N. V. Prokof’ev and B. V. Svistunov, Bold diagrammatic Monte Carlo technique: When the signproblem is welcome, Phys. Rev. Lett. 99, 250201 (2007). → pages 96, 97, 101, 104[186] N. V. Prokof’ev and B. V. Svistunov, Bold diagrammatic Monte Carlo: A generic sign-problemtolerant technique for polaron models and possibly interacting many-body problems, Phys. Rev. B77, 125101 (2008). → pages 96, 97, 101, 104[187] R. H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys.Rev. Lett. 58, 86 (1987). → page 97[188] M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation (Oxford UniversityPress, Oxford; New York 2010). → page 97[189] H. Gould and J. Tobochnik, Overcoming critical slowing down, Computers in Physics 3, 82 (1989).→ page 97[190] U. Wolff, Critical slowing down, Nucl. Phys. B Proc. Suppl. 17, 93 (1990). → page 97[191] N. V. Prokof’ev and B. V. Svistunov, Polaron problem by diagrammatic quantum Monte Carlo,Phys. Rev. Lett. 81, 2514 (1998). → page 97[192] G. Bighin, (private communication, 2019). → page 97[193] M. C. Kennedy and A. O’Hagan, Bayesian calibration of computer models, J. Royal Stat. Soc. B 63,425 (2001). → pages 98, 99, 100[194] K. Van Houcke, F. Werner, E. Kozik, N. V. Prokof’ev, B. V. Svistunov, M. J. H. Ku, A. T. Sommer,L. W. Cheuk, A. Schirotzek, and M. W. Zwierlein, Feynman diagrams versus Fermi-gas Feynmanemulator, Nat. Phys. 8, 366 (2012). → page 101[195] A. S. Mishchenko, N. Nagaosa, and N. V. Prokof’ev, Diagrammatic Monte Carlo method formany-polaron problems, Phys. Rev. Lett. 113, 166402 (2014). → page 101[196] I. S. Tupitsyn and N. V. Prokof’ev, Stability of Dirac liquids with strong coulomb interaction, Phys.Rev. Lett. 118, 026403 (2017). → page 101[197] A. L. Fetter and J. D. Walecka, Quantum theory of many-particle systems (Dover Publications,Mineola, N.Y 2003). → page 101122[198] L. D. Landau and E. M. Lifshitz, Quantum Mechanics: Non-Relativistic Theory, number 3 in Courseof Theoretical Physics (Elsevier, Singapore 2007), 3rd ed., rev. and enl., authorized Engl. reprintedition. → page 101[199] G. Bighin, T. Tscherbul, and M. Lemeshko, Diagrammatic Monte Carlo approach to angularmomentum in quantum many-particle systems, Phys. Rev. Lett. 121, 165301 (2018). → page 104123Appendix ASupplementary material for Andersonlocalization for particles with long-rangehoppingA.1 Outline of the Derivation for Equation 2.26Here, I give an outline of the derivation of Eqn. 2.26, the eigenvalues of the Hamiltonian (2.24) for an ideal3D cubic lattice with isotropic hopping of arbitrary range (2.22). We begin with the Hamiltonian (2.24) andevaluateHˆIα |k1k2k3〉=1√VN−12∑i, j,h,l,m,n,α,β ,δ=−N−12ei2piN (k1α+k2β+k3δ )ζrα(1−δilδ jmδhn)cˆ†lmncˆi jhcˆ†αβδ |0〉 , (A.1)where|k1k2k3〉= cˆ†k1k2k3 |0〉 (A.2)=1√VN−12∑α,β ,δ=−N−12ei2piN (k1α+k2β+k3δ )cˆ†αβδ |0〉 , (A.3)124V ≡ N3, r ≡√(l− i)2+(m− j)2+(n−h)2, and |0〉 is the vacuum state. Using the commutation relations[cˆlmn, cˆ†i jσ]≡ δliδm jδσn (A.4)[cˆlmn, cˆi jσ]=[cˆ†lmn, cˆ†i jσ]≡ 0, (A.5)and evaluating the sums over α,β ,δ in Eqn. A.1, we getHˆIα |k1k2k3〉=ζ√VN−12∑i, j,h,l,m,n=−N−12ei2piN (k1i+k2 j+k3h)1rα(1−δilδ jmδhn)cˆ†lmn |0〉 (A.6)Then, to evaluate the eigenvalues (2.26), we multiply from the left by 〈k1k2k3| and obtain〈k1k2k3|HˆIα |k1k2k3〉=ζVN−12∑i, j,h,l,m,n,α,β ,δ=−N−121rα(1−δilδ jmδhn)ei2piN (k1(i−α)+k2( j−β )+k3(h−δ )) 〈0|cαβδ cˆ†lmn|0〉(A.7)=ζVN−12∑i, j,h,l,m,n=−N−121rα(1−δilδ jmδhn)ei2piN (k1(i−l)+k2( j−m)+k3(h−n)) (A.8).Note that one can perform this reordering of the sum:N−12∑i,l=−N−12h(i− l) =N−1∑r1=−(N−1)(N−|r1|)h(r1), (A.9)where r1 = i− l, h is some function of i− l, and (N−|r1|) counts the number of times a particular value ofr1 appears in the double sum on the left-hand side.125We perform a similar reordering for Eqn. A.8 and obtain〈k1k2k3|HˆIα |k1k2k3〉=ζVN−1∑r1=−(N−1)(N−|r1|)N−12∑j,h,m,n=−N−121rα(1−δr10δ jmδhn)ei2piN (k1r1+k2( j−m)+k3(h−n))(A.10)=ζVN−12∑j,h,m,n=−N−12ei2piN (k2( j−m)+k3(h−n))N−1∑r1=−(N−1)(N−|r1|) 1rα(1−δr10δ jmδhn)ei2piN (k1r1)(A.11)=ζVN−12∑j,h,m,n=−N−12ei2piN (k2( j−m)+k3(h−n)) f ( j−m,h−n) , (A.12)where r1 = i− l, now r =√r21 +(m− j)2+(n−h)2, andf ( j−m,h−n)≡N−1∑r1=−(N−1)(N−|r1|) 1rα(1−δr10δ jmδhn)ei2piN (k1r1) (A.13)= N(1−δ( j−m)0δ(h−n)0)((m− j)2+(n−h)2) α2+2N−1∑r1=1(N− r1)cos(2piN k1r1)rα, (A.14)where the first term of the second line corresponds to r1 = 0 and the second term arises as δr10 = 0 (r1 6= 0in the second term).Repeating the process for r2 ≡ j−m, we get〈k1k2k3|HˆIα |k1k2k3〉=ζVN−12∑h,n=−N−12ei2piN k3(h−n)g(h−n) , (A.15)whereg(h−n)≡N−1∑r2=−(N−1)(N−|r2|) f (r2,h−n)ei 2piN (k2r2) (A.16)= N f (0,h−n)+2N−1∑r2=1(N− r2) f (r2,h−n)cos(2piNk2r2). (A.17)126Finally, repeating the process again for r3 ≡ h−n, we get〈k1k2k3|HˆIα |k1k2k3〉=ζV[Ng(0)+2N−1∑r3=1(N− r3)g(r3)cos(2piNk3r3)]. (A.18)We then haveg(0) = N f (0,0)+2N−1∑r2=1(N− r2) f (r2,0)cos(2piNk2r2)(A.19)= 2NN−1∑r1=1(N− r1)cos(2piN k1r1)rα1+2NN−1∑r2=1(N− r2)cos(2piN k2r2)rα2+4N−1∑r1=1N−1∑r2=1(N− r1)(N− r2)cos(2piN k1r1)+ cos(2piN k2r2)(r21 + r22) α2(A.20)asf (r2,0) = N(1−δr20)rα2+2N−1∑r1=1(N− r1)cos(2piN k1r1)(r21 + r22) α2and (A.21)f (0,0) = 2N−1∑r1=1(N− r1)cos(2piN k1r1)rα1. (A.22)To obtain the right hand side of Eqn. A.22, note that (1−δ( j−m)0δ(h−n)0)((m− j)2+(n−h)2) α2= 0 when j = m and h = n. This isbecause, in the Hamiltonian (2.24), the term (1−δilδ jmδhn)rα is understood to be zero when i = l, j = m, andh = n (i.e. the particle can hop anywhere except to its original site).Then,g(r3) = N f (0,r3)+2N−1∑r2=1(N− r2) f (r2,r3)cos(2piNk2r2)(A.23)=N2rα3+2NN−1∑r1=1(N− r1)cos(2piN k1r1)(r21 + r23) α2+2NN−1∑r2=1(N− r2)cos(2piN k2r2)(r22 + r23) α2+4N−1∑r1=1N−1∑r2=1(N− r1)(N− r2)cos(2piN k1r1)cos(2piN k2r2)rα, (A.24)127asf (r2,r3) =N(r22 + r23) α2+2N−1∑r1=1(N− r1)cos(2piN k1r1)rαand (A.25)f (0,r3) =Nrα3++2N−1∑r1=1(N− r1)cos(2piN k1r1)(r21 + r23) α2, (A.26)where the fact that r3 > 0 was taken into account (see last term of Eqn. A.18). We then substitute Eqns. A.24and A.20 into Eqn. A.18. Finally, by expanding and collecting terms, we obtain Eqn. 2.26 after relabellingr1,r2, and r3 to l,m, and n, respectively.128Appendix BSupplementary material for molecularhyperfine interferometryB.1 Schro¨dinger Equation for Eigenstate CoefficientsUsing Eqn. 3.19, the time-independent Schro¨dinger equation isHˆ |ER˜〉= E |ER˜〉 (B.1)Kˆ |ER˜〉=(E− HˆR)|ER˜〉∑R∫dx ΦER˜R (x)Kˆ |xR〉=∑R∫dx ΦER˜R (x)(E− HˆR)|xR〉∑R∫dx ΦER˜R (x)〈x0R0| Kˆ |xR〉=∑R∫dx ΦER˜R (x)〈x0R0|(E− HˆR)|xR〉 (B.2)where Hˆ is the total Hamiltonian (3.2) in the current region, Kˆ ≡ pˆ22m , we use HˆR as a shorthand for HˆR(~Bloc),and the last line was multiplied by 〈x0R0|.129The different terms can be evaluated as〈x0R0| Kˆ |xR〉= δR0R 〈x0|∫dkh¯2k22m|k〉〈k|x〉=∫ dk2piδR0Rh¯2k22meik(x0−x), (B.3)〈x0R0|E |xR〉= δR0Rδ (x− x0)E, (B.4)〈x0R0| HˆR |xR〉= δ (x− x0)HRR0R, (B.5)where |k〉 is a momentum state with wavenumber k, and m is the mass of the molecule. The additionalfactor of (2pi)−1 in Eqn. (B.3) comes from 〈x|k〉 ≡ (2pi)− 12 eikx. After inserting these three equations intoEqn. (B.2) and evaluating most of the sums, we obtain∫dx∫ dk2pih¯2k22meik(x0−x)ΦER˜R0 (x) =ΦER˜R0 (x0)E−∑RHRR0RΦER˜R (x0) (B.6)Noting that k2eikx0 = − ∂ 2∂x20 eikx0 and∫ dk2pi eik(x0−x) = δ (x0− x), we obtain Eqn. (3.20) after the relabellingx0→ x.B.2 Coefficient Relations Across a DiscontinuitySince ΦER˜R (x) ∈C1(x) for a specific value of R and given Eqn. (3.26), we get the defining equations for thecontinuity of the wavefunction aslimx→0−ΦER˜R+ (x) = limx→0+ΦER˜R+ (x)limx→0−∑R−ΦER˜R− (x)S∗R−R+ = limx→0+ΦER˜R+ (x)limx→0−∑R−S∗R−R+(AR−eikR−x+BR−e−ikR−x)= limx→0+AR+eikR+x+BR+e−ikR+x (Eqn. (3.22))AR+ +BR+ =∑R−S∗R−R+ (AR−+BR−) , (B.7)where S∗R−R+ ≡ 〈R+|R−〉, kR± ≡√2m(E−ER±)h¯ , and ER± ≡ 〈R±|HˆR(~B(0±))|R±〉. There are NR such equations,one for each value of R+.130Correspondingly, the defining equations for the continuity of the first derivative of the coefficients arelimx→0−∂∂xΦER˜R+ (x) = limx→0+∂∂xΦER˜R+ (x)limx→0−∑R−∂∂xΦER˜R− (x)S∗R−R+ = limx→0+∂∂xΦER˜R+ (x)limx→0−∑R−S∗R−R+∂∂x(AR−eikR−x+BR−e−ikR−x)= limx→0+∂∂x(AR+eikR+x+BR+e−ikR+x)(Eqn. (3.22))limx→0−∑R−S∗R−R+ ikR−(AR−eikR−x−BR−e−ikR−x)= limx→0+ikR+(AR+eikR+x−BR+e−ikR+x)AR+−BR+ =∑R−S∗R−R+kR−kR+(AR−−BR−) , (B.8)Solving Eqns. (B.7) and (B.8) for the coefficients AR+ and BR+ , we obtain Eqns. (3.27) and (3.28).B.3 Computational Parameters Used for the Application toortho-HydrogenWe take the mean velocity v0 = 1436.14ms−1 and the velocity spread to be 4 %FWHM. When performingthe integral of Eqn. (3.18), we take a k-space grid spacing ∆k = 1×104 cm−1 and integrate from −7σk to+7σk, where σk is the Gaussian width in momentum space as defined in Section 3.4. For the magnetic fieldprofile and the angles between the two branches of the apparatus, see Figure 3.3. The relative probabilitiesused for the state selector probabilities PR0 and the detector coefficients cRD are given in Table B.1.Where applicable, the parameters above were chosen to match those in the supplementary informationof Godsi et al. [44].Table B.1: Relative probabilities of the state selector ηmImJ and the detector κmImJ . The state selectorprobabilities PR0 are calculated as PR0 = PmImJ ≡ ηmImJ/∑mImJ ηmImJ and the detector coefficientscRD are calculated as cRD = cmImJ ≡ κmImJ/∑mImJ κmImJ .mI 1 1 1 0 0 0 -1 -1 -1mJ 1 0 -1 1 0 -1 1 0 -1ηmImJ 1.0000 0.9755 0.7901 0.1465 0.1111 0.0738 0.0343 0.0299 0.0258κmImJ 1.00 0.96 0.93 0.53 0.42 0.37 0.21 0.19 0.16131
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Interference in complex quantum systems : from localization...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Interference in complex quantum systems : from localization in high-dimensional lattices to surface spin… Cantin, Joshua Tyler 2019
pdf
Page Metadata
Item Metadata
Title | Interference in complex quantum systems : from localization in high-dimensional lattices to surface spin echo with molecules |
Creator |
Cantin, Joshua Tyler |
Publisher | University of British Columbia |
Date Issued | 2019 |
Description | This thesis considers (1) novel manifestations and applications of quantum interference in complex systems and (2) development of new approaches to study complex quantum systems. First, I examine a general model of particles with long-range hopping amplitudes. For at least 30 years, it has been widely accepted that these particles do not undergo Anderson localization in 3D lattices. We show that these particles do undergo Anderson localization in 3D lattices if their hopping amplitudes are isotropic. In contrast, particles with anisotropic long-range hopping appear to follow the widely held expectations. We show these results by demonstrating that cooperative shielding extends to 3D cubic lattices with isotropic long-range hopping, but not with anisotropic long-range hopping and by computing the scaling behaviour of participation ratios and energy level statistics. Secondly, I develop a fully quantum mechanical model of molecular surface spin-echo experiments, which study surface properties and dynamics by scattering molecules off the sample surface. This model, based on the transfer matrix method, incorporates molecular hyperfine degrees of freedom, allows for the efficient calculation of the experimental signal given a molecule-surface scattering matrix, and permits us to begin addressing the inverse scattering problem. This fully quantum model is required to properly describe these experiments as the semi-classical methods used to describe experiments using helium-3 atoms do not take magnetic-field induced momentum changes into account. We apply our method to the case of ortho-hydrogen and then apply Bayesian optimization to determine the molecule-surface scattering matrix elements from a calculated signal, for a scattering matrix defined by three parameters. Our work sets the stage for Bayesian optimization to solve the inverse scattering problem for these experiments. Finally, I propose using Bayesian model calibration to improve the convergence of Monte Carlo calculations in regions where the sign problem or critical slowing down are an issue. Specifically, Bayesian model calibration would correct poorly converged Monte Carlo calculations with the information from a small number of well-converged Monte Carlo calculations. As a simple proof of principle demonstration, we apply Bayesian model calibration to a diagrammatic Monte Carlo calculation of the scattering length of a spherical potential barrier. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2019-09-04 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0380808 |
URI | http://hdl.handle.net/2429/71600 |
Degree |
Doctor of Philosophy - PhD |
Program |
Chemistry |
Affiliation |
Science, Faculty of Chemistry, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2019-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2019_november_cantin_joshua.pdf [ 3.86MB ]
- Metadata
- JSON: 24-1.0380808.json
- JSON-LD: 24-1.0380808-ld.json
- RDF/XML (Pretty): 24-1.0380808-rdf.xml
- RDF/JSON: 24-1.0380808-rdf.json
- Turtle: 24-1.0380808-turtle.txt
- N-Triples: 24-1.0380808-rdf-ntriples.txt
- Original Record: 24-1.0380808-source.json
- Full Text
- 24-1.0380808-fulltext.txt
- Citation
- 24-1.0380808.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0380808/manifest