Topological Nematic Spin Liquid on the Square-Kagome Lattice

TL;DR

This paper investigates the ground states of the square-kagome lattice, revealing a topological nematic quantum spin liquid phase with unique properties, supported by mean-field theory and dynamical structure factor calculations.

Contribution

It introduces a new topological nematic quantum spin liquid phase in the square-kagome lattice and characterizes its stability and experimental signatures.

Findings

Discovery of a topological nematic quantum spin liquid phase.

Identification of incommensurate magnetic orders.

Provision of dynamical structure factors for experimental detection.

Abstract

The ground state of the spin$-1/2$ kagome antiferromagnet remains uncertain despite decades of active research. Here we step aside from this debated question to address the ground-state nature of a related, and potentially just as rich, system made of corner-sharing triangles: the square-kagome lattice (SKL). Our work is motivated by the recent synthesis of a distorted SKL compound mentioned in [Morita & Tohyama, J. Phys. Soc. Japan 87, 043704 (2018)]. We have studied its spin$-1/2$ $J_{1}$-$J_{2}$ phase diagram with an unrestricted Schwinger boson mean-field theory (SBMFT). We show that, in addition of agreeing with previous studies, three original phases appear: two incommensurate orders and a topological quantum spin liquid with weak nematicity. The topological order is characterized by fluxes on specific gauge-invariant quantities and the phase is stable under anisotropic perturbations relevant for experiments. Finally, we provide dynamical structure factors of the reported phases that could be observed in inelastic neutron scattering.

Equations

$$\mathcal{H}=\sum _{⟨i,j⟩}{J}_{ij}{\hat{\mathbf{S}}}_i\cdot {\hat{\mathbf{S}}}_j,$$

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

$$\hat{\mathcal{H}}$$

Full text

I Introduction

The search of novel topological phases is one of the most active field in Condensed Matter, often accompanied by the emergence of fractionalized excitations Laughlin (1983); Wen (1991); Moessner and Sondhi (2001); Kitaev (2006); Castelnovo et al. (2012). The characteristic absence of a local order parameter has made frustrated magnetism a successful playground for topological phases, since frustration inherently hinders order. Everyone knows the textbook example of a simple antiferromagnetic triangle. But despite decades of research, this minimal brick of frustration has not unveiled all its secrets. The ground state (GS) of the quantum kagome antiferromagnet, made of corner-sharing triangles, is still under debate between a gapless Dirac spin liquid Hermele et al. (2008); Iqbal et al. (2011); Läuchli et al. (2016); Liao et al. (2017); He et al. (2017); Ralko et al. (2018), a gapped $\mathbb{Z}_{2}$ spin liquid Yan et al. (2011); Jiang et al. (2012); Depenbrock et al. (2012), or something else Läuchli et al. (2016); Ralko et al. (2018). In this context, it is appealing to step away from the cumbersome kagome problem, and to consider an alternate, and potentially just as rich, network of this elementary brick of frustration.

The square-kagome lattice (SKL) offers such a possibility Siddharthan and Georges (2001) [Fig. 1]. While the kagome lattice supports only smallest loops of 6 sites, the SKL is paved with two types of small loops of length 4 and 8. On one hand, this asymmetry suggests a more localized quantum entanglement, and thus a possibly more tractable treatment. On the other hand, the SKL is also a rare and promising example where quantum resonance is not governed by the shortest loops Ralko and Rousochatzakis (2015); Ralko et al. (2018). Adding to that the recent announcement of a spin$-1/2$ material with distorted square-kagome geometry Morita and Tohyama (2018); Hasegawa et al. (2018); Fujihala et al. (2017), and the theoretical proposal to realise the SKL in optical lattices Glaetzle et al. (2014), it has become a timely question to understand the quantum GS of this model, moving beyond previous works on exact diagonalization (ED) Richter et al. (2009); Nakano and Sakai (2013); Rousochatzakis et al. (2013); Morita and Tohyama (2018) and perturbative methods Siddharthan and Georges (2001); Tomczak and Richter (2003); Ralko and Rousochatzakis (2015).

