Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Axion term in topological insulators with broken time reversal and parity Reis, Itamar 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2014_spring_reis_itamar.pdf [ 949.84kB ]
JSON: 24-1.0166899.json
JSON-LD: 24-1.0166899-ld.json
RDF/XML (Pretty): 24-1.0166899-rdf.xml
RDF/JSON: 24-1.0166899-rdf.json
Turtle: 24-1.0166899-turtle.txt
N-Triples: 24-1.0166899-rdf-ntriples.txt
Original Record: 24-1.0166899-source.json
Full Text

Full Text

Axion term in topological insulators with brokentime reversal and paritybyItamar ReisB.Sc, University of Tel-Aviv, 2010a thesis submitted in partial fulfillmentof the requirements for the degree ofMaster of Scienceinthe faculty of graduate and postdoctoralstudies(Physics)The University Of British Columbia(Vancouver)March 2014c© Itamar Reis, 2014AbstractThe main subject of this work is the axion term in the effective electromag-netic action of topological insulators, which is responsible for the specialelectromagnetic properties of these materials. The axion term is charac-terized by a parameter θ, which can only take the values of 0, for regularinsulators, or pi, for topological insulators, respecting at least one of the timereversal or parity symmetries. A non zero axion term leads to a variety ofmeasurable phenomena, generally referred to as the magneto-electric effects.We focus our interest on the value the parameter θ takes for a topologicalinsulator, when both time reversal and parity are broken. In this case θ nolonger must be quantized to 0 or pi. We use a lattice model for a topologicalinsulator, and introduce a symmetry breaking term in the Hamiltonian.We numerically find the value of θ in this case using calculations of themagnitude of various magneto-electric effects. The results are compared tothe theoretical prediction. We find that θ is no longer quantized when aspecific symmetry breaking term is introduced.iiPrefaceThis dissertation is unpublished work, ultimately based on results from anumerical calculation. These results are compared with known theoreticalpredictions. The code for the calculation was written by the author, I. Reis.The diagrammatic theoretical calculation found in chapter 2 was carried outby my supervisor M. Franz.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . viiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Axion term in a system with broken symmetries . . . . . . 63 Numerical probing of the θ parameter in the axion term . 123.1 Numerical model for a topological insulator . . . . . . . . . . 123.1.1 The spectrum of the model . . . . . . . . . . . . . . . 133.2 Topological phases . . . . . . . . . . . . . . . . . . . . . . . . 143.2.1 Using θ . . . . . . . . . . . . . . . . . . . . . . . . . . 143.2.2 Using ν0 . . . . . . . . . . . . . . . . . . . . . . . . . . 163.3 P and T breaking term . . . . . . . . . . . . . . . . . . . . . 173.3.1 Relation to field theory Lagrangian . . . . . . . . . . . 193.4 Probing θ using magneto-electric effects . . . . . . . . . . . . 213.4.1 m5 domain wall . . . . . . . . . . . . . . . . . . . . . . 21iv3.4.2 Flux insertion . . . . . . . . . . . . . . . . . . . . . . . 243.4.3 Magnetic monopole . . . . . . . . . . . . . . . . . . . 273.4.4 Usage of pi4 rotation symmetry in the numerical calcu-lation . . . . . . . . . . . . . . . . . . . . . . . . . . . 293.5 Physical interpretation of non zero m5 . . . . . . . . . . . . . 304 Calculation of the θ parameter from band structure . . . . 314.1 Low energy contribution . . . . . . . . . . . . . . . . . . . . . 344.2 Full Brillouin zone contribution . . . . . . . . . . . . . . . . . 365 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 375.1 m5 domain wall . . . . . . . . . . . . . . . . . . . . . . . . . . 375.2 θ from band structure . . . . . . . . . . . . . . . . . . . . . . 425.3 Comparison of different numerical calculations . . . . . . . . 456 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52vList of FiguresFigure 2.1 First Feynman diagram for the EM propagator . . . . . . 8Figure 2.2 First Feynman diagram for the EM current . . . . . . . . 9Figure 3.1 m5 domain wall set up . . . . . . . . . . . . . . . . . . . . 22Figure 3.2 Flux insertion setup diagram . . . . . . . . . . . . . . . . 25Figure 3.3 Charge on flux tube Vs. magnetic flux . . . . . . . . . . . 26Figure 5.1 Charge distribution for the m5 domain wall setup . . . . 38Figure 5.2 ∆θ Vs. , numerical results . . . . . . . . . . . . . . . . . 40Figure 5.3 ∆θ Vs. , field theory prediction . . . . . . . . . . . . . . 41Figure 5.4 Comparison between the theoretical prediction and thenumerical results . . . . . . . . . . . . . . . . . . . . . . . 41Figure 5.5 Charge distribution for the m5 domain wall setup, smallband gap . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 5.6 ∆θ Vs. , full BZ prediction . . . . . . . . . . . . . . . . . 44Figure 5.7 θ Vs. m5, different calculation methods . . . . . . . . . . 45Figure 5.8 Charge density for the m5 domain wall setup, using agradual domain wall . . . . . . . . . . . . . . . . . . . . . 46Figure 5.9 Charge density for the flux insertion setup . . . . . . . . . 47Figure 5.10 Charge density in the magnetic monopole setup . . . . . . 48Figure 5.11 Charge accumulated up to distance R from a magneticmonopole Vs. R . . . . . . . . . . . . . . . . . . . . . . . 49viAcknowledgmentsI would like to thank my supervisor M. Franz for introducing me to thesubject and for his advice and comments throughout the work on this thesis.Thanks also goes to D. Marchand for his help with computational issues andto M. Vazifeh for suggestions during the weekly group meetings.viiDedicationDedicated to B., Z., B. and S.viiiChapter 1IntroductionA special property of topological insulators is a unique effective electromag-netic action (ref. [1] and [2]). Like all insulators, the effective electromag-netic action inside topological insulators contains the regular partS0 =12∫d4x(E2 −1µB2), (1.1)but unlike regular insulators, for topological insulators there is an extraterm, called the axion term. That is, for topological insulators the effectiveelectromagnetic action is S = S0 + Sθ, withSθ =∫d4xθ2pie22piE ·B =∫d4xθαpiE ·B. (1.2)Here α, the fine structure constant, isα =e24pi. (1.3)We work in natural units~ = c = 1. (1.4)1The two Maxwell’s equations that are altered due to the axion term areGauss’ and Ampere’s laws∇ ·E = ρ−e24pi2∇θ ·B∇×B =∂E∂t+ j +e24pi2(∇θ ×E +∂θ∂tB).(1.5)The name ’axion’ comes from the name of a hypothetical elementaryparticle, suggested in ref. [3] in 1977 and later discussed in ref [4] and [5],in relation to the problem of strong charge parity, in the field of Quantumchromodynamics. The effect of such a particle on electromagnetism in vac-uum was found in 1987 by F. Wilczek (ref. [6]). In his paper he showsthat the existence of the axion particle will create an additional term in theeffective electromagnetic action, equal to the above axion term.In the years past since their prediction, axions were not experimentallydiscovered. Searches still continue today, the axion being a possible darkmatter composite (ref [7], [8] and [9]). Although the axion itself was notfound yet, its effective electromagnetic action, the axion term, does appearin nature, in topological insulators. The difference between the two cases isthe physical interpretation of θ. For axions it is the field of the particle, andfor topological insulators θ is a parameter describing the material.At first look, the axion term might seem to break both time reversaland parity symmetries. We know that under parity: E→ −E and B→ B,and so Sθ → −Sθ. Under time reversal we have a similar story: E → Eand B → −B, and again Sθ → −Sθ. What saves the day is the fact thatthe axion term is invariant to changes of 2pi in the θ parameter. That is,the partition function and all physical quantities are not effected by thetransformation θ → θ + 2pin, where n is an integer. This is shown in ref.[10]. This fact tells us that θ = pi is equivalent to θ = −pi, which means thatwith this value of θ, the axion term respects both P and T . We see that wehave exactly two options that respect the symmetries, θ = pi and the trivialθ = 0. For regular insulators of course θ = 0, and for topological insulatorswe have θ = pi.2On a surface between topological insulators and regular insulators thevalue of θ needs to change from pi to 0. Although P is generally brokenon surfaces, T should still be respected. This fact does not allow θ to bechanged continuously. The result is a conducting surface state for which θis not defined, and as such allows the transition between 0 and pi. The exis-tence and the special properties of this surface state are the most interestingfeatures of topological insulators.The original theoretical proposal of topological insulators (ref. [11] and[12]) predict the existence of the metallic surface state. This surface statewas also the smoking gun to be later experimentally discovered (ref. [13],[14] and [15]) in materials such as Bi2Se3 and Bi2Te3, to prove the existenceof the topological insulator state of matter.A different option for the surface between a topological insulator and aregular insulator appears when T is artificially broken on the surface, forexample with a ferromagnetic film. In this case θ can change continuouslyand we can get the properties of the surface state from the axion term. Thequantity of interest is the electric current on the surfacejµ =∂L∂Aµ. (1.6)To get the current we write the axion term in a different form:Sθ =∫d4xLθ =∫d4xθ2pie216piµναβFµνFαβ. (1.7)The current isjµ =e24pi2µναβ∂νθ∂αAβ. (1.8)From this current we can get the characteristics of the surface. For examplea half quantized hall conductance, taking the z axis to be perpendicular tothe surface we havejx =12e2hEy →σxy =12e2h.(1.9)3Where this is the total current flowing in the region in which θ(z) changesfrom pi to 0.We saw that one can distinguish topological insulators from regular insu-lators by the value of θ. A different way to characterize insulators is with aZ2 topological invariant as introduced in a series of papers by Fu, Kane andMele (ref. [11], [12] and [16]). The motivation behind the definition of thistopological invariant, called ν0, is coming from the metallic surface states,as opposed to the definition of θ coming from the electromagnetic response.The definition is such that a material with non-zero ν0 will have metallicsurface states. Although defined very differently, there is an equivalencebetween θ and ν0, as proved in ref. [17],(−1)ν0 = eiθ. (1.10)Later we present the ways to calculate the topological invariants and do thecalculation for a lattice model.Next we turn to the case of a system with broken T and P symmetries.In this case the argument described above, for the quantization of the θparameter, does not apply. In rest of the work we discuss the value θ takesfor a topological insulator with broken symmetries.Some discussions of this problem are found in the literature. Ref. [18]suggested the idea of a dynamical θ field. To have a dynamical θ field thefirst thing to do is remove the quantization condition from θ, this is doneby breaking T and P. In the paper a calculation of θ using band structuremethods suggested that there is a correction to θ coming from a specifickind of symmetry breaking term. In this case the form of the correction isnot discussed.An explicit formula for the relation between θ and the magnitude ofthe symmetry breaking term can be found in ref. [19]. As will be discussedshortly, many properties of topological insulators can be described using fieldtheory approach. In the paper this field theory description of topologicalinsulators is used to find the form θ takes when the symmetries are broken.The same relation is obtained using a calculation similar to the one found4in the well known paper by Goldstone and Wilczek in ref. [20].In the next section we present the field theory calculation for θ. Inlater sections we test this prediction using a lattice model for a topologicalinsulator, with a symmetry breaking term. We numerically find the valueof θ in this case using calculations of the magnitude of magneto-electriceffects, such as the Witten effect (ref. [21] and [22]) and charge accumulatedin the flux insertion setup (ref. [23]). The results are also compared to aband structure calculation of θ, using the same lattice model. For the bandstructure calculation we follow formulas from ref. [1] and [18].5Chapter 2Axion term in a system withbroken symmetriesThe electrons occupying the low energy states in a topological insulatorbehave as Dirac fermions and can be described by the Dirac LagrangianL = Ψ¯γµ (i∂µ − eAµ) Ψ +mΨ¯Ψ. (2.1)Many properties of topological insulators can be derived by only consideringthis low energy theory (ref. [1]). This fact is unique when compared tothe much more complicated theories usually needed to understand materialproperties.The appearance of the axion term in the effective electromagnetic actionis one of the properties we can understand starting from the Lagrangianabove. Starting from this Lagrangian we do a chiral rotation,Ψ→ exp(iβ2γ5)Ψ, (2.2)whereγ5 = iγ0γ1γ2γ3, {γµ, γ5} = 0, γ25 = 1. (2.3)6Due to the following propertyΨ¯ = Ψ†γ0 → Ψ† exp(−iβ2γ5)γ0= Ψ†γ0 exp(iβ2γ5)= Ψ¯ exp(iβ2γ5),(2.4)the only term in the Lagrangian effected by this transformation is the masstermmΨ¯Ψ→ mΨ¯ exp(iβ2γ5)Ψ= m cosβΨ¯Ψ + im sinβΨ¯γ5Ψ.(2.5)The γ5 term breaks both P and T . The other effect of the transformationis coming from the path integral measure, using Fujikawa’s method in ref.[24] we haveDΨ¯DΨ→ DΨ¯DΨ exp(∫d4xβ2pie216piµναβFµνFαβ), (2.6)that is, this transformation creates an axion term in effective electromagneticaction. The θ parameter is changed by θ → θ + β. If β = pi we getmΨ¯Ψ → −mΨ¯Ψ,. We see that when the mass changes sign, θ changesbetween 0 and pi. We note that θ is not determined for a specific sign of themass. We can only say that it is changing between 0 and pi (that is, from 0to pi, or from pi to 0) when the sign of the mass changes. When β 6= 0, pi weget a non zero P and T breaking term in the Dirac Lagrangian, which allowsa P and T breaking axion term in the effective electromagnetic action.Next we review a more simple way to get the same result, the correctionto θ due to P and T breaking. We start again from the Dirac LagrangianL = Ψ¯γµ (i∂µ − eAµ) Ψ +mΨ¯Ψ +m5Ψ¯iγ5Ψ. (2.7)7Figure 2.1: First Feynman diagram with fermionic contribution to theelectromagnetic field propagator.We want to check what would be the impact of adding a P and T breakingterm on the value of θ. The formal way to do this is to integrate out thefermions. To do this one calculates the contribution to the electromagneticfield propagator, coming from the interaction with the fermions. The firstFeynman diagram contributing to this is shown in Fig. 2.1. This diagramis UV-divergent with a result that depend on the regulator used.Here we want to use a shortcut following ref. [20]. We use the electriccurrent generated by the axion term, the ’axion current’Sθ =∫d4xLθ =∫d4xθ2pie216piµναβFµνFαβjµ =∂L∂Aµ=e24pi2µναβ∂νθ∂αAβ.We also evaluate the electric current starting from the electrons Lagrangianwith a P and T breaking term, then match the result to the form of theaxion current above. From this procedure we can see what is the relationbetween θ and the m5. We evaluatejµ =∂L∂Aµ= eΨ¯γµΨ, (2.8)in terms of the external fields m5(x) and Aµ(x). This is a perturbativecalculation for the expectation value of the current < Ψ¯γµΨ > using Feyn-man diagrams. m5 is treated as a perturbation. The leading contribution8Figure 2.2: First Feynman diagram contributing to the electromag-netic current. This electromagnetic current is to be comparedwith the ’axion current’, the ’axion current’ being the currentproduced by a non-zero axion term in the effective electromag-netic action.that gives an axion current form comes from the diagram shown in Fig.2.2. There is also a diagram with m5 and Aβ interchanged giving the samecontribution. This diagram gives< jα > = em5(−p− q)Aβ(q)×∫d4k(2pi)4Tr[γα(/k −m)−1eγβ(/k + /q −m)−1γ5(/k − /p−m)−1]= e2m5(−p− q)Aβ(q)Zαβ(p, q).(2.9)The fermion propagator is1/k −m=/k +mk2 −m2, (2.10)and the slash notation is used/k = kµγµ. (2.11)We expand Zαβ(p, q) to first order in p, q. As that is the order that appears9in the axion current we expect to getZαβ(p, q) = pµqν∂2Zαβ∂pµ∂qν∣∣∣∣p,q=0. (2.12)We use∂∂pµ(/k − /p−m)−1∣∣∣∣p=0= −γµk2 −m2+2kµ(/k +m)(k2 −m2)2. (2.13)When the calculation is carried out all the contributions coming from thesecond term are found to add up to zero. We are left with only the firsttermZαβ(p, q) = pµqν∫d4k(2pi)41(k2 −m2)3Tr[γα(/k +m)γβγνγ5(−1)γµ].(2.14)As a trace of an odd number of γ matrices is zero, we see that from the(/k +m) factor only the m contributes. We useTr[γαγβγµγνγ5]= 4iαβµν , (2.15)to getZαβ(p, q) = mpµqν4iαβµν∫d4k(2pi)41(k2 −m2)3= mpµqν4iαβµνI(m2).(2.16)The integral is carried out using analytic continuation to the complex planek → iK,I(m2) = −iΩ4(2pi)4∫ ∞0dKK3(K2 +m2)3=−i32m2pi2, (2.17)where Ω4 is the solid angle in 4 dimensions. We integrate over p and q toget the current in real space, we also add the factor of two coming from thesecond diagram< jα >=e24pi2αβµν∂µm5m∂νAβ. (2.18)Remembering that we are working in the limit m5  m and with a constant10m this result is equal to< jα >=e24pi2αβµν∂µ arctan(m5m)∂νAβ. (2.19)When compared to the axion current coming from the axion term in theelectromagnetic action, this leads toθ(m5) = arctan(m5m)+ Const. (2.20)We know that when m5 = 0, θ = 0, pi and thus the constant is the aboveresult can be 0 or pi.In ref. [19] the same result is obtained by a calculation of the electromag-netic field propagator. We note that it might look like there is a differenceof a factor of two between the result above to the result from ref. [19]. Thisis resolved due to the fact that in ref. [19] the result is for the change in θ ona surface between a region with +m,+m5 and a region with −m,+m5. Inthis case if the correction to θ is indeed arctan(m5m), then on such a surfacewe would get a correction to the change in θ of 2 arctan(m5m).11Chapter 3Numerical probing of the θparameter in the axion termWe want to test the theoretical prediction numerically using a lattice modelfor a topological insulator. In order to numerically probe the value of θ onecan make use of the magneto-electric effects generated by the axion term.All the effects we used can be derived from the axion currentjµ =e24pi2µναβ∂νθ∂αAβ,and specifically from the charge densityj0 = ρ =e24pi2∇θ ·B. (3.1)assuming θ does not depend on time, which is the case here. We numericallycalculate this charge density in different setups and from the result deducethe value of θ.3.1 Numerical model for a topological insulatorFor the numerical calculation we use the method of exact diagonalizationfor a simple topological insulator model used in ref. [23] and [21]. Firstwe introduce the model respecting time reversal and parity, later we add a12symmetry breaking term. We take a cubic lattice with two orbitals per site,denoted c and d. The Hamiltonian consists of two parts, H = HSO + Hcd.The first part is the spin orbit couplingHSO = iλ∑j,µΨ†jτzσµΨj+µ + h.c. (3.2)Here Ψj = (cj,↑, cj,↓, dj,↑, dj,↓)T where c†j , cj are creation and annihilationoperators for the c orbital in the j site. τi are the Pauli matrices in orbitalspace, σi are in spin space. µ = x, y, z. The second partHcd = ∑jΨ†jτxΨj − t∑j,µΨµj τxΨj+µ + h.c, (3.3)contains a spin dependent hopping term. In momentum space this Hamil-tonian takes the form H =∑k Ψ†kHkΨk, whereHk = −2λ∑µτzσµ sin kµ + τxmk (3.4)andmk = − 2t∑µcos kµ. (3.5)This can be written asHk =4∑a=1Γada(k), (3.6)withΓa = (τzσx, τzσy, τzσz, τx) (3.7)andda(k) = (−2λ sin kx,−2λ sin ky,−2λ sin kz,mk) . (3.8)3.1.1 The spectrum of the modelThis Hamiltonian is easily diagonalized due to the fact that it is composedof anti commuting matrices. We can square the Hamiltonian to get H2k =134λ2∑µ sin2 kµ +m2k, and from this we get the spectrumEk = ±√4λ2∑µsin2 kµ +m2k. (3.9)There are two bands, each band is doubly degenerate. This is due to the factthat there are four states per unit cell which gives four states per momentumk. The gap between the bands can close at the following 8 points, called ΓpointsΓ1 = (0, 0, 0) (3.10)Γ2,3,4 = (pi, 0, 0) , (0, pi, 0) , (0, 0, pi) (3.11)Γ5,6,7 = (pi, pi, 0) , (0, pi, pi) , (pi, 0, pi) (3.12)Γ8 = (pi, pi, pi) . (3.13)For Γ1 the gap will close for  = 6t, for Γ2,3,4 at  = 2t , for Γ5,6,7 at  = −2tand for Γ8 at  = −6t.The Γ points are also called time reversal invariant points, as at thosepoints Hk is invariant under time reversal (and also under parity). Thisis due to the fact that at those points k = −k, and as explained later Hktransforms under time reversal to H−k.3.2 Topological phases3.2.1 Using θOne way to find the topological phases of this model is by using the bandstructure formula for θ (ref. [1])θ =14pi∫BZd3kijk Tr[Ai∂jAk + i23AiAjAk], (3.14)14where Aαβi is the Berry connection,Aαβi (k) = −i〈αk∣∣∣∣∂∂ki∣∣∣∣βk〉. (3.15)α, β are band labels, |αk > are eigenfunctions of Hk. The trace is overoccupied states. This formula for θ is generally hard to evaluate, but for theDirac Hamiltonian in use hereHk =4∑a=1Γada(k)Γa = (τzσx, τzσy, τzσz, τx)da(k) = (−2λ sin kx,−2λ sin ky,−2λ sin kz,mk)it is possible to get an analytical expression (ref. [21]). This is done byfinding the eigenstates of the Hamiltonian. The expression isθ =12pi∫BZd3kαβµν1|d(k)|4dα∂kxdβ∂kydµ∂kzdν . (3.16)Following ref. [21], from the structure of the integral we see that in the limitof small mass, mk  λ, the contribution is coming only from the Γ points.In this limit for each Dirac point we expand da(k) around k = Γ to getθΓ =12pi∫d3kvΓ,xvΓ,yvΓ,zmΓ(4λ2k2 +m2Γ)2 , (3.17)wherevΓ,i = −2λ (cos Γx, cos Γy, cos Γz) (3.18)are the Dirac velocities. Due to the fact that the entire contribution is fromsmall k we can take the integration limits to infinity and the integral givesθΓ = −pi2sgn (vΓ,xvΓ,yvΓ,zmΓ) . (3.19)15The total θ isθ =∑ΓθΓ. (3.20)This result is obtained in the limit of mk  λ, but it remains valid whenthe model parameters are continuously changed, as long as mk does not gothrough zero at any of the Γ points. The reason is that θ can not changeunless the gap closes, as θ is quantized to be either 0 or pi for any insulator(with P or T symmetries). When the model parameters are changed con-tinuously θ can not change unless a metallic state occurs, that is when thegap close.The sign of the Dirac velocities only depends on λ, the sign of the massesdepends on the relation between  and t. We assume λ, t > 0. In this casewe can writeθ =pi2(−sgn (+ 6t) + 3sgn (+ 2t)− 3sgn (− 2t) + sgn (− 6t)) . (3.21)The result for different values of  is < −6t −6t <  < −2t −2t <  < 2t 2t <  < 6t 6t < θ 0 pi 0 (2pi) pi 03.2.2 Using ν0A different but equivalent way of finding the topological phases of the modelis by calculating the Z2 invariant. For Hamiltonians with both T and Psymmetries this can be done relatively easily, following ref. [12]. At theΓ points we have Kramer degenerate pairs of states which share the parityeigenvalue due to [T ,P] = 0. The value of ν0 is calculated from the parityeigenvalues at the Γ points. For every Γ point, we take the product of theparity eigenvalues over all the occupied bands, to get δΓδΓ =∏mξm(Γ), (3.22)16where m is the index of the occupied Kramer pairs and ξm = ±1 is theparity eigenvalue of the pair. The Z2 invariant is(−1)ν0 =∏ΓδΓ. (3.23)For the 4×4 Hamiltonian in use here we have, at half filling, one occupiedKramer pair for each Γ point. The Hamiltonian at the Γ points isHΓ = τxmΓ, (3.24)and as explained in the next section, the parity operator (at the Γ points)isP = τx. (3.25)In this case, for the negative energy states, the parity eigenvalue is ′+′ if mΓis negative, and ′−′ if mΓ is positive, that isξΓ = −sgn(mΓ), (3.26)and(−1)ν0 =∏Γsgn(mΓ). (3.27)To get he final result we start from a small  where all the mΓs are negativeand ν0 = 0. When we cross  = −6t, m(pi,pi,pi) changes sign and we getν0 = 1. We continue is this manner to get < −6t −6t <  < −2t −2t <  < 2t 2t <  < 6t 6t < ν0 0 1 0 1 03.3 P and T breaking termIn the numerical calculation we will want to break T and P in the bulk. Todo this we need to know what operator form these symmetries take for themodel in use here. Time Reversal for a single spin 12 particle is implementedby T = Kσy where K is complex conjugate and σy is in spin space, T 2 = −1.17T changes the direction of time and thus for momentum k → −k and forposition x→ x.We want to get a condition on the 4 × 4 Hk for the entire Hamilto-nian H =∑k Ψ†kHkΨk to be invariant under time reversal. We use thetransformation rule for the creation and annihilation operatorsT cjT−1 = cjT ckT−1 = T∑jexp (ijk)cjT −1 =∑jexp (−ijk)T cjT −1 = c−k,(3.28)where the minus sign in the exponent is due to the complex conjugate prop-erty. The transformation of the entire Hamiltonian isT(∑kΨ†kHkΨk)T −1 =∑kΨ†−kT HkT−1Ψ−k =∑kΨ†−kσyH?kσyΨ−k=∑kΨ†kHkΨk,(3.29)and we get the condition for the 4× 4 Hk in the modelσyH?kσy = H−k. (3.30)Parity is implemented by switching between the two orbitals as well as in-verting position and momentum; k→ −k , x→ −x. We getτxHkτx = H−k. (3.31)We want to find a term that breaks those two conditions and also anticommutes with all other four matrices in the Hamiltonian. There is onlyone matrix which satisfies those conditions, τy. We add to our Hamiltoniana term∆Hk = m5τy. (3.32)This is the symmetry breaking term we will use in the numerical cal-culation. We note that this is not the only way to break the symmetries.18Looking at the structure of the modelΓa = (τzσx, τzσy, τzσz, τx, τy) ,da(k) = (−2λ sin kx,−2λ sin ky,−2λ sin kz,mk, 0) ,we can add the following symmetry breaking terms, expanded around a Γpoint∆da = (m1,m2,m3, 0,m5) . (3.33)The low energy calculation suggested that the m5 term will change the valueof θ, this is not the case for m1, m2 and m3. In the rest of the work we willfocus on just the m5 term.3.3.1 Relation to field theory LagrangianThe relation between the model and the Lagrangian used in the field theorycalculation comes from the low energy behavior the model takes near the Γpoints. The low energy expansion of the model have the same form at everyΓ point but with different parametersHk,Γ =∑ıτzσivΓ,iki + τxmΓ + τym5, (3.34)where i = x, y, z andmΓ1 = − 6t ; mΓ2,3,4 = − 2t ; mΓ5,6,7 = + 2t ; mΓ8 = + 6t (3.35)are the masses at the different Γ points. The Dirac velocities arevΓ,i = −2λ (cos Γx, cos Γy, cos Γz) .We want to use this expansion and get the parameters of the Dirac La-grangian,L = Ψ¯γµ (i∂µ − eAµ) Ψ +MΨ¯Ψ +M5Ψ¯iγ5Ψ,19that is, we want to see what is the relation between m5, m and the Diracvelocities from the lattice model, to M5 and M from the Dirac Lagrangian.For this purpose we need the Hamiltonian corresponding to the Dirac La-grangian. We useH =∂L∂(∂0Ψ)∂0Ψ +∂L∂(∂0Ψ†)∂0Ψ† − L, (3.36)with∂L∂(∂0Ψ)= iΨ†∂L∂(∂0Ψ†)= 0,(3.37)to getH = Ψ†(∑i−iγ0γi∇i +Mγ0 +M5iγ0γ5)Ψ, (3.38)and in momentum spaceH = Ψ†(∑iγ0γiki +Mγ0 +M5iγ0γ5)Ψ. (3.39)By changing integration variables in∫d3k, we can divide the lattice Hamil-tonian Hk,Γ by 2λ (assuming λ is positive). For this work we are onlyinterested in the ratio M5M of the Dirac Lagrangian parameters, which doesnot depend on λ. This ratio is what appears in the theoretical predictionwe want to test. There is one thing to be careful about, the sign of M5 inthe Dirac Lagrangian. This is not determined completely by the sign of m5in the lattice Hamiltonian, and in fact depends on the Dirac velocities. Thiscomes from the fact that we need to define the matrices γµ in the Lagrangiandifferently for each Γ pointγµ = (τx,− cos Γxτxτzσx,− cos Γyτxτzσy,− cos Γzτxτzσz) . (3.40)20We haveM5γ5 = −im5τxτy = m5τz, (3.41)but γ5 also needs to satisfyγ5 = iγ0γ1γ2γ3 = −τxi cos Γx cos Γy cos Γziτxτz = τz cos Γx cos Γy cos Γz.(3.42)Those two equations gives the bottom lineM5 = m5 cos Γx cos Γy cos Γz. (3.43)This tells us that different Γ points will contribute to the correction to θwith a different sign for the same m5 in the lattice model. When we look atthe list of Γ points above it’s easy to see that for Γ1, M5 = m5, for Γ2,3,4,M5 = −m5 , for Γ5,6,7, M5 = m5 and for Γ8, M5 = −m5.3.4 Probing θ using magneto-electric effectsTo probe θ numerically we use the magneto-electric effects generated by theaxion term. The magnitude of those effects depends on the value of θ, andby numerically calculating this magnitude we can extract θ.3.4.1 m5 domain wallOne way to probe θ is the following. We add a P and T breaking term∑j Ψ†jm5τyΨj , and we use a non-constant m5 inside the bulk of the topo-logical insulator. We will take m5(x) to be zero in one half of the sampleand some non-zero value in the other half. We also add a constant magneticfield in the x direction B = Bxˆ, A = Byzˆ. The setup is shown in Fig 3.1.We want to look at the charge density induced in the region where m5changes. We remember the formula for the charge densityj0 = ρ =e24pi2∇θ ·B.For a single Γ point the theoretical prediction says θ = Const+arctan(m5/m),21Figure 3.1: m5 domain wall setup. m5 = 0 in half of the topologicalinsulator, and non zero in the other half. A magnetic fieldperpendicular to change inm5 will result in charge accumulatingon the domain wall.where m is the mass at the Γ point. The charge we expect to get (integratedalong x in the region where θ changes) isQ =∫dxρ =e24pi2B∫dx∇xθ =e24pi2∆θB, (3.44)where∆θ = arctan(m5m0)(3.45)is the change in θ across the domain wall.The implementation of the magnetic field in the calculation is done usingthe Pierels substitution, that is, all hopping terms get a factorexp(i2piΦ0∫ jiA · dl), (3.46)where the hopping is between site i and site j andΦ0 =he=2pie(3.47)22is the magnetic flux quantum in natural units. Zeeman terms are not in-cluded.For the vector potential in use, only the hopping terms in the z directionget a factor, which is the followingexp (ieBy). (3.48)We note that the Hamiltonian remains translational invariant in the zdirection after adding the varying m5 and the magnetic field. In this case kzremains an eigenvalue, and instead of diagonalizing the Hamiltonian for a3D lattice , we are left with a kz dependent Hamiltonian on the xy plane. Wediagonalize this Hamiltonian for Lz values of kz ∈ [0, 2pi] which simulatesa lattice with Lz layers in the z direction. The Hamiltonian we need todiagonalize is the result of doing a Fourier transform only in the z direction,to the real space Hamiltonian. H =∑kz Hkz and Hkz = HkzSO +Hkzcd +Hkzz ,whereHkzSO = iλ∑j,µΨ†j,kzτzσµΨj+µ,kz + h.c (3.49)Hkzcd = ∑jΨ†j,kzτxΨj,kz − t∑j,µΨµj,kzτxΨj+µ,kz + h.c (3.50)Hkzz =∑jΨ†j,kz (−2λτzσz sin (kz + eBy)− 2tτx cos (kz + eBy)) Ψj,kz(3.51)∆Hkz =∑jm5(j)Ψ†j,kzτyΨj,kz . (3.52)µ = x, y and j is only in the xy plane.If we work with a L×L lattice on the xy plane, then for each kz we havea 4L2 × 4L2 Hamiltonian. We diagonalize and get the eigenvectors ψn,kz ,n = 1 . . . 4L2. For half filling we take the 2L2 eigenvectors with the negativeeigenvalues. This is due to our expectation that for small perturbations(small m5 and B) we still get two negative and two positive bands like inthe original Hamiltonian, and for each kz separately half the eigenvalues will23be negative and half will be positive. The charge density is calculated fromρ (x, y) =∑n,kz ,l(|ψn,kz (x, y, l)|2 −12). (3.53)l = 1, 2, 3, 4 is for summing over the orbitals and the spins. The −12 is due tothe fact that we want to see the change from the situation with no magneticfield and no m5, in which we have two electrons in each site (at half filling).We have open boundary conditions on the y axis, that is, the sites at y =0 and y = L do not have hopping terms to their left and right respectively.We must use open boundary conditions on the y axis due to the vectorpotential dependence on y. For the x axis we can use either open or periodicboundary conditions. For periodic boundary conditions we have hoppingterms between the sites at x = 0 and the sites at x = L.For the x axis we expect different charge distribution depending on theboundary conditions. For open boundary conditions, the surfaces at x = 0and x = L are surfaces between topological insulator and vacuum so thereshould be charge accumulated there corresponding to the jump in θ, fromthe value it takes inside the topological insulator, pi plus a correction due tom5, and θ = 0 in the vacuum. For periodic boundary conditions, the jumpin θ can only come from different m5 between the x = 0 surface and thex = L surface.3.4.2 Flux insertionA different setup that allows a calculation of θ is the flux insertion. Fol-lowing ref [23] we consider a topological insulator with a thin tube carryingmagnetic flux. The topological insulator is coated by a magnetic film thatbreaks time reversal on the surface, in order to create a gap for the surfacestates. The setup is shown in the Fig. 3.2, taken from ref [23]. In this setup,charge with a magnitude proportional to the value of θ will be accumulatedon the surface of the topological insulator, in the vicinity of the flux tube.We consider the magnetic flux inside the tube being raised from 0 tosome value ηΦ0 = η 2pie . We use Faraday’s law to get the electric field in24Figure 3.2: Flux insertion setup diagram, from ref. [23]. On the topand bottom surfaces of the topological insulator we have ferro-magnetic coating. The yellow arrow is the magnetic flux, the redarrow illustrates the electric field generated when the magneticflux is raised from zero.cylindrical coordinates (ρ, φ, z), and take the magnetic field to be in the zdirection∇×E = −∂B∂t→1ρ∂(ρEφ)∂ρ= −∂Bz∂t,(3.54)andEφ(ρ) = −∂Bz∂tR22ρ, (3.55)where R is the radius of the tube. On the z surfaces of the topological insu-lator we have Hall conductance σxy = e24pi , due to the magnetic coating thatopens a gap on the surface. This result follows from the electromagnetic ac-tion discussed in the introduction, and also from the microscopic description25Figure 3.3: Charge accumulated on the flux tube Vs. the magneticflux, from ref. [23]. The setup is shown in Fig. 3.2. Themagnetic flux in the tube is Φ = ηΦ0, when Φ0 is the magneticflux quantum.of the surface state (ref. [12]). The Hall conductance creates a currentjρ(ρ) =e24piEφ(ρ) = −σxy∂Bz∂tR22ρ. (3.56)The magnetic field is Bz(t) =Φ(t)piR2 . The charge accumulated near the fluxtube isQ =∫ T0dt∫ 2pi0ρdφ jρ(ρ) =∫ T0dt∫ 2pi0ρdφ∂Φ(t)∂tσxy12piρ= ηΦ0σxy = η2piee24pi= ηe2.(3.57)This is the charge in the case of a T and P invariant topological insulator,where θ = pi. In our case, when θ is not quantized at pi we will have adifferent hall conductance σxy = θe24pi2 and the charge will beQ = ηeθ2pi. (3.58)In order to use this setup we set η in the range between 0 and 0.5. From ref[23], the dependence of the charge on η is shown in Fig. 3.3.For η = 1 we expect to get Q = 0, as an integer flux quantum can26be removed by a gauge transformation. The effect at η = 0.5, called thewormhole effect, is the topic of ref [23]. They find that the flux tube turnsconducting for this value of η. A value of θ different than pi is expected tokeep the overall shape the same, but change the slope of the lines between0 and 0.5, and between 0.5 and 1.The vector potential needed for the Pierels substitution in the latticemodel for this setup is A = Aφφˆ withAφ =ηΦ0ρ2piR2ρ < RAφ =ηΦ02piRρ > R.(3.59)We consider a very thin tube for which all the sites are outside of the tube.In this setup we need a T breaking term on the surface. For the i surfacewe add the following term to the Hamiltonian,∆HSi = Ω∑j∈i−surfΨ†jσiΨj , (3.60)where Ω is the strength of the surface magnetization.We use open boundary conditions on the z axis and can use either openor periodic boundary conditions on the x and y axes. For periodic boundaryconditions on x and y we put additional flux tubes in the locations nxLxˆ+nyLyˆ from the original tube, nx, ny being integers. This is necessary in orderto make the system periodic.3.4.3 Magnetic monopoleIn this setup we place a magnetic monopole inside a topological insulator,following ref. [21]. With this setup we have what is called the Witten effect(ref [22]). A charge with a magnitude proportional to θ is accumulated closeto the magnetic monopole.To see the effect of the magnetic monopole we look at Maxwell’s equa-tions in presence of the axion term. A non zero axion term will add extra27terms to Gauss’ law and to Ampere’s law∇ ·E = ρ−e24pi2∇θ ·B∇×B =∂E∂t+ j +e24pi2(∇θ ×E +∂θ∂tB),where ρ and j are the electric charge density and electric current. The twoother equations remain unchanged. For a unit magnetic monopole at theorigin we have∇ ·B = Φ0δ(r)B =Φ04pir2rˆ.(3.61)The divergence of the magnetic field at the origin will destroy the topologicalorder there and we will have θ(r = 0) = 0. Further away from the magneticmonopole the original value of θ will be restored. We think of θ as a functionof the radius with θ(r = 0) = 0 and θ(r = R) = pi for large enough R. Wedo a volume integration of the second term on the right hand side of themodified Gauss’ law, to get an effective electric charge bound to the magneticmonopoleQ = −4pi∫ R0r2e24pi2∂θ∂rΦ04pir2= −e2Φ04pi2∫ R0∂θ∂r= −e2pi.(3.62)In our case, when the symmetries are broken and originally inside the topo-logical insulator θ 6= pi we have Q = − θ2pie.For the Pierels substitution, the vector potential used in the lattice modelis (where (θ, φ) are the spherical angles)A = −Φ0(1 + cos θ)∇φ. (3.63)283.4.4 Usage of pi4 rotation symmetry in the numericalcalculationIn both the flux insertion and the magnetic monopole setups we have a pi4rotation symmetry around the z axis. We can use the symmetry to make theexact diagonalization procedure more efficient. We note that in both caseswe choose the vector potential to respect this symmetry. We start with aL3 lattice and (4L3) × (4L3) Hamiltonian. Using the symmetry we blockdiagonalize the Hamiltonian into four L3 × L3 parts. The diagonalizationof those four L3 × L3 matrices is faster then one (4L3)× (4L3) matrix andallows us to use larger values of L.This is done by writing the Hamiltonian in the basis of eigenstates ofthe pi4 rotation operator. Thepi4 rotation operator mixes the sites in groupsof four. In each of those groups, starting from one site, one can get to theother three by consequent rotations of pi4 around the z axis. In the subspaceof only one of those groups, we can write the operator0 1 0 00 0 1 00 0 0 11 0 0 0. (3.64)We can think of the 1’s as being 14×4 in orbital and spin space. For eacheigenstate of the matrix above we construct 4L4 × L × L eigenstates of thefull pi4 rotation operator.As this is a unitary operator that commutes with the Hamiltonian, eigen-states with different eigenvalues will have zero Hamiltonian matrix element,and so when we write the Hamiltonian in this new basis it will be blockdiagonal.Even after exploiting the symmetry to make the calculation more ef-ficient, the flux tube and magnetic monopole setups remain much slowerthan the m5 domain wall. For this reason we will mostly work with the m5domain wall setup, and use the other two to confirm the results.293.5 Physical interpretation of non zero m5As discussed, topological insulators can be described by the theory of asingle Dirac fermion. That is, by the LagrangianL = Ψ¯γµ (i∂µ − eAµ) Ψ +mΨ¯Ψ.When we add a m5 term to this low energy theory we get the same result forθ for any topological insulator model considered. But for different models,the physical interpretation of the m5 term can be different. For the latticemodel in use here we have∆Hk = m5τy,where τy is in orbital space. While in this case the interpretation of thesymmetry breaking term is not clear, in other topological insulator models,that do describe real solids, there is a physical interpretation. In ref. [18] canbe found a discussion regarding the physical interpretation of a m5 term, fora lattice model describing the low energy bands of the topological insulatorsBi2Te3, Bi2Se3 and Sb2Te3. They find that the m5 term represents in thiscase a staggered Zeeman field, pointing in the opposite direction for thetwo sub-lattices in their model. What can create such a Zeeman field isanti-ferromagnetic order.30Chapter 4Calculation of the θparameter from bandstructureWe can get the value for θ in the model with the P and T breaking termin another way. To get the topological phases of the model without thesymmetry breaking term we used the band structure formula for θθ =14pi∫BZd3kijk Tr[Ai∂jAk + i23AiAjAk].We want to do a similar calculation to see what is the effect of breakingP and T . We can use the results from ref. [1] and ref. [18], for a DiracHamiltonian of the formHk =5∑a=1Γada(k),where Γa are anti commuting matrices. In this case we have∂θ(ξ)∂ξ=34pi∫d3k1|d|5abcdeda∂kzdb∂kydc∂kxdd∂ξde. (4.1)31where ξ is some parameter in the Hamiltonian. We remember that for themodel used in this work we haveΓa = (τzσx, τzσy, τzσz, τx, τy)da(k) = (−2λ sin kx,−2λ sin ky,−2λ sin kz,mk,m5) .This result is obtained by first calculating the Chern number for a DiracHamiltonian in 4 + 1 dimensions (ref. [1], equation (64))C2 =38pi2∫d4Kabcdeda∂Kzdb∂Kydc∂Kxdd∂Kwde|d(K)|5. (4.2)The physical meaning of the Chern number is of less importance here. TheChern number was also shown to be (ref. [1], equation (54))C2 =132pi2∫d4KijklTr[FijFkl], (4.3)whereFij = ∂iAj − ∂jAi, (4.4)andAαβi = −i < αK |∂∂Ki|βK > . (4.5)|βK > is the β eigenstate of HK. Here K = (Kx,Ky,Kz,Kw) is a fourdimensional vector. A dimensional reduction from 4 + 1 to 3 + 1 dimensionscan be performed by taking the momentum in the extra dimension, Kw, tobe a parameter. We denote this parameter with ξ. From the results of thedimensional reduction we see (ref. [1], equation (76))∂θ(ξ)∂ξ=14pi∫d3kijkTr[FξiFjk]. (4.6)32When this is compared to the Chern number formulas, we get the result for∂θ(ξ)∂ξ in terms of the Hamiltonian parameters∂θ(ξ)∂ξ=34pi∫d3k1|d(k)|5abcdeda∂kzdb∂kydc∂kxdd∂ξde. (4.7)The next step is to get a formula for θ, we follow Ref. [18]. We look at thefollowing parameterizationda(k, λ) = (d1(k), d2(k), d3(k), d4(k) + ξ, d5(k)) . (4.8)The special role of d4 here comes from the fact that Γ4 is the only P and Tinvariant matrix in the Hamiltonian, the a = 1, 2, 3 terms are only invariantwith the sin k factors. In any Dirac Hamiltonian we will have only one matrixthat is invariant under P and T (ref. [25]). With this parameterization wehave∂θ(ξ)∂ξ=34pi∫d3kabcd4da∂kzdb∂kydc∂kxdd(|dξ=0|2 − d24 + (d4 + ξ)2) 52, (4.9)andθ(ξ = 0) = θ(ξ =∞)−∫ ∞0dξ∂θ(ξ)∂ξ. (4.10)First we look at the constant θ(ξ =∞). The Hamiltonian at ξ =∞ isH(ξ →∞) = ξΓ4. (4.11)Thanks to the specific choice of parameterization we get a Hamiltonian thatis T symmetric and clearly topologically trivial, so we have θ(ξ = ∞) = 0.Next we evaluate the integral∫ ∞0dξ34pi∫d3kabcd41(|dξ=0|2 − d24 + (d4 + ξ)2) 52da∂kzdb∂kydc∂kxdd.(4.12)33We use ∫ ∞0dx1(a2 − b2 + (b+ x)2)52=(b+ x)(3a2 − b2 + 4bx+ 2x2)3(a2 − b2)2(a2 + x(2b+ x))32∣∣∣∣x=∞x=0=23(a2 − b2)2−b(3a2 − b2)3(a2 − b2)2a3=b(b2 − 2ab+ a2)− 4a2b+ 2a3 + 2ab23(a− b)2(a+ b)2a3=2a+ b3(a+ b)2a3,(4.13)and get the formula for θ for Dirac Hamiltonians with non zero d5θ =14pi∫BZd3k2|d|+ d4(|d|+ d4)2 |d|3ijkldi∂kxdj∂kydk∂kzdl, (4.14)where i, j, k, l = 1, 2, 3, 5. We evaluate the integral twice. First we look atthe contribution from momenta close to the Γ points, this could be doneexplicitly. We also evaluate the entire integral numerically.4.1 Low energy contributionFor a single Γ point we haveda(Γ) = (vx,Γkx, vy,Γky, vz,Γkz,mΓ,m5) ,and we can writeijkldi∂kxdj∂kydk∂kzdl = m5vx,Γvy,Γvz,Γ. (4.15)We take the integration limits to cover the entire k space (which is the sameas field theory calculations). In this case we can get rid of the −2λ factors inthe Dirac velocities by changing variables ki → −2λki. This gives an overallminus sign from changing ∞ to −∞ in the integration limits for each ki.In the end of the day we have an overall sign factor coming from the Dirac34velocities −sgn (vx,Γvy,Γvz,Γ). Without writing this factor we haveθΓ = m5∫ ∞0dkk22√k2 +m20 +m25 +m0(√k2 +m20 +m25 +m0)2√k2 +m20 +m253, (4.16)where this is after the integration over dΩ. This is the contribution fromone Γ point, and m0 = mΓ. We writem = m0 + im5, (4.17)andm = |m| exp iβ. (4.18)We can change integration variables,tanα =k|m|, (4.19)to getθ = sinβ∫ pi20dα2√tan2 α+ 1 + cosβ(√tan2 α+ 1 + cosβ)2√tan2 α+ 13tan2 αcos2 α= sinβ∫ pi20dα2 + cosα cosβ(1 + cosα cosβ)2sin2 α.(4.20)This integral givessinβ2 arctan((cosβ−1) tan α2√1−cos2 β)√1− cos2 β−sinα cosαcosβ cosα+ 1∣∣∣∣α=pi2α=0= 2 arctan(cosβ − 1sinβ)= 2 arctan(−2 sin2 β22 sin β2 cosβ2)= −β,(4.21)and the final result is∆θΓ = βΓ sgn (vx,Γvy,Γvz,Γ) . (4.22)35That is, the contribution to θ coming from each Γ point is the phase of thecomplex mass m = m0 + im5 at the Γ point, with a sign determined by theDirac velocities.4.2 Full Brillouin zone contributionWe want to check whether or not there is a contribution from momentafar from the Γ points. To do this we can evaluate the entire θ integralnumerically. We use a 3D grid for kx, ky, kz and evaluateθ =−2λ3pi∫BZd3k2√−2λ∑i=x,y,z sin ki2 +m2k +m25 +m2k(√. . .+m2k)2√. . .3×m5 cos kx cos ky cos kz,(4.23)where mk = − 2t∑i=x,y,z cos ki and in this case m5 is constant.36Chapter 5Results5.1 m5 domain wallWe start with the m5 setup shown in Fig. 3.1, for which m5 = 0 on half thesample and takes some non zero value on the other half. The charge densitycalculated numerically for this setup generally have a shape as shown in Fig.5.1.We note that the charge density on edges of the sample in the x directionis expected as θ jumps there from pi (up to a correction of order m5m0 ) to 0.The different sign on the two edges is due to different relative directionbetween the jump of θ and the magnetic field. The charge accumulated onthe edges of the sample in the y direction is explained by the fact that thissurface is conducting.We are interested in the charge accumulated due to the jump in m5. Toget this charge we sum over the charge density across the domain wall for aspecific yQ =∑xρ(x, y). (5.1)The range of x over which we sum depends on the sample size and the m5profile in use. It is determined by looking at the charge density distribution,and by comparing the result obtained with different ranges to the theoreticalprediction. The result for this charge is reliable only if it does not vary37Figure 5.1: Three images of the charge distribution for the m5 domainwall setup. All three images are of the same charge density.The first one (top left) is of the entire 14 × 14 xy plane, thesecond one (top right) is of the entire sample in the y directionbut without the edges on the x direction. The third image(bottom) is without the edges on both axes. Without edgesmeans without two sites on each edge. We are mainly interestedin the third image in which we can see clearly the charge densitydue to the jump in m5. m5 = 0.05 on half the sample andm5 = 0 on the other half, L = 14 and we use open boundaryconditions on both the x and y axes. Other model parametersare 2pia2Φ0B = eB = 0.2 (a = 1 is the lattice constant), t = λ = 1, = 4, Lz = 160.38significantly when we change the x range by 1− 2 lattice sites. In this casewe can be sure that there is no overlap between the charge on the m5 domainwall and the charge on the topological insulator surface. We do this sum fora y value at the center of the sample. Again, a reliable result does not varysignificantly with y.We present the results in terms of ∆θ instead of charge∆θ =4pi2LzBQ. (5.2)This is coming from the axion charge density formulaρ =14pi2∆θB.The factor of Lz is because Lz values of kz corresponds to a lattice with Lzsites along the z direction. Our result for the charge density on the xy planeis effectively the sum of the charge density along z.We will look at the dependence of ∆θ on the masses at the Γ points. Forthe eight Γ points in our model the masses arem1 = − 6t ; m2,3,4 = − 2t ; m5,6,7 = + 2t ; m8 = + 6t.We control these masses by changing . Fig. 5.2 is showing the numericalresults for ∆θ, for different values of .The theoretical prediction is that each Γ point will contribute± arctan(m5m0)to θ, where m0 is the mass at the Γ point. The sign is determined by theDirac velocities at the Γ point. The prediction for ∆θ due to P and Tbreaking m5 is∆θ =− arctan(m5+ 6)+ 3 arctan(m5+ 2)− 3 arctan(m5− 2)+ arctan(m5− 6).(5.3)For m5 = 0.05 this theoretical prediction is shown in Fig. 5.3. Next wewant to compare the numerical data to the theoretical prediction. This is39Figure 5.2: Numerical results for the change in θ due to a non zerom5 Vs. .  controls the masses at the Γ points. m5 = 0.05. Inthis case Lz = 100, all other parameters are the same as before(eB = 0.2, m5 = 0.05, t = λ = 1, L = 14 in the rest of the workwe use t = λ = 1, L = 14 unless stated otherwise).done in Fig. 5.4.We see that the general features are similar but the agreement is notexact. It is clear that there is a correction to the value of θ due to m5 andit looks like the theoretical formula is a good prediction for this correction.At the values of  in which the gap closes for one or more of the Γ points,we do not expect to get good results from the numerical calculation. Thereason is that when the gap is small there is a problem with the localizationof the charge on the domain wall. This can be seen in Fig. 5.5. Here  = 5.4,other parameters are the same as in the image of the charge density for  = 4shown in Fig. 5.1. We see that the charge that was localized on the surfacebefore is now spread out on the sample. The charge on the m5 domain wallis hard to see here.40Figure 5.3: Theoretical prediction for the change in θ due to a nonzero m5 Vs.  which controls the masses at the Γ points. m5 =0.05.Figure 5.4: Comparison between the theoretical prediction and thenumerical results. m5 = 0.05.41Figure 5.5: Charge distribution for the m5 domain wall setup, for thecase of small band gap. In this case the charge is not localizedon the domain wall. To be compared with Fig. 5.1We want to find the cause of the discrepancy between the theoreticaland the numerical result for values of  in which the gap is large, to do thiswe turn to the band structure evaluation of θ.5.2 θ from band structureWhen we use band structure formulas to evaluate θ we get the followingintegralθ =−2λ3pi∫BZd3k2√−2λ∑i=x,y,z sin ki2 +m2k +m25 +m2k(√. . .+m2k)2√. . .3m5 cos kx cos ky cos kz.42mk =  − 2t∑i=x,y,z cos ki. We saw that the low energy contribution tothis integral is the same as the field theory prediction. Here we look atthe numerical evaluation of the full Brillouin zone integral. In Fig. 5.6we compare the full Brillouin zone prediction, the low energy prediction,and the results from the lattice model. For a range of  we calculate bothpredictions for ∆θ, the difference in θ between the case with m5 = 0.05 andthe case with m5 = 0. For the lattice model what is shown is the jump in θon a surface between m5 = 0 and m5 = 0.05. That is, the results from theprevious section.We see that other than in the vicinity of  = 6 and  = 2 there is an exactagreement between the lattice model and the full Brillouin zone prediction,as expected. We see that there is a deviation from the low energy prediction.The fact that we see the same deviation in two different calculations suggeststhat it is a real deviation and that the low energy should be thought of asan approximation.We can also look at the dependence of the full θ on m5. For the latticemodel result, to get the full θ, we add pi to the ∆θ we get from the numericalcalculation. This is done in Fig. 5.7. Again we see a good agreement betweenthe full Brillouin zone prediction to the lattice calculation with a deviationfrom the low energy prediction.We conclude that the low energy prediction, is an approximation to theform θ takes in the presence of an im5γ5 term. In the next section we presentmore results obtained also from the flux insertion and monopole setup toconfirm this conclusion.43Figure 5.6: ∆θ Vs. , Comparison between the numerical results andthe low energy and full Brillouin zone predictions. Numericalresults in red, low energy prediction in green, full Brillouin zoneprediction in blue. m5 = 0.05. ∆θ is the difference in θ betweenthe case with m5 = 0.05 and the case with m5 = 0. Bottomfigure is zoomed out.44Figure 5.7: θ Vs. m5, Low energy prediction in green, full Brillouinzone prediction in blue, lattice model results in red.  = 4.5.3 Comparison of different numericalcalculationsIn order to make sure that the deviation of the numerical result in the m5domain wall setup from the low energy prediction is not a numerical error,we compare it to results from the magnetic monopole and flux insertionsetups. We also try a gradual m5 domain wall. This is in order to check ifthe deviation from the low energy prediction, seen in the last section, is aresult of the abrupt change in m5 in the original m5 domain wall setup.First we try gradual domain wall and periodic boundary conditions, inthe m5 domain wall setup. We use the following m5 profilem5(x) =12m5(tanh(x− L4ξ)− tanh(x− 3L4ξ)). (5.4)We calculate the charge density and the charge accumulated on the domainwall for ξ = 1, 2, 3. The charge density in this setup is shown in Fig. 5.8.45Figure 5.8: Charge density for the m5 domain wall setup, using agradual domain wall in which m5 changes gradually from m5 =0 on the surfaces to m5 = 0.1 in the middle of the sample.ξ defined above controls the width of the domain wall. Hereξ = 3 and Lx = 50. Periodic boundary conditions on the xaxis. Other model parameters are  = 4, B = 0.001 Ly = 8,Lz = 100.The change in θ due to m5 is ∆θ = 4pi2BLzQ, Q is the sum of the chargedensity from x = 0 to x = Lx2 for some y. We get∆θξ=1,Lx=40 = 0.214∆θξ=2,Lx=40 = 0.220∆θξ=3,Lx=50 = 0.213.(5.5)(for ξ = 3 we use Lx = 50, as with Lx = 40 the charge distributions fromthe two domain walls starts to overlap). In the original setup (m5 zero onhalf the sample) we had∆θm5,0 = 0.215. (5.6)Next we look at the flux insertion setup depicted in Fig. 3.2, with a constantm5 = 0.1 in the bulk, open boundary conditions, L = 16, Ω = 1 (additional46Figure 5.9: Charge density for the flux insertion setup. The top sur-face of the sample is shown. m5 = 0.1 in the entire sample.L = 16. Φ = ηΦ0 is the magnetic flux inside the tube, η = 0.1.Ω = 1 is the magnitude of the time reversal breaking on thesurface. To get to charge accumulated on the flux tube we sumthe charge on the entire upper half of the sample.time reversal breaking on the surface is needed here, as we use small m5),and all other model parameters the same. The charge density on the topsurface is shown in Fig. 5.9.We calculate the charge QFI in the entire upper half of the sample (z >L2 ). The lower half of the sample will have an opposite charge. To get thechange in θ from the case with m5 = 0 we use∆θFI =2piQFIη− pi. (5.7)Φ = ηΦ0 is the flux inside the flux tube. We get∆θFI = 0.223. (5.8)For the magnetic monopole setup we use m5 = 0.1, open boundaryconditions, L = 16 and all other model parameters the same. The magneticmonopole is located at the middle of the sample. Fig. 5.10 shows the charge47Figure 5.10: Charge density in the magnetic monopole setup. Thez = L2 surface is shown. m5 = 0.1 in the entire sample. L = 16.The magnetic monopole is of unit magnetic charge.density for the z = L2 plane. The charge accumulated up to distance R fromthe monopole is shown in Fig. 5.11.We take the charge accumulated on the magnetic monopole to be oneof the points on the plateau from Fig. 5.11. We note that the charge onthe entire sample is zero, here the negative counter part to the charge nearthe monopole is accumulated on the surface. Similarly to the flux insertionsetup, to get the change in θ from the case with m5 = 0 we use∆θmonopole = 2piQmonopole − pi, (5.9)to get∆θmonopole = 0.219. (5.10)The full Brillouin zone prediction with the same model parameters gives∆θFull BZ = 0.217. (5.11)48Figure 5.11: Charge accumulated up to distance R from a unit mag-netic monopole Vs. R. m5 = 0.1,  = 4, L = 16. We takethe charge accumulated on the monopole from the plateau ofR ' 4 to R ' 6. For large R the charge goes to 0 as the totalcharge in the sample is zero (that is, the same as in the casewith m5 = 0 or without the magnetic monopole).The low energy prediction is∆θLow E = arctan(0.110)− 3 arctan(0.16)+3 arctan(0.12)− arctan(0.1−2)= 0.159.(5.12)We sum the results in the following tablenumerical calculations Full BZ prediction Low E prediction∆θ 0.213 - 0.223 0.217 0.159We see that all the different numerical calculations agree with each other,and with the full Brillouin zone prediction, and that there is a deviation fromthe low energy prediction. We get similar behavior for different model pa-rameters. We also tried using different boundary conditions where possible,as explained in the model description section, the results does not changesignificantly. This supports the conclusion that the low energy prediction isan approximation.49Chapter 6ConclusionsThe main goal of this work was to numerically calculate θ in a system withbroken P and T , and compare the result with the theoretical predictioncoming from various field theory and band structure calculations. Differentlow energy calculations yield the simple result, saying that the contributionto θ from each Γ point is the phase of the complex mass, m = m0 + im5, atthe Γ point (up to a constant being either 0 of pi).The numerical results show a good but not exact agreement with theprediction coming from low energy calculations, and an exact agreementwith a calculation taking into account the full Brillouin zone.We do not have an analytical form for the exact form θ takes. For DiracHamiltonian models this exact form can be evaluated numerically relativelyeasily using the integralθ =14pi∫BZd3k2|d|+ d4(|d|+ d4)2 |d|3ijkldi∂kxdj∂kydk∂kzdl.with the simple low energy result usually being an approximation. Thespecial role d4 takes here is assuming that in the Dirac Hamiltonian, Γ4 isinvariant under time reversal and parity.As mentioned, for a number of topological insulator models describingreal materials, m5 is related to anti-ferromagnetic order. In those cases weconclude that anti-ferromagnetic order will change the magnitude of the50magneto-electric effects. This change could be calculated or approximatedusing the methods described above.51Bibliography[1] Xiao-Liang Qi, Taylor L. Hughes, and Shou-Cheng Zhang.Topological field theory of time-reversal invariant insulators. PRB,78(195424), 2008. → pages 1, 5, 6, 14, 31, 32[2] A. M. Essin, J. E. Moore, and D. Vanderbilt. Magnetoelectricpolarizability and axion electrodynamics in crystalline insulators.PRL, 102(146805), 2009. → pages 1[3] R. D. Peccei and H. R. Quinn. Cp conservation in the presence ofpseudoparticles. PRL, 38(1440), 1977. → pages 2[4] S. Weinber. A new light boson? PRL, 40(223), 1978. → pages 2[5] F.Wilczek. Problem of strong p and t invariance in the presence ofinstantons. PRL, 40(279), 1978. → pages 2[6] F. Wilczek. Two applications of axion electrodynamics. PRL, 58(18),1987. → pages 2[7] L. D. Duffy and K. van Bibber. Axions as dark matter particles. Newjournal of physics, 11(105008), 2009. → pages 2[8] L. D. Duffy, P. Sikivie, D. B. Tanner, S. J. Asztalos, C. Hagmann,D. Kinion, L. J. Rosenberg, K. van Bibber, D. B. Yu, , and R. F.Bradley. High resolution search for dark-matter axions. PRD,74(012006), 2006. → pages 2[9] W. Beiglbck R. Beig and W. Domcke. Axions: Theory, Cosmology,and Experimental Searches. Springer, New York, 2008. → pages 2[10] M.M. Vazifeh and M. Franz. Quantization and 2pi periodicity of theaxion action in topological insulators. PRB, 82(233103), 2010. →pages 252[11] L. Fu, C. L. Kane, and E. J. Mele. Topological insulators in threedimensions. PRL, 98(106803), 2007. → pages 3, 4[12] L. Fu and C. L. Kane. Topological insulators with inversionsymmetry. PRB, 76(045302), 2007. → pages 3, 4, 16, 26[13] D. Hsieh, D. Qian, Y. Xia L. Wray, Y. S. Hor, R. J. Cava, and M. Z.Hasan. A topological dirac insulator in a quantum spin hall phase.Nature (London), 452(970), 2008. → pages 3[14] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil,D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan. Observation of alarge-gap topological-insulator class with a single dirac cone on thesurface. Nature Physics, 5(398), 2009. → pages 3[15] Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo, X. L. Qi,H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S. C. Zhang, I. R. Fisher,Z. Hussain, and Z.-X. Shen. Experimental realization of athree-dimensional topological insulator, bi2te3. Science, 325(178),2009. → pages 3[16] L. Fu and C. L. Kane. Time reversal polarization and a z2 adiabaticspin pump. PRB, 74(195312), 2006. → pages 4[17] Zhong Wang, Xiao-Liang Qi, and Shou-Cheng Zhang. Equivalenttopological invariants of topological insulators. New journal ofphysics, 12(065007), 2010. → pages 4[18] Rundong Li, Jing Wang, Xiao-Liang Qi, and Shou-Cheng Zhang.Dynamical axion field in topological magnetic insulators. NaturePhysics, (1534), 2010. → pages 4, 5, 30, 31, 33[19] Michael Mulligan and F. J. Burnell. Topological insulators avoid theparity anomaly. PRB, 88(085104), 2013. → pages 4, 11[20] J. Goldstone and F. Wilczek. Fractional quantum numbers onsolitons. PRL, 47(14), 1981. → pages 5, 8[21] G. Rosenberg and M. Franz. Witten effect in a crystalline topologicalinsulator. PRB, 82(035105), 2010. → pages 5, 12, 15, 27[22] E. Wittten. PRL, 86(283), 1979. → pages 5, 2753[23] G. Rosenberg, H.-M. Guo, and M. Franz. Wormhole effect in a strongtopological insulator. PRB, 82(041104(R)), 2010. → pages 5, 12, 24,25, 26, 27[24] K. Fujikawa and H. Suzuki. Path Integrals and Quantum Anomalies.Oxford University Press, 2004. → pages 7[25] M. E. Peskin and D. V. Schroeder. An Introduction to Quantum FieldTheory. Perseus Press, Cambridge, MA, 1995. → pages 3354


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items