Quantum spin liquid and magnetic order in a two-dimensional non-symmorphic lattice: considering the distorted Kagome lattice of Volborthite

TL;DR

This paper investigates quantum spin liquids and magnetic orders in a non-symmorphic lattice model inspired by Volborthite, revealing potential phases and signatures relevant for experimental detection.

Contribution

It classifies possible spin liquid and magnetic phases in a non-symmorphic lattice using PSG and Schwinger boson methods, specific to Volborthite's anisotropic structure.

Findings

Quantum spin liquids connect to incommensurate spiral order.

Spin liquids exhibit symmetry fractionalization signatures.

Magnetic phases include incommensurate spiral and commensurate spin density wave.

Abstract

The Kagome-lattice-based material, Volborthite, $\mathrm{Cu_3 V_2 O_7 (OH)_2 \cdot 2 H_2 O}$, has been considered as a promising platform for discovery of unusual quantum ground states due to the frustrated nature of spin interactions. Here we explore possible quantum spin liquid and magnetically ordered phases in a two-dimensional non-symmorphic lattice described by $p2gg$ layer space group, which is consistent with the spatial anisotropy of the spin model derived from density functional theory (DFT) for Volborthite. Using the projective symmetry group (PSG) analysis and Schwinger boson mean field theory, we classify possible spin liquid phases with bosonic spinons and investigate magnetically ordered phases connected to such states. It is shown, in general, that only translationally invariant mean field states are allowed in two-dimensional non-symmorphic lattices, which simplifies the classification considerably. The mean field phase diagram of the DFT-derived spin model is studied and it is found that possible quantum spin liquid phases are connected to two types of magnetically ordered phases, a coplanar incommensurate $(q,0)$ spiral order as the ground state and a closely competing coplanar commensurate $(\pi,\pi)$ spin density wave order. In addition, periodicity enhancement of the two-spinon continuum, a signature of symmetry fractionalization, is found in the spin liquid phases connected to the $(\pi,\pi)$ spin density wave order. We discuss relevance of these results to recent and future experiments on Volborthite.

Tables

Table 1. Table 1: Critical bosonic density κ c subscript 𝜅 𝑐 \kappa_{c} of different spin liquid states ( p 2 , p 3 ) subscript 𝑝 2 subscript 𝑝 3 \left(p_{2},p_{3}\right) .

$(p_2,p_3)$	${\kappa }_c$
$(0,0)$	$0.327$
$(1,0)$	$0.351$
$(0,1)$	$0.362$
$(1,1)$	$0.368$

Table 2. Table 2: Comparison of total energy E 𝐸 E (scaled by the total number of sites N 𝑁 N ) of four spin liquid states ( p 2 , p 3 ) subscript 𝑝 2 subscript 𝑝 3 (p_{2},p_{3}) at several κ < κ c 𝜅 subscript 𝜅 𝑐 \kappa<\kappa_{c} . The energy shown here is in unit of 10 − 2 ​ J superscript 10 2 𝐽 10^{-2}J , where J 𝐽 J is the strongest coupling constant. In low κ 𝜅 \kappa limit, the states with same p 2 subscript 𝑝 2 p_{2} has exactly the same spectrum and energy. The energy begins to differ around κ = 0.25 𝜅 0.25 \kappa=0.25 . ( 1 , 0 ) 1 0 (1,0) is identified to be the most energetically favorable state.

$\kappa \backslash (p_2,p_3)$	$(0,0)$	$(0,1)$	$(1,0)$	$(1,1)$
$0.05$	-0.951587	-0.951587	-0.953141	-0.953141
$0.10$	-1.973095	-1.973095	-1.984894	-1.984894
$0.15$	-3.064663	-3.064663	-3.101634	-3.101634
$0.20$	-4.226422	-4.226422	-4.307044	-4.307044
$0.25$	-5.458552	-5.458552	-5.603122	-5.603122
$0.26$	-5.713968	-5.714021	-5.873327	-5.873327
$0.27$	-5.972844	-5.973016	-6.147242	-6.147242
$0.28$	-6.235185	-6.235556	-6.425370	-6.425271
$0.29$	-6.500999	-6.501673	-6.707870	-6.707606
$0.30$	-6.770304	-6.771429	-6.994841	-6.994180
$0.31$	-7.043122	-7.044775	-7.286310	-7.285185
$0.32$	-7.319378	-7.321825	-7.582275	-7.580489
$0.33$	$-$	-7.602579	-7.882738	-7.880159
$0.34$	$-$	-7.887169	-8.187831	-8.184325
$0.35$	$-$	-8.175661	$-$	-8.492923

Equations

$$H=\sum _{ij}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

$${\mathbf{S}}_i=\frac{1}{2}\sum _{\alpha \beta }{b}_{i\alpha }^{†}{\mathbf{\sigma }}_{\alpha \beta }{b}_{i\beta },$$

$$\left[b_{i\alpha },b_{j\beta }^{†}\right]$$

$$\left[b_{i\alpha },b_{j\beta }\right]$$

$${\hat{n}}_i=\sum _{\alpha }{b}_{i\alpha }^{†}{b}_{i\alpha }.$$

$${\mathbf{S}}_i^2$$

$$=\frac{{\hat{n}}_i}{2}\left(\frac{{\hat{n}}_i}{2}+1\right),$$

$${\hat{n}}_i=2S\equiv \kappa .$$

$$⟨{\hat{n}}_i⟩=\kappa .$$

$${\hat{A}}_{ij}$$

$${\hat{B}}_{ij}$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_j=:{\hat{B}}_{ij}^{†}{\hat{B}}_{ij}:-{\hat{A}}_{ij}^{†}{\hat{A}}_{ij}.$$

$$H=\sum _{ij}{J}_{ij}\left(:{\hat{B}}_{ij}^{†}{\hat{B}}_{ij}:-{\hat{A}}_{ij}^{†}{\hat{A}}_{ij}\right)-\sum _i{\mu }_i\left({\hat{n}}_i-\kappa \right).$$

$$H_{\mathrm{MF}}=$$

$$-\left(A_{ij}^{\ast }{\hat{A}}_{ij}+h.c.-{\left|A_{ij}\right|}^2\right)]$$

$$-\sum _i{\mu }_i\left({\hat{n}}_i-\kappa \right).$$

$$\frac{\partial ⟨H_{\mathrm{MF}}⟩}{\partial {\mathcal{O}}_{ij}}=0⟺{\mathcal{O}}_{ij}=⟨{\hat{\mathcal{O}}}_{ij}⟩,{\mathcal{O}}_{ij}=A_{ij},B_{ij},$$

$$\frac{\partial ⟨H_{\mathrm{MF}}⟩}{\partial {\mu }_i}=0⟺\kappa =⟨{\hat{n}}_i⟩.$$

$$:{\hat{B}}_{ij}^{†}{\hat{B}}_{ij}:+{\hat{A}}_{ij}^{†}{\hat{A}}_{ij}=\frac{1}{4}{\hat{n}}_i\left({\hat{n}}_j-{\delta }_{ij}\right),$$

$${\mathbf{S}}_i\cdot {\mathbf{S}}_j$$

$$=S^2-2{\hat{A}}_{ij}^{†}{\hat{A}}_{ij},$$

$$H_{\mathrm{MF}}=$$

$$+\sum _{J_{ij}<0}{J}_{ij}\left(B_{ij}^{\ast }{\hat{B}}_{ij}+h.c.-{\left|B_{ij}\right|}^2\right)$$

$$-\sum _i{\mu }_i\left({\hat{n}}_i-\kappa \right)$$

$$\frac{1}{4}\sum _{s=1,2,4,5}⟨{\hat{n}}_{i,s}⟩$$

$$\frac{1}{2}\sum _{s=3,6}⟨{\hat{n}}_{i,s}⟩$$

$$G:$$

$$A_{ij}⟶e^{i\left(\varphi (i)+\varphi (j)\right)}{A}_{ij},$$

$$B_{ij}⟶e^{i\left(-\varphi (i)+\varphi (j)\right)}{B}_{ij}.$$

$$X:b_{i\alpha }⟶b_{X(i)\alpha }.$$

$$A_{ij}{b}_{i\alpha }^{†}{b}_{j\beta }^{†}$$

$$\stackrel{G_X}{⟶}{A}_{ij}{e}^{-i\left({\varphi }_X(X(i))+{\varphi }_X(X(j))\right)}{b}_{X(i)\alpha }^{†}{b}_{X(j)\beta }^{†}.$$

$$A_{ij}{e}^{-i\left({\varphi }_X(X(i))+{\varphi }_X(X(j))\right)}=A_{X(i)X(j)},$$

$$B_{ij}{e}^{-i\left(-{\varphi }_X(X(i))+{\varphi }_X(X(j))\right)}=B_{X(i)X(j)}.$$

