From 6c09cade95e20c65bef6abb49b06dd06cfe9559c Mon Sep 17 00:00:00 2001 From: Tom Hodson Date: Wed, 25 Jan 2023 16:04:40 +0000 Subject: [PATCH] Update thesis --- _thesis/0_Preface/0.1_Abstract.html | 2 +- _thesis/1_Introduction/1_Intro.html | 64 +- _thesis/2_Background/2.1_FK_Model.html | 38 +- _thesis/2_Background/2.2_HKM_Model.html | 46 +- _thesis/2_Background/2.4_Disorder.html | 12 +- .../3.1_LRFK_Model.html | 8 +- .../3.2_LRFK_Methods.html | 29 +- .../3.3_LRFK_Results.html | 26 +- .../4.1_AMK_Model.html | 26 +- .../4.2_AMK_Methods.html | 6 +- .../4.3_AMK_Results.html | 18 +- _thesis/5_Conclusion/5_Conclusion.html | 2 +- .../A.1.2_Fermion_Free_Energy.html | 14 +- .../A.1_Particle_Hole_Symmetry.html | 12 +- .../A.2_Markov_Chain_Monte_Carlo.html | 46 +- .../6_Appendices/A.3_Lattice_Generation.html | 4 +- .../6_Appendices/A.4_Lattice_Colouring.html | 6 +- _thesis/6_Appendices/A.5_The_Projector.html | 18 +- _thesis/toc.html | 2 +- .../lattice_construction_animated.gif | Bin 842418 -> 831984 bytes .../lattice_construction_animated.svg | 3416 +++++++++-------- .../stokes_theorem/stokes_theorem.svg | 783 ++-- .../binder_cumulants/binder_cumulants.svg | 2384 +++++------- 23 files changed, 3279 insertions(+), 3683 deletions(-) diff --git a/_thesis/0_Preface/0.1_Abstract.html b/_thesis/0_Preface/0.1_Abstract.html index c3bf04a..c317ebf 100644 --- a/_thesis/0_Preface/0.1_Abstract.html +++ b/_thesis/0_Preface/0.1_Abstract.html @@ -51,7 +51,7 @@ image: - @@ -137,14 +136,14 @@ H_{\mathrm{FK}} = & -\;t \sum_{\langle i,j \rangle} c^\dagger_{i}c_{j} + \;U

Spin-orbit coupling is a relativistic effect that, very roughly, corresponds to the fact that in the frame of reference of a moving electron the electric field of nearby nuclei looks like a magnetic field to which the electron spin couples. This couples the spatial and spin parts of the electron wavefunction. The lattice structure can therefore influence the form of the spin-spin interactions, leading to spatial anisotropy in the effective interactions. This spatial anisotropy can frustrate an MI state  [60,61] leading to more exotic ground states than the AFM order we have seen so far. As with the Hubbard model, interaction effects are only strong or weak in comparison to the bandwidth or hopping integral \(t\). Hence, we will see strong frustration in materials with strong spin-orbit coupling \(\lambda\) relative to their bandwidth \(t\).

-

In certain transition metal based compounds, such as those based on Iridium and Ruthenium, the lattice structure, strong spin-orbit coupling and narrow bandwidths lead to effective spin-\(\tfrac{1}{2}\) Mott insulating states with strongly anisotropic spin-spin couplings. These transition metal compounds, known as Kitaev materials, draw their name from the celebrated Kitaev Honeycomb (KH) model which is expected to model their low temperature behaviour  [56,62–65].

-

At this point, we can sketch out a phase diagram like that of fig. 4. When both electron-electron interactions \(U\) and spin-orbit couplings \(\lambda\) are small relative to the bandwidth \(t\), we recover standard band theory of band insulators and metals. In the upper left, we have the simple Mott insulating state as described by the Hubbard model. In the lower right, strong spin-orbit coupling gives rise to topological insulators characterised by symmetry protected edge modes and non-zero Chern number. Kitaev materials occur in the region where strong electron-electron interaction and spin-orbit coupling interact. See  [66] for a more expansive version of this diagram.

-

The KH model  [67] was one of the first exactly solvable spin models with a QSL ground state. It is defined on the 2D honeycomb lattice and provides an exactly solvable model that can be reduced to a free fermion problem via a mapping to Majorana fermions. This yields an extensive number of static \(\mathbb Z_2\) fluxes tied to an emergent gauge field. The model is remarkable not only for its QSL ground state, but also for its fractionalised excitations with non-trivial braiding statistics. It has a rich phase diagram hosting gapless, Abelian and non-Abelian phases  [68] and a finite temperature phase transition to a thermal metal state  [69]. It has been proposed that its non-Abelian excitations could be used to support robust topological quantum computing  [70–72].

+

In certain transition metal based compounds, such as those based on iridium and ruthenium, the lattice structure, strong spin-orbit coupling and narrow bandwidths lead to effective spin-\(\tfrac{1}{2}\) Mott insulating states with strongly anisotropic spin-spin couplings. These transition metal compounds, known as Kitaev materials, draw their name from the celebrated Kitaev Honeycomb (KH) model which is expected to model their low temperature behaviour  [56,62–65].

+

At this point, we can sketch out a phase diagram like that of fig. 4. When both electron-electron interactions \(U\) and spin-orbit couplings \(\lambda\) are small relative to the bandwidth \(t\), we recover standard band theory of band insulators and metals. In the upper left, we have the simple Mott insulating state as described by the Hubbard model. In the lower right, strong spin-orbit coupling gives rise to topological insulators characterised by symmetry protected edge modes and non-zero Chern number. Kitaev materials occur in the region where strong electron-electron interaction and spin-orbit coupling interact. See ref.  [66] for a more expansive version of this diagram.

+

The KH model  [67] was one of the early exactly solvable spin models with a QSL ground state. It is defined on the 2D honeycomb lattice and provides an exactly solvable model that can be reduced to a free fermion problem via a mapping to Majorana fermions. This yields an extensive number of static \(\mathbb Z_2\) fluxes tied to an emergent gauge field. The model is remarkable not only for its QSL ground state, but also for its fractionalised excitations with non-trivial braiding statistics. It has a rich phase diagram hosting gapless, Abelian and non-Abelian phases  [68] and a finite temperature phase transition to a thermal metal state  [69]. It has been proposed that its non-Abelian excitations could be used to support robust topological quantum computing  [70–72].

The KH and FK models have quite a bit of conceptual overlap. They can both be seen as models of spinless fermions coupled to a classical Ising background field. This is what makes them exactly solvable. At finite temperatures, fluctuations in their background fields provide an effective disorder potential for the fermionic sector, so both models can be studied at finite temperature with Markov chain Monte Carlo methods  [69,73].

As Kitaev points out in his original paper, the KH model remains solvable on any trivalent \(z=3\) graph which can be three-edge-coloured. Indeed, many generalisations of the model exist  [74–78]. Notably, the Yao-Kivelson model  [79] introduces triangular plaquettes to the honeycomb lattice leading to spontaneous chiral symmetry breaking. These extensions all retain translation symmetry. This is likely because edge-colouring, finding the ground state and understanding the QSL properties are much harder without it  [80,81]. Undeterred, this gap led us to wonder what might happen if we remove translation symmetry from the Kitaev model. This would be a model of a trivalent, highly bond anisotropic but otherwise amorphous material.

-

Amorphous materials do not have long-range lattice regularities but may have short-range regularities in the lattice structure, such as fixed coordination number \(z\) as in some covalent compounds. The best examples are amorphous Silicon and Germanium with \(z=4\) which are used to make thin-film solar cells  [82,83]. Recently, it has been shown that topological insulating (TI) phases can exist in amorphous systems. Amorphous TIs are characterised by similar protected edge states to their translation invariant cousins and generalised topological bulk invariants  [84–90]. However, research on amorphous electronic systems has mostly focused on non-interacting systems with a few exceptions, for example, to account for the observation of superconductivity  [91–95] in amorphous materials or very recently to understand the effect of strong electron repulsion in TIs  [96].

+

Amorphous materials do not have long-range lattice regularities but may have short-range regularities in the lattice structure, such as fixed coordination number \(z\) as in some covalent compounds. The best examples are amorphous silicon and germanium with \(z=4\) which are used to make thin-film solar cells  [82,83]. Recently, it has been shown that topological insulating (TI) phases can exist in amorphous systems. Amorphous TIs are characterised by similar protected edge states to their translation invariant cousins and generalised topological bulk invariants  [84–90]. However, research on amorphous electronic systems has mostly focused on non-interacting systems with a few exceptions, for example, to account for the observation of superconductivity  [91–95] in amorphous materials or very recently to understand the effect of strong electron repulsion in TIs  [96].

Amorphous magnetic systems have been investigated since the 1960s, mostly through the adaptation of theoretical tools developed for disordered systems  [97–100] and with numerical methods  [101,102]. Research on classical Heisenberg and Ising models accounts for the observed behaviour of ferromagnetism, disordered antiferromagnetism and widely observed spin glass behaviour  [103]. However, the role of spin-anisotropic interactions and quantum effects in amorphous magnets has not been addressed.

-

In chapter 4, I will address the question of whether frustrated magnetic interactions on amorphous lattices can give rise to genuine quantum phases, i.e. to long-range entangled QSL  [49,104–106]. We will find that the answer is yes. I will introduce the amorphous Kitaev model, a generalisation of the KH model to random lattices with fixed coordination number three. I will show that this model is a solvable, amorphous, chiral spin liquid. As with the Yao-Kivelson model  [79], the amorphous Kitaev model retains its exact solubility but the presence of plaquettes with an odd number of sides leads to a spontaneous breaking of time reversal symmetry. I will confirm prior observations that the form of the ground state is relatively simple  [77,107] and unearth a rich phase diagram displaying Abelian as well as a non-Abelian chiral spin liquid phases. Furthermore, I will show that the system undergoes a finite-temperature phase transition to a thermal metal state and discuss possible experimental realisations.

+

In chapter 4, I will address the question of whether frustrated magnetic interactions on amorphous lattices can give rise to genuine quantum phases, i.e., to long-range entangled QSL  [49,104–106]. We will find that the answer is yes. I will introduce the amorphous Kitaev model, a generalisation of the KH model to random lattices with fixed coordination number three. I will show that this model is a solvable, amorphous, chiral spin liquid. As with the Yao-Kivelson model  [79], the amorphous Kitaev model retains its exact solubility but the presence of plaquettes with an odd number of sides leads to a spontaneous breaking of time reversal symmetry. I will confirm prior observations that the form of the ground state is relatively simple  [77,107] and unearth a rich phase diagram displaying Abelian as well as a non-Abelian chiral spin liquid phases. Furthermore, I will show that the system undergoes a finite-temperature phase transition to a thermal metal state and discuss possible experimental realisations.

The next chapter, Chapter 2, will introduce some necessary background to the FK model, the KH model, and disorder and localisation.

Next Chapter: 2 Background

@@ -477,7 +476,8 @@ H_{\mathrm{FK}} = & -\;t \sum_{\langle i,j \rangle} c^\dagger_{i}c_{j} + \;U
    -

  1. In fact the FQH state can be described within a Ginzburg-Landau like paradigm but it requires a non-local order parameter  [51,52]. Phases like this are said to possess topological order defined by the fact that they cannot be smoothly deformed into a product state  [53,54].↩︎

    2. +

  2. The Pauli exclusion principle is special in that it can be treated much more simply that other interactions, this relates to the fact that it is a hard constraint.↩︎

    3. +

  3. In fact the FQH state can be described within a Ginzburg-Landau like paradigm but it requires a non-local order parameter  [51,52]. Phases like this are said to possess topological order defined by the fact that they cannot be smoothly deformed into a product state  [53,54].↩︎

diff --git a/_thesis/2_Background/2.1_FK_Model.html b/_thesis/2_Background/2.1_FK_Model.html index 567a0c3..bd8c336 100644 --- a/_thesis/2_Background/2.1_FK_Model.html +++ b/_thesis/2_Background/2.1_FK_Model.html @@ -83,32 +83,32 @@ image:

The Model

-Figure 1: The dispersion (upper row) and density of states (lower row) obtained from a cubic lattice model H = \sum_{i} V_i c^\dagger_{i}c_{i} - t \sum_{\langle i,j\rangle} c^\dagger_{i}c_{j} in 1D. (a) With no external potential. (b) With a static charge density wave background V_i = (-1)^i (c) A static charge density wave background with 2% binary disorder. The top row shows the analytic dispersion in orange compared with the integral of the DOS in dotted black. - +Figure 1: The dispersion (upper row) and density of states (lower row) obtained from a cubic lattice model H = \sum_{i} V_i c^\dagger_{i}c^{\phantom{\dagger}}_{i} - t \sum_{\langle i,j\rangle} c^\dagger_{i}c^{\phantom{\dagger}}_{j} in 1D. (a) With no external potential. (b) With a static charge density wave background V_i = (-1)^i (c) A static charge density wave background with 2% binary disorder. The top row shows the analytic dispersion in orange compared with the integral of the DOS in dotted black. +

The Falicov-Kimball (FK) model is one of the simplest models of the correlated electron problem. It captures the essence of the interaction between itinerant and localised electrons. It was originally introduced to explain the metal-insulator transition in f-electron systems. However, in its long history, the FK model has been interpreted variously as a model of electrons and ions, binary alloys or crystal formation  [1–4]. In terms of immobile fermions \(d_i\) and light fermions \(c_i\) and with chemical potential fixed at half-filling, the model reads

\[\begin{aligned} -H_{\mathrm{FK}} = & \;U \sum_{i} (d^\dagger_{i}d_{i} - \tfrac{1}{2})\;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) -\;t \sum_{\langle i,j\rangle} c^\dagger_{i}c_{j}.\\ +H_{\mathrm{FK}} = & \;U \sum_{i} (d^\dagger_{i}d_{i} - \tfrac{1}{2})\;(c^\dagger_{i}c^{\phantom{\dagger}}_{i} - \tfrac{1}{2}) -\;t \sum_{\langle i,j\rangle} c^\dagger_{i}c^{\phantom{\dagger}}_{j}.\\ \end{aligned}\]

-

Here we will only discuss the hypercubic lattices, i.e the chain, the square lattice, the cubic lattice and so on. The connection to the Hubbard model is that we have relabelled the up and down spin electron states and removed the hopping term for one spin state. This is equivalent to taking the limit of infinite mass ratio  [5].

-

Like other exactly solvable models  [6], the FK model possesses extensively many conserved degrees of freedom \([d^\dagger_{i}d_{i}, H] = 0\). Similarly, the Kitaev model contains an extensive number of conserved fluxes. So in both models, the Hilbert space breaks up into a set of sectors in which these operators take a definite value. Crucially, this reduces the interaction terms in the model from being quartic in fermion operators to quadratic. This is what makes the two models exactly solvable, in contrast to the Hubbard model. For the FK model the interaction term \((d^\dagger_{i}d_{i} - \tfrac{1}{2})\;(c^\dagger_{i}c_{i} - \tfrac{1}{2})\) becomes quadratic when \(d^\dagger_{i}d_{i}\) is replaced with one of its eigenvalues \(\{0,1\}\). The same thing happens in the Kitaev model, though after first applying a clever transformation which we will discuss later.

-

Due to Pauli exclusion, maximum filling occurs when each lattice site is fully occupied, \(\langle n_c + n_d \rangle = 2\). Here we will focus on the half filled case \(\langle n_c + n_d \rangle = 1\). The ground state phenomenology as the model is doped away from the half-filled state can be rich  [7,8] but the half-filled point has symmetries that make it particularly interesting. From this point on, we will only consider the half-filled point.

-

At half-filling and on bipartite lattices, the FK model is particle-hole symmetric. That is, the Hamiltonian anticommutes with the particle hole operator \(\mathcal{P}H\mathcal{P}^{-1} = -H\). As a consequence, the energy spectrum is symmetric about \(E = 0\), which is the Fermi energy. The particle hole operator corresponds to the substitution \(c^\dagger_i \rightarrow \epsilon_i c_i, d^\dagger_i \rightarrow d_i\) where \(\epsilon_i = +1\) for the A sublattice and \(-1\) for the B sublattice  [4]. The absence of a hopping term for the heavy electrons means they do not need the factor of \(\epsilon_i\) but they would need it in the corresponding Hubbard model. See appendix A.1 for a full derivation of the particle-hole symmetry.

-

We will later add a long-range interaction between the localised electrons. At that point we will replace the immobile fermions with a classical Ising field \(S_i = 1 - 2d^\dagger_id_i = \pm1\) which I will refer to as the spins.

+

Here we will only discuss the hypercubic lattices, i.e., the chain, the square lattice, the cubic lattice and so on. The connection to the Hubbard model is that we have relabelled the up and down spin electron states and removed the hopping term for one spin state. This is equivalent to taking the limit of infinite mass ratio  [5].

+

Like other exactly solvable models  [6], the FK model possesses extensively many conserved degrees of freedom \([d^\dagger_{i}d_{i}, H] = 0\). Similarly, the Kitaev model contains an extensive number of conserved fluxes. So in both models, the Hilbert space breaks up into a set of sectors in which these operators take a definite value. Crucially, this reduces the interaction terms in the model from being quartic in fermion operators to quadratic. This is what makes the two models exactly solvable, in contrast to the Hubbard model. For the FK model the interaction term \((d^\dagger_{i}d_{i} - \tfrac{1}{2})\;(c^\dagger_{i}c^{\phantom{\dagger}}_{i} - \tfrac{1}{2})\) becomes quadratic when \(d^\dagger_{i}d_{i}\) is replaced with one of its eigenvalues \(\{0,1\}\). The same thing happens in the Kitaev model, though after first applying a clever transformation which we will discuss later.

+

Due to Pauli exclusion, maximum filling occurs when each lattice site is fully occupied, \(\langle n_c + n_d \rangle = 2\). Here we will focus on the half-filled case \(\langle n_c + n_d \rangle = 1\). The ground state phenomenology as the model is doped away from the half-filled state can be rich  [7,8] but the half-filled point has symmetries that make it particularly interesting. From this point on, we will only consider the half-filled point.

+

At half-filling and on bipartite lattices, the FK model is particle-hole symmetric. That is, the Hamiltonian anticommutes with the particle hole operator \(\mathcal{P}H\mathcal{P}^{-1} = -H\). As a consequence, the energy spectrum is symmetric about \(E = 0\), which is the Fermi energy. Fig. 1 shows this in action, in the presence of a periodic potential a gap in the energy spectrum opens symmetrically about \(E = 0\). The particle hole operator corresponds to the substitution \(c^\dagger_i \rightarrow \epsilon_i c_i, d^\dagger_i \rightarrow d_i\) where \(\epsilon_i = +1\) for the A sublattice and \(-1\) for the B sublattice  [4]. The absence of a hopping term for the heavy electrons means they do not need the factor of \(\epsilon_i\) but they would need it in the corresponding Hubbard model. See appendix A.1 for a full derivation of the particle-hole symmetry.

+

We will later add a long-range interaction between the localised electrons. At that point we will replace the immobile fermions with a classical Ising field \(S_i = \tfrac{1}{2}(1 - 2d^\dagger_id_i) = \pm\tfrac{1}{2}\) which I will refer to as the spins.

\[\begin{aligned} -H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) -\;t \sum_{\langle i,j\rangle} c^\dagger_{i}c_{j}.\\ +H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c^{\phantom{\dagger}}_{i} - \tfrac{1}{2}) -\;t \sum_{\langle i,j\rangle} c^\dagger_{i}c^{\phantom{\dagger}}_{j}.\\ \end{aligned}\]

The FK model can be solved exactly with dynamic mean field theory in the infinite dimensional limit  [9–12]. In lower dimensional systems it has radically different behaviour, as we shall see.

Phase Diagrams

-Figure 2: Schematic Phase diagram of the Falicov-Kimball model in dimensions greater than two. At low temperature the classical fermions (spins) settle into an ordered charge density wave state (antiferromagnetic state). The schematic diagram for the Hubbard model is the same. Reproduced from  [9,13]. - +Figure 2: Schematic Phase diagram of the Falicov-Kimball model in dimensions two or more At low temperature the classical fermions (spins) settle into an ordered charge density wave state (antiferromagnetic state). The schematic diagram for the Hubbard model is the same. Reproduced from  [9,13]. +

In dimensions greater than one, the FK model exhibits a phase transition at some \(U\) dependent critical temperature \(T_c(U)\) to a low temperature ordered phase  [14]. In terms of the heavy electrons this corresponds to them occupying only one of the two sublattices A and B, known as a Charge Density Wave (CDW) phase. In terms of spins, this is an antiferromagnetic phase.

In the disordered region above \(T_c(U)\), there are two insulating phases. For weak interactions \(U << t\), thermal fluctuations in the spins act as an effective disorder potential for the fermions. This causes them to localise, giving rise to an Anderson insulating (AI) phase  [15] which we will discuss more in section 2.3. For strong interactions \(U >> t\), the spins are not ordered. Nevertheless, their interaction with the electrons opens a gap, leading to a Mott insulator analogous to that of the Hubbard model  [16]. The presence of an interaction driven phase like the Mott insulator in an exactly solvable model is part of what makes the FK model such an interesting system.

-

By contrast, in the 1D FK model there is no Finite-Temperature Phase Transition (FTPT) to an ordered CDW phase  [17]. Indeed, dimensionality is crucial for the physics of both localisation and FTPTs. In 1D, disorder generally dominates: even the weakest disorder exponentially localises all single particle eigenstates. In the 1D FK model, this means that the whole spectrum is localised at all finite temperatures  [18–20]. Although at low temperatures, the localisation length may be so large that the states appear extended in finite sized systems  [13]. Only longer-range correlations of the disorder potential can potentially induce localisation-delocalisation transitions in 1D  [21–23]

+

By contrast, in the 1D FK model there is no Finite-Temperature Phase Transition (FTPT) to an ordered CDW phase  [17]. Indeed, dimensionality is crucial for the physics of both localisation and FTPTs. In 1D, disorder generally dominates: even the weakest disorder exponentially localises all single particle eigenstates. In the 1D FK model, this means that the whole spectrum is localised at all finite temperatures  [18–20]. Although at low temperatures, the localisation length may be so large that the states appear extended in finite sized systems  [13]. Only longer-range correlations of the disorder potential can potentially induce localisation-delocalisation transitions in 1D  [21–23].

The absence of finite temperature ordered phases in 1D systems is a general feature. It can be understood as a consequence of the fact that domain walls are energetically cheap in 1D. Thermodynamically, short-range interactions just cannot overcome the entropy of thermal defects in 1D. However, the addition of longer range interactions can overcome this  [24,25].

The absence of an FTPT in the short ranged FK chain is far from obvious because the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction mediated by the fermions  [26–29] decays as \(r^{-1}\) in 1D  [30]. This could, in principle, induce the necessary long-range interactions for the classical Ising background to order at low temperatures  [24,31]. However, Kennedy and Lieb established rigorously that, at half-filling, a CDW phase only exists at \(T = 0\) for the 1D FK model  [25].

The 1D FK model has been studied numerically, perturbatively in interaction strength \(U\) and in the continuum limit  [32]. The main results are that for attractive \(U > U_c\) the system forms electron spin bound state ‘atoms’ which repel one another  [33] and that the ground state phase diagram has a fractal structure as a function of electron filling, a devil’s staircase  [34,35].

@@ -117,7 +117,7 @@ H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) -\

