Topological Nematic Spin Liquid on the Square-Kagome Lattice

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.


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 [1][2][3][4][5]. 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 [6][7][8][9][10][11], a gapped Z 2 spin liquid [12][13][14], or something else [8,11]. 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 [15] [ 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 [11,16]. Adding to that the recent announcement of a spin−1/2 material with distorted square-kagome geometry [17][18][19], and the theoretical proposal to realise the SKL in optical lattices [20], it has become a timely question to understand the quantum GS of this model, moving beyond previous works on exact diagonalization (ED) [17,[21][22][23] and perturbative methods [15,16,24].
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 , relevant to experiments, will be considered at the end of the paper. Square-kagome lattice with its two nonequivalent sets of sites, u and v, respectively 4 (blue) and 2 (red) per unit cell, and the (e1, e2) lattice translation vectors (left). Spin coupling J1 links neighboring u sites on square plaquettes while J2 connects (u, v) sites on octagons. Two symmetry equivalent paths in the Brillouin zone connecting high symmetry points (right).

The Hamiltonian then reads
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 [25] 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. 23. Then, as x closes to 1 from both sides, we report the appearance of two gapless incommensurate orders, I1 and I2. Eventually when x ∼ 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 [17] and provide inelastic neutron scattering signatures of the reported phases. Ferri. They extend until energy level crossings (located at the dashed vertical lines) with the plaquette phase at x 0.37 and the ferrimagnet at x 2.04 through a first order transition scenario. This is corroborated by σ collapsing down to zero at the transition point between I1 and the plaquette phase.

II. SCHWINGER BOSON MEAN-FIELD THEORY
SBMFT has proven successful in a variety of models: to analytically study antiferromagnetic GS [26], to classify spin liquids using projective symmetry groups [27][28][29] and even to establish phase diagrams [30][31][32]. The method consists in introducing two Schwinger bosons,b † i↑ andb † i↓ , on each site i of the lattice through the following spin mappingŜ where σ are the Pauli matrices and α, β = ↑, ↓. For this mapping to be correct, an on-site boson constraint b † i↑b i↑ +b † i↓b i↓ = 2S has to be fulfilled. This constraint being very hard to implement exactly, it is thus imposed on average by introducing Lagrange multipliers {µ i }. Two SU(2)-invariant bond operators [33], the singlet operator on the oriented bond (i → j)Â ij = (b i↑bj↓ −b i↓bj↑ )/2 and the boson hopping operatorB ij = (b † i↑b j↑ +b † i↓b j↓ )/2, are introduced. They are respectively favored in disordered and ordered phases. Following a mean-field decoupling of these operators, Eq. 1 becomes: where A ij = Φ 0 |Â ij |Φ 0 and B ij = Φ 0 |B ij |Φ 0 are mean-field parameters computed in the gapped boson vacuum |Φ 0 at T = 0. For a system with n u sites per unit cell, this GS is obtained by diagonalization of the (2n u ) × (2n u ) q-dependent Hamiltonian matrix expressed in the Fourier space -after performing a Choleski decomposition [32,34,35] -on a Brillouin zone of linear size l containing l × l momenta (n u × l × l lattice sites). If one considers translationally invariant solutions, only one Lagrange multiplier per sublattice is needed, {µ s } s=0,..,nu−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 = ( √ 3 − 1)/2 < 1/2 so that our theory naturally fulfils the sum rule i Ŝ 0 ·Ŝ i = N S(S + 1) without the extra 3/2 factor appearing for S = 1/2 [36,37].
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 {µ 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 |Φ 0 . Finally, a new Hamiltonian is constructed and the selfconsistent 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 [30,31] 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 [23]. 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 → 0, this staggered state remains energetically favored up to x 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 → 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 σ 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 σ, but a linear extrapolation [gray line in Fig. 2] confirms that σ reaches zero when x 0.38.
For x 1, the GS is a Lieb ferrimagnet [17,23,38]. 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 [23,39]. 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 [40]. 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 [40] 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 ≈ 2.04, but the presence of a small intervening phase -hinted at but not described in Refs. [17,23] -just below the ferrimagnet cannot be ruled out because of convergence difficulties at large system sizes for x ∼ 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 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 ∆ 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 Q vector is obtained [ Fig. 2] which is not accessible by ED with periodic boundary conditions [17,23] because of the size limitation (∼ 30 sites). Note however that in principle such incommensurability could be approached by ED considering twisted boundary conditions [41]. The transition points of these two gapless incommensurate phases with the topological state are defined by the opening of the gap at x 0.84 and x 1.27 respectively. Dashed arrows indicate an extra phase of π on the corresponding mean-field parameter, Aij (middle) and Bij (right). The blue shaded loop represents the smallest WL with non-trivial flux ΦA = π and the red shaded loops the 4 smallest with ΦB = π. Note that in the Aij WL case (middle), π-rotation also leads to a flux of ΦA = π.

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 ∈ [0.84, 1.27], as explained below.
Within SBMFT, gauge invariant quantities called Wilson loops (WLs) [27,42] 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 nontrivial gauge-invariant flux of π emerges. On Fig. 3, this flux corresponds respectively to the phase Φ A of A ij (−A * jk )A kl (−A * lm )A mn (−A * ni ) (middle) and the phase Φ B of B ij B jk B ki (right). Note that Φ A is invariant under a π-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 [43], it is possible to evidence the topological character of the phase. A fourfold topological degeneracy is obtained as expected for a 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 π/2−rotational symmetry. This nematicity can be measured by considering a directional order parameter [31] where ... averages over all J 2 bonds connecting a v site and one of its four u d nearest-neighbors. Ψ takes a small but non-zero value; e.g. Ψ = 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 [23], later challenged by a resonating 6-site loop (R6L) state [16] 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∆E TSL = 0.0004 < ∆E PVBC = 0.0171 < ∆E R6L = 0.0281 at x = 1. In addition, the amplitudes of S i · S j 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 ≈ 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. [44].

VI. ROLE OF THE ANISOTROPY
As mentioned in the introduction, a square-kagome compound has recently been announced in Refs. [17,19]. 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 [ Fig. 5 left]. While the full phase diagram of this model is beyond the scope of the present work (see Ref. [17] for its exact-diagonalization contour), the TNSL is shown to be stable over an extended region of the (J 2 , J 2 ) plane [ Fig. 5 right] making it a promising candidate for future materials.  Fig. 2. We have only tested the stability of our TNSL phase against the disproportion between J2 and J 2 . We have tracked down the points (red circles) at which the nematicity is lost (within our mesh precision symbolized by the thick red lines), indicating the boundaries of the TNSL.

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(q, ω) [ Fig. 4] (see [32,35] for technical details). Note that we adopt the same normalization as in [35], setting the maximum of S(q, ω) to unity. At small x, S(q, ω) 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 [44] emerging naturally from an antiferromagnetic insulator, possibly connected to a quantum spin liquid -see the vanishing nematic order parameter at x ≈ 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 [17,19] and the SKL appears to be realisable in optical lattices [20]. 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 Q soft or the gap ∆. 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 Q soft and ∆ 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 Q soft minimizing the dispersion relation. We note that the thermodynamic limit of ∆ 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 The topological nature of the Topological Nematic Spin Liquid (TNSL) can be evidenced by following the flux insertion procedure proposed in Bieri et al [43]. 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 [43]. For lattice gauge fields though, the flux insertion has to follow certain constraints, and a Z 2 quantum spin liquid only authorizes the insertion of a flux of π.
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 [43]. This corresponds to two operators,T 1 andT 2 verifyingT 1,2 O i,j = ∓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 [43], 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 π per cut line, any even l lattice size will have a trivial flux in this construction, while odd l can have non-trivial π-flux inserted. (0,1) The four topological sectors obtained on the 24site unit cell system with odd values of the system size l. Thick dashed bonds correspond to a phase of −1 on the meanfield parameters Aij. One can go from one gauge definition to another one by reversing all signs of the parameters that cross the winding cut lines (red and blue) by applyingT1,2 operator. The topological sector is defined by the presence or not of a π flux along a winding WL making the complete turn of the torus in the natural directions of the lattice (see text and Fig. 8). Now, starting from the TNSL described in the main text and applying operatorsT 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 nonlocal and reflects the topological nature of the phase. It is displayed in unit of π 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 0 or lπ. It is thus necessarily trivial (mod 2π) 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 → ∞ (light gray line). This is the fourfold topological degeneracy [45] within the ergodicity breaking of the nematic order. This degeneracy is enhanced to eightfold when including the rotation symmetry of the two nematic states.  9. Evidence of the topological degeneracy at the thermodynamic limit. Energy of the 4 Ansätze defined in Fig. 7 versus 1/l 2 . The thermodynamic limit is obtained for a 6site unit-cell system of size l = 900 where no fluctuations can be observed. For l even (blue points), the 4 Ansätze are exactly degenerate and their winding WL fluxes are (0, 0). For l odd, the four topological sectors (0, 0) (red), (1, 1) (black), (0, 1) and (1, 0) (both gray) gradually converge to the same thermodynamic limit as 1/l 2 → 0.