$${\varphi }_X(\mathbf{r})⟶\varphi (\mathbf{r})+{\varphi }_X(\mathbf{r})-\varphi (X^{-1}(\mathbf{r})).$$

$$T_x^{-1}{T}_y^{-1}{T}_x{T}_y=I.$$

Full text

1. I Introduction

2. II Lattice Structure and Spin Model

3. II.1 Unit Cell and Space Group

4. II.2 Microscopic Spin Model

5. III Schwinger Boson Mean Field Theory

6. IV Projective Symmetry Group Analysis

7. IV.1 Algebraic PSG

8. IV.2 Mean Field Ansatz

9. V Emergent Quantum Phases and Their Physical Properties

10. V.1 Spin Liquid States

11. V.2 Magnetically Ordered States

12. V.3 Periodicity Enhancement of Two-Spinon Continuum

13. VI Discussion

14. A Space Group and Algebraic Identities of Non-Symmorphic Kagome Lattice

15. B Solution to Bosonic PSG

16. B.1 Time Reversal Symmetry in Bosonic PSG

17. C Mean Field Calculation Techniques

18. C.1 Bogoliubov Transformation

19. C.2 Simulated Annealing

20. D Linear Approximation Scheme

21. E Energetics of Spin Liquid States

I Introduction

The interplay between lattice structure and spin exchange interactions is the defining characteristic of frustrated quantum magnets that hold the promise for exotic quantum ground states such as quantum spin liquid phases.Balents (2010) A prominent example of frustrated magnets is the Kagome lattice system with localized $S=1/2$ moments. There exist several materials that host $S=1/2$ local moments on various versions of the Kagome lattice. For example, a number of experimentsMendels et al. (2007); Helton et al. (2007) observe signatures of a possible spin liquid ground state in the isotropic Kagome antiferromagnet,Shores et al. (2005) Herbertsmithite $\mathrm{ZnCu_{3}(OH)_{6}Cl_{2}}$. Another material, Volborthite $\mathrm{Cu_{3}V_{2}O_{7}(OH)_{2}\cdot 2H_{2}O}$, has a distorted Kagome lattice Hiroi et al. (2001); Yoshida et al. (2012) and a magnetic order appears below 1K. Clearly, the difference in lattice structure plays an important role in the determination of the quantum ground state.

Volborthite, however, shows a rather complex response to an external magnetic field and its full magnetic phase diagram is still under investigation.Yoshida et al. (2009a, b); Okamoto et al. (2011); Ishikawa et al. (2015); Yoshida et al. The small energy scale for the magnetic order and unusual response to an external magnetic field make it difficult to identify the magnetic order below $1\,\mathrm{K}$ and may also suggest that a number of exotic quantum ground states may arise in this system. On the other hand, a recent thermal conductivity measurement finds signatures of entropy-carrying charge-neutral excitations above $1\,\mathrm{K}$, which suggests that the phase above $1\,\mathrm{K}$ may be connected to a putative quantum spin liquid with spinon quasiparticles.Watanabe et al. Over the years, Volborthite has also inspired several theoretical studies.Yavors’kii et al. (2007); Wang et al. (2007); Schnyder et al. (2008); Kaneko et al. (2010); Parker and Balents

In this work, we investigate possible quantum spin liquid and magnetically ordered phases on the distorted Kagome lattice, in view of the spatially anisotropic spin exchange model proposed earlier by the density functional theory (DFT) computations.Janson et al. (2016) The distorted Kagome lattice corresponding to the so-called $J-J^{\prime}-J_{1}-J_{2}$ model represents the non-symmorphic layer space group, $p2gg$, that possesses glide symmetry, the combination of reflection and a fractional translation. Motivated by the experiments and the DFT result, we first classify possible spin liquid phases with bosonic spinons for the $p2gg$ space group via the projective symmetry group (PSG) analysis Wen (2002) of the Schwinger boson Sachdev (1992); Wang and Vishwanath (2006); Messio et al. (2013) mean field states. It is shown that only the translationally invariant mean field states are possible in two-dimensional non-symmorphic lattices. As a result, there are only eight possible spin liquid phases characterized by three $\mathbb{Z}_{2}$ variables when both of the spatial and time reversal symmetries are taken into account. Using these results, we investigate the mean field ground states of the $J-J^{\prime}-J_{1}-J_{2}$ model. We identify the stable quantum spin liquid states and obtain magnetically ordered phases that arise from these spin liquid states by condensing the bosonic spinons.

It is found that the spin liquid ground state of this model is connected to a coplanar incommensurate $(q,0)$ spiral order in the semiclassical limit while a different, but highly competing, spin liquid state is related to a coplanar commensurate $(\pi,\pi)$ spin density wave phase, where the amplitude of the spin density varies from site to site. The corresponding classical model is also studied by simulated annealing and we find that the classical magnetic order is indeed consistent with the $(q,0)$ spiral order for the model parameters determined by the DFT computations. This correspondence is natural as the length of classical spin is fixed and hence the competing spin density wave order is simply not possible in the classical model. Interestingly, the spin liquid state related to the $(\pi,\pi)$ spin density wave order exhibits periodicity enhancement, namely an extra periodicity beyond the normal periodicity given by the lattice structure, in the two-spinon continuum. Such periodicity enhancement is a signature of the so-called symmetry fractionalization,Essin and Hermele (2013, 2014); Lu et al. which originates from the existence of spatial inversion and time reversal symmetry.

We argue that it may be necessary to consider both spin liquid states described above on equal footing when we apply these results to experiments. This is because the mean field energetics may not be accurate enough to determine the true ground state and the model parameters determined by the DFT computations may also allow some variations. In fact, both of the incommensurate $(q.0)$ spiral order and the commensurate $(\pi,\pi)$ spin density wave order may be compatible with the existing nuclear magnetic resonance (NMR) experimental data on Volborthite below 1K. An important question is whether the paramagnetic state above 1K can be regarded as a continuation of certain quantum spin liquid state, as suggested by the recent thermal conductivity measurement. In view of this possibility, it will be interesting to perform a neutron scattering experiment on Volborthite above 1K and look for signatures of the two-spinon continuum and especially periodicity enhancement of such continuum of excitations.

The rest of this paper is organized as follows. In Section II, we introduce the lattice structure that is compatible with the spin model derived from DFT computations. We explain the non-symmorphic nature of this lattice and some details of the DFT-derived spin model. In Section III, we outline the Schwinger boson mean field theory. In Section IV, we perform the PSG analysis of quantum spin liquid phases with bosonic spinons and construct relevant mean field ansatz for the spin liquid states. In Section V, we analyze the energetics of the mean field spin liquid states via Schwinger boson mean field theory and the two-spinon continuum in each spin liquid state. We investigate the magnetically ordered states obtained by condensing bosonic spinons in the spin liquid phases. We then compare these results with the simulated annealing study of the corresponding classical model. In Section VI, we discuss relevance of our results to the existing and future experiments on Volborthite.

II Lattice Structure and Spin Model

In Volborthite $\mathrm{Cu_{3}V_{2}O_{7}(OH)_{2}\cdot 2H_{2}O}$, two layers of distorted Kagome lattice, each consisting of edge-sharing $\mathrm{CuO_{6}}$ octahedra, are separated by non-magnetic $\mathrm{V_{2}O_{7}}$ pillars and $\mathrm{H_{2}O}$ molecules.Yoshida et al. (2012) It is reasonable to assume that the interaction between different Kagome layers is negligible and we will therefore focus on just a single Kagome layer. The localized $S=1/2$ moment at each site on the Kagome lattice is carried by $\mathrm{Cu^{2+}}$ ions.Yoshida et al. (2012) Moreover, there are two crystallographically distinct $\mathrm{Cu^{2+}}$ sites, which suggests two different magnetically active orbitals.Janson et al. (2016) As shown below, the structure of the distorted Kagome layer is described by the plane crystallographic group $p2gg$, whose space group is discussed in Section II.1. This non-symmorphic version of the Kagome lattice (non-symmorphic Kagome lattice hereafter) possesses glide symmetry. We introduce the microscopic spin model of Volborthite, derived from a recent density functional theory (DFT) calculation Janson et al. (2016) in Section II.2.

II.1 Unit Cell and Space Group

The non-symmorphic Kagome lattice has six sites (or sublattices) per unit cell (FIG. 1). Denote the lattice constant along $x$- and $y$-direction by $b$ and $a$ respectively. Then, the coordinate of a generic site has the form $(xb,ya,s)$, which we simply write as $(x,y,s)$, where $x,y\in\mathbb{Z}$ specify the unit cell to which the site belongs, and $s=1,\ldots,6$ indexes the sublattice. The space group of non-symmorphic Kagome lattice is generated by $\pi$-rotation $C_{2}$ and glide $h$, which consists of reflection and half lattice translation. In general, a non-symmorphic operation combines a point group operation (e.g. rotation and reflection) with a fractional lattice translation, which cannot be rewritten in terms of point group operations and full lattice translations by switching to another coordinate system.Parameswaran et al. (2013) The non-symmorphic symmetry has important implication on the translational invariance of mean field ansatzes, which is discussed in Section IV.2.

