ABSTRACT Title of dissertation: QUANTUM MANY-BODY PHENOMENA IN ULTRA COLD ATOMS IN OPTICAL LATTICES Anzi Hu, Doctor of Philosophy, 2011 Dissertation directed by: Professor Bei-Lok Hu Department of Physics Two models are discussed here to illustrate the quantum many-body phenomena in mixtures of ultra-cold atoms in optical lattices. The first model describes a mixture of two species of bosonic atoms of equal masses in optical lattices and the second describes a mixture of heavy bosonic atoms and light fermionic atoms in optical lattices. For both models, we assume the trap is present and use parameters typical in experiment. For the first model, the discussion is aimed at providing a thorough description of the collective behavior of the binary mixture in various interaction regions, with emphasis on two many-body phenomena, pairing and anti-pairing, as a result of the inter-species inter- action. The pairing leads to a new type of superfluid order, called the paired superfluid (PSF) and the anti-pairing leads to another type of superfluid order, called the counter- flow superfluid (CFSF). In addition, we discuss the coexistence of charge density wave order with the three superfluid orders in the strong interaction region. We use both Lut- tinger liquid theory and the time evolving block decimation (TEBD) method to study this model in one dimension. The discussion is organized in three parts: the phase diagram and the correlation functions; the noise correlation functions; and the transport proper- ties. Two phase diagrams are constructed to map the different orders in the parameter space. The correlation functions, include noise correlations, are carefully examined for the determination of the orders and for possible detection methods. In the end, the transport properties of the PSF and CFSF orders are studied through the dipole oscillation induced by trap displacement. For the second model, examining a mixture of heavy bosons and light fermions, the discussion is oriented toward determining the thermal properties of the mixture for attrac- tive inter-species interactions. This work is motivated by experiments creating artificial molecules through optical and magnetic control of ultra-cold atoms. We use the strong coupling (SC) expansion method to evaluate the density profile, the onsite inter-species correlations, the density fluctuations and the entropy per particle. Analytical expressions are derived for all the quantities above as well as the partition function. To benchmark the accuracy, the SC calculations are compared with inhomogeneous dynamical mean field theory (IDMFT) and Monte Carlo (MC) simulation. From the calculations, we find that 1) the efficiency of creating pre-formed molecules is significantly increased by confining the mixtures onto optical lattices; 2) the temperature of the mixtures in optical lattices can be reliably estimated through the density gradient and the density fluctuations. QUANTUM MANY-BODY PHENOMENA IN ULTRA-COLD ATOMS IN OPTICAL LATTICES by Anzi Hu Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2011 Advisory Committee Professor Bei-Lok Hu, Chair Professor Carl Williams Professor James Freericks (Georgetown University) Professor Eite Tiesinga Professor Victor Galitski Professor Ramani Duraiswami, Dean?s representative c Copyright by Anzi Hu 2011 to my parents ii Acknowledgments Five years of living and studying in Maryland has been quite a journey. I am grateful to all the people who helped and befriended me along the way and will always cherish the memories of my time in graduate school. I am fortunate enough to meet Dr. Bei-Lok Hu in the very beginning of the graduate school. I am thankful for the many interesting discussions we had, his introduction to Dr. Carl Williams and other scientists at National Institute of Stan- dards and Technology (NIST), his support for me to pursue research in the field of ultra-cold atoms and his kind reminder that as a theorist, it is always important to keep an eye on the bigger picture and think deeper. Besides research, I am grateful for his sharing many classical music resources in the area, which has become an essential part of my life in the past few years. I am equally fortunate to be under the guidance of Dr. Carl Williams and am grateful to him for giving me the opportunity to work at NIST and interact with so many brilliant scientists. I would like to thank Dr. Williams for his introduc- tion to computational physics and specifically the time evolving block decimation (TEBD) method, which has lead to fruitful research and helped my understanding of many-body effects in reduced dimensions, and his introduction to Dr. Maciej Ma?ska of University of Silesia and Dr. Jim Freericks of Georgetown University. I am also very grateful for being able to walk into Dr. Williams? office whenever I have a question, despite his busy schedule. I cannot imagine being in today?s position without his guidance and advices. During my stay at NIST, I have met and interacted with many brilliant scien- tists. Among them, I would like to first thank Dr. Ludwig Mathey for his guidance iii and collaboration on studies of one dimensional Bose mixtures. I am grateful for his introduction to the rich many-body phenomena in one dimension and to the ultra-cold atom experiments, for answering so many of my naive questions and for helping me through the scientific paper writing process. It has been a pleasure working with him. I would like to thank Dr. Ippei Danshita for his assistance on the TEBD code and very helpful discussions. I would like to thank Dr. Eite Tiesinga for his advice and discussions, particularly for his involvement in the research on one dimensional Bose mixtures and his careful editing of our paper. I also would like thank Dr. Charles Clark for giving me the opportunity of working with Dr. Mathey and Dr. Danshita in his division, spending time editing our paper on the noise correlations and his kind invitation to the QIBEC talks as well as many other events. I would like thank Dr. Paul Julienne for the helpful discussions on atomic systems in optical lattices. While working on the Bose-Fermi mixtures, I was fortunate to meet Dr. Maciej Ma?ska and Dr. Jim Freericks. I am especially grateful for Dr. Ma?ska?s hospitality during my stay in Poland and for his teachings on the Monte Carlo method and on numerical simulation in general. I would like to thank Dr. Jim Freericks for his in- troduction to the strong coupling expansion method. I have benefited greatly from his thorough and solid knowledge in condensed matter physics and his careful attitude towards research. Besides research, I have enjoyed many graduate courses. I still feel quite fortu- nate for learning statistical mechanics from Dr. Michael Fisher. His teaching was iv not just about the knowledge but also about scientific thinking and questioning and his homework problems often felt like research projects, which can be chal- lenging sometimes, but always stimulating. I also want to thank Dr. Victor Galitski for his class on solid state physics, where I got my first glimpse of renormalization group theory, and many other important concepts and methods in condensed mat- ter physics. Although I have not worked on any research related with relativity, Dr. Ted Jacobson?s introduction to general relativity was so fresh and enjoyable that I was completely absorbed in the curved space during that course. I am also very happy that I chose to take the scientific computing class by Dr. Ramani Du- raiswami at the end of my coursework Dr. Duraiswami?s class introduced me to the basic computer science knowledge that turned out to be very helpful to my research and has since greatly influenced the way I write programming codes. I owe my deepest thanks to my family - especially my mother, who has been staying with me for the past half year. She is always my best friend and best sup- porter. She and my dad have always respected and supported me even when they disagree with my choices. Words cannot express my gratitude for their uncondi- tional love. In retrospective, so many people inside and outside school have helped me along my five years in Maryland. Without them, this thesis would never be possi- ble. I can?t believe myself how lucky I have been. v Contents List of Figures ix List of Tables xxiv 1 Introduction 1 1.1 Second quantization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.1 Bose Hubbard model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.2 Multi-component Hamiltonian . . . . . . . . . . . . . . . . . . . . . . 7 1.2 Phase diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2.1 Off diagonal long-range order . . . . . . . . . . . . . . . . . . . . . . 9 1.2.2 Order parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.3 Superfluid to Mott insulator phase transition . . . . . . . . . . . . . 12 1.3 Tomonaga-Luttinger liquid theory . . . . . . . . . . . . . . . . . . . . . . . . 13 1.4 Detection Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.1 Time-of-flight measurement . . . . . . . . . . . . . . . . . . . . . . . . 17 1.4.2 Noise correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.4.3 In situ measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2 Time evolving block decimation method 21 2.1 Basic TEBD method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1.1 State decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.1.2 Single-site and two-site operations . . . . . . . . . . . . . . . . . . . . 25 2.1.3 Time evolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 2.1.4 Observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2 Extended TEBD method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.1 Number conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.2.2 One dimensional two-component Hamiltonian . . . . . . . . . . . . 35 2.3 Efficiency and accuracy estimation . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4 Example: one-dimensional Bose Hubbard model . . . . . . . . . . . . . . . . 37 vi 2.4.1 Homogeneous System . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4.2 Harmonic trap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4.3 Accuracy estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3 One dimensional binary Bose mixtures in optical lattices (I): paired and counter- flow superfluidity 46 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2 Two-component Bose Hubbard model . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Luttinger liquid approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.4 Numerical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4.1 Phase diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.4.1.1 Phase diagram at a fixed hopping parameter . . . . . . . . 62 3.4.1.2 Phase diagram at half-filling . . . . . . . . . . . . . . . . . . 63 3.4.2 Realization and detection . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4 One dimensional binary Bose mixtures in optical lattices (II): noise correlation 73 4.1 Noise correlations and quasi-orders . . . . . . . . . . . . . . . . . . . . . . . 74 4.2 Luttinger liquid approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.3 Numerical approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.3.1 Homogeneous systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.3.2 Trapped systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.3.3 Determination of Luttinger parameters from experimental data . . . 91 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5 One dimensional binary Bose mixtures in optical lattices (III): transport proper- ties 93 5.1 Dynamics after a trap displacement . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.1 Model and parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.2 Simulation of the real time evolution . . . . . . . . . . . . . . . . . . . 95 vii 5.2 Result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 5.2.1 Brief displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2.2 Constant displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6 Two dimensional Bose-Fermi mixtures in optical lattices: pre-formed molecules and novel thermometry 105 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.2 Fermi-Bose Falicov-Kimball model . . . . . . . . . . . . . . . . . . . . . . . . 108 6.3 Strong coupling expansion formalism . . . . . . . . . . . . . . . . . . . . . . 112 6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 6.4.1 Comparison with the IDMFT and MC calculations . . . . . . . . . . . 119 6.4.2 Finite-size effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 6.5 Thermometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 6.5.1 Temperature and Density fluctuations . . . . . . . . . . . . . . . . . . 128 6.5.2 Fluctuation calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 6.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 7 Conclusion and further discussions 136 Bibliography 139 viii List of Figures 2.1 Number fluctuations j and density distributionhnjifor U=t = 1(a), 4(b), 10(c), 20(d). The total number of particle is 50. The density distribution is uniformly 1 except for at the boundary. In (a) and (b), the system is in the SF state, this is reflected by the large number fluctuation. In (c) and (d), the system is in the MI state and the number fluctuations are much reduced. . . 39 2.2 One-body density matrix Gij = hbyibjifor U=t = 1(a), 4(b), 10(c) and 20(d). The parameters are the same as in Fig. 2.1. In (a) and (b) the system is in the SF state. In (c) and (d) the system is in the MI state. In (a), Gij decays very slowly and the correlation function is forced to zero by the boundary condition at the edge of the lattice. In (b), Gij decays faster than in (a) but still qualifies as algebraic decay (see also Fig. 2.3) with a larger decay power. In (c) and (d), Gijbecomes short-ranged and the decay behavior changes to exponential. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.3 Algebraic and exponential fit for Gij. Here we choose i = 25, j = 25;:::50. The functionG25;j represent the decay behavior of the off-diagonal one body density matrix. Through the fit, one can obtain the Luttinger parameter K, K = 9:4 in (a) and K = 6:4 in (b). The decay length X0 can be found for the exponential decay, X0 = 1:5 in (c) and X0 = 0:75 in (d). . . . . . . . . . . . . 41 2.4 Density distribution and number fluctuations in a trapped system. In (a), the system is in the SF state with U=t = 1 In (b), the interaction is strong, U=t = 20. The density forms a plateau near the center of the trap. This is where the MI state is formed. The fluctuations are reduced at the plateau area as the result of the MI order. The SF state still exists at the edge of the trap. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 ix 2.5 One-body density matrix in a trapped system. The parameters in (a) and (b) are the same as in Fig. 2.4. In (a), the system is in the SF state. The density matrix is reduced dramatically to zero near the edge of the cloud. In (b), the off-diagonal decay is exponential in the MI region. The decay changes to algebraic as a result of the SF state at the edge. . . . . . . . . . . . . . . . . . 43 2.6 Vector [l] as a function of l and . The decomposition rank is 50 and the size of the system is 50 lattice sites. In (a) and (b), the value of decreases to around 10 4 for = 50 at the center. In (c) and (d), the value of decays rapidly as increases and is less than 10 9 at = 50. . . . . . . . . . . . . . 44 2.7 Decay of [25] in SF and MI states. Here, [25]is chosen because it relates with the Schmidt rank at the center of the lattice, which usually has the slowest decay with regard to . The vector [25] decays much faster in the MI region than in the SF region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.8 Behavior of [j] for a trapped system. In (a), [l] corresponds to the SF state of Fig. 2.4 (a) and in (b), [l] corresponds to the MI-SF mixed state of Fig. 2.4 (b). The behavior of reflects the orders that exists in the trap. . . . . . . . 45 3.1 Sketch of a condensate of pairs. Atoms of each species (red/green) pair together and form a paired superfluid (PSF) state. . . . . . . . . . . . . . . . 47 3.2 Sketch of a condensate of anti-pairs. Here, atoms of one species are strongly anti-correlated with atoms of the other species, creating a counterflow su- perfluid (CFSF) state. These composite bosons can also be thought of as a pair of one atom of one species and one hole of the other species. . . . . . . 48 3.3 Phase diagram of a bosonic mixture at non-unit and non-half-filling. For at- tractive interactions U12 and K < 2 the system can form a paired superfluid state, in the regime labeled PSF and PSF(CDW). This phase can coexist with CDW order for weaker interactions. For large repulsive (attractive) interac- tionsU12 the system phase separates (PS) (collapses (CL)). For the remaining regime the system shows single particle superfluidity (SF). This can coexist with CDW order, resulting in a quasi-supersolid (SS) regime. . . . . . . . . . 53 x 3.4 Phase diagram of a bosonic mixture at half-filling. In addition to the phases that appear in Fig. 3.3, the system now develops a counterflow superfluid (CFSF) phase, which can coexist with CDW order. . . . . . . . . . . . . . . . 55 3.5 Correlation functions RA, RS, and G on a logarithmic scale as a function of distanceji jj. The indexiis 40, the center of the 80 lattice sites. The squares are the numerical data. The blue lines are exponential fits to the data and red dotted lines are algebraic fits. Note that the scale of the vertical axis of the graphs differs by orders of magnitude. In (a), we show an example for the paired superfluid phase at =0.3, t = 0:02U, and U12 = 0:16U. RA decays exponentially andRS decays algebraically. The single-particle corre- lation function decays exponentially, implying the absence of single-particle superfluidity. In (b), we show an example for the counterflow superfluid phase at = 0:5, t = 0:02U, and U12 = 0:2U. The anti-pair correlation function RA decays algebraically, while the pair correlation function decays exponentially. Single-particle superfluidity is again absent. The algebraic fits deviate from the data aroundji jj 40, due to the boundary condi- tions of our numerical calculations. . . . . . . . . . . . . . . . . . . . . . . . . 58 xi 3.6 (a) KS and KA as a function of U12 as extracted from the fit of the corre- lation functions, RS and RA. The filling is 0.7 and t=U is 0.02. Around U12=U 0:06, the anti-pair correlation function changes from algebraic to exponential decay. This corresponds to the transition from the PSF to SF phase. When RA decays exponentially, KA is formally set to zero. For Ka + Ks . 2, the system has CDW order. Error bars are one standard de- viation uncertainties obtained from the power-law fit to the numerical data. (b) Comparison of KA obtained from our RG and TEBD calculations. The red square connected by lines are the TEBD results while all other lines are determined from the RG flow with flow parameter l = 3;4;7, and 10, where l is defined in Eq. 3.10. The error bars are as in panel (a). The PSF-to-SF transition obtained from TEBD is around U12=U t 0:06, while the RG cal- culation shows that for l = 10, the transition occurs near U12=U t 0:01. We interpret the regime between U12=U t 0:06 and 0:01 the cross-over region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 xii 3.7 Phase diagram for a homogeneous system with 80 sites and the hopping parameter t = 0:02U as a function of filling and inter-species interac- tion U12=U. The horizontal axis shows three disconnected regions in U12=U. The solid lines are the estimated phase boundaries based on the TEBD re- sults and the dotted line is the PSF-to-SF phase boundary predicted by our RG calculation (see Eq. 3.20). For attractive interaction U12 . 0:06U, the system forms a paired-superfluid (PSF). The state collapses(CL) for U12 . 0:7U. For U12 & 0:06 and U12 . U the system shows single-particle su- perfluidity (SF). The system phase-separates (PS) for U12 & 1 and forms two single-particle superfluids (SF). Open circles are the points whereKS+KA < 2 and charge density wave (CDW) order coexists with a superfluid phase (SF,PSF, or CFSF). At half and unit filling there exist special phases. For re- pulsive interactionU12 & 0:08U and half-filling, the system forms a counter- flow superfluid (CFSF). For unit filling, we find a Mott-Insulator (MI) phase for interactionsjU12j. U. Finally, in the PS region at half- and unit-filling, the system forms two individual MI states. . . . . . . . . . . . . . . . . . . . 61 3.8 Phase diagram at half-filling as a function of U12=U and t=U. The solid lines are estimated phase boundaries from the TEBD calculation and the dotted lines are the phase boundaries predicted by the RG calculation (see Eqs. 3.20 and 3.24). For large repulsive interaction, the system phase separates (PS) and for large attractive interaction, the system collapses (CL). For moder- ate interactions and for t=U . 0:2, the system shows paired superfluidity (PSF) on the attractive side and counterflow superfluidity on the repulsive side. Both PSF and CFSF can coexist with charge density wave (CDW) order when t. 0:1U. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 xiii 3.9 Density distribution of a trapped system for t = 0:02U. (a) Attractive inter- action U12. The trap frequency is = 1 10 5U and the number of atoms is 20 for each species. For attractive interactions, the density distributions of the two species are identical. For U12 = 0:01U (curve I) the system is superfluid. For U12 = 0:11U (curve II) and U12 = 0:21U (curve III), the system is in the paired superfluid (PSF) state. As U12 becomes more neg- ative the distribution gradually shrinks in size. (b) Repulsive interaction U12 = 0:2U with = 8 10 5U and 30 atoms of each species. The red and green curves correspond to the species, respectively. The density distribu- tion has a ?plateau? with half-filling in the center of the trap. The system is in a counter-flow superfluid (CFSF) state. The two species have weak inter- locked density modulations around half filling. . . . . . . . . . . . . . . . . 66 3.10 Density distribution after a time-of-flight expansion. We assume 87Rb atoms and use an expansion time of 0.03s. The hopping energy is t = 0:02U. Panel (a): For attractive interaction U12, we show the TOF expansion of a SF state at U12 = 0:01U (red line) and of a PSF state at U12 = 0:21U (green line). The two curves correspond to the expansion of the densities shown as curve I and III in Fig. 3.9(a) The trap frequency is = 1 10 5U. Panel (b): For repulsive interaction, we show a TOF expansion of a SF state at U12 = 0:01U and of a CFSF state at U12 = 0:21U. The trap frequency is = 8 10 5U. . 67 3.11 Density distribution of molecules after time-of-flight expansion of state III in Fig. 3.9(a). The expansion time is 0.03s. We assume two hyperfine states of 87Rb. These are converted into Feshbach molecules at T = 0 via a fast ramp across a resonance. We assume a complete conversion. The strongly peaked interference pattern of molecules indicates the presence of a quasi- condensate of pairs. For comparison, we also show the TOF expansion of atoms in the PSF phase for the same parameters. The broad Lorentzian dis- tribution demonstrates the absence of single-particle SF. . . . . . . . . . . . 71 xiv 3.12 Structure factor at filling = 0:3. For U12 = 0:01U the system is in the SF regime (dashed line) and for U12 = 0:07U the system is in the PSF regime (continuous line). Cusps atjkj = 2 only occur for U12 = 0:07U indicating the coexistence of CDW with PSF order. . . . . . . . . . . . . . . . 71 3.13 Structure factor S+(k) (blue line) after applying a =2 pulse in the CFSF phase. The quasi-condensate of anti-pairs generates an algebraic peak atk = 0. The cusp also appear in the Fourier transform of the anti-pair correlation function Ra(i;j) =hby1;ib2;iby2;jb1;ji(red dashed line). . . . . . . . . . . . . . . 72 4.1 Noise correlations in different phases, derived from Luttinger liquid theory. In the MI state (column (a)), -function like correlations along k = k0 in G11(k;k0) are visible, whereasG12 nearly vanishes. In the SF state (column (b)), with Luttinger parameters KA = 1:03 and KS = 0:96, we find various contributions in G11, especially - function along k = k0. In Fig. 4.2, we show the contour plots forG11 andG12 for the same state, where we can see the negative correlations at k = 0 and k0 = 0, as well as pairing correlation along k = k0, which is similar to the single-species result in Ref. [35]. G12 shows similar features, but the bunching contribution is an algebraic peak, rather than a -function. In (c) we show an example for the PSF phase, with KA = 0:01 and KS ? 1:3, in (d) an example for the CFSF phase, with KS = 0:01 and KA ? 1:2. In the PSF state, the inter-species correlation G12(k;k0) has strong correlations along k = k0, a reflection of pairing. In the CFSF state, the peak is formed along k = k0 direction, an indication of the anti-pairing (particle-hole) formation in the CFSF state. . . . . . . . . . . 77 xv 4.2 Noise correlations,G11(k;k0) (a) andG12(k;k0) (b), in the SF state. The data used here are the same as those used in Fig. 4.1 (b). We create the non- linear gray scales by plotting functions tanh(500G11) (a) and tanh(105G12) (b), in order to magnify the details around k = k0 = 0. The labels of the color-bar reflect the values ofG11(k;k0) andG12(k;k0). In these plots, we can clearly see the negative correlations between the quasi-condensate (k = 0 and k0 = 0) and the higher momentum states at the quantum depletion, as well as the anti-pairing correlation along k = k0 and the pairing correlation along k = k0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.3 Noise correlations for a homogeneous system of 40 lattice sites, calculated with the TEBD method. The frames (a) ? (d) correspond to the examples (a) ? (d) in Table 4.1. In (a), we show the noise correlations of a MI state. In the plot ofG11(k;k0), there is a strong correlation along the direction k = k0, whereas the noise correlation functionG12(k;k0) essentially vanishes. In (b), we show the noise correlations of a SF state. Here, we can see the peak around k = k0 corresponding to the - function bunching peak predicted by Luttinger Liquid theory (see also Fig. 4.1 (b)). ForG12, we find negative value at k = k0 = 0, which is different from the Luttinger Liquid result (Fig. 4.1 (b)). Other structures predicted by Luttinger Liquid theory can be seen in Fig. 4.4, whereG12 andG11 are plotted in a non-linear color scale to magnify the structures around k = k0 = 0. In (c) and (d), we show the noise correlations of the PSF and CFSF state, respectively. In the PSF state (c), the inter-species correlationG12(k;k0) has strong correlations along k = k0, a consequence of pairing (see also Fig. 4.1 (c)). In the CFSF state (d), the peak is formed along the direction k = k0, an indication of anti-pairing in the CFSF state (see also Fig. 4.1). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 xvi 4.4 Noise correlations,G11(k;k0) (a) andG12(k;k0) (b), in the SF state of a homo- geneous system. The values ofG11(k;k0) andG12(k;k0) are exactly the same as in Fig. 4.3 (b). We create non-linear gray scales by plotting tanh(10G11) and tanh(200G12) in linear scales. The labels of the color-bar reflects the val- ues ofG11(k;k0) andG12(k;k0). The features aroundk = k0 = 0 are magnified as a result of the non-linear scale. In (a), we find the features predicted by LL calculations (4.2 (a)). In addition, we can see a weak correlation at around k = k0 2kF , wherekF = =aL = 0:5 =aL. This is where a strong correla- tion (cusps) will develop if CDW order is present. This feature can also been shown in Luttinger Liquid calculations at around k t 0 and k0 t 2kF (Eq. 4.23). In (b), we find that the structures along k = k0 is similar with the ones in Luttinger Liquid calculations, however, the structures along k = k0 is negative, different from the Luttinger Liquid predictions (see also Fig. 4.2). The difference may be understood as a result of different boundary condi- tions used for the finite-size calculations: the numerical calculations use a ?hard-wall? boundary condition, whereas the Luttinger Liquid calculations assume a periodic boundary condition. . . . . . . . . . . . . . . . . . . . . . 84 4.5 Left : The correlation function C11(q), as defined in Eq. 4.2; right: the struc- ture factor S(k) for a quasi-supersolid state. The parameters are given in example (5) of Table. 4.1. The Luttinger parameters are KA 1:4 and KS 0:57. The filling fraction is = 0:2, hence the "Fermi wave vector" kF is 0:2. At momentum 2kF, both quantities develop cusps, indicating the presence of CDW order. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 xvii 4.6 Noise correlations in the trapped system. The system size is 80 sites and t=U = 0:02. In (a), the system is in the SF state. The particle number of each species is 30, the trap frequency 8 10 5U and U12=U = 0:01. In (b), the particle number of each species is 40, the trap frequency 1 10 4U and U12=U = 0:11. The system has both MI and PSF orders. The MI state forms a plateau at unit-filling at the center of the trap and the PSF is formed at the edge. The PSF state at the edge causes the small peak along the k = k0 direction, similar to the one in (c). However, this peak is at a much smaller amplitude than the one shown in (c), where the whole system is a PSF state. In (c), the particle number of each species is 20, the trap frequency 1 10 5U and U12=U = 0:11. The whole system is in the PSF state. A strong pairing correlation is formed along k = k0 direction. In (d), the particle number of each species is 30, the trap frequency is 8 10 5U and U12=U = 0:2. The system has both CFSF and SF order. The CFSF order forms a plateau at half-filling at the center of the trap and the SF state towards the edges of the trap. The CFSF order causes a strong anti-pairing (particle-hole) correlation along k = k0direction. At the same time, the SF order adds to the "dips" along k = 0 and k0 = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4.7 CorrelationsC12(q) andD12(q) for the states that are described in Fig. 4.7. In (a), we show the behavior of C12(q) (Eq. 4.2) in SF, MI, PSF and CFSF states. The strong anti-pairing (particle-hole) correlations in the CFSF state gives a strong signal around q = 0 in C12(q). This strong signal is also unique to the CFSF state and therefore can be used to detect to the CFSF order. In (b), we show the behavior of D12(q) (Eq. 4.3) in SF, MI, PSF and CFSF states. The strong pairing correlation in the PSF state is the reason for the high peak around q = 0 in C12(q): This suggests measuring C12(q) is a good way of detecting the PSF order. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 xviii 4.8 Noise correlation C11(q) and structure factor S(k) in a PSF/CDW state. The system size is 80 sites and there are 20 particles of each species ( = 0:25). The trap frequency is = 10 5U , the hopping t = 0:02U and the inter- species interaction is U12 = 0:11U. The density at the center of the trap is roughly 0:45 per site and the cusps are developed around 0:9 . The inhomogeneity of a trapped system means that the ?Fermi wave vector? kF is no longer , where is the average filling of the system. Instead, kF can be evaluated as ncenter, where ncenter is the density at the center of the trap, 90 5.1 Dipole oscillation in the CFSF state. Note that time t is in the unit of ~=J for all plots. At t = 0, the density distributions of both species have the same plateau at half-filling and then the trap of species 1 is briefly perturbed by a displacement of one lattice site [see Fig. 5.2 (a)]. In (a) and (b), we show the density distributions of species 1 and 2 as a function of time after the brief displacement. Species 1 and 2 have exactly opposite oscillation. In (c), we show the time evolution of the total density of both species. There is no oscillation in the total density and the plateau at the center is fixed at one. In (d), we show the time evolution of the relative density between the two species. The relative density reflects the oscillatory motion in species 1 and species 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.2 (a) Trap displacement as a function of time. Only the trap of species 1 is displaced at t = 0. The displacement is one lattice site and is removed at t = 1J=~. (b)-(d) Dipole oscillation of the center of mass in SF, PSF and CFSF states induced by the displacement of species 1?s trap. In the SF state (b), only species 1 are excited at t = 0 as a result of displacement. In PSF (c) and CFSF (d) states, both species are excited at t = 0 as a result of their pairing and anti-pairing order respectively. The red(green) color stands for species 1(2) for all plots in this chapter. . . . . . . . . . . . . . . . . . . . . . . 98 xix 5.3 (a) Trap displacement as a function of time. Here, two traps are displaced by 1 lattice site at t = 0 and brought back to the original place at t = J=~. (b) Dipole oscillation in the PSF state after the same displacement. (c) Trap dis- placement as a function of time. Here, two traps are displaced by 1 lattice site from the center and brought back to the center at t = J=~. (d) Dipole os- cillation of center of mass in the CFSF state after the opposite displacement. The parameters used are shown in Table. 5.1. Note that there is no dipole oscillation for the CFSF state under the same displacement and for the PSF state under the opposite trap displacement. . . . . . . . . . . . . . . . . . . 100 5.4 (a) Trap displacement as a function of time. Both traps are displaced by 1 lattice site at t = 0. In (b), the system is at SF state. Both species show the heavily damped dipole oscillation. In (c), the system at the PSF state and the two species show identical heavily damped dipole oscillations. In (d), the system is in a CFSF state. The dipole oscillation is overdamped in this case and the center of mass stays near the original equilibrium position before the displacement. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 5.5 (a) Trap displacement as a function of time. The traps are displaced in the opposite direction by 1 lattice site at t = 0. (b) Dipole oscillation of the center of mass induced by the trap displacement. In this case, the PSF state prohibit the separation of the two species and the center of mass stays at the original equilibrium position of the trap before the displacement. For both CFSF and SF states, the center of mass of each species moves apart from each other under the displacement. In the SF case, the oscillation converges to the new equilibrium position of the traps after the displacement. But, in the CFSF case, each species keep moving away from each other, until they reach the edge of the CFSF state. In (c)-(e), we show the density distribution at different times for the CFSF state. . . . . . . . . . . . . . . . . . . . . . . . 102 xx 6.1 (Color on-line) Efficiency E as a function of temperature calculated by the SC (red cross), IDMFT (blue triangle) and MC (green circle) methods. The interaction parameters, Ubb and Ubf, are shown in each plot. In (a), the SC calculation differs from the other two methods for T < 1t=kB. For this re- gion, the SC expansion formulas derived here are no longer accurate. In (b)- (f), all three methods give almost identical results. These calculations also show that almost 100% efficiency is reached for relatively strong attraction, Ubf 6t, at low temperature, T 1t=kB. In (b), we blow-up the area inside the black square in (a), which corresponds to the low temperature region, where T = 0:2t=kB to 0:5t=kB. In this re- gion, we find that T1 shows large relative fluctuations (deviation) from the mean value and the mean value of T1 differs relatively greater from T. The extracted temperature T2 however still shows excellent agreement with the input temperature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 xxiii List of Tables 3.1 Definitions of MI, SF, CFSF and PSF orders in terms of the large x behavior of the correlation functions RS(x), RA(x), and G(x), Rn;a(x). A: algebraic decay of the form x ; E: exponential decay of the form e x. A correlation function is said to exhibit quasi-order when it is subject to algebraic decay with < 2. In this system, the algebraic decay for RS, RA and G always has < 2, while Rn;a can have 2. CDW quasi-order exists only when Rn;a is described by < 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.1 The parameters used in the numerical examples and the Luttinger parame- ters extracted from the algebraic fit of correlation functions, RA andRS. The Luttinger parameters are set to zero when the correlations decay exponen- tially. The hopping parameter t is 0:02U for all cases. The parameters are chosen to represent different orders that can exist in this system. . . . . . . 82 5.1 Parameters used to represent the SF, PSF and CFSF states. The hopping parameter is set atU=J = 8 for all cases. The particle number of each species is the same and denoted as Na. stands for the particle number of either species. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 xxiv Chapter 1 Introduction Incredibly rich quantum many-body effects are observed (and predicted) in ultra cold atoms [1, 2, 3]. Study in this area now involves methods and concepts from atomic physics, condensed matter physics, quantum information theory and quantum chemistry. With the development of the technology for cooling and manipulating atoms, dilute gases of alkali or alkaline-earth atoms can be cooled to temperatures below micro-Kelvin. The extremely low temperature leads to the realization of quantum degenerate gases such as the Bose- Einstein condensate [4, 5] and the Fermi degenerate gas [6]. With the ability to adjust the interaction strength through the Feshbach resonance [7] and the possibility of changing the dimensionality with optical potentials and of generating optical lattices [8], the dilute gases can display many-body phenomena that are characteristic of strongly correlated sys- tems, which typically exist only in the dense and strongly interacting quantum liquids of condensed matter or nuclear physics. In principle, an accurate quantum description of such many-body phenomena requires solving the many-body Schr?dinger equation with the inclusion of the scattering interac- tion and the interaction with the optical potential. In practice, the Schr?dinger equation of such large systems is impossible to solve directly. However, the resemblance between ultra-cold atoms in optical lattices and a conventional condensed matter system is often made, with atoms in analogy to electrons and the periodic optical lattice to the periodic crystalline lattice [9, 10]. Such resemblance leads to suggestions that the many-body phe- nomena for ultra-cold atoms can be studied through the methods used in the condensed matter theory, such as non-relativistic field operators and the second quantization of the many-body Hamiltonian. Describing ultra-cold atoms in optical lattices with second quan- tized many-body Hamiltonians has lead to the successful prediction of the quantum phase transition between the superfluid (SF) and the Mott insulator (MI) states as the amplitude 1 of the optical lattice increases [10]. This SF to MI transition has been observed in three [11], quasi-two [12] and quasi-one dimensions [13] in ultra-cold atom experiments. In this thesis, the discussion is focused on two types of ultra-cold atom mixtures in optical lattices in reduced dimensions. The first mixture consists of two species of bosonic atoms of equal mass. We assume that the system is in one dimension, which can be studied efficiently by the Luttinger liquid theory and the time evolving block decimation (TEBD) method. The discussion for this system is carried on in four parts: the introduction to the TEBD method in Chapter 2; the phase diagram and correlation functions in Chapter 3, the noise correlation functions in Chapter 4, and the transport properties in Chapter 5. In Chapter 6, we discuss the second model, which consists of heavy bosonic and light fermionic atoms. The discussion is focused on the pre-formed molecules in such systems. For the rest of this chapter, we would like to give a brief introduction to the basic concepts and methods that are most relevant to discussions in the later chapters. 2 Publications in the PhD work: Chapter 3: ? Counterflow and paired superfluidity in one-dimensional Bose mixtures in optical lat- tices, Anzi Hu, L. Mathey, Ippei Danshita, Eite Tiesinga, Carl J. Williams, and Charles W. Clark, Phys. Rev. A 80, 023619 (2009). Chapter 4: ? Noise correlations of one-dimensional Bose mixtures in optical lattices, Anzi Hu, L. Mathey, Carl J. Williams, and Charles W. Clark, Phys. Rev. A 81, 063602 (2010). Chapter 5: ? Detecting paired and counterflow superfluidity via dipole oscillations, Anzi Hu, L. Mathey, Ippei Danshita, Carl J. Williams, and Charles W. Clark, in preparation. Chapter 6 ? Improving the efficiency of ultracold dipolar molecule formation by first loading onto an optical lattice, J. K. Freericks, M. M. Ma?ska, Anzi Hu, Thomas M. Hanna, C. J. Williams, P. S. Julienne, and R. Lema?nski, Phys. Rev. A 81, 011605 (2010). ? Efficiency for preforming molecules from mixtures of light Fermi and heavy Bose atoms in optical lattices: the strong-coupling-expansion method, Anzi Hu, J. K. Freericks, M. M. Ma?ska, C. J. Williams, in preparation. 3 1.1 Second quantization Theoretical studies on interacting many-body systems almost always start with a second quantized Hamiltonian. Here, we discuss the general procedure of the second quantiza- tion (for detailed discussion, see for example [14]) and then derive several of the many- body models relevant to ultra-cold atoms in optical lattices. Creation and annihilation operators: In general, we can assume that there exists a quantum- mechanical basis that describes the number of particles occupying each state in a complete set of single-particle states and introduce the vector state as jn1;n2;:::;n1i (1.1) where the notation means that there aren1 particles in the eigenstate 1, n2 in the eigenstate 2, etc. We also assume that this basis is complete and orthonormal. The creation and annihilation operators are introduced as operators that create or annihilate one particle from the corresponding eigenstate. For bosonic particles, the creation and annihilation operators follow the same rule as those for harmonic oscillators: h br;bys i = r;s; [br;bs] = h byr;bys i = 0: (1.2) For fermionic particles, the creation and annihilation operators satisfy the anticommuta- tion rules, n cr;cys o = r;s;fcr;csg= n cyr;cys o = 0; (1.3) where the anticommutator is defined by the following relationship, fA;Bg AB +BA: (1.4) Second quantization: For almost all systems with two-body interactions, the Hamilto- 4 nian takes the form H = NX k=1 T(xk) + 12 NX k6=l=1 U(xk;xl): (1.5) Here T(xk) is the single-particle energy, which includes the kinetic energy, the potential, etc. And U(xk;xl) is the energy due to interaction between particles. The second quanti- zation of the Hamiltonian is introduced by applying the field operator, ^ (x) = X j j(x)aj; (1.6) onto the Hamiltonian ^H = Z d3x^ y(x)T(x) ^ (x) + 12 Z Z d3xd3x0^ y(x) ^ y(x0)U(x;x0) (x0) (x): (1.7) Here, the field operator aj can be either bosonic or fermionic and the function j(x) is the single-particle wave function. Note that the relationship represented in Eq. (1.7) is inde- pendent of the specific representation of field operators. For field operators represented by annihilation operators, the equation above leads to a discretized Hamiltonian which we will derive in the following section. A different way of representing the field operators is discussed in the section of Luttinger liquid theory. 1.1.1 Bose Hubbard model As suggested in Ref. [10], the Bose Hubbard model can describe a system of bosonic atoms in the optical lattices. In this case, the field operator ^ (x) is constructed to describe the behavior of atoms in the lowest energy band of the energy spectrum of a single atom in an optical lattice, ^ (x) = X j w(x xj)bj (1.8) where w(x xj) is the Wannier function of the lowest band and bj corresponds to the bosonic annihilation operator at lattice site j. The Hamiltonian for the bosonic system is 5 written as, H = Z d3x^ y(x) ~2mr2 +V0(x) +VT(x) ^ (x) +12 4 as~ 2 m Z d3x Z d3x0^ y(x) ^ y(x0) (x x0) ^ (x0) ^ (x): (1.9) Here, the first term corresponds to the kinetic energy, p2=2m. The function V0(x) is the periodic potential created by the optical lattice and Vt(x) is the external trapping potential. In the simplest case, the optical lattice potential can be written as standing waves of wave vector k = 2 = , with the wavelength of the laser light, V0(x) = 3X j=1 Vj0 sin2(kxj): (1.10) The interaction potential is approximated by a on-site interaction with scattering length as and atomic mass m. After the second quantization, the Hamiltonian is written as, H = X hj;j0i tj;j0byjbj0 + X j jnj + 12U X j nj(nj 1) (1.11) where nj is the number operator at site j, nj = byjbj . The parameter tj;j0 is the hopping parameter between site j and j0, tj;j0 = Z d3xw (x xj) ~ 2 2mr 2 +V0(x) w(x xj0): (1.12) For the simplest case, the sites j and j0 are adjacent sites. The parameter j describes the energy offset of each lattice site, j = Z d3xVT(x)jw(x xj)j2 VT(xj): (1.13) 6 And the parameter U describes the interaction between two atoms at site j, U = 4 as~ 2 m Z d3xjw(x)j4: (1.14) Typically, we consider repulsive interactions, as > 0, and therefore U is always positive. This model can be further simplified if we assume the energy offset due to the external trapping potential is negligible and the system is homogeneous. In this case, and when a grand canonical ensemble is considered, the energy offset term jnj is replaced by the chemical potential term, nj. 1.1.2 Multi-component Hamiltonian For ultra-cold atom systems with internal degrees of freedom, one can derive multi-component Hubbard-type models similar to the way the Bose Hubbard model is derived. The internal degrees of freedom can be the result of mixed species, mixed atomic states, and/or mixed energy bands. Here, we consider mixed species as an example and introduce to denote different species. The bosonic and fermionic field operators are now written as a sum of different species, ^ b(x) = X j; w (x xj)bj; ; (1.15) and ^ f(x) = X j; w (x xj)cj; : (1.16) Here bj; is the bosonic annihilation operator of species at site j and cj; is the fermionic annihilation operator of species at site j. The function w is the single-particle Wannier function of species . The quantization is carried out in a similar way except for the Hamiltonian is now a matrix in internal degrees of freedom. If we assume there are no cross-species correlations, the Hamiltonian can be written as a sum of single-species Hamiltonians, H . The param- eters in each H depends on each species? atomic properties and their Wannier functions. 7 For a given lattice, it is possible that different species experience very different hopping and interaction strengths, as a result of different atomic masses and scattering lengths. If we assume that there are cross-species interactions, i.e. the interaction part of the Hamiltonian is not diagonal, one additional term appears in the Hamiltonian, U 0n ;jn 0;j: (1.17) Here n ;j is either the bosonic or fermionic number operator of species at site j. For a Bose-Fermi mixture, n ;j and n 0;j represent the bosonic and fermionic number operator respectively. If we assume that the species are associated with each other, i.e. the single particle Hamiltonian is not diagonal, two additional terms can occur in the Hamiltonian: one cor- responds to transferring one species to another at one site, t ; 0;jay ;ja 0;j +h:c:; (1.18) and the other corresponds to transferring one species at site j to another species at site j0, t ;j; 0;j0ay ;jaa0;j0 +h:c:: (1.19) Note that such off-diagonal terms can only occur within bosonic or fermionic mixtures, not Bose-Fermi mixtures. In this thesis, the first model we describe is a two-component Bose Hubbard model in which both components have the same hopping and the same repulsive intra-species interactions. The inter-species interaction in this model can be either positive or negative [see Eq. (3.1)]. For the second model, we consider tunable bosonic intra-species interac- tions and tunable inter-species interactions and neglect the quantum effect of the hopping of the heavy bosons. This leads to the spinless Bose-Fermi Falicov Kimball model [see Eqs. (6.1-6.3) ]. 8 1.2 Phase diagram 1.2.1 Off diagonal long-range order First discussed by Penrose in Ref. [15], the quantum coherence in Bose-Einstein conden- sates and the superfluidity in liquid Helium can be characterized by the non-decaying behavior of the single-particle density matrix in the coordinate space representation, i.e. 1(x;x0) =h^ y(x) ^ (x)i constant; (1.20) as jx x0j!1. Off diagonal long-range order (ODLRO) is introduced to describe this non-decaying behavior. A more general discussion on ODLRO is given in Ref. [16]. The discussion in Ref. [16] turns out to be relevant not just to superfluidity, but also the pairing and counter-flow effects to be discussed in later chapters. Here, we consider a many-particle system with fixed number of particles whose den- sity matrix is denoted as , Tr( ) = 1; (1.21) We define the reduced density matrices 1, 2, 3.... as hij 1jji = Tr(aj ayi) hijj 1jkli = Tr(akal ayjayi) (1.22) etc: ; wherei,j,... represent single particle states andai,aj,..., are the corresponding annihilation operators. The density matrices 1, 2, and 3 stands for one-, two- and three- particle re- duced density matrix. Ref. [16], it argues that ODLRO can occur in groups of particles that are composed of bosons and an even number of fermions. In this generalized statement, the Bose Einstein condensate and the superfluid correspond to the ODLRO of m = 1. The superconductivity in electronic systems corresponds to the ODLRO of a special group of m = 2, where the group is composed of one spin-up and one spin-down fermion. 9 1.2.2 Order parameter Also originating from the theory for dilute Bose gases is mean-field theory introduced by Bogoliubov (1947). The key point in mean-field theory is the separation of the condensate contribution from the bosonic field operator. In other words, we can introduce a complex function (x), (x) =h^ (x)i: (1.23) This function (x) is a classical field (not an operator). It has the meaning of an order pa- rameter and characterizes the off-diagonal long-range behavior of the one-particle density matrix, h^ y(x) ^ (x0)i! (x) (x0); (1.24) for largejx x0j. It immediate follows that ifh^ y(x) ^ (x0)i goes to a constant asjx x0j goes to infinity, the order parameter (x) should also approach a constant. In the mean- field theory, the non-zero order parameter is used as an indication of the existence of its corresponding ODLRO. The concept of the order parameter can be applied to study other long-range orders, such as superconductivity. In the BCS theory [17], the order parameter can be written as X k hc k;1ck;2i: (1.25) Here, the index 1 stands for spin down electrons and 2 stands for spin up electrons. If we define the field, ^ (x) = X k e ikxc k; ; (1.26) the order parameter can be written as the expectation value of the pair, h^ 1(x) ^ 2(x)i: (1.27) Combining Yang?s discussion on ODLRO and the order parameter introduced in the 10 mean-field theory, we find that in general one can define an operator to represent to the group of m particles, whose reduced density matrix presents an ODLRO, Om = a1a2:::am: (1.28) Here O is an operator of the product consisting of the annihilation operators, a1, a2, ..., am. For the superfluid and the Bose Einstein condensate of identical bosonic particles, this operator is the bosonic field, OSF = ^ (x): (1.29) In superconductivity, the operator is the pair of one spin-up and one spin-down fermion, ^ 1(x) and ^ 2(x) , OBCS = ^ 1(x) ^ 2(x): (1.30) The correlation function of Om, hOym(x)Om(x0)i, is the m-particle reduced density matrix and the non-decaying behavior of the correlation function signals the existence of the cor- responding ODLRO. For many quantum many-body models, ODLRO occurs only below a certain critical temperature, when the thermal fluctuation is low. This is the case for the Bose-Einstein condensate, superfluid and superconductivity. ODLRO may also disappear at zero tem- perature, as a result of the competition between the kinetic energy, the interaction and the external potential. It is also possible that for a given model, several different ODLROs can exist in different parameter regions and sometimes, different ODLROs can coexist, such as the case for supersolid[18]. How to determine the existence of an ODLRO therefore becomes very important, because the ODLRO characterizes the properties of the state. A parameter map can be drawn to show the parameter regions of different ODLRO and the border where the transition occurs. Such a parameter map is often called the phase dia- gram, which is often used to study the phase transitions for a given model. 11 1.2.3 Superfluid to Mott insulator phase transition The existence of two phases in the Bose Hubbard model can be understood from the two extreme cases: 1) the non-interacting case, where U = 0; 2) the extremely strong repulsive interaction case, where the hopping is completely suppressed, i.e. t = 0. In the case U = 0, the model is quadratic in the operators bj and byj. Transforming the Hamiltonian into the momentum space and diagonalizing the Hamiltonian. We find that the energy spectrum is E = X q [ 2tcos(q )]nq: (1.31) and the ground state of the model corresponds to the state where all the atoms occupy the q = 0 state of the lowest band, j Ni(U = 0) = 1pN 0 @ 1p NL X j byj 1 A N j0i; (1.32) where N is the total number of particles and NLis the total number of lattice site. In the case t = 0, the Hamiltonian is diagonal in the Fock state basis, Qj j ji, where j corresponds to the occupation number at site j. The energy at each site j is a function of j as Ej = j +U j( j 1): (1.33) For the case whereN=NL = 1, the ground state of the system is a MI state with one particle at each site, j N=NLi(t = 0) = NY j=1 byjj0i: (1.34) If we define the order parameter for the SF order as haji, it is easy to show that haji is non-zero in the SF state of Eq. (1.32) and zero in the MI state of Eq. (1.34). These two extreme cases also imply that the phase transition between the SF and the MI states is the result of the competition between the kinetic energy (represented by t), which trys to delocalize the particles and reduce the phase fluctuations, and the combination of the inter- 12 actions and the periodic potential (represented by U), which tries to localize the particles and reduce the density fluctuations. Depending on the values of U and t, either the SF or MI state has lower energy and is energetically favored at very low temperature. The competition between the kinetic energy and the interactions also plays an important role for the pairing and anti-pairing ordering in the mixtures, which is going to be discussed in Chapter 3. For more rigorous studies on the SF to MI phase transition, one can use many differ- ent methods, such as mean-field theory [19, 20] which studies the behavior of the order parameterhaji, the strong-coupling expansion method [22] which studies the energy gap between the ground state in the atomic limit (t = 0), j i and the state with one added particle, Piayij iand one hole, Piaij i, or quantum Monte Carlo simulation which studies either the correlation function or the density fluctuations [21]. Generally speaking, the SF-MI phase diagram as a function of andt=U usually contains loops for the MI state, where the density is fixed at some integer number and is not susceptible to a small change of . This indicates the incompressibility of the MI state. 1.3 Tomonaga-Luttinger liquid theory Quasi-one-dimensional systems can be achieved by loading cold atoms onto a two-dimensional lattice, where the atoms are confined by the cigar-shaped one-dimensional tubes created by the lattice. In these arrays of one-dimensional tubes, interesting quantum many-body phe- nomena have been observed, such as the Tonks-Girardeau gas [23] and the one-dimensional SF to MI transition. The many-body physics in one dimension has attracted much inter- est because it is in many senses quite different from the physics in higher dimensions. In one-dimensional systems, the quantum fluctuations play an much enhanced role and there is no true long-range order even at zero temperature. Instead, one can define quasi- long range order through the asymptotic behavior of various correlation functions. The- oretically, these correlation functions and other observables are calculated based on the Tomonaga-Luttinger liquid theory. In this subsection, we briefly describe how Tomonaga- Luttinger liquid theory is applied to one-dimensional bosonic atoms as an example. More 13 detailed discussion about Tomonaga-Luttinger liquid theory and its application in cold atom systems can be found in [24, 25]. We start the discussion by considering the effective Hamiltonian for a homogeneous one-dimensional system of bosonic atoms in lattices, H = Z dx ~ 2 2m@x ^ y(x)@x ^ (x) +V(x) (x) +g2 Z dx^ y(x) ^ y(x) ^ (x) ^ (x); (1.35) where V(x) is the lattice potential in one-dimensional system, V(x) = V0 sin2(kx) and g is the on-site interaction strength. In Tomonaga-Luttinger liquid theory, the Hamiltonian is rewritten in terms of collective variables that describe the density and phase fluctua- tions. The advantage of the density-phase presentation is that the interaction term be- comes quadratic in this presentation, therefore can be treated non-perturbatively. To first illustrate the quadratic form of the Hamiltonian, we consider the case of one-dimensional interacting systems where the lattice is absent, V(x) = 0. The field operator in the density-phase representation up to the leading order is written as ^ (x) = ( 0 + 1 @x ) 1=2X p ei2p( 0x+ (x))ei (x): (1.36) Here and are the long wavelength density and phase fields that obey the standard commutation relation [@x (x); (y)] = i (x y): (1.37) The variable 0 is the average density. The field @x corresponds to the density fluctuation and is assumed to be small compared with 0. Using Eq. 1.36 and retaining only the terms that can become dominant at the low temperature, the Hamiltonian is rewritten as, H0 = ~v2 Z L 0 dx K (@x )2 + 1K (@x )2 (1.38) Here the parameters K and v are dimensionless and depend on both the density 0 and 14 the interaction strength. This Hamiltonian leads to the action, S=~ = 12 K Z dx Z d 1 v(@ ) 2 + (@x )2 (1.39) and the partition function is expressed via functional integrals, Z = Z D (x; ) Z D (x; )eS=~; (1.40) where = @x . The single-particle Green?s function and the density correlations are then calculated based on the partition function and the leading orders in the two functions are written as h^ y(x) ^ (0)i= A x 1 2K +:::: (1.41) h (x) (0)i= 0 + K2 y 2 x2 (y2 +x2)2 +A3 cos(2 0x) 1 x 2K +:::: (1.42) Here Ajis a prefactor. The parameters y and are related to the interaction parameters and the temperature. It is interesting to note that the asymptotic behavior of the correlation functions is solely determined by the parameter K. For non-interacting systems, K = 1and we recover the non-decaying single-particle Green?s function, which indicates the existence of the ODLRO and the system is condensed in the lowest momentum state. As the repulsion increase (K decreases), the correlation function decays faster and the system has less and less tendency towards superfluidity. For purely local interactions,K = 1 is the minimum, corresponding to the infinitely large local interaction. In this limit, the system is the Tonks-Girardeau gas [26]. Now we consider the effect of adding an optical lattice. Here the existence of a lattice introduces an underlying periodicity for the density, characterized by the lattice constant aL and the quasi-momentum, q = 2 =aL. In the interaction process, the conservation of momentum now takes the form k1 + k2 k3 k4 = q. This can be understood as the particles transferring momentum back and forth with the lattice. Such process is known 15 as the umklapp process. Now in this case, the term ei2p( 0x) in Eq. 1.36 can become non- oscillating for p6= 0 when p 0x becomes an integer. This is possible in the lattice system, because the position x is discretized as jaL , where j stands for the j th site in the lattice. When p 0 = na 1L is satisfied, the corresponding term in the Hamiltonian becomes non- trivial and it leads a new term in the Hamiltonian, HL/g0 Z dxcos [2p (x)]: (1.43) For the Bose-Hubbard model with onsite interaction,HL is non-trivial only for the casep = 1 and it corresponds to an integer number of atoms per site, i.e. 0a = n, with n being any integer. The parameter g0 depends on the optical lattice potential, g0 0V0(V0= )n0 1, where is the chemical potential and n0 is the number of bosons per site. This Hamiltonian takes the form of the sine-Gordon model and can no longer be com- puted exactly. Treating HL perturbatively and use the renormalization procedure, one can find that the parameter K is now determined by flow equations. The behavior of the flow equations reveals that the value ofK can either converge to a fixed pointK , in which case HL is irrelevant, or continue to grow, in which case the system goes to the strong coupling region. Along the separatrix between these two regions, K converges to the fixed point Kc. For the case where HL is irrelevant, the Hamiltonian is still written in the quadratic form and the asymptotic behavior of correlation functions can still be written in the form of Eq. (1.41) with K replaced by the fixed point K . For the case where the system goes to the strong coupling region, HL can be expanded around = 0, assuming the density is localized at each lattice site, and the second order expansion of HL leads to a ?mass? term in the Hamiltonian. The asymptotic behavior becomes exponential. For the Bose Hubbard model, the separatrix between these two regions at commensu- rate fillings corresponds to the fixed point Kc = 2. If jK 2j> 0 , the system is in the SF region. Otherwise, the system is the MI region. It is also worth noting that the sine- Gordon problem is a two-dimensional (one space one imaginary time) problem and the SF-MI transition described here falls into the Berezinskii-Kosterlitz-Thouless (BKT) tran- sition [27]. The SF to MI transition in one dimension has also been studied numerically 16 [29, 28] 1.4 Detection Methods 1.4.1 Time-of-flight measurement The time-of-flight (TOF) measurement in the context of ultra-cold atom experiments means measuring the absorption image of the density distribution of atoms after the free expan- sion when all confinement, including the trap and the optical lattices, is turned off [1]. It is one of the most frequently used experimental methods to detect the coherence properties of ultra-cold atoms. Here we briefly explain the mathematical description of the time-of- flight measurement and the reason why the time-of-flight image can reflect the momentum distribution before the flight. Let?s assume at t = 0 (in this section t stands for time, not hopping), both the trap and lattices are turned off and all confined atoms are released and start to expand. Because the gas is usually very dilute, we also assume the expansion is ballistic. Ballistic expansion in free space can be described by the time evolution of the field operator, ^ (p;t) 1 (2 )3 Z d3pe i~p2t=2M ^ (p); (1.44) where ^ (p) = 1 (2 )3=2 Z d3re ip r ^ (r) = Z d3r X j w(r Rj)bje ip r = X j ~w(p)e ip Rjbj: (1.45) Here, M is the atomic mass, bj is the annihilation operator at site j, w(r) is the Wan- nier function and ~w(p) is the Fourier transform of the Wannier function. At time t after the expansion, the density distribution at position x at time t is the expectation value of 17 ^ y(x;t) ^ (x;t) h^n(x;t)i = h^ y(x;t) ^ (x;t)i = 1(2 )3 Z d3p1 Z d3p2e i(p1 p2) xh^ y(p1;t) ^ (p2;t)i ? (M=~t)3jw(k)j2G(k); (1.46) where k is related to x by k = Mx=~t, andG(k) characterizes the coherence properties of the initial many-body state, G(k) = X j;k eik (Rj Rk)hbyjbki: (1.47) If we redefine k and Rj as the momentum and the position vectors in the lattice system as k = 2 n=NLaL, and Rj = jaL, where aL is the lattice constant, NL is the total num- ber of lattice sites and n = 0; 1;::M, the function G(k) above is exactly the momentum distribution of the atoms in the lattice. The relationship between the density distribution after the time-of-flight and the mo- mentum distribution before the time-of-flight often leads to the simplified relationship h^n(x;t)iTOF h^n(k)itrap; (1.48) where the position x is related to the momentum k before the expansion with k = Mx=~t. This relationship is usually sufficient to estimate the main features in the time-of-flight measurement. 1.4.2 Noise correlation In the time-of-flight measurement, each experimental image corresponds to a single real- ization of the density, not the expectation value. Moreover, each pixel in the image records on average a substantial number N of atoms, while in each single realization of an exper- iment, the number of atoms on the pixel exhibits shot-noise fluctuations of relative order 18 1=pN [1]. These shot-to-shot fluctuations can be characterized by the density-density correlation functionG (x;x0), G (x;x0;t) = h^n (x;t)^n (x0;t)i h^n (x;t)ih^n (x0;t)i: (1.49) Here and are the species indices for a system of multiple species. For the system of one species, and can be omitted. Similar to the relationship between the time-of-flight density distribution and the momentum distribution before the expansion, the density- density correlation (or noise correlation) function after the time-of-flight is related to the density-density correlation in the momentum space before the expansion, G (x;x0;t) h^n (k)^n (k0)itrap h^n (k)itraph^n (k0)itrap; (1.50) where, again, k = Mx=~t. The noise correlation function, G (x;x0;t) , can be particu- larly useful for detecting certain strongly correlated orders that can not be revealed in the averaged density distribution image. This is because the noise correlation includes contri- butions from higher-order correlations as h^n (k;t)^n (k0;t)i X j;k eik (Rj Rj0)+ik0 (Rk Rk0)hby jb j0by kb k0i: (1.51) This relationship between the noise correlation function and the higher-order correla- tion functions was first explained in Ref. [30]. Since then, such analysis has been used to experimentally demonstrate the SF to MI phase transition [31, 32] and the formation of fermionic pairs [33]. Theoretically, it has been shown that the noise correlation functions can also be used to detect the charge density wave order [34, 35, 36, 37, 38]. Noise correla- tion functions are calculated in one- and two-component Fermi systems [36, 38], as well as bosonic systems in the hard-core limit [40]. 19 1.4.3 In situ measurement Although the time-of-flight image and its density fluctuations can reveal rich information about the long-range order, the measurement can only extract information from the mo- mentum space of the lattice system. The direct observation of atoms and their density distribution in situ has been a challenge for many years. Recent progress in this area has been made: the site-resolved optical imaging of single atoms has been demonstrated in lattices with large spacing [41] and in sparsely populated one-dimensional arrays [42]; imaging of 2D arrays of ?tubes? with large occupations with an electron microscope [43] and optical imaging systems [44]. Most recently, the density profile up to a single site resolution was measured experimentally for a thin layer of atoms in a two dimensional lattice [45] and the high-resolution measurement of the density pro- file in both SF and MI states in two dimension have been also achieved [46, 47]. The in situ measurements make it possible to directly study the observables in real space and real time. It makes it possible to directly study the onsite number fluctuations. It provides new ways of studying the phase diagram. For example, the variance of the local chemical potential can be controlled through the trapping potential. In the trapped system, the density develops a plateau corresponding to integer filling fractions in the MI state and the width of the plateau is related to the width of the MI loops in the phase diagram [47]. It is also possible to directly observe the dynamics. In Refs. [45, 46], the atomic tunneling between the lattice site is directly recorded and the tunneling rate is measured for different lattice depths. 20 Chapter 2 Time evolving block decimation method 2.1 Basic TEBD method 2.1.1 State decomposition The time evolving block decimation (TEBD) method is originally developed in the context of quantum information, as a method to simulate quantum computations that involve only a limited amount of entanglement [48]. In the context of quantum information, this method describes a chain of qubits coupled to their neighboring qubits. The total Hilbert space H is a product of the Hilbert space of each qubit, which are two-dimensional vector spaces with the qubit states, j0i and j1i, as the orthonormal basis. The TEBD method is soon extended to simulate a wide range of one dimensional quantum many-body models outside quantum information and have shown to be an effective tool to study the static and dynamical properties of quantum many-body models [50]. Here, we discuss how the TEBD method is applied to simulate one-dimensional Hubbard- type models. The Hubbard-type models are similar with qubit systems in the sense that its Hilbert space H can be decomposed as a product of local Hilbert spaces, H = Ml=1Hl: (2.1) Here, l refers to the lth lattice site and M is the number of sites. The local Hilbert space at site l, Hl , has a local dimension of d, d = 2 for fermionic systems and d is infinite for bosons. The orthonormal basis in H is expressed as the tensor product of the Fock state on 21 each site,jj1ijj2i jjMi. Any statej iin H is represented as a superposition of the basis, j i = dX j1;j2;:::;jM=1 cj1;j2;:::;jMjj1ijj2i jjMi; (2.2) where cj1;j2;:::;jM is the superposition coefficient. The dimension of the superposition coef- ficients grows like dM and are very expensive to calculate for a large system. For example, in a system of 100 lattice sites with a local dimension of 2, the dimension of the Hilbert space is as large as dM = 2100 1033. It is necessary to introduce some kind of truncation that can greatly reduce the coefficient space, yet effectively simulate the system. In the TEBD method, this truncation is inspired by the measure of entanglement. Ac- cording to the Schmidt decomposition (SD), the decomposition of a statej iwith respect to the bipartite splitting of the system at site l, i.e. [1;:::;l 1;l] : [l + 1;l + 2;:::;M] is written as j i = lX l=1 [l] lj A lij B li: (2.3) Here, the vector [l] l and l are the coefficients and rank of the Schmidt decomposition. The statesj A liform an orthonormal basis for subsystem A of lattice sites [1;:::;l 1;l]. The statesj B liform an orthonormal basis for subsystem B of lattice sites [l + 1;l + 2;:::;M]. The Schmidt coefficients [l] l satisfy X l j lj2 = 1 (2.4) and h A lj i= [l] lj B li: (2.5) The Schmidt coefficients are closely related to the eigenvalues of the reduced density ma- 22 trices of subsystems, A = TrB(j ih j) = X l [l] l 2j A lih A lj; (2.6) and B = TrA(j ih j) = X l [l] l 2j B lih B lj; (2.7) where TrB denotes the partial trace over subsystem B and TrA the partial trace over sub- system A. The Schmidt rank l is often regarded as a measure of entanglement in quantum in- formation theory [51, 48]. Large values of l usually means more entangled subsystems A and B. For l = 1, the state j i can be written as a tensor product of j Ai and j Bi, which means subsystems A and B are separable (not entangled). For all possible bipartite splittings, the maximum possible value of l is dM=2, when the bipartite splitting is at the center of the system. State decomposition: Based on the Schmidt decomposition, the TEBD method introduce the following way to decompose the statej iof Eq. 2.2, cj1;j2;:::;jM = X 1=1 X 2=1 X M 1=1 [1]j1 1 [1] 1 [2]j2 1 2 [2] 2 [M 2] M 2 [M 1]jM 1 M 2 M 1 [M 1] M 1 [M]jM M 1 : (2.8) This decomposition is composed of M tensors { [1],..., [M]} and M 1 vectors { [1],..., [M 1]}. The index jl in [l] corresponds to the local basis statejjliand take values in {1, ..., d 1} and the indices l in s and s takes values in {1,..., }. Expect for [1] and [M], all s are rank-3 tensor with the dimension of d . The matrices [1] and [M] have the dimension of d. The vector [l] has elements, which are ordered in a decreasing order, i.e. [l]1 [l]2 ::: [l] . The total number of parameters in all s and s is (M 2)d 2 + (M 1) + 2d , which scales like M(d 2 + ) for large M. For an arbitrary state j i; can scale exponentially with M and the total number of 23 parameters after the decomposition is comparable to dM. However, for states where scales linearly with M, this decomposition requires only ploy(M) parameters and is much more efficient than the direct diagonalization. This means that the state decomposition in the TEBD method is efficient for a less entangled state, where the Schmidt rank dM=2. Fortunately, for certain one dimensional many-body models, the ground state and low energy excitation states satisfy such condition. In [49], it is shown that the Schmidt coefficients [l] decay roughly exponentially with index for the ground and low energy excitation states, [l] exp( K ): (2.9) This state decomposition directly gives the Schmidt decomposition. For a state j i represented by Eqs. 2.2 and 2.8, the Schmidt decomposition at site l is given by, j i= X l=1 [l] lj A lij B li; (2.10) where the statesj A liandj B liare given by, j A li = X 1=1 X 2=1 X l 1=1 [1]j1 1 [1] 1 [2]j2 1 2 [l 1] l 1 [l]jl l 1 ljj1ijj2i:::jjli; (2.11) and j B li = X l+1=1 X l+2=1 X M 1=1 [l+1]j1 l l+1 [l+1] l+1 [l+2]j2 l+1 l+2 [M 1] M 1 [M]jM M 1jjl+1ijjl+2i:::jjMi: (2.12) Besides its relationship with Schmidt decomposition, the state decomposition also pro- vides an efficient way for local operations. In the following section, we discuss how local 24 operations are carried out for the decomposed state. 2.1.2 Single-site and two-site operations For most one-dimensional Hubbard-type models, the Hamiltonian is written as a sum of only single-site and two-site operations, H = MX l=1 Hl = MX l=1 (K[1]l +K[2]l ); (2.13) where K[1]l is some single-site operator at site l and K[2] is some two-site operator that involves usually nearest-neighbor or next-nearest-neighbor sites. Here, we only discuss the single-site and double-site operations, but the method can be easily extended to the three- or higher-site operations. However, as we shall see in the discussion that the com- putational cost grows rapidly for higher-site operations. It is reasonable to expect that the TEBD method is most efficient for models with localized operations. Single-site operation: For a single-site operator O at site l, O = X jl;kl Ojlkljjlihklj; (2.14) only [l] needs to be updated, and the new state is given by [l]jl l 1 l = dX kl=1 Ojlk0 l [l]k1 l 1 l: (2.15) Nearest-neighbor operation: For a operator Q of two neighboring sites, l and l + 1, Q = X Qjljl+1klkl+1jjlijjl+1ihkl+1jhklj; (2.16) we need to update the tensors [l], [l+1] and [l]. To explain the update, we first rewrite the statej iinto four state, the states before site l, on site l, on site l+ 1 and after site l+ 1 : 1) an orthonormal basis for sites [1:::l 1] represented byj i, 25 j i = X 1=1 X 2=1 X l 1=1 [1]j1 1 [1] 1 [2]j2 1 2 [l 2] l 1 [l 1]jl 1 l 1 jj1ijj2i:::jjl 1i; (2.17) 2) the Fock state basis for site l,jji; 3) the Fock state basis for site l + 1,jki; 4) an orthonormal basis for sites [l + 2;:::;M], represented byj i, j i = X l+2=1 X l+3=1 X M 1=1 [l+2]j1 l+2 [l+2] l+2 [l+3]j2 l+2 l+3 [M 1] M 1 [M]jM M 1jjl+1ijjl+2i:::jjMi: (2.18) Together, the statej iis written as j i= X ; ; =1 dX j;k=1 [l 1] [l]j [l] [l+1]k [l+1] j ijjijkij i: (2.19) If we introduce a rank-4 tensor jk = X =1 [l 1] [l]j [l] [l+1]k [l+1] ; (2.20) the expression forj iis further simplified as j i= X ; =1 jk j ijjijkij i: (2.21) Now when the two-site operator Q is applied, Q is acting on the rank-4 tensor . Let e be the updated and renormalized after applying Q, 26 e jk = P j0;k0Q jk j0k0 j0k0 Pj0;k0Qjkj0k0 j0k0 : (2.22) To update [l], [l+1] and [l] one just need to find e [l], e [l] and e [l+1] that satisfy, e jk = X =1 [l 1] e [l]j e [l] e [l+1]k [l+1] : (2.23) The new tensors e [l],e [l] and e [l+1]can be found through the singular value decomposition (SVD) of e . Specifically, we construct a d d matrix T, T(j 1) + ;(k 1) + = e jk (2.24) and through the single value decomposition we obtain, T = U V ; (2.25) where U is an d d unitary matrix, is an d d diagonal matrix with non-negative real numbers on the diagonal, and V is the conjugate transpose of V, an d d unitary matrix. The SVD is performed numerically with only the first columns of U, the first diagonal elements of (in decreasing order) and the first rows of V as the output result, T(j 1) + ;(k 1) + = P =1 U(j 1) + ; V ;(k 1) + : (2.26) The values of e [l]j ;e [l] ; e [l+1] are determined as, e [l]j = U(j 1) + ; [l 1] ;e [l] = ; e [l+1] = V ;(k 1) + [l+1] : (2.27) Next-nearest-neighbor operator: In the TEBD simulation, the next-neighboring operator cannot directly be apply to its corresponding sites, instead the site in between should also be included. To apply one next-neighboring operator as a two-site operator while main- 27 taining the tensor product sequence in Eq. 2.8, totally three nearest-neighbor operations are needed. To illustrate the three operations, let?s assume V be an operator acting on site l and l + 2, V = dX kl;kl+2=1 Vjljl+2klkl+2jjlijjl+2ihkl+2jhklj: (2.28) The three steps of applying V on the statej iare: 1) Exchange the state at site l + 1 with the state at site l + 2 by applying the swapping operator at the site l + 1. The swapping operator at an arbitrary site l is defined as Q[l]swap = X j0l;jl;j0l+1;jl+1 jl;j0l+1 jl;j0l+1jjlijjl+1ihj0l+1jhj0lj; (2.29) where x;y is the Dirac delta function (x y). Since the swapping operator is a nearest- neighbor operator, it is directly applied to the state. 2) Apply V as a nearest-neighbor operator on the swapped state (note now the state at site l + 1 corresponds to the state at site l + 2 before the swapping), V = dX kl;kl+2=1 Vjljl+2klkl+2jjlijjl+1ihkl+1jhklj: (2.30) 3) Exchange the updated state at site l + 1 with the state at site l + 2 again by applying the swapping operator Q[l+1]swap. After the three steps, we effectively applyV onto the state at sitel and at sitel+2. Note that the stat at site l+ 1 is also updated as a result of swapping. For further discussion, see also Ref. [52]. 2.1.3 Time evolution In the TEBD method, both the ground state and the dynamic calculation are based on the application of the time evolution operator, or the propagator, which is decomposed into a series of local operations. This can be done because the Hamiltonian consists only local operations, i.e. 28 H = X l Hl;l+s; (2.31) where s is the number of sites around site l. If s = 1, the Hamiltonian involves only single-site and nearest-neighbor operations. If s = 2, the Hamiltonian involves single- site, nearest-neighbor and next-nearest-neighbor operations. To illustrate the basic steps involved in the time-evolution simulation, we will first explain the case of s = 1, where the Hamiltonian involves only single-site and nearest-neighbor operation. Let the time propagator be exp( iHt) with the Plank constant ~ = 1. Because Hl;l+1 does not necessarily commute with each other, the time propagator can not be written directly as a product of exp( iHl;l+1t), exp( iHt) = exp( i X l Hl;l+1t)6= Y l exp( iHl;l+1t): (2.32) But notice that Hl;l+1 and Hl+2;l+3 always commute with each other. If we divide the Hamiltonian into the even site part Heven and the odd site part Hodd, H = Heven +Hodd = X evenl Hl;l+1 + X oddl Hl;l+1: (2.33) every Hl;l+1 commutes with each other within Heven and Hodd. Using the second order Suzuki-Trotter expansion [54] we can rewrite the propagator for a very short time t as follows, e iH t = e iHodd t=2e iHeven te iHodd t=2 +O( t3): = Y oddl e iHl;l+1 t=2 Y evenl e iHl;l+1 t Y oddl e iHl;l+1 t=2 +O( t3): (2.34) In this way, exp( iH t) is decomposed into a product of local operations e iHl;l+1 t , each of which can be applied onto the statej ithrough a nearest-neighbor operation. For our calculation, t is on the order of 0:01 and the error introduced by this expansion is on the 29 order of 10 6. To further reduce the error, one use higher-order Suzuki-Trotter expansions [54]. Imaginary-time evolution: the imaginary time propagation method states that the ground state is the state that the imaginary-time evolution of a superposition of eigenstates con- verges into as the imaginary time, = it, goes to infinite, j igs = lim !1 exp( H )j intikexp( H )j intik : (2.35) In the TEBD method, the initial superposition of eigenstates can either be the input state, or the state after applying the propagator with a very large t 10. We assume such states are a superposition of many eigenstates, including the ground state, of the system. 2.1.4 Observables To numerically study the static and dynamical properties, it is essential to calculate various quantities, such as the total energy , the density distribution and varies correlation func- tions. Generally speaking, there are two ways for calculating the quantities in the TEBD method: 1) through calculating reduced density matrix and 2) through the inner tensor product. For quantities that involve only single-site or local sites, it is often more efficient to calculate it through the reduced density matrix. For quantities that involves many sites or far-apart sites, the inner tensor product of the whole state is more efficient. As an ex- ample, we explain the calculation of the energy E, the density distribution (x) and the correlation functions. The calculation of other quantities can be easily extended from the examples. Single-site observables: for the density distribution and other single-site observables, they are calculated through the reduced density matrix at the corresponding site. Let l be the reduced density of state at site l, l is directly calculated from [l], [l]and [l+1], l = X j; j+1 [l] j [l]jl j j+1 [l+1] j+1( [l] j [l]j0l j j+1 [l+1] j+1) jjlihj0lj: (2.36) 30 The density n(x) at site l is then calculated as n(x) =h j^nlj i= Tr(^nl l); (2.37) where ^nl is the number operator at the site l. Two-site observables: For the total energy E and other observables that involves only neighboring site operators, the calculation is based on the reduced two-site density matrix, l;l+1, which can be written in terms of the rank-4 tensor in Eq. 2.20 , l;l+1 = X l; l+2 jljl+1 lal+2 j 0 lj 0 l+1 l l+1 jjlijjl+1ihj0l+1jhj0lj: (2.38) The total energy E is a sum of local two-site operations over all lattice sites, E = X l El = X l h jHl;l+1j i= Tr(Hl;l+1 l;l+1): (2.39) For an observable involving sitesl andl+k, wherek> 1 , the reduced two-site density of state can be found as a inner product of the s and s betweenl andl+k. For very large k; this calculation can be numerically expensive. Correlation functions: for the calculation of the correlation functionhOyl1Ol2i, one first ap- ply the single-site operator Ol2 ontoj iwhich updates the state toj 1i(non-normalized) and then apply Oyl1 ontoj 1iand update the state toj 2i. The correlation function is then the tensor inner product betweenj iandj 2i, h jOyl2Ol1j i=h jOyl1j 1i=h j 2i; (2.40) where j 1i= Ol1j i;j 2i= Oyl2j i: (2.41) This method can be easily adapted fornth order correlation functions,h jOylnOln 1::::Oyl2Ol1j i, by applying operators in the sequence of Ol1 , Oyl2, .., Oln 1, ; Oylnontoj i, 31 h jOylnOln 1::::Oyl2Ol1j i = h jOylnOln 1:::Oyl2 (Ol1j i) = h jOylnOln 1::: h Oyl2 (Ol1j i) i ::: = h j n OylnOln 1::: h Oyl2 (Ol1j i) io : (2.42) In our calculation, this method have been applied to calculate the pairing, anti-pairing and the four point correlations. Now, we have explained all basic components of the TEBD method and in the next subsection, we will discuss some of additional methods to improve the efficiency and extend the application. 2.2 Extended TEBD method 2.2.1 Number conservation In the TEBD method, the computational cost is mainly determined by the size of the trun- cated Hilbert space. The truncated Hilbert space can be further reduced by taking into account the existence of conserved quantities in the studied systems. In our work, we are interested in the ultra-cold atoms in a confined lattice system, it is reasonable to assume that the total number of atoms in the system is conserved. Here we briefly describe how the number conservation can be implemented in the TEBD method. In the Bose-Hubbard model, the orthonormal basis for the local Hilbert space at an ar- bitrary site l is the Fock statejjli, where jl corresponds to the number of particles at site l. Let the total number to be fixed be atN, the number conservation require that the or- thonormal basis for the system,jj1ijj2i:::jjMi, satisfyPljl =N. This is easy to satisfy for the input state in the TEBD method. The question is how to preserve the number conser- vation through the time propagation. Now because the propagation is decomposed into a series of local two-site operations (see Eq. 2.34), the number conservation is preserved as long as the two-site operation conserves the total number of particles. 32 To implement the number conservation in the two-site operation, it is obvious that the local two-site Hilbert space should only contain the states whose particle number equals to the total particle numberN minus the particle number of the rest of the lattice. Now the problem here is that the local states are linked with the states on the left and right through tensor products. The numbers of particles on the right and left are not numbers, but a vector. To illustrate this problem, let us first label the lattice sites from the left to right with numbers 1, 2, 3,....M and consider a two-site operation at an arbitrary site l. The statej i is written as product of statesj i,jji,jki,j ias in Eq. (2.21). Let NL[l]( ) be the total number of particles in statej iand let NR[l+1]( ) be the total number of particles in statej i. The vector NL[l] contains all the possible total number of the particles on the left of siteland the vectorNR[l+1] contains all the possible total number of the particles on the right of site l + 1. The statesjjiandjkimust satisfy j +k +NL[l]( ) +NR[l+1]( ) =N: (2.43) Now let the total number of particles in statej iand statejjibe NJ( ;j) j + NL[l]( ) and the total number of particles in statej iand statejkibe NK( ;k) = k +NR[l]( ). The biggest and smallest possible value in N1 are NJ;M = max [NJ( ;j)] = min h max(NL[l]) + (d 1);N min(NR[l+1]) i ; (2.44) NJ;m = min [NJ( ;j)] = max h min(NL[l]);N max(NR[l+1]) i ; (2.45) where d is the local dimension for all sites. Similarly, we find the biggest and smallest possible values of N2 are NK;M = max [NK( ;k)] = min h max(NR[l]) + (d 1);N min(NL[l+1]) i ; (2.46) NK;m = min [NK( ;k)] = max h min(NR[l]);N max(NL[l+1]) i : (2.47) There are totally (NJ;M NJ;m+1) possible values for the total particle number for state 33 j i and state jjicombined. We therefore reconstruct e into (NJ;M NJ;m + 1) matrices, B[v], where v = 1;2::;(NJ;M NJ;m + 1). The matrix element B[v]mn is equal to e jmkn m n where fjm; m;kn; ngbe themth combination ofj and that satisfyj+NL[l] = v+NJ;m 1 and the nth combination of k and that satisfy k+NR[l] =N v NJ;m + 1. We also store the indicesfjm; m;kn; ngin matrices LL[v] and LR[v] as LL[v]m;1 = jm; LL[v]m;2 = m; LR[v]n;1 = kn; LR[v]n;2 = n: (2.48) The SVD of B[v] yields a series of unitary matrices U[v], V [v] and the diagonal matrices S[v] whose diagonal elements are the single values for B[v], B[v] = U[v]S[v]V [v] : (2.49) Now, we have a series of matrices S[v] whose diagonal elements are the values fore [l]. But we cannot assign the values in S[l] to e [l] arbitrarily, because the elements in e [l] must be in the decreasing order. In other words, the value of e [l] can only be th largest value in all S[l]s. Its corresponding column vectors in U[v] and V [v] should be used to update e [l]j and e [l]k . To satisfy this requirement, e [l] , e [l]j and e [l+1]k are determined in the following way. First, let us assume S[y]x is the th largest element in all S[v]. We then immediately have e [l] = S[y]x : (2.50) The xth column of matrix U[y] is then used to update [l]j , [l]j = U [y]mx [l 1] (2.51) where the indices j and are determined by LL[y] , 34 j = LL[y]m;1; = LL[y]m;2: (2.52) Similarly, the xth row of matrix V [y] is used to update [l+1]k , [l+1]k = V [y]xn [l+1] ; (2.53) where the indices k and are determined by LR[y], k = LR[y]n;1; = LR[y]n;2: (2.54) For the elements ine [l],e [l] ande [l+1] that are not updated in the above procedure, they are all set to zero. Tensorse [l],e [l] ande [l+1] are the updated [l], [l] and [l+1] and the number conserved two-site operation is completed. Repeat the procedures for each application of exp( iHl;l+1 t), the total number is conserved during the whole time propagation opera- tion. 2.2.2 One dimensional two-component Hamiltonian The straightforward way of presenting the two-component system is to double the local Hilbert space dimension, but this is not favorable because of the high computational cost for increasing the local dimension. A more efficient way is to map the two-component Bose-Hubbard model (Eq. 3.1) onto the one-component Hamiltonian with next-nearest- neighbor hoppings and nearest-neighbor interactions, H = t 2N 2X l=1 (bylbl+2 + h:c:) +U12 X oddl nlnl+1 +U2 2NX l=1 nl(nl 1); (2.55) 35 whereN is the number of sites in the original two-species Hamiltonian. In this one-species Hamiltonian, there are 2N sites, each of which is indexed byl. The odd siteslcorrespond to species 1 and the even sites to species 2. Hopping between neighboring sites tbya;iba;i+1 in Eq. 3.1 is mapped onto a next-nearest-neighbor hopping tbylbl+2 in Eq. 2.55. Similarly, the inter-species onsite-interactionU12n1;in2;i is mapped onto the nearest-neighbor interaction U12nlnl+1. This type of mapping has been successfully applied to treat the two-legged Bose-Hubbard model [53]. The reason why the mapping reduces the cost is explained as follows. The cost in the TEBD method scales as Md3 3. For the two-species system with N sites, M = N, and with a local dimension of D for each species, i.e. d = D2 for two species, the cost scales as MD6 3. On the other hand, for the mapped Hamiltonian with 2N sites, M = 2N, and with a local dimension ofD, the cost only scales as 2ND3 . In our calculation, we setd = 3 or 5. The mapping makes the computation five to ten times faster. To apply the time propagation of the mapped Hamiltonian, we need to split the Hamil- tonian into three parts as H = Hint +Hoddhop +Hevenhop , where Hint = NX m=1 [U12n2m 1n2m +Un2m 1(n2m 1 1) +Un2m(n2m 1)]; (2.56) Hoddhop = t X oddm (by2m 1b2m+1 +by2mb2m+2 + h:c:); Hevenhop = t X even m (by2m 1b2m+1 +by2mb2m+2 + h:c:): Subsequently, we use the second-order Suzuki-Trotter expansion to decompose e i ^H as e iH = e iHint =2e iHoddhop =2e iHevenhop e iHoddhop =2 e iHint =2 +O( 3); (2.57) Each of the operatorse iHint =2,e iHoddhop =2, ande i ^Hevenhop can be decomposed into a product of two-site operators, which can be efficiently applied to the matrix product state j i. 36 We use swapping techniques to apply the next-nearest-neighbor operators e iHoddhop =2 and e i ^Hevenhop . 2.3 Efficiency and accuracy estimation Based on our previous discussion on the TEBD method, it is easy to see that the complexity of the TEBD method is proportional to M 3d2, where M is the system-size and 3d2 is the complexity of a two-site operation. For a given system-size and a local dimension, the computational cost of the TEBD method depends largely on the value of . To reduce the computational cost, it is ideal to set to a small value. From the viewpoint of accuracy, larger is favorable because it includes more entan- gled states into the decomposition. In fact, for a system of system-size M and fixed local dimension d, the TEBD method is exact if is fixed at the maximum value of dM=2. But this value of makes the TEBD method even less efficient than direct diagonalization. The optimal value of is crucial for achieving the desired accuracy with minimum computa- tional cost. Usually, we determine the accuracy by studying the dependence of on . Numerically, we find that decreases very fast as increases. For the value of chosen in our calculation, the value ofj jis less than 10 12. For models involving bosonic particles, where a cut-off is needed for the local dimen- sion d, the accuracy of the calculation also depends on the value of d. The accuracy can be estimated by compare the calculations with local dimensions of d and d + 1. When d is large enough for the state considered, there is very small difference in the calculation results between d and d+ 1. 2.4 Example: one-dimensional Bose Hubbard model To demonstrate the capabilities of the TEBD method, we would like to discuss the compu- tation of the ground state of one dimensional Bose-Hubbard model in different parameter regions. As explained in Chap. 1, the model is written in terms of the bosonic creation and annihilation operators byi and bj for particles on site j, 37 H = t M 1X j=1 (byjbj+1 +h:c:) + MX j=1 Unj(nj 1) MX j=1 (j j c)2 n j: (2.58) Here we assume a closed lattice system of M lattice sites and the lattice constant a. The operator nj is the number operator on site j, nj = byjbj; The parameter t is the hopping parameter, U is the onsite interaction parameter. As a result of the implementation of the total number conservation, the chemical potential is set to zero. Instead, we use N to present the total number of the particles in the system. The term (j jc)2 corresponding to the harmonic trap that is located at the center of the lattice, jc = (M + 1)=2. As is discuss in the previous chapter, there are two phases that exist in the Bose- Hubbard model: Mott insulator (MI) and superfluid (SF). In this section, we will discuss the TEBD calculation of the ground state in the MI and SF regions in one dimension. We will show the calculation of the one body density matrix Gi;j = hbyibji , the momentum distribution n(k) and the on-site number fluctuations j = p(nj hnji)2. 2.4.1 Homogeneous System when the trap frequency = 0, the model corresponding to a homogeneous lattice system with hard-wall boundary conditions. In the SF region, the number fluctuations is very high. The asymptotic behavior of the one-body reduced density matrix, Gij = hbyibji, is algebraic, Gij ji jj 1=K, where K is the Luttinger parameter (Chapter 1). In the MI region, the number fluctuation is low and the Green?s function decays rapidly. The asymptotic behavior of Gij is exponential, Gij exp( ji jj), where is the decay constant. As an example, we choose a system of 50 lattice sites and totally 50 particles. The local dimension for each site d is 5, which corresponds to the maximum occupation number of 4. The parameter is set to 50. In Fig. 2.1, we show the behavior of the number fluctuation, and density distribution, hnji, as a function of location j in the SF and MI region. The density is at unit filling for both cases, while the fluctuation in the SF region is much higher than the one in the MI 38 10 20 30 40 500 0.2 0.4 0.6 0.8 1 1.2 j Paricle Number 10 20 30 40 500 0.2 0.4 0.6 0.8 1 1.2 j Particle Number 10 20 30 40 500 0.2 0.4 0.6 0.8 1 1.2 j Particle Number 10 20 30 40 500 0.2 0.4 0.6 0.8 1 1.2 j Particle Number ? nj? D j (a) (b) (c) (d) Figure 2.1: Number fluctuations j and density distribution hnji for U=t = 1(a), 4(b), 10(c), 20(d). The total number of particle is 50. The density distribution is uniformly 1 except for at the boundary. In (a) and (b), the system is in the SF state, this is reflected by the large number fluctuation. In (c) and (d), the system is in the MI state and the number fluctuations are much reduced. region. In Fig. 2.2 and 2.3, we discuss the behavior of the one body reduced density matrix for U=t =1, 2, 10 and 20. For U=t = 1 and 2, the system is in the SF state. In both cases, Gij decays algebraically. For U=t = 2, Gij has a larger decay power. This is the result of the increased U=t in (b) compared with (a). As U=t continues to increase, Gij becomes short-ranged and the decay behavior changes to exponential. This is the case for U=t = 10 and 20, where the system is in the MI state. In Fig. 2.3, we study the decay behavior of Gij by considering i = 25 and j = 25;26;::50. We can use the power law and exponential fit function to extract the decay parameters. From the decay power, one can determine the Luttinger parameter K. We find that K decreases as U=t increases till around 2, where the transition to MI happens and the decay behavior changes to the exponential decay. In the MI region, Gij decays very fast from the diagonal j = i. This decay can approximated by a exponential function, exp( x). The decay constant becomes bigger for stronger U=t. 39 10 20 30 40 50 1020 3040 500 0.5 1 1.5 10 20 30 40 50 1020 3040 500 0.5 1 1.5 10 20 30 40 50102030 4050 0 0.5 1 1.5 10 20 30 40 5010 2030 4050 0 0.5 1 1.5 Gi,j Gi,jGi,j Gi,j j j j j i i i i (a) (b) (c) (d) Figure 2.2: One-body density matrix Gij =hbyibjifor U=t = 1(a), 4(b), 10(c) and 20(d). The parameters are the same as in Fig. 2.1. In (a) and (b) the system is in the SF state. In (c) and (d) the system is in the MI state. In (a), Gij decays very slowly and the correlation function is forced to zero by the boundary condition at the edge of the lattice. In (b), Gij decays faster than in (a) but still qualifies as algebraic decay (see also Fig. 2.3) with a larger decay power. In (c) and (d), Gijbecomes short-ranged and the decay behavior changes to exponential. 40 0 5 10 15 20 250.5 0.6 0.7 0.8 0.9 11.1 |i?j| G(|i?j|) 0 5 10 15 20 2510 ?10 10?5 100 |i?j| G(|i?j|) 0 5 10 15 20 25 10?10 10?5 100 |i?j| G(|i?j|) 0 5 10 15 20 250.1 0.3 0.5 0.70.9 1.1 |i?j| G(|i?j|) exp(?x/X0) |x|(?1/K) data (a) (d)(c) (b) Figure 2.3: Algebraic and exponential fit for Gij. Here we choose i = 25, j = 25;:::50. The function G25;j represent the decay behavior of the off-diagonal one body density matrix. Through the fit, one can obtain the Luttinger parameter K, K = 9:4 in (a) and K = 6:4 in (b). The decay length X0 can be found for the exponential decay, X0 = 1:5 in (c) and X0 = 0:75 in (d). 2.4.2 Harmonic trap The trap introduces a local potential variance on each site. To understand the effect of the local potential, one can use the local density approximation and treat local potential variance as the chemical potential at sitej, j = (j jc)2. The state at sitej is then the state of the homogeneous system with chemical potential of (j jc)2 and the same U and t. In the SF-MI phase diagram, the MI state exists within lobes of and t=U. In each of the MI lobe, the density is fixed as a integer filling as changes. This means that the local compressibility @n=@ j should be zero. In other words, in MI state, the density is held at an integer filling despite of the trap. The density in the SF state on the other hand will vary over lattice sites and reflect the trap potential. In Figs. 2.4 and 2.5, we show the number fluctuations, the density distribution and one body density matrix for the trapped system at U=t = 1 and U=t = 20. The parameters for this calculation is = 0:016 and N = 25. The trap frequency is strong enough to confine the particle within the trapping potential and hard-wall boundary conditions can 41 10 20 30 40 500 0.5 1 1.5 2 j number of particle 10 20 30 40 500 0.5 1 1.5 2 j number of particle D j ? nj? (a) (b) Figure 2.4: Density distribution and number fluctuations in a trapped system. In (a), the system is in the SF state with U=t = 1 In (b), the interaction is strong, U=t = 20. The density forms a plateau near the center of the trap. This is where the MI state is formed. The fluctuations are reduced at the plateau area as the result of the MI order. The SF state still exists at the edge of the trap. be neglected. ForU=t = 1, the whole system is in a SF state. The number fluctuation is high [Fig. 2.4 (a)] and the functionGi;j shows long-range correlation [Fig. 2.4 (b)]. ForU=t = 20, the MI state is at the center of trap, while the SF state is in the edge. The incompressibility of the MI state leads to the plateau at the center of the trap in the density distribution [Fig. 2.4 (b)]. The difference between the MI and SF states can be seen behavior of the number fluctuation in Fig. 2.4 (b) and in the one-body density matrix Gi;j in Fig. 2.5 (b). 2.4.3 Accuracy estimation As is discussed in section 2.3, the truncation parameter plays an important role of de- termine the accuracy of the TEBD method. In Fig. 2.6 and 2.7, we show examples of [l] as a function of l and for the SF and MI states discussed here. We find that in the SF region, [l] decays to the order of 10 4 when approaches to = 50 and in the MI region, [l] decays to the order of 10 10 as approaches = 50. The decay becomes faster and faster as U=t increases. In the limit of U=t!1, the ground of the model is the Fock state, which is the basis vector state for the TEBD method. Even for the SF region, where the ground state is closer to the many-body state of the lowest single-particle momentum, the value of [l] still decays fast (see Fig. 2.7). Throughout the whole parameter region, we find that reasonable accuracy can be achieved with M , where M is the systemsize. 42 10 20 30 40 501020 3040 500 0.5 1 1.5 2 i j 10 20 30 40 5020 40 0 0.5 1 ij (a)G i,j Gi,j (b) Figure 2.5: One-body density matrix in a trapped system. The parameters in (a) and (b) are the same as in Fig. 2.4. In (a), the system is in the SF state. The density matrix is reduced dramatically to zero near the edge of the cloud. In (b), the off-diagonal decay is exponential in the MI region. The decay changes to algebraic as a result of the SF state at the edge. The computational cost for one parameter set is around 3 CPU hours for this system-size, which is much more efficient than the direct diagonalization. In Fig. 2.8, we show the behavior of [l] in a trapped system. In Fig. 2.8(a), the whole system is in a SF state. The value of [l] decays faster and the edge of the trap and the slowest at the center of trap. A slight different behavior is found in Fig. 2.8 (b), where [l] near the center decays much faster than that at the edge. This is the result of a MI state near the center. Overall, we found reasonable accuracy for this calculation as well. 43 10 20 30 40 501 1020 3040 5010?5 100 10 20 30 40 50 1 1020 3040 50 10?10 10?5 100 10 20 30 40 50 1 1 10 2030 405010?10 10?5 100 1 10 20 30 40 50 10 2030 4050 1 10?10 10?5 100 ?[l]? l ?[l]?l ?[l]? l ?[l]? l l l l l ?l ?l ?l ?l (a) (b) (c) (d) Figure 2.6: Vector [l] as a function of l and . The decomposition rank is 50 and the size of the system is 50 lattice sites. In (a) and (b), the value of decreases to around 10 4 for = 50 at the center. In (c) and (d), the value of decays rapidly as increases and is less than 10 9 at = 50. 5 10 15 20 25 30 35 40 45 50110?10 10?8 10?6 10?4 10?2 100 ? ?[25] ? U/t=1 U/t=4U/t=10 U/t=20 Figure 2.7: Decay of [25] in SF and MI states. Here, [25]is chosen because it relates with the Schmidt rank at the center of the lattice, which usually has the slowest decay with regard to . The vector [25] decays much faster in the MI region than in the SF region. 44 0102030 405001020 30405010?10 10?5 100 l? ?[l] ? 0102030 405001020 30405010?10 10?5 100 l? ?[l] ? (b)(a) Figure 2.8: Behavior of [j] for a trapped system. In (a), [l] corresponds to the SF state of Fig. 2.4 (a) and in (b), [l] corresponds to the MI-SF mixed state of Fig. 2.4 (b). The behavior of reflects the orders that exists in the trap. 45 Chapter 3 One dimensional binary Bose mixtures in optical lattices (I): paired and counter-flow superfluidity 3.1 Introduction Bose-Einstein condensation [4] is a fascinating many-body phenomenon, which demon- strates the significance of quantum statistics at low temperature. Identical bosons can oc- cupy the same single particle state and are in fact more likely to do so than classical parti- cles. At a critical temperature, a gas of bosons undergoes a phase transition towards a state in which a macroscopic fraction of the particles occupy the lowest energy state, creating a condensate. Such a state was realized in ultra-cold atom systems in [5], demonstrating that the technology of cooling and manipulating atoms had reached a level of control with which novel states of matter could be generated and studied. In the case of a Fermi gas, the Pauli exclusion principle prevents such a phenomenon to occur, because no single particle state can be more than singly occupied. However, the phenomenon of condensation can still occur in Fermi systems via a different mechanism: fermions can form pairs to create composite bosons. The bosonic particles then form a condensate of pairs. Conventional superconductors, for example, were understood as a condensate of electron pairs [17]. In ultra-cold atoms, fermionic condensates of this type were created in [55]. Interestingly, this mechanism of condensation of pairs is not limited to fermionic sys- tems but can occur in bosonic systems as well. In fermionic systems, formation of Bosonic pairs necessarily occurs before condensation. In bosonic systems this mechanism can be favored energetically, and will typically be in competition with single particle condensa- tion. In [56, 57], two types of composite bosons were predicted for a binary Bose mixture in 46 Figure 3.1: Sketch of a condensate of pairs. Atoms of each species (red/green) pair together and form a paired superfluid (PSF) state. a optical lattice: pairs and anti-pairs. For attractive mutual interactions, a bosonic mixture can form pairs of atoms which then form a paired superfluid (PSF) state, as is visualized in Fig. 3.1. For repulsive interactions, at special fillings, the atoms can form anti-pairs, which can be interpreted as pairs of one atom of one species and one hole of the other species. These anti-pairs can then generate a counterflow superfluid (CFSF) state, visualized in Fig. 3.2. Most of their simulations were performed for two dimensional systems. In one-dimensional gases quantum phases have quasi-long range order (QLRO), rather than true long range order. QLRO of an operator O(x) is defined as follows: The correla- tion function R(x) =hOy(x)O(0)ifalls off algebraically as R(x) jxj 2 asjxj!1with > 0. Various order parameters O(x) will be defined in the text. In contrast in higher dimensional bosonic systems correlation functions can have true long range order, where correlation functions approach a finite value. Power-law scaling in a 1D optical lattice has been observed in [23]. They observed the Tonks-Girardeau regime of strongly interacting bosons. In this chapter we consider a two-component Bose mixture held in an optical lattice that only allows atoms to hop in one spatial dimension. We ask the question of how the super- fluid as well as other phases or orders can be realized. We assume that the two species of the mixture have the same filling , restricted to the range 0 < 1. The phase diagram of these mixtures is determined using Tomonaga-Luttinger liquid theory [24], which gives the universal phase diagram in terms of a few effective parameters. Based on the universal 47 Figure 3.2: Sketch of a condensate of anti-pairs. Here, atoms of one species are strongly anti-correlated with atoms of the other species, creating a counterflow superfluid (CFSF) state. These composite bosons can also be thought of as a pair of one atom of one species and one hole of the other species. phase diagram, we generate the numerical phase diagram using the time-evolving block decimation (TEBD) method discuss in Chapter 2. With these two approaches we find that CFSF can exist for = 1=2 (half-filling) and repulsive interaction, whereas PSF can exist for < 1 and attractive interaction (see also [58]). We also find that charge density wave (CDW) quasi-order can coexist with both PSF and CFSF, as well as single particle superfluidity (SF). The regimes in which CDW and SF quasi-order coexist constitute a quasi-supersolid phase [39, 18]. Similarly, the regimes where CDW and PSF quasi-order coexist is a quasi-supersolid of pairs and in the case of CFSF, a quasi-supersolid of anti-pairs. Previous work has predicted coexistence of CDW and PSF for 1D Bose mixtures [59, 39] and bi-layer 2D lattice bosons with long-range inter- actions [60], and that of CDW and CFSF for 1D Bose-Fermi mixtures [18, 61]. We then address the question whether PSF and CFSF can be realized and detected in experiment. To simulate the effect of a global trap, we numerically study a mixture confined by a harmonic trap and find that PSF and CFSF can indeed exist in such trapped systems. Their existence can be detected through various measurements. The PSF phase can be detected by using a Feshbach ramp, similar to what has been used in BEC-BCS experiments [55], which generates a quasi-condensate signal in the resulting molecules. The CFSF phase can be detected by applying a =2 pulse followed by Bragg spectroscopy. 48 This generates a quasi-condensate signal in the structure factor. Time-of-flight expansion can also be used to show the absence of single particle superfluidity in PSF and CFSF. Measuring the structure factor via Bragg spectroscopy can be one way of detecting CDW order. This chapter is organized as follows: in Section 3.2, we introduce the model that is used to describe the system; in Section 3.3, we use Tomonaga-Luttinger liquid theory to derive the phase diagram. The numerical results are discussed in Section 3.4. Specifically, phase diagrams of the homogeneous system are presented in Section 3.4.1, and the realization and detection of PSF and CFSF are discussed in Section 3.4.2. We conclude in Section 3.5. 3.2 Two-component Bose Hubbard model Ultra-cold bosonic atoms in optical lattices can be well described by Bose Hubbard mod- els [10]. Here, we consider a mixture of two types of atoms confined to a one-dimensional lattice system. The Hamiltonian of such a system is given by H = t X a=1;2 N 1X i=1 (bya;iba;i+1 +h:c:) +U12 NX i=1 n1;in2;i +U2 X a=1;2 NX i=1 na;i(na;i 1): (3.1) We denote the different types of atoms with index a = 1;2, and the lattice site with index i. We assume that the two species have equal particle density 1, the same intra-species interaction U > 0 and hopping parameter t > 0. The inter-species interaction is given by U12. The operators bya;i and ba;i are the creation and annihilation operators for atoms of type a and site i and na;i = bya;iba;i are the number operators. 3.3 Luttinger liquid approach The universal behavior of this system can be found within a Tomonaga-Luttinger liquid description [24]. First, we switch to a continuum description, ba;i!ba(x), and express the 49 operators ba(x) through a bosonization identity, according to Haldane [25, 62]: ba(x) = [n+ a(x)]1=2 X m e2mi a(x)ei a(x); (3.2) where the real-space density of each species is n = =aL and aL is the lattice constant. The lattice sites are at positions x = iaL. This expression is a phase-density representation of the Bose operators, in which the square root of the density operator has been written in an intricate way. The fields 1;2(x) describe the small amplitude and the long wave length density fluctuations. The fields 1;2(x) are given by 1;2(x) = nx+ 1;2(x), where 1;2(x) = Rxdy 1;2(y). The fields 1;2(x) describe the phase, and are conjugate to the density fluctuations 1;2(x). The contact interactions between the densities in 3.1 written in Haldane?s representa- tion generate an infinite series of terms that contain exp(2m1i( nx+ 1)+2m2i( nx+ 2)), where m1 and m2 are some integers. A term of this form can only drive a phase transi- tion, if the oscillatory part 2 m1nx + 2 m2nx vanishes for all lattice sites. This leads to the requirement m1 +m2 = m3, with m3 another integer [59]. As a further requirement, small integers m1 and m2 are necessary, because the scaling dimension of the term scales quadratically in m1 and m2. For the range 0 < 1, we find that there are three different cases: unit-filling ( = 1), half-filling ( = 1=2), and non-commensurate filling ( 6= 1 and 6= 1=2). It can be checked, using renormalization group arguments as below, that higher forms of commensurability do not generate new phases, but that either phase separation or collapse is reached first. Our numerical findings are consistent with this. Non-commensurate filling. The action of the system, assuming a short-range spatial cut- off r0, at non-commensurate filling is given by [24, 62, 39]: S = Z d2r[ X j=1;2 1 2 K (@v j)2 + (@x j)2 +U12aL 2v~ @x 1@x 2 + 2g (2 r 0)2 cos(2 1 2 2)] (3.3) 50 The first line of the action is characterized by a Luttinger parameter K and a velocity v, contained in r = (v ;x). This part of the action, without the coupling between the two fields a(x), generates a linear dispersion ! = vjkj, where . v should therefore be inter- preted as the phonon velocity. The Luttinger parameterK is a measure of the intra-species interaction U. We will be interested in the regime U &t, in which we have approximately [63] K 1 + 8tU sin : (3.4) The velocity v can also be related to the parameters of the underlying Hubbard model by v vF(1 8t cos =U) (3.5) where vF is the ?Fermi velocity? of an identical system of fermions, vF = 2(aLt=~) sin , and kF is the ?Fermi wave vector?, kF = n. Here, ~ is the Planck constant. The two fields a(x) are coupled by the inter-species interaction. The interaction term U12n1n2 in the underlying Hubbard model generates both the term containing @x 1@x 2, as well as the backscattering term [24, 59] containing cos(2 1 2 2). The action S is only well-defined with a short-range cut-off r0. It is proportional to 1=n. At this scale, g is approximately given by g = U12aL=(v~): (3.6) We diagonalize the quadratic part of the action by switching to the symmetric and antisymmetric combinations S=A = 1p2( 1 2). For the two sectors we find KS=A = (1=K2 U12aL=(v~ K)) 1=2 (3.7) as effective Luttinger parameters. To lowest order inU12 this givesKS=A K U12aLK2=(2 v~). The effective velocities are vS=A = vp1 U12aLK=( v~). Collapse (phase separation) of 51 the superfluid phase is when vS=A is imaginary. We note that KS diverges when collapse (CL) is approached, and that KA diverges as the system approaches phase separation (PS). The anti-symmetric sector contains the nonlinear backscattering term cos(2p2 A). To study its effect, we use an RG approach. We renormalize the short-range cut-off r0 to a slightly larger value, and correct for it at one-loop order. The resulting flow equations are given by [24]: dg dl = (2 2KA)g (3.8) dKA dl = g2 2 2K 3 A (3.9) The flow parameter l is given by l = loge r0 0 r0 ; (3.10) where r00 is the new cut-off that has been created in the RG process. The flow equations 3.8 and 3.9 have two qualitatively different fixed points: Either g diverges, which in turn renormalizes KA to zero, or g is renormalized to zero for finite KA = K A. In the latter case, the action S is quadratic in S and A. For the parameter KS, we use the bare value given in Eq. 3.7. As mentioned in the introduction, we can determine the phase diagram by studying the long-range scaling behavior of correlation functions,hOy(x)O(y)i, of various order pa- rameters O(x). In particular, the single-particle superfluid order parameter is OSF = ba(x) with a = 1;2. The CDW order is related to the 2kF wavevector component of the density operator, OCDW = na. PSF is described by OPSF = b1(x)b2(x), and CFSF by OCFSF = by1(x)b2(x). In the homogeneous system, it suffices to study G(x) = hbya(x)ba(0)i;a = 1;2 (3.11) Rn;a(x) = hna(x)na(0)i;a = 1;2 (3.12) RS(x) = hby1(x)by2(x)b1(0)b2(0)i (3.13) 52 U123 2 1 0 1 2 3 1 1.5 2 SF CL PSK PSF SS(CDW) PSF aL/vh Figure 3.3: Phase diagram of a bosonic mixture at non-unit and non-half-filling. For at- tractive interactions U12 and K < 2 the system can form a paired superfluid state, in the regime labeled PSF and PSF(CDW). This phase can coexist with CDW order for weaker in- teractions. For large repulsive (attractive) interactions U12 the system phase separates (PS) (collapses (CL)). For the remaining regime the system shows single particle superfluidity (SF). This can coexist with CDW order, resulting in a quasi-supersolid (SS) regime. RA(x) = hby1(x)b2(x)b1(0)by2(0)i: (3.14) We find that away from collapse (CL) and phase separation (PS), the correlation functions scale either algebraically or exponentially. For algebraic scaling, we have G(x) jxj SF 2; SF = 2 1=(4KS) 1=(4KA) (3.15) Rn;a(x) cos(2kFx)jxj CDW 2; CDW = 2 KS KA (3.16) RS(x) jxj PSF 2; PSF = 2 1=KS (3.17) RA(x) jxj CFSF 2; CFSF = 2 1=KA: (3.18) where the scaling exponents O are determined by KS and KA after the RG flow. For the case that g diverges in Eqs. 3.8 and 3.9 and KA is undefined, these expressions can still 53 be used. We set KA to zero, and find that CDW and PSF are well defined. Hence Rn;a and RS still show algebraic scaling. On the other hand, CFSF and SF become 1and G and RA scale exponentially. We can identify regimes where different scaling exponents are positive based on the relationship between the scaling exponents and KS=A after the flow. This determines the different quasi-long range orders that are present. The resulting phase diagram is shown in Fig. 3.3, as a function K and U12aL=(v~), as appearing in the action in Eq. 3.3. These two parameters determine the initial values of the flow equations through equations 3.7 and 3.6. We can estimate the phase boundary between PSF and SF. For small U12aL=(v~) this boundary is near the point KA = 1 and g = 0. For that limit, Eq. (3.9) can be linearized to dKA dl = g2 2 2 (3.19) and the expression A = 2(1 KA)2 g2 =4 becomes an invariant of the flow. From the properties of the RG flow of a Berezinskii-Kosterlitz-Thouless transition (see e.g. [24, 27]), the phase boundary is given by A = 0 and g < 0. Using the expressions of KA and v in terms of the Hubbard parameters, we estimate the critical interaction U12 for PSF to occur at U12 U c = 32 t2 U2 sin 2( ): (3.20) The phase boundary between supersolid (SS) and SF has been derived in Ref. [39]. Half-filling. In the case of half-filling, another non-linear term has to be introduced in the action Suk = 2guk(2 r 0)2 Z d2rcos(2 1 + 2 2): (3.21) This term describes Umklapp scattering. At the initial cut-off r0 1=n, guk is approxi- mately given by U12aL=v. In addition to the RG flow in the antisymmetric sector we now 54 3 2 1 0 1 2 31 1.5 2 PSSFCL CFSFPSF (CDW) (CDW) PSF CFS F K U12- -- aL/vh Figure 3.4: Phase diagram of a bosonic mixture at half-filling. In addition to the phases that appear in Fig. 3.3, the system now develops a counterflow superfluid (CFSF) phase, which can coexist with CDW order. also have dguk dl = (2 2KS)guk (3.22) dKS dl = g2uk 2 2K 3 S (3.23) in the symmetric sector. Proceeding along the same lines as for the non-commensurate case, we find the phase diagram shown in Fig. 3.4. We estimate the SF-CFSF phase boundary in the same way as the PSF-SF boundary. We find U12 U c = 32 t2 U2 sin 2( ): (3.24) Unit-filling. At unit-filling we have to introduce a term of the form S1 = 2g1(2 r 0)2 Z d2r(cos(2 1) + cos(2 2)): (3.25) 55 The resulting RG flow for this system is given by dguk dl = (2 2KS)guk + 3 g21(KA KS) 2 (3.26) dg dl = (2 2KA)g + 3 g21(KS KA) 2 (3.27) dg1 dl = (2 KS +KA 2 + 3 gukKS +g KA )g1 (3.28) dKA dl = g2 2 2K 3 A g21 16 2 (KS +KA)K 2 A (3.29) dKS dl = g2uk 2 2K 3 S g21 16 2 (KS +KA)K 2 S (3.30) where 3 is some non-universal parameter [64]. The behavior of this set of equations de- pends strongly on the initial value of g1. For small values of g1, four phases can be stable: Single-particle superfluidity, CFSF, PSF and a Mott phase. For large values only single- particle SF and MI are stable. We determine with our numerical approach, that the Hub- bard model falls into the second category, i.e. there is only a single-particle SF and a Mott state at unit-filling. Having established the universal behavior of the system from Tomonaga-Luttinger liq- uid theory, we now want to connect the phase diagram with the parameters in the Hub- bard model. The expressions 3.4 and 3.5, which relate the Luttinger parameter K and the velocity v to microscopic parameters of the Hubbard model, are only approximate, no full analytic expression is known. In addition, only some phase boundaries are predicted re- liably, because we use perturbative RG in the g . We expect that the analytic calculation only predicts the general structure of the phase diagram, as well as the decay behavior of the correlation functions. To obtain the phase diagram in terms of the parameters in the Hubbard model, we need to use numerical methods. The next section describes the numerical determination of the phase diagram. 3.4 Numerical approach We use the time-evolving-block-decimation (TEBD) method to obtain an approximate ground state solution. We considerN lattice sites with hard-wall boundary conditions and express 56 RS(x) RA(x) G(x) Rn;a(x) MI E E E A SF A A A A CFSF E A E A PSF A E E A CDW A or E A or E A or E A ( < 2) Table 3.1: Definitions of MI, SF, CFSF and PSF orders in terms of the largexbehavior of the correlation functions RS(x), RA(x), and G(x), Rn;a(x). A: algebraic decay of the form x ; E: exponential decay of the form e x. A correlation function is said to exhibit quasi-order when it is subject to algebraic decay with < 2. In this system, the algebraic decay for RS, RA and G always has < 2, while Rn;a can have 2. CDW quasi-order exists only when Rn;a is described by < 2. the Hubbard parameters in units of the intra-species interaction U. The number of sites N is equal to 80, unless otherwise noted. In our numerical analysis, we limit the particle number on each site and each species to two for filling 6 0:8 and four otherwise. Once we obtain the ground state, we calculate the energy, density distributions, correlation func- tions, and the structure factor to identify the quasi-long range order and other properties of the ground state. For example, to determine whether a SF, PSF, or CFSF is present, we study the decay behavior of the correlation functions, G(x), RA(x), and RS(x), defined in Eqs. 3.11, 3.14, and 3.13, respectively. If both RA and RS decay algebraically, the system is in a single- particle superfluid (SF) state. If both are exponential, the system is in a Mott insulator(MI) state. If RS or RA decays algebraically, the system is in the PSF or CFSF state, respectively. These relationships are summarized in Table 3.1. In Fig. 3.5(a) and (b), we show the decay behavior of the correlation functions in the PSF and CFSF phase, respectively. As the Hamiltonian is discrete, the correlation func- tions are calculated as discrete functions: G(i;j) =hbya;iba;ji, RS(i;j) =hby1;iby2;ib1;jb2;ji, and RA(i;j) = hby1;ib2;ib1;jby2;ji. For the PSF phase, RA(i;j) decays exponentially, while RS de- cays algebraically. It is also worthwhile to notice that the single-particle Green?s function decays exponentially, implying the absence of single-particle superfluidity. For the CFSF phase, RA decays algebraically while RS decays exponentially. Single-particle superfluid- ity is again absent. Behavior ofKS andKA: We study the decay behavior ofRS andRA in more detail. Using 57 0 204010?15 10?10 10?5 100 R A(i,j) 0 204010?2 10?1 100 R S(i,j) 0 204010?15 10?10 10?5 100 G(i,j) |i?j| |i?j| |i?j| (a) 0 204010?2 10?1 100 R A(i,j) 0 204010?15 10?10 10?5 100 R S(i,j) 0 204010?15 10?10 10?5 100 G(i,j) |i?j||i?j| |i?j| (b) Figure 3.5: Correlation functions RA, RS, and G on a logarithmic scale as a function of distance ji jj. The index i is 40, the center of the 80 lattice sites. The squares are the numerical data. The blue lines are exponential fits to the data and red dotted lines are algebraic fits. Note that the scale of the vertical axis of the graphs differs by orders of magnitude. In (a), we show an example for the paired superfluid phase at =0.3,t = 0:02U, and U12 = 0:16U. RA decays exponentially and RS decays algebraically. The single- particle correlation function decays exponentially, implying the absence of single-particle superfluidity. In (b), we show an example for the counterflow superfluid phase at = 0:5, t = 0:02U, and U12 = 0:2U. The anti-pair correlation function RA decays algebraically, while the pair correlation function decays exponentially. Single-particle superfluidity is again absent. The algebraic fits deviate from the data around ji jj 40, due to the boundary conditions of our numerical calculations. 58 the fit function, c ji jj 2, wherecand are the fitting parameters, we obtain the power- law exponent and, hence, the Luttinger parameters KS and KA based on Eqs. 3.17 and 3.18. In Fig. 3.6(a), we show these KS and KA as a function of U12, for non-commensurate filling. A Luttinger parameter is formally set to zero when its correlation function decays exponentially. For U12 < 0:06U, RA decays exponentially, while for U12 > 0:06U, RA decays alge- braically, and KA. increases as U12 increases. The system undergoes a PSF to SF transition at U12 = 0:06U. On the other hand, KS decreases monotonically for U12 > 0:6U. For U12 < 0:6U the numerics failed to converge to a homogeneous state. This indicates that the system collapses, and we therefore cannot extract a Luttinger liquid parameter. We can observe charge density wave (CDW) order for a range of U12=U in Fig. 3.6. According to Eq. 3.16, this order exists when KS +KA < 2. In fact, it co-exists with the SF, PSF or CFSF order. At half-filling, KS will go to zero at a critical, positive value of U12. This indicates the transition from the SF to CFSF phase. Finite-size effect: The behavior ofKA=S stated above is affected by the size of the system. Finite size effects can ?smooth out? a sudden change in KA=S at the phase transition. This effect can be estimated from the RG flow calculation by integrating Eqs. 3.8 and 3.9 out to a finite value l rather than to infinity. In Fig. 3.6(b), we show an example of a finite-l RG calculation in the vicinity of the PSF-to-SF transition. We see that as l increases, KA dramatically changes for the attractive U12. In the limit of l !1, KA becomes discon- tinuous and ?jumps? from 0 to 1 at U12 t 0:01U. This is where the PSF-to-SF transition occurs. This transition is a Berezinskii-Kosterlitz-Thouless transition [27, 24]. In order to compare the RG result with our TEBD result, we associate the system size N with the flow parameter l, based on the relation in Eq. 3.10. The cut-off r0 is the lattice constant aL and r00 = NaL. For N = 80 we have l = 4:4 and we find that the RG and TEBD are in good agreement. The regime between U12=U 0:06 and 0:01 is a cross-over regime due to the finite size of the system. Collapse and phase separation: For largejU12j, the system approaches collapse or phase separation. According to Tomonaga-Luttinger liquid theory, KS !1 as the system ap- proaches collapse and KA !1 as the system approaches phase separation. As seen in 59 ?0.8?0.6?0.4?0.20 0.20.40.60 1 2 3 4 U12/U K A/S KAK S SFCL PSF?0.06 (a) ?0.10 0.1 0.2 0.3 0.40 0.5 1 1.5 2 U12/U K A l=3l=4 l=7l=10 Numerical (b) U12??0.06U U12??0.01U Figure 3.6: (a) KS and KA as a function of U12 as extracted from the fit of the correlation functions, RS and RA. The filling is 0.7 and t=U is 0.02. Around U12=U 0:06, the anti- pair correlation function changes from algebraic to exponential decay. This corresponds to the transition from the PSF to SF phase. When RA decays exponentially, KA is formally set to zero. ForKa+Ks . 2, the system has CDW order. Error bars are one standard deviation uncertainties obtained from the power-law fit to the numerical data. (b) Comparison ofKA obtained from our RG and TEBD calculations. The red square connected by lines are the TEBD results while all other lines are determined from the RG flow with flow parameter l = 3;4;7, and 10, where l is defined in Eq. 3.10. The error bars are as in panel (a). The PSF- to-SF transition obtained from TEBD is around U12=U t 0:06, while the RG calculation shows that for l = 10, the transition occurs near U12=U t 0:01. We interpret the regime between U12=U t 0:06 and 0:01 the cross-over region. 60 0.20.1 0.30 0.2 0.4 0.6 0.8 1 U12/U 0.60.811.2 0 0.2 0.4 0.6 0.8 1 ?0.2 0?0.1?0.3?1.2?1?0.8?0.60 0.2 0.4 0.6 0.8 1 ? SF CL PSCFSFPSF MI Figure 3.7: Phase diagram for a homogeneous system with 80 sites and the hopping param- eter t = 0:02U as a function of filling and inter-species interaction U12=U. The horizontal axis shows three disconnected regions in U12=U. The solid lines are the estimated phase boundaries based on the TEBD results and the dotted line is the PSF-to-SF phase bound- ary predicted by our RG calculation (see Eq. 3.20). For attractive interactionU12 . 0:06U, the system forms a paired-superfluid (PSF). The state collapses(CL) for U12 . 0:7U. For U12 & 0:06 and U12 .U the system shows single-particle superfluidity (SF). The system phase-separates (PS) for U12 & 1 and forms two single-particle superfluids (SF). Open cir- cles are the points whereKS+KA < 2 and charge density wave (CDW) order coexists with a superfluid phase (SF,PSF, or CFSF). At half and unit filling there exist special phases. For repulsive interaction U12 & 0:08U and half-filling, the system forms a counterflow super- fluid (CFSF). For unit filling, we find a Mott-Insulator (MI) phase for interactionsjU12j.U. Finally, in the PS region at half- and unit-filling, the system forms two individual MI states. 61 Fig. 3.6, we indeed find such a tendency in our TEBD calculations. For U12 > 0:8U (not shown), KA increases rapidly to values around 10, indicating a possible phase separation. ForU12 < 0:6U, due to the slow decay of the correlation functionRS and the finite-size of our system, we are unable to extract an accurateKS from the numerical result. On the other hand, we observe a peaked density distribution for U12 < 0:6U, indicating a collapse. In the phase separation regime, G(x) has algebraic decay except for = 0:5 or 1, where it has exponential decay. An algebraic decay implies two spatially-separated single-species su- perfluids while the exponential decay implies two spatially-separated Mott insulators.[65]. 3.4.1 Phase diagram We study the phase diagram as a function of filling and parameters of the Hubbard Hamiltonian. Assuming a positive U, the system can be fully characterized in terms of , t=U, and U12=U. Our results are shown in Fig. 3.7 for a fixed hopping parameter and in Fig. 3.8 for half filling. 3.4.1.1 Phase diagram at a fixed hopping parameter In Fig. 3.7 we show the phase diagram for filling fractions between 0 and 1 and the inter- action U12=U between -1.1 and 1.1. The symbols correspond to numerical data points at which the phases have been characterized. Different markers represent the different or- ders. The orders are determined from the decay behavior of the three correlation functions RA, RS, and G. For weak attractive inter-species interaction, 0:06 < U12=U < 0, the system is in a SF state. As U12 grows more attractive, paired superfluidity (PSF) occurs. The critical U12 is largest, 0:08U, at half-filling and gradually decreases away from half-filling. This phase boundary differs from that predicted by our RG calculation (Eq. 3.20), plotted as the dotted line in Fig. 3.7. This discrepancy is the result of the finite-size effect discussed in Fig. 3.6(b). In the SF to PSF cross-over regime, charge density wave (CDW) order can coexist. According to the phase diagram Fig. 3.3, for attractive interaction, CDW order can co-exist only with PSF order. In our numerical work, we observed the CDW order 62 slightly outside the numerical phase boundary of PSF but within the RG phase boundary of PSF. The sub-regime where CDW and PSF coexist ends when U12=U . 0:4. When the inter-species attraction is comparable to the intra-species repulsion, U12 . U, the system collapses (CL) and no long-range order is present. For repulsive inter-species interaction and U12 < U, the system is in a SF state for all non-commensurate fillings. Within the SF regime, there is a smaller parameter re- gion where CDW order coexist with the SF order. This subregime is a quasi-supersolid regime. The boundary between a normal superfluid and a quasi-supersolid is estimated by RG calculation in Ref. [39]. At half-filling, counterflow superfluidity (CFSF) occurs when 0:08 . U12=U . 1. Within the CFSF regime, the CDW order can coexist, forming a quasi-supersolid of anti-pairs. It also worthwhile to point out that at half-filling, CDW order only exists within the PSF and CFSF regimes. At unit filling, our numerical results do not show evidence of PSF or CFSF for any U12. We find a Mott insulator (MI) state forjU12j 0:14U. We can compare this phase diagram with the half-filling phase diagram in Fig. 3.4 obtained from Tomonaga-Luttinger liquid theory. Especially, we can compare the loca- tion of the phase boundary between SF and PSF(CFSF). To do so, we plot the RG phase boundaries, described by Eqs. 3.20 and 3.24, onto our phase diagram. The area near the two boundaries is interpreted as the cross-over regime where finite-size effects modify the phase boundary. 63 0.5 1 0 0.05 0.1 0.15 0.2 0.25 ?1 ?0.5 0 0 0.05 0.1 0.15 0.2 0.25 t/U PSF PSSFCL U12 / U CFSF Figure 3.8: Phase diagram at half-filling as a function of U12=U and t=U. The solid lines are estimated phase boundaries from the TEBD calculation and the dotted lines are the phase boundaries predicted by the RG calculation (see Eqs. 3.20 and 3.24). For large repulsive interaction, the system phase separates (PS) and for large attractive interaction, the system collapses (CL). For moderate interactions and for t=U . 0:2, the system shows paired superfluidity (PSF) on the attractive side and counterflow superfluidity on the repulsive side. Both PSF and CFSF can coexist with charge density wave (CDW) order when t . 0:1U. 64 3.4.2 Realization and detection Having established the phase diagram for the homogeneous system, we now discuss how to realize and detect the PSF and CFSF phases. First, we need to modify the Hubbard Hamiltonian in Eq. 3.1 because in any ultra-cold atom experiment an additional trapping potential is present. We add a harmonic potential, (j jc)2(n1;j + n2;j), where j is the site index and jc is the index at the center of the system. The TEBD method is used to find the ground state. We consider a system of 80 lattice sites and adjust the total number of particles and the trap frequency so that the number of particles is negligible at the edge of the lattice. We again determine the orders of the system by studying the correlation functions in Table 3.1. We find that, in spite of the presence of the trap, the correlation functions still show exponential or algebraic scaling away from the edge of the lattice. In fact, a correlation function can have different decay behavior in different parts of the trap. We also find that SF, PSF, and CFSF still exist. The remainder of this article focuses on experimental signatures that distinguish between these orders by calculating the density distribution, the time-of-flight image after an expansion, or the structure factor for Bragg spectroscopy. Density distribution: We find that in a trapped system PSF and CFSF can only exist when the density distribution satisfies certain conditions. For PSF, the density of each species at the center of the trap, ncenter, must be less than one atom per site or equivalently per lattice constant aL. (The density is largest at the center.) For CFSF, ncenter must satisfy ncenteraL = 1=2. Once such conditions are satisfied, the critical value of U12 for PSF and CFSF is close to the one for a homogeneous system (See Figs. 3.7 and 3.8). In Fig. 3.9(a) we show density distributions for three attractive interactions U12 and a hopping parameter equal to the one used for Fig. 3.7. For all attractive interactions, the density distributions of each species are the same. For more attractive inter-species interaction, the density distribution concentrates near the center of the trap. There is no discontinuous change in the density distribution when the system goes from SF to PSF. In Fig. 3.9(b) we show the density distribution for U12 = 0:2U. In this case in the center of the trap, where the density distribution is constant or has a ?plateau?, the system is in a 65 0 20 40 60 800 0.2 0.4 0.6 n(x) (per site) (a) IIIII I x (site) 0 20 40 60 800 0.2 0.4 0.6 n(x) (per site) species 1species 2 (b) x (site) Figure 3.9: Density distribution of a trapped system for t = 0:02U. (a) Attractive interac- tion U12. The trap frequency is = 1 10 5U and the number of atoms is 20 for each species. For attractive interactions, the density distributions of the two species are iden- tical. For U12 = 0:01U (curve I) the system is superfluid. For U12 = 0:11U (curve II) and U12 = 0:21U (curve III), the system is in the paired superfluid (PSF) state. As U12 be- comes more negative the distribution gradually shrinks in size. (b) Repulsive interaction U12 = 0:2U with = 8 10 5U and 30 atoms of each species. The red and green curves correspond to the species, respectively. The density distribution has a ?plateau? with half- filling in the center of the trap. The system is in a counter-flow superfluid (CFSF) state. The two species have weak interlocked density modulations around half filling. CFSF state. The ?plateau? is at half-filling consistent with predictions from a local density approximation and noting that in Fig. 3.7 CFSF only occurs at = 1=2. Towards the edge, where the density is decreasing sharply, it is in a SF state. The plateau implies that the system is incompressible in the center. Time of flight measurement: A widely used measurement technique in the field of ultra- cold atoms is measuring the density of atoms after a time-of-flight (TOF) expansion. The 1D optical lattice potential and the harmonic trap are abruptly turned off at time T = 0 and the atoms expand freely afterwards. We calculate the density at time T, according to na(x;T) =hcya(x;T)ca(x;T)i (3.31) with a = 1;2. The operators ca(x;T) are related to the lattice operator ba;j according to ba(x;T) = NX j=1 w(x rj;T)ba;j; (3.32) 66 ?0.5 0 0.50 0.5 1 n(x) (mm ?1 ) x (mm) SF PSF (a) ?0.5 0 0.50 0.5 1 n(x) (mm ?1 ) x (mm) CFSF (b)SF Figure 3.10: Density distribution after a time-of-flight expansion. We assume 87Rb atoms and use an expansion time of 0.03s. The hopping energy is t = 0:02U. Panel (a): For attractive interaction U12, we show the TOF expansion of a SF state at U12 = 0:01U (red line) and of a PSF state at U12 = 0:21U (green line). The two curves correspond to the expansion of the densities shown as curve I and III in Fig. 3.9(a) The trap frequency is = 1 10 5U. Panel (b): For repulsive interaction, we show a TOF expansion of a SF state at U12 = 0:01U and of a CFSF state at U12 = 0:21U. The trap frequency is = 8 10 5U. wherew(x;T) = q d=p2 (T)2 exp( x2=(4 (T)2)) describes the free expansion from the initial Gaussian wave-function of an atom in a lattice site and (T)2 = d2 +iT~=(2m). The parameterdis the width of the initial Gaussian state andmis the atomic mass. The density distribution na(x;T) is then given by na(x;T) = NX j1;j2=1 w (x rj1;T)w(x rj2;T)G(j1;j2); where G(j1;j2) is the single-particle Green?s function. In Fig. 3.10 we show examples of TOF expansions of PSF, CFSF, and SF order. For the SF phase, we find a strongly peaked interference pattern, reflecting the single-particle quasi-long range order. For both PSF and CFSF phases, the TOF density shows a broad Lorentzian distribution, which is due to the exponential decay of the single-particle Green?s function. Feshbach ramp: In order to detect the superfluidity of pairs, we consider applying a Feshbach ramp to pairwise project the atoms onto molecules formed by one atom from each species, which is similar to detection of fermionic pairs in the BCS regime [55]. In those experiments, a fast ramp across a Feshbach resonance was used, followed by a time- of-flight expansion. The density distribution of the molecules showed the superfluidity of 67 fermionic pairs. We propose a similar detection for bosonic pairs in PSF. To give a simple estimate of a TOF image after a Feshbach ramp, we imagine that bosons of different species on the same lattice site are converted into molecules. This leads to the replacement b1;jb2;j !Mj, where Mj is the molecule annihilation operator. A TOF density of the molecules at position x and time T is given by nM(x;T) = NX j1;j2=1 w (x rj1;T)w(x rj2;T)Rs(j1;j2): (3.33) In the expanding wave functionw(x;T), the massmis replaced by the mass of the molecule. We assume the same initial width d. In a more realistic estimate, the conversion efficiency to molecules would not be 100%, but approximately given by the square of the overlap of the molecular wave function and the single-atom wave functions. This leads to a re- duced signal. The spatial dependence, however, remains the same. In Fig. 3.11, we see an example of the density of molecules after TOF and, for comparison, the atomic den- sity after TOF for the PSF state. The strongly peaked molecular distribution indicates the quasi-condensate of the bosonic pairs. The single-atom density is a broad Lorentzian dis- tribution, indicating the absence of single-particle SF. Bragg spectroscopy: To detect the presence of CDW order, one can use Brag spectroscopy [66, 67]. The quantity that is measured in those experiments is either the dynamic or static structure factor. Here we calculate the static structure factor Sa(k) for species a = 1;2. It is defined as Sa(k) = 1N X j1;j2 e ikaL(j1 j2)(hna(j1)na(j2)i hna(j1)ihna(j2)i): (3.34) For wave vectors k near twice the ?Fermi wavevector? kF, the structure factor S(k) jjkj 2kFj1 CDW with CDW = 2 KS KA [24]. In our system, KS +KA is always larger than 1 and, thus, 1 CDW is positive. Consequently, the structure factor does not diverge. In the CDW regime with KS + KA < 2 the power 1 CDW, however, is less than one. 68 This gives S(k) cusps at 2kF when CDW quasi-long range order is present. In Fig. 3.12 we show examples of S(k) for a case with and without CDW. Bragg Spectroscopy preceded by a /2 pulse: To detect CFSF order, we propose the fol- lowing detection method. It applies to the case that the mixture is composed of atoms in different internal states rather than different atomic species. First, we apply a =2 pulse, which transfers the atoms into the superpositions b1=2;i ! b ;i = (b1;i b2;i)=p2. We then measure the structure factor, which now corresponds to the Fourier transform of the density correlations Rn (i;j) = hn ;in ;ji hn ;iihn ;ji. In terms of the original b1=2;i operators these density correlations are given by Rn (i;j) = 14h(n1;i +n2;i)(n1;j +n2;j)i 14(hn1;ii+hn2;ii)(hn1;ji+hn2;ji) +12hby1;ib2;iby2;jb1;ji (3.35) The last term in the above equation is the correlation functionRa(i;j) of the order parame- ter of CFSF, b1;jby2;j. In Fig. 3.13, we show the structure factor S+(k), the Fourier transform of Eq. 3.35, as well as the Fourier transform of Ra(i;j). Both S+(k) and the Fourier trans- form of Ra(i;j) have a cusp around k = 0. The cusp is due to the long-range correlations of the anti-pairs in the CFSF. The two functions are nearly identical near k = 0, indicating that the momentum distribution of anti-pairs can be measured by determining the struc- ture factor S+(k) . 3.5 Summary We have studied ground state properties of one-dimensional Bose mixtures in an optical lattice using both Tomonaga-Luttinger liquid theory and the time-evolving block decima- tion method. We first discussed the zero-temperature phase diagram in a homogeneous system at different filling fractions and different parameter regimes. We have shown that 1D Bose mixtures in an optical lattice can have quasi-long range orders that include su- perfluid, paired superfluid (PSF), counterflow superfluid (CFSF), and Mott insulator. We 69 also found that each type of superfluid order can coexist with charge density wave (CDW) order and that in both PSF and CFSF phases single particle superfluidity (SF) is absent. In addition, we discussed ways of realizing and detecting these phases experimentally. We propose using a Feshbach ramp to probe the momentum distribution of pairs in the PSF, which shows signatures of the quasi-condensate of pairs. To detect the CFSF for a mixture composed of two atomic hyperfine states, we propose to measure the static struc- ture factor by using Bragg spectroscopy preceded by a =2 pulse. A sharp peak in the structure factor was shown to be dominated by the contribution from the momentum dis- tribution of anti-pairs in the CFSF phase. Finally, we suggest to detect CDW order with Bragg spectroscopy. 70 ?0.5 0 0.50 0.5 1 n(x) (mm ?1 ) x (mm) Single?Particle Pair Figure 3.11: Density distribution of molecules after time-of-flight expansion of state III in Fig. 3.9(a). The expansion time is 0.03s. We assume two hyperfine states of 87Rb. These are converted into Feshbach molecules atT = 0 via a fast ramp across a resonance. We assume a complete conversion. The strongly peaked interference pattern of molecules indicates the presence of a quasi-condensate of pairs. For comparison, we also show the TOF expansion of atoms in the PSF phase for the same parameters. The broad Lorentzian distribution demonstrates the absence of single-particle SF. 0 0 10 20 S(k) SF PSF(CDW) pi/2 pi?pi ?pi/2 kaL Figure 3.12: Structure factor at filling = 0:3. For U12 = 0:01U the system is in the SF regime (dashed line) and for U12 = 0:07U the system is in the PSF regime (continuous line). Cusps atjkj = 2 only occur for U12 = 0:07U indicating the coexistence of CDW with PSF order. 71 00 20 40 60 S +(k) S+(k) FT of Ra ?pi/2?pi pi/2 pikaL Figure 3.13: Structure factor S+(k) (blue line) after applying a =2 pulse in the CFSF phase. The quasi-condensate of anti-pairs generates an algebraic peak at k = 0. The cusp also appear in the Fourier transform of the anti-pair correlation function Ra(i;j) = hby1;ib2;iby2;jb1;ji(red dashed line). 72 Chapter 4 One dimensional binary Bose mixtures in optical lattices (II): noise correlation In Chapter 1, we have discussed the study of noise correlations as a way of probing ultra- cold atom systems. The noise correlations are measured as the spatial correlations of the density in the fully expanded atomic cloud after turning off the trap. If atomic interactions during the expansion can be ignored, one can use the noise correlation measurements to infer momentum space correlations in the initial state. Such an analysis has been used to demonstrate the phase transition between superfluid (SF) and Mott insulator (MI) states [31, 32], as well as the formation of fermionic pairs [33]. In Refs. [38, 36], noise correlations were shown to be an effective probe of such orders in 1D Fermi systems, for both one- and two-component systems. Similar studies have been done for 1D bosonic systems, either in the hard-core limit [40] or using Luttinger liquid (LL) theory [35]. In Ref. [35], the signature of condensates and quasi-condensates was discussed in detail. Noise correlations can also be used to study the phases of binary bosonic mixtures. In the previous chapter, we established the phase diagram of a binary mixture exhibiting SF, PSF, CFSF and MI orders, and we showed that each of the superfluid orders can coexist with the CDW order. We also show that because the PSF and CFSF orders are the result of inter-species pairing, they do not provide a signature in the momentum distributions of the individual atomic species. In this paper, we show that noise correlation measurements provide distinctive signals of both the PSF and CFSF orders. Ref. [71] shows that noise correlations characteristic of the PSF/CFSF orders can be observed even in a system of only four atoms. Here we calculate the noise correlation spectra from the ground state obtain by the TEBD method, which is supported by analytical calculations based on LL theory. We make appropriate comparisons between results for homogeneous and trapped systems. 73 To evaluate the noise correlations, we first assume ballistic expansion and long expan- sion time and define the noise correlations as the density correlations in momentum space, Gaa0(k;k0) = hna;kna0;k0i hna;kihna0;k0i (4.1) where a, a0 are species indices (a;a0 = 1;2), k, k0 are momenta, and na;k and na0;k0 are the occupation operators in momentum space. We also consider the derived quantitiesCaa0(q) and Daa0(q), defined as Caa0(2q) = Z dk hna(k +q)na0(k q)ihn a(k +q)ihna0(k q)i ; (4.2) and Daa0(2q) = Z dk hna(k +q)na0(q k)ihn a(k +q)ihna0(q k)i : (4.3) Each of these quantities can capture the main features of the noise correlations for partic- ular types of order and can be directly measured in experiments [32]. We will present all our results first in the form ofGaa0(k;k0) and then use Caa0(q) and Daa0(q) to highlight the key features. This chapter is organized as follows: in Sec. 4.1, we first summarize the different quasi- orders discussed in the previous chapter and introduce the noise correlations in such mix- ture; in Sec. 4.2, we show our LL results and predict the generic signature of the noise correlations for different orders. In Sec. 4.3, we present our numerical calculation of noise correlations for both homogeneous and trapped systems. In Sec. 4.4 we summarize the main features of the noise correlations for 1D binary Bose mixtures. 4.1 Noise correlations and quasi-orders In previous chapter, we discuss the phase diagram for the binary Bose mixtures in opti- cal lattices and show that the single-species superfluid (SF) has the order parameter ba(x) (a = 1;2) and its corresponding correlation functionG(x) =hbya(x)ba(0)i; the paired super- 74 fluid (PSF) has the order parameter b1(x)b2(x) and its corresponding correlation function RS(x) = hby1(x)by2(x)b1(0)b2(0)i; the counter-flow superfluid (CFSF) has the order param- eter b1(x)by2(x) and its corresponding correlation function RA(x) = hby1(x)b2(x)b1(0)by2(0)i: The CDW order parameter is na(x) (a = 1;2) and the corresponding correlation function Rn;a(x) = hna(x)na(0)i. The asymptotic behavior of the correlation functions at large x is listed in Table. 3.1. We calculate the noise correlations from the four-point correlation function: Gaa0(k;k0) = NX j1;;j2;j3;j4=1 Laa0(j1;j2;j3;j4)ei[kj12+k0j34] hna(k)ihna0(k0)i; (4.4) (4.5) where j12 j1 j2, j34 j3 j4 andLaa0 is the four-point correlation function, Laa0(j1;j2;j3;j4) =hbya;j1ba;j2bya0;j3ba0j4i: (4.6) It is easy to see that the correlation functions RS, RA and Rn;a are the special cases ofLaa0, L12(j1;j2;j1;j2) = RS(j1;j2); L12(j1;j2;j2;j1) = RA(j1;j2); (4.7) Laa(j1;j2;j2;j1) = Rn;a(j1;j2) +na;j1: The noise correlationG12, therefore, contains the Fourier transform of RS and RA, gS(k;k0) = X j1;j2 RS(j1;j2)ei(k+k0)(j1 j2) (4.8) 75 and gA(k;k0) = X j1;j2 RA(j1;j2)ei(k k0)(j1 j2) (4.9) andGaa contains the Fourier transform of Rn;a, gn;a = X j1;j2 Rn;a(j1;j2)ei(k k0)(j1 j2): (4.10) If RS(j1;j2) decays asjj1 j2j 1=KS, we find that gs scales asjk + k0j 1=KS. Similarly, if RA decays asjj1 j2j 1=KA, gA scales asjk k0j 1=KA. For the PSF phase, we find that gs(k;k0) is the dominant term ofG12(k;k0) with a strong peak around k = k0. This peak is the signal of the PSF order. Similarly, for the CFSF phase, we find that the functiongA(k;k0) becomes dominant around k = k0 inG12(k;k0). The peak around k = k0 is the signal of the CFSF order. These remarks are made to give the reader an intuitive interpretation of the relationship between the noise correlations and the long-range orders. In the following section, we will explain the calculation of the noise correlations via LL theory and show that the features mentioned above are indeed reflected in the LL calculation results. 4.2 Luttinger liquid approach In this section we determine the generic behavior of the noise correlations using a Luttinger liquid approach. This formalism has been applied to one-dimensional Fermi systems in Ref. [36], and additionally to single-species bosonic systems in Ref. [35], where a detailed description of these calculations was given. Here, we use an analogous derivation for the case of a bosonic mixture. We outline key steps of the derivation, but refer the reader to Ref. [35] for a detailed description of the method. As described in Chapter 3, we use the Luttinger liquid theory to construct the field operator for both species [Eq. (3.2)] and obtain the action of Eq. 3.3. We calculate the noise correlations at the Gaussian fixed point, corresponding to the SF phase. Here, the system separates into symmetric and anti-symmetric degrees of freedom, defined as S;A = ( 1 2)=p2 and S;A = ( 1 2)=p2. The action can be written either in terms of the 76 000 0.51G11(k,k?) pi ?pipi?pik?aL kaL 00?5 05 10 ?pipipi G11(k,k?) ?pi kaLk?aL 000 0.51G11(k,k?) pi?pik?aL pi ?pikaL 00 00.5 1G11(k,k?) kaL ?pi pi pi ?pik?aL 00?5 0510 15x 10?3 ?pi G12(k,k?) k?aL kaL?pi pipi (a) MI 00?50 510 15x 10?3 ?pi G12(k,k?) pi?pik?aL pi kaL (b) SF 00 0.20.4 0.6 0 G12(k,k?) pi?pik?aL kaL?pipi (c) PSF 00 00.2 0.40.6G12(k,k?) kaL?pik?aL pipi ?pi (d) CFSF Figure 4.1: Noise correlations in different phases, derived from Luttinger liquid theory. In the MI state (column (a)), -function like correlations along k = k0 inG11(k;k0) are visible, whereasG12 nearly vanishes. In the SF state (column (b)), with Luttinger parametersKA = 1:03 and KS = 0:96, we find various contributions in G11, especially - function along k = k0. In Fig. 4.2, we show the contour plots for G11 and G12 for the same state, where we can see the negative correlations at k = 0 and k0 = 0, as well as pairing correlation along k = k0, which is similar to the single-species result in Ref. [35]. G12 shows similar features, but the bunching contribution is an algebraic peak, rather than a -function. In (c) we show an example for the PSF phase, with KA = 0:01 and KS ? 1:3, in (d) an example for the CFSF phase, with KS = 0:01 and KA ? 1:2. In the PSF state, the inter-species correlationG12(k;k0) has strong correlations along k = k0, a reflection of pairing. In the CFSF state, the peak is formed along k = k0 direction, an indication of the anti-pairing (particle-hole) formation in the CFSF state. kaL k?a L 0 0 ?5 ?10 1 5p ?p p ? 10?3 ?p (a) kaL k?a L 0 0 ?2.5 ?0.50 0.5 2.5p(b) ?p ?p p ? 10?5 Figure 4.2: Noise correlations,G11(k;k0) (a) andG12(k;k0) (b), in the SF state. The data used here are the same as those used in Fig. 4.1 (b). We create the non-linear gray scales by plotting functions tanh(500G11) (a) and tanh(105G12) (b), in order to magnify the details around k = k0 = 0. The labels of the color-bar reflect the values ofG11(k;k0) andG12(k;k0). In these plots, we can clearly see the negative correlations between the quasi-condensate (k = 0 and k0 = 0) and the higher momentum states at the quantum depletion, as well as the anti-pairing correlation along k = k0 and the pairing correlation along k = k0. 77 phase fields S = X j=S;A Z d2rj hKj 2 [(@vj j) 2 + (@x j)2] i ; (4.11) or in terms of the fields S;A S = X j=S;A Z d2rj h 1 2 Kj [(@vj j) 2 + (@x j)2] i : (4.12) The velocities vS;A are the phonon velocities of the symmetric/anti-symmetric modes, and rS;A = (vS;A ;x). The parametersKS;A are the Luttinger parameters of the symmetric/anti- symmetric sector. To calculate the noise correlations away from the SF regime, we take the limitsKS;A!0 to describe the phases in which either or bothRS=A have short-ranged cor- relations (exponential decay). This approximation corresponds to the limit that the length scale of the exponential decay is much smaller than any other length scale of the system. Our calculation could be extended in a straightforward way to include a finite decay length of the exponential decay. We start out by calculatinghna;kifor small momentum k 0, for which the Bose oper- ators are given by ba pnei a. Forhna;kiwe find: hna;ki n Z dx12eikx12e 12h( a(2) a(1))2i; (4.13) where a(1) refers to a(x1), and similarly for a(2), and x12 = x1 x2. The correlation functionh( a(2) a(1))2ican be rewritten in terms of correlation functions for S;A. Using the Gaussian action above, we find h( S=A(2) S=A(1))2i = 12K S=A log r 20 +x212 r20 ; (4.14) 78 where r0 is a short-range cut-off. With that we find hnki n Z dx12eikx12F(x12); (4.15) where F(x) = r2 0 r20+x2 g : (4.16) The exponent g is given by g = 1=8KS + 1=8KA. Next we evaluate the expectation value hnknk0ialong the same lines. We obtain: hn1;kn1;k0 i n2 Z eikx12+ik0x34F(x12)F(x34)A; (4.17) where A = (r2 0 +x214)(r20 +x223) (r20 +x213)(r20 +x224) h ; (4.18) and Eq. 4.17 is a volume integral over the spatial variables x12, x23, x34. The exponent h is given by h = 1=8KS 1=8KA. We combine these expressions to get the correlation functionG11(k;k0): G11(k;k0) n2 Z eikx12+ik0x34F(x12)F(x34)(A 1): (4.19) ForG12(k;k0) we proceed analogously, and find h = 1=8KS + 1=8KA. For the finite-size systems that we treat here, we evaluate these integrals numerically, by choosing a finite length L of the system, and by replacing each spatial variable x by (L=2 ) sin(2 x=L) (see Ref. [24]). To compare with the TEBD calculations of a homogeneous system in the next section, we choose the values of the Luttinger parameters, KA and KS to those obtained from the TEBD calculations and they are listed in Table. 4.1. 79 In Fig. 4.1 (b) and 4.2, we show an example for the SF regime. In the upper panel of Fig. 4.1 (b), we show G11(k;k0), in the lower panel G12(k;k0). The Luttinger parameters are KA = 1:03 and KS = 0:96. The ratio L=r0 was chosen as L=r0 = 20, corresponding to the particle number of each species in the numerical example. The shape of G11(k;k0) is the same as the noise correlation function for a single bosonic SF, which was discussed in Ref. [35]. It has the characteristic features of a superfluid: positive correlations along k = k0, which indicates pairing correlations; negative correlations for the axes k = 0 and k0 = 0, indicating the negative correlations between the quasi-condensate and the higher momenta due to pair fluctuations; and bunching correlations along k = k0. ForG12(k;k0) we find qualitatively a similar shape, with the main difference, that the bunching along k = k0 does not have a -function contribution, but only algebraic terms. We note that for a system of two non-interacting species, i.e. for U12 = 0, KS = KA, andG12 vanishes. In Figs. 4.1 (c) and (d) we show the noise correlations for the PSF and the CFSF phase, respectively. For the PSF example, the Luttinger parameters are KS = 1:3 and KA = 0:01. For L=r0 we again pick L=r0 = 20. For the CFSF phase, the parameters are KA = 1:2 and KS = 0:01. In the PSF regime, we find a strong pairing signature in G12, similar to the pairing signature in Fermi mixtures [36, 35]. In the CFSF example, an strong anti-pairing signature is found inG12. We can obtain the functional form of these signatures in the limit L!1, by applying similar arguments to what has been given in Ref. [35]. In the PSF region, We rewrite the noise correlation integral in terms of z = (x12 x34)=2, h+ = (x14 + x23)=2 and h = (x14 x23)=2. We then note that forG12 and forKA!0, the exponenth = 1=8KS+1=8KA diverges. This enforces the integrand to be negligible away from z, h+ 0. Thus the integral evaluates to G12 jk +k0j 1=KS: (4.20) This is the shape that would be approached in an infinite system by the noise correlations shown in Fig. 4.1 (c), lower panel. The deviation from the pure power law is due to the finite size of the system. With similar arguments one can show that in the CFSF regime the 80 inter-species noise correlation approaches G12 jk k0j 1=KA; (4.21) for L!1. Again, the deviation from a pure power law is due to the finite size of the sys- tem. Furthermore, one can show that for both PSF and CFSF orders,G11(k;k0) approaches (k k0), in the limit of infinite size. Equations 4.20 and 4.21 show that there is a simple re- lationship betweenG12 andKS in the PSF regime andG12 andKA in the CFSF regime. This suggests that a careful measurement ofG12 can be used to extract the value of the Luttinger parameters appropriate to the system. This is further confirmed by our numerical calcu- lations for a trapped system, where we show that the algebraic relationship described by Eqs. 4.20 and 4.21 remains valid in the presence of a harmonic trap. We discuss prospects for experimental determination of Luttinger parameters in Sec. 4.3.3. The MI result in Fig. 4.1 (a) is obtained by setting both KA and KS to 0.01. In this case, the ground state closely approximates a simple product of MI states of each species. Thus, G11 approaches a -function, whereasG12 nearly vanishes. Next we calculate the noise correlations for the case k 0 and k0 2kF, where kF is the Fermi wavevector defined above. Essentially the same calculation can be done for k0 2kF, and k0 0 and k 2kF. na;k is still given by the expression (4.15), but na;k0 now needs to be calculated with the operator representation b(x) =pnexp(2i (x)) exp(i (x)). With that we find hnq0+2kFi n Z dx12eiq0x12F0(x12); (4.22) whereF0(x12) has the same form as before but with an exponent g0 = 1=8KS + 1=8KA + (KS +KA)=2. The noise correlations take the form G11(k;q0) n2 Z eikx12+iq0x34F(x12)F0(x34)(A 1): (4.23) 81 Parameter setting Order Luttinger Parameters (a) U12=U = 0:01, = 1 MI KA = KS = 0 (b) U12=U = 0:01, = 0:5 SF KA?1:03, KS?0:96 (c) U12=U = 0:11, = 0:5 PSF KA = 0, KS?1:3 (d) U12=U = 0:11, = 0:5 CFSF KA?1:2, KS = 0 (e) U12=U = 0:26, = 0:2 SF & CDW KA?1:4;KS?0:57 Table 4.1: The parameters used in the numerical examples and the Luttinger parameters extracted from the algebraic fit of correlation functions, RA and RS. The Luttinger param- eters are set to zero when the correlations decay exponentially. The hopping parameter t is 0:02U for all cases. The parameters are chosen to represent different orders that can exist in this system. We therefore note that around the points k 0 and k0 2kF, and k 2kF and k0 0 the integrand is multiplied by a contribution that is of the form of the integrand of the static structure factor S(q) Z eiqx12 r2 0 r20 +x212 (KS+KA)=2 ; (4.24) which can create cusps in the noise correlation when the system is in the CDW regime. These cusps are found in our numerical calculations and are discussed in the next section. 4.3 Numerical approach The calculation of noise correlations is based on the ground state generated by the TEBD method. We set the Schmidt rank = 100 and the local dimension d = 5. We use imaginary-time propagation to generate the ground state. After obtaining the ground state, we calculate the correlation functions, RA, RS, Rn;a and G, and determine the quasi- long range order present in the system based on the relationship shown in Table. 3.1. Furthermore, we can extract the value of the Luttinger parameters, KA and KS, from the numerically calculated correlation functions. We use these parameters in a Luttinger Liq- uid calculation to compare the numerical and the analytical results. The main challenge of determining the noise correlation functions is the high compu- tational cost of calculating the four-point function, Laa0(j1;j2;j3;j4) (Eq. 4.6), which is estimated to scale as 3d3N4. For the system sizes used in this paper, we use parallel com- 82 000 12G11(k,k?) ?pi pipi ?pik?aL ka L 00 ?5 10 05 ?pi G11(k,k?) k?aL kaLpipi ?pi 000 12 3 pi?pi pi ?pi G11(k,k?) k?aL kaL 000 12 3 pi G11(k,k?) ?pi ?pik?aL kaLpi 00?0.05 00.05 ?pi pipi ?pi G12(k,k?) k?aL kaL (a) MI 00?0.05 00.05 ?pi ?pipipi G12(k,k?) k?aL kaL (b) SF 000 0.20.4 pi?pi G12(k,k?) k?aL pi ?pikaL (c) PSF 00 00.10.2 0.30.4G12(k,k?) kaLpi?pi ?pipik?aL (d) CFSF Figure 4.3: Noise correlations for a homogeneous system of 40 lattice sites, calculated with the TEBD method. The frames (a) ? (d) correspond to the examples (a) ? (d) in Table 4.1. In (a), we show the noise correlations of a MI state. In the plot ofG11(k;k0), there is a strong correlation along the direction k = k0, whereas the noise correlation function G12(k;k0) essentially vanishes. In (b), we show the noise correlations of a SF state. Here, we can see the peak around k = k0 corresponding to the - function bunching peak predicted by Luttinger Liquid theory (see also Fig. 4.1 (b)). ForG12, we find negative value atk = k0 = 0, which is different from the Luttinger Liquid result (Fig. 4.1 (b)). Other structures predicted by Luttinger Liquid theory can be seen in Fig. 4.4, whereG12 andG11 are plotted in a non- linear color scale to magnify the structures around k = k0 = 0. In (c) and (d), we show the noise correlations of the PSF and CFSF state, respectively. In the PSF state (c), the inter-species correlationG12(k;k0) has strong correlations along k = k0, a consequence of pairing (see also Fig. 4.1 (c)). In the CFSF state (d), the peak is formed along the direction k = k0, an indication of anti-pairing in the CFSF state (see also Fig. 4.1). puting algorithms to speed up the calculation by parallelization the computation of Laa0 along the indices ji. 4.3.1 Homogeneous systems In this section we discuss the numerical results for noise correlations of a homogeneous system of 40 lattice sites, subject to the hard-wall or ?open? boundary condition, in which the wave function is required to vanish on the fictitious sites of index 0 and N + 1 im- plied by Eq. 3.1. We consider five parameter sets listed in Table 4.1, representing different regimes of the phase diagram of 1D Bose mixtures. 83 k?a L kaL0 0 ?0.25 ?0.05 0 0.05 0.25pi pi (a) ?pi?pi kaL k?a L 0 0 ?12.5 ?2.50 2.5 12.5p p (b) ? 10?3 ?p ?p Figure 4.4: Noise correlations, G11(k;k0) (a) and G12(k;k0) (b), in the SF state of a homo- geneous system. The values ofG11(k;k0) andG12(k;k0) are exactly the same as in Fig. 4.3 (b). We create non-linear gray scales by plotting tanh(10G11) and tanh(200G12) in linear scales. The labels of the color-bar reflects the values of G11(k;k0) and G12(k;k0). The fea- tures around k = k0 = 0 are magnified as a result of the non-linear scale. In (a), we find the features predicted by LL calculations (4.2 (a)). In addition, we can see a weak corre- lation at around k = k0 2kF , where kF = =aL = 0:5 =aL. This is where a strong correlation (cusps) will develop if CDW order is present. This feature can also been shown in Luttinger Liquid calculations at around k t 0 and k0 t 2kF (Eq. 4.23). In (b), we find that the structures along k = k0 is similar with the ones in Luttinger Liquid calculations, however, the structures along k = k0 is negative, different from the Luttinger Liquid predictions (see also Fig. 4.2). The difference may be understood as a result of different boundary conditions used for the finite-size calculations: the numerical calculations use a ?hard-wall? boundary condition, whereas the Luttinger Liquid calculations assume a periodic boundary condition. 84 Superfluid and Mott insulator For the Hamiltonian of Eq. 3.1, in the non-interacting case, U12 = 0, SF and MI are the only two possible orders. In the interacting case, SF and MI orders are still encountered, when the inter-species interaction is weak. For the Hamiltonian of Eq. 3.1 with t U, the MI state exists for anyjU12j.U, until the occurrence of collapse (U12 . U) or phase separation (U12 & U). The SF state however exists only whenjU12j U. In either of SF or MI phases, the quasi-order is formed in each individual species and the cross-species correlation is weak. For the MI state (Fig. 4.3 (a)) we find that G11(k0;k) shows strong correlations along the direction k0 = k, in agreement with the LL theory result shown in Fig. 4.1 (a). We also find that the correlations along k0 = k are not uniform and that the peak along k0 = k resembles a Lorentzian distribution in k imposed upon a constant. This Lorentzian is due to the characteristic scale of the correlation functions. This contribution was ignored in the before-mentioned approximation in the LL calculation, but could be included in a straightforward manner. The cross-species noise correlation,G12(k;k0), on the other hand, is essentially zero, indicating the absence of cross-species correlations in the MI state. For the SF state (Fig. 4.3 (b) and Fig. 4.4), we consider the case where there is weak repulsion between the two species (Table 4.1 (2)). The Luttinger parameters are KA = 1:03 and KS ? 0:96, which were extracted from the correlation functions RS and RA by numerical fitting. From the upper panel in 4.3 (b), we see that G11 has the characteristic features of a quasi-condensate [35]: the positive correlations along k = k0, which indicate pairing; the negative correlations between k = 0 and finite k0, as well as between k0 = 0 and finite k; and a -function like correlation along k = k0, corresponding to bosonic bunching. The lower panel in Fig. 4.3 (b), we see G12, which shows similar features, except for the -function along k = k0, which is "softened" into a power-law divergence and a slight negative value at k = k0 = 0. For a system of two non-interacting superfluids, i.e. U12 = 0, we have KS = KA, andG12 = 0. 85 0.5 1 1.5 2 2.5 qaL C 11(q) ?pi 0 pi ?0.4pi 0.4pi (a) 0 5 10 kaL S(k) 0.4pi?0.4pi pi?pi 0 (b) Figure 4.5: Left : The correlation function C11(q), as defined in Eq. 4.2; right: the structure factor S(k) for a quasi-supersolid state. The parameters are given in example (5) of Table. 4.1. The Luttinger parameters are KA 1:4 and KS 0:57. The filling fraction is = 0:2, hence the "Fermi wave vector" kF is 0:2. At momentum 2kF, both quantities develop cusps, indicating the presence of CDW order. Paired superfluid and counter-flow superfluid We now discuss the noise correlations of the PSF and CFSF states. The noise correlation G12(k;k0) is particularly important for these two phases, because it can verify the existence of PSF and CFSF orders. Unlike SF and MI states, PSF and CFSF states are characterized by order parameters that contain both species and therefore cannot be reflected in any single-species observables, such as the single-particle Green?s function,Ga(x) or the single- particle momentum distribution [70]. The noise correlation function G12(k;k0) measures the correlations between the momentum occupancies of the two species, and thus provides a direct probe of these orders. We have shown in the previous section, that the peak along k = k0 inG12(k;k0) indicates the PSF order and that alongk = k0 indicates the CFSF order. These features are verified in our numerical calculation ofG12 from the ground state. In Fig. 4.3 (c), we show the noise correlations in the PSF state. The parameters are listed in (c) of Tab. 4.1. The existence of PSF order is - as usual - determined by the behavior of the RS(x) and RA(x). RA(x) decays exponentially and RS(x) algebraically with Luttinger parameter KS ? 1:3. For the noise correlation function G12(k;k0), we find that a peak is formed along k = k0, which is a consequence of the pairing correlations. In Fig. 4.3 (d), we show our numerical results for the CFSF example (d) in Table 4.1. Based on the behavior of the RS(x) and RA(x) we verify that the system is in a CFSF state with KA?1:2, and an 86 exponentially decaying RS(x). ForG12 we find that a peak is formed along the diagonal direction, as a result of correlations of anti-pairs (b1by2). These findings are consistent with the predictions of LL theory (see Fig. 4.1). We note thatG12(k;k0) is enhanced in magnitude in the PSF and the CFSF phase compared to the MI and the SF phase, with a strongly altered functional form. Charge density wave In certain parameter regimes of the phase diagram, charge density wave (CDW) order can coexist with each of the three superfluid orders, SF, PSF and CFSF. In Sect. 4.2, we use LL theory to show that CDW order can be reflected in the functionG11 and that the behavior of G11 aroundk = k0 2kF resembles the structure factorS(k). The reason for the resemblance can be understood in a simple way, by recalling the definition of the structure factor S(k) = 1N X j1;j2 e ik(j1 j2)(hnj1nj2i hnj1ihnj2i): (4.25) As mentioned in Sect. 4.1 the density correlation function is "contained" in the noise corre- lations and the term , NX j1;j2=1 hby1;j1b1;j2by1;j2b1;j1iei[k(j1 j2)+k0(j2 j1)]; is part of the full sum that needs to be taken for G11. This term can also be written as a function of the density operator, nj = byjbj , as NX j1;j2=1 hnj1nj2ie i(k0 k)(j1 j2) + (k k0)hnj1ie i(k0 k)j1 This shows that G11 and S(k) (Eq. 4.25) have the same Fourier transform of the density correlation function. If S(k) develops cusps at 2kF, where kF = [70], when CDW order is present, we expectG11 to have similar cusps at k = k0 2kF. In Fig. 4.5, we show one example of a quasi-supersolid (SS) state [39], where CDW order coexists with SF order. 87 The parameters are listed in (e) of Table 4.1. In the plot, the correlation function C11(q), an integration ofG11(k;k0) along the direction k = k0 (Eq. 4.2), is compared with the structure factor S(k) of the same state. In both functions, we can see cusps appearing at 2kF. 4.3.2 Trapped systems We now discuss how the different types of order are affected by the presence of a trapping potential. To simulate the effect of a trap, we add a harmonic potential, (j jc)2(n1;j+n2;j) to the Hamiltonian in Eq. 3.1, where j is the site index and jc is the index at the center of the system. We then use the TEBD method to calculate the ground state. We also increase the system size to 80 lattice sites, and choose the total number of particles and the trap frequency to ensure that the boundary effect is negligible. One interesting feature of a trapped system is that different orders can coexist in the trap. A well-known example is the MI plateau at the center of the trap surrounded by a SF at the edge [72]. For repulsive inter-species interaction, we find coexistence of a CFSF plateau with a SF at its edge and a MI plateau with PSF at the edges for attractive inter- species interactions [70]. Despite the potential complication of coexistence of orders, we find clear signals for the pairing correlations of the PSF phase and the anti-pairing correla- tions of the CFSF phase. In Fig. 4.6, we show the behavior ofG12(k;k0) in four different cases, where the orders at the center of the trap are SF, MI, PSF and CFSF respectively. We find that the general behavior of the noise correlation in a trap is very similar to its homogeneous counterpart. In Fig. 4.6 (c) and (d),G12 shows clearly the feature of pairing correlations in the PSF state and the anti-pairing correlations in the CFSF state. In addition, we see some minor features attributed to the coexisting orders. In the case of CFSF in a trapped system, we can see the "dip" along k = 0 and k0 = 0 because of the coexistence with the SF order. On the other hand, in the case of a MI in a trap, we can see pairing correlations as a result of the residual PSF state at the edges. This pairing signal is much smaller than when the whole system is in the PSF state. To show that the peaks along k = k0 and k = k0 inG12 are detectable in experiments, 88 00 ?5 0 5x 10?3 pi pi ?pi G12(k,k?) ?pi kaLk?aL (a) MI 00 ?0.1 0 0.1 ?pi pipi ?pi G12(k,k?) kaLk?aL (b) SF 00 ?0.1?0.05 00.05 0.1G12(k,k?) pipi ?pi ?pik?a L kaL (c) PSF 00 ?0.1 0 0.1 0.2G12(k,k?) ?pi pi ?pipik?aL kaL (d) CFSF Figure 4.6: Noise correlations in the trapped system. The system size is 80 sites and t=U = 0:02. In (a), the system is in the SF state. The particle number of each species is 30, the trap frequency 8 10 5U and U12=U = 0:01. In (b), the particle number of each species is 40, the trap frequency 1 10 4U andU12=U = 0:11. The system has both MI and PSF orders. The MI state forms a plateau at unit-filling at the center of the trap and the PSF is formed at the edge. The PSF state at the edge causes the small peak along the k = k0 direction, similar to the one in (c). However, this peak is at a much smaller amplitude than the one shown in (c), where the whole system is a PSF state. In (c), the particle number of each species is 20, the trap frequency 1 10 5U and U12=U = 0:11. The whole system is in the PSF state. A strong pairing correlation is formed along k = k0 direction. In (d), the particle number of each species is 30, the trap frequency is 8 10 5U andU12=U = 0:2. The system has both CFSF and SF order. The CFSF order forms a plateau at half-filling at the center of the trap and the SF state towards the edges of the trap. The CFSF order causes a strong anti-pairing (particle-hole) correlation along k = k0direction. At the same time, the SF order adds to the "dips" along k = 0 and k0 = 0. 89 Figure 4.7: Correlations C12(q) and D12(q) for the states that are described in Fig. 4.7. In (a), we show the behavior of C12(q) (Eq. 4.2) in SF, MI, PSF and CFSF states. The strong anti-pairing (particle-hole) correlations in the CFSF state gives a strong signal aroundq = 0 in C12(q). This strong signal is also unique to the CFSF state and therefore can be used to detect to the CFSF order. In (b), we show the behavior of D12(q) (Eq. 4.3) in SF, MI, PSF and CFSF states. The strong pairing correlation in the PSF state is the reason for the high peak around q = 0 in C12(q): This suggests measuring C12(q) is a good way of detecting the PSF order. 0.8 1 1.2 0.9 1.1 qaL C 11(q) (a) pi?pi 0 0 5 10 15 20 25 S(k) kaL (b) ?pi 0 pi Figure 4.8: Noise correlationC11(q) and structure factorS(k) in a PSF/CDW state. The sys- tem size is 80 sites and there are 20 particles of each species ( = 0:25). The trap frequency is = 10 5U , the hopping t = 0:02U and the inter-species interaction is U12 = 0:11U. The density at the center of the trap is roughly 0:45 per site and the cusps are developed around 0:9 . The inhomogeneity of a trapped system means that the ?Fermi wave vec- tor? kF is no longer , where is the average filling of the system. Instead, kF can be evaluated as ncenter, where ncenter is the density at the center of the trap, 90 we also calculate C12(q) (Eq. 4.2) and D12(q) (Eq. 4.3) for the four states. In the correlation C12(q) (Fig. 4.7 (a)), a high peak at q = 0 only appears in the case of the CFSF state. This peak corresponds to the peak in G12 along k = k0 in the CFSF state and is a reflection of the anti-pair correlation in the CFSF state. Similarly, in D12(q), the high peak at q = 0 only appears in the PSF state, as a result of the strong pairing correlations in the PSF state. A similar measurement has been performed for fermionic mixtures to detect the pairing of fermions [33]. In addition to the PSF and CFSF order, we also look for the signal of CDW order in G11(k;k0) in the trapped system. In a trapped system, the CDW order is more difficult to establish especially in the PSF and SF states, because the varying local density makes the "Fermi wave vector" n a spatially varying quantity. However, we can still see weakened cusps forming at the momentum roughly corresponding to 2 ncenter, where ncenter is the density at the center of the trap. This may indicate that in the trapped system, the CDW order in PSF and SF states has a wave vector corresponding to the density at the center of the trap. For the CFSF state, because the system has a plateau at half filling, the wave vector 2kF is =aL. Compared to the homogeneous case, this feature is slightly diminished due to the effect of the coexisting SF state in the trapped system. In Fig. 4.8, we show one case where CDW order coexists with PSF order in a trap. The system size is 80 sites and there are 20 particles of each species. The trap frequency = 10 5U , t = 0:02U and U12 = 0:11U. The density at the center of the trap is roughly 0:45 per site. The cusps are developed around 0:9 , which is roughly 2 ncenter. 4.3.3 Determination of Luttinger parameters from experimental data Another important question is whether we can use the noise correlationG12 to measure the Luttinger parameters,KS andKA, in the PSF and CFSF regimes. The LL calculation shows that as the system size approaches infinity, the noise correlationG12 approaches a power law decay with the power 1=KS in the PSF regime and with 1=KA in the CFSF regime (see Eqs. 4.20 and 4.21). In our numerical results for C12(q) and D12(q), we indeed find that the decay from the peak at q = 0 satisfies the algebraic decay. To find out the power 91 of the algebraic decay, we fit the function C12(q) in the PSF regime and D12(q) in the CFSF regime with the fitting function, F(q) = Ajsin(2q)j 1=K +B; (4.26) where B is the minimum value of C12(q) or D12(q) and A and K are the fitting parameters. In the PSF case (U12=U = 0:11), we find that K is 1:3 0:1. This is indeed very close to the value of KS, which is estimated at 1:4 0:1 obtained by the algebraic fit of RS. In the CFSF case (U12=U = 0:2), we find that K is roughly 1:48 0:1, while the value of KA extracted from the algebraic fit of RA is also at 1:48 0:12. Because of the singularity at q = 0, a reasonable values of K can be obtained by a simple algebraic decay function, Aq 1=K+B, around smallq. This shows that even in a trapped system, one can still assume a algebraical relationship predicted in the LL theory (Eqs. 4.20 and 4.21) and estimate the values of the Luttinger parameters by studying the power of the decay from the peak at q = 0. 4.4 Summary What?s most interesting in the behavior of the noise correlation is that it provides a complete set of signatures for all the phases that can possibly exist in the 1D binary Bose mixtures. These includes distinctive signatures for SF, MI, PSF, CFSF and CDW orders. In particu- lar, we show the direct measurement of pairing order and the anti-pairing order through noise correlations, which can not be detected through the averaged time-of-flight density distribution. These features still can be observed even for a harmonically trapped system and we show that the modifications induced by the inhomogeneity understood in terms of the results for the homogeneous system. Even with the presence of a trap, the noise correlations in the PSF/CFSF regime still have distinctive peak structure that are being re- lated with the Luttinger parameters, KS/ KA . This suggests that one can use the noise correlation to estimate the Luttinger parameters in these two regimes. 92 Chapter 5 One dimensional binary Bose mixtures in optical lattices (III): transport properties In [73], the transport property of an one-dimensional single-species Bose gas in a harmonic trap and a periodic potential is studied through the dipole oscillation of the center of mass of the atoms excited by suddenly displacing the total harmonic trap of the gas. Depending on the depth of the periodic potential, the dipole oscillation can change from undamped when the periodic potential is zero, to heavily damped when the periodic potential is rel- atively weak and the gas is in the strongly interacting superfluid (SF) state, and to over- damped when the periodic potential is strong. In [74, 75], the dipole oscillation of 1D Bose gas is studied numerically and have found good agreement with the experiment. Here, we consider the transport properties for one-dimensional binary Bose mixtures. As studied in the two previous chapters, there are two additional phases in the mixtures: the paired superfluid and the counter-flow superfluid. In addition to their distinctive sig- nals in the noise correlation (Chap. 4), we will show that the pairing (anti-pairing) mech- anism in PSF (CFSF) states leads to very different transport properties. For the PSF state, because the particles are bounded as pairs, two species tend to behave the same way un- der perturbations, even when the kinetic energy injected into the two species are different. This is reflected as the suppression of the relative motion and the enhancement of the total motion of the particles. The opposite case can be found for the CFSF state, where because of the anti-pairing correlation, different species tend to avoid each other and move in ex- actly the opposite directions. The total motion in this state is strongly suppressed and the relative motion much enhanced. It also important to notice that the pairing and anti-pairing ordering can exist in other systems as well. If we assume that each particle carries the same amount of charge and the two species are distinguished by up and down spins, the total (relative) motion discussed 93 here corresponds to the charge (spin) mobility of the system. If we consider the case that the two species are distinguished by their location (like the case atoms localed on different lattice planes interacting through dipolar interactions [60]), the total (relative) motion here corresponds to the correlated (anti-correlated) motion between planes. 5.1 Dynamics after a trap displacement 5.1.1 Model and parameters We consider two independently controllable traps for the two species and let a be the trapping frequency of species a and ja;c(t) are the location of the trap minimum, which we allow to be time dependent. This corresponds to a time-dependent term in the Hamilto- nian (Eq. 3.1), NX j=1 X a=1;2 a(j ja;c(t))2na;j: (5.1) Here, na;j is the number operators andj is the lattice site with indexj. For this calculation, we assume a system of 80 lattice sites and both traps are initially centered at the trap, ja;c(t = 0) = 40:5. The displacement is the location of the trap with regard to the center of the lattice, Da(t) = ja;c(t) ja;c(0): (5.2) Both Da(t) and ja;c(t) are written in units of the lattice site. The parameter J is used to replace the hopping parameter t in the Hamiltonian of Eq. 3.1 to avoid the confusion with time t. We consider four parameter sets as listed in Table 5.1. The sets correspond to the SF, PSF and CFSF states respectively. The hopping parameters are rather large in the PSF and CFSF region. This choice is made in order to exclude the influence of the CDW order that coexists with the PSF and CFSF order for smaller hoppings. The inter-species interactions are based on the phase diagram in Figs. 3.7 and 3.8. The particle number and the trap frequency of either species are chosen to make sure that the density at the center of the 94 SF PSF CFSF U12=J -0.1 -5.6 5.6 a=J 2 10 3 2 10 4 8 10 3 Na 20 20 21 Na(Fig. 5.3) 10 10 21 Table 5.1: Parameters used to represent the SF, PSF and CFSF states. The hopping param- eter is set at U=J = 8 for all cases. The particle number of each species is the same and denoted as Na. stands for the particle number of either species. trap can not be greater than 1 and the density is far from the edge of the lattice to minimize the boundary effects. The particle number and trap frequency chosen for the CFSF state also ensure that the density at the center can be higher than 0:5. We also consider one additional low density case for the PSF state in order to present the general feature of the PSF for all different fillings. 5.1.2 Simulation of the real time evolution The initial state for the time evolution is the ground state of the system generated by the imaginary time propagation of the TEBD method. Details about the method can be found in Chap. 2. The trap is displaced at t = 0 and the real-time propagation with the TEBD method is used. Even although the Hamiltonian is time dependent, it is a close system and we can assume that for the propagation step from t to t+ t, where t is sufficiently small, the Hamiltonian is time independent, j (t+ t)i= e i~H(t) tj (t)i: (5.3) This real-time propagator is applied on to the state in the same way as the imaginary-time propagation (see Chap. 2 for details), except that the time step now is real. Similar to the imaginary-time propagation calculation, we use the second order Trotter-expansion to decompose the propagator for a small time step t. Here, t is set as tJ=~ = 0:05, where ~ is the Plank constant and is set to 1 in the simulation. The error introduced the expansion is estimated as O ( tJ=~)3 10 4. The real time propagator is a unitary transformation and the total number and total energy are constant during the propagation. Because of 95 the Trotter-expansion, however, errors are introduced. We use the number conservation algorithm (explained in Chap. 2) in the propagation to ensure that the total number is conserved. The total energy is monitored during the propagation and we find it stays closely to a constant, with fluctuations on the order of the expansion error, throughout the whole propagation. In addition to the total energy, several other observables, such as the density and corre- lations, are calculated along the propagation. In this chapter, we focus on the behavior of the density distributionhna;j(t)ias a function of time. To capture the collective behavior of each species, we calculated the center of mass for each species, defined as xc;a(t) = X j jhna;j(t)i Na : (5.4) At t = 0, both centers of mass are at the center of the trap, xc;a(t = 0) = 40:5. After the displacement, the center of mass starts oscillating, which is a reflection of the collective response to the displacement of the corresponding species. 5.2 Result In Chap. 3, we have discussed our numerical calculation of the density distributions in the SF, PSF and CFSF states in a trap (see Fig. 3.9). We have shown that both PSF and SF are compressible, while the CFSF state is incompressible. For a trapped system, the PSF order exists for the whole system. The CFSF state, however, exists only near the center of trap where a plateau at half-filling is formed and a SF state coexists with the CFSF state at the edge of the trap. To probe the dynamical response of the PSF and CFSF state to an external perturbation, we consider two different type of the displacement. The first is a very brief displacement and the trap is brought back to the center. The second is a small displacement and the traps are placed in the new position for the rest of the time-evolution. Physically speaking, the brief displacement is like inject a small impulse and the constant displacement is like placing the system in a new potential configuration. 96 0 20 40 60 80 02040 6080100 0 0.2 0.4 0.6 (b) jt0 20 40 60 80 02040 6080100 0 0.2 0.4 0.6 jt 0 5002040 6080100 00.2 0.40.6 0.81 (c) jt 0 50020406080100 ?0.4?0.2 00.2 0.4 (d) jt (a) n1,j(t)+n2,j(t) n1,j(t)?n2,j(t) n1,j(t) n2,j(t) Figure 5.1: Dipole oscillation in the CFSF state. Note that time t is in the unit of ~=J for all plots. At t = 0, the density distributions of both species have the same plateau at half- filling and then the trap of species 1 is briefly perturbed by a displacement of one lattice site [see Fig. 5.2 (a)]. In (a) and (b), we show the density distributions of species 1 and 2 as a function of time after the brief displacement. Species 1 and 2 have exactly opposite oscillation. In (c), we show the time evolution of the total density of both species. There is no oscillation in the total density and the plateau at the center is fixed at one. In (d), we show the time evolution of the relative density between the two species. The relative density reflects the oscillatory motion in species 1 and species 2. 5.2.1 Brief displacement In this section, we would like to start the discussion with a compelling example for the counter-flow property in the CFSF phase (In Fig. 5.1). Here, only the trap of species 1 is displaced for a very brief time t (see Fig. 5.2(a)). The species 1 response to the displace- ment with oscillations that can be observed in its density distributionhn1;ji. Interestingly, atoms of species 2 instantaneously start to move in exactly the opposite way, even though their trap has not been displaced. What?s more interesting is that the distribution of the total density in the CFSF state stays unchanged throughout [Fig. 5.1 (c)] and the oscilla- tory motion is completely captured in the relative motion [Fig. 5.1 (d)]. This mechanism of moving with fixed total density can be understood intuitively from the anti-pairing (or particle-hole) order. The particle-hole pairing requires that the local density of individual 97 0 100 200 300 40039.5 40 40.5 41 41.5 t x C 0 100 200 300 40040.4 40.45 40.5 40.55 t x C 0 200 400100 300 40.539.5 38.537.5 41.542.5 43.5 t x C 0 10010 ?1 0 1 t Displacement species 1 species 2 (a) (b) (c) (d) Figure 5.2: (a) Trap displacement as a function of time. Only the trap of species 1 is dis- placed at t = 0. The displacement is one lattice site and is removed at t = 1J=~. (b)-(d) Dipole oscillation of the center of mass in SF, PSF and CFSF states induced by the displace- ment of species 1?s trap. In the SF state (b), only species 1 are excited at t = 0 as a result of displacement. In PSF (c) and CFSF (d) states, both species are excited at t = 0 as a result of their pairing and anti-pairing order respectively. The red(green) color stands for species 1(2) for all plots in this chapter. 98 species satisfynj;1 = 1 nj;2 , where 1 nj;2 represents the density of the hole of species 2 for at site j. This relationship does not restrict the behavior of the individual species, rather it puts a condition on the collective behavior of the two species. Although the density distribution gives an intuitive picture of the oscillation, the am- plitude of the oscillation can be too small to be visible in the density distribution. To study the oscillatory motion systemically, we consider the motion of the center of mass. In Fig. 5.2, we show the time-evolution of the center of mass for individual species in the SF, PSF and CFSF states. The trap displacement considered here is the same as the one in the pre- vious paragraph. In Fig. 5.2 (b) we show the response of a SF state. Here, we see the excitation of only species 1 in the beginning. As a result of weakly attraction between the two species interact, the oscillatory motion of species 2 is induced gradually in a later time. In Fig. 5.2 (c), we show the response of the PSF state. In this case, due to the pairing or- der, both species start moving instantaneously and the time evolution of the two centers of mass are identical. The energy injected into species 1 is transformed into the total motion of both species. In Fig. 5.2 (d), we consider the CFSF case, the density distribution of which is discussed in the previous section. From the center of mass motion, we can clearly see that the oscillatory motion of species 1 is perfectly matched by an counter-flowing motion of species 2 and the total center of mass xc;1(t) +xc;2(t) is kept almost constant. From Fig. 5.2 c) and d) it is apparent that a characteristic feature of the transport in the CFSF (PSF) phase is that the total (relative) density oscillations are strongly suppressed. If we adopt a pseudo-spin language in which species 1 corresponds to spin-up and species 2 to spin-down, we can say that in the PSF state, the perturbation of the system leads to excitations in the charge sector of the system and in the CFSF state, the perturbation leads to excitation in the spin sector. If we apply the opposite type of impact for the two phases, i.e. total motion excitation for the CFSF state and relative motion excitation for the PSF state, the motion of the atoms will be entirely inhibited. To illustrate such effects, we consider displacing the traps in the same and opposite directions. The same displacements only excite the total motion. In this case, we expect only the PSF state be excited. For the opposite displacement, all energy is injected in the relative motion. Therefore only the CFSF state is excited. This is indeed what we find in our simulation. In Fig. 5.3, we show 99 0 50 100 150 200 250 300 350 400 45040.2 40.3 40.4 40.5 40.6 40.7 t x C 0 50 100 150 200 250 300 350 400 45034.5 36.5 38.5 40.5 42.5 44.5 46.5 t x C 0 10 100 ?1 0 1 tDisplacement 0 10010?2 ?1 0 1 2 t Displacement (b) (c) (d) (a) Figure 5.3: (a) Trap displacement as a function of time. Here, two traps are displaced by 1 lattice site at t = 0 and brought back to the original place at t = J=~. (b) Dipole oscillation in the PSF state after the same displacement. (c) Trap displacement as a function of time. Here, two traps are displaced by 1 lattice site from the center and brought back to the center at t = J=~. (d) Dipole oscillation of center of mass in the CFSF state after the opposite displacement. The parameters used are shown in Table. 5.1. Note that there is no dipole oscillation for the CFSF state under the same displacement and for the PSF state under the opposite trap displacement. 100 0 100 200 300 400 500 600 700 800 90038.5 39 39.5 40 40.5 41 t x C 0 400 800 1200 160038.5 39 39.5 40 40.5 41 t x C 0 100 200 300 400 500 60038.5 3939.5 4040.5 41 t x C 0?2 ?1 0 1 2 tDisplacement (c) (a) (b) (d) Figure 5.4: (a) Trap displacement as a function of time. Both traps are displaced by 1 lattice site at t = 0. In (b), the system is at SF state. Both species show the heavily damped dipole oscillation. In (c), the system at the PSF state and the two species show identical heavily damped dipole oscillations. In (d), the system is in a CFSF state. The dipole oscillation is overdamped in this case and the center of mass stays near the original equilibrium position before the displacement. an example of the excitation of the total motion of the PSF state and the relative motion of the CFSF state under the same and opposite oscillation respectively. The inhibition of motion is also found in the constant displacement, which is discussed in the following section. 5.2.2 Constant displacement Besides a brief, temporary displacement of the trap to impart an impulse on the atoms, we also consider displacing the trap permanently. To demonstrate the CFSF and PSF orders, we consider simultaneous displacements in the same and in opposite directions. In Fig. 101 0 100 200 300 400 500 600 700 800 900 ?8 ?6 ?4 ?5 0 2 4 6 8 10 ?10 t x C 0 20 40 60 800 0.5 1 t=55 0 20 40 60 800 0.5 1 t=100 0 20 40 60 800 0.5 1 t=590 0?2 ?1 0 1 2 t displacement (c) (d) (e) (b) (a) PSF SFCFSF Figure 5.5: (a) Trap displacement as a function of time. The traps are displaced in the opposite direction by 1 lattice site at t = 0. (b) Dipole oscillation of the center of mass induced by the trap displacement. In this case, the PSF state prohibit the separation of the two species and the center of mass stays at the original equilibrium position of the trap before the displacement. For both CFSF and SF states, the center of mass of each species moves apart from each other under the displacement. In the SF case, the oscillation converges to the new equilibrium position of the traps after the displacement. But, in the CFSF case, each species keep moving away from each other, until they reach the edge of the CFSF state. In (c)-(e), we show the density distribution at different times for the CFSF state. 102 5.4 we consider the time-evolution under a constant displacement of the same amplitude. The response of three phases is markedly different. Both the SF and the PSF state exhibit oscillatory motion around the new minimum of the trap. However, the frequencies differ not only because of the different trapping frequencies but also because of the higher ef- fective mass of the pairs in the PSF state. The motion of the CFSF state is almost entirely inhibited. The atoms stay essentially at the original location and the transport of the total density is suppressed. The small oscillations that are visible in the time evolution of xc;1(t) and xc;2(t) are due to oscillations of the superfluid fraction at the edges. In Fig. 5.5 we consider displacing the traps in opposite directions. For the SF state, the two atomic species oscillate around the new trap minimum. The PSF state shows essen- tially no response, because the formed pairs are not broken up by the force imparted by the trap displacement. The most dramatic response is in the CFSF state. Here, the displace- ment of one lattice site leads to roughly ten-lattice-site relative displacement of the center of mass of each species. The small displacement in opposite directions is enhanced by the counter-propagating character of anti-pairing. To illustrate the response of the state fur- ther, we show density distributions of each of the two species as well as the total density. The total density distribution is again essentially fixed throughout the entire evolution. It is worth noting that the separation between the two species is limited by the size of the trap. It is reasonable to expect that the equilibrium position for the two species after the trap displacement will always be at the opposite edges of the CFSF plateau. If the one di- mensional lattice is connected like a ring, this may lead to persistent counter-propagating flow. 5.3 Summary The different responses of the PSF and CFSF states to a potential shift (Figs. 5.4 and 5.5), as well as to impulses (Fig. 5.2 and 5.3) vivid illustrate the defining features of the PSF and CFSF order, where the cross-species correlation leads to novel transport properties of the system. Under certain perturbations, the suppression of the individuality leads to similar transport properties as a Mott insulator state of a single-species system. However, for other 103 perturbations, the PSF and CFSF states manifests their superfluid transport properties. Because such a selective excitation pattern is purely the result of pairing and anti-pairing ordering, it is reasonable to expect that very similar transport properties can be found for all the systems where the pairing (or anti-pairing) order is dominant. We therefore consider the transport properties based on the total/relative motion be generically true for PSF/CFSF order and is not limited for 1D binary Bose mixtures. Because the trap displacement can be easily controlled in the experiment, the transport properties predicted here should be realizable in the experiments once the PSF/CFSF order is present in the system. By increasing the displacement amplitude in the opposite way for the PSF state, it could potentially break pair ordering and therefore recover the individual oscillation. The displacement amplitude then becomes an indicator of the binding energy of the pairs. Similarly, large displacements in the same direction can break the anti-pairs in the CFSF state and provide information about the energy spectrum in the CFSF state. 104 Chapter 6 Two dimensional Bose-Fermi mixtures in optical lattices: pre-formed molecules and novel thermometry 6.1 Introduction In Chaps. 3-5, we focus on discuss the new orders induced by the inter-species interaction in a mixture of bosonic atoms. Here, we would like to consider another different systems and the work discussed here is mostly motivated by the growing interests of ultra-cold polar molecules [76], as they have the promise for being a new state of quantum degenerate matter, with unique properties. In order to have a large dipole moment, the polar molecules must be in their rovibrational ground state, where further cooling can ultimately lead to quantum degenerate dipolar matter [77]. Such polar molecules can have long-range, anisotropic or three-body interactions [78], which may lead to novel quantum phases [58, 80] and new applications in quantum information science [81]. In most ultra-cold polar molecule experiments, one starts with a mixture of ultra-cold gases of atoms of different species, for example various isotopic combinations of K and Rb [85, 84, 86, 83, 82]. These atoms can form a weakly bound state through a magnetic field sweep over the Feshbach resonance [87, 86]. To create molecules with significantly higher dipolar moments, the loosely bound Feshbach molecules are coherently transferred to a ground state with very high efficiency through stimulated Raman adiabatic passage (STIRAP) [88, 89, 90, 91]. Although the rate of transferring a Feshbach molecule to the ground state is very high, the overall efficiency for forming dipolar molecules is still low due to the low efficiency of forming the loosely bound Feshbach molecules during the field sweep. In Ref. [90], the fermionic 40K and the bosonic87Rb atoms are trapped by an optical trap. The efficiency to form the Feshbach molecule depends on the phase-space density of the two species. But, 105 because the Fermi cloud stops shrinking once it reaches the quantum degenerate regime, and the Bose cloud continues to shrink as it Bose condenses, this phase space density is low at low temperature, and never reaches appreciable sizes at higher temperatures, as the clouds become more diffuse. On the other hand, if the mixture is first loaded onto an optical lattice, the motion of the atoms can be more strongly confined, and it is possible to create a large area where exactly one atom of each species sits at the same lattice site, leading to a reduced three body loss [87] and almost unit efficiency [92] for pre-forming the molecules. When mixtures of 40K and 87Rb are loaded into an optical lattice, the atoms of each species are influenced by the optical lattice differently[?]. With the same optical lattice depth, the heavy atoms usually have much lower tunneling rate than the light atoms be- cause of their significantly larger mass. In Ref. [92], it was shown that for sufficient lat- tice depths, the hopping rate of Rb is more than an order of magnitude less than that of K. It is therefore reasonable to ignore the quantum effects of the tunneling of the heavy bosonic atoms while allowing the light fermionic atoms to hop between nearest neighbors (a classical effect of the motion of the Rb atoms is taken into account by averaging over all energetically favorable distributions of Rb atoms). Such systems can be described by the Fermi-Bose Falicov-Kimball model [93, 94, 95]. Using this model, we quantitatively deter- mine the probability of having exactly one atom of each species per lattice site in order to optimize the formation of dipolar molecules. For the Falicov-Kimball model, the phenomena of pre-forming molecules has been dis- cussed for Fermi-Fermi mixtures or Fermi-hard-core-Bose mixtures [96] on a homogeneous lattice and Fermi-Fermi mixtures in a harmonic trap [97]. In previous work [92], we consid- ered the Fermi-soft-core-Bose mixtures in a harmonic trap and determined the efficiency for pre-forming molecules as the probability to have exactly one atom of each species per site. We used inhomogeneous dynamical mean-field theory (IDMFT) and Monte Carlo (MC) techniques to calculate the efficiency as well as the density profile and the entropy per particle. Both of these methods have advantages and disadvantages. The IDMFT approach is approximate for two-dimensional systems, but it can calculate both the effi- ciency and the entropy per particle. The MC method is numerically exact after it reaches 106 thermal equilibrium, but it can not calculate the contributions to entropy coming from the heavy particles. Both methods require large computational times to calculate properties of a trapped system of reasonable size. Using these methods, we have shown that the effi- ciency is significantly increased by first loading onto an optical lattice before forming the molecules and near unit efficiency can be achieved with parameters that are realistic for current experiments. The efficiency of pre-formed molecules is also likely to be affected by the heating (the temperature increase) induced by loading onto an optical lattice [1, 99, 98]. Considering that thermal fluctuations generally destroy the ordering and the localization of the parti- cles, it is reasonable to expect that the efficiency of having exactly one Rb atom and one K atom per site should be reduced if the temperature becomes too high. On the other hand, if the temperature is low enough, the presence of the lattice significantly increases the ef- ficiency, almost to unity in the case of deep lattices. The temperature of the lattice system, however, remains difficult to measure in experiment [104, 100, 102, 103, 101]. Instead, it is often assumed that the process of loading atoms onto optical lattices is adiabatic and therefore the total entropy of the system is conserved [98, 106, 1, 99, 105]. Determined based on the thermal properties of the gas before adding lattices, the entropy per particle is then used as an effective temperature scale for the lattice system [107, 1]. There are also several proposals for directly determining the temperature for systems of bosonic atoms [109, 108, 110], fermionic atoms [111] or the magnetic systems [112]. In Ref. [104], a general thermometry scheme is derived based on the fluctuation-dissipation theorem. Through quantum MC simulation, this proposal is shown to be applicable to the non-interacting fermionic systems [104] and interacting bosonic systems [113]. In our paper, we discuss light-Fermi-heavy-Bose mixtures in optical lattices based on the strong-coupling (SC) expansion method (perturbation theory in the hopping). The cal- culation is oriented to develop an efficient way of estimating the efficiency of pre-forming molecules for a given experimental system. With the SC expansion method, we obtain analytical expressions for the efficiency of pre-forming molecules, the entropy per parti- cle and the local charge compressibilities. The behavior of the efficiency is studied both as a function of entropy per particle and temperature. The determination of temperature 107 is further studied by applying the thermometry proposal in Ref. [104] to the Fermi-Bose mixture. To benchmark the accuracy, we compare the SC calculation with the IDMFT and MC calculations for all parameters considered. Overall, we find excellent agreement be- tween the three methods. Such agreement even extends to the low temperature region when the interaction is strong enough. This is particularly useful, given the fact that the SC expansion calculation is significantly faster than the IDMFT and MC calculations. Such a speedup makes it possible to consider much larger lattice sizes to eliminate the bound- ary effects, to scan the large parameter space for optimal parameter regions for pre-formed molecules and to estimate the density fluctuations and other properties. 6.2 Fermi-Bose Falicov-Kimball model For mixtures of heavy bosons and light fermions, such as 87Rb/40K mixtures, the hopping parameter for the heavy bosons (87Rb) is usually more than an order of magnitude less than the hopping parameter for the light fermions (40K) when one takes reasonable lattice depths (greater than 15 Rb recoil energies) [92]. In this case, we can ignore the quantum- mechanical effects of the hopping of the heavy bosons and describe such mixtures with the Fermi-Bose Falicov Kimball model in the presence of a trap potential. The Hamiltonian of this model is written as H = H0 +Hh = X j H0j +Hh; (6.1) with H0j = (Vj f)fyjfj +Ubffyjfjbyjbj + (Vj b)byjbj + 12Ubbbyjbj(byjbj 1); (6.2) and Hh = X jj0 tjj0fyjfj0: (6.3) 108 Here, j, j0 label the sites of a two-dimensional square lattice, with a lattice constant, a. The symbols fyj and fj denote the creation and annihilation operators for the fermions at lattice site j, respectively. The symbols byj and bj denote the creation and annihilation operators for the bosons at lattice site j, respectively. The fermionic operators satisfy the canonical anticommutation relation ffj;fyj0g = j;j0 and the bosonic operators satisfy the canonical commutation relation [bj;byj0] = j;j0. The quantity Vj is the trap potential, which is assumed to be a simple harmonic-oscillator potential centered at the center of the finite lattice. We assume that the jth site has the coordinate (xj;yj), so that Vj can be written as Vj = t ~ 2ta 2 x2j +y2j ; (6.4) where is the trap frequency. The quantity f is the chemical potential for fermions and b is the chemical potential for bosons. Combining the trap potential and the chemical potentials, we can define an effective position dependent local chemical potential for the fermions, f;j = f Vj, and for the bosons, b;j = b Vj. Ubf is the interaction en- ergy between fermions and bosons and Ubb is the interaction energy between the soft-core bosons. The symbol tjj0 is the hopping energy for fermions to hop from site j0 to site j. We consider a general tjj0 for the formal developments in the earlier part of the next sec- tion, but later specialize to the case of nearest-neighbor hopping with amplitude t, which we will take to be the energy unit. We also set the lattice constant, a equal to one. The efficiency for pre-forming molecules is defined as the averaged joint probability of having exactly one boson and exactly one fermion on a lattice site, E = 1N X j h^Pj1;1i= 1N X j Tr ^ Pj1;1e H ; (6.5) with = (kBT) 1 the inverse temperature. We define the projection operator ^Pj1;1 for having exactly one boson and one fermion at site j; ^Pj1;1 =jnb;j = 1;nf;j = 1ihnb;j = 1;nf;j = 1j; (6.6) 109 and N is the smaller value in the total numbers of bosons and fermions, Nb and Nf. In our case, we assume equal number of bosons and fermions, therefore N = Nb = Nf. In general, one can obtainE directly from Eq. (6.5) for a readily diagonalized Hamilto- nian. In our case, the efficiencyE is derived by distinguishing the contribution from terms corresponding to nb;j = 1 in the expression for the density of fermions. We assume that the density of bosons and fermions at site j, b;j and f;j, can be both be written as a series in terms of the occupation number of bosons at site j in the following way, b;j =hbyjbji= X nb;j Wj(nb;j)nb;j; (6.7) and f;j =hfyjfji= X nb;j Wj(nb;j)enf;j(nb;j): (6.8) Here nb;j is the occupation number of bosons on site j, nb;j = 0;1;:::. The coefficient Wj(nb;j) is the probability of having exactlynb;j bosons at sitej and the coefficientenf;j(nb;j) is the probability for having one fermion on sitej for the occupation numbernb;j. The joint probability of having exactly one boson and one fermion at site j can be written as Ej =Wj(nb;j = 1)enf;j(nb;j = 1); (6.9) and the efficiencyE is the average ofEj over all sites, E = P jEj N = P jWj(nb;j = 1)enf;j(nb;j = 1) N : (6.10) It can be shown that the expression for the efficiency obtained in this way is the same as from Eq. (6.5). Now, the efficiency is obtained directly from the density of bosons and fermions, which can be easily derived from the partition functionZ by taking derivatives 110 with respect to the appropriate chemical potentials, b;j = 1 @ln(Z)@ b;j ; (6.11) and f;j = 1 @ln(Z)@ f;j : (6.12) To study the behavior of the efficiency as a function of the entropy per particle, we evaluate the entropy per particle by dividing the total entropy by the total number of particles, s = S=(Nb +Nf) = kBln(Z) kB@ln(Z)@ =(Nb +Nf): (6.13) It is worthwhile noticing that the formalism development in this section is based on the grand-canonical ensemble. This ensemble is appropriate because we assume that in the lattice system both the energy and the number of particles fluctuate. This may seem in contradiction with the use of the entropy per particle as an effective temperature scale, because strictly speaking entropy is used as a parameter only for the micro-canonical en- semble. This contradiction is resolved because the entropy per particle is assumed to be conserved during the process of turning on the optical lattice. It is a conserved quantity when comparing the systems before and after turning on the optical lattice, which is par- ticularly useful from the experimental point of view, since the experiments often start with- out the optical lattices. For the lattice system itself, assuming it is in thermal equilibrium, it is more reasonable to consider it with the grand-canonical ensemble, allowing the energy and number fluctuations. The difference between the different ensembles of course is not problematic if we assume the system is large enough to be in the thermodynamical limit, where all three ensembles are equivalent. 111 6.3 Strong coupling expansion formalism In this section, we explain the SC expansion formalism. We first discuss the evaluation of the partition functionZ, approximated by the second-order expansion in terms of the hopping, Hh. From the partition function, we derive the expressions for the density of fermions shown in Eqs. (6.32) to (6.34), the density of bosons in Eqs. (6.35) to (6.37), the efficiency in Eqs. (6.38) to (6.40) and the entropy per particle in Eqs. (6.44) to (6.47). For readers who prefer to see the final expressions, we suggest skipping the following deriva- tion and referring to the equations listed above for the corresponding quantities. The evaluation of the partition function in the SC approach starts with the exact so- lution of the atomic Hamiltonian H0. Hence, we use an interaction picture with respect to H0, where for any operator A, we define the (imaginary) time-dependent operator A( ) = e H0Ae H0. The partition function is written using the standard relation, Z = Tr e H = Tr e H0U( ;0) : (6.14) Here, U( ;0) = T exp hR 0 Hh( )d i is the evolution operator with T being the time- ordering operator for imaginary times. Expanding the exponential inU( ;0) up to second order in Hh( ) and evaluating the resulting traces with respect to equilibrium ensembles of H0, we have U( ;0) ? 1 + + 12 Z 0 d 1 Z 0 d 2T Hh( 1)Hh( 2): (6.15) Here, we note that the first order correction to the partition function vanishes because the hopping connects different sites. Substituting Eq. (6.15) into Eq. (6.14), we can write the partition function as, Z =Z(0)(1 +Z(2)); (6.16) 112 whereZ(0) is the partition function in the atomic limit (t = 0), Z(0) = Tr e H0 ; (6.17) andZ(2) corresponds to the second-order term in the expansion ofU divided byZ(0), Z(2) = 12Z(0) Tr e H0 Z 0 Z 0 d 1d 2T Hh( 1)Hh( 2) : (6.18) To simplify the notation, we introduce f;j(nb;j) to represent the negative of the fermionic part of the Hamiltonian H0j [Eq. (6.2)] when there is a fermion at site j, f;j(nb;j) f Vj Ubfnb;j; (6.19) and b;j(nb;j) for the negative of the bosonic part of the Hamiltonian H0j, b;j(nb;j) ( b Vj)nb;j Ubbnb;j(nb;j 1)=2: (6.20) Both f;j and b;j depend on the number of bosons at site j. The effective fugacities for bosonic and fermionic particles can then be written as the exponential of f;j and b;j respectively, f;j(nb) = exp [ f;j(nb;j)]; (6.21) and b;j(nb) = exp [ b;j(nb;j)]: (6.22) The atomic partition functionZ(0) can then be written in terms of the effective fugacities as, Z(0) = jZ(0)j ; (6.23) 113 whereZ(0)j is the atomic partition function at site j, Z(0)j = X nb;j b;j(nb;j)(1 + f;j(nb;j)): (6.24) Now we evaluate the second term in the partition function,Z(2) of Eq. (6.31). To satisfy the total number conservation, only terms with j = k0 and j0 = k in Hh( 1)Hh( 2) are non- zero after the trace andZ(2) is reduced into a sum of products of the fermionic annihilation and creation operators at the same site, Z(2) = 12 Z 0 Z 0 d 1d 2 X jk tjktkj Tr h T e H0jfyj( 1)fj( 2) i =Z(0)j Tr h T e H0kfk( 1)fyk( 2) i =Z(0)k : (6.25) Using the cyclic permutation relationship of the trace, the products can be represented by the local atomic Green?s function, Gjj( ) = Tr h T e H0jfj( )fyj(0) i =Z(0)j ; (6.26) andZ(2) is expressed as integrations of the atomic Green?s functions in terms of their rela- tive times, Z(2) = 12 Z 0 Z 0 d 1d 2 X jk tjktkjGkk( 1 2)Gjj( 2 1): (6.27) Solving the Heisenberg equation of motion for the annihilation operator fj( ), @fj( ) @ = e H0 [H0;fj]e H0 (6.28) one easily finds the expression for the annihilation operatorfj( ) in the interaction picture, fj( ) = e f;j(nb;j) fj(0); (6.29) 114 Substituting Eq. (6.29) into Eq. (6.26), we obtain the atomic Green?s function in terms of the effective fugacities as, Gjj( ) = 8 >>< >>: Pnb b;j(nb)Z(0) j e f;j; > 0 P nb b;j(nb) f;j(nb) Z(0)j e f;j; < 0 (6.30) We now perform the integration over 1 and 2 inZ(2) and obtain the final expression forZ(2), Z(2) = 12 X jk tjktkj X nb;jnb;k b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k ( f;j(nb;j) f;k(nb;k)) f;j(nb;j) f;k(nb;j) : (6.31) Note that the partition function we derived here is not limited to the case of nearest- neighbor hopping with a uniform hopping parameter t. Eq. (6.31) can be applied to de- scribe hopping between arbitrary sites j and k and the hopping parameter tjk can vary over different sites of the lattice. Observables are evaluated by taking appropriate derivatives of the partition function. In calculating the derivatives, we truncate all final expressions to include only terms through the order of t2jk. Also note that because sites j and k are different sites, we do not normally have denominators equal to zero in Eq. (6.31), but in any case, the formulas are always finite as can be verified by l?H?pital?s rule. During numerical calculations of the observ- ables, the denominator, f;j(nb;j) f;k(nb;j), can become too small and cause numerical errors. In our calculations, we use the Taylor expansion in terms of f;j(nb;j) f;k(nb;j) around zero when the absolute value of f;j(nb;j) f;k(nb;j) is less than 10 5. The density distribution is evaluated by taking the derivative of the partition function with respect to the appropriate local chemical potential [Eqs. (6.12) and (6.11)]. For the den- sity of fermions at sitej, the expression constitutes two terms corresponding to derivatives fromZ(0) andZ(2), f;j = 1 @In(Z)@ f;j = (0)f;j + (2)f;j; (6.32) 115 where (0)f;j is the density of fermions in the atomic limit, (0)f;j = 1 @In Z(0)j @ f = P nb;j f;j(nb;j) b;j(nb;j) Z(0)j ; (6.33) and (2)f;j is the total contribution to the density at site j from particles hopping from all possible sites, (2)f;j = X k tjktkj X nb;jnb;k ( b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k 2 4 (1 (0) j;k) f;j(nb;j) + (0) j;k f;k(nb;k) f;j(nb;j) f;k(nb;j) + f;k(nb;j) f;j(nb;k)[ f;j(nb;j) f;k(nb;j)]2 #) : (6.34) Similarly, the density of bosons at site j is written as a sum of the atomic density and the hopping contribution as, b;j = 1 @ln(Z)@ b;j = (0)b;j + (2)b;j; (6.35) where (0)b;j = 1 @In Z(0)j @ b;j = P nb;j nb;j b;j(nb;j) [1 + f;j(nb;j)] Z(0)j ; (6.36) and (2)b;j = X k tjktkj X nb;jnb;k " (nb;j (0)b;j) b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k [ f;j(nb;j) f;k(nb;k)] f;j(nb;j) f;k(nb;j) : (6.37) The expression for the efficiency is obtained from the density distributions of the fermions and bosons. Similar to the expression for the densities, the efficiency consists of two terms, one corresponding to the atomic limit and one corresponding to the contributions from the 116 hopping, Ej =E(0)j +E(2)j ; (6.38) where E(0)j = b;j(nb;j) f;j(nb;j)jnb;j=1 Z(0)j ; (6.39) and E(2)j = X k X nb;j;nb;k b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k 8 >< >: b;j(n0b;j) f;j(n0b;j) Z(0)j ! n0b;j=1 [ f;j(nb;j) f;k(nb;k)] f;j(nb;j) f;k(nb;j) + nb;j;1 f;j(nb;j) f;j(nb;j) f;k(nb;j) + f;k(nb;j) f;j(nb;k)[ f;j(nb;j) f;k(nb;j)]2 #) : (6.40) For the trapped system, the local chemical potential j includes both the global chem- ical potential and the trapping potential Vj. The derivatives with regard to the local chemical potential or the chemical potential leads to different physical quantities. For the Fermi-Bose mixture considered here, the cross-derivatives should also be evaluated. Specifically, the derivative with regard to the global chemical potential ( b + f) corre- sponds to the total number fluctuations, = @ 2 lnZ @2( b + f) = h ^ Nf + ^Nb 2 i h^Nf + ^Nbi2: (6.41) Here we define the total number operators, ^Nf = Pjfyjfj and ^Nb = Pjbyjbj. The global compressibility is introduced as the response of the local density to the change of the global 117 chemical potentials, gj = @ 2 lnZ @( f;j + b;j)@( b + f) = h h fyjfj +byjbj ^ Nf + ^Nb i hfyjfj +byjbjih ^Nf + ^Nbi i : (6.42) And the local compressibility, or the onsite number fluctuation, is determined from the derivatives with regard to the local chemical potential, lj = @ 2 lnZ @2( b;j + f;j) = h fyjfj +byjbj 2 i hfyjfj +byjbji2 : (6.43) Both the global and local compressibilities are derivatives of the density distributions and can be obtained from the density expressions above. Finally, we obtain the expression for the entropy per particle defined in Eq. (6.13) by averaging the total entropy of the system and we again write the entropy per particle in terms of the atomic limit expression and the contributions from the hopping, s = 1N X j S(0)j + 1N X j S(2)j : (6.44) Here S(0)j is the entropy at site j in the atomic limit, S(0)j =kB = ln Z(0)j j; (6.45) where the parameter j corresponds to the onsite energy at site j in the atomic limit, j = @ln(Z (0)) @ = 1 Z(0)j X nb;j f b;j(nb;j) b;j(nb;j) [1 + f;j(nb;j)] 118 + f;j(nb;j) b;j(nb;j) f;j(nb;j)g: (6.46) The averaged contributions from the hopping at site j is S(2)j ; S(2)j =kB = ln(1 +Z(2)) @ln(1 +Z (2)) @ = 2 2 X k X nb;jnb;k b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k [ f;j(nb;j) f;k(nb;j)] f[ f;j(nb;j) f;k(nb;k)] [ b;j(nb;j) + b;k(nb;k) j k] + f;j(nb;j) f;j(nb;j) f;k(nb;k) f;k(nb;j)g: (6.47) This ends the discussion on the derivation of the SC expansion method formulas. In general, the expressions obtained above are accurate in the case when the hopping is much smaller than interaction strength and the temperature is very high ( t is small). In this pa- rameter region, the SC method can evaluate physical quantities, like the density distribu- tion, efficiency, compressibility and entropy, very efficiently. The total number of particles is fixed by varying the chemical potentials, b and f. To maximize the efficiency and reduce three body loss, we consider the low density region with attractive interspecies interactions and repulsive bosonic interactions. For other strong-coupling regions, the for- mulas developed above are equally applicable but not further discussed in this paper. 6.4 Results 6.4.1 Comparison with the IDMFT and MC calculations For a perturbative method like the SC expansion method, it is always necessary to de- termine the parameter regions where the approximation is valid. Here, we use the pre- vious results obtained from IDMFT and MC methods [92] as a reference to determine the accuracy of the SC calculation. It is also worthwhile to notice that the three meth- ods require substantially different computational times. The SC calculation usually takes 119 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (a) U bb=5.7t,Ubf=?2t SC MC IDMFT 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (b) Ubb=5.7t,Ubf=?6t 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (c) Ubb=5.7t,Ubf=?10t 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (d) U bb=11.5t,Ubf=?8t 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (e) U bb=11.5t,Ubf=?12t 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (f) Ubb=11.5t,Ubf=?16t Figure 6.1: (Color on-line) EfficiencyE as a function of temperature calculated by the SC (red cross), IDMFT (blue triangle) and MC (green circle) methods. The interaction param- eters, Ubb and Ubf, are shown in each plot. In (a), the SC calculation differs from the other two methods for T < 1t=kB. For this region, the SC expansion formulas derived here are no longer accurate. In (b)-(f), all three methods give almost identical results. These cal- culations also show that almost 100% efficiency is reached for relatively strong attraction, Ubf 6t, at low temperature, T 1t=kB, the SC calculations agree nicely with the other methods even for a relatively weak cross-species interaction, Ubf = 2t. The difference between the SC calculation and the other two methods can be under- stood from the fact that the SC method is a perturbative method based on the atomic limit of the Hamiltonian, t = 0 and that the properties derived from the SC expansion are dom- inated by the atomic-limit behavior with relatively small corrections from the hopping. In the atomic limit, bosons and fermions are completely localized and the only density fluc- tuations are due to thermal fluctuations. For the low density case considered here, the bosons always form a plateau of unit filling at the center of trap at low temperature and the fermions are attracted by the bosons one by one and form an almost identical plateau. The efficiency therefore always converges to unity as temperature deceases. In Fig. 6.1, we indeed find the efficiency from the SC calculation always goes to one at low tempera- tures. The convergence to unity is also true for the IDMFT and MC calculations for all the parameters except for Ubf=t = 2 and Ubb=t = 5:7. That?s where the SC calculation differs from the IDMFT and MC calculation. It is reasonable to assume that the SC calculation can be applied to the region where the ground state of the system is a localized, Mott insulator like state. 121 10?1 100 1010 1 2 3 kBT/t entropy per particle (k B) Ubb=11.5t,Ubf=?16t SCIDMFT Figure 6.2: (Color on-line) (a) Entropy per particle as a function of temperature T. The SC calculation is marked with red crosses and the IDMFT calculation by the blue line. We find excellent agreement between the SC calculation and the IDMFT calculation. 0 0.5 1 1.5 2 2.50 0.2 0.4 0.6 0.8 1 Entropy per particle (kB) Efficiency (b) Ubb=5.7t Ubf=?10t (SC) Ubb=?6t (SC) Ubf=?10t (IDMFT) Ubf=?6t (IDFMT) 0 0.5 1 1.5 2 2.50 0.2 0.4 0.6 0.8 1 Entropy per particle (kB) Efficiency (a) Ubb=11.5t Ubf=?16t (SC) Ubf=?8t (SC) Ubf=?16t (IDFMT) Ubf=?8t (IDMFT) Figure 6.3: (Color on-line) Efficiency as a function of entropy per particle for different interaction parameters. Note here that we didn?t include the case of Ubf = 2t, because it is already shown in Fig. 6.1 that the SC calculation is not accurate for low temperatures in this case. In (a) and (b), we consider two different bosonic interaction strengths and five different inter-species interaction strengths. For all parameters, the efficiency reaches 100% when the entropy is very low. For an intermediate entropy, with an entropy per particle around 1kB, the efficiency is around 80%. 122 The SC calculation of the entropy per particle is also compared with the IDMFT and MC calculations for all the parameters using Eqs. (6.44-6.47). The conclusion of the comparison is similar with the efficiency calculation, that the SC calculation is accurate except forUbf = 2t. In Fig. 6.2 we use one example, Ubb = 11:5t and Ubf = 16t, to represent all the cases where the SC calculation agrees with the IDMFT calculation. As the temperature increases, the entropy per particle starts to saturate at around 2:3kB. In the next section, we will show that this saturation is actually the result of finite-size effects. In Fig. 6.3, we show the behavior of the efficiency as a function of the entropy per particle. This figure can be compared with Fig. 2 in Ref. [92], where the IDMFT calculation is discussed. We verify the findings from the previous work that for strongly attractive inter-species interactions, an efficiency of 100% can be achieved at low temperature (low entropy) region. For an entropy per particle around 1kB, a 80% efficiency can still be reached. This efficiency is much higher than what has been achieved in experiment [90]. In the following discussion on the SC calculation result, we no longer consider the case of Ubf = 2t. This is also based on the consideration that the interaction of Ubf = 2t is too weak to achieve the desired high efficiency of pre-formed molecules and therefore is not in the parameter region of the main interest in this paper. 6.4.2 Finite-size effects In our calculations, we always assume a hard-wall boundary condition at the edge of the lattice. In experiments, however, the atoms are confined only by the trapping potential. This additional confinement imposed by the boundary condition can potentially affect the accuracy of our calculation. This finite-size effect can be neglected if the system is so large that the atoms trapped by the trapping potential almost never reach the edge the system. This, however, is not always the case for the 50 50 lattice. This problem is difficult to address with the IDMFT and MC methods, because of the high computational costs. The SC method, on the other hand, can calculate much larger systems for a fraction of the cost. In this section, we discuss our calculation for different lattice sizes and discuss finite size effects for different lattice sizes. To benchmark the SC calculations, the trap frequency 123 0 5 10 15 20 25 300 0.2 0.4 0.6 0.8 1 radial distance ? f(r) (a) T=1 t/kB 0 25 50 75 1001251500 0.2 0.4 0.6 0.8 1 radial distance ? f(r) (b) T=20 t/kB N=50 N=100N=200N=300 Figure 6.4: (Color on-line) Finite-size effects on the radial density profile. We assume a two- dimensional N N square lattice with hard-wall boundary conditions. The dotted lines indicate the boundaries of different lattices.The interaction parameters are Ubf = 16t, Ubb = 11:5t. We use the density distribution of the fermionic particle to represent the general dependence of density on the lattice size. In (a), we consider the case of low tem- perature, T = t=kB. Here, the density distribution is concentrated at the center of the trap and there is no difference between different lattice sizes. In (b), we consider the case of high temperature at T = 20t=kB. Here, the density is confined mainly by the size of the lattice. For N = 50, the density is confined at the edge of the lattice, r = 25. For N = 100, the density is again confined at the edge, r = 50. For both N = 200 and 300, the density goes to zero before reaching the edge of the lattice and the two distributions overlap with each other. We estimate that finite-size effects are eliminated for the 300 300 square lattice for the trap frequency and number of particles considered here. 124 10?1 100 1010 0.2 0.4 0.6 0.8 1 kBT/t Efficiency (b) 10?1 100 1010 1 2 3 4 5 kBT/t Entropy per particle (k B) (a) N=50 N=100N=200N=300 Figure 6.5: (Color on-line) Finite-size effects on the entropy per particle and the efficiency. We assume a two-dimensional N N square lattice with hard-wall boundary conditions. The interaction parameters are Ubf = 16t, Ubb = 11:5t. In (a), we show the behavior of the entropy per particle as a function of temperature for different system sizes. We see the entropy is significantly affected by the finite size when the lattice is smaller than around 200 200. The finite-size effect is not noticeable at lower temperature (T < 1t=kB). 125 0 1 2 3 4 50 0.2 0.4 0.6 0.8 1 Entropy per particle (kB) Efficiency (a) Ubb=5.7t Ubf=?10t Ubf=?6t 0 1 2 3 4 50 0.2 0.4 0.6 0.8 1 Entropy per particle (kB) Efficiency (b) Ubb=11.5t Ubf=?8t Ubf=?12t Ubf=?16t Figure 6.6: (Color on-line) Efficiency as a function of entropy per particle for a 300 300 square lattice system. We consider 625 atoms for each species. Compared with Fig. 6.3, the efficiency is significantly higher for the same value of the entropy per particle in the 300 300 lattice system when the entropy per particle is large. On the other hand, the behavior is similar in both lattice systems when the entropy per particle is less than 1kB. The unit efficiency is reached roughly when the entropy per particle is less than 0.5kB: 126 and the total number of particles are fixed for all different lattice sizes. We assume the largest lattice sizes are sufficient to neglect the finite-size effects. In Fig. 6.4, we show the density profile as a function of the lattice size at two temperatures, T = t=kB (a) and T = 20t=kB (b). Here Fig. 6.4 (a) represents the scaling behavior in the low temperature region, where there is no significant difference between different lattice sizes and Fig. 6.4 (b) represents the scaling behavior in the high-temperature region, where the system of small lattice size is highly affected by the boundary effect. Note that the horizontal axes are different scales in the two panels. The parameters used in the plots are Ubf = 16t and Ubb = 11:5t. We find similar behavior of the density profile for all the other parameters. In Fig. 6.5 (a), we show entropy per particle as a function of temperature at different lattice sizes. In this plot, we find that for small lattices, the entropy per particle becomes saturated at high temperature, while for large lattices it keep increasing as the temperature increases. The saturation is understood as the result of the finite-size effects. When the temperature is high, atoms tend to expand to a larger area in the trap, which leads to a large cloud size and higher entropy. When atoms expand to the edge of the lattice, the possible occupied sites are now constrained and the entropy stays similar even though the temperature increases, hence the saturation. When the lattice is sufficiently large, atoms can freely expand as the temperature increases and the entropy keeps increasing. The confinement of the atomic cloud in high temperature also affects the efficiency calculation. In Fig. 6.5 (b), we find that the efficiency saturates to a higher value for smaller lattices. This is because the confinement increases the density overlap between the two species. In the low temperature region, the atoms are close to unit filling at the center of the trap and the efficiency is similar for all difference lattice sizes. We find that a lattice of 300 300 sites is sufficient to eliminate the finite size effects for our parameter regions. Hence, we use this lattice size for the efficiency and entropy per particle calculations. In fig. 6.6, we show the result for the efficiency as a function of the entropy per particle. We estimate the calculation result from the 50 50 lattice is accurate when the temperature is around or below T = 1:25t=kB. 127 6.5 Thermometry 6.5.1 Temperature and Density fluctuations Based on the fluctuation-dissipation theorem, the compressibility can be related to the den- sity fluctuations as [104], = @ (r)@ = 1k BT h (r)Ni h (r)ihNi; (6.48) where (r) is the radial density profile, is the chemical potential and N is the total num- ber of particles. For a system with a spherically symmetric harmonic trapping potential, Vtr2, the local chemical potential at a radial distancer is Vtr2. Within the local density approximation, the trapping potential is interpreted as a variance in the chemical potential and the compressibility in the trapped system can be re-written as a function of the density gradient, @ (r) @ = 1 2Vtr @ (r) @r : (6.49) These two equations lead to a simple relationship between the density gradient and the density fluctuations in the trapped system, 12V tr @ (r) @r = 1 kBT (h (r)Ni h (r)ihNi): (6.50) Based on this relationship, one can determine the temperature from the independently measured density gradient and density fluctuations. For a two dimensional system, a simplified relationship can be found by integrating the above equation over all the two dimensional plane, Vt (0) = 1 kBT hN2i hNi2 : (6.51) Here, (0) stands for the density at the center of the trap. With the development of in situ measurements, it is now possible to measure the den- 128 0 10 20 30 40 50?0.2 0 0.2 0.4 0.6 (c) radial distance density fluctuation 0 10 20 30 40 50?0.2 00.2 0.40.6 0.81 (b) radial distance density fluctuation 0 10 20 30 40 50?1 ?0.5 0 0.5 1 1.5 (a) radial distance density fluctuation 0 10 20 30 40 50?0.2 0 0.2 0.4 0.6 (d) radial distance density fluctuation (kBT)?1(?f+?b) ?(2Vtr)?1?(?f+?b)/?r Figure 6.7: (Color on-line) Density fluctuations averaged over different numbers of sam- ples. The density fluctuations shown here are the total density fluctuations divided by the input temperature, T = 2tk 1B . All the fluctuations are compared with (2Vtr) 1@( b + f)=@r. According to Eq. (6.50), these two quantities should equal to each other. In (a)-(c), the total number of configuration generated is 2 105, with a different sampling strategy. In (a), one sample is taken at every 103 configurations, which gives a total of 200 samples to average over. The statistical error in this case is very large. In (b), one sample is taken at every 100 configurations, which gives a total of 2000 samples. The statistical error is reduced compared with (a). In (c), the total number of samples is 2 104. The statistical error is the smallest among (a) to (c). In (d), a total of 2 106 configurations are generated and 2 104 samples are taken at every 100 configurations. sity gradient and the fluctuations [111, 47] in experiment and this thermometry scheme has shown promise to be a reliable way of estimating the temperature [104, 113]. Here we test this method for the Bose-Fermi mixtures and Eqs. (6.50 and 6.51) are extended to mix- tures by considering the density as the total density of both species and the total number as the total number of both species. With the SC method, we calculate the density gradient directly from the density profile expressions. To simulate the fluctuations measured in the experiments, we use a simplified MC simulation explained in the next section. 6.5.2 Fluctuation calculation The MC simulation method generates a large collection of states (or configurations) that satisfies the thermal equilibrium criteria. Such collection of states constitutes a thermal 129 ensemble. In the ensemble, each state (or configuration) gives one density distribution, analogous to one single shot image of the density in the experiment. By averaging over all configurations, one obtains the averaged distribution of particles. Deviations between different configurations are the fluctuations. In our simplified MC method, we use the SC method to determine the density distribution for a given temperature and then use the probability as a reference for configuration generation. The ensemble of configurations is decided to be large enough if it can reproduce the input probabilities. Determining the joint probability : the joint probability, Pjn;m, is the joint probability of having n bosons and m fermions at site j. For m = 1, the joint probability of having n bosons and 1 fermion at site j can be found from the fermionic density distribution, similar to the calculation of the local efficiencyEj (indeed,Ej = Pj1;1), Pjn;1 = Pj(0)n;1 Pj(1)n;1 +Pj(2)n;1 ; (6.52) where we again write the probability as a sum of the probability in the atomic limit, Pj(0)n;1 = b;j(n) f;j(n) Z(0)j ; (6.53) and the contributions from the hopping, Pj(1)n;1 = X k b;j(n) f;j(n) Z(0)j X nb;j;nb;k " b;j(nb;j) b;k(nb;k) Z(0)j Z(0)k f;j(nb;j) f;k(nb;k) f;j(nb;j) f;k(nb;k) (6.54) Pj(2)n;1 = X k X nb;k b;j(n) b;k(nb;k) Z(0)j Z(0)k f;j(n) f;j(n) f;k(nb;k) 130 + f;k(n) f;j(nb;k)[ f;j(n) f;k(nb;k)]2 # : (6.55) Once the joint probability Pjn;1 is determined, the complementary probability Pjn;0 is found based on the relationship in the atomic limit, X n " Pj(0)n;1 + b;j(n) Z(0)j # = 1: (6.56) Taking into the account the hopping contributions, we can write Pjn;0 as, Pjn;0 = b;j(n) Z(0)j +Pj(1)n;1 Pj(2)n;1 : (6.57) We assume each lattice site is independent and the joint probability at site j is sufficient to determine the density distribution at site j. The joint probabilities are evaluated for all the lattice sites and stored in a table before the MC procedure. Simulation procedure: we use a random number generator to generate configurations with reference to the joint probability table. Specifically the simulation includes the fol- lowing steps: 1) Create a table for the values of ePjn;m corresponding to the sum of the joint probability of having up to n bosons and up to m fermions at site j = 1, i.e. ePjn;m = nX k=0 mX l=0 Pjk;l: (6.58) 2) Generate a random number x between 0 and 1. 3) Find the smallest ePjn0;m0 that is larger than x. The number of bosons and fermions at site j is then equal to n0 and m0. 4) Repeat steps (2) and (3) to another site, j = 2, until all the lattice sites are considered. Store the configuration. 5) Repeat steps (2)-(4)N times to generateN configurations. 131 To avoid auto-correlation between adjacent configurations, we choose every otherM 1 configurations as samples. The total number of samples is then Ns =N=M. Averaging over all the samples, we obtain the fermionic and bosonic part of the density fluctuation as f(b)(r) =h f(b)(r)(Nf +Nb)i h f(b)(r)ihNf +Nbi; (6.59) and the total density fluctuation is the sum of f and b. The total number fluctuation is defined as =h(Nf +Nb)2i hNf +Nbi2: (6.60) Here the bracket stands for the averaging over all samples in analogy to the experimental measurement of the fluctuations. 6.5.3 Results The fluctuation calculation is carried out for a 300 300 lattice with all five sets of param- eters. Overall we find very similar behavior for all the parameters and we use parameters Ubb = 11:5t and Ubf = 16t as an example. In our simulation, the fluctuations between different configurations are from both the random number generator and the thermal fluc- tuations. The difference between them is that the thermal fluctuations are independent of ensemble sizes and the sampling size. We find that the correct thermal fluctuation calcula- tion requires a large number of samples ( 104) and large ensemble sizes ( 106). Because of the similarity between the simulation and experimental measurement, this may suggest that a large number of shots are needed in the experiments to obtain the correct thermal fluctuations. Note that we consider here the results from a single plane as one shot, not the averaged results over many planes as reported in Ref. [113]. In Fig. 6.7, we discuss the sampling effects by comparing the fluctuations obtained from different samples with the compressibility calculated from the density gradient. When the number of samples are small, the fluctuations are largely random deviations from the average value. In Fig. 6.7 (a), the fluctuations can be equally positive and negative, which does not even satisfy the condition that the total fluctuations should be always positive. As 132 0 5 10 15 200 5 10 15 20 input temperature (t/kB) Fit temperature (t/k B) T1 T2 T 0.10.20.30.40.50.6?0.2 0 0.2 0.4 0.6 input temperature (t/kB) Fit temperature (t/k B) Figure 6.8: (Color on-line) Extracted temperature as a function of the input temperature. The fluctuations are obtained from 2 104 samples out of 2 106 configurations. The value of T1 is the mean of T1(r) averaged over 12 < r=a < 25 and the error-bar for T1 is the standard deviation in T1(r) [Eq. (6.61)]. The value of T2 is obtained through Eq. (6.62). The input temperature T is drawn as a straight blue line in both plots. In (a), we show our result for the full range of the input temperature, from T = 0:2t=kB to 20t=kB. In this plot, we find very good overall agreement of T1 and T2 with the input temperature for the temperature range considered, particularly for T > 1t=kB. In (b), we blow-up the area inside the black square in (a), which corresponds to the low temperature region, where T = 0:2t=kB to 0:5t=kB. In this region, we find that T1 shows large relative fluctuations (deviation) from the mean value and the mean value of T1 differs relatively greater from T. The extracted temperature T2 however still shows excellent agreement with the input temperature. 133 the number of samples grows, the random noise starts to be averaged out and the fluctua- tions start to agree with the fluctuation-dissipation theorem. In Fig. 6.7 (c), the fluctuations agree very nicely with the relationship predicted by Eq. (6.50). To show that 104 samples are sufficient, we consider an even larger ensemble, with 2 106 configurations [Fig. 6.7(d)] and find that the two ensembles produce almost identical results. This shows that the fluc- tuation calculations obtained in this way are independent of the ensemble size and should correspond to the thermal fluctuations of the system. With Eqs. (6.50) and (6.51), we define two extracted temperatures, T1 and T2. Let T1(r) be the extracted temperature obtained in terms of the density fluctuations and the density gradient, kBT01(r) = f(r) + b(r)(2V tr) 1@( f(r) + b(r))=@r ; 12