In this paper, we focus on the spin-$1/2$ Heisenberg model on the SKL. The asymmetry of the lattice suggests to treat two nonequivalent bonds $J_{1}$ and $J_{2}$, as displayed in Fig. 1. A further anisotropy $J_{2}^{\prime}$, relevant to experiments, will be considered at the end of the paper. The Hamiltonian then reads

[TABLE]

and its zero-temperature phase diagram [Fig. 2] is studied by means of an unrestricted Schwinger-boson mean-field theory (SBMFT). SBMFT is able not only to treat on an equal footing quantum spin liquids and magnetically ordered phases Auerbach (1994) but also to reach very large sizes, a rare feature for a quantum approach. Defining the parameter $x=J_{2}/J_{1}>0$, we show that for small and large values, our SBMFT is consistent with the plaquette and ferrimagnetic phases reported in Ref. Rousochatzakis et al., 2013. Then, as $x$ closes to 1 from both sides, we report the appearance of two gapless incommensurate orders, I1 and I2. Eventually when $x\sim 1$, they give rise to a topological nematic spin liquid (TNSL). Finally, to make contact with materials, we discuss the robustness of the TNSL under further anisotropy Morita and Tohyama (2018) and provide inelastic neutron scattering signatures of the reported phases.

II Schwinger boson mean-field theory

SBMFT has proven successful in a variety of models: to analytically study antiferromagnetic GS Sachdev (1992), to classify spin liquids using projective symmetry groups Wang and Vishwanath (2006); Wang (2010); Messio et al. (2013) and even to establish phase diagrams Cabra et al. (2011); Zhang and Lamas (2013); Merino and Ralko (2018). The method consists in introducing two Schwinger bosons, $\hat{b}^{\dagger}_{i\uparrow}$ and $\hat{b}^{\dagger}_{i\downarrow}$, on each site $i$ of the lattice through the following spin mapping

[TABLE]

where $\boldsymbol{\sigma}$ are the Pauli matrices and $\alpha,\beta=\ \uparrow,\downarrow$. For this mapping to be correct, an on-site boson constraint $\hat{b}^{\dagger}_{i\uparrow}\hat{b}_{i\uparrow}+\hat{b}^{\dagger}_{i\downarrow}\hat{b}_{i\downarrow}=2S$ has to be fulfilled. This constraint being very hard to implement exactly, it is thus imposed on average by introducing Lagrange multipliers $\{\mu_{i}\}$. Two SU(2)-invariant bond operators Flint and Coleman (2009), the singlet operator on the oriented bond $(i\rightarrow j)$ $\hat{A}_{ij}=(\hat{b}_{i\uparrow}\hat{b}_{j\downarrow}-\hat{b}_{i\downarrow}\hat{b}_{j\uparrow})/2$ and the boson hopping operator $\hat{B}_{ij}=(\hat{b}^{\dagger}_{i\uparrow}\hat{b}_{j\uparrow}+\hat{b}^{\dagger}_{i\downarrow}\hat{b}_{j\downarrow})/2$, are introduced. They are respectively favored in disordered and ordered phases. Following a mean-field decoupling of these operators, Eq. 1 becomes:

[TABLE]

where $A_{ij}=\langle\Phi_{0}|\hat{A}_{ij}|\Phi_{0}\rangle$ and $B_{ij}=\langle\Phi_{0}|\hat{B}_{ij}|\Phi_{0}\rangle$ are mean-field parameters computed in the gapped boson vacuum $|\Phi_{0}\rangle$ at $T=0$. For a system with $n_{u}$ sites per unit cell, this GS is obtained by diagonalization of the $(2n_{u})\times(2n_{u})$ $q$-dependent Hamiltonian matrix expressed in the Fourier space – after performing a Choleski decomposition Toth and Lake (2015); Halimeh and Punk (2016); Merino and Ralko (2018) – on a Brillouin zone of linear size $l$ containing $l\times l$ momenta ($n_{u}\times l\times l$ lattice sites). If one considers translationally invariant solutions, only one Lagrange multiplier per sublattice is needed, $\{\mu_{s}\}_{s=0,..,n_{u}-1}$, as well as $4n_{u}$ mean-field parameters, as shown in Fig. 3 (left) for $n_{u}=6$. No difference was noted on the phase diagram when using bigger unit cells up to $n_{u}=24$. We choose the convention $S=(\sqrt{3}-1)/2<1/2$ so that our theory naturally fulfils the sum rule $\sum_{i}\langle\hat{\mathbf{S}}_{0}\cdot\hat{\mathbf{S}}_{i}\rangle=NS(S+1)$ without the extra $3/2$ factor appearing for $S=1/2$ Messio et al. (2017); Bauer and Fjærestad (2017).