For convenience of subsequent analysis, we also consider the lattice translations $T_{x}$ and $T_{y}$ along two independent directions $\hat{\mathbf{x}}$ and $\hat{\mathbf{y}}$. We fix the center of rotation at the center of hexagon in the $(0,0)$ unit cell, and the glide axis to the horizontal line passing through $(0,0,4)$ and $(0,0,5)$. $h$ is therefore the reflection about the glide axis followed by translation by $b\hat{\mathbf{x}}/2$. In Appendix A we show explicitly how a site with coordinates $(x,y,s)$ transforms under the space group operations.

We refer to the direction along sublattices $\ldots-1-2-1-2-\ldots$ (equivalently $\ldots-4-5-4-5-\ldots$ as they are related by symmetry) as the chain direction. The sites between adjacent chains, which have either sublattice index $3$ or $6$, are known as interstitial sites.

II.2 Microscopic Spin Model

According to the DFT calculation by Janson et al,Janson et al. (2016) the interaction between $S=1/2$ local moments on the Volborthite Kagome layer is described by the following Heisenberg model

[TABLE]

with four leading exchange interactions. They are the first and second nearest neighbour couplings along the chain direction, $J_{1}$ and $J_{2}$, and two different couplings between the chain and interstitial spins, $J$ and $J^{\prime}$, as shown in FIG. 1. The ratio between these coupling constants is given by $J:J^{\prime}:J_{1}:J_{2}=1:-0.2:-0.5:0.2$, with negative (positive) sign indicating ferromagnetic (antiferromagnetic) interaction.

III Schwinger Boson Mean Field Theory

In this section, we outline the Schwinger boson mean field approach Sachdev (1992); Wang and Vishwanath (2006); Messio et al. (2013) for the generic Heisenberg Hamiltonian (1). We start with Schwinger boson representation of the spin operator,

[TABLE]

where $b_{i\alpha}^{\dagger}$ ($b_{i\alpha}$) creates (annihilates) a bosonic spinon with spin $\alpha\in\{\uparrow,\downarrow\}$ at site $i$. Here $\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli matrices. The bosonic operators obey commutation relation,

[TABLE]

The total number of Schwinger bosons at site $i$ is represented by the number operator

[TABLE]

On the other hand, for localized moments of spin $S$, the total spin operator reads Messio et al. (2013)

[TABLE]

where the second equality follows from (2) and (4). This imposes a constraint on the number operator,

[TABLE]

Following Refs. Wang and Vishwanath, 2006; Sachdev, 1992, we allow $\kappa$ to be any positive real number, i.e. $S$ is not restricted to $1/2$. At the mean field level, the constraint is replaced by its ground state expectation value,Wang and Vishwanath (2006)

[TABLE]

Next, we define the bond operators Wang and Vishwanath (2006)

[TABLE]

with $\epsilon_{\alpha\beta}$ being the antisymmetric tensor. $\hat{A}_{ij}$ ($\hat{B}_{ij}$) is known as singlet pairing (hopping) channel. It can be checked explicitly that both $\hat{A}_{ij}$ and $\hat{B}_{ij}$ are invariant under global $SU(2)$ spin rotation. In terms of bond operators, the dot product of spin operators can be written as

[TABLE]

Heisenberg Hamiltonian (1) then becomes

[TABLE]

We have imposed the contraint (6) on all sites by introducing Lagrange multipliers (chemical potentials) $\mu_{i}$ in the Hamiltonian. Here $:\;:$ denotes normal ordering.

The Hamiltonian (10) is quartic in bosonic operators $b$. We perform mean field decoupling on $\hat{B}_{ij}^{\dagger}\hat{B}_{ij}$ and $\hat{A}_{ij}^{\dagger}\hat{A}_{ij}$ to obtain a Hamiltonian quadratic in $b$.Wang and Vishwanath (2006) Such a decoupling preserves the global $SU(2)$ spin rotation symmetry, which is present in the original Heisenberg Hamiltonian (1). Then, we replace $\langle\hat{A}_{ij}\rangle$ and $\langle\hat{B}_{ij}\rangle$ by complex-valued variational parameters $A_{ij}$ and $B_{ij}$ respectively. The mean field Hamiltonian reads

[TABLE]

Extremizing the expectation value of mean field Hamiltonian with respect to the variational parameters yields the following self-consistent equations,Wang and Vishwanath (2006)

[TABLE]

and the chemical potential yields the constraint (7),

[TABLE]

In the mean field decoupling scheme above, however, $A_{ij}$ ($B_{ij}$) channel would present an unbounded free energy for $J_{ij}<0$ ($J_{ij}>0$). In order to obtain a controlled mean field solution, we use the following identity Messio et al. (2013)

[TABLE]

which leads to two variants of (9),

[TABLE]

where we have assumed the constraint (6) and $i\neq j$. Therefore, for $J_{ij}>0$ ($J_{ij}<0$), we write the spin scalar product as (15b) ((15a)), which contains only $A_{ij}$ ($B_{ij}$). The multiplicative constant $2$ and the additive constant $\pm S^{2}$ can be further dropped without qualitatively changing the theory. Mean field decoupling then leads to the following Hamiltonian

[TABLE]

while the self-consistent equations are given by the same expression (12) as before. In addition, we assume that symmetry-related sites have the same chemical potential. Since there exist two chemically inequivalent sites, we denote the chemical potential at the chain sites ($s=1,2,4,5$) by $\mu_{1}$, and the interstitial sites ($s=3,6$) by $\mu_{2}$. The constraint (7) is then split into two,

[TABLE]

where $i$ now labels the unit cells rather than the individual sites. Notice that $\kappa$ is still the same at both sites since all spins have the same amplitude.

In practice, the self-consistent equations (12) are solved iteratively in momentum space. The procedure and techniques of mean field calculation are discussed in Appendix C.

IV Projective Symmetry Group Analysis

We first classify possible $\mathbb{Z}_{2}$ spin liquid phases with bosonic spinons for the two-dimensional non-symmorphic lattice described by $p2gg$ space group. Notice that this classification is independent of the spin Hamiltonian. If a specific spin model is chosen, one can investigate the spin liquid ground state that minimizes the mean field energy. We use the projective symmetry group (PSG) analysis,Wen (2002); Wang and Vishwanath (2006); Messio et al. (2013) namely we classify possible spin liquid mean field ansatzes that are invariant under the lattice symmetry transformations followed by a gauge transformation. Recall that the global spin rotation symmetry is automatically enforced by the choice of mean field parameters $A_{ij}$ and $B_{ij}$. Hence, in the following, we consider the space group and time reversal symmetry for the PSG analysis.

First, we observe that the mean field Hamiltonian (16) has $U(1)$ gauge redundancy. It is invariant under the following local $U(1)$ transformation Wang and Vishwanath (2006); Messio et al. (2013)

[TABLE]

Notice that the gauge transformed bosonic operators $\tilde{b}_{i\alpha}\equiv e^{i\phi(i)}b_{i\alpha}$ leave the representation (2) of spin invariant, satisfy the same commutation relation (3) and constraint (6), so that it describes the same physical spin. This implies that mean field ansatzes differing by a gauge transformation would correspond to the same physical state.

Now let us consider a space group element $X$ acting on the bosonic operator via

[TABLE]

Instead of considering the action on the bosonic operator, we can equivalently view $X$ as acting on the ansatz. We want $H_{\mathrm{MF}}$ to respect every symmetry operation $X$, such that the physics it describes is invariant under the transformation (19). Therefore, under the action of $X$, the modified ansatz is equivalent to the original one up to a gauge transformation, because they result in the same physical state. The set of all compound operations consisting of symmetry and gauge transformations that leave the ansatz invariant is defined as the projective symmetry group (PSG).Wen (2002) To illustrate this, consider for example the term $A_{ij}\hat{A}_{ij}^{\dagger}\sim A_{ij}b_{i\alpha}^{\dagger}b_{j\beta}^{\dagger}$ in the mean field Hamiltonian (16). Applying the symmetry transformation $X$ followed by the associated gauge transformation $G_{X}(\mathbf{r})=\exp\left[i\phi_{X}(\mathbf{r})\right]$ on such a term, we get Hwang et al. (2013)

[TABLE]

With the sites $i$ and $j$ being summed over, the invariance of the mean field ansatz and thus $H_{\mathrm{MF}}$ under $G_{X}X$ requires that

[TABLE]

or $A_{ij}=\exp\left[i\phi_{X}(X(i))+i\phi_{X}(X(j))\right]A_{X(i)X(j)}$ as (18b) and (19) might intuitively suggest. It can be similarly shown that

[TABLE]