Long-Ranged Ising model

The suppression of phase transitions is a common phenomenon in 1D systems and the Ising model serves as the canonical illustration of this. In terms of classical spins \(S_i = \pm 1\) the standard Ising model reads

-

\[H_{\mathrm{I}} = \sum_{\langle ij \rangle} S_i S_j\]

+

\[H_{\mathrm{I}} = \sum_{\langle ij \rangle} S_i S_j.\]

Like the FK model, the Ising model shows an FTPT to an ordered state only in 2D and above. This can be understood via Peierls’ argument  [24,25] to be a consequence of the low energy penalty for domain walls in 1D systems.

Figure 3: Domain walls in the 1D Ising model cost finite energy because they affect only one interaction. In the Long-Range Ising (LRI) model it depends on how the interactions decay with distance. @@ -125,16 +125,16 @@ H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) -\

Following Peierls’ argument, consider the difference in free energy \(\Delta F = \Delta E - T\Delta S\) between an ordered state and a state with single domain wall as in fig. 3. If this value is negative, it implies that the ordered state is unstable with respect to domain wall defects, and they will thus proliferate, destroying the ordered phase. If we consider the scaling of the two terms with system size \(L\), we see that short range interactions produce a constant energy penalty \(\Delta E\) for a domain wall. In contrast, the number of such single domain wall states scales linearly with system size so the entropy is \(\propto \ln L\). Thus, the entropic contribution dominates (eventually) in the thermodynamic limit and no finite temperature order is possible. In 2D and above, the energy penalty of a large domain wall scales like \(L^{d-1}\) which is why they can support ordered phases. This argument does not quite apply to the FK model because of the aforementioned RKKY interaction. Instead, this argument will give us insight into how to recover an ordered phase in the 1D FK model.

-Figure 4: The thermodynamic behaviour of the long-range Ising model H_{\mathrm{LRI}} = J \sum_{i\neq j} |i - j|^{-\alpha} S_i S_j as the exponent of the interaction \alpha is varied. In my simulations I stick to a value of \alpha = \tfrac{5}{4} the complexity of non-universal critical exponents. - +Figure 4: The thermodynamic behaviour of the long-range Ising model H_{\mathrm{LRI}} = J \sum_{i\neq j} |i - j|^{-\alpha} S_i S_j as the exponent of the interaction \alpha is varied. In my simulations I stick to a value of \alpha = \tfrac{5}{4} to avoid the complexity of non-universal critical exponents that arise above \alpha = \tfrac{3}{2} +

In contrast, the LRI model \(H_{\mathrm{LRI}}\) can have an FTPT in 1D.

-

\[H_{\mathrm{LRI}} = \sum_{ij} J(|i-j|) S_i S_j = J \sum_{i\neq j} |i - j|^{-\alpha} S_i S_j\]

-

Renormalisation group analyses show that the LRI model has an ordered phase in 1D for \(1 < \alpha < 2\)   [36]. Peierls’ argument can be extended  [31] to long-range interactions to provide intuition for why this is the case. Let’s consider again the energy difference between the ordered state \(|\ldots\uparrow\uparrow\uparrow\uparrow\ldots\rangle\) and a domain wall state \(|\ldots\uparrow\uparrow\downarrow\downarrow\ldots\rangle\). In the case of the LRI model, careful counting shows that this energy penalty is \[\Delta E \propto \sum_{n=1}^{\infty} n J(n)\qquad{(1)}\]

+

\[H_{\mathrm{LRI}} = \sum_{ij} J(|i-j|) S_i S_j = J \sum_{i\neq j} |i - j|^{-\alpha} S_i S_j.\]

+

Renormalisation group analyses show that the LRI model has an ordered phase in 1D for \(1 < \alpha < 2\)   [36]. Peierls’ argument can be extended  [31] to long-range interactions to provide intuition for why this is the case. Let’s consider again the energy difference between the ordered state \(|\ldots\uparrow\uparrow\uparrow\uparrow\ldots\rangle\) and a domain wall state \(|\ldots\uparrow\uparrow\downarrow\downarrow\ldots\rangle\). In the case of the LRI model, careful counting shows that this energy penalty is \[\Delta E \propto \sum_{n=1}^{\infty} n J(n),\qquad{(1)}\]

because each interaction between spins separated across the domain by a bond length \(n\) can be drawn between \(n\) equivalent pairs of sites. The behaviour then depends crucially on how eq. 1 scales with system size. Ruelle proved rigorously for a very general class of 1D systems that if \(\Delta E\) or its many-body generalisation converges to a constant in the thermodynamic limit then the free energy is analytic  [37]. This rules out a finite order phase transition, though not one of the Kosterlitz-Thouless type. Dyson also proved this, though with a slightly different condition on \(J(n)\)  [36].

-

With a power law form for \(J(n)\), there are a few cases to consider: For \(\alpha = 0\) i.e infinite range interactions, the Ising model is exactly solvable and mean field theory is exact  [38]. This limit is the same as the infinite dimensional limit. For \(\alpha \leq 1\) we have very slowly decaying interactions. \(\Delta E\) does not converge as a function of system size so the Hamiltonian is non-extensive, a topic not without some considerable controversy  [39–41] that I will not consider further here. For \(1 < \alpha < 2\), we get a phase transition to an ordered state at a finite temperature. This is what we want! At the special point \(\alpha = 2\), the energy of domain walls diverges logarithmically. This turns out to be a Kostelitz-Thouless transition  [31]. Finally, for \(\alpha > 2\) we have very quickly decaying interactions and domain walls again have a finite energy penalty. Hence, Peirels’ argument holds and there is no phase transition.

+

With a power law form for \(J(n)\), there are a few cases to consider: For \(\alpha = 0\) i.e., infinite range interactions, the Ising model is exactly solvable and mean field theory is exact  [38]. This limit is the same as the infinite dimensional limit. For \(\alpha \leq 1\) we have very slowly decaying interactions. \(\Delta E\) does not converge as a function of system size so the Hamiltonian is non-extensive, a topic not without some considerable controversy  [39–41] that I will not consider further here. For \(1 < \alpha < 2\), we get a phase transition to an ordered state at a finite temperature. This is what we want! At the special point \(\alpha = 2\), the energy of domain walls diverges logarithmically. This turns out to be a Kostelitz-Thouless transition  [31]. Finally, for \(\alpha > 2\) we have very quickly decaying interactions and domain walls again have a finite energy penalty. Hence, Peirels’ argument holds and there is no phase transition.

One final complexity is that for \(\tfrac{3}{2} < \alpha < 2\) renormalisation group methods show that the critical point has non-universal critical exponents that depend on \(\alpha\)   [42]. To avoid this potential confounding factor we will park ourselves at \(\alpha = 1.25\) when we apply these ideas to the FK model.

-

Were we to extend this to arbitrary dimension \(d\), we would find that thermodynamics properties, generally both \(d\) and \(\alpha\), long-range interactions can modify the ‘effective dimension’ of thermodynamic systems  [43].

+

Were we to extend this to arbitrary dimension \(d\), we would find that, in general, both \(d\) and \(\alpha\) affect the thermodynamic properties of the model. Long-range interactions essentially modify the ‘effective dimension’ of thermodynamic system  [43].

Next Section: The Kitaev Honeycomb Model

diff --git a/_thesis/2_Background/2.2_HKM_Model.html b/_thesis/2_Background/2.2_HKM_Model.html index 2f1b41f..881c1d4 100644 --- a/_thesis/2_Background/2.2_HKM_Model.html +++ b/_thesis/2_Background/2.2_HKM_Model.html @@ -84,15 +84,15 @@ image:

The Kitaev Honeycomb Model

-Figure 1: (a) The Kitaev honeycomb model is defined on a honeycomb lattice. The special feature of the honeycomb lattice that makes the model solvable is that each vertex is joined by exactly three bonds, i.e. the lattice is trivalent. One of three labels is assigned to each bond \{x,\;y,\;z\} here represented by colour. (b) After transforming to the Majorana representation we get an emergent gauge degree of freedom u_{jk} = \pm 1 that lives on each bond, the bond variables. These are antisymmetric, u_{jk} = -u_{kj}, so we represent them graphically with arrows on each bond that point in the direction that u_{jk} = +1. (c) The Majorana transformation can be visualised as breaking each spin into four Majoranas b_i^x,\;b_i^y,\;b_i^z,\;c_i. The x, y and z Majoranas then pair along the bonds forming conserved \mathbb{Z}_2 bond operators u_{jk} = \langle i b_i^\alpha b_j^\alpha \rangle. The remaining c_i operators form an effective quadratic Hamiltonian H = \frac{i}{4} \sum_{\langle i,j\rangle_\alpha} 2J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j. - +Figure 1: (a) The Kitaev honeycomb model is defined on a honeycomb lattice. The special feature of the honeycomb lattice that makes the model solvable is that each vertex is joined by exactly three bonds, i.e., the lattice is trivalent. One of three labels is assigned to each bond \{x,\;y,\;z\} here represented by colour. (b) After transforming to the Majorana representation we get an emergent gauge degree of freedom u_{jk} = \pm 1 that lives on each bond, the bond variables. These are antisymmetric, u_{jk} = -u_{kj}, so we represent them graphically with arrows on each bond that point in the direction that u_{jk} = +1. (c) The Majorana transformation, discussed later in the main text, can be visualised as breaking each spin into four Majoranas b_i^x,\;b_i^y,\;b_i^z,\;c_i. The b_i^x,\;b_i^y and b_i^z Majoranas then pair along the bonds forming conserved \mathbb{Z}_2 bond operators u_{jk} = \langle i b_i^\alpha b_j^\alpha \rangle. The remaining c_i operators form an effective quadratic Hamiltonian H = \frac{i}{4} \sum_{\langle i,j\rangle_\alpha} 2J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j. +

The Spin Hamiltonian

The Kitaev Honeycomb (KH) model is an exactly solvable model of interacting spin\(-1/2\) spins on the vertices of a honeycomb lattice. Each bond in the lattice is assigned a label \(\alpha \in \{ x, y, z\}\) and couples two spins along the \(\alpha\) axis. See fig. 1 for a diagram of the setup.

This gives us the Hamiltonian \[ H = - \sum_{\langle j,k\rangle_\alpha} J^{\alpha}\sigma_j^{\alpha}\sigma_k^{\alpha}, \qquad{(1)}\] where \(\sigma^\alpha_j\) is the \(\alpha\) component of a Pauli matrix acting on site \(j\) and \(\langle j,k\rangle_\alpha\) is a pair of nearest-neighbour indices connected by an \(\alpha\)-bond with exchange coupling \(J^\alpha\). Kitaev introduced this model in his seminal 2006 paper  [1].

-

The KH can arise as the result of strong spin-orbit couplings in, for example, the transition metal based compounds  [2–6]. The model is highly frustrated: each spin would like to align along a different direction with each of its three neighbours but this cannot be achieved even classically  [7,8]. This frustration leads the model to have a Quantum Sping Liquid (QSL) ground state, a complex many-body state with a high degree of entanglement but no long-range magnetic order even at zero temperature. While the possibility of a QSL ground state was suggested much earlier  [9], the KH model was the first exactly solvable models of the QSL state. The KH model has a rich ground state phase diagram with gapless and gapped phases, the latter supporting fractionalised quasiparticles with both Abelian and non-Abelian quasiparticle excitations. Anyons have been the subject of much attention because, among other reasons, they can be braided through spacetime to achieve noise tolerant quantum computations  [10]. At finite temperature the KH model undergoes a phase transition to a thermal metal state  [11]. The KH model can be solved exactly via a mapping to Majorana fermions. This mapping yields an extensive number of static \(\mathbb Z_2\) fluxes tied to an emergent gauge field with the remaining fermions are governed by a free fermion hamiltonian.

-

This section will go over the standard model in detail, first discussing the spin model, then detailing the transformation to a Majorana hamiltonian that allows a full solution while enlarging the Hamiltonian. I will discuss the properties of the emergent gauge fields and the projector. The next section will discuss anyons, topology and the Chern number, using the KH model as an explicit example.I will then discuss the ground state found via Lieb’s theorem as well as work on generalisations of the ground state to other lattices. Finally I will present the phase diagram.

+

The KH can arise as the result of strong spin-orbit couplings in, for example, the transition metal based compounds  [2–6]. The model is highly frustrated: each spin would like to align along a different direction with each of its three neighbours but this cannot be achieved even classically  [7,8]. This frustration leads the model to have a Quantum Spin Liquid (QSL) ground state, a complex many-body state with a high degree of entanglement but no long-range magnetic order even at zero temperature. While the possibility of a QSL ground state was suggested much earlier  [9], the KH model was the first exactly solvable models of the QSL state. The KH model has a rich ground state phase diagram with gapless and gapped phases, the latter supporting fractionalised quasiparticles with both Abelian and non-Abelian quasiparticle excitations. Anyons have been the subject of much attention because, among other reasons, they can be braided through spacetime to achieve noise tolerant quantum computations  [10]. At finite temperature the KH model undergoes a phase transition to a thermal metal state  [11]. The KH model can be solved exactly via a mapping to Majorana fermions. This mapping yields an extensive number of static \(\mathbb Z_2\) fluxes tied to an emergent gauge field with the remaining fermions are governed by a free fermion hamiltonian.

+

This section will go over the standard model in detail, first discussing the spin model, then detailing the transformation to a Majorana hamiltonian that allows a full solution while enlarging the Hamiltonian. I will discuss the properties of the emergent gauge fields and the projector. The next section will discuss anyons, topology and the Chern number, using the KH model as an explicit example. I will then discuss the ground state found via Lieb’s theorem as well as work on generalisations of the ground state to other lattices. Finally, I will present the phase diagram.

The Spin Model

@@ -103,7 +103,7 @@ image:

As discussed in the introduction, spin hamiltonians like that of the KH model arise in electronic systems as the result the balance of multiple effects  [5]. For instance, in certain transition metal systems with \(d^5\) valence electrons, crystal field and spin-orbit couplings conspire to shift and split the \(d\) orbitals into moments with spin \(j = 1/2\) and \(j = 3/2\). Of these, the bandwidth \(t\) of the \(j= 1/2\) band is small, meaning that even relatively meagre electron correlations (such as those induced by the \(U\) term in the Hubbard model) can lead to the opening of a Mott gap. From there we have a \(j = 1/2\) Mott insulator whose effective spin-spin interactions are again shaped by the lattice geometry and spin-orbit coupling leading some materials to have strong bond-directional Ising-type interactions  [12,13]. In the KH model the bond directionality refers to the fact that the coupling axis \(\alpha\) in terms like \(\sigma_j^{\alpha}\sigma_k^{\alpha}\) is strongly bond dependent.

In the spin hamiltonian eq. 1, we can already tease out a set of conserved fluxes that will be key to the model’s solution. These fluxes are the expectations of Wilson loop operators

\[\hat{W}_p = \prod_{\langle i,j\rangle_\alpha\; \in\; p} \sigma_i^{\alpha}\sigma_j^{\alpha},\]

-

the products of bonds winding around a closed path \(p\) on the lattice. These operators commute with the Hamiltonian and so have no time dynamics. The winding direction does not matter so long as it is fixed. By convention we will always use clockwise. Each closed path on the lattice is associated with a flux. The number of conserved quantities grows linearly with system size and is thus extensive. This is a common property for exactly solvable systems and can be compared to the heavy electrons present in the Falicov-Kimball model. The square of two loop operators is one so any contractible loop can be expressed as a product of loops around plaquettes of the lattice, as in fig. 3. For the honeycomb lattice, the plaquettes are the hexagons. The expectations of \(\hat{W}_p\) through each plaquette, the fluxes, are therefore enough to describe the whole flux sector. We will focus on these fluxes, denoting them by \(\phi_i\). Once we have made the mapping to the Majorana Hamiltonian, I will explain how these fluxes can be connected to an emergent \(B\) field which makes their interpretation as fluxes clear.

+

the products of bonds winding around a closed path \(p\) on the lattice. These operators commute with the Hamiltonian and so have no time dynamics. The winding direction does not matter so long as it is fixed. By convention, we will always use clockwise. Each closed path on the lattice is associated with a flux. The number of conserved quantities grows linearly with system size and is thus extensive. This is a common property for exactly solvable systems and can be compared to the heavy electrons present in the Falicov-Kimball model. The square of two loop operators is one so any contractible loop can be expressed as a product of loops around plaquettes of the lattice, as in fig. 3. For the honeycomb lattice, the plaquettes are the hexagons. The expectations of \(\hat{W}_p\) through each plaquette, the fluxes, are therefore enough to describe the whole flux sector. We will focus on these fluxes, denoting them by \(\phi_i\). Once we have made the mapping to the Majorana Hamiltonian, I will explain how these fluxes can be connected to an emergent \(B\) field which makes their interpretation as fluxes clear.

Figure 3: In the KH model, Wilson loop operators \hat{W}_p = \prod_{\langle i,j\rangle_\alpha\; \in\; p} \sigma_i^{\alpha}\sigma_j^{\alpha} can be composed via multiplication to produce arbitrary contractible loops. As a consequence, we need only to keep track of the value of the flux through each plaquette \phi_i. This relationship between the u_{ij} around a region and the fluxes inside it is evocative of Stokes’ theorem from classical electromagnetism. Indeed, it turns out to be closely related as we shall see later. @@ -116,40 +116,40 @@ image:

Majorana fermions are something like ‘half of a complex fermion’ and are their own antiparticle. From a set of \(N\) fermionic creation \(f_i^\dagger\) and anhilation \(f_i\) operators, we can construct \(2N\) Majorana operators \(c_m\). We can do this construction in multiple ways subject to only mild constraints required to keep the overall commutations relations correct  [1]. Majorana operators square to one, but otherwise have standard fermionic anti-commutation relations.

\(N\) spins can be mapped to \(N\) fermions with the well-known Jordan-Wigner transformation and indeed this approach can be used to solve the Kitaev model  [15]. Here I will introduce the method Kitaev used in the original paper as this forms the basis for the results that will be presented in this thesis. Rather than mapping to \(N\) fermions, Kitaev maps to \(4N\) Majoranas, effectively \(2N\) fermions. In contrast to the Jordan-Wigner transformation which makes fermions out of strings of spin operators in order to correctly produce fermionic commutation relations, the Kitaev transformation maps each spin locally to four Majoranas. The downside is that this enlarges the Hilbert space from \(2^N\) to \(4^N\). We will have to employ a projector \(\hat{P}\) to come back down to the physical Hilbert space later. As everything is local, I will drop the site indices \(ijk\) in expressions that refer to only a single site.

The mapping is defined in terms of four Majoranas per site \(b_i^x,\;b_i^y,\;b_i^z,\;c_i\) such that

-

\[\tilde{\sigma}^x = i b^x c,\; \tilde{\sigma}^y = i b^y c,\; \tilde{\sigma}^z = i b^z c\qquad{(2)}\]

-

The tildes on the spin operators \(\tilde{\sigma_i^\alpha}\) emphasise that they live in this new extended Hilbert space and are only equivalent to the original spin operators after applying a projector \(\hat{P}\). The form of the projection operator can be understood in a few ways. From a group-theoretic perspective, before projection, the operators \(\{\tilde{\sigma}^x, \tilde{\sigma}^y, \tilde{\sigma}^z\}\) form a representation of the gamma group \(G_{3,0}\). The gamma groups \(G_{p,q}\) have \(p\) generators that square to the identity and \(q\) that square (roughly) to \(-1\). The generators otherwise obey standard anticommutation relations. The well known gamma matrices \(\{\gamma^0, \gamma^1, \gamma^2, \gamma^3\}\) represent \(G_{1,3}\) the quaternions \(G_{0,3}\) and the Pauli matrices \(G_{3,0}\).

-

The Pauli matrices, however, have the additional property that the chiral element \(\sigma^x \sigma^y \sigma^z = \pm i\) is not fully determined by the group properties of \(G_{3,0}\), but it is equal to \(i\) in the Pauli matrices. Therefore, to fully reproduce the algebra of the Pauli matrices, we must project into the subspace where \(\tilde{\sigma}^x \tilde{\sigma}^y \tilde{\sigma}^z = +i\). The chiral element of the gamma matrices for instance \(\gamma_5 = i\gamma^0 \gamma^1 \gamma^2 \gamma^3\) is of central importance in quantum field theory. See  [16] for more discussion of this group theoretic view.

+

\[\tilde{\sigma}^x = i b^x c,\; \tilde{\sigma}^y = i b^y c,\; \tilde{\sigma}^z = i b^z c.\qquad{(2)}\]

+

The tildes on the spin operators \(\tilde{\sigma_i}^\alpha\) emphasise that they live in this new extended Hilbert space and are only equivalent to the original spin operators after applying a projector \(\hat{P}\). The form of the projection operator can be understood in a few ways. From a group-theoretic perspective, before projection, the operators \(\{\tilde{\sigma}^x, \tilde{\sigma}^y, \tilde{\sigma}^z\}\) form a representation of the gamma group \(G_{3,0}\). The gamma groups \(G_{p,q}\) have \(p\) generators that square to the identity and \(q\) that square (roughly) to \(-1\). The generators otherwise obey standard anticommutation relations. The well-known gamma matrices \(\{\gamma^0, \gamma^1, \gamma^2, \gamma^3\}\) represent \(G_{1,3}\) the quaternions \(G_{0,3}\) and the Pauli matrices \(G_{3,0}\).

+

The Pauli matrices, however, have the additional property that the chiral element \(\sigma^x \sigma^y \sigma^z = \pm i\) is not fully determined by the group properties of \(G_{3,0}\), but it is equal to \(i\) in the Pauli matrices. Therefore, to fully reproduce the algebra of the Pauli matrices, we must project into the subspace where \(\tilde{\sigma}^x \tilde{\sigma}^y \tilde{\sigma}^z = +i\). The chiral element of the gamma matrices for instance \(\gamma_5 = i\gamma^0 \gamma^1 \gamma^2 \gamma^3\) is of central importance in quantum field theory. See ref.  [16] for more discussion of this group theoretic view.

So the projector must project onto the subspace where \(\tilde \sigma^x \tilde \sigma^y \tilde \sigma^z = i\). If we work this through, we find that in general \(\tilde \sigma^x \tilde \sigma^y \tilde\sigma^z = iD\) where \(D = b^x b^y b^z c\) must be the identity for every site. In other words, we can only work with physical states \(|\phi\rangle\) that satisfy \(D_i|\phi\rangle = |\phi\rangle\) for all sites \(i\). From this we construct an on-site projector \(P_i = \frac{1 + D_i}{2}\) and the overall projector is simply \(P = \prod_i P_i\).

Another way to see what this is doing physically is to explicitly construct the two intermediate fermionic operators \(f\) and \(g\) that give rise to these four Majoranas. Denoting a fermion state by \(|n_f, n_g\rangle\) the Hilbert space is the set \(\{|00\rangle,|01\rangle,|10\rangle,|11\rangle\}\). We can map these to Majoranas with, for example, this definition

\[\begin{aligned} b^x = (f + f^\dagger),\;\;& b^y = -i(f - f^\dagger),\\ -b^z = (g + g^\dagger),\;\;& c = -i(g - g^\dagger), +b^z = (g + g^\dagger),\;\;& c = -i(g - g^\dagger). \end{aligned}\]

-