In order to compute the value of the mean-field parameters minimizing the energy of the system, and describing the GS for different values of $x$, we have applied the following method. First, we start from any given Ansatz $\{A_{ij},B_{ij}\}$. Then, we search for the set of Lagrange multipliers $\{\mu_{s}\}$ fulfilling the constraint by using a least square minimization. After diagonalizing the Hamiltonian, we compute the new set of parameters $\{A_{ij},B_{ij}\}$ from $|\Phi_{0}\rangle$. Finally, a new Hamiltonian is constructed and the self-consistent procedure is repeated until convergence. We have used an arbitrary tolerance of at least $10^{-9}$ on the mean-field parameters and $10^{-12}$ on the energy. This unrestricted Cabra et al. (2011); Zhang and Lamas (2013) method allows to find solutions without imposing any ad-hoc Ansatz to be miminized, that could be too restrictive and overlook important emergent features such as our reported TNSL. Note that this algorithm is derivative-free in contrast to usual SBMFT approaches and naturally deals with complex mean-field parameters, crucial in the study of gauge-invariant flux induced topological properties.

III Plaquette and ferrimagnet

For small values of $x$, our GS is composed of resonating singlets on the square plaquettes formed by the $u$ sites. At the mean-field level, $v$ spins are uncorrelated leading to a large degeneracy which is difficult to treat within our approach. However, it has been shown that virtual plaquette excitations induce an effective coupling between $v$ sites, responsible for a staggered valence-bond order Rousochatzakis et al. (2013). Such perturbative effects are beyond SBMFT, but can be mimicked by introducing a coupling $J_{p}$ between $v$ sites across each plaquette, which allows us to overcome convergence issues induced by a large degeneracy of the mean-field solution at these parameters. Taking the limit $J_{p}\to 0$, this staggered state remains energetically favored up to $x\simeq 0.37(1)$, point at which a first order phase transition to the I1 incommensurate order occurs. Note that the finite gap induced by $J_{p}$ vanishes as $J_{p}\to 0$ [Fig. 2]. To confirm this boundary, we have used the fact that the plaquette phase is made of local modes – and thus with a flat dispersion relation – to construct a flatness parameter $\sigma$ as the average of standard deviations of the bands. Close to the boundary, convergence issues due to Bose condensations and degeneracies make it difficult to measure $\sigma$, but a linear extrapolation [gray line in Fig. 2] confirms that $\sigma$ reaches zero when $x\simeq 0.38$.

For $x\gg 1$, the GS is a Lieb ferrimagnet Lieb (1989); Rousochatzakis et al. (2013); Morita and Tohyama (2018). Spins on sites $u$ and $v$ are pointing in opposite directions resulting in a finite magnetization of $M=1/3$. Classically, this phase extends to $x=2$ (resp. 1) for Heisenberg (resp. Ising) spins Rousochatzakis et al. (2013); Pohle et al. (2016). Our SBMFT is unfortunately not designed for describing configurations outside the $S_{z}=0$ sector since it cannot fulfill the boson constraint even on average Feldner et al. (2011). To estimate the boundary between I2 and the ferrimagnetic phase, we have measured the level crossing between the extrapolation of the I2 energy measured by SBMFT Feldner et al. (2011) and the classical energy of the ferrimagnetic phase given by $(4S^{2}J_{1}-8S^{2}J_{2})/6$ (per site). They cross around a coherent value of $x\approx 2.04$, but the presence of a small intervening phase – hinted at but not described in Refs. Rousochatzakis et al. (2013); Morita and Tohyama (2018) – just below the ferrimagnet cannot be ruled out because of convergence difficulties at large system sizes for $x\sim 1.5-2.0$.

IV Incommensurate magnetic orders