Then, the collection of all $G_{X}X$ that leave the ansatz invariant is called the PSG. This definition also includes $G_{\mathcal{T}}\mathcal{T}$ for the antiunitary time reversal $\mathcal{T}$. We discuss how to treat time reversal symmetry explicitly in bosonic PSG in Appendix B.1.

Suppose that $G_{X}X\in\mathrm{PSG}$ for a space group element $X$. If we apply a gauge transformation on the mean field ansatz, say $A_{ij}\longrightarrow GA_{ij}$, then $G_{X}\longrightarrow GG_{X}XG^{-1}X^{-1}$ such that $G_{X}X$ is still an element of the PSG.Wang and Vishwanath (2006) The corresponding change in the phase of the bosonic operator is given by

[TABLE]

Due to the antiunitarity of $\mathcal{T}$, $G_{\mathcal{T}}$ transforms in a different manner (54), which is explained in Appendix B.1.

The PSG elements of the form $G_{I}I$ with $I$ being the identity element of the space group form a subgroup of PSG called the invariant gauge group (IGG). Messio et al. (2013) Alternatively, we can view IGG as the set of all pure gauge transformations that leave the ansatz invariant. For the mean field Hamiltonian (16) in which both $A_{ij}$ and $B_{ij}$ are present, the IGG is just $\mathbb{Z}_{2}=\{-1,1\}$. Wang and Vishwanath (2006)

IV.1 Algebraic PSG

The algebraic relations among space group elements constrain the possible forms of gauge transformations $G_{X}$. Wang and Vishwanath (2006) For instance, consider the string of translation operators that equals to identity,

[TABLE]

Suppose $G_{T_{x}}T_{x},G_{T_{y}}T_{y}\in\mathrm{PSG}$, then we must have

[TABLE]

or, in terms of phases,

[TABLE]

The PSG in which the gauge transformations $G_{X}$ satisfy algebraic constraints such as (24) is called algebraic PSG. For a given lattice, there is only a finite number of independent algebraic identities such as (23). These identities can be found by inspecting how two distinct space group elements commute, which are listed in Appendix A in the case of the non-symmorphic Kagome lattice. Representing the gauge transformations by their phases, the final solution of algebraic PSG is given by

[TABLE]

with

[TABLE]

The three independent $\mathbb{Z}_{2}$ variables $p_{2},p_{3},p_{13}\in\{0,1\}$ lead to $2^{3}=8$ distinct bosonic spin liquid states. Detailed derivation of the algebraic PSG can be found in Appendix B. We remark that, while $p_{2}$ and $p_{3}$ arise entirely from the space group considerations, $p_{13}$ is introduced only when time reversal symmetry/invariance is explicitly enforced. If only the spatial symmetry is enforced, the spin liquid phases are simply classified by $(p_{2},p_{3})$, which would in principle allow both time reversal invariant as well as time reversal breaking spin liquid states. In the mean field analysis of spin liquid phases in Section V, we will study possible spin liquid states by requiring only the spatial symmetry. However, we will find that all of the stable mean field solutions, which belong to the $(p_{2},p_{3})$ classification, satisfy time reversal symmetry and they are related to the time reversal invariant $(p_{1},p_{2},p_{13}=0)$ states.

IV.2 Mean Field Ansatz

Considering the spin model described in Section II.2 and the generic mean field Hamiltonian (16), we can see that there are altogether four independent pairing and hopping amplitudes per unit cell, which we denote by $A$, $B^{\prime}$, $B_{1}$ and $A_{2}$, depending on which exchange coupling they are associated with. All other amplitudes can be generated from these by lattice symmetry transformations. For example, let us fix $A=A_{(0,0,3)\longrightarrow(0,0,4)}$. The amplitude $A_{\left(1,0,3\right)\longrightarrow(0,0,2)}$ is related to $A_{(0,0,3)\longrightarrow(0,0,4)}$ by $C_{2}$ and the PSG ensures that the mean field ansatz is invariant under $G_{C_{2}}C_{2}$. Using (20),

[TABLE]

If time reversal symmetry is considered, the corresponding phase $\phi_{\mathcal{T}}$ further constrains the complex phase of $A$, $B^{\prime}$, $B_{1}$ and $A_{2}$ (see Appendix B.1).

Notice that the phase variables related to the lattice translations, as shown in (25a) and (25b), are trivial. To construct $H_{\mathrm{MF}}$, it is therefore sufficient to determine various relations between mean field amplitudes such as (26) in the $\left(0,0\right)$ unit cell, as all other unit cells have the same relations via (20) and (21) with $X=T_{x},T_{y}$. In other words, the mean field ansatz does not go beyond the physical unit cell. This is a consequence of the non-symmorphic symmetry of the lattice, which we explain as follows. The algebraic identities (36) and (38) impose the constraints ${\tilde{T}_{x}}^{-1}\tilde{h}^{2}=\pm 1\equiv\eta_{h}$ and $\tilde{h}^{-1}\tilde{T}_{y}\tilde{h}\tilde{T}_{y}=\pm 1\equiv\eta_{hT_{y}}$ on PSG, where we have used the abbreviation $\tilde{X}=G_{X}X$. Then,

[TABLE]

(27) dictates that the mean field ansatz constructed from PSG (see Section IV.2) can never enlarge the physical unit cell, because it would require ${\tilde{T}_{x}}^{-1}{\tilde{T}_{y}}^{-1}\tilde{T}_{x}\tilde{T}_{y}=-1$ if the contrary were true. This argument also holds for three other non-symmorphic plane crystallographic groups $pg$, $p2mg$ and $p4gm$. Such feature is characteristically different from the isotropic Kagome lattice, where, for example, the glide symmetry is absent. The algebraic PSG of the isotropic Kagome lattice allows certain spin liquid states such as the $\pi$-flux state (in the terminology of Ref. Wang and Vishwanath, 2006) that is given by a mean field ansatz that enlarges the physical unit cell.

V Emergent Quantum Phases and Their Physical Properties

As mentioned earlier, we consider the spin liquid phases labeled by $(p_{2},p_{3})$, which are constrained from the spatial symmetries. We determine the ground state of the spin model derived in the DFT computation, for different values of $2S=\kappa$. In Schwinger boson mean field theory, the condensation of spinons at certain wavevector arises as the number of bosons per site, $\kappa$, increases and then exceeds a critical value $\kappa_{c}$.Wang and Vishwanath (2006) This corresponds to the transition from a given spin liquid state to a magnetically ordered state. Notice that increasing the “spin” magnitude $S=\kappa/2$ reduces quantum fluctuations. The magnetically ordered phases that arise right after the transition can be obtained by analyzing the critical eigenmode near the transition. We identify such magnetically ordered phases obtained from different spin liquid phases. We also compare these results to the simulated annealing result of the classical model, which is equivalent to $\kappa\longrightarrow\infty$ limit in the Schwinger boson formulation.

In order to better characterize the spin liquid states, we compute the dispersion of the lower boundary of two-spinon continuum. As shown below, the spin liquid phases labeled by $p_{3}=1$ exhibit periodicity enhancement, namely the two-spinon continuum are invariant under translation by $\mathbf{k}=(\pm\pi,\pm\pi)$ in momentum space, which leads to spectral doubling in the Brillouin zone. We make the connection between this phenomenon and the notion of symmetry fractionalization Essin and Hermele (2013, 2014) introduced in previous works.

V.1 Spin Liquid States

Consider the spin liquid states labeled by $(p_{2},p_{3})$, where $p_{2},p_{3}\in\{0,1\}$. The critical $\kappa=\kappa_{c}$, where the spinons condense and the minimum of spinon dispersion $\omega_{\mathrm{min}}$ touches zero energy, is shown in TABLE 1. When $\kappa>\kappa_{c}$, a magnetic order arises. The value of $\kappa_{c}$ is obtained by computing $\omega_{\rm min}$ at various values of $\kappa$ and then making a linear extrapolation, as shown in FIG. 6 in Appendix D. The values of $\kappa_{c}$ fall between 0.3 and 0.4.

We compare the total energy $E=\langle H_{\mathrm{MF}}\rangle$ of four spin liquid states at several values of $\kappa<\kappa_{c}$, as shown in TABLE 2 in Appendix E. At $\kappa$ far below $\kappa_{c}$, the spin liquid states with the same $p_{2}$ have exactly the same energy, which happens because each of the hopping amplitudes $B_{1}$ and $B^{\prime}$ vanishes while each of the pairing amplitudes $A_{2}$ and $A$ has the same magnitude. As $\kappa$ approaches $\kappa_{c}$, the situation changes and all of these amplitudes become finite, which lifts the degeneracy. We find that $(p_{2},p_{3})=(1,0)$ is the most energetically favorable state among four spin liquid phases. Nevertheless, it should be noticed that the energy of $(1,1)$ spin liquid state is particularly close to $(1,0)$. We will pay particular attention to $(1,0)$ and $(1,1)$ spin liquid states as they are closely competing phases.