Working through the algebra, we see that the operator \(D = b^x b^y b^z c\) is equal to the fermion parity \(D = -(2n_f - 1)(2n_g - 1) = \pm1\) where \(n_f,\; n_g\) are the number operators. So setting \(D = 1\) everywhere is equivalent to restricting to the \(\{|01\rangle,|10\rangle\}\) though we could equally well have used \(D = -1\).

-

Expanding the product \(\prod_i P_i\) out, we find that the projector corresponds to a symmetrisation over \(\{u_{ij}\}\) states within a flux sector and overall fermion parity \(\prod_i D_i\), see  [17] or appendix A.6 for the full derivation. The significance of this is that an arbitrary many-body state can be made to have non-zero overlap with the physical subspace via the addition or removal of just a single fermion. This implies that, in the thermodynamic limit, the projection step is not generally necessary to extract physical results

+

Working through the algebra, we see that the operator \(D = b^x b^y b^z c\) is equal to the fermion parity \(D = -(2n_f - 1)(2n_g - 1) = \pm1\) where \(n_f,\; n_g\) are the number operators. So setting \(D = 1\) everywhere is equivalent to restricting to the states \(\{|01\rangle\) and \(|10\rangle\}\), though we could equally well have used \(D = -1\).

+

Expanding the product \(\prod_i P_i\) out, we find that the projector corresponds to a symmetrisation over \(\{u_{ij}\}\) states within a flux sector and overall fermion parity \(\prod_i D_i\), see ref.  [17] or appendix A.6 for the full derivation. The significance of this is that an arbitrary many-body state can be made to have non-zero overlap with the physical subspace via the addition or removal of just a single fermion. This implies that, in the thermodynamic limit, the projection step is not generally necessary to extract physical results

We can now rewrite the spin hamiltonian in Majorana form with the caveat that they are only strictly equivalent after projection. The Ising interactions \(\sigma_j^{\alpha}\sigma_k^{\alpha}\) decouple into the form \(-i (i b^\alpha_i b^\alpha_j) c_i c_j\). We factor out the bond operators \(\hat{u}_{ij} = i b^\alpha_i b^\alpha_j\) which are Hermitian and, remarkably, commute with the Hamiltonian and each other.

\[\begin{aligned} \tilde{H} &= - \sum_{\langle i,j\rangle_\alpha} J^{\alpha}\tilde{\sigma}_i^{\alpha}\tilde{\sigma}_j^{\alpha}\\ - &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} \hat{u}_{ij} \hat{c}_i \hat{c}_j + &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} \hat{u}_{ij} \hat{c}_i \hat{c}_j. \end{aligned}\]

-

The bond operators \(\hat{u}_{ij}\) square to one so have eigenvalues \(\pm1\). As they’re conserved we will work in their eigenbasis and take off the hats in the Hamiltonian.

+

The bond operators \(\hat{u}_{ij}\) square to one so have eigenvalues \(\pm1\). As they are conserved we will work in their eigenbasis and take off the hats in the Hamiltonian.

\[\begin{aligned} -H &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j +H &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c}_j. \end{aligned}\qquad{(3)}\]

The Fermion Problem

We now have a quadratic Hamiltonian, eq. 3, coupled to a classical field \(u_{ij}\). What follows is relatively standard theory for quadratic Hamiltonians  [18].