Between the plaquette and ferrimagnetic phases, two ordered phases with Bose-Einstein condensation of the Schwinger bosons emerge at wave vector $\mathbf{Q}$ [Fig. 2]. By spontaneous symmetry breaking, the ordering wavevector takes the form $(0,Q)$ or $(Q,0)$. To extract the values of $Q$ and the energy gap $\Delta$ at the thermodynamic limit, we have considered the asymptotic mean-field Ansatz at very large system size up to $62424$ sites, before minimizing the dispersion relation in the continuum (see Appendix A). A gapless phase with incommensurate ${\bf Q}$ vector is obtained [Fig. 2] which is not accessible by ED with periodic boundary conditions Rousochatzakis et al. (2013); Morita and Tohyama (2018) because of the size limitation ($\sim 30$ sites). Note however that in principle such incommensurability could be approached by ED considering twisted boundary conditions Sindzingre et al. (2010). The transition points of these two gapless incommensurate phases with the topological state are defined by the opening of the gap at $x\simeq 0.84$ and $x\simeq 1.27$ respectively.

V Topological nematic phase

But the most remarkable outcome of this unrestricted SBMFT analysis is certainly the presence of an extended topological nematic spin liquid (TNSL) for $x\in[0.84,1.27]$, as explained below.

Within SBMFT, gauge invariant quantities called Wilson loops (WLs) Wang and Vishwanath (2006); Tchernyshyov et al. (2006) can be constructed from $A_{ij}$ and $B_{ij}$ parameters. We find two types of elementary plaquettes of 6 and 3 sites at which a non-trivial gauge-invariant flux of $\pi$ emerges. On Fig. 3, this flux corresponds respectively to the phase $\Phi_{A}$ of $A_{ij}(-A^{*}_{jk})A_{kl}(-A^{*}_{lm})A_{mn}(-A^{*}_{ni})$ (middle) and the phase $\Phi_{B}$ of $B_{ij}B_{jk}B_{ki}$ (right). Note that $\Phi_{A}$ is invariant under a $\pi$-rotation of the WL [see Fig. 3 (middle)]. Such non-zero fluxes suggest the possibility of non-trivial order and help in discriminating different phases otherwise similar by their symmetries. By defining two winding WLs around the lattice torus Bieri et al. (2016), it is possible to evidence the topological character of the phase. A fourfold topological degeneracy is obtained as expected for a $\mathbb{Z}_{2}$ spin-liquid state (see Appendix B for details). This degeneracy is enhanced to eightfold when including the broken rotation symmetry of the nematicity (see below). It is important to note that the physical properties and local WLs fluxes remain unchanged in each topological sector.

The gap of the TNSL reflects the presence of singlets on all bonds and indicates an absence of long-range dipolar order. However, while being homogeneous on the $J_{1}$ bonds (resonating plaquettes), singlets are weakly inhomogeneous on the $J_{2}$ ones. The latter on the SKL form a decorated square lattice, with horizontal and vertical zig-zag lines. In the TNSL, the singlets on these horizontal lines have a different amplitudes than on the vertical ones, breaking $\pi/2-$rotational symmetry. This nematicity can be measured by considering a directional order parameter Zhang and Lamas (2013) defined as $\Psi=|\sum_{d=1,4}(-1)^{d}\langle\hat{\mathbf{S}}_{\rm v}\cdot\hat{\mathbf{S}}_{\textrm{u}_{d}}\rangle|,$ where $\langle...\rangle$ averages over all $J_{2}$ bonds connecting a $v$ site and one of its four $u_{d}$ nearest-neighbors. $\Psi$ takes a small but non-zero value; e.g. $\Psi=0.0014$ at $x=1$.

However, since the mean-field decoupling prevents access to singlet-singlet correlations, it is relevant to compare the TNSL with other possible scenarios. First, the weak nematicity of the TNSL motivates to consider a similar topological Anstaz but where rotational symmetry is restored, labeled TSL. Then, making contact with previous works, a “pinwheel” valence-bond crystal (PVBC) had also been proposed at $x=1$ Rousochatzakis et al. (2013), later challenged by a resonating 6-site loop (R6L) state Ralko and Rousochatzakis (2015) when going beyond the nearest-neighbor valence bond basis. Interestingly, the SBMFT counterparts of these 3 states are metastable solutions of our algorithm, with excitation energies: magenta$\Delta E_{\textrm{TSL}}=0.0004$ $<\Delta E_{\textrm{PVBC}}=0.0171<\Delta E_{\textrm{R6L}}=0.0281$ at $x=1$. In addition, the amplitudes of $\langle{\bf S_{i}}\cdot{\bf S_{j}}\rangle$ are similar between $J_{1}$ and $J_{2}$ bonds for both the TNSL and TSL Ansätze going down to 0.1% difference at $x=1$. Interestingly at the point $x\approx 1.143$, the nematicity is lost and the ground state of our system thus becomes a quantum spin liquid. Those energy differences and similarities in singlet amplitudes are not in favor of a valence bond crystal (PVBC or R6L) but are consistent with a quantum spin liquid with hidden nematic order, as defined in Ref. Grover et al. (2010).