To gain further insight about these spin liquid phases, we compute the dispersion $\Omega_{\mathbf{k}}$ of the lower boundary of two-spinon continuum, which is given by

[TABLE]

where $\omega_{\mathbf{q}}$ is the one-spinon dispersion. We plot $\Omega_{\mathbf{k}}$ of each spin liquid state with $\kappa$ close to $\kappa_{c}$, for $-2\pi\leq k_{x},k_{y}\leq 2\pi$, as shown in FIG 2. Darker regions indicate lower energy.

Notice that $\Omega_{\mathbf{k}}$ with the same $p_{3}$ have similar profiles. For $p_{3}=0$, the minimum of the two-spinon continuum occurs at $\mathbf{k}=(\pm q,0)$ (as well as at the zone center), which suggests that $(q,0)$ magnetic order would develop beyond $\kappa>\kappa_{c}$. On the other hand, for $p_{3}=1$, the minimum occurs at $(\pi,\pi)$ (as well as at the zone center), which indicates that $(\pi,\pi)$ magnetic order would arise beyond $\kappa_{c}$. Details of the two-spinon continuum and the associated magnetic ordering patterns will be discussed in the next two sections.

Finally, it can be explicitly checked that four spin liquid solutions labelled by $(p_{2},p_{3})$ are time-reversal invariant by computing the fluxWang and Vishwanath (2006); Messio et al. (2013) piercing through the length-6 hexagon and three independent length-8 rhombi on the non-symmorphic Kagome lattice. We confirm that these fluxes are always [math] or $\pi$ in our solution, which is a necessary condition for the time reversal invariance. It can be checked that these solutions are indeed related to $(p_{2},p_{3},p_{13}=0)$ states in the full classification of the time reversal invariant spin liquid phases. Hence the flux counting is consistent with the PSG analysis.

V.2 Magnetically Ordered States

When the spin magnitude or the bosonic density $\kappa=2S$ reaches its critical value $\kappa_{c}$, the excitation spectrum becomes gapless and spinons condense at particular wavevectors $\mathbf{k}_{c}$, where $\omega_{\mathbf{k}_{c}}=0$, causing a phase transition from a spin liquid to a long range magnetically ordered state. The spinor form of the spinon operator $\Psi_{\mathbf{k}}=(b_{\mathbf{k}\uparrow},b_{-\mathbf{k}\downarrow}^{\dagger})^{\mathrm{T}}$ (see (68) in Appendix C), gains a finite expectation value at these $\mathbf{k}$ points, which is proportional to the respective critical eigenvector that becomes soft at $\mathbf{k}_{c}$. Analyzing these eigenmodes, we can determine the real space ordering patterns.Sachdev (1992); Wang and Vishwanath (2006) Below we focus on the two most energetically favorable spin liquid phases, namely $(p_{2},p_{3})=(1,0)$ and $(1,1)$. It is found that $\mathbf{k}_{c}=\pm(0.44\pi,0)$ for the state $(1,0)$, while $\mathbf{k}_{c}=\pm(\pi/2,\pi/2)$ for $(1,1)$.

More precisely, we first determine the expectation value of the real space spinor form $x_{i}\equiv(\langle b_{i\uparrow}\rangle,\langle b_{i\downarrow}\rangle)^{\mathrm{T}}$ via the Fourier transformation of $\langle\Psi_{\mathbf{k}}\rangle=(\langle b_{\mathbf{k}\uparrow}\rangle,\langle b_{-\mathbf{k}\downarrow}^{\dagger}\rangle)^{\mathrm{T}}$, which is dominated by the contribution at $\mathbf{k}=\mathbf{k}_{c}$. The real space spin configuration is then obtained from

[TABLE]

which we plot in FIG. 3a and 3b for states $(1,0)$ and $(1,1)$ respectively.

In the case of the (1,0) spin liquid state, a coplanar incommensurate $(q,0)$ spiral order with $q=0.88\pi$ develops as $\kappa$ reaches the critical value. Here the spins rotate by $0.88\pi$ under lattice translation $T_{x}$, while they do not change under $T_{y}$. The amplitude of spin is the same among the chain (interstitial) sites, but it is larger at the chain sites than at the interstitial sites. On the other hand, a coplanar commensurate $(\pi,\pi)$ spin density wave order develops if one starts from the $(1,1)$ spin liquid state. Here the spins rotate by $\pi$ under lattice translation $T_{x}$ or $T_{y}$. The amplitude of spin is the same among the chain (interstitial) sites, but it is larger at the chain sites than at the interstitial sites. Notice that in both configurations, the spins which interact by the dominant antiferromagnetic coupling $J$ are anti-aligned.

To elucidate the nature of these magnetically ordered phases, we investigate the ground state of the classical model, where we treat the spins in the Heisenberg model (1) as three-component vectors of fixed length, $\mathbf{S}_{i}=(S^{x}_{i},S^{y}_{i},S^{z}_{i})$ and $\lvert\mathbf{S}_{i}\rvert=1$. We use the same exchange interactions given in Section II.2. Simulated annealing is employed to obtain the ground state spin configuration on a lattice with $64\times 32$ unit cells, which turns out to be a coplanar $(q,0)$ spiral order with $q=3\pi/4$, as shown in FIG. 4.

Notice that the ground state spin configuration in the classical Heisenberg model agrees quite well with that arising from the most energetically favorable spin liquid state $(1,0)$, except the numerical value of $q$ in the ordering wave vector and the uniformity of spin amplitude. Hence it is natural to conclude that the $(q,0)$ spiral order derived from the (1,0) spin liquid state is continuously connected to the classical limit. It is, therefore, natural to call this state a $(q,0)$ spiral order. On the other hand, the spin density wave order arising from the $(1,1)$ spin liquid state has no classical analog. Hence the emergence of this state is a purely quantum effect as this state can only appear as the ground state if the amplitude of spins varies over different sites. This is the reason why we call this state a $(\pi,\pi)$ spin density wave. Since the local moments of Volborthite carry $S=1/2$, this material is not close to the classical limit, which suggests that the $(q,0)$ spiral and $(\pi,\pi)$ spin density wave states may be competing magnetic orders below $1\,\mathrm{K}$ in Volborthite.

V.3 Periodicity Enhancement of Two-Spinon Continuum

The dispersion of the lower boundary of two-spinon continuum is shown in FIG. 2 for different spin liquid phases. It can be seen that the spectra with $p_{3}=1$ exhibit periodicity enhancement, namely the translation by $\mathbf{k}=(\pm\pi,\pm\pi)$ in momentum space leaves the spectra invariant, which is beyond the periodicity allowed by the lattice translational symmetry. Such an enhanced periodicity, which leads to the spectral doubling,Essin and Hermele (2014) is a consequence of the spatial inversion $C_{2}$ and time reversal $\mathcal{T}$ symmetries. To prove this, we first observe that a physical spin operator, which is a bilinear form of two spinon operators, must transform trivially under the symmetry operations that amount to identity. However, each spinon can transform projectively under the same symmetry operations.Lu et al. It can gain a phase of $\pm\pi$ which is determined by the solution of algebraic PSG. In particular, the $\mathbb{Z}_{2}$ variable $p_{3}$ characterizes how a spinon transforms under the algebraic identities (34) and (35), i.e. $\tilde{C}_{2}\tilde{T}_{x}\tilde{C}_{2}^{-1}\tilde{T}_{x}=(-1)^{p_{3}}$ and $\tilde{C}_{2}\tilde{T}_{y}{\tilde{C}_{2}}^{-1}\tilde{T}_{y}=(-1)^{p_{3}}$ (see Appendix B), where $\tilde{X}\equiv G_{X}X$.

Let us denote a one-spinon momentum eigenstate by $\lvert\mathbf{q}\rangle$, which is necessarily an energy eigenstate with eigenvalue $\omega_{\mathbf{q}}$. Since momentum is the generator of translation, we get $\tilde{T}_{a}\lvert\mathbf{q}\rangle=\exp(iq_{a})\lvert\mathbf{q}\rangle$ where $a=x,y$. Consider the state $\lvert\mathbf{q}^{\prime}\rangle=\tilde{\mathcal{T}}\tilde{C}_{2}\lvert\mathbf{q}\rangle$, which is degenerate with $\lvert\mathbf{q}\rangle$ although in general $\mathbf{q}^{\prime}\neq\mathbf{q}$. To see how $\mathbf{q}^{\prime}$ is related to $\mathbf{q}$, we apply the one-spinon translation operator $\tilde{T}_{a}$,

[TABLE]

where we have used the fact that translation commutes with time reversal at one-spinon level, i.e. $\tilde{T}_{a}\tilde{\mathcal{T}}=\tilde{\mathcal{T}}\tilde{T}_{a}$, from our solution of the algebraic PSG (see Appendix B). (30) tells us that $\mathbf{q}^{\prime}=\mathbf{q}+p_{3}(\pi,\pi)$ while $\omega_{\mathbf{q}^{\prime}}=\omega_{\mathbf{q}}$. Since the two-spinon energy takes the form $\Omega_{\mathbf{p}+\mathbf{q}}=\omega_{\mathbf{p}}+\omega_{\mathbf{q}}$, we have