Because of the antisymmetry \(J^{\alpha} u_{ij}\), the eigenvalues of eq. 3 come in pairs \(\pm \epsilon_m\). We organise the eigenmodes of \(H\) into pairs, such that \(b_m\) and \(b_m'\) have energies \(\epsilon_m\) and \(-\epsilon_m\). The transformation \(Q\)

-

\[(c_1, c_2... c_{2N}) Q = (b_1, b_1', b_2, b_2' ... b_{N}, b_{N}')\]

+

\[(c_1, c_2... c_{2N}) Q = (b_1, b_1', b_2, b_2' ... b_{N}, b_{N}'),\]

puts the Hamiltonian into normal mode form

\[H = \frac{i}{2} \sum_m \epsilon_m b_m b_m'.\]

The determinant of \(Q\) appears when evaluating the projector explicitly, otherwise, the \(b_m\) are merely an intermediate step. From them, we form fermionic operators

-

\[ f_i = \tfrac{1}{2} (b_m + ib_m')\]

+

\[ f_i = \tfrac{1}{2} (b_m + ib_m'),\]

with their associated number operators \(n_i = f^\dagger_i f_i\). These let us write the Hamiltonian neatly as

\[ H = \sum_m \epsilon_m (n_m - \tfrac{1}{2}).\]

The energy of the ground state \(|n_m = 0\rangle\) of the many-body system at fixed \(\{u_{ij}\}\) is

-

\[E_{0} = -\frac{1}{2}\sum_m \epsilon_m \]

+

\[E_{0} = -\frac{1}{2}\sum_m \epsilon_m, \]

and we can construct any state from a particular choice of \(n_m = 0,1\). If we only care about the ground state energy \(E_{0}\), it is possible to skip forming the fermionic operators. The eigenvalues obtained directly from diagonalising \(J^{\alpha} u_{ij}\) come in \(\pm \epsilon_m\) pairs. We can take half the absolute value of the set to recover \(\sum_m \epsilon_m\) directly.

@@ -160,10 +160,10 @@ H &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c} Figure 4: A honeycomb lattice with edges in grey, along with its dual, the triangle lattice in red. The vertices of the dual lattice are the faces of the original lattice and, hence, are the locations of the vortices. (Left) The action of the gauge operator D_j at a vertex is to flip the value of the three u_{jk} variables (black lines) surrounding site j. The corresponding edges of the dual lattice (red lines) form a closed triangle. (Middle) Composing multiple adjacent D_j operators produces a large, closed dual loop or multiple disconnected dual loops. Dual loops are not directed like Wilson loops. (Right) A non-contractable loop which cannot be produced by composing D_j operators. All three operators can be thought of as the action of a vortex-vortex pair that is created, one of them is transported around the loop, and then the two annihilate again. Note that every plaquette has an even number of u_{ij}s flipped on its edge. Therefore, all retain the same flux \phi_i. -

In addition, the bond operators form a highly degenerate description of the system. The operators \(D_i = b^x_i b^y_i b^z_i c_i\) commute with \(H\) forming a set of local symmetries. The action of \(D_i\) on a state is to flip the values of the three \(u_{ij}\) bonds that connect to site \(i\). This changes the bond configuration \(\{u_{ij}\}\) but leaves the flux configuration \(\{\phi_i\}\) unchanged. Physically, we interpret \(u_{ij}\) as a gauge field with a high degree of degeneracy and \(\{D_i\}\) as the set of gauge symmetries. The Majorana bond operators \(u_{ij}\) are an emergent, classical, \(\mathbb{Z}_2\) gauge field! The flux configuration \(\{\phi_i\}\) is what encodes physical information about the system without all the gauge degeneracy.

+

In addition, the bond operators form a highly degenerate description of the system. The operators \(D_i = b^x_i b^y_i b^z_i c_i\) commute with \(H\) forming a set of local symmetries. The action of \(D_i\) on a state is to flip the values of the three \(u_{ij}\) bonds that connect to site \(i\). This changes the bond configuration \(\{u_{ij}\}\) but leaves the flux configuration \(\{\phi_i\}\) unchanged. Physically, we interpret \(u_{ij}\) as a gauge field with a high degree of degeneracy and \(\{D_i\}\) as the set of gauge symmetries. Products of the gauge symmetries correspond to closed loops on the dual lattice, see fig. 4. The Majorana bond operators \(u_{ij}\) are an emergent, classical, \(\mathbb{Z}_2\) gauge field! The flux configuration \(\{\phi_i\}\) is what encodes physical information about the system without all the gauge degeneracy.

The ground state of the KH model is the flux configuration where all fluxes are one \(\{\phi_i = +1\; \forall \; i\}\). This can be proven via Lieb’s theorem  [19] which gives the lowest energy magnetic flux configuration for a system of electrons hopping in a magnetic field. Kitaev remarks in his original paper that he was not initially aware of the relevance of Lieb’s 1994 result. This is not surprising because at first glance the two models seem quite different. Yet the connection is quite instructive for understanding the KH and its generalisations.

Lieb discussed a model of mobile electrons

-

\[H = \sum_{ij} t_{ij} c^\dagger_i c_j\]

+

\[H = \sum_{ij} t_{ij} c^\dagger_i c_j,\]

where the hopping terms \(t_{ij} = |t_{ij}|\exp(i\theta_{ij})\) incorporate Aharanhov-Bohm (AB) phases  [20] \(\theta_{ij}\). The AB phases model the effect of a slowly varying magnetic field on the electrons through the integral of the magnetic vector potential

\[\theta_{ij} = \int_i^j \vec{A} \cdot d\vec{l},\]

a Peierls substitution  [21]. If we map the Majorana form of the Kitaev model to Lieb’s model, we see that our \(t_{ij} = i J^\alpha u_{ij}\). The \(i u_{ij} = \pm i\) correspond to AB phases \(\theta_{ij} = \pi/2\) or \(3\pi/2\) along each bond.

@@ -177,7 +177,7 @@ H &= i \sum_{\langle i,j\rangle_\alpha} J^{\alpha} u_{ij} \hat{c}_i \hat{c} &= \prod_{\mathcal{P}_i} \exp{(i\theta_{jk}})\\ &= \exp \left( i \sum_{\mathcal{P}_i} \theta_{jk} \right)\\ &= \exp \left( i \oint_{\mathcal{P}_i} \vec{A} \cdot d\vec{l} \right)\\ - &= \exp \left( i Q_i \right) + &= \exp \left( i Q_i \right). \end{aligned}\qquad{(5)}\]

Thus, we can interpret the fluxes \(\phi_i\) as the exponential of magnetic fluxes \(Q_m\) of some fictitious gauge field \(\vec{A}\) and the bond operators as \(i u_{ij} = \exp i \int_i^j \vec{A} \cdot d\vec{l}\). In this analogy to classical electromagnetism, the sets \(\{u_{ij}\}\) that correspond to the same \(\{\phi_i\}\) are all gauge equivalent as we have already seen via other means. The fact that fluxes can be written as products of bond operators and composed is a consequence of eq. 5. If the lattice contains odd plaquettes, as in the Yao-Kivelson model  [26], the complex fluxes that appear are a sign that chiral symmetry has been broken.

In full, Lieb’s theorem states that the ground state has magnetic flux \(Q_i = \sum_{\mathcal{P}_i}\theta_{ij} = \pi \; (\mathrm{mod} \;2\pi)\) for plaquettes with \(0 \; (\mathrm{mod}\;4)\) sides and \(0 \; (\mathrm{mod}\;2\pi)\) for plaquettes with \(2 \; (\mathrm{mod}\;4)\) sides. In terms of our fluxes, this means \(\phi = -1\) for squares, \(\phi = 1\) for hexagons and so on.

@@ -195,7 +195,7 @@ A honeycomb lattice (in black) along with its dual (in red). (Left) The product Figure 6: Wilson loops that wind the major or minor diameters of the torus measure flux winding through the hole of the doughnut/torus (red) or through the filling (green). If they made doughnuts which had both had a jam filling and a hole, this analogy would be a lot easier to make  [27]. -

A final but important point to mention is that the local fluxes \(\phi_i\) are not quite all there is. We’ve seen that products of \(\phi_i\) can be used to construct the flux associated with arbitrary contractible loops. On the plane contractible loops are all there is. However, on the torus we can construct two global fluxes \(\Phi_x\) and \(\Phi_y\) which correspond to paths tracing the major and minor axes. The four sectors spanned by the \(\pm1\) values of these fluxes are gapped away from one another, but only by virtual tunnelling processes so the gap decays exponentially with system size  [1]. Physically \(\Phi_x\) and \(\Phi_y\) could be thought of as measuring the flux that threads through the hole of the doughnut. In general, surfaces with genus \(g\) have \(g\) ‘handles’ and \(2g\) of these global fluxes. At first glance, it may seem that these fluxes would not have much relevance to physical realisations of the Kitaev model which are likely to have a planar geometry. However these fluxes are closely linked to topology and the existence of anyonic quasiparticle excitations in the model, which we will discuss next.

+

A final but important point to mention is that the local fluxes \(\phi_i\) are not quite all there is. We’ve seen that products of \(\phi_i\) can be used to construct the flux associated with arbitrary contractible loops. On the plane contractible loops are all there is. However, on the torus we can construct two global fluxes \(\Phi_x\) and \(\Phi_y\) which correspond to paths tracing the major and minor axes. The four sectors spanned by the \(\pm1\) values of these fluxes are gapped away from one another, but only by virtual tunnelling processes so the gap decays exponentially with system size  [1]. Physically \(\Phi_x\) and \(\Phi_y\) could be thought of as measuring the flux that threads through the hole of the doughnut. In general, surfaces with genus \(g\) have \(g\) ‘handles’ and \(2g\) of these global fluxes. At first glance, it may seem that these fluxes would not have much relevance to physical realisations of the Kitaev model which are likely to have a planar geometry. However, these fluxes are closely linked to topology and the existence of anyonic quasiparticle excitations in the model, which we will discuss next.

Anyons, Topology and the Chern number

@@ -203,7 +203,7 @@ A honeycomb lattice (in black) along with its dual (in red). (Left) The product Figure 7: Worldlines of particles in 2D can become tangled or braided with one another. -

To discuss different ground state phases of the KH model, we must first review the topic of anyons and topology. The standard argument for the existence of Fermions and Bosons goes like this: the quantum state of a system must pick up a factor of \(\pm1\) if two identical particles are swapped. Only \(\pm1\) are allowed since swapping twice must correspond to the identity. This argument works in 3D for states without topological degeneracy, which seems to be true of the real world, but condensed matter systems are subject to no such constraints.

+

To discuss different ground state phases of the KH model, we must first review the topic of anyons and topology. The standard argument for the existence of fermions and bosons goes like this: the quantum state of a system must pick up a factor of \(\pm1\) if two identical particles are swapped. Only \(\pm1\) are allowed since swapping twice must correspond to the identity. This argument works in 3D for states without topological degeneracy, which seems to be true of the real world, but condensed matter systems are subject to no such constraints.

In gapped condensed matter systems, all equal time correlators decay exponentially with distance  [28]. Put another way, gapped systems support quasiparticles with a definite location in space and finite extent. As such it is meaningful to consider what would happen to the overall quantum state if we were to adiabatically carry out a series of swaps as described above. This is known as braiding. Recently, braiding in topological systems has attracted interest because of proposals to use ground state degeneracy to implement both passively fault tolerant and actively stabilised quantum computations  [29–31].

First we realise that, in 2D, swapping identical particles twice is not topologically equivalent to the identity, see fig. 7. Instead it corresponds to encircling one particle around the other. This means we can in general pick up any complex phase \(e^{i\theta}\) upon exchange, hence the name any-ons. Those that pick up a complex phase are known as Abelian anyons because complex multiplication commutes and hence the group of braiding operations on them forms an Abelian group.

The KH model has a topologically degenerate ground state with sectors labelled by the values of the topological fluxes \((\Phi_x\), \(\Phi_y)\). Consider the operation in which a quasiparticle pair is created from the ground state, transported around one of the non-contractible loops, and then annihilated together, call them \(\mathcal{T}_{x}\) and \(\mathcal{T}_{y}\). These operations move the system around within the ground state manifold and they need not commute. This leads to non-Abelian anyons. As Kitaev pointed out, these operations are not specific to the torus: the operation \(\mathcal{T}_{x}\mathcal{T}_{y}\mathcal{T}_{x}^{-1}\mathcal{T}_{y}^{-1}\) corresponds to an operation in which none of the particles crosses the torus, rather one simply winds around the other. Hence, these effects are relevant even for the planar case!

@@ -225,7 +225,7 @@ This all works the same way for the amorphous lattice but the diagram is a lot m

Setting the overall energy scale with the constraint \(J_x + J_y + J_z = 1\) yields a triangular phase diagram. In each of the corners one of the spin-coupling directions dominates, \(|J_\alpha > |J_\beta| + |J_\gamma|\), yielding three equivalent \(A_\alpha\) phases while the central triangle around \(J_x = J_y = J_z\) is called the B phase. Both phases support two kinds of quasiparticles, fermions and \(\mathbb{Z}_2\)-vortices. In the A phases, the vortices have bosonic statistics with respect to themselves but act like fermions with respect to the fermions, hence they are Abelian anyons. This phase has the same anyonic structure as the Toric code  [29]. The B phase can be described as a semi-metal of the Majorana fermions  [5]. Since the B phase is gapless, the quasiparticles aren’t localised and so don’t have braiding statistics.

An external magnetic can be used to break chiral symmetry. The lowest order term which breaks chiral symmetry but retains the solvability of the model is the three spin term \[ -\sum_{(i,j,k)} \sigma_i^{\alpha} \sigma_j^{\beta} \sigma_k^{\gamma} +\sum_{(i,j,k)} \sigma_i^{\alpha} \sigma_j^{\beta} \sigma_k^{\gamma}, \] where the sum \((i,j,k)\) runs over consecutive indices around plaquettes. The addition of this to the spin model leads to two bond terms in the corresponding Majorana model. The effect of breaking chiral symmetry is to open a gap in the B phase. The vortices of the gapped B phase are non-Abelian anyons. This phase has the same anyonic exchange statistics as \(p_x + ip_y\) superconductor  [38], the Moore-Read state for the \(\nu = 5/2\) fractional quantum Hall state  [39] and many other systems  [40–44]. Collectively these systems have attracted interest as possible physical realisations for quantum computers whose operations are based on braiding operations.

At finite temperatures, recent work has shown that the KH model undergoes a transition to a thermal metal phase. Vortex disorder causes the fermion gap to fill up and the DOS has a characteristic logarithmic divergence at zero energy which can be understood from random matrix theory  [11].

To surmise, the KH model is remarkable because it combines three key properties. First, the form of the Hamiltonian can plausibly be realised by a real material. Candidate materials, such as \(\alpha\mathrm{-RuCl}_3\), are known to have sufficiently strong spin-orbit coupling and the correct lattice structure to behave according to the KH model with small corrections  [5,45]. Second, its ground state is the canonical example of the long sought after QSL state, its dynamical spin-spin correlation functions are zero beyond nearest neighbour separation  [46]. Its excitations are anyons, particles that can only exist in 2D that break the normal fermion/boson dichotomy. Third, and perhaps most importantly, this model is a rare many-body interacting quantum system that can be treated analytically. It is exactly solvable. We can explicitly write down its many-body ground states in terms of single particle states  [1]. The solubility of the KH model, like the FK model, comes about because the model has extensively many conserved degrees of freedom. These conserved quantities can be factored out as classical degrees of freedom, leaving behind a non-interacting quantum model that is easy to solve.

diff --git a/_thesis/2_Background/2.4_Disorder.html b/_thesis/2_Background/2.4_Disorder.html index 430ba78..1cfde97 100644 --- a/_thesis/2_Background/2.4_Disorder.html +++ b/_thesis/2_Background/2.4_Disorder.html @@ -88,7 +88,7 @@ H = -t\sum_{\langle jk \rangle} c^\dagger_j c_k + \sum_j V_j c_j^\dagger c_j.

Localisation phenomena are strongly dimension dependent. In 3D the scaling theory of localisation  [5,6] shows that Anderson localisation is a critical phenomenon with critical exponents both for how the conductivity vanishes with energy when approaching the mobility edge and for how the localisation length increases below it. By contrast, in 1D disorder generally dominates. Even the weakest disorder exponentially localises all single particle eigenstates in the 1D Anderson model. Only long-range spatial correlations of the disorder potential can induce delocalisation  [7–12].

Later localisation was found in disordered interacting many-body systems:

\[ -H = -t\sum_{\langle jk \rangle} c^\dagger_j c_k + \sum_j V_j c_j^\dagger c_j + U\sum_{jk} n_j n_k +H = -t\sum_{\langle jk \rangle} c^\dagger_j c_k + \sum_j V_j c_j^\dagger c_j + U\sum_{jk} n_j n_k. \] Here, in contrast to the Anderson model, localisation phenomena are robust to weak perturbations of the Hamiltonian. This is called many-body localisation  [13,14].

Both many-body localisation and Anderson localisation depend crucially on the presence of quenched disorder. Quenched disorder takes the form a static background field drawn from an arbitrary probability distribution to which the model is coupled. Disorder may also be introduced into the initial state of the system rather than the Hamiltonian. This has led to ongoing interest in the possibility of disorder-free localisation where the disorder is instead annealed. In this scenario, the disorder necessary to generate localisation is generated entirely from the thermal fluctuations of the model.

The concept of disorder-free localisation was first proposed in the context of Helium mixtures  [15] and then extended to heavy-light mixtures in which multiple species with large mass ratios interact. The idea is that the heavier particles act as an effective disorder potential for the lighter ones, inducing localisation. Two such models  [16,17] instead find that the models thermalise exponentially slowly in system size, which Ref.  [16] dubs Quasi-MBL.

@@ -96,7 +96,7 @@ H = -t\sum_{\langle jk \rangle} c^\dagger_j c_k + \sum_j V_j c_j^\dagger c_j + U

In Chapter 3 we will consider a generalised FK model in 1D and study how the disorder generated near a 1D thermodynamic phase transition interacts with localisation physics.

Topological Disorder

-

So far we have considered disorder as a static or dynamic field coupled to a model defined on a translation invariant lattice. Another kind of disordered system that is worthy of study are amorphous systems. Amorphous systems have disordered bond connectivity, so called topological disorder. As discussed in the introduction these include amorphous semiconductors such as amorphous Germanium and Silicon   [20–23]. While materials do not have long-range lattice structure they can enforce local constraints such as the approximate coordination number \(z = 4\) of silicon.

+

So far we have considered disorder as a static or dynamic field coupled to a model defined on a translation invariant lattice. Another kind of disordered system that is worthy of study are amorphous systems. Amorphous systems have disordered bond connectivity, so called topological disorder. As discussed in the introduction these include amorphous semiconductors such as amorphous germanium and silicon   [20–23]. While materials do not have long-range lattice structure they can enforce local constraints such as the approximate coordination number \(z = 4\) of silicon.

Topological disorder can be qualitatively different from other disordered systems. Disordered graphs are constrained by fixed coordination number and the Euler equation. A standard method for generating such graphs with coordination number \(d+1\) is Voronoi tessellation  [24,25]. The Harris  [26] and the Imry-Mar  [27] criteria are key results on the effect of disorder on thermodynamic phase transitions. The Harris criterion signals when disorder will affect the universality of a thermodynamic critical point while the Imry-Ma criterion simply forbids the formation of long-range ordered states in \(d \leq 2\) dimensions in the presence of disorder. Both these criteria are modified for the case of topological disorder. This is because the Euler equation and vertex degree constraints lead to strong anti-correlations which mean that topological disorder is effectively weaker than standard disorder in 2D  [28,29]. This does not apply to 3D Voronoi lattices where the Euler equation contains an extra volume term and so is effectively a weaker constraint. It is worth exploring how QSLs and disorder interact. The KH model has been studied subject to both flux  [30] and bond  [31] disorder. In some instances it seems that disorder can even promote the formation of a QSL ground state  [32]. It has also been shown that the KH model exhibits disorder-free localisation after a quantum quench  [33].

In chapter 4 we will put the Kitaev model onto 2D Voronoi lattices and show that much of the rich character of the model is preserved despite the lack of long-range order.

@@ -111,18 +111,18 @@ H = -t\sum_{\langle jk \rangle} c^\dagger_j c_k + \sum_j V_j c_j^\dagger c_j + U

For low dimensional systems with quenched disorder, transfer matrix methods can be used to directly extract the localisation length. These work by turning the time independent Schrödinger equation \(\hat{H}|\psi\rangle = E|\psi\rangle\) into a matrix equation linking the amplitude of \(\psi\) on each \(d-1\) dimensional slice of the system to the next and looking at average properties of this transmission matrix. This method is less useful for systems like the FK model where the disorder as a whole must be sampled from the thermodynamic ensemble.

A more versatile method is based on the Inverse Participation Ratio (IPR). The IPR is defined for a normalised wave function \(\psi_i = \psi(x_i), \sum_i |\psi_i|^2 = 1\) as its fourth moment  [6]:

\[ -P^{-1} = \sum_i |\psi_i|^4 +P^{-1} = \sum_i |\psi_i|^4. \]

-

The name derive from the fact that this operator acts as a measure of the volume where the wavefunction is significantly different from zero. They can alternatively be thought of as providing a measure of the average diameter \(R\) from \(R = P^{1/d}\). See fig. 1 for the distinction between \(R\) and \(\lambda\).

+

The name derives from the fact that this operator acts as a measure of the volume where the wavefunction is significantly different from zero. They can alternatively be thought of as providing a measure of the average diameter \(R\) from \(R = P^{1/d}\). See fig. 1 for the distinction between \(R\) and \(\lambda\).

For localised states, the inverse participation ratio \(P^{-1}\) is independent of system size while for plane wave states in \(d\) dimensions \(P^{-1} = L^{-d}\). States may also be intermediate between localised and extended, described by their fractal dimensionality \(d > d* > 0\):

\[ -P(L)^{-1} \sim L^{-d*} +P(L)^{-1} \sim L^{-d*}. \]

Such intermediate states tend to appear as critical phenomena near mobility edges  [34]. For finite size systems, these relations only hold once the system size \(L\) is much greater than the localisation length. When the localisation length is comparable to the system size, the states still contribute to transport, this is the aforementioned weak localisation effect  [35,36].

In the following two chapters I will use an energy resolved IPR \[ \begin{aligned} DOS(\omega) &= \sum_n \delta(\omega - \epsilon_n)\\ -IPR(\omega) &= DOS(\omega)^{-1} \sum_{n,i} \delta(\omega - \epsilon_n) |\psi_{n,i}|^4 +IPR(\omega) &= DOS(\omega)^{-1} \sum_{n,i} \delta(\omega - \epsilon_n) |\psi_{n,i}|^4, \end{aligned} \] where \(\psi_{n,i}\) is the wavefunction corresponding to the energy \(\epsilon_n\) at the ith site. In practice I bin the IPRs into finely spaced bins in energy space and use the mean IPR within each bin.

diff --git a/_thesis/3_Long_Range_Falicov_Kimball/3.1_LRFK_Model.html b/_thesis/3_Long_Range_Falicov_Kimball/3.1_LRFK_Model.html index c5ee1f3..86f7dac 100644 --- a/_thesis/3_Long_Range_Falicov_Kimball/3.1_LRFK_Model.html +++ b/_thesis/3_Long_Range_Falicov_Kimball/3.1_LRFK_Model.html @@ -98,11 +98,11 @@ image:

We interpret the FK model as a model of spinless fermions, \(c^\dagger_{i}\), hopping on a 1D lattice against a classical Ising spin background, \(S_i \in {\pm \frac{1}{2}}\). The fermions couple to the spins via an onsite interaction with strength \(U\) which we supplement by a long-range interaction, \[ J_{ij} = 4\kappa J\; (-1)^{|i-j|} |i-j|^{-\alpha}, \]

-

between the spins. The additional coupling is very similar to that of the long-range Ising (LRI) model. It stabilises the Antiferromagnetic (AFM) order of the Ising spins which promotes the finite temperature CDW phase of the fermionic sector.

-

The hopping strength of the electrons, \(t = 1\), sets the overall energy scale and we concentrate throughout on the particle-hole symmetric point at zero chemical potential and half filling  [6].

+

between the spins, see fig. 1. The additional coupling is very similar to that of the long-range Ising (LRI) model. It stabilises the antiferromagnetic (AFM) order of the Ising spins which promotes the finite temperature CDW phase of the fermionic sector.

+

The hopping strength of the electrons, \(t = 1\), sets the overall energy scale and we concentrate throughout on the particle-hole symmetric point at zero chemical potential and half-filling  [6].

\[\begin{aligned} -H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) -\;t \sum_{i} (c^\dagger_{i}c_{i+1} + \textit{h.c.)}\\ -& + \sum_{i, j}^{N} J_{ij} S_i S_j +H_{\mathrm{FK}} = & \;U \sum_{i} S_i\;(c^\dagger_{i}c^{\phantom{\dagger}}_{i} - \tfrac{1}{2}) -\;t \sum_{i} (c^\dagger_{i}c^{\phantom{\dagger}}_{i+1} + \textit{h.c.)}\\ +& + \sum_{i, j}^{N} J_{ij} S_i S_j. \label{eq:HFK}\end{aligned}\]

Without proper normalisation, the long-range coupling would render the critical temperature strongly system size dependent for small system sizes. Within a mean field approximation, the critical temperature scales with the effective coupling to all the neighbours of each site, which for a system with \(N\) sites is \(\sum_{i=1}^{N} i^{-\alpha}\). Hence, the normalisation \(\kappa^{-1} = \sum_{i=1}^{N} i^{-\alpha}\), renders the critical temperature independent of system size in the mean field approximation. This greatly improves the finite size behaviour of the model.

Taking the limit \(U = 0\) decouples the spins from the fermions, which gives a spin sector governed by a classical long-range Ising model. Note, the transformation of the spins \(S_i \to (-1)^{i} S_i\) maps the AFM model to the FM one. As discussed in the background section, Peierls’ classic argument can be extended to long-range couplings to show that, for the 1D LRI model, a power law decay of \(\alpha < 2\) is required for a FTPT. This is because the energy of defect domain scales with the system size when the interactions are long-range and can overcome the entropic contribution. A renormalisation group analysis supports this finding and shows that the critical exponents are only universal for \(\alpha \leq 3/2\)  [7–9]. In the following, we choose \(\alpha = 5/4\) to avoid the additional complexity of non-universal critical points.

diff --git a/_thesis/3_Long_Range_Falicov_Kimball/3.2_LRFK_Methods.html b/_thesis/3_Long_Range_Falicov_Kimball/3.2_LRFK_Methods.html index 414f713..ba5ff9c 100644 --- a/_thesis/3_Long_Range_Falicov_Kimball/3.2_LRFK_Methods.html +++ b/_thesis/3_Long_Range_Falicov_Kimball/3.2_LRFK_Methods.html @@ -80,51 +80,52 @@ image:

Thermodynamics of the LRFK Model

Figure 1: Two Markov Chain Monte Carlo (MCMC) walks starting from the CDW state for a system with N = 100 sites and 10,000 MCMC steps but at a temperature close to but above the ordered state (left column) and much higher than it (right column). In this simulation, only a single spin can be flipped per step according to the Metropolis-Hastings Algorithm. The staggered magnetisation m = N^{-1} \sum_i (-1)^i \; S_i order parameter is plotted below. At both temperatures the thermal average of m is zero, while the initial state has m = 1. The higher temperature allows the MCMC to converge more quickly and to fluctuate about the mean with a shorter autocorrelation time. t = 1, \alpha = 1.25, T = {2.5,5}, J = U = 5 - +
-

A classical Markov Chain Monte Carlo (MCMC) method allows us to solve our LRFK model efficiently, yielding unbiased estimates of thermal expectation values.

+

A classical Markov Chain Monte Carlo (MCMC) method allows us to solve our LRFK model efficiently, yielding unbiased estimates of thermal expectation values, see fig. 1.

Since the spin configurations are classical, the LRFK Hamiltonian can be split into a classical spin part \(H_s\) and an operator valued part \(H_c\).

\[\begin{aligned} H_s& = - \frac{U}{2}S_i + \sum_{i, j}^{N} J_{ij} S_i S_j \\ -H_c& = \sum_i U S_i c^\dagger_{i}c_{i} -t(c^\dagger_{i}c_{i+1} + c^\dagger_{i+1}c_{i}) \end{aligned}\]

+H_c& = \sum_i U S_i c^\dagger_{i}c^{\phantom{\dagger}}_{i} -t(c^\dagger_{i}c^{\phantom{\dagger}}_{i+1} + c^\dagger_{i+1}c^{\phantom{\dagger}}_{i}). \end{aligned}\]

The partition function can then be written as a sum over spin configurations, \(\vec{S} = (S_0, S_1...S_{N-1})\):

\[\begin{aligned} \mathcal{Z} = \mathrm{Tr} e^{-\beta H}= \sum_{\vec{S}} e^{-\beta H_s} \mathrm{Tr}_c e^{-\beta H_c} .\end{aligned}\]

-

The contribution of \(H_c\) to the grand canonical partition function can be obtained by performing the sum over eigenstate occupation numbers giving \(-\beta F_c[\vec{S}] = \sum_k \ln{(1 + e^{- \beta \epsilon_k})}\) where \({\epsilon_k[\vec{S}]}\) are the eigenvalues of the matrix representation of \(H_c\) determined through exact diagonalisation. This gives a partition function containing a classical energy which corresponds to the long-range interaction of the spins, and a free energy which corresponds to the quantum subsystem.

+

The contribution of \(H_c\) to the grand canonical partition function can be obtained by performing the sum over eigenstate occupation numbers giving \(-\beta F_c[\vec{S}] = \sum_k \ln{(1 + e^{- \beta \epsilon_k})}\) where \({\epsilon_k[\vec{S}]}\) are the eigenvalues of the matrix representation of \(H_c\) determined through exact diagonalisation. This gives a partition function containing a classical energy which corresponds to the long-range interaction of the spins, and a free energy which corresponds to the quantum subsystem

\[\begin{aligned} -\mathcal{Z} = \sum_{\vec{S}} e^{-\beta H_S[\vec{S}] - \beta F_c[\vec{S}]} = \sum_{\vec{S}} e^{-\beta E[\vec{S}]}\end{aligned}\]

+\mathcal{Z} = \sum_{\vec{S}} e^{-\beta H_S[\vec{S}] - \beta F_c[\vec{S}]} = \sum_{\vec{S}} e^{-\beta E[\vec{S}]}.\end{aligned}\]

Markov Chain Monte Carlo and Emergent Disorder

Classical MCMC defines a weighted random walk over the spin states \((\vec{S}_0, \vec{S}_1, \vec{S}_2, ...)\), such that the likelihood of visiting a particular state converges to its Boltzmann probability \(p(\vec{S}) = \mathcal{Z}^{-1} e^{-\beta E}\). Hence, any observable can be estimated as a mean over the states visited by the walk  [4–6],

\[\begin{aligned} \label{eq:thermal_expectation} -\langle O \rangle & = \sum_{\vec{S}} p(\vec{S}) \langle O \rangle\\ - & = \sum_{i = 0}^{M} \langle O\rangle \pm \mathcal{O}(M^{-\tfrac{1}{2}}) +\langle O \rangle & = \sum_{\vec{S}} p(\vec{S}) \langle O \rangle_{\vec{S}}\\ + & = \sum_{i = 0}^{M} \langle O\rangle_{\vec{S}_i} \pm \mathcal{O}(M^{-\tfrac{1}{2}}), \end{aligned}\]

-

where the former sum runs over the entire state space while the later runs over all the state visited by a particular MCMC run.

+

where the former sum runs over the entire state space while the latter runs over all the states visited by a particular MCMC run,

\[\begin{aligned} -\langle O \rangle_{\vec{S}}& = \sum_{\nu} n_F(\epsilon_{\nu}) \langle O \rangle{\nu} +\langle O \rangle_{\vec{S}}& = \sum_{\nu} n_F(\epsilon_{\nu}) \langle O \rangle{\nu}, \end{aligned}\]

-

Where \(\nu\) runs over the eigenstates of \(H_c\) for a particular spin configuration and \(n_F(\epsilon) = \left(e^{-\beta\epsilon} + 1\right)^{-1}\) is the Fermi function.

+

where \(\nu\) runs over the eigenstates of \(H_c\) for a particular spin configuration and \(n_F(\epsilon) = \left(e^{-\beta\epsilon} + 1\right)^{-1}\) is the Fermi function.

The choice of the transition function for MCMC is under-determined as one only needs to satisfy a set of balance conditions for which there are many solutions  [7]. Here, we incorporate a modification to the standard Metropolis-Hastings algorithm  [8] gleaned from Krauth  [9].

The standard algorithm decomposes the transition probability into \(\mathcal{T}(a \to b) = p(a \to b)\mathcal{A}(a \to b)\). Here, \(p\) is the proposal distribution, that we can directly sample from, while \(\mathcal{A}\) is the acceptance probability. The standard Metropolis-Hastings choice is

\[\mathcal{A}(a \to b) = \min\left(1, \frac{p(b\to a)}{p(a\to b)} e^{-\beta \Delta E}\right)\;,\]

with \(\Delta E = E_b - E_a\). The walk then proceeds by sampling a state \(b\) from \(p\) and moving to \(b\) with probability \(\mathcal{A}(a \to b)\). The latter operation is typically implemented by performing a transition if a uniform random sample from the unit interval is less than \(\mathcal{A}(a \to b)\) and otherwise repeating the current state as the next step in the random walk. The proposal distribution is often symmetric, so it does not appear in \(\mathcal{A}\). Here, we flip a small number of sites in \(b\) at random to generate proposals, which is a symmetric proposal.

-

In our computations  [10], we employ a modification to this algorithm based on the observation that the free energy of the FK system is composed of a classical part which is much quicker to compute than the quantum part. Hence, we can obtain a computational speed up by first considering the value of the classical energy difference \(\Delta H_s\) and rejecting the transition if the former is too high. We only compute the quantum energy difference \(\Delta F_c\) if the transition is accepted. We then perform a second rejection sampling step based upon it. This corresponds to two nested comparisons with the majority of the work only occurring if the first test passes. This modified scheme has the acceptance function \[\mathcal{A}(a \to b) = \min\left(1, e^{-\beta \Delta H_s}\right)\min\left(1, e^{-\beta \Delta F_c}\right)\;.\]

+

In our computations  [10], we employ a modification to this algorithm based on the observation that the free energy of the FK system is composed of a classical part which is much quicker to compute than the quantum part. Hence, we can obtain a computational speed up by first considering the value of the classical energy difference \(\Delta H_s\) and rejecting the transition if the former is too high. We only compute the quantum energy difference \(\Delta F_c\) if the transition is accepted. We then perform a second rejection sampling step based upon it. This corresponds to two nested comparisons with the majority of the work only occurring if the first test passes. This modified scheme has the acceptance function

+

\[\mathcal{A}(a \to b) = \min\left(1, e^{-\beta \Delta H_s}\right)\min\left(1, e^{-\beta \Delta F_c}\right).\]

For the model parameters used, we find that with our new scheme the matrix diagonalisation is skipped around 30% of the time at \(T = 2.5\) and up to 80% at \(T = 1.5\). We observe that for \(N = 50\), the matrix diagonalisation, if it occurs, occupies around 60% of the total computation time for a single step. This rises to 90% at N = 300 and further increases for larger N. We therefore get the greatest speedup for large system sizes at low temperature where many prospective transitions are rejected at the classical stage and the matrix computation takes up the greatest fraction of the total computation time. The upshot is that we find a speedup of up to a factor of 10 at the cost of very little extra algorithmic complexity.

Our two-step method should be distinguished from the more common method for speeding up MCMC, which is to add asymmetry to the proposal distribution to make it as similar as possible to \(\min\left(1, e^{-\beta \Delta E}\right)\). This reduces the number of rejected states, which brings the algorithm closer in efficiency to a direct sampling method. However, it comes at the expense of requiring a way to directly sample from this complex distribution. This is a problem which MCMC was employed to solve in the first place. For example, recent work trains restricted Boltzmann machines (RBMs) to generate samples for the proposal distribution of the FK model  [11]. The RBMs are chosen as a parametrisation of the proposal distribution that can be efficiently sampled from, while offering sufficient flexibility that they can be adjusted to match the target distribution. Our proposed method is considerably simpler and does not require training while still reaping some of the benefits of reduced computation.

Scaling

-Figure 2: (Left) The order parameters, \langle m^2 \rangle(solid) and 1 - f (dashed) describing the onset of the charge density wave phase of the LRFK model at low temperature with staggered magnetisation m = N^{-1} \sum_i (-1)^i S_i and fermionic order parameter f = 2 N^{-1}|\sum_i (-1)^i \; \langle c^\dagger_{i}c_{i}| \rangle . (Right) The crossing of the Binder cumulant, B = \langle m^4 \rangle / \langle m^2 \rangle^2, with system size provides a diagnostic that the phase transition is not a finite size effect, it is used to estimate the critical lines shown in the phase diagram later. All plots use system sizes N = [10,20,30,50,70,110,160,250] and lines are coloured from N = 10 in dark blue to N = 250 in yellow. The parameter values U = 5,\;J = 5,\;\alpha = 1.25 except where explicitly mentioned. - +Figure 2: (Left) The order parameters, \langle m^2 \rangle (solid) and 1 - f (dashed) describing the onset of the charge density wave phase of the LRFK model at low temperature with staggered magnetisation m = N^{-1} \sum_i (-1)^i S_i and fermionic order parameter f = 2 N^{-1}|\sum_i (-1)^i \; \langle c^\dagger_{i}c^{\phantom{\dagger}}_{i}| \rangle . (Right) The crossing of the Binder cumulant, B = \langle m^4 \rangle / \langle m^2 \rangle^2, with system size provides a diagnostic that the phase transition is not a finite size effect, it is used to estimate the critical lines shown in the phase diagram later. All plots use system sizes N = [10,20,30,50,70,110,160,250] and lines are coloured from N = 10 in dark blue to N = 250 in yellow. The parameter values U = 5,\;J = 5,\;\alpha = 1.25 except where explicitly mentioned. +

To improve the scaling of finite size effects, we make the replacement \(|i - j|^{-\alpha} \rightarrow |f(i - j)|^{-\alpha}\), in both \(J_{ij}\) and \(\kappa\), where \(f(x) = \frac{N}{\pi}\sin \frac{\pi x}{N}\). \(f\) is smooth across the circular boundary and its effect diminished for larger systems  [12]. We only consider even system sizes given that odd system sizes are not commensurate with a CDW state.

To identify critical points, I use the Binder cumulant \(U_B\) defined by

\[ -U_B = 1 - \frac{\langle\mu_4\rangle}{3\langle\mu_2\rangle^2} +U_B = 1 - \frac{\langle\mu_4\rangle}{3\langle\mu_2\rangle^2}, \]

where \(\mu_n = \langle(m - \langle m\rangle)^n\rangle\) are the central moments of the order parameter \(m = \sum_i (-1)^i (2n_i - 1) / N\). The Binder cumulant evaluated against temperature is a diagnostic for the existence of a phase transition. If multiple such curves are plotted for different system sizes, a crossing indicates the location of a critical point while the lines do not cross for systems that don’t have a phase transition in the thermodynamic limit  [2,3].

Next Section: Results

diff --git a/_thesis/3_Long_Range_Falicov_Kimball/3.3_LRFK_Results.html b/_thesis/3_Long_Range_Falicov_Kimball/3.3_LRFK_Results.html index 3404723..bdfd749 100644 --- a/_thesis/3_Long_Range_Falicov_Kimball/3.3_LRFK_Results.html +++ b/_thesis/3_Long_Range_Falicov_Kimball/3.3_LRFK_Results.html @@ -78,7 +78,7 @@ image:

Looking at the results of our MCMC simulations, we find a rich phase diagram with a CDW FTPT and interaction-tuned Anderson versus Mott localised phases similar to the 2D FK model  [1]. We explore the localisation properties of the fermionic sector and find that the localisation lengths vary dramatically across the phases and for different energies. The results at moderate system sizes indicate the coexistence of localised and delocalised states within the CDW phase. We then introduce a model of uncorrelated binary disorder on a CDW background. This disorder model gives quantitatively similar behaviour to the LRFK model but we are able to simulate it on much larger systems. For these larger system sizes, we find that all states are eventually localised with a localisation length which diverges towards zero temperature indicating that the results at moderate system size suggestive of coexistence were due to weak localisation effects.

Figure 1: Phase diagrams of the 1D long-range Falicov-Kimball model. (Left) The TJ plane at U = 5: the CDW ordered phase is separated from a disordered Mott insulating phase by a critical temperature T_c, linear in J. (Right) On the TU plane at J = 5, the disordered phase is split into two. At large U there is a Mott insulator phase characterised by the presence of a gap at E=0 in the single particle energy spectrum. At small U there is an Anderson phase characterised by the absence of a gap. U_c is independent of temperature indicating that the FTPT is primarily driven by the long-range coupling term in J. At U = 0 the fermions are decoupled from the spins forming a simple Fermi gas. - +

Phase Diagram

@@ -87,15 +87,15 @@ image:

The CDW transition temperature is largely independent from the strength of the interaction \(U\). This demonstrates that the phase transition is driven by the long-range term \(J\) with little effect from the coupling to the fermions \(U\). The physics of the spin sector in the LRFK model mimics that of the LRI model and is not significantly altered by the presence of the fermions. In 2D the transition to the CDW phase is mediated by an RKYY-like interaction  [3]. However, this is insufficient to stabilise long-range order in 1D. That the critical temperature \(T_c\) does not depend on \(U\) in our model further confirms this.

The main order parameter for this model is the staggered magnetisation \(m = N^{-1} \sum_i (-1)^i S_i\) that signals the onset of a CDW phase at low temperature. However, my main interest concerns the additional structure of the fermionic sector in the high temperature phase. Following Ref.  [1], we can distinguish between the Mott and Anderson insulating phases. The Mott insulator is characterised by a gapped DOS in the absence of a CDW, instead the gap is driven entirely by interaction effect. Thus, the opening of a gap for large \(U\) is distinct from the gap-opening induced by the translational symmetry breaking in the CDW state below \(T_c\). The Anderson phase is gapless but, as we explain below, shows localised fermionic eigenstates. It therefore has an insulating character.

-Figure 2: Energy resolved DOS(\omega) against system size N in all three phases. The Charge Density Wave (CDW) phase is shown in both the high and low U regime for completeness. The top left panel shows the Anderson phase at U = 2 and high T = 2.5, this phase is gapless, but does not conduct due to Anderson localisation. In the lower left pane at U = 2 and low T = 1.5, CDW order sets in, allowing the single particle eigenstates to become extended but opening a gap in their band structure. In the upper right panel at U = 5 and high T = 2.5, the states are localised by disorder and an interaction driven gap opens. This is a Mott insulator. Finally the CDW phase at U = 5 and T = 1.5 is qualitatively similar to the lower left panel except that the gap scales with U. For all the plots J = 5,\;\alpha = 1.25. - +Figure 2: Energy resolved DOS(\omega) against system size N in all three phases. The Charge Density Wave (CDW) phase is shown in both the high and low U regime for completeness. The top left panel shows the Anderson phase at U = 2 and high T = 2.5, this phase is gapless, but does not conduct due to Anderson localisation. In the lower left pane at U = 2 and low T = 1.5, CDW order sets in, allowing the single particle eigenstates to become extended but opening a gap in their band structure. In the upper right panel at U = 5 and high T = 2.5, the states are localised by disorder and an interaction driven gap opens. This is a Mott insulator. Finally, the CDW phase at U = 5 and T = 1.5 is qualitatively similar to the lower left panel except that the gap scales with U. For all the plots J = 5,\;\alpha = 1.25. +

Localisation Properties

-Figure 3: The IPR(\omega) scaling with N at fixed energy for each phase and for points both in the gap (\omega_0) and in a band (\omega_1). The slope of the line yields the scaling exponent \tau defined by \mathrm{IPR} \propto N^{-\tau}). \tau close to zero implies that the states at that energy are localised while \tau = -d corresponds to extended states where d is the system dimension. All but the bands of the charge density wave phase are approximately localised with \tau is very close to zero. The bands in the charge density wave phase are localised with long localisation lengths at finite temperatures that extend to infinity as the temperature approaches zero. For all the plots J = 5,\;\alpha = 1.25. The measured \tau_0,\tau_1 for each figure are: (a) (0.06\pm0.01, 0.02\pm0.01 (b) 0.04\pm0.02, 0.00\pm0.01 (c) 0.05\pm0.03, 0.30\pm0.03 (d) 0.06\pm0.04, 0.15\pm0.05 We show later that the apparent slight scaling of the IPR with system size in the localised cases can be explained by finite size effects due to the changing defect density with system size rather than due to delocalisation of the states. - +Figure 3: The IPR(\omega) scaling with N at fixed energy for each phase and for points both in the gap (\omega_0) and in a band (\omega_1). The slope of the line yields the scaling exponent \tau defined by \mathrm{IPR} \propto N^{-\tau}. \tau close to zero implies that the states at that energy are localised while \tau = -d corresponds to extended states where d is the system dimension. All but the bands of the charge density wave phase are approximately localised with \tau is very close to zero. The bands in the charge density wave phase are localised with long localisation lengths at finite temperatures that extend to infinity as the temperature approaches zero. For all the plots J = 5,\;\alpha = 1.25. The measured \tau_0,\tau_1 for each figure are: (a) 0.06\pm0.01, 0.02\pm0.01 (b) 0.04\pm0.02, 0.00\pm0.01 (c) 0.05\pm0.03, 0.30\pm0.03 (d) 0.06\pm0.04, 0.15\pm0.05 We show later that the apparent slight scaling of the IPR with system size in the localised cases can be explained by finite size effects due to the changing defect density with system size rather than due to delocalisation of the states. +

The MCMC formulation suggests viewing the spin configurations as a form of annealed binary disorder whose probability distribution is given by the Boltzmann weight \(e^{-\beta H_S[\vec{S}] - \beta F_c[\vec{S}]}\). This makes apparent the link to the study of disordered systems and Anderson localisation. These systems are typically studied by defining the probability distribution for the quenched disorder potential externally. Here, by contrast, we have a translation invariant system with disorder as a natural consequence of the Ising background field conserved by the dynamics.

@@ -104,14 +104,14 @@ image:

In the limits of zero and infinite temperature, our model becomes a simple tight-binding model for the fermions. At zero temperature, the spin background is in one of the two translation invariant AFM ground states with two gapped fermionic CDW bands at energies

\[E_{\pm} = \pm\sqrt{\frac{1}{4}U^2 + 2t^2(1 + \cos ka)^2}\;.\]

-

At infinite temperature, all the spin configurations become equally likely and the fermionic model reduces to one of binary uncorrelated disorder in which all eigenstates are Anderson localised  [4]. An Anderson localised state, centred around \(r_0\), has magnitude that drops exponentially over some localisation length \(\xi\) i.e \(|\psi(r)|^2 \sim \exp{-|r - r_0|/\xi}\). Calculating \(\xi\) directly is numerically demanding. Therefore, we determine if a given state is localised via the energy-resolved IPR and the DOS defined as

+

At infinite temperature, all the spin configurations become equally likely and the fermionic model reduces to one of binary uncorrelated disorder in which all eigenstates are Anderson localised  [4]. An Anderson localised state, centred around \(r_0\), has magnitude that drops exponentially over some localisation length \(\xi\) i.e., \(|\psi(r)|^2 \sim \exp{-|r - r_0|/\xi}\). Calculating \(\xi\) directly is numerically demanding. Therefore, we determine if a given state is localised via the energy-resolved IPR and the DOS defined as

\[\begin{aligned} \mathrm{DOS}(\vec{S}, \omega)& = N^{-1} \sum_{i} \delta(\epsilon_i - \omega)\\ -\mathrm{IPR}(\vec{S}, \omega)& = \; N^{-1} \mathrm{DOS}(\vec{S}, \omega)^{-1} \sum_{i,j} \delta(\epsilon_i - \omega)\;\psi^{4}_{i,j}\end{aligned}\]

+\mathrm{IPR}(\vec{S}, \omega)& = \; N^{-1} \mathrm{DOS}(\vec{S}, \omega)^{-1} \sum_{i,j} \delta(\epsilon_i - \omega)\;\psi^{4}_{i,j},\end{aligned}\]

where \(\epsilon_i\) and \(\psi_{i,j}\) are the \(i\)th energy level and \(j\)th element of the corresponding eigenfunction, both dependent on the background spin configuration \(\vec{S}\).

-Figure 5: The DOS (a) and scaling exponent \tau (b) as a function of energy for the CDW phase to the gapless Anderson insulating phase at U=2. Regions where the DOS is close to zero are shown in white. The scaling exponent \tau is obtained from fits to IPR(N) = A N^{-\lambda} for a range of system sizes. J = 5,\;\alpha = 1.25 - +Figure 5: The DOS (a) and scaling exponent \tau (b) as a function of energy for the CDW phase to the gapless Anderson insulating phase at U=2. Regions where the DOS is close to zero are shown in white. The scaling exponent \tau is obtained from fits to \mathrm{IPR}(N) = A N^{-\lambda} for a range of system sizes. J = 5,\;\alpha = 1.25 +

The scaling of the IPR with system size \(\mathrm{IPR} \propto N^{-\tau}\) depends on the localisation properties of states at that energy. For delocalised states, e.g. Bloch waves, \(\tau\) is the physical dimension. For fully localised states \(\tau\) goes to zero in the thermodynamic limit. However, for special types of disorder, such as binary disorder, the localisation lengths can be large, comparable to the system size. This can make it difficult to extract the correct scaling. An additional complication arises from the fact that the scaling exponent may display intermediate behaviours for correlated disorder and in the vicinity of a localisation-delocalisation transition  [5,6]. The thermal defects of the CDW phase lead to a binary disorder potential with a finite correlation length. The key question for our system is then: How is the \(T=0\) CDW phase with fully delocalised fermionic states connected to the fully localised phase at high temperatures?

@@ -124,11 +124,11 @@ image:

In the CDW phases, at \(U=2\) and \(U=5\), we find that states within the gapped CDW bands, e.g. at \(\omega_1\), have scaling exponents \(\tau = 0.30\pm0.03\) and \(\tau = 0.15\pm0.05\), respectively. This surprising finding suggests that the CDW bands are partially delocalised with multi-fractal behaviour of the wavefunctions  [6]. This phenomenon would be unexpected in a 1D model as they generally do not support delocalisation in the presence of disorder except as the result of correlations in the emergent disorder potential  [7,8]. However, we show later, by comparison to an uncorrelated Anderson model, that these nonzero exponents are a finite size effect and the states are localised with a finite \(\xi\) similar to the system size. This is an example of weak localisation. As a result, the IPR does not scale correctly until the system size has grown much larger than \(\xi\). Fig. 7 shows that the scaling of the IPR in the CDW phase does flatten out eventually.

-

Next, we use the DOS and the scaling exponent \(\tau\) to explore the localisation properties over the energy-temperature plane in fig. 5 and fig. 4. Gapped areas are shown in white, which highlights the distinction between the gapped Mott phase and the ungapped Anderson phase. In-gap states appear just below the critical point, smoothly filling the bandgap in the Anderson phase and forming islands in the Mott phase. As in the finite  [9] and infinite dimensional  [10] cases, the in-gap states merge and are pushed to lower energy for decreasing U as the \(T=0\) CDW gap closes. Intuitively, the presence of in-gap states can be understood as a result of domain wall fluctuations away from the AFM ordered background. These domain walls act as local potentials for impurity-like bound states  [9].

-

To understand the localisation properties we can compare the behaviour of our model with that of a simpler Anderson disorder model in which the spins are replaced by a CDW background with uncorrelated binary defect potentials. This is defined by replacing the spin degree of freedom in the FK model \(S_i = \pm \tfrac{1}{2}\) with a disorder potential \(d_i = \pm \tfrac{1}{2}\) controlled by a defect density \(\rho\) such that \(d_i = -\tfrac{1}{2}\) with probability \(\rho/2\) and \(d_i = \tfrac{1}{2}\) otherwise. \(\rho/2\) is used rather than \(\rho\) so that the disorder potential takes on the zero temperature CDW ground state at \(\rho = 0\) and becomes a random choice over spin states at \(\rho = 1\) i.e the infinite temperature limit.

+

Next, we use the DOS and the scaling exponent \(\tau\) to explore the localisation properties over the energy-temperature plane in fig. 5 and fig. 4. Gapped areas are shown in white, which highlights the distinction between the gapped Mott phase and the ungapped Anderson phase. In-gap states appear just below the critical point, smoothly filling the bandgap in the Anderson phase and forming islands in the Mott phase. As in the finite  [9] and infinite dimensional  [10] cases, the in-gap states merge and are pushed to lower energy for decreasing \(U\) as the \(T=0\) CDW gap closes. Intuitively, the presence of in-gap states can be understood as a result of domain wall fluctuations away from the AFM ordered background. These domain walls act as local potentials for impurity-like bound states  [9].

+

To understand the localisation properties we can compare the behaviour of our model with that of a simpler Anderson disorder model in which the spins are replaced by a CDW background with uncorrelated binary defect potentials. This is defined by replacing the spin degree of freedom in the FK model \(S_i = \pm \tfrac{1}{2}\) with a disorder potential \(d_i = \pm \tfrac{1}{2}\) controlled by a defect density \(\rho\) such that \(d_i = -\tfrac{1}{2}\) with probability \(\rho/2\) and \(d_i = \tfrac{1}{2}\) otherwise. \(\rho/2\) is used rather than \(\rho\) so that the disorder potential takes on the zero temperature CDW ground state at \(\rho = 0\) and becomes a random choice over spin states at \(\rho = 1\) i.e., the infinite temperature limit.

\[\begin{aligned} -H_{\mathrm{DM}} = & \;U \sum_{i} (-1)^i \; d_i \;(c^\dagger_{i}c_{i} - \tfrac{1}{2}) \\ -& -\;t \sum_{i} c^\dagger_{i}c_{i+1} + c^\dagger_{i+1}c_{i} +H_{\mathrm{DM}} = & \;U \sum_{i} (-1)^i \; d_i \;(c^\dagger_{i}c^{\phantom{\dagger}}_{i} - \tfrac{1}{2}) \\ +& -\;t \sum_{i} c^\dagger_{i}c^{\phantom{\dagger}}_{i+1} + c^\dagger_{i+1}c^{\phantom{\dagger}}_{i}, \end{aligned}\]

Figures -fig. 6 and -fig. 7 compare the FK model to the disorder model at different system sizes, matching the defect densities of the disorder model to the FK model at \(N = 270\) above and below the CDW transition. We find very good, even quantitative, agreement between the FK and disorder models. This suggests that correlations in the spin sector do not play a significant role.

As we can sample directly from the disorder model, rather than through MCMC, the samples are uncorrelated. Hence we can evaluate much larger system sizes with the disorder model. This enables us to pin down the correct localisation effects. In particular, what appear to be delocalised states for small system sizes eventually turn out to be states with large localisation length. The localisation length diverges towards the ordered zero temperature CDW state. The interplay of interactions, which here produces a binary potential, and localisation can be very intricate and the added advantage of a 1D model is that we can explore very large system sizes.

diff --git a/_thesis/4_Amorphous_Kitaev_Model/4.1_AMK_Model.html b/_thesis/4_Amorphous_Kitaev_Model/4.1_AMK_Model.html index 52feda8..bc82df2 100644 --- a/_thesis/4_Amorphous_Kitaev_Model/4.1_AMK_Model.html +++ b/_thesis/4_Amorphous_Kitaev_Model/4.1_AMK_Model.html @@ -87,13 +87,13 @@ image:

The material in this chapter expands on work presented in

  [1] Cassella, G., D’Ornellas, P., Hodson, T., Natori, W. M., & Knolle, J. (2022). An exact chiral amorphous spin liquid. arXiv preprint arXiv:2208.08246.

All the code is available online as a Python package called Koala  [2].

-

This was a joint project of Gino, Peru and me with advice and guidance from Willian and Johannes, all authors of the above. The project grew out of an interest the three of us had in studying amorphous systems, coupled with Johannes’ expertise on the Kitaev model. The idea to use Voronoi partitions came from  [3] and Gino did the implementation of this. The idea and implementation of the edge colouring using SAT solvers and the mapping from flux sector to bond sector using A* search were both entirely my work. Peru produced the numerical evidence for the ground state and implemented the local markers. Gino and I did much of the rest of the programming for Koala collaboratively, often pair programming, this included the phase diagram, edge mode and finite temperature analyses as well as the derivation of the projector in the amorphous case.

+

This was a joint project of Gino, Peru and me with advice and guidance from Willian and Johannes, all authors of the above. The project grew out of an interest the three of us had in studying amorphous systems, coupled with Johannes’ expertise on the Kitaev model. The idea to use Voronoi partitions came from ref.  [3] and Gino did the implementation of this. The idea and implementation of the edge colouring using SAT solvers and the mapping from flux sector to bond sector using A* search were both entirely my work. Peru produced the numerical evidence for the ground state and implemented the local markers. Gino and I did much of the rest of the programming for Koala collaboratively, often pair programming, this included the phase diagram, edge mode and finite temperature analyses as well as the derivation of the projector in the amorphous case.

Chapter Summary

-

In this chapter, I will first define the amorphous Kitaev (AK) model and discuss the construction of amorphous lattices. Second, in the methods section I will discuss the details of voronisation and graph colouring. Finally I will present and interpret the results obtained.

+

In this chapter, I will first define the amorphous Kitaev (AK) model and discuss the construction of amorphous lattices. Second, in the methods section I will discuss the details of voronisation and graph colouring. Finally, I will present and interpret the results obtained.

From its introduction it was known that the Kitaev Honeycomb (KH) model is solvable on any trivalent lattice. Consequently, it has been generalised to many such lattices  [4–7] but so far none that entirely lack translation symmetry. Here we will do just that.

-

Amorphous lattices are characterised by local constraints but no long-range order. They arise, for instance, in amorphous semiconductors like Silicon and Germanium  [8,9]. Recent work has shown that topological insulating (TI) phases, characterised by protected edge states and topological bulk invariants, can exist in amorphous systems  [10–16]. TI phases, however, arise in non-interacting systems. In this context, we might ask whether Quantum Spin Liquid (QSL) systems and the Kitaev Honeycomb (KH) model, in particular, could be realised on amorphous lattices. The phases of the KH model have many similarities with TIs but differ in that the KH model is an interacting system. In general, research on amorphous electronic systems has been focused mainly on non-interacting systems with the exception of amorphous superconductivity  [17–21] or very recent work looking to understand the effect of strong electron repulsion in TIs  [22].

+

Amorphous lattices are characterised by local constraints but no long-range order. They arise, for instance, in amorphous semiconductors like silicon and germanium  [8,9]. Recent work has shown that topological insulating (TI) phases, characterised by protected edge states and topological bulk invariants, can exist in amorphous systems  [10–16]. TI phases, however, arise in non-interacting systems. In this context, we might ask whether Quantum Spin Liquid (QSL) systems and the Kitaev Honeycomb (KH) model, in particular, could be realised on amorphous lattices. The phases of the KH model have many similarities with TIs but differ in that the KH model is an interacting system. In general, research on amorphous electronic systems has been focused mainly on non-interacting systems with the exception of amorphous superconductivity  [17–21] or very recent work looking to understand the effect of strong electron repulsion in TIs  [22].

The KH model is a magnetic system. Magnetism in amorphous systems has been investigated since the 1960s, mostly through the adaptation of theoretical tools developed for disordered systems  [23–26]. This is not always ideal, we have already seen that the topological disorder of amorphous lattices can be qualitatively different from standard bond or site disorder, especially in 2D  [27,28]. Research focused on classical Heisenberg and Ising models has accounted for the observed behaviour of ferromagnetism, disordered antiferromagnetism and widely observed spin glass behaviour  [29]. However, the role of the spin-anisotropic interactions and quantum effects that we see in the KH model has not been addressed in amorphous magnets. It is an open question whether frustrated magnetic interactions on amorphous lattices can give rise to genuine quantum phases such as QSLs  [30–33]. This chapter will answer that question by demonstrating that the Kitaev model on amorphous lattices leads to a kind of QSL called a chiral spin liquid.

In this section, I will discuss how to generalise the KH to an amorphous lattice. The methods section discusses how to generate amorphous lattices using Voronoi partitions of the plane  [10,12], colour them using a SAT solver and how to map back and forth between gauge field configurations and flux configurations. In the results section, I will show extensive numerical evidence that the AK model follows the simple generalisation to Lieb’s theorem  [34] found by other works  [4–7]. I then map out the phase diagram of the AK model and show that the chiral phase around the symmetric point (\(J_x = J_y = J_z\)) is gapped and non-Abelian. We use a quantised local Chern number \(\nu\)  [10,35] as well as the presence of protected chiral Majorana edge modes to determine this. Finally, I look at the role of finite temperature fluctuations and show that the proliferation of flux excitations leads to an Anderson transition, similar to that of the Falicov-Kimball model, to a thermal metal phase  [36–38]. Finally, I consider possible physical realisations of the AK model and other generalisations.

@@ -101,19 +101,19 @@ image:

The Model

-Figure 1: (a) The standard Kitaev model is defined on a honeycomb lattice. The special feature of the honeycomb lattice that makes the model solvable is that each vertex is joined by exactly three bonds, i.e. the lattice is trivalent. One of three labels is assigned to each (b). We represent the antisymmetric gauge degree of freedom u_{jk} = \pm 1 with arrows that point in the direction u_{jk} = +1 (c). The Majorana transformation can be visualised as breaking each spin into four Majoranas which then pair along the bonds. The pairs of x,y and z Majoranas become part of the classical \mathbb{Z}_2 gauge field u_{ij}. This leaves a single Majorana c_i per site. - +Figure 1: (a) The standard Kitaev model is defined on a honeycomb lattice. The special feature of the honeycomb lattice that makes the model solvable is that each vertex is joined by exactly three bonds, i.e., the lattice is trivalent. One of three labels is assigned to each (b). We represent the antisymmetric gauge degree of freedom u_{jk} = \pm 1 with arrows that point in the direction u_{jk} = +1 (c). The Majorana transformation can be visualised as breaking each spin into four Majoranas which then pair along the bonds. Pairs of b_i^x,\;b_i^y and b_i^z Majoranas become part of the classical \mathbb{Z}_2 gauge field u_{ij}. This leaves a single Majorana c_i per site. +
-

The KH model is solvable on any lattice which satisfies two properties: it must be trivalent and it must three-edge-colourable. The first property means every vertex must have three edges attached to it  [39,40]. 2D Voronoi lattices are a well studied class of amorphous trivalent lattices  [10,12,41]. Given a set of seed points, the Voronoi partition divides the plane into basins, based on which seed point is closest by some metric, usually the euclidean metric. The basins of each seed point form the plaquettes of the resulting lattices, while the boundaries become the edges. The Voronoi partition exists in arbitrary dimension \(d\) and produces lattices with degree \(d+1\) except for degenerate cases with measure zero  [42,43]. Voronoi lattices in 2D are trivalent so lend themselves naturally to the Kitaev model.

-

Other methods of lattice generation exist. One can connect randomly placed sites based on proximity  [11] or create simplices from random sites  [44]. However these methods do not present a natural way to restrict the vertex degree to a constant. The most unbiased way to select trivalent graphs would be to sample uniformly from the space of possible trivalent graphs. There has been some work on how to do this using a Markov Chain Monte Carlo approach  [45]. However, it does not guarantee that the resulting graph is planar, which is necessary to be able to three-edge-colour the lattice, our second constraint.

+

The KH model is solvable on any lattice which satisfies two properties: it must be trivalent and it must three-edge-colourable. The first property means every vertex must have three edges attached to it  [39,40]. 2D Voronoi lattices are a well-studied class of amorphous trivalent lattices  [10,12,41]. Given a set of seed points, the Voronoi partition divides the plane into basins, based on which seed point is closest by some metric, usually the euclidean metric. The basins of each seed point form the plaquettes of the resulting lattices, while the boundaries become the edges. The Voronoi partition exists in arbitrary dimension \(d\) and produces lattices with degree \(d+1\) except for degenerate cases with measure zero  [42,43]. Voronoi lattices in 2D are trivalent so lend themselves naturally to the Kitaev model.

+

Other methods of lattice generation exist. One can connect randomly placed sites based on proximity  [11] or create simplices from random sites  [44]. However, these methods do not present a natural way to restrict the vertex degree to a constant. The most unbiased way to select trivalent graphs would be to sample uniformly from the space of possible trivalent graphs. There has been some work on how to do this using a Markov Chain Monte Carlo approach  [45]. However, it does not guarantee that the resulting graph is planar, which is necessary to be able to three-edge-colour the lattice, our second constraint.

The second constraint, three-edge-colourability, requires that we must be able to assign labels to each bond \(\{x,y,z\}\) such that no two edges with the same label meet at a vertex. Such an assignment is known as a three-edge-colouring. For translation invariant models we need only find a solution for the unit cell. This problem is usually small enough that this can be done by hand or using symmetry. For amorphous lattices, the difficulty is that, to the best of my knowledge, the problem of edge-colouring these lattices in general is in NP. To find colourings in practice, we will employ a standard method from the computer science literature for finding solutions of NP problems called a SAT solver, this is discussed in more detail in the methods secton.

-

We find that for large lattices there are many valid colourings. In the isotropic case \(J^\alpha = 1\) the colouring has no physical significance as the definition of the four Majoranas at a site is arbitrary. In the anisotropic case this symmetry is broken at the local level but we nevertheless expect the lattices to exhibit a self averaging behaviour in larger systems such that the choice of colouring doesn’t matter.

+

We find that for large lattices there are many valid colourings. In the isotropic case \(J^\alpha = 1\) the colouring has no physical significance as the definition of the four Majoranas at a site is arbitrary. In the anisotropic case this symmetry is broken at the local level but we nevertheless expect the lattices to exhibit a self-averaging behaviour in larger systems such that the choice of colouring doesn’t matter.

Figure 2: (Bond Sector) A state in the bond sector is specified by assigning \pm 1 to each edge of the lattice. However, this description has a substantial gauge degeneracy. To remove it, we decompose each state into the product of three kinds of objects: (Flux Sector) The main physically relevant quantities. Only a small number of bonds need to be flipped (compared to some arbitrary fixed reference) to reconstruct the flux sector. Here, the edges are chosen from a spanning tree of the dual lattice, so there are no loops. (Gauge Field) The ‘loopiness’ of the bond sector is in this part. This is a network of loops that can always be written as a product of the gauge operators D_j. (Topological Sector) Finally, there are two loops that have no effect on the vortex sector, nor can they be constructed from gauge symmetries D_j. These can be thought of as two fluxes \Phi_{x/y} that thread through the major and minor axes of the torus. Measuring \Phi_{x/y} amounts to constructing Wilson loops around the axes of the torus. We can flip the value of \Phi_{x} by transporting a vortex pair around the torus in the y direction, as shown here. In each of the three figures on the right, black bonds correspond to those that must be flipped, while red line are those same edges on the dual lattice. Composing the three objects together gives back the original bond sector on the left.
-

On a lattice with the above properties, the solution for the KH model laid out in section 2.2 remains applicable to our AK model. See fig. 1 for an example lattice generated by our method. The main differences are twofold. Firstly, the lattices are no longer bipartite in general and therefore contain plaquettes with an odd number of sides which enclose flux \(\pm i\). This leads the AK model to have a ground state with spontaneously broken chiral symmetry  [7,46–52]. This is analogous to the behaviour of the original Kitaev model in response to a magnetic field. One ground state is related to the other by globally inverting the imaginary \(\phi_i\) fluxes  [47]. Secondly, as the model is no longer translationally invariant, Lieb’s theorem for the ground state flux sector no longer applies. However as discussed in the background, a simple generalisation of Lieb’s theorem has been shown numerically to be applicable to many generalised Kitaev models  [4–7]. This generalisation states that the ground state flux configuration depends on the number of sides of each plaquette as

-

\[\phi = -(\pm i)^{n_{\mathrm{sides}}}\qquad{(1)}\]

+

On a lattice with the above properties, the solution for the KH model laid out in section 2.2 remains applicable to our AK model. See fig. 1 for an example lattice generated by our method. The main differences are twofold. Firstly, the lattices are no longer bipartite in general and therefore contain plaquettes with an odd number of sides which enclose flux \(\pm i\). This leads the AK model to have a ground state with spontaneously broken chiral symmetry  [7,46–52]. This is analogous to the behaviour of the original Kitaev model in response to a magnetic field. One ground state is related to the other by globally inverting the imaginary \(\phi_i\) fluxes  [47]. Secondly, as the model is no longer translationally invariant, Lieb’s theorem for the ground state flux sector no longer applies. However, as discussed in the background, a simple generalisation of Lieb’s theorem has been shown numerically to be applicable to many generalised Kitaev models  [4–7]. This generalisation states that the ground state flux configuration depends on the number of sides of each plaquette as

+

\[\phi = -(\pm i)^{n_{\mathrm{sides}}},\qquad{(1)}\]

with a twofold global chiral degeneracy (picking either \(+i\) or \(-i\) in eq. 1).

To verify numerically that Lieb’s theorem generalises to the AK model, the obvious approach would be via exhaustive checking of flux configurations. However, this is problematic because the number of states to check scales exponentially with system size. We side-step this by gluing together two methods, we first work with lattices small enough that we can fully enumerate their flux sectors but tile them to reduce finite size effects. We then show that the effect of tiling scales away with system size.

@@ -136,11 +136,11 @@ The same as @fig:flood_fill but for the amorphous lattice.

The Euler Equation

Euler’s equation provides a convenient way to understand how the states of the AK model factorise into flux sectors, gauge sectors and topological sectors as in fig. 2. The Euler equation states that if we embed a lattice with \(B\) bonds, \(P\) plaquettes and \(V\) vertices onto a closed surface of genus \(g\), (\(0\) for the sphere, \(1\) for the torus) then

-

\[B = P + V + 2 - 2g\]

+

\[B = P + V + 2 - 2g.\]

For the case of the torus where \(g = 1\), we can rearrange this and exponentiate it to read:

-

\[2^B = 2^{P-1}\cdot 2^{V-1} \cdot 2^2\]

+

\[2^B = 2^{P-1}\cdot 2^{V-1} \cdot 2^2.\]

There are \(2^B\) configurations of the bond variables \(\{u_{ij}\}\). Each of these configurations can be uniquely decomposed into a flux sector, a gauge sector and a topological sector, see fig. 2. Each of the \(P\) plaquette operators \(\phi_i\) takes two values but vortices are created in pairs so there are \(2^{P-1}\) vortex sectors in total. There are \(2^{V-1}\) gauge symmetries formed from the \(V\) symmetry operators \(D_i\) because \(\prod_{j} D_j = \mathbb{I}\) is enforced by the projector. Finally, the two topological fluxes \(\Phi_x\) and \(\Phi_y\) account for the last factor of \(2^2\).

-

In a trivalent lattice, there are three bonds for every 2 vertices. Substituting \(3V = 2B\) into Euler’s equation tells us that any trivalent lattice on the torus with \(N\) plaquettes has \(2N\) vertices and \(3N\) bonds. Since each bond is part of two plaquettes this implies that the mean number of sides of a plaquette is exactly six and that odd sides plaquettes must come in pairs.

+

In a trivalent lattice, there are three bonds for every 2 vertices. Substituting \(3V = 2B\) into Euler’s equation tells us that any trivalent lattice on the torus with \(N\) plaquettes has \(2N\) vertices and \(3N\) bonds. Since each bond is part of two plaquettes this implies that the mean number of sides of a plaquette is exactly six and that odd sided plaquettes must come in pairs.

Next Section: Methods

diff --git a/_thesis/4_Amorphous_Kitaev_Model/4.2_AMK_Methods.html b/_thesis/4_Amorphous_Kitaev_Model/4.2_AMK_Methods.html index 8d1d921..4ad277b 100644 --- a/_thesis/4_Amorphous_Kitaev_Model/4.2_AMK_Methods.html +++ b/_thesis/4_Amorphous_Kitaev_Model/4.2_AMK_Methods.html @@ -92,10 +92,10 @@ image:

To be solvable, the AK model requires that each edge in the lattice be assigned a label \(x\), \(y\) or \(z\), such that each vertex has exactly one edge of each type connected to it, a three-edge-colouring. This problem must be distinguished from that considered by the famous four-colour theorem  [9]. The four-colour theorem is concerned with assigning colours to the vertices of planar graphs, such that no vertices that share an edge have the same colour.

-

For a graph of maximum degree \(\Delta\), \(\Delta + 1\) colours are always enough to edge-colour it. An \(\mathcal{O}(mn)\) algorithm exists to do this for a graph with \(m\) edges and \(n\) vertices  [10]. Graphs with \(\Delta = 3\) are known as cubic graphs. Cubic graphs can be four-edge-coloured in linear time  [11]. However we need a three-edge-colouring of our cubic graphs, which turns out to be more difficult. Cubic, planar, bridgeless graphs can be three-edge-coloured if and only if they can be four-face-coloured  [12]. Bridges are edges that connect otherwise disconnected components. An \(\mathcal{O}(n^2)\) algorithm exists for these  [13]. However, it is not clear whether this extends to cubic, toroidal bridgeless graphs.

-

A four-face-colouring is equivalent to a four-vertex-colouring of the dual graph, see appendix A.4. So if we could find a four-vertex-colouring of the dual graph we would be done. However vertex-colouring a toroidal graph may require up to seven colours  [14]! The complete graph of seven vertices \(K_7\) is a good example of a toroidal graph that requires seven colours.

+

For a graph of maximum degree \(\Delta\), \(\Delta + 1\) colours are always enough to edge-colour it. An \(\mathcal{O}(mn)\) algorithm exists to do this for a graph with \(m\) edges and \(n\) vertices  [10]. Graphs with \(\Delta = 3\) are known as cubic graphs. Cubic graphs can be four-edge-coloured in linear time  [11]. However, we need a three-edge-colouring of our cubic graphs, which turns out to be more difficult. Cubic, planar, bridgeless graphs can be three-edge-coloured if and only if they can be four-face-coloured  [12]. Bridges are edges that connect otherwise disconnected components. An \(\mathcal{O}(n^2)\) algorithm exists for these  [13]. However, it is not clear whether this extends to cubic, toroidal bridgeless graphs.

+

A four-face-colouring is equivalent to a four-vertex-colouring of the dual graph, see appendix A.4. So if we could find a four-vertex-colouring of the dual graph we would be done. However, vertex-colouring a toroidal graph may require up to seven colours  [14]! The complete graph of seven vertices \(K_7\) is a good example of a toroidal graph that requires seven colours.

Luckily, some problems are easier in practice. Three-edge-colouring cubic toroidal graphs is one of those things. To find colourings, we use a Boolean Satisfiability Solver or SAT solver. A SAT problem is a set of statements about a set of boolean variables \([x_1, x_2\ldots]\), such as “\(x_1\) or not \(x_3\) is true”. A solution to a SAT problem is an assignment \(x_i \in {0,1}\) that satisfies all the statements  [15]. General purpose, high performance programs for solving SAT problems have been an area of active research for decades  [16]. Such programs are useful because, by the Cook-Levin theorem  [17,18], any NP problem can be encoded (in polynomial time) as an instance of a SAT problem. This property is what makes SAT one of the subset of NP problems called NP-Complete. It is a relatively standard technique in the computer science community to solve NP problems by first transforming them to SAT instances and then using an off-the-shelf SAT solver. The output of this can then be mapped back to the original problem domain.

-

Whether graph colouring problems are in NP or P seems to depend delicately on the class of graphs considered, the maximum degree and the number of colours used. It is therefore possible that a polynomial time algorithm may exist for our problem. However using a SAT solver turns out to be fast enough in practice that it is by no means the rate limiting step for generating and solving instances of the AK model. In appendix A.4 I detail the specifics of how I mapped edge-colouring problems to SAT instances and show a breakdown of where the computational effort is spent, the majority being on matrix diagonalisation.

+

Whether graph colouring problems are in NP or P seems to depend delicately on the class of graphs considered, the maximum degree and the number of colours used. It is therefore possible that a polynomial time algorithm may exist for our problem. However, using a SAT solver turns out to be fast enough in practice that it is by no means the rate limiting step for generating and solving instances of the AK model. In appendix A.4 I detail the specifics of how I mapped edge-colouring problems to SAT instances and show a breakdown of where the computational effort is spent, the majority being on matrix diagonalisation.

Mapping between flux sectors and bond sectors

diff --git a/_thesis/4_Amorphous_Kitaev_Model/4.3_AMK_Results.html b/_thesis/4_Amorphous_Kitaev_Model/4.3_AMK_Results.html index 0f53ced..c0be5cd 100644 --- a/_thesis/4_Amorphous_Kitaev_Model/4.3_AMK_Results.html +++ b/_thesis/4_Amorphous_Kitaev_Model/4.3_AMK_Results.html @@ -94,28 +94,28 @@ image:

Results

-

This section contains our results on the AK model, we first look at how we checked numerically that Lieb’s theorem generalises to our model. Next we compute the ground state diagram and look at the two phases that arise there. We then use a local Chern marker and the presence of edge modes to characterise these phases as having Abelian or non-Abelian statistics. Finally we look at the finite temperature behaviour of the model.

+

This section contains our results on the AK model, we first look at how we checked numerically that Lieb’s theorem generalises to our model. Next we compute the ground state diagram and look at the two phases that arise there. We then use a local Chern marker and the presence of edge modes to characterise these phases as having Abelian or non-Abelian statistics. Finally, we look at the finite temperature behaviour of the model.

The Ground State Flux Sector

We will check that Lieb’s theorem generalises to our model by enumerating all the flux sectors of many separate amorphous lattice realisations. However, we have two seemingly irreconcilable problems. Finite size effects have a large energetic contribution for small systems  [1] so we would like to perform our analysis for very large lattices. For an amorphous system with \(N\) plaquettes, \(2N\) edges and \(3N\) vertices we have \(2^{N-1}\) flux sectors to check and diagonalisation scales with \(\mathcal{O}(N^3)\). That exponential scaling makes it difficult to work with lattices much larger than \(16\) plaquettes with the resources.

To get around this, we instead look at periodic systems with amorphous unit cells. For a similarly sized periodic system with \(A\) unit cells and \(B\) plaquettes in each unit cell where \(N \sim AB\) things get much better. We can use Bloch’s theorem to diagonalise this system in about \(\mathcal{O}(A B^3)\) operations, and more importantly there are only \(2^{B-1}\) flux sectors to check. We fully enumerated the flux sectors of \(\sim\) 25,000 periodic systems with disordered unit cells of up to \(B = 16\) plaquettes and \(A = 100\) unit cells. However, showing that our guess is correct for periodic systems with disordered unit cells is not quite convincing on its own as we have effectively removed longer-range disorder from our lattices.

-

The second part of the argument is to show that the energetic effect of introducing periodicity scales away as we go to larger system sizes and has already diminished to a small enough value at 16 plaquettes, which is indeed what we find. From this, we argue that the results for small periodic systems generalise to large amorphous systems. In the isotropic case (\(J^\alpha = 1\)), Lieb’s theorem correctly predicts the ground state flux sector for all of the lattices we tested. For the toric code phase (\(J^x = J^y = 0.25, J^z = 1\)) all but around (\(\sim 0.5 \%\)) lattices had ground states conforming to our conjecture. In these cases, the energy difference between the true ground state and our prediction was on the order of \(10^{-6} J\).

-

The spin Kitaev Hamiltonian is real and therefore has time reversal symmetry. However in the ground state the flux \(\phi_p\) through any plaquette with an odd number of sides has imaginary eigenvalues \(\pm i\). Thus, states with a fixed flux sector spontaneously break time reversal symmetry. Kiteav noted this in his original paper but it was first explored in a concrete model by Yao and Kivelson for a translation invariant Kitaev model with odd sided plaquettes  [2].

+

The second part of the argument is to show that the energetic effect of introducing periodicity scales away as we go to larger system sizes and has already diminished to a small enough value at 16 plaquettes, which is indeed what we find. From this, we argue that the results for small periodic systems generalise to large amorphous systems. In the isotropic case (\(J^\alpha = 1\)), Lieb’s theorem correctly predicts the ground state flux sector for all of the lattices we tested. For the toric code phase (\(J^x = J^y = 0.25, J^z = 1\)) all but around \(\sim 0.5 \%\) of lattices had ground states conforming to our conjecture. In these cases, the energy difference between the true ground state and our prediction was on the order of \(10^{-6}]\;J\).

+

The spin Kitaev Hamiltonian is real and therefore has time reversal symmetry. However, in the ground state the flux \(\phi_p\) through any plaquette with an odd number of sides has imaginary eigenvalues \(\pm i\). Thus, states with a fixed flux sector spontaneously break time reversal symmetry. Kiteav noted this in his original paper but it was first explored in a concrete model by Yao and Kivelson for a translation invariant Kitaev model with odd sided plaquettes  [2].

Flux sectors come in degenerate pairs, where time reversal is equivalent to inverting the flux through every odd plaquette, a general feature for lattices with odd plaquettes  [3,4]. This spontaneously broken symmetry serves a role analogous to the external magnetic field in the original honeycomb model, leading the AK model to have a non-Abelian anyonic phase without an external magnetic field.

Ground State Phase Diagram

The triangular \(J_x, J_y, J_z\) phase diagram of this family of models arises from setting the energy scale with \(J_x + J_y + J_z = 1\). The intersection of this plane and the unit cube is what yields the equilateral triangles seen in diagrams like fig. 1. The KH model has an Abelian gapped phase in the anisotropic region (the A phase) and is gapless in the isotropic region. The introduction of a magnetic field breaks the chiral symmetry, leading to the isotropic region becoming a gapped, non-Abelian phase, the B phase.

-

Similar to the KH model with a magnetic field, we find that the amorphous model is only gapless along critical lines, see fig. 1 (Left). Interestingly, in finite size systems the gap closing exists in only one of the four topological sectors though the sectors become degenerate in the thermodynamic limit. Nevertheless, this could be a useful way to define the (0, 0) topological flux sector for the amorphous model which otherwise has no natural way to choose it.

-

In the honeycomb model, the phase boundaries are located on the straight lines \(|J^x| = |J^y| \;+ \;|J^z|\) and permutations of \(x,y,z\). These are shown as dotted lines in fig. 1 (Right). We find that on the amorphous lattice these boundaries exhibit an inward curvature, similar to honeycomb Kitaev models with flux or bond disorder  [5–10].

+

Similar to the KH model with a magnetic field, we find that the amorphous model is only gapless along critical lines, see the left panel of fig. 1 (left panel). Interestingly, in finite size systems the gap closing exists in only one of the four topological sectors though the sectors become degenerate in the thermodynamic limit. Nevertheless, this could be a useful way to define the (0, 0) topological flux sector for the amorphous model which otherwise has no natural way to choose it.

+

In the honeycomb model, the phase boundaries are located on the straight lines \(|J^x| = |J^y| \;+ \;|J^z|\) and permutations of \(x,y,z\). These are shown as dotted lines in fig. 1 (right panel). We find that on the amorphous lattice these boundaries exhibit an inward curvature, similar to honeycomb Kitaev models with flux or bond disorder  [5–10].

Figure 1: The phase diagram of the model can be characterised by an equilateral triangle whose corners indicate points where J_\alpha = 1, J_\beta = J_\gamma = 0 while the centre denotes J_x = J_y = J_z. (Center) To compute critical lines efficiently in this space, we evaluate the order parameter of interest along rays shooting from the corners of the phase diagram. The ray highlighted in red defines the values of J used for the left figure. (Left) The fermion gap as a function of J for an amorphous system with 20 plaquettes, where the x axis is the position on the red line in the central figure from 0 to 1. For finite size systems the four topological sectors are not degenerate and only one of them (in green) has a true gap closing. (Right) The Abelian A_\alpha phases of the model and the non-Abelian B phase separated by critical lines where the fermion gap closes. Later we will show that the Chern number \nu changes from 0 to \pm 1 from the A phases to the B phase. Indeed, the gap must close in order for the Chern number to change  [11]. - +

Abelian or non-Abelian statistics of the Gapped Phase

The two phases of the amorphous model are gapped as we can see from the finite size scaling of fig. 4. The next question is: do these phases support excitations with trivial, Abelian or non-Abelian statistics? To answer that we turn to Chern numbers  [12–14]. As discussed earlier the Chern number is a quantity intimately linked to both the topological properties and the anyonic statistics of a model. Here we will make use of the fact that the Abelian/non-Abelian character of a model is linked to its Chern number  [1]. The Chern number is only defined for the translation invariant case so we instead use a family of real space generalisations of the Chern number that work for amorphous systems called local topological markers  [15–17].

-

There are many possible choices here, indeed Kitaev defines one in his original paper on the KH model  [1]. Here we use the crosshair marker of  [18] because it works well on smaller systems. We calculate the projector \(P = \sum_i |\psi_i\rangle \langle \psi_i|\) onto the occupied fermion eigenstates of the system in open boundary conditions. The projector encodes local information about the occupied eigenstates of the system and in gapped systems it is exponentially localised  [19]. The name crosshair comes from the fact that the marker is defined with respect to a particular point \((x_0, y_0)\) by step functions in x and y

+

There are many possible choices here, indeed Kitaev defines one in his original paper on the KH model  [1]. Here we use the crosshair marker of  [18] because it works well on smaller systems. We calculate the projector \(P = \sum_i |\psi_i\rangle \langle \psi_i|\) onto the occupied fermion eigenstates of the system in open boundary conditions. The projector encodes local information about the occupied eigenstates of the system and in gapped systems it is exponentially localised  [19]. The name crosshair comes from the fact that the marker is defined with respect to a particular point \((x_0, y_0)\) by step functions in \(x\) and \(y\)

\[\begin{aligned} \nu (x, y) = 4\pi \; \Im\; \mathrm{Tr}_{\mathrm{B}} \left ( @@ -161,7 +161,7 @@ image:

Discussion and Conclusion

-

In this chapter we have looked at an extension of the KH model to amorphous lattices with coordination number three. We discussed a method to construct arbitrary trivalent lattices using Voronoi partitions, how to embed them onto the torus and how to edge-colour them using a SAT solver. We showed numerically that the ground state flux sector of the model is given by a simple extension of Lieb’s theorem. The model has two gapped QSL phases. The two phases support excitations with different anyonic statistics, Abelian and non-Abelian, distinguished using a real-space generalisation of the Chern number  [18]. The presence of odd-sided plaquettes in the model resulted in spontaneous breaking of time reversal symmetry, leading to the emergence of a chiral spin liquid phase. Finally we showed evidence that the amorphous system undergoes an Anderson transition to a thermal metal phase, driven by the proliferation of vortices with increasing temperature. The AK model is an exactly solvable model of the chiral QSL state, one of the first models to exhibit a topologically non-trivial quantum many-body phase on an amorphous lattice. As such this study provides a number of future lines of research.

+

In this chapter we have looked at an extension of the KH model to amorphous lattices with coordination number three. We discussed a method to construct arbitrary trivalent lattices using Voronoi partitions, how to embed them onto the torus and how to edge-colour them using a SAT solver. We showed numerically that the ground state flux sector of the model is given by a simple extension of Lieb’s theorem. The model has two gapped QSL phases. The two phases support excitations with different anyonic statistics, Abelian and non-Abelian, distinguished using a real-space generalisation of the Chern number  [18]. The presence of odd-sided plaquettes in the model resulted in spontaneous breaking of time reversal symmetry, leading to the emergence of a chiral spin liquid phase. Finally, we showed evidence that the amorphous system undergoes an Anderson transition to a thermal metal phase, driven by the proliferation of vortices with increasing temperature. The AK model is an exactly solvable model of the chiral QSL state, one of the first models to exhibit a topologically non-trivial quantum many-body phase on an amorphous lattice. As such this study provides a number of future lines of research.

Experimental Realisations and Signatures

We should consider whether a physical amorphous system that supports a QSL ground state could exist. The search for translation invariant Kitaev systems is already motivated by the prospect of a physically realised QSL state, Majorana fermions and direct access to a system with emergent \(\mathbb{Z}_2\) gauge physics  [29]. In addition to all this, an amorphous Kitaev model would provide the possibility of exploring the CSL state and potentially very different routes to a physical realisation. One route would be to ask if any crystalline Kitaev material candidates can be heated and rapidly quenched  [30–32] to produce amorphous analogues that might preserve enough of their local structure to support a QSL state.

@@ -175,7 +175,7 @@ image:

The KH model can be extended to 3D either on trivalent lattices  [54,55] or it can be generalised to an exactly solvable spin-\(\tfrac{3}{2}\) model on 3D four-coordinate lattices  [2,57–68]. In  [57], the 2D square lattice with 4 bond types (\(J_w, J_x, J_y, J_z\)) is considered. Since Voronoi partitions in 3D produce lattices of degree four, one interesting generalisation of this work would be to look at the spin-\(\tfrac{3}{2}\) Kitaev model on amorphous lattices.

We did not perform a full Markov Chain Monte Carlo (MCMC) simulation of the AK model at finite temperature but the possible extension to a 3D model with an FTPT would motivate this full analysis. This MCMC simulation would be a numerically challenging task but poses no conceptual barriers  [23,25,27]. Doing this would, first, allow one to look for possible violations of the Harris criterion  [69] for the Ising transition of the flux sector. Recall that topological disorder in 2D has radically different properties to that of other kinds of disorder due to the constraints imposed by the Euler equation and maintaining coordination number which allows it to violate otherwise quite general rules like the Harris criterion  [70,71]. Second, incorporating the projector in addition to MCMC would allow for a full investigation of whether the effect of topological degeneracy is apparent at finite temperatures, this is done for the KH model in  [23].

Next, one could investigate whether a QSL phase may exist for other models defined on amorphous lattices with a view to more realistic prospects of observation. Do the properties of the Kitaev-Heisenberg model generalise from the honeycomb to the amorphous case?  [41,43,72–74] Alternatively we might look at other lattice construction techniques. For instance we could construct lattices by linking close points  [75] or create simplices from random sites  [76]. Lattices constructed using these methods would likely have a large number of lattice defects where \(z \neq 3\) in the bulk, leading to many localised Majorana zero modes.

-

We found a small number of lattices for which Lieb’s theorem did not correctly predict the true ground state flux sector. I see two possibilities for what could cause this. Firstly it could be a finite size effect that is amplified by certain rare lattice configurations. It would be interesting to try to elucidate what lattice features are present when Lieb’s theorem fails. Alternatively, it might be telling that the ground state conjecture failed in the toric code A phase where the couplings are anisotropic. We showed that the colouring does not matter in the B phase. However an avenue that I did not explore was whether the particular choice of colouring for a lattice affects the physical properties in the toric code A phase. It is possible that some property of the particular colouring chosen is what leads to these rare failures of Lieb’s theorem.

+

We found a small number of lattices for which Lieb’s theorem did not correctly predict the true ground state flux sector. I see two possibilities for what could cause this. Firstly, it could be a finite size effect that is amplified by certain rare lattice configurations. It would be interesting to try to elucidate what lattice features are present when Lieb’s theorem fails. Alternatively, it might be telling that the ground state conjecture failed in the toric code A phase where the couplings are anisotropic. We showed that the colouring does not matter in the B phase. However, an avenue that I did not explore was whether the particular choice of colouring for a lattice affects the physical properties in the toric code A phase. It is possible that some property of the particular colouring chosen is what leads to these rare failures of Lieb’s theorem.

Overall, there has been surprisingly little research on amorphous quantum many-body phases despite there being plenty of material candidates. I expect the exact chiral amorphous spin liquid to find many generalisations to realistic amorphous quantum magnets.

Next Chapter: 5 Conclusion

diff --git a/_thesis/5_Conclusion/5_Conclusion.html b/_thesis/5_Conclusion/5_Conclusion.html index 3e2e2c7..ccc456c 100644 --- a/_thesis/5_Conclusion/5_Conclusion.html +++ b/_thesis/5_Conclusion/5_Conclusion.html @@ -70,7 +70,7 @@ image: Figure 1: An example of a Metal Organic Framework (MOF). MOFs are large synthetic organic ligands coordinated with metal ions that can cross link to form crystals. It has been demonstrated that these materials can sustain Kitaev-like interactions  [5]. They can also be made into amorphous materials  [6]. This image from  [7] is in the public domain. -

Unlike the 2D FK model and 1D LRFK models, the KH and AK models don’t have a Finite-Temperature Phase Transition (FTPT). They immediately disorder at any finite temperature  [3]. However, generalisations of the KH model to 3D do in general have an FTPT. Indeed, the role of dimensionality has been a key theme in this work. Both localisation and thermodynamic phenomena depend crucially on dimensionality, with thermodynamic order generally suppressed and localisation effects strengthened in low dimensions. The graph theory that underpins the KH and AK models itself also changes strongly with dimension. Voronisation in 2D produces trivalent lattices, on which the spin-\(1/2\) AK model is exactly solvable. Meanwhile in 3D, voronisation gives us \(z=4\) lattices upon which a spin-\(3/2\) generalisation to the KH model is exactly solvable  [8–10]. Similarly, planar and toroidal graphs are a uniquely 2D construct. Satisfying planarity imposes constraints on the connectivity of planar graphs. This leads amorphous planar graphs to have strong anti-correlations which can violate otherwise robust bounds like the Harris criterion  [11]. Contrast this with Anderson localisation in 1D where only longer range correlations in the disorder can produce surprising effects  [12–17].

+

Unlike the 2D FK model and 1D LRFK models, the KH and AK models don’t have a finite-temperature phase transition (FTPT). They immediately disorder at any finite temperature  [3]. However, generalisations of the KH model to 3D do in general have an FTPT. Indeed, the role of dimensionality has been a key theme in this work. Both localisation and thermodynamic phenomena depend crucially on dimensionality, with thermodynamic order generally suppressed and localisation effects strengthened in low dimensions. The graph theory that underpins the KH and AK models itself also changes strongly with dimension. Voronisation in 2D produces trivalent lattices, on which the spin-\(1/2\) AK model is exactly solvable. Meanwhile in 3D, voronisation gives us \(z=4\) lattices upon which a spin-\(3/2\) generalisation to the KH model is exactly solvable  [8–10]. Similarly, planar and toroidal graphs are a uniquely 2D construct. Satisfying planarity imposes constraints on the connectivity of planar graphs. This leads amorphous planar graphs to have strong anti-correlations which can violate otherwise robust bounds like the Harris criterion  [11]. Contrast this with Anderson localisation in 1D where only longer range correlations in the disorder can produce surprising effects  [12–17].

Looking towards future work, the LRFK model provides multiple possible routes. One interesting idea is to park the model at a thermal critical point. This would generate a scale-free disorder potential, which could potentially lead to complex localisation physics not often seen in 1D  [12,18]. Like other solvable models of disorder-free localisation, the LRFK model should also exhibit intriguing out-of-equilibrium physics, such as slow entanglement dynamics. This could be used to help understand these phenomena in more generic interacting systems  [19]. There is also the rich ground state phenomenology of the FK model as a function of filling  [20], such as the devil’s staircase  [21] as well as superconductor like states  [22]. Could the LRFK model stabilise these at finite temperature? Finally, a topological variant of the LRFK model akin to the Su-Schrieffer-Heeger (SSH) model could be an interesting way to probe the interplay of topological bound states and thermal domain wall defects.

Looking at the AK model, we discussed whether its physics might be realisable in amorphous versions of known KH candidate materials  [23]. Alternatively, we might be able to engineer them in synthetic materials, such as Metal Organic Frameworks (MOFs). Work on MOFs has already explored the possibility of both Kitaev-like interactions and amorphous lattices  [5,6]. It is an open question whether the superexchange couplings that generate Kitaev interactions could survive the transition to an amorphous lattice. If the interactions do survive, there will likely be many defects of different kinds present in the resulting material. These might take the form of dangling bonds, vertex degree \(> 3\) or violations of the colouring conditions. I therefore hope that future work examines how robust the QSL and CSL states of the AK model are to these kinds of disorder. This is a difficult task, as many of these classes of defects would break the integrability of the AK model that we relied on to make the work computationally feasible  [24–28]. While we are considering models with defects, we might consider alternate lattice construction techniques  [29,30]. Lattices constructed using these methods would have defects but may have other desirable properties when compared to Voronoi lattices.

In terms of experimental signatures, we discussed the quantised thermal Hall effect  [31–34], local probes such as spin-polarised scanning tunnelling microscopy  [35–37], and longitudinal heat transport signatures  [38]. One possible difficulty is that the introduction of topological disorder may dilute some of these signatures. On the brighter side, topological disorder may also suppress competing interactions that would otherwise induce magnetic ordering. This could potentially widen the class of materials that could host a QSL or CSL ground state.

diff --git a/_thesis/6_Appendices/A.1.2_Fermion_Free_Energy.html b/_thesis/6_Appendices/A.1.2_Fermion_Free_Energy.html index 2696279..227e10a 100644 --- a/_thesis/6_Appendices/A.1.2_Fermion_Free_Energy.html +++ b/_thesis/6_Appendices/A.1.2_Fermion_Free_Energy.html @@ -70,30 +70,30 @@ image:

Evaluation of the Fermion Free Energy

There are \(2^N\) possible configurations of the spins in the LRFK model. In the language of ions and electrons (immobile and mobile species), we define \(n^k_i\) to be the occupation of the \(i\)th site of the \(k\)th configuration. The quantum part of the free energy can then be defined through the quantum partition function \(\mathcal{Z}^k\) associated with each state \(n^k_i\):

\[\begin{aligned} -F^k &= -1/\beta \ln{\mathcal{Z}^k} \\ +F^k &= -1/\beta \ln{\mathcal{Z}^k}, \\ \end{aligned}\]

-

Such that the overall partition function is:

+

such that the overall partition function is:

\[\begin{aligned} \mathcal{Z} &= \sum_k e^{- \beta H^k} Z^k \\ -&= \sum_k e^{-\beta (H^k + F^k)} \\ +&= \sum_k e^{-\beta (H^k + F^k)}. \\ \end{aligned}\]

-

Because fermions are limited to occupation numbers of 0 or 1, \(Z^k\) simplifies nicely. If \(m^j_i = \{0,1\}\) is defined as the occupation of the level with energy \(\epsilon^k_i\) then the partition function is a sum over all the occupation states labelled by \(j\):

+

Fermions are limited to occupation numbers of 0 or 1, so \(Z^k\) simplifies nicely. If \(m^j_i = \{0,1\}\) is defined as the occupation of the level with energy \(\epsilon^k_i\) then the partition function is a sum over all the occupation states labelled by \(j\):

\[\begin{aligned} Z^k &= \mathrm{Tr} e^{-\beta F^k} = \sum_j e^{-\beta \sum_i m^j_i \epsilon^k_i}\\ &= \sum_j \prod_i e^{- \beta m^j_i \epsilon^k_i}= \prod_i \sum_j e^{- \beta m^j_i \epsilon^k_i}\\ &= \prod_i (1 + e^{- \beta \epsilon^k_i})\\ -F^k &= -1/\beta \sum_k \ln{(1 + e^{- \beta \epsilon^k_i})} +F^k &= -1/\beta \sum_k \ln{(1 + e^{- \beta \epsilon^k_i})}. \end{aligned}\]

Observables can then be calculated from the partition function, for examples the occupation numbers:

\[\begin{aligned} \langle N \rangle &= \frac{1}{\beta} \frac{1}{Z} \frac{\partial Z}{\partial \mu} = - \frac{\partial F}{\partial \mu}\\ &= \frac{1}{\beta} \frac{1}{Z} \frac{\partial}{\partial \mu} \sum_k e^{-\beta (H^k + F^k)}\\ - &= 1/Z \sum_k (N^k_{\mathrm{ion}} + N^k_{\mathrm{electron}}) e^{-\beta (H^k + F^k)}\\ + &= 1/Z \sum_k (N^k_{\mathrm{ion}} + N^k_{\mathrm{electron}}) e^{-\beta (H^k + F^k)},\\ \end{aligned}\]

with the definitions:

\[\begin{aligned} N^k_{\mathrm{ion}} &= - \frac{\partial H^k}{\partial \mu} = \sum_i n^k_i\\ -N^k_{\mathrm{electron}} &= - \frac{\partial F^k}{\partial \mu} = \sum_i \left(1 + e^{\beta \epsilon^k_i}\right)^{-1}\\ +N^k_{\mathrm{electron}} &= - \frac{\partial F^k}{\partial \mu} = \sum_i \left(1 + e^{\beta \epsilon^k_i}\right)^{-1}.\\ \end{aligned}\]

Next Section: Particle-Hole Symmetry

diff --git a/_thesis/6_Appendices/A.1_Particle_Hole_Symmetry.html b/_thesis/6_Appendices/A.1_Particle_Hole_Symmetry.html index 812db38..11e1df2 100644 --- a/_thesis/6_Appendices/A.1_Particle_Hole_Symmetry.html +++ b/_thesis/6_Appendices/A.1_Particle_Hole_Symmetry.html @@ -66,16 +66,16 @@ image:

Particle-Hole Symmetry

The Hubbard and FK models on a bipartite lattice have particle-hole (PH) symmetry \(\mathcal{P}^\dagger H \mathcal{P} = - H\), accordingly they have symmetric energy spectra. The associated symmetry operator \(\mathcal{P}\) exchanges creation and annihilation operators along with a sign change between the two sublattices. In the language of the Hubbard model of electrons \(c_{\alpha,i}\) with spin \(\alpha\) at site \(i\) the particle hole operator corresponds to the substitution of new fermion operators \(d^\dagger_{\alpha,i}\) and number operators \(m_{\alpha,i}\) where

-

\[d^\dagger_{\alpha,i} = \epsilon_i c_{\alpha,i}\] \[m_{\alpha,i} = d^\dagger_{\alpha,i}d_{\alpha,i}\]

+

\[d^\dagger_{\alpha,i} = \epsilon_i c_{\alpha,i}\] \[m_{\alpha,i} = d^\dagger_{\alpha,i}d_{\alpha,i},\]

the lattices must be bipartite because to make this work we set \(\epsilon_i = +1\) for the A sublattice and \(-1\) for the even sublattice  [1].

The entirely filled state \(\ket{\Omega} = \sum_{\alpha,i} c^\dagger_{\alpha,i} \ket{0}\) becomes the new vacuum state \[d_{i\sigma} \ket{\Omega} = (-1)^i c^\dagger_{i\sigma} \sum_{j\rho} c^\dagger_{j\rho} \ket{0} = 0.\]

-

The number operator \(m_{\alpha,i} = 0,1\) now counts holes rather than electrons \[ m_{\alpha,i} = c_{\alpha,i} c^\dagger_{\alpha,i} = 1 - c^\dagger_{\alpha,i} c_{\alpha,i}.\]

-

In the case of nearest neighbour hopping on a bipartite lattice this transformation also leaves the hopping term unchanged because \(\epsilon_i \epsilon_j = -1\) when \(i\) and \(j\) are on different sublattices: \[ d^\dagger_{\alpha,i} d_{\alpha,j} = \epsilon_i \epsilon_j c_{\alpha,i} c^\dagger_{\alpha,j} = c^\dagger_{\alpha,i} c_{\alpha,j} \]

+

The number operator \(m_{\alpha,i} = 0,1\) now counts holes rather than electrons \[ m_{\alpha,i} = c^{\phantom{\dagger}}_{\alpha,i} c^\dagger_{\alpha,i} = 1 - c^\dagger_{\alpha,i} c^{\phantom{\dagger}}_{\alpha,i}.\]

+

In the case of nearest neighbour hopping on a bipartite lattice this transformation also leaves the hopping term unchanged because \(\epsilon_i \epsilon_j = -1\) when \(i\) and \(j\) are on different sublattices: \[ d^\dagger_{\alpha,i} d_{\alpha,j} = \epsilon_i \epsilon_j c^{\phantom{\dagger}}_{\alpha,i} c^\dagger_{\alpha,j} = c^\dagger_{\alpha,i} c^{\phantom{\dagger}}_{\alpha,j}. \]

Defining the particle density \(\rho\) as the number of fermions per site: \[ - \rho = \frac{1}{N} \sum_i \left( n_{i \uparrow} + n_{i \downarrow} \right) + \rho = \frac{1}{N} \sum_i \left( n_{i \uparrow} + n_{i \downarrow} \right). \]

-

The PH symmetry maps the Hamiltonian to itself with the sign of the chemical potential reversed and the density inverted about half filling: \[ \text{PH} : H(t, U, \mu) \rightarrow H(t, U, -\mu) \] \[ \rho \rightarrow 2 - \rho \]

-

The Hamiltonian is symmetric under PH at \(\mu = 0\) and so must all the observables, hence half filling \(\rho = 1\) occurs here. This symmetry and known observable acts as a useful test for the numerical calculations.

+

The PH symmetry maps the Hamiltonian to itself with the sign of the chemical potential reversed and the density inverted about half-filling: \[ \text{PH} : H(t, U, \mu) \rightarrow H(t, U, -\mu) \] \[ \rho \rightarrow 2 - \rho. \]

+

The Hamiltonian is symmetric under PH at \(\mu = 0\) and so must all the observables, hence half-filling \(\rho = 1\) occurs here. This symmetry and known observable acts as a useful test for the numerical calculations.

Next Section: Markov Chain Monte Carlo

diff --git a/_thesis/6_Appendices/A.2_Markov_Chain_Monte_Carlo.html b/_thesis/6_Appendices/A.2_Markov_Chain_Monte_Carlo.html index 2f69b80..307a9a8 100644 --- a/_thesis/6_Appendices/A.2_Markov_Chain_Monte_Carlo.html +++ b/_thesis/6_Appendices/A.2_Markov_Chain_Monte_Carlo.html @@ -92,21 +92,21 @@ image:

Markov Chain Monte Carlo (MCMC) is a useful method whenever we have a probability distribution that we want to sample from but there is not direct sampling way to do so.

Direct Random Sampling

-

In almost any computer simulation the ultimate source of randomness is a stream of (close to) uniform, uncorrelated bits generated from a pseudo random number generator. A direct sampling method takes such a source and outputs uncorrelated samples from the target distribution. The fact they’re uncorrelated is key as we’ll see later. Examples of direct sampling methods range from the trivial: take n random bits to generate integers uniformly between 0 and \(2^n\) to more complex methods such as inverse transform sampling and rejection sampling  [1].

+

In almost any computer simulation the ultimate source of randomness is a stream of (close to) uniform, uncorrelated bits generated from a pseudo random number generator. A direct sampling method takes such a source and outputs uncorrelated samples from the target distribution. The fact they are uncorrelated is key as we’ll see later. Examples of direct sampling methods range from the trivial: take n random bits to generate integers uniformly between 0 and \(2^n\) to more complex methods such as inverse transform sampling and rejection sampling  [1].

In physics the distribution we usually want to sample from is the Boltzmann probability over states of the system \(S\): \[ \begin{aligned} -p(S) &= \frac{1}{\mathcal{Z}} e^{-\beta H(S)} \\ +p(S) &= \frac{1}{\mathcal{Z}} e^{-\beta H(S)}, \\ \end{aligned} \] where \(\mathcal{Z} = \sum_S e^{-\beta H(S)}\) is the normalisation factor and ubiquitous partition function. In principle we could directly sample from this, for a discrete system there are finitely many choices. We could calculate the probability of each one and assign each a region of the unit interval which we could then sample uniformly from.

-

However if we actually try to do this we will run into two problems, we can’t calculate \(\mathcal{Z}\) for any reasonably sized systems because the state space grows exponentially with system size. Even if we could calculate \(\mathcal{Z}\), sampling from an exponentially large number of options quickly become tricky. This kind of problem happens in many other disciplines too, particularly when fitting statistical models using Bayesian inference  [2].

+

However, if we actually try to do this we will run into two problems, we can’t calculate \(\mathcal{Z}\) for any reasonably sized systems because the state space grows exponentially with system size. Even if we could calculate \(\mathcal{Z}\), sampling from an exponentially large number of options quickly become tricky. This kind of problem happens in many other disciplines too, particularly when fitting statistical models using Bayesian inference  [2].

MCMC Sampling

-

So what can we do? Well it turns out that if we’re willing to give up in the requirement that the samples be uncorrelated then we can use MCMC instead.

+

So what can we do? Well it turns out that if we are willing to give up in the requirement that the samples be uncorrelated then we can use MCMC instead.

MCMC defines a weighted random walk over the states \((S_0, S_1, S_2, ...)\), such that in the long time limit, states are visited according to their probability \(p(S)\).  [3–5].   [6]

-

\[\lim_{i\to\infty} p(S_i) = P(S)\]

+

\[\lim_{i\to\infty} p(S_i) = P(S).\]

In a physics context this lets us evaluate any observable with a mean over the states visited by the walk. \[\begin{aligned} -\langle O \rangle & = \sum_{S} p(S) \langle O \rangle_{S} = \sum_{i = 0}^{M} \langle O\rangle_{S_i} + \mathcal{O}(\tfrac{1}{\sqrt{M}})\\ +\langle O \rangle & = \sum_{S} p(S) \langle O \rangle_{S} = \sum_{i = 0}^{M} \langle O\rangle_{S_i} + \mathcal{O}(\tfrac{1}{\sqrt{M}}).\\ \end{aligned}\]

The samples in the random walk are correlated so the samples effectively contain less information than \(N\) independent samples would. As a consequence the variance is larger than the \(\langle O^2 \rangle - \langle O\rangle^2\) form it would have if the estimates were uncorrelated. Methods of estimating the true variance of \(\langle O \rangle\) and decided how many steps are needed will be considered later.

@@ -128,17 +128,17 @@ p(S) &= \frac{1}{\mathcal{Z}} e^{-\beta H(S)} \\

We can quite easily write down the properties that \(\mathcal{T}\) must have in order to yield the correct target distribution. Since we must transition somewhere at each step, we first have the normalisation condition that \[\sum\limits_S \mathcal{T}(S' \rightarrow S) = 1.\]

Second, let us move to an ensemble view, where rather than individual walkers and states, we think about the probability distribution of many walkers at each step. If we start all the walkers in the same place the initial distribution will be a delta function and as we step the walkers will wander around, giving us a sequence of probability distributions \(\{p_0(S), p_1(S), p_2(S)\ldots\}\). For discrete spaces we can write the action of the transition function on \(p_i\) as a matrix equation

\[\begin{aligned} -p_{i+1}(S) &= \sum_{S' \in \{S\}} p_i(S') \mathcal{T}(S' \rightarrow S) +p_{i+1}(S) &= \sum_{S' \in \{S\}} p_i(S') \mathcal{T}(S' \rightarrow S). \end{aligned}\]

This equation is essentially just stating that total probability mass is conserved as our walkers flow around the state space.

In order that \(p_i\) converges to our target distribution \(p\) in the long time limit, we need the target distribution to be a fixed point of the transition function

\[\begin{aligned} -P(S) &= \sum_{S'} P(S') \mathcal{T}(S' \rightarrow S) +P(S) &= \sum_{S'} P(S') \mathcal{T}(S' \rightarrow S). \end{aligned} \] Along with some more technical considerations such as ergodicity which won’t be considered here, global balance suffices to ensure that a MCMC method is correct  [7].

A sufficient but not necessary condition for global balance to hold is called detailed balance:

\[ -P(S) \mathcal{T}(S \rightarrow S') = P(S') \mathcal{T}(S' \rightarrow S) +P(S) \mathcal{T}(S \rightarrow S') = P(S') \mathcal{T}(S' \rightarrow S). \]

In practice most algorithms are constructed to satisfy detailed rather than global balance, though there are arguments that the relaxed requirements of global balance can lead to faster algorithms  [8].

The goal of MCMC is then to choose \(\mathcal{T}\) so that it has the desired thermal distribution \(P(S)\) as its fixed point and converges quickly onto it. This boils down to requiring that the matrix representation of \(T_{ij} = \mathcal{T}(S_i \to S_j)\) has an eigenvector with entries \(P_i = P(S_i)\) with eigenvalue 1 and all other eigenvalues with magnitude less than one. The convergence time depends on the magnitude of the second largest eigenvalue.

@@ -146,33 +146,33 @@ P(S) \mathcal{T}(S \rightarrow S') = P(S') \mathcal{T}(S' \rightarro

The Metropolis-Hastings Algorithm

-

The Metropolis-Hastings algorithm breaks the transition function into a proposal distribution \(q(S \to S')\) and an acceptance function \(\mathcal{A}(S \to S')\). \(q\) must be a function we can directly sample from, and in many cases takes the form of flipping some number of spins in \(S\), i.e if we’re flipping a single random spin in the spin chain, \(q(S \to S')\) is the uniform distribution on states reachable by one spin flip from \(S\). This also gives the symmetry property that \(q(S \to S') = q(S' \to S)\).

+

The Metropolis-Hastings algorithm breaks the transition function into a proposal distribution \(q(S \to S')\) and an acceptance function \(\mathcal{A}(S \to S')\). \(q\) must be a function we can directly sample from, and in many cases takes the form of flipping some number of spins in \(S\), i.e., if we are flipping a single random spin in the spin chain, \(q(S \to S')\) is the uniform distribution on states reachable by one spin flip from \(S\). This also gives the symmetry property that \(q(S \to S') = q(S' \to S)\).

The proposal \(S'\) is then accepted or rejected with an acceptance probability \(\mathcal{A}(S \to S')\), if the proposal is rejected then \(S_{i+1} = S_{i}\). Hence:

-

\[\mathcal{T}(S\to S') = q(S\to S')\mathcal{A}(S \to S')\]

+

\[\mathcal{T}(S\to S') = q(S\to S')\mathcal{A}(S \to S').\]

The Metropolis-Hasting algorithm is a slight extension of the original Metropolis algorithm which allows for non-symmetric proposal distributions \(q(S\to S') \neq q(S'\to S)\). It can be derived starting from detailed balance  [6]:

\[ -P(S)\mathcal{T}(S \to S') = P(S')\mathcal{T}(S' \to S) +P(S)\mathcal{T}(S \to S') = P(S')\mathcal{T}(S' \to S), \]

inserting the proposal and acceptance function

\[ -P(S)q(S \to S')\mathcal{A}(S \to S') = P(S')q(S' \to S)\mathcal{A}(S' \to S) +P(S)q(S \to S')\mathcal{A}(S \to S') = P(S')q(S' \to S)\mathcal{A}(S' \to S), \]

rearranging gives us a condition on the acceptance function in terms of the target distribution and the proposal distribution which can be thought of as inputs to the algorithm

\[ -\frac{\mathcal{A}(S \to S')}{\mathcal{A}(S' \to S)} = \frac{P(S')q(S' \to S)}{P(S)q(S \to S')} = f(S, S') +\frac{\mathcal{A}(S \to S')}{\mathcal{A}(S' \to S)} = \frac{P(S')q(S' \to S)}{P(S)q(S \to S')} = f(S, S'). \]

The Metropolis-Hastings algorithm is the choice

\[ \begin{aligned} \label{eq:mh} -\mathcal{A}(S \to S') = \min\left(1, f(S,S')\right) +\mathcal{A}(S \to S') = \min\left(1, f(S,S')\right), \end{aligned} \] for the acceptance function. The proposal distribution is left as a free choice.

-

Noting that \(f(S,S') = 1/f(S',S)\), we can see that the MH algorithm satifies detailed balance by considering the two cases \(f(S,S') > 1\) and \(f(S,S') < 1\).

+

Noting that \(f(S,S') = 1/f(S',S)\), we can see that the MH algorithm satisfies detailed balance by considering the two cases \(f(S,S') > 1\) and \(f(S,S') < 1\).

By choosing the proposal distribution such that \(f(S,S')\) is as close as possible to one, the rate of rejections can be reduced and the algorithm sped up. This can be challenging though, as getting \(f(S,S')\) close to 1 would imply that we can directly sample from a distribution very close to the target distribution. As MCMC is usually applied to problems for which there is virtually no hope of sampling directly from the target distribution, it’s rare that one can do so approximately.

When the proposal distribution is symmetric as ours is, it cancels out in the expression for the acceptance function and the Metropolis-Hastings algorithm is simply the choice

\[ -\mathcal{A}(S \to S') = \min\left(1, e^{-\beta\;\Delta F}\right) +\mathcal{A}(S \to S') = \min\left(1, e^{-\beta\;\Delta F}\right), \]

where \(F\) is the overall free energy of the system, including both the quantum and classical sector.

To implement the acceptance function in practice we pick a random number in the unit interval and accept if it is less than \(e^{-\beta\;\Delta F}\):

@@ -207,11 +207,11 @@ P(S)q(S \to S')\mathcal{A}(S \to S') = P(S')q(S' \to S)\mathcal{ states[i] = current_state

As discussed in the main text, for the model parameters used, we find that with our new scheme the matrix diagonalisation is skipped around 30% of the time at \(T = 2.5\) and up to 80% at \(T = 1.5\). We observe that for \(N = 50\), the matrix diagonalisation, if it occurs, occupies around 60% of the total computation time for a single step. This rises to 90% at N = 300 and further increases for larger N. We therefore get the greatest speedup for large system sizes at low temperature where many prospective transitions are rejected at the classical stage and the matrix computation takes up the greatest fraction of the total computation time. The upshot is that we find a speedup of up to a factor of 10 at the cost of very little extra algorithmic complexity.

-

This modified scheme has the acceptance function \[\mathcal{A}(a \to b) = \min\left(1, e^{-\beta \Delta H_s}\right)\min\left(1, e^{-\beta \Delta F_c}\right)\;.\]

+

This modified scheme has the acceptance function \[\mathcal{A}(a \to b) = \min\left(1, e^{-\beta \Delta H_s}\right)\min\left(1, e^{-\beta \Delta F_c}\right).\]

We can show that this satisfies the detailed balance equations as follows. Defining \(r_c = e^{-\beta H_c}\) and \(r_q = e^{-\beta F_q}\) our target distribution is \(\pi(a) = r_c r_q\). This method has \(\mathcal{T}(a\to b) = q(a\to b)\mathcal{A}(a \to b)\) with symmetric \(p(a \to b) = \pi(b \to a)\) and \(\mathcal{A} = \min\left(1, r_c\right) \min\left(1, r_q\right)\)

-

Substituting this into the detailed balance equation gives: \[\mathcal{T}(a \to b)/\mathcal{T}(b \to a) = \pi(b)/\pi(a) = r_c r_q\]

+

Substituting this into the detailed balance equation gives: \[\mathcal{T}(a \to b)/\mathcal{T}(b \to a) = \pi(b)/\pi(a) = r_c r_q.\]

Taking the LHS and substituting in our transition function: \[\begin{aligned} -\mathcal{T}(a \to b)/\mathcal{T}(b \to a) = \frac{\min\left(1, r_c\right) \min\left(1, r_q\right)}{ \min\left(1, 1/r_c\right) \min\left(1, 1/r_q\right)}\end{aligned}\]

+\mathcal{T}(a \to b)/\mathcal{T}(b \to a) = \frac{\min\left(1, r_c\right) \min\left(1, r_q\right)}{ \min\left(1, 1/r_c\right) \min\left(1, 1/r_q\right)},\end{aligned}\]

which simplifies to \(r_c r_q\) as \(\min(1,r)/\min(1,1/r) = r\) for \(r > 0\).

Autocorrelation Time

@@ -219,9 +219,9 @@ P(S)q(S \to S')\mathcal{A}(S \to S') = P(S')q(S' \to S)\mathcal{ Figure 1: (Upper) 10 MCMC chains starting from the same initial state for a system with N = 150 sites and 3000 MCMC steps. At each MCMC step, n spins are flipped where n is drawn from Uniform(1,N) and this is repeated N^2/100 times. The simulations therefore have the potential to necessitate 10*N^2 matrix diagonalisations for each 100 MCMC steps. (Lower) The normalised autocorrelation (\langle m_i m_{i-j}\rangle - \langle m_i\rangle \langle m_i \rangle) / Var(m_i)) averaged over i. It can be seen that even with each MCMC step already being composed of many individual flip attempts, the autocorrelation is still non negligible and must be taken into account in the statistics. t = 1, \alpha = 1.25, T = 2.2, J = U = 5 -

At this stage one might think we’re done. We can indeed draw independent samples from our target Boltzmann distribution by starting from some arbitrary initial state and doing \(k\) steps to arrive at a sample. These are not, however, independent samples. In fig. 1 it is already clear that the samples of the order parameter \(m\) have some autocorrelation because only a few spins are flipped each step. Even when the number of spins flipped per step is increased that it can be an important effect near the phase transition. Let’s define the autocorrelation time \(\tau(O)\) informally as the number of MCMC samples of some observable O that are statistically equal to one independent sample or equivalently as the number of MCMC steps after which the samples are correlated below some cut-off, see  [9]. The autocorrelation time is generally shorter than the convergence time so it therefore makes sense from an efficiency standpoint to run a single walker for many MCMC steps rather than to run a huge ensemble for \(k\) steps each.

+

At this stage one might think we are done. We can indeed draw independent samples from our target Boltzmann distribution by starting from some arbitrary initial state and doing \(k\) steps to arrive at a sample. These are not, however, independent samples. In fig. 1 it is already clear that the samples of the order parameter \(m\) have some autocorrelation because only a few spins are flipped each step. Even when the number of spins flipped per step is increased that it can be an important effect near the phase transition. Let’s define the autocorrelation time \(\tau(O)\) informally as the number of MCMC samples of some observable O that are statistically equal to one independent sample or equivalently as the number of MCMC steps after which the samples are correlated below some cut-off, see ref.  [9]. The autocorrelation time is generally shorter than the convergence time so it therefore makes sense from an efficiency standpoint to run a single walker for many MCMC steps rather than to run a huge ensemble for \(k\) steps each.

Once the random walk has been carried out for many steps, the expectation values of \(O\) can be estimated from the MCMC samples \(S_i\): \[ - \langle O \rangle = \sum_{i = 0}^{N} O(S_i) + \mathcal{O}(\frac{1}{\sqrt{N}}) + \langle O \rangle = \sum_{i = 0}^{N} O(S_i) + \mathcal{O}(\frac{1}{\sqrt{N}}). \]

The samples are correlated so the N of them effectively contains less information than \(N\) independent samples would, in fact roughly \(N/\tau\) effective samples. As a consequence the variance is larger than the \(\langle O^2 \rangle - \langle O \rangle ^2\) form it would have if the estimates were uncorrelated. There are many methods in the literature for estimating the true variance of \(\langle O \rangle\) and deciding how many steps are needed but my approach has been to run a small number of parallel chains, which are independent, in order to estimate the statistical error produced. This is a slightly less computationally efficient because it requires throwing away those \(k\) steps generated before convergence multiple times but it is conceptually simple.

@@ -239,7 +239,7 @@ P(S)q(S \to S')\mathcal{A}(S \to S') = P(S')q(S' \to S)\mathcal{
  • Choosing n from Uniform(1, N) and then flipping n sites for some fixed N.
  • Attempting to tune the proposal distribution for each parameter regime.
    • -

    Fro fig. 2 we see that even at moderately high temperatures \(T > T_c\) flipping one or two sites is the best choice. However for some simulations at very high temperature flipping more spins is warranted.

    +

    Fro fig. 2 we see that even at moderately high temperatures \(T > T_c\) flipping one or two sites is the best choice. However, for some simulations at very high temperature flipping more spins is warranted.

    Next Section: Lattice Generation

    diff --git a/_thesis/6_Appendices/A.3_Lattice_Generation.html b/_thesis/6_Appendices/A.3_Lattice_Generation.html index 7b94432..fd84d74 100644 --- a/_thesis/6_Appendices/A.3_Lattice_Generation.html +++ b/_thesis/6_Appendices/A.3_Lattice_Generation.html @@ -78,7 +78,7 @@ image:

    Three key pieces of information allow us to represent amorphous lattices. The majority of the graph connectivity is encoded by an ordered list of edges \((i,j)\). These are ordered to represent both directed and undirected graphs. This is useful for defining the sign of bond operators \(u_{ij} = - u_{ji}\).

    Information about the embedding of the lattice onto the torus is encoded into a point on the unit square associated with each vertex. The torus is unwrapped onto the square by defining an arbitrary pair of cuts along the major and minor axes. For simplicity, we take these axes to be the lines \(x = 0\) and \(y = 0\). We can wrap the unit square back up into a torus by identifying the lines \(x = 0\) with \(x = 1\) and \(y = 0\) with \(y = 1\).

    Finally, we need to encode the topology of the graph. This is necessary because, if we are simply given an edge \((i, j)\) we do not know how the edge gets from vertex i to vertex j. One method would be taking the shortest path, but it could also ‘go the long way around’ by crossing one of the cuts. To encode this information, we store an additional vector \(\vec{r}\) associated with each edge. \(r_i^x = 0\) means that edge i does not cross the x. \(r_i^x = +1\) (\(-1\)) means it crossed the cut in a positive (negative) sense.

    -

    This description of the lattice has a very nice relationship to Bloch’s theorem. Applying Bloch’s theorem to a periodic lattice essentially means wrappping the unit cell onto a torus. Variations that happen at longer length scales than the size of the unit cell are captured by the crystal momentum. The crystal momentum inserts a phase factor \(e^{i \vec{q}\cdot\vec{r}}\) onto bonds that cross to adjacent unit cells. The vector \(\vec{r}\) is exactly what we use to encode the topology of our lattices.

    +

    This description of the lattice has a very nice relationship to Bloch’s theorem. Applying Bloch’s theorem to a periodic lattice essentially means wrapping the unit cell onto a torus. Variations that happen at longer length scales than the size of the unit cell are captured by the crystal momentum. The crystal momentum inserts a phase factor \(e^{i \vec{q}\cdot\vec{r}}\) onto bonds that cross to adjacent unit cells. The vector \(\vec{r}\) is exactly what we use to encode the topology of our lattices.

    Figure 1: Bloch’s theorem can be thought of as transforming from a periodic Hamiltonian on the plane to the unit cell defined on a torus. In addition we get some phase factors e^{i\vec{k}\cdot\vec{r}} associated with bonds that cross unit cells that depend on the sense in which they do so \vec{r} = (\pm1, \pm1). Representing graphs on the torus turns out to require a similar idea, we unwrap the torus to one unit cell and keep track of which bonds cross the cell boundaries. @@ -91,7 +91,7 @@ image:

    In the main text we discuss the problem of three-edge-colouring, assigning one of three labels to each edge of a graph such that no edges of the same label meet at a vertex. To solve this in practice I use a solver called MiniSAT  [1]. Like most modern SAT solvers, MiniSAT requires the input problem to be specified in Conjunctive Normal Form (CNF). CNF requires that the constraints be encoded as a set of clauses of the form

    -

    \[x_1 \;\textrm{or}\; -x_3 \;\textrm{or}\; x_5\]

    +

    \[x_1 \;\textrm{or}\; -x_3 \;\textrm{or}\; x_5,\]

    that contain logical ORs of some subset of the variables where any of the variables may also be logically NOT’d, which we represent by negation here. A solution of the problem is one that makes all the clauses simultaneously true.

    I encode the edge colouring problem by assigning \(3B\) boolean variables to each of the \(B\) edges of the graph, \(x_{i\alpha}\) where \(x_{i\alpha} = 1\) indicates that edge \(i\) has colour \(\alpha\). For edge colouring graphs we need two types of constraints: 1. Each edge is exactly one colour. 2. No neighbouring edges are the same colour.

    The first constraint is a product of doing this mapping to boolean variables. The solver does not know anything about the structure of the problem unless it is encoded into the variables. Let’s say we have three variables that correspond to particular edge being red \(r\), green \(g\) or blue \(b\). To require that exactly one of the variables be true, we can enforce that no pair of variables be true: -(r and b) -(r and g) -(b and g)

    diff --git a/_thesis/6_Appendices/A.4_Lattice_Colouring.html b/_thesis/6_Appendices/A.4_Lattice_Colouring.html index b6d00fd..557e11b 100644 --- a/_thesis/6_Appendices/A.4_Lattice_Colouring.html +++ b/_thesis/6_Appendices/A.4_Lattice_Colouring.html @@ -74,9 +74,9 @@ image:

    A four-face-colouring can be converted into a three-edge-colouring quite easily: 1. Assume the faces of G can be four-coloured with labels (0,1,2,3) 2. Label each edge of G according to \(i + j \;\textrm{mod}\; 3\) where i and j are the labels of the face adjacent to that edge. For each edge label there are two face label pairs that do not share any face labels. i,e the edge label \(0\) can come about either from faces \(0 + 3\) or \(1 + 2\).

    Explicitly, the mapping from face labels to edge labels is:

    \[\begin{aligned} -0 + 3 \;\mathrm{or}\; 1 + 2 &= 0 \;\mathrm{mod}\; 3\\ -0 + 1 \;\mathrm{or}\; 2 + 3 &= 1 \;\mathrm{mod}\; 3\\ -0 + 2 \;\mathrm{or}\;1 + 3 &= 2 \;\mathrm{mod}\; 3\\ +0 + 3 \;\mathrm{or}\; 1 + 2 &= 0 \;\mathrm{mod}\; 3,\\ +0 + 1 \;\mathrm{or}\; 2 + 3 &= 1 \;\mathrm{mod}\; 3,\\ +0 + 2 \;\mathrm{or}\;1 + 3 &= 2 \;\mathrm{mod}\; 3.\\ \end{aligned} \]

      diff --git a/_thesis/6_Appendices/A.5_The_Projector.html b/_thesis/6_Appendices/A.5_The_Projector.html index ab8b0df..f4182c2 100644 --- a/_thesis/6_Appendices/A.5_The_Projector.html +++ b/_thesis/6_Appendices/A.5_The_Projector.html @@ -71,19 +71,19 @@ image:
      The relationship between the different Hilbert spaces used in the solution. __needs updating__
      -->

      The physical states are defined as those for which \(D_i |\phi\rangle = |\phi\rangle\) for all \(D_i\). Since \(D_i\) has eigenvalues \(\pm1\), the quantity \(\tfrac{(1+D_i)}{2}\) has eigenvalue \(1\) for physical states and \(0\) for extended states so is the local projector onto the physical subspace.

      -

      Therefore, the global projector is \[ \mathcal{P} = \prod_{i=1}^{2N} \left( \frac{1 + D_i}{2}\right)\]

      -

      for a toroidal trivalent lattice with \(N\) plaquettes \(2N\) vertices and \(3N\) edges. As discussed earlier, the product over \((1 + D_j)\) can also be thought of as the sum of all possible subsets \(\{i\}\) of the \(D_j\) operators, which is the set of all possible gauge symmetry operations.

      -

      \[ \mathcal{P} = \frac{1}{2^{2N}} \sum_{\{i\}} \prod_{i\in\{i\}} D_i\]

      +

      Therefore, the global projector is \[ \mathcal{P} = \prod_{i=1}^{2N} \left( \frac{1 + D_i}{2}\right),\]

      +

      for a toroidal trivalent lattice with \(N\) plaquettes \(2N\) vertices and \(3N\) edges. As discussed earlier, the product over \((1 + D_j)\) can also be thought of as the sum of all possible subsets \(\{i\}\) of the \(D_j\) operators, which is the set of all possible gauge symmetry operations

      +

      \[ \mathcal{P} = \frac{1}{2^{2N}} \sum_{\{i\}} \prod_{i\in\{i\}} D_i.\]

      Since the gauge operators \(D_j\) commute and square to one, we can define the complement operator \(C = \prod_{i=1}^{2N} D_i\) and see that it takes each set of \(\prod_{i \in \{i\}} D_j\) operators and gives us the complement of that set. We will shortly see why \(C\) is the identity in the physical subspace, as noted earlier.

      We use the complement operator to rewrite the projector as a sum over half the subsets of \(\{i\}\) - referred to as \(\Lambda\). The complement operator deals with the other half

      -

      \[ \mathcal{P} = \left( \frac{1}{2^{2N-1}} \sum_{\Lambda} \prod_{i\in\{i\}} D_i\right) \left(\frac{1 + \prod_i^{2N} D_i}{2}\right) = \mathcal{S} \cdot \mathcal{P}_0\]

      +

      \[ \mathcal{P} = \left( \frac{1}{2^{2N-1}} \sum_{\Lambda} \prod_{i\in\{i\}} D_i\right) \left(\frac{1 + \prod_i^{2N} D_i}{2}\right) = \mathcal{S} \cdot \mathcal{P}_0.\]

      To compute \(\mathcal{P}_0\), the main quantity needed is the product of the local projectors \(D_i\) \[\prod_i^{2N} D_i = \prod_i^{2N} b^x_i b^y_i b^z_i c_i \] for a toroidal trivalent lattice with \(N\) plaquettes \(2N\) vertices and \(3N\) edges.

      -

      First, we reorder the operators by bond type. This does not require any information about the underlying lattice.

      -

      \[\prod_i^{2N} D_i = \prod_i^{2N} b^x_i \prod_i^{2N} b^y_i \prod_i^{2N} b^z_i \prod_i^{2N} c_i\]

      -

      The product over \(c_i\) operators reduces to a determinant of the Q matrix and the fermion parity, see  [1]. The only difference from the honeycomb case is that we cannot explicitly compute the factors \(p_x,p_y,p_z = \pm\;1\) that arise from reordering the b operators such that pairs of vertices linked by the corresponding bonds are adjacent.

      -

      \[\prod_i^{2N} b^\alpha_i = p_\alpha \prod_{(i,j)}b^\alpha_i b^\alpha_j\]

      +

      First, we reorder the operators by bond type. This does not require any information about the underlying lattice,

      +

      \[\prod_i^{2N} D_i = \prod_i^{2N} b^x_i \prod_i^{2N} b^y_i \prod_i^{2N} b^z_i \prod_i^{2N} c_i.\]

      +

      The product over \(c_i\) operators reduces to a determinant of the Q matrix and the fermion parity, see ref.  [1]. The only difference from the honeycomb case is that we cannot explicitly compute the factors \(p_x,p_y,p_z = \pm\;1\) that arise from reordering the b operators such that pairs of vertices linked by the corresponding bonds are adjacent,

      +

      \[\prod_i^{2N} b^\alpha_i = p_\alpha \prod_{(i,j)}b^\alpha_i b^\alpha_j.\]

      However, they are simply the parity of the permutation from one ordering to the other and can be computed in linear time with a cycle decomposition  [2].

      -

      We find that \[\mathcal{P}_0 = 1 + p_x\;p_y\;p_z\; \hat{\pi} \; \mathrm{det}(Q^u) \; \prod_{\{i,j\}} -iu_{ij}\]

      +

      We find that \[\mathcal{P}_0 = 1 + p_x\;p_y\;p_z\; \hat{\pi} \; \mathrm{det}(Q^u) \; \prod_{\{i,j\}} -iu_{ij},\]

      where \(p_x\;p_y\;p_z = \pm 1\) are lattice structure factors and \(\mathrm{det}(Q^u)\) is the determinant of the matrix mentioned earlier that maps \(c_i\) operators to normal mode operators \(b'_i, b''_i\). These depend only on the lattice structure.

      \(\hat{\pi} = \prod{i}^{N} (1 - 2\hat{n}_i)\) is the parity of the particular many-body state determined by fermionic occupation numbers \(n_i\). As discussed in  [1], \(\hat{\pi}\) is gauge invariant in the sense that \([\hat{\pi}, D_i] = 0\).

      This implies that \(det(Q^u) \prod -i u_{ij}\) is also a gauge invariant quantity. In translation invariant models this quantity which can be related to the parity of the number of vortex pairs in the system  [3].

      diff --git a/_thesis/toc.html b/_thesis/toc.html index df03653..001ec0f 100644 --- a/_thesis/toc.html +++ b/_thesis/toc.html @@ -1,7 +1,7 @@