VI Role of the anisotropy

As mentioned in the introduction, a square-kagome compound has recently been announced in Refs. Fujihala et al. (2017); Morita and Tohyama (2018). While its parametrisation is still speculative, the authors have pointed out that treating an additional distortion between nonequivalent $J_{2}$ bonds could be essential. We have then challenged our TNSL phase against such a distorsion parametrized by two couplings $J_{2}$ and $J_{2}^{\prime}$ [Fig. 5 left]. While the full phase diagram of this model is beyond the scope of the present work (see Ref. Morita and Tohyama (2018) for its exact-diagonalization contour), the TNSL is shown to be stable over an extended region of the $(J_{2},J_{2}^{\prime})$ plane [Fig. 5 right] making it a promising candidate for future materials.

VII Dynamical structure factors

A powerful experimental probe of magnetic phases is inelastic neutron scattering, readily available by our method via the spin-spin dynamical structure factors $S(\mathbf{q},\omega)$ [Fig. 4] (see Halimeh and Punk (2016); Merino and Ralko (2018) for technical details). Note that we adopt the same normalization as in Halimeh and Punk (2016), setting the maximum of $S(\mathbf{q},\omega)$ to unity. At small $x$, $S(\mathbf{q},\omega)$ reveals flat bands consistent with the localized resonating square plaquettes. The three other phases I1, I2 and TNSL have quite similar dynamical responses, with a clear broken symmetry indicated by a spectral weight imbalance between the two paths in the Brillouin zone. This similarity, together with the smooth evolution of $Q$ across I1 and I2 in Fig. 2 and the continuous gap closing, suggest that TNSL originates from the quantum melting of the incommensurate orders. We stress that a few key differences (gap, intensities) would make them distinguishable in experiments.

VIII Outlook

The present paper opens the path for several promising directions. We propose indeed a rare example of a topological nematic spin liquid Grover et al. (2010) emerging naturally from an antiferromagnetic insulator, possibly connected to a quantum spin liquid – see the vanishing nematic order parameter at $x\approx 1.143$. Alternative methods, such as DMRG, would be welcome to confirm our scenario. Keeping in mind the kagome debate, is the gap of the TNSL robust beyond mean-field ? Can we further describe the quantum melting of the incommensurate phases ? These questions are particularly relevant since the TNSL is stable against realistic microscopic perturbations Fujihala et al. (2017); Morita and Tohyama (2018) and the SKL appears to be realisable in optical lattices Glaetzle et al. (2014). Square-Kagome (artificial) materials offer thus an exciting direction for exotic topological quantum states.

Acknowledgements – We thank Owen Benton, Jaime Merino, Laura Messio, Ioannis Rousochatzakis and Nic Shannon for fruitful discussions. The authors thank anonymous referees for insightful questions that helped in improving the manuscript. TL & AR (resp. LJ) acknowledge hospitality from the LOMA in Bordeaux (resp. the Néel Institute in Grenoble). This work was supported by the “Agence Nationale de la Recherche” under Grant No. ANR-18-CE30-0011-01 and travel grants from the GdR MEETICC.

Appendix A Reaching the thermodynamic limit

Due to the incommensurate properties of the phases of the model, strong size effects appear. It is thus quite delicate to properly extract the thermodynamic limit of some quantities such as the ordering wavevector $\mathbf{Q}_{\text{soft}}$ or the gap $\Delta$. Thankfully, the mean-field parameters, that are integrated quantities (summed over the full BZ), are found to be much more stable in function of the size. For example, in the TNSL phase, no significant variations can be noticed for sizes larger than $l=78$ for $n_{u}=6$. We have thus exploited this advantage to extract the thermodynamic limit of the mean-field parameters instead of the physical quantities.