[TABLE]

When $p_{3}=0$, this equation is trivially satisfied. When $p_{3}=1$, we have $\Omega_{\mathbf{k}+(\pi,\pi)}=\Omega_{\mathbf{k}}$, which is the periodicity enhancement we observe in, for example, FIG. 5. We have thus shown that the enhanced periodicity arises from the symmetry that combines spatial inversion $C_{2}$ and time reversal $\mathcal{T}$. This class of enhanced periodicity has been discussed in Ref. Essin and Hermele, 2014 and explained through the language of symmetry fractionalization,Essin and Hermele (2013) where the symmetry action on a composite physical operator/state (e.g. spin) can be represented by the product of individual symmetry action on each constituent spinon operator/state. Hence our result can also be interpreted as a consequence of symmetry fractionalization.

VI Discussion

In this work, we investigate possible quantum spin liquid and magnetically ordered phases in a two-dimensional non-symmorphic lattice, motivated by the experiments on Volborthite, $\mathrm{Cu_{3}V_{2}O_{7}(OH)_{2}\cdot 2H_{2}O}$, and the DFT computations. In Volborthite, a magnetic ordering Yoshida et al. (2009b); Ishikawa et al. (2015); Yoshida et al. occurs below $1\,\mathrm{K}$ while a recent thermal conductivity measurement Watanabe et al. found some signatures of possible spin liquid behavior right above $1\,\mathrm{K}$. This may suggest that the paramagnetic state above $1\,\mathrm{K}$ may be proximate to a putative quantum spin liquid that is intimately related to the magnetic order below $1\,\mathrm{K}$. The nature of the magnetic order below $1\,\mathrm{K}$, however, is not fully understood. Earlier nuclear magnetic resonance (NMR) measurements suggest that the magnetic order could be a helical order or a spin density wave.

To address these issues, we use the Schwinger boson mean field theory Sachdev (1992); Wang and Vishwanath (2006); Messio et al. (2013) and PSG Wen (2002) to understand the connection between possible magnetic order and quantum spin liquid phases on equal footing, in the spin model Janson et al. (2016) derived from the DFT computation. The DFT-derived model suggests that the underlying lattice has non-symmorphic symmetry, described by $p2gg$ planar space group. We first show that mean field spin liquid states in two-dimensional non-symmorphic lattices do not enlarge the lattice unit cell because of the glide symmetry, the combination of refection and a fractional translation. We analyze the resulting mean field ground states and find that there are two competing quantum spin liquid states that lead to a $(q,0)$ spiral magnetic order and a $(\pi,\pi)$ spin density wave, respectively, upon the spinon condensation. This suggests that the $(q,0)$ spiral order and $(\pi,\pi)$ spin density wave would be competing magnetically ordered state below $1\,\mathrm{K}$, and two spin liquid phases mentioned above are candidate paramagnetic states above $1\,\mathrm{K}$. Notice that both of the $(q,0)$ spiral order and $(\pi,\pi)$ spin density wave would be consistent with earlier NMR experiments. It is also found that the spin liquid phase related to the $(\pi,\pi)$ spin density wave shows periodicity enhancement of the two-spinon continuum, a signature of the so-called symmetry fractionalization.Essin and Hermele (2013, 2014) Hence a future neutron scattering experiment may be able to find such signatures in case that this spin liquid state is relevant to the paramagnetic state above $1\,\mathrm{K}$.

Finally, the spin liquid states with bosonic spinons in two dimensions are necessarily gapped phases as a magnetic ordering will occur if the spin-carrying bosons become gapless. While a small excitation gap at finite temperature in the spin liquid phases with bosonic spinons may explain the large thermal conductivity discovered in Volborthite, it will also be useful to consider fermionic versions of these spin liquid phases, which naturally allows gapless spin liquid states and a large thermal conductivity. This would be an excellent topic of future study.

Acknowledgements.

We thank Sopheak Sorn for initial collaboration and Robert Schaffer for useful discussion. This work was supported by the NSERC of Canada and the Center for Quantum Materials at the University of Toronto. T. M. was supported by Advanced Leading Graduate Course for Photon Science (ALPS). Some of the computations were performed on the GPC supercomputer at the SciNet HPC Consortium.Loken et al. (2010) SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.

Appendix A Space Group and Algebraic Identities of Non-Symmorphic Kagome Lattice

Below we list the action of space group elements, $T_{x}$, $T_{y}$, $C_{2}$ and $h$, on coordinate $(x,y,s)$.

[TABLE]

The algebraic relations among $T_{x}$, $T_{y}$, $C_{2}$ and $h$ are

[TABLE]

If time reversal symmetry is considered, we further have

[TABLE]

(32) $-$ (41) constrain the possible form of gauge transformation $G_{X}$ associated with symmetry operation $X=T_{x},T_{y},C_{2},h,\mathcal{T}$ such that $G_{X}X\in\mathrm{PSG}$. Note that (36) and (37) are not independent of each other, but we are going to use both forms to find the solution to algebraic PSG in Appendix B.

Appendix B Solution to Bosonic PSG

Here we determine the form of gauge transformations $G_{X},X=T_{x},T_{y},C_{2},h,\mathcal{T}$ that are consistent with the algebraic constraints imposed by the identities (32) $-$ (41). Relevant discussion can be found in Section IV.1. We represent the gauge transformation by its phase,

[TABLE]

First and foremost, we follow the usual convention Wang and Vishwanath (2006) of fixing

[TABLE]

The algebraic identity (32) leads to

[TABLE]

where $p_{n}=0$ or $1\,\mathrm{mod}\,2$. Simplifiying yields

[TABLE]

Using gauge freedom, we can fix $\phi_{T_{y}}(0,y,s)=0$. Therefore,

[TABLE]

Before we proceed to the calculation of rotation and glide, we can eliminate additional parameters in (36) and (39) by a global $\mathbb{Z}_{2}$ gauge. Schematically, (39) gives

[TABLE]

where $\ldots$ denotes some undetermined phases of $C_{2}$ and $h$. The trick is to add a global phase of $p_{8}\pi$ to $\phi_{T_{y}}(x,y,s)$, which does no harm because the IGG is $\mathbb{Z}_{2}$. Now $p_{8}\pi$ appears on both sides of (43), so they cancel out each other. Since all other identities involve an even number of $T_{y}$, without loss of generality we can set $p_{8}=0\,\mathrm{mod}\,2$. Schematically, (36) gives

[TABLE]

where $\ldots$ denotes some undetermined phases of $h$. Similar to previous case, we add a global phase of $p_{5}\pi$ to $\phi_{T_{x}}(x,y,s)$. (44) then forces $p_{5}=0\,\mathrm{mod}\,2$.

From (34) and (35), we have

[TABLE]

From (33), we have

[TABLE]

Consider the sublattice dependent gauge transformation defined by

[TABLE]

$\phi_{T_{x}}$ and $\phi_{T_{y}}$ are invariant under $G_{1}$, while $\phi_{C_{2}}$ is modified by

[TABLE]

Therefore, by an appropriate choice of $\phi_{1}$, we can fix

[TABLE]

Furthermore, with the gauge transformations

[TABLE]

we can fix

[TABLE]

From (37) and (38), we have

[TABLE]

From (36), starting from $(x,y,1)$, we have

[TABLE]

Since R.H.S. has no coordinate dependence, this forces $p_{1}=0\,\mathrm{mod}\,2$. From (36), starting from $(x,y,4)$,

[TABLE]

From (36), starting from $(x,y,5)$, we have

[TABLE]

The two equations above implies that $p_{6}=0\,\mathrm{mod}\,2$. From (39), starting from $(x,y,3)$, we have

[TABLE]

From (39), starting from $(x,y,6)$, we have

[TABLE]

The two condition above implies that $p_{7}=p_{4}\,\mathrm{mod}\,2$. Multiplying the last equation by 2 gives

[TABLE]

Therefore $p_{4}=p_{3}\,\mathrm{mod}\,2$. To eliminate $p_{4}$ in favor of $p_{3}$ in $\phi_{C_{2}}(0,0,6)$, we have to treat division by 2 carefully since there may be an additional phase of $\pi$,

[TABLE]

where $q=0$ or $1\,\mathrm{mod}\,2$. $q$ can be expressed in terms of $p_{2}$ and $p_{3}$ by

[TABLE]

or $q=\left(p_{2}+p_{3}\right)\,\mathrm{mod}\,2$. Therefore,

[TABLE]

From (39), starting from $(x,y,1)$,

[TABLE]

From (39), starting from $(x,y,2)$,

[TABLE]

Evaluating (39) starting from $(x,y,s=4,5)$ generates the same two equations above.

From (36), starting from $(x,y,3)$,

[TABLE]

Now, we should determine the form of $\phi_{h}(0,0,s)$. Using the sublattice dependent gauge transformation defined by

[TABLE]

we first fix

[TABLE]

We are then left with the system of linear equations from previous calculations

[TABLE]

Using a similar sublattice dependent gauge transformation, we can fix

[TABLE]

which then implies

[TABLE]

The final solution to algebraic PSG involving only spatial symmetries is given by

[TABLE]

with

[TABLE]

which is described by only two $\mathbb{Z}_{2}$ variables $p_{2}$ and $p_{3}$.

B.1 Time Reversal Symmetry in Bosonic PSG

Conventionally, time reversal symmetry is enforced by flux argument in Schwinger boson mean field theory, unlike the fermionic approach Wen (2002); Lu et al. (2011); Huang et al. which involves time reversal symmetry directly in PSG calculation together with the space group. Since Lu et al proposed a scheme of unifying bosonic and fermionic theories of spin liquid through vison PSG,Lu et al. there have been several attempts Chatterjee et al. ; Schaffer et al. ; Yang and Wang (2016) to treat time reversal symmetry in bosonic PSG on equal footing with fermionic PSG. Below we clarify such an approach and solve the bosonic PSG from algebraic identities (40) and (41) involving time reversal $\mathcal{T}$, so that we can establish connection between bosonic and fermionic spin liquid states in the future.

In analogy to fermionic PSG, we can define the action of $G_{X}X\in\mathrm{PSG}$ on the pairing ansatz $A_{ij}$ for a space group element $X$ as

[TABLE]

while the hopping ansatz $B_{ij}$ follows a similar transformation rule. Note that (49) is consistent with (20) and the derivation of algebraic PSG above. As argued in Ref. Wen, 2002, time reversal operator complex conjugates the ansatz

[TABLE]

In fermionic PSG, when the ansatz is complex conjugated by time reversal, we have the freedom to perform $SU(2)$ gauge transformation $i\tau^{2}$ to trade the complex conjugation for a minus sign.Wen (2002) However, such an $SU(2)$ gauge redundancy is not present in the bosonic Hamiltonian (16), so the complex conjugation (50) cannot be removed. If the system has time reversal symmetry, then $A_{ij}$ and $A_{ij}^{*}$ describe the same physical state, and thus by the argument in Section IV they must be equal up to a $U(1)$ gauge transformation, say

[TABLE]

Recall that PSG is defined as the group of compound operators $G_{X}X$ that leave the mean field ansatz invariant. For time reversal, the PSG element $G_{\mathcal{T}}\mathcal{T}$ is constructed such that, acting $\mathcal{T}$ first complex conjugates the ansatz, acting $G_{\mathcal{T}}$ next brings the complex conjugated ansatz back to the original one. Using (51),

[TABLE]

or simply (c.f. (49))

[TABLE]

The action of $\left(G_{\mathcal{T}}\mathcal{T}\right)^{-1}$ is naturally first multiply the inverse phase factor to the ansatz and then complex conjugate everything. The net effect of $\left(G_{\mathcal{T}}\mathcal{T}\right)^{-1}$ is thus same as $G_{\mathcal{T}}\mathcal{T}$,

[TABLE]

Suppose that $G_{X}X\in\mathrm{PSG}$ for a space group element $X$. Applying a gauge transformation on the mean field ansatz, $A_{ij}\longrightarrow GA_{ij}$, $G_{X}$ changes as (22) such that $G_{X}X$ is still a PSG element. However, since the time reversal operator complex conjugates the ansatz, it is easy to see that for the same gauge transformation above, $G_{\mathcal{T}}$ changes by

[TABLE]

such that $G_{\mathcal{T}}\mathcal{T}$ is still a PSG element.

The time reversal satisfies the algebraic identities (40) and (41), which constrain the form of $G_{\mathcal{T}}$ by

[TABLE]

and

[TABLE]

respectively, where $X=T_{x},T_{y},C_{2},h$ and $p_{X\mathcal{T}},p_{\mathcal{T}}\in\{0,1\}$. Applying (55) to the mean field ansatz step by step,

[TABLE]

To obtain the next line, we have used the fact that the previous line transforms as $A_{ij}$ by the definition of PSG, regardless of how complicated the expression is. This leads to the consistent condition

[TABLE]

On the other hand, applying (56) to the mean field ansatz yields

[TABLE]

where we have manually inserted a factor of $-1$ on both sides of $A_{ij}$ in the last line, because $\mathcal{T}^{2}$ acting on a $S=1/2$ object should produce a minus sign. This leads to

[TABLE]

or $p_{\mathcal{T}}=1$.

With the gauge $G_{X}$ associated with $X=T_{x},T_{y},C_{2},h$ fixed as (45) $-$ (48), we now use (58) to solve for $G_{\mathcal{T}}$. $X=T_{x},T_{y}$ gives

[TABLE]

$X=C_{2}$ gives

[TABLE]

which implies $p_{11}=p_{10}\,\mathrm{mod}\,2$ and $p_{12}=\left(p_{2}+p_{3}+p_{10}\right)\,\mathrm{mod}\,2$. Finally, $X=h$ gives

[TABLE]

which implies $p_{10}=0\,\mathrm{mod}\,2$. We are left with five equations to determine the form of $\phi_{\mathcal{T}}\left(0,0,s\right)$,

[TABLE]

Consider the sublattice dependent gauge transformation defined by

[TABLE]

which does not affect the phase $\phi_{X}$ associated with space group elements $X$. It changes $\phi_{\mathcal{T}}$ by

[TABLE]

By an appropriate choice of $\phi_{5}$, we can fix $\phi_{\mathcal{T}}\left(0,0,1\right)=0$ and $\phi_{\mathcal{T}}\left(0,0,5\right)=\left(p_{2}+p_{3}\right)\pi$ from (61). (63) and (65) then implies $\phi_{\mathcal{T}}\left(0,0,2\right)=p_{13}\pi$ and $\phi_{\mathcal{T}}\left(0,0,4\right)=\left(p_{2}+p_{3}+p_{13}\right)\pi$. Furthermore, with a similar sublattice dependent gauge transformation defined by

[TABLE]

we can fix $\phi_{\mathcal{T}}\left(0,0,3\right)=0$ and $\phi_{\mathcal{T}}\left(0,0,6\right)=p_{13}\pi$ from (64). In conclusion,

[TABLE]

with

[TABLE]

Solving the bosonic PSG involving time reversal symmetry in such a manner introduce an independent $\mathbb{Z}_{2}$ variable $p_{13}$ on top of $p_{2}$ and $p_{3}$, which arise from consideration of spatial symmetries only. There are in total $2^{3}=8$ possible bosonic spin liquid states labeled by $\left(p_{2},p_{3},p_{13}\right)$ that respect the space group of non-symmorphic Kagome lattice and time reversal symmetry. From (51), we see that $2\,\mathrm{Arg}\{A_{ij}\}=-\phi_{\mathcal{T}}(i)-\phi_{\mathcal{T}}(j)$, therefore time reversal symmetry restricts an ansatz to be either real or imaginary by (66).

Appendix C Mean Field Calculation Techniques

The flow of a mean field iteration is outlined as follows. First, we fix the bosonic density $\kappa$. Define time step $t$ such that it is initially [math] and increases by $1$ after an iteration is completed. At time $t=0$, we choose some random value for the mean field parameters $A_{ij}$ and $B_{ij}$. At time $t>0$, we solve for $\mu_{1}$ and $\mu_{2}$ from the constraints (17a) and (17b) with the value of $A_{ij}$ and $B_{ij}$ from time $t-1$. Call the solutions $\mu_{1}(t)$ and $\mu_{2}(t)$. Then, we evaluate the expectation values $\langle\hat{A}_{ij}\rangle$ and $\langle\hat{B}_{ij}\rangle$ with the set of inputs $\{\mu_{1}(t),\mu_{2}(t),A_{ij}(t-1),B_{ij}(t-1)\}$, which outputs $A_{ij}(t)$ and $B_{ij}(t)$ respectively. The mean field iteration is now completed and $t$ increases by $1$. If the value of mean field parameters converge upon a sufficiently large number of iterations, then a mean field solution is obtained.

In practice, mean field theory is solved by iteration of self-consistent equations (12) in momentum space. We perform Fourier transformation

[TABLE]

to obtain the mean field Hamiltonian (16), which is written in terms of the coupling constants $J,J^{\prime},J_{1},J_{2}$ and the mean field parameters $A,B^{\prime},B_{1},A_{2}$, in momentum space. In (67), $i$ labels the individual unit cells, $s=1,\ldots,6$ indexes the sublattice and $N$ is the total number of unit cells (the total number of sites is therefore $N_{s}=6N$). Define the $2n$-component vector spinon field Wang and Vishwanath (2006)