Minimizing the dispersion relation for these fixed Ansätze in the continuum of the thermodynamic limit gives access to $\mathbf{Q}_{\text{soft}}$ and $\Delta$ simultaneously, and circumventes the convergence issues. To illustrate the efficiency of the procedure, the value of the gap for different system sizes up to 4860000 sites is displayed on Fig. 6 (blue dots) and compared with the extracted valuemagenta in the thermodynamic limit (gray line). The oscillations reflect the incommensurability of $\mathbf{Q}_{\text{soft}}$ minimizing the dispersion relation. We note that the thermodynamic limit of $\Delta$ corresponds to the minimum of the oscillations.

This method has proven very useful to determine the precise location of the gap opening in function of $J_{2}$. However, numerical convergence is quite difficult beyond $l=102$ in the gapless phases. As a result, the gap closes in a small and finite area in Fourier space instead of a single Bragg point. Our finite-size scalings have confirmed that this area converge to a single point at the thermodynamic limit.

Appendix B Topological degeneracy

The topological nature of the Topological Nematic Spin Liquid (TNSL) can be evidenced by following the flux insertion procedure proposed in Bieri et al Bieri et al. (2016). The idea is that within the parton construction of quantum spin liquids, piercing the torus by additional gauge fluxes leads to locally indistinguishable degenerate states Bieri et al. (2016). For lattice gauge fields though, the flux insertion has to follow certain constraints, and a ${\mathbb{Z}}_{2}$ quantum spin liquid only authorizes the insertion of a flux of $\pi$.

To construct these states following the constraint, one has to define “cut” lines winding around the lattice in the two natural directions as depicted in Fig. 7, and to change the sign of any mean-field parameter crossing the line Bieri et al. (2016). This corresponds to two operators, $\hat{T}_{1}$ and $\hat{T}_{2}$ verifying $\hat{T}_{1,2}{\cal O}_{i,j}=\mp{\cal O}_{i,j}$ whether the bond $(i,j)$ crosses a cut-line or not. This construction is similar to the one used to determine topological order in quantum dimer models Bieri et al. (2016), and is easy to apply in real space where a single cut line can be defined. In our Fourier representation however, defining a cut line on the unit-cell implies $l$ copies in real space. As a result, working with an enlarged unit-cell of $n_{u}=24$ sites as depicted in Fig. 7, we can generate $l$ cut lines in each direction for a system of linear size $l$. We clearly see that since the flux insertion can only be of $\pi$ per cut line, any even $l$ lattice size will have a trivial flux in this construction, while odd $l$ can have non-trivial $\pi$-flux inserted.

Now, starting from the TNSL described in the main text and applying operators $\hat{T}_{1,2}$ on it, 4 different Ansätze can be constructed, as shown in Fig. 7. For each of them, we have calculated through the SBMFT the energy and we have verified that the local Wilson loops defined in the text as well as the physical parameters were unchanged. In addition, we have defined two winding Wilson loops based on $A_{ij}$ parameters going in the two lattice directions as depicted in Fig. 8. The obtained fluxes are non-local and reflects the topological nature of the phase. It is displayed in unit of $\pi$ for each Ansatz in Fig. 7.

As mentioned previously, for a given system size $l$, the total flux insertion of a winding WL can only be [math] or $l\pi$. It is thus necessarily trivial (mod $2\pi$) for even $l$; the four solutions of Fig. 7 are degenerate but do not belong to different topological sectors. Hence, within this approach, one needs to consider odd values of $l$ to establish the topological degeneracy. For finite and odd $l$ system size, the energies are split in three groups, TS(0,0), TS(1,1) and TS(1,0)/(0,1), as displayed in Fig. 9. Finite size scaling shows that these energies converge to the same ground state energy as the one obtained for even values of $l$, defining a thermodynamic limit as $l\to\infty$ (light gray line). This is the fourfold topological degeneracy Misguich (2012) within the ergodicity breaking of the nematic order. This degeneracy is enhanced to eightfold when including the rotation symmetry of the two nematic states.