[TABLE]

where we have used the abbreviation $\vec{b}_{\mathbf{k}\uparrow}=(b_{\mathbf{k},\uparrow,1},\ldots,b_{\mathbf{k},\uparrow,n})^{\mathrm{T}}$ and $\vec{b}_{-\mathbf{k}\downarrow}^{\dagger}=(b_{-\mathbf{k},\downarrow,1}^{\dagger},\ldots,b_{-\mathbf{k},\downarrow,n}^{\dagger})^{\mathrm{T}}$, with $n=6$ being the total number of sublattices. The mean field Hamiltonian can be expressed as

[TABLE]

$\mathrm{D}_{\mathbf{k}}$ is a $2n\times 2n$ matrix that depends on the coupling constants, mean field parameters and chemical potentials. Evaluating the expectation value of bond and number operator, as in (12) and (13), requires diagonalizing the Fourier transformed mean field Hamiltonian by Bogoliubov transformation, which is discussed in Appendix C.1.

If the parameter space is small, we can just assign arbitrary initial values to the mean field parameters and iterate the self-consistent equations until they converge. Different set of initial values should be tried to avoid identifying solution that leads to local but not global minimum of the energy. Unfortunately, since we have $8$ independent mean field parameters in this particular problem, our parameter space is considerably large (a generic complex-valued parameter counts twice because the real and imaginary part are independent). Conventional mean field iteration is no longer an ideal primary tool to tackle such problem, as it becomes difficult to check whether a solution corresponds to the global minimum, and slow convergence is expected. We need a more efficient way to explore the solution space. This is done through simulated annealing (SA), which is discussed in Appendix C.2.

SA is basically a method of probabilistic search, which suggest an approximate solution to the global minimum. We can however refine the solution from SA by inputing it as initial condition for mean field iteration. Also, for bosonic densities $\kappa$ close to the critical value $\kappa_{c}$, tiny changes can easily yield numerically insensible result, which our algorithm rejects by construct, rendering the exploration of solution space difficult. Therefore, mean field iteration is still important. Our general strategy is to employ SA for several $\kappa$ sufficiently lower than $\kappa_{c}$, and then perform mean field iteration with solutions from SA as initial condition. We tune up $\kappa$ gradually and iterate the self-consistent equations, with the solution from previous mean field iteration as initial condition, until we reach $\kappa_{c}$.

C.1 Bogoliubov Transformation

Bogoliubov transformation allows us to diagonalize the Fourier transform matrix $\mathrm{D}_{\mathbf{k}}$ by introducing a new pair of annihilation and creation operators, $\gamma_{\mathbf{k}}$ and $\gamma_{\mathbf{k}}^{\dagger}$, which satisfy the original bosonic commutation relation (3). The $b$ operators are linear combinations of them,Wang and Vishwanath (2006); Sachdev (1992)

[TABLE]

where $\Gamma_{\mathbf{k}}$ is the vector spinon field (68) of $\gamma$ operators. The $2n\times 2n$ matrix $\mathrm{M}_{\mathbf{k}}$ is refered to as Bogoliubov transformation matrix. Define the $2n\times 2n$ diagonal matrix in which the first $n$ nonzero entries equal to $1$ and the remaining equal to $-1$,

[TABLE]

From the commutation relation (3) among the $b$ operators, we have the following identity Sachdev (1992)

[TABLE]

(72) is also true for $\Gamma_{\mathbf{k}}$ since Bogoliubov transformation has to preserve the commutation relation. Then,

[TABLE]

or simply Sachdev (1992)

[TABLE]

Diagonalization of the Fourier transform matrix $\mathrm{D}_{\mathbf{k}}$ requires

[TABLE]

where $\mathcal{E}_{\mathbf{k}}=\mathrm{diag}\,(\omega_{\mathbf{k},\uparrow,1},\ldots,\omega_{\mathbf{k},\uparrow,n},\omega_{-\mathbf{k},\downarrow,1},\ldots,\omega_{-\mathbf{k},\downarrow,n})$ is the eigenvalue matrix. This implies Sachdev (1992)

[TABLE]

Combining (73) and (74) gives Sachdev (1992)

[TABLE]

Therefore, we are really diagonalizing $\tau_{3}\mathrm{D}_{\mathbf{k}}$ instead of merely $\mathrm{D}_{\mathbf{k}}$ in the usual linear algebra sense, and the eigenvalue matrix is $\tau_{3}\mathcal{E}_{\mathbf{k}}$. The expectation value of $H_{\mathrm{MF}}$ and other operators (e.g. bond operators in (12)) is evaluated with respect to the ground state of $\gamma$ spinon defined by $\gamma_{\mathbf{k},\alpha,s}\lvert 0\rangle=0$ for every $\mathbf{k}$, $\alpha$ and $s$. If we express these operators in terms of $\gamma$ operators, their expectation values are given by summing certain matrix elements of $\mathrm{M}_{\mathbf{k}}$ over $\mathbf{k}$.

The diagonalization (75) is done numerically for complicated systems like Kagome lattice in which $\mathrm{M}_{\mathbf{k}}$ does not admit a nice analytical expression. We recommend Refs. Messio et al., 2013; Copla, 1978 to readers who are interested in the details of computer algorithm of Bogoliubov transformation.

C.2 Simulated Annealing

Since we have a large parameter space spanned by four independent complex-valued mean field amplitudes $A,B^{\prime},B_{1},A_{2}$, we implement simulated annealing (SA) to determine the solution that minimize the mean field energy $\langle H_{\mathrm{MF}}\rangle$. We follow the basic procedure outlined in Ref. Corana et al., 1987.

We would like to minimize a multivariable cost function $f\left(x_{1},\ldots,x_{n}\right)$ that is bounded below. First, we choose an arbitrary initial configuration $(x_{1}^{0},\ldots,x_{n}^{0})\in\mathbb{R}^{n}$ and evaluate $f_{0}\equiv f(x_{1}^{0},\ldots,x_{n}^{0})$. Then, we explore the parameter space by perturbing the variables probabilistically, for instance $x_{1}^{0}\longrightarrow x_{1}^{0}+\delta$ with $\delta$ being a random perturbation. Let $\epsilon$ be a positive constant that characterizes the magnitude of $\delta$, and $p\in[0,1]$ be a random number. We require the random perturbation $\delta$ to be bounded by $\epsilon$, which can be done by choosing $\delta=(-1+2p)\epsilon$. We evaluate $f_{1}\equiv f(x_{1}^{0}+\delta,x_{2}^{0},\ldots,x_{n}^{0})$ and compare $f_{1}$ with $f_{0}$. If $f_{1}\leq f_{0}$, we accept the change completely. If $f_{1}>f_{0}$, we accept the change with probability $\exp\left(-\left(f_{1}-f_{0}\right)/T\right)$, where $T>0$ is a parameter analogous to temperature in Boltzmann factor. This acceptance rule is known as Metropolis criterion. After that, we move on to perturb the remaining variables $x_{2}^{0},\ldots,x_{n}^{0}$, one at a time, in similar fashion. If every variable is perturbed once, we say that a cycle is completed. Carrying out a sufficiently large number of cycles while gradually decreasing the temperature, we can (hopefully) find a configuration $\left(y_{1},\ldots,y_{n}\right)$ that well estimates the global minimum.

The temperature $T$ in the acceptance probability $\exp\left(-\left(f-f_{0}\right)/T\right)$ plays a crucial role in the exploration of parameter space and the convergence of solution. For large $T$ compared to the typical change in cost function, $\Delta f=\overline{f-f_{0}}$, fluctuations in the variables $x_{i}$ occur more often such that a larger portion of the solution space is explored. For small $T$, solutions that yield higher cost tend to be rejected. We usually start with large $T$ and decrease it in subsequent cycles.

In our system, the cost function is just the mean field energy,

[TABLE]

with $\kappa$ fixed at some value. The chemical potentials $\mu_{1}$ and $\mu_{2}$ are not freely varying parameters because they are solved self-consistently by (17a) and (17b) for a given set of mean field amplitudes. Since each of the mean field parameters is complex-valued, we have to perturb the real and imaginary parts independently, and the solution space is effectively $\mathbb{R}^{8}$.

Convergent solution is found upon implementation of SA. As a check, we input the solution from SA as initial condition for mean field iteration, and verify that it is also the convergent solution of self-consistent equations (12).

Appendix D Linear Approximation Scheme

FIG. 6 shows the linear approximation scheme we use to determine the critical bosonic density $\kappa_{c}$ from mean field theory. Related discussion can be found in Section V.1.

Appendix E Energetics of Spin Liquid States

TABLE 2 lists the energy of four spin liquid states $(p_{2},p_{3})$ that respect the space group of non-symmorphic Kagome lattice at several bosonic densities $\kappa$ below the critical value $\kappa_{c}$. Related discussion can be found in Section V.1.