Z2 flux binding to higher-spin impurities in the Kitaev spin liquid

Z2 flux binding to higher-spin impurities in the Kitaev spin liquid

Introduction

Magnetic impurities in strongly correlated electron systems are valuable tools for probing hidden physical phenomena. A well-known example is the Kondo effect, where the scattering of conducting electrons by magnetic impurities in metals or quantum dots reveals crucial insights into the low-energy physics of both the bulk material and the impurities1,2,3,4,5. This effect, characterized by the screening of the magnetic impurity by the conduction electrons leading to the formation of a Kondo singlet, has profound implications for understanding many-body interactions and has been extensively studied in both theoretical and experimental contexts.

The study of magnetic impurities extends to a variety of systems beyond conventional metals. In topological insulators, for instance, magnetic impurities can break time-reversal symmetry, leading to the opening of a gap at the Dirac point on the surface states and potentially inducing novel magnetic phases6,7. These systems provide a rich playground for exploring the interplay between magnetism and topology, with potential applications in spintronics and quantum computing.

In low-dimensional spin systems such as quantum spin liquids (QSLs), the introduction of impurities can reveal even more exotic phenomena8,9,10,11. QSLs, which are characterized by a lack of conventional magnetic order even at zero temperature due to strong quantum fluctuations, offer a unique environment where impurities can induce localized excitations and modify the emergent gauge fields.

In the context of the Kitaev spin liquid (KSL) model12 -a paradigmatic example of a two-dimensional QSL with fractionalized excitations and emergent gauge fields—introducing impurities, whether magnetic (spin-S sites) or non-magnetic (vacancies)13,14,15,16,17,18,19,20,21,22, can lead to various novel phenomena. These include localized bound states15,19,23,24, flux binding effects13,14,17,18,22, and modifications in the system’s topological nature20,22. Such effects not only provide new insights into impurity physics in QSLs but also enhance our understanding of their overall behavior. Due to the localized nature of the low-energy fractionalized excitations in the above phenomena, characterization of KSL based on the signatures in scanning tunneling microscopy (STM) has been proposed in various theoretical works21,23,24,25,26,27,28,29,30,31,32,33. In particular, the above proposals show that the inelastic tunneling spectroscopy can access the real-space spin-spin correlation function of the Kitaev model, which is very sensitive to defects, open edges, and local flux structures.

So far, two types of local impurities have been relatively well-studied in the Kitaev model. The first type is vacancies. It has been demonstrated that vacancies in the Kitaev model lead to almost zero-energy localized bound states and flux-binding effects13,14,15,18,19,34, which can potentially be probed by thermodynamics35,36 and STM21,23,24. The second type is spin-S impurities, which are coupled to KSL at a given site via Kondo coupling. The studies of Kondo impurities16,17,37 have highlighted several remarkable properties of the Kondo effect in the Kitaev model. In the presence of a spin-1/2 Kondo impurity, the fluxes in the three plaquettes adjacent to the impurity site are no longer individually conserved. However, their product (the flux in the impurity plaquette) and all outer fluxes remain conserved17. Furthermore, a topological transition occurs from the zero-flux state to a bound-flux state attached to the impurity site as a function of Kondo coupling16,17.

In this work, we theoretically investigate the behavior of Simp = 1 and Simp = 3/2 impurities in KSL by means of numerical exact diagonalization (ED) and density matrix renormalization group (DMRG) methods as well as phenomenological models in the Majorana representation. We will focus on the case with a single magnetic impurity, although the case of multiple impurities can be straightforwardly extended.

We show that the behavior of KSL with a magnetic impurity strongly depends on whether the impurity has a half-integer or integer spin. This dependence, which we demonstrate by considering two cases of magnetic impurities, Simp = 1 and Simp = 3/2, echoes the recent findings by Ma38, which shows that the nature of the spin-SZ2 KSL differs based on whether the spin is integer or half-integer. While Ma’s work introduces 2S-flavor Majorana representation for the pure spin-S KSL and identifies the Z2 gauge fluxes as conserved quantities, we investigate the behavior of a spin-S impurity embedded within a spin-1/2 Kitaev spin liquid, focusing on the impurity magnetization and the ground-state flux sector as functions of varying impurity coupling strength and spin size. This mixed-spin KSL system, though less explored in the literature, holds significant potential for realization in Kitaev materials containing magnetic impurities.

Similarly to vacancies, the magnetic impurities can bind Z2-fluxes in the lattice. We show that varying the coupling of spin-S impurity with the surrounding spin-1/2 KSL can drive a phase transition between bound-flux and zero-flux sectors. Furthermore, the point at which this transition occurs depends on the magnitude of Simp. This is significant because, in the presence of a time-reversal symmetry-breaking magnetic field, each Z2 flux can bind a Majorana zero mode, thereby realizing an Ising anyon governed by non-Abelian statistics12. Therefore, demonstrating that magnetic impurities can trap Z2 fluxes provides a pathway to realizing Ising anyons in these systems.

Results

Model

To describe the Kitaev model with magnetic impurities, we begin with the Hamiltonian:

$$H=-Jsumlimits_{mathop {j,k notin {{Lambda}},}limits_{langle {jk}rangle_{mu }}}{S}_{j}^{mu }{S}_{k}^{mu }-gsumlimits_{mathop {{j}{in} {{Lambda }},{k}{notin} {{Lambda }},}limits_{langle jkrangle_{mu }} }{tilde{S}}_{j}^{mu }{S}_{k}^{mu }quad (mu =x,y,z).$$
(1)

Here, Sμ = σμ/2 represents the μ-component of the spin-1/2 operator, with σ denoting the Pauli matrices. ({tilde{S}}^{mu }) represents the μ-component of the spin-S operator of the impurity. The set of impurity sites is denoted by Λ. J > 0 denotes the ferromagnetic coupling strength between spin-1/2 operators away from the impurity, while g > 0 denotes the ferromagnetic coupling strength between the impurity spins and the spin-1/2 operators of the original Kitaev model. We mainly consider the case that a single impurity is not located on the edges.

Since the interaction on each honeycomb bond remains Kitaev-like, we can define a triple-plaquette flux operator in the vicinity of the magnetic impurity as:

$${W}_{I}equiv {2}^{12}{S}_{1}^{x}{S}_{2}^{x}{S}_{3}^{y}{S}_{4}^{z}{S}_{5}^{z}{S}_{6}^{z}{S}_{7}^{x}{S}_{8}^{y}{S}_{9}^{y}{S}_{10}^{y}{S}_{11}^{z}{S}_{12}^{x},$$
(2)

where the product is taken over all bonds forming a 12-site plaquette around the impurity site labeled as “0” (see Fig. 1a). This operator captures the flux configuration within the three plaquettes surrounding the magnetic impurity and allows us to analyze the effect of the impurity on the local flux dynamics. A Z2 flux operator at a plaquette p in the bulk, where there is no impurity spin, is symbolically expressed as ({W}_{p}={2}^{6}{prod }_{jin p}{S}_{j}^{mu }) in the spin-1/2 basis12. Both Wp and WI commute with the Hamiltonian H and with each other, taking ± 1 as their eigenvalues, respectively. Consequently, the total Hilbert space is divided into individual flux sector subspaces: ({mathcal{L}},{ = bigoplus }_{{w}_{{p}_{1}},{w}_{{p}_{2}},cdots ,{w}_{I}},{{mathcal{L}}}_{{w}_{{p}_{1}},{w}_{{p}_{2}},cdots ,{w}_{I}}). We refer to the sector with wI = − 1 as the bound-flux sector and the sector with wI = + 1 as the zero-flux sector, when all plaquettes not in the vicinity of the impurity have wp = 1, as introduced in a vacancy case13.

Fig. 1: Integer/half-integer dependence of magnetic impurity.
Z2 flux binding to higher-spin impurities in the Kitaev spin liquid

a A spin-S impurity embedded in the Kitaev spin liquid. The inset shows site labels around the impurity site (0) and the plaquette labels used in Eq. (5). b, c Typical expectation values of internal flux operators ({{mathbb{W}}}_{a}) including the impurity site (0) in the bound-flux sector (left) and the zero-flux sector (right) for Simp = 1 and Simp = 3/2, respectively. d Local magnetization of the impurity spin embedded in spin-1/2 Kitaev spin liquid, calculated via spin-space exact diagonalization in a 24-site cluster.

Full size image

Internal plaquette operators

Around the impurity, the three adjacent plaquette operators must incorporate the higher-spin operator ({tilde{S}}_{0}^{alpha }), such that they are constructed by the unitary operators of π-rotation, as introduced by Baskaran et al.39. We first quote the basic properties of these unitary operators40:

$${tilde{R}}_{j}^{alpha }={e}^{ipi {tilde{S}}_{j}^{alpha }},quad {left({tilde{R}}_{j}^{alpha }right)}^{2}={(-1)}^{2{tilde{S}}_{j}},quad {tilde{R}}_{j}^{alpha }{tilde{R}}_{j}^{beta }={(-1)}^{2{tilde{S}}_{j}}{tilde{R}}_{j}^{beta }{tilde{R}}_{j}^{alpha },quad {tilde{R}}_{j}^{alpha }{tilde{R}}_{j}^{beta }={tilde{R}}_{j}^{gamma },$$
(3)

where αβγ, and (α, β, γ) (x, y, z) obey cyclic permutations. The above relations are defined on the same site j, otherwise the operators simply commute on different sites. For a general hexagonal plaquette with arbitrary spin on the corners, one can define the plaquette operator as

$${{mathbb{W}}}_{p}=prod _{jin p}{tilde{R}}_{j}^{{alpha }_{j}},quad {({{mathbb{W}}}_{p})}^{2}=prod _{jin p}{(-1)}^{2{tilde{S}}_{j}}.$$
(4)

Specifically, for the three plaquettes shown in Fig. 1a, we can write down (the tilde is dropped for spin-1/2 sites for clarity)

$$begin{array}{rcl}&&{{mathbb{W}}}_{{I}_{1}}={R}_{6}^{z},{R}_{7}^{x},{R}_{8}^{y},{R}_{9}^{z},{tilde{R}}_{0}^{x},{R}_{5}^{y},\ &&{{mathbb{W}}}_{{I}_{2}}={R}_{4}^{z},{R}_{5}^{x},{tilde{R}}_{0}^{y},{R}_{1}^{z},{R}_{2}^{x},{R}_{3}^{y},\ &&{{mathbb{W}}}_{{I}_{3}}={tilde{R}}_{0}^{z},{R}_{9}^{x},{R}_{10}^{y},{R}_{11}^{z},{R}_{12}^{x},{R}_{1}^{y},end{array}$$
(5)

and then it is straightforward to show that ({{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}={(-1)}^{2{S}_{{rm{imp}}}+2{S}_{5}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{1}}). Therefore, one can see that the mutual commutation relation between the internal plaquette operators depends on the size of the spins shared by the two hexagons, which is consistent with the previous work on the mixed-spin Kitaev model41. In our impurity model, only the impurity site can have a general spin size and all the rest are spin-1/2. This leads to the important properties of the internal plaquette operators:

$$left{begin{array}{ll}left[{{mathbb{W}}}_{a},,{{mathbb{W}}}_{b}right]=0,quad {({{mathbb{W}}}_{a})}^{2}=+1,quad &{mathrm{if}},,{S}_{{rm{imp}}},,{text{is}}, {text{a}}, {text{half}}hbox{-}{text{integer}},,\ left{{{mathbb{W}}}_{a},,{{mathbb{W}}}_{b}right}=0,,,,{({{mathbb{W}}}_{a})}^{2}=-1,quad &{mathrm{if}},,{S}_{{rm{imp}}},,{text{is}}, {text{an}},{text{integer}},,end{array}right.$$
(6)

where ab and (a, b) (I1, I2, I3).

Note that a triple-plaquette operator and internal plaquette operators are directly related. They satisfy the following relation:

$${{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{3}}={(-1)}^{2{S}_{{rm{imp}}}}{W}_{I}.$$
(7)

For a half-integer spin impurity, the conserved value wI (= ± 1) can be decomposed into the product of ({w}_{{I}_{1}},{w}_{{I}_{2}},) and ({w}_{{I}_{3}}), each taking ± 1. Thus, in the case of the Simp = 3/2, the triple-plaquette operator WI and three internal plaquette operators ({{mathbb{W}}}_{a}) are related by ({W}_{I}=-{{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{3}}). There are four possible configurations for (({w}_{{I}_{1}},{w}_{{I}_{2}},{w}_{{I}_{3}})=(+1,+1,+1),(+1,-1,-1),(-1,+1,-1)) and (− 1, − 1, + 1) that can realize the bound-flux sector wI = −1 (see the left panel of Fig. 1c). Similarly, there are four degenerate internal flux configurations that can realize the zero-flux sector (see the right panel of Fig. 1c). In contrast, for Simp = 1, the internal ({{mathbb{W}}}_{a}) always takes on a purely imaginary value with some real coefficient (cin {mathbb{R}}). This makes it independent of wI (see Fig. 1b). These distinct behaviors for different impurity spin values are clearly observed in the DMRG calculations, which we will discuss more later.

Integer/half-integer dependence of the impurity magnetization

Based on the algebra of the internal plaquette operators (6), one can derive an interesting integer/half-integer effect on the magnetization of the impurity spin. Notice that each internal plaquette operator contains only one of the three spin components of the impurity. This leads to the commutation/anticommutation relations between the internal plaquette operators and impurity spin operators:

$$begin{array}{rcl}&&left[{{mathbb{W}}}_{{I}_{1}},{tilde{S}}_{0}^{x}right]=left{{{mathbb{W}}}_{{I}_{1}},{tilde{S}}_{0}^{y}right}=left{{{mathbb{W}}}_{{I}_{1}},{tilde{S}}_{0}^{z}right}=0\ &&left[{{mathbb{W}}}_{{I}_{2}},{tilde{S}}_{0}^{y}right]=left{{{mathbb{W}}}_{{I}_{2}},{tilde{S}}_{0}^{z}right}=left{{{mathbb{W}}}_{{I}_{2}},{tilde{S}}_{0}^{x}right}=0\ &&left[{{mathbb{W}}}_{{I}_{3}},{tilde{S}}_{0}^{z}right]=left{{{mathbb{W}}}_{{I}_{3}},{tilde{S}}_{0}^{x}right}=left{{{mathbb{W}}}_{{I}_{3}},{tilde{S}}_{0}^{y}right}=0.end{array}$$
(8)

From Eq. (4), it follows that if the Simp is a half-integer, the square of the internal plaquette operators is one (({{mathbb{W}}}_{a}^{2}=+1)) and they all mutually commute ((left[{{mathbb{W}}}_{a},{{mathbb{W}}}_{b}right]=0)). Using these properties, we can demonstrate that:

$$leftlangle {tilde{S}}_{0}^{z}rightrangle =leftlangle {tilde{S}}_{0}^{z},{{mathbb{W}}}_{{I}_{1}}^{2}rightrangle =-leftlangle {{mathbb{W}}}_{{I}_{1}},{tilde{S}}_{0}^{z},{{mathbb{W}}}_{{I}_{1}}rightrangle =-{({w}_{{I}_{1}})}^{2}leftlangle {tilde{S}}_{0}^{z}rightrangle =-leftlangle {tilde{S}}_{0}^{z}rightrangle ,$$
(9)

which implies (langle {tilde{S}}_{0}^{z}rangle =0). The derivation is applicable to the other two components, (langle {tilde{S}}_{0}^{x}rangle =langle {tilde{S}}_{0}^{y}rangle =0), because one can always find an internal plaquette operator that anticommutes with the impurity spin component. The same argument applies to the three neighboring components (({S}_{1}^{x}), ({S}_{9}^{y}), and ({S}_{5}^{z})), as well as to all other spin-1/2 operators on the lattice. In contrast, for an integer-spin impurity, ({w}_{{I}_{a}}) takes pure imaginary values, not quantized to Z2, leading to the possibility of a non-zero local magnetization (langle {tilde{S}}_{0}^{z}rangle, ne, 0). Similarly, the μ-component (μ = x, y, z) of spin-1/2 neighboring with the impurity spin on μ-bond can have a nonzero spin moment; the other two components should be zero.

This spin-size dependence of the impurity magnetization can be confirmed through numerical exact diagonalization, as illustrated in Fig. 1d. Here we perform ED calculation for the spin Hamiltonian in a 24-site cluster on a cylinder geometry with a single impurity to evaluate the z-component of the magnetic impurity’s local spin moment. For any finite g, the local spin moment is numerically zero for both Simp = 1/2 and 3/2 cases, which aligns precisely with our analytical findings. Conversely, in the Simp = 1 case, the local moment takes finite values in the range [− 1, 1] at random due to local dynamics induced by the anticommutation relations of internal plaquette operators. The points plotted in Fig. 1d represent the results obtained in a single run. Averaging ({tilde{S}}_{0}^{z}) across multiple independent runs would yield a mean value approaching zero, highlighting a clear qualitative difference from the half-integer impurity case.

Reentrance effect of the bound-flux sector in zero and finite magnetic fields

Here, we focus on the ground-state flux configuration in the presence of the impurity. Using DMRG within the ITensorsGPU.jl package42, we numerically determine the ground-state flux configurations as a function of g/J in a 48-site cylinder with a single spin-S impurity located in the bulk [see Fig. 2a]. Figure 2b shows the transitions between the bound-flux and zero-flux sectors. Starting from the bound-flux sector at g = 0, the system undergoes a flux-sector transition under a weak but nonzero coupling g, leading to the zero-flux sector. This behavior is qualitatively consistent across all Simp cases, indicating the instability of the bound-flux sector in a “quasivacancy” problem as reported in ref.18,34. However, the transition point g1 depends on the spin size of the impurity, shifting to the left as Simp increases. Additionally, for the Simp = 3/2 impurity case, the system undergoes a second flux-sector transition at the transition point g2, resulting in the appearance of the reentrant bound-flux sector. Note again that, in both flux sectors of the Simp = 3/2 impurity case, a triple-plaquette operator WI and three internal plaquette operators ({{mathbb{W}}}_{a}) are related as ({W}_{I}=-{{mathbb{W}}}_{{I}_{1}}{{mathbb{W}}}_{{I}_{2}}{{mathbb{W}}}_{{I}_{3}}). Thus, each flux sector is fourfold degenerate.

Fig. 2: DMRG results of flux-sector transitions.
figure 2

a The bound-flux sector of a 48-site cluster used in calculations. Numbers of every plaquette in the bulk shows 〈Wp〉, while the number next to the impurity site shows 〈WI〉. b Flux-sector transition calculated by DMRG with a single spin-S impurity embedded in spin-1/2 Kitaev spin liquid. c–e Phase diagrams of 〈WI〉 as a function of the magnetic field strength h/J and the impurity-bulk coupling strength g/J for the cases Simp = 1/2, 1, and 3/2, respectively. We use the same color scale for all three subfigures: the yellow (gray) region corresponds to the bound-flux (zero-flux) sector. White dot lines denote the edges of each pixel.

Full size image

A few remarks are in order. First, the flux gap exhibits a position dependence of the impurity, possibly due to the edge effect on the flux gap as described by Feng et al.43. Specifically, the flux gap decreases in the cylinder geometry as the impurity approaches one of the edges. This positional dependence might contribute to the energy difference between the two flux sectors. For Simp = 3/2, this effect is mild and only changes the transition points g1 and g2. In contrast, for Simp = 1, we observe a significant qualitative change in the flux-sector transitions in the strong g region due to the position of the impurity, even on the 48-site cylinder. This change is influenced not only by the specifics of the numerical conditions – such as cluster shape, size, and boundary conditions, – but also by the dynamical properties around the impurity, particularly the lack of quantization of internal flux operators ({{mathbb{W}}}_{a}), as described at Eq. (6). Second, the two-impurity case exhibits behavior qualitatively similar to the single-impurity case. Specifically, every impurity spin with Simp = 3/2 binds the Z2 flux in a wide parameter range, while for Simp = 1, there is a strong position dependence on the flux-sector transitions. Third, our main findings are not specific for the 48-site cylindrical cluster. We have confirmed the same trends in other finite-size cylindrical clusters with different shapes. Numerical evidence supporting these arguments can be found in Supplementary information.

We also examined the stability of the bound-flux sector under a uniform magnetic field numerically using DMRG in the same finite-size cluster. For this analysis, we considered the Hamiltonian ({H}_{{rm{total}}}=H+h{sum }_{j,mu }{S}_{j}^{mu }), with the field applied in the [111] direction in the spin basis. We obtained three phase diagrams of 〈WI〉 as a function of g and h, shown in Fig. 2c–e, corresponding to the cases of Simp = 1/2, 1, and 3/2, respectively.

In the cases of Simp = 1/2 and 1, the bound-flux sector at nonzero g is very fragile in the presence of the external field, and this fragility is independent of the impurity’s position (Supplementary Information). In contrast, the reentrant bound-flux sector for a spin-3/2 impurity exhibits some stability for finite field strengths. This stability also depends on the coupling strength g: as g increases, the bound-flux sector tends to withstand stronger fields. This behavior can be understood by noting that the energy difference between the bound-flux and zero-flux sectors under zero magnetic field increases monotonically with g/J in the reentrant flux sector regime. This increase in energy difference provides the bound-flux sector with greater stability in the presence of magnetic fields.

Triple-plaquette analysis: Majorana representation and effective-coupling model

In this section, we aim to understand better the previous numerical findings, including the spin-size-dependent flux-sector transitions and the reentrance of the bound-flux sector. Similar phenomena have been studied in the site-diluted KSL, where a π-flux can be trapped by a true vacancy or quasivacancy13,18,21,24,34. These studies were based on the exact solution of KSL using the Majorana representation for spin-1/212. In these studies, some of the authors have shown that the low-energy modes introduced by quasivacancies are highly localized, allowing for a clear distinction in the energy spectrum between zero- and bound-flux sectors18,24. Recently, on the other hand, the Majorana parton construction was generalized to study uniform spin-S KSLs38. The essential idea of this construction is recognizing spin-S operator as 2S flavors of spin-1/2s with constraints. This motivates us to re-examine our findings in the mixed-spin KSL using the fermionic approach and the multiple-flavor representation, which we outlined in the Methods. To this end, we focus on a triple-plaquette cluster, as shown in the inset of Fig. 1a, to demonstrate that even the system on the minimal cluster can capture the key findings from the previous section.

In the following, we will discuss two fermionic models in the triple-plaquette cluster: the Simp = 3/2 impurity model with four-body Majorana interaction term, and the effective-coupling model for a general Simp. In a general case, the Hamiltonian (1) with a magnetic impurity in Majorana fermion representation reads:

$$H=frac{J}{4}sumlimits_{mathop {j,knotin {{Lambda }},}limits_{{langle jkrangle }_{mu }} },{u}_{jk}^{mu },i{c}_{j}{c}_{k}+frac{g}{4}sumlimits_{mathop {jin {{Lambda }},knotin {{Lambda }},}limits_{{langle jkrangle }_{mu }} }mathop{sum }limits_{a=1}^{2tilde{S}}(i{gamma }_{aj}^{mu }{b}_{k}^{mu })(i{gamma }_{aj}^{0}{c}_{k}),$$
(10)

where Λ = 0 denotes the impurity site in the triple-plaquette cluster (see, again, the inset of Fig. 1a). Here, we introduce bμ and c Majorana operators for the bulk spin-1/2 operator and γ-Majorana operators for the impurity spin as used in ref. 38. The Z2 gauge field on all the 12 edges is conserved and has eigenvalues ({u}_{jk}^{mu }=pm 1). The triple-plaquette flux, wI, is thus determined by ({w}_{I}=prod {u}_{jk}^{mu }). The second term in Eq. (10) represents the transformed impurity coupling, which becomes a four-Majorana interaction. Since (i{gamma }_{aj}^{mu }{b}_{k}^{mu }) does not commute with H, this quartic term can not be trivially rewritten as a quadratic form. In our ED calculation, this four-Majorana interaction is treated as it is by preparing 2M Majorana fermion operators composed of (M(in {mathbb{N}})) complex fermion operators in a binary number basis.

We visualize in Fig. 3a the triple-plaquette cluster with Simp = 3/2 in the Majorana representation. In this case, there are 27 Majorana fermions in the whole 13-spin-site system, including 12 from c-Majoranas of spin-1/2s, 12 from γ-Majoranas of the impurity site, and 3 from b-Majoranas at nearest-neighbor sites of the impurity.

Fig. 3: Fermionic representation and analysis for the triple-plaquette cluster.
figure 3

a Majorana representation for the spin-3/2 impurity case in the triple-plaquette cluster, following the parton construction by Ma38. The remaining sites are spin-1/2 degrees of freedom, described using Kitaev’s original parton construction12. b Energy difference, ΔE, between the bound-flux and zero-flux sectors for the Simp = 3/2 case. Magenta (orange) curve: ΔE computed with (without) projection operators. Green line: wI calculated in the spin Hamiltonian (1). c ΔE calculated in the effective-coupling model, Eq (18). The coefficients are given in the main text. From Simp = 1/2 to higher impurity spins, the substitution (Delta to sqrt{2{S}_{{rm{imp}}}}Delta) is applied.

Full size image

Using ED for the many-body Majorana Hamiltonian (10) under projections for the spin-3/2 impurity, we calculate the energy difference between the bound-flux (wI = − 1) and zero-flux (wI = 1) states of the triple-plaquette cluster. This energy difference, defined as

$$Delta Eequiv {E}_{{rm{bound}}}-{E}_{{rm{zero}}},$$
(11)

is evaluated as a function of g/J (see the magenta curve in Fig. 3b). When ΔE < 0 (ΔE > 0), the system realizes the bound-flux (zero-flux) sector at a given g/J, which is illustrated as a yellow (gray) shaded area in the figure. We observe two key behaviors: (i) the initial transition from the bound-flux sector to the zero-flux sector at the transition point g1, and (ii) the reentrant transition back to the bound-flux sector at the second transition point g2. For g > g2, ΔE decreases monotonically, indicating the stability of the reentrant bound-flux sector. Moreover, these two transition points match exactly with the flux-sector transition points obtained in ED of the spin-basis Hamiltonian. The calculated wI in spin-basis ED is shown as the green line in the Fig. 3b. It is also worth mentioning here that the correct behavior of flux-sector transition calculated in the Majorana representation is only obtained under the presence of projection operators Pa and ({P}_{tilde{S} = 3/2}) (see Methods for details). Without these projection operators, as shown by the orange curve in Fig. 3b), the reentrant bound-flux sector does not appear, emphasizing the need for proper constraints on the four-body Majorana terms.

The presented numerical results confirm that even a minimal cluster with a Simp = 3/2 impurity can qualitatively capture the flux-sector transitions and the reentrant effect, provided a proper treatment of projection. It reinforces the localized picture of the flux-binding effect by a site defect.

To further extract the essential ingredient that leads to a reentrant bound-flux sector, we propose a heuristic approach based on the effective coupling between the c-Majorana eigenmodes of the 12-site plaquette and the Majorana zero modes (({gamma }_{a}^{0})) introduced by the impurity. In the first constituent, the L = 12 plaquette is considered as a fermion hopping problem on a ring, where L denotes the number of sites on the ring44. The difference between the zero- and π-flux sectors is translated into the difference between the periodic and antiperiodic boundary conditions (PBC and APBC) of the ring (see Fig. 4a). Here, the π-flux sector corresponds to the bound-flux sector we have discussed. If the nearest-neighbor hopping strength is J/4 (corresponds to the Kitaev coupling in the spin Hamiltonian), the free-fermion Hamiltonian on the ring reads:

$${{mathcal{H}}}_{{rm{PBC/APBC}}}^{{rm{ring}}}=frac{J}{4}mathop{sum }limits_{j=1}^{L-1}left({a}_{j}^{dagger }{a}_{j+1}+h.,c.right)pm frac{J}{4}left({a}_{L}^{dagger }{a}_{1}+h.,c.right),$$
(12)

with the corresponding energy eigenvalues

$$left{begin{array}{l}{epsilon }_{{rm{PBC}}}=-frac{J}{2}cos left(frac{2pi n}{L}right),quad \ {epsilon }_{{rm{APBC}}}=-frac{J}{2}cos left[frac{(2n+1)pi }{L}right],quad end{array}right.$$
(13)

where n is the integers from 0 to L − 1. Since the ground-state energy of the system is calculated by the sum of all negative-energy modes, it is straightforward to conclude that L = 4n + 2 favors the zero-flux sector (periodic boundary) and L = 4n favors the π-flux sector (antiperiodic boundary). This is consistent with the prediction by Lieb’s theorem of flux configuration45, even only a single plaquette being considered here46. Note that in this L = 4n case, ϵPBC contains two zero modes while ϵAPBC contains no zero modes under zero magnetic fields. This is the key difference that helps us understand the flux-sector transition in the presence of impurity.

Fig. 4: Eigenstates of the L = 12 fermionic ring.
figure 4

a The energy spectrum of the fermionic ring with periodic or antiperiodic boundary. The former corresponds to the zero-flux sector of WI, while the latter corresponds to the bound-flux sector. b The eigenfunctions for some selected eigenmodes. The sites in red can directly hop to the additional impurity fermion. However, if the sum over amplitudes on the red sites vanishes, the effective coupling will be zero.

Full size image

To tackle the quasivacancy or impurity problem with Simp = 1/2, we add an additional fermion that only directly couples to three sites of the L = 12 ring with effective coupling strength Δ. This Δ can vary for different eigenmodes of the ring coupled to the impurity fermions, but in general, it is proportional to g. Therefore, the effective-coupling Hamiltonian is a tight-binding matrix between some of the eigenmodes of the ring and the impurity fermions. By the symmetry of the wavefunction, only one zero mode (ϵ = 0) and one particle-hole pair (ϵ = ± αJ) of the periodic ring (i.e., zero-flux sector) can hybridize with the impurity (see Methods). This leads to a simple tight-binding matrix:

$$left.begin{array}{cc}&{{mathcal{H}}}_{{rm{PBC}}}^{{rm{eff}}}=left(begin{array}{cccc}0&0&0&{Delta }_{1}\ 0&-alpha J&0&{Delta }_{2}\ 0&0&alpha J&{Delta }_{2}\ {Delta }_{1}&{Delta }_{2}&{Delta }_{2}&0end{array}right)end{array}right.,$$
(14)

where Δi ~ g is the effective coupling strength and we assume Δi ≡ Δ for simplicity. This Hamiltonian gives the eigenvalues of the effective-coupling model when Δ 1:

$${epsilon }_{{rm{PBC}}}^{{prime} }approx pm Delta ,,,,pm left(alpha J-frac{{Delta }^{2}}{alpha J}right),$$
(15)

where the prime is added to effective-coupling model eigenvalues ({epsilon }^{{prime} }) in order to be distinguished from the hopping-ring model eigenvalues ϵ. One can simply understand the above results by the perturbation theory. If the eigenmode of the L = 12 ring is a zero-energy mode, the effective coupling results in a degenerate perturbation theory with energy correction Δ. On the other hand, if the eigenmodes of the ring have finite energy, the non-degenerate perturbation theory gives rise to second-order corrections. Therefore, we can conclude that the ground-state total energy, which is the sum of all negative eigenvalues, has the general expressions (in the unit of J):

$${E}_{{rm{PBC}}}^{{rm{eff}}}approx -{A}_{1}Delta -{B}_{1}{Delta }^{2}-{C}_{1},$$
(16)

where A1, B1, and C1 are positive constants. In the bound-flux case, there is no zero-energy eigenmode in the plaquette model, so the ground-state energy expression does not include the linear term in Δ when Δ 1:

$${E}_{{rm{APBC}}}^{{rm{eff}}}approx -{B}_{2}{Delta }^{2}-{C}_{2}.$$
(17)

In the effective-coupling model, the ground-state energy difference defined in Eq. (11) translates into (Delta E={E}_{{rm{APBC}}}^{{rm{eff}}}-{E}_{{rm{PBC}}}^{{rm{eff}}}), which is the measure of the flux-sector transition. It can be modeled as

$$Delta Eapprox {A}_{1}Delta +({B}_{1}-{B}_{2}){Delta }^{2}+({C}_{1}-{C}_{2}),$$
(18)

and the coefficients are fitted by the exact-diagonalization result for the 13-site fermion-hopping model, which gives A1 ≈ 1.163, (B1B2) ≈ − 1.079, and (C1C2) ≈ − 0.261. Apparently, one can see that when Δ 1, ΔE is a concave quadratic function and provides the possibility of first (bound-to-zero) and second (zero-to-bound) flux-sector transitions. However, if the predicted second transition happens at the strength ({Delta }_{2}^{* }) beyond the validity of the quadratic approximation, it may not be seen in the exact diagonalization result. This is exactly what happens for the Simp = 1/2 case shown in the Fig. 3c, where the exact diagonalization curve starts to deviate from the quadratic curve at Δ ≈ 0.5J.

The next important question is how to incorporate the higher-spin impurity. Here, we conjecture that the essential ingredient is the additional ({gamma }_{a}^{0}) Majorana fermions that provide more entries of the tight-binding matrix. For example, for Simp = 3/2, one introduces three additional zero modes and considers the perturbative effects on the zero mode (ϵ0 = 0) and the finite-energy modes (ϵα = ± αJ) separately:

$$begin{array}{lll},{{mathcal{H}}}_{0}=left(begin{array}{cccc}0,&Delta ,&Delta ,&Delta \ Delta ,&0,&0,&0\ Delta ,&0,&0,&0\ Delta ,&0,&0,&0end{array}right),,,{epsilon }_{0}^{{prime} }=0,pm sqrt{3}Delta \\ ,{{mathcal{H}}}_{alpha }=left(begin{array}{ccccc}-alpha J,&0,&Delta ,&Delta ,&Delta \ 0,&alpha J,&Delta ,&Delta ,&Delta \ Delta ,&Delta ,&0,&0,&0\ Delta ,&Delta ,&0,&0,&0\ Delta ,&Delta ,&0,&0,&0end{array}right),,,{epsilon }_{alpha }^{{prime} }approx 0,pm left[alpha J+frac{3{Delta }^{2}}{alpha J}+{mathcal{O}}left(frac{{Delta }^{4}}{{(alpha J)}^{3}}right)right]end{array}$$
(19)

This indicates that the higher-spin effect can be incorporated by making the substitution (Delta to sqrt{2{S}_{{rm{imp}}}}Delta). In Fig. 3c, we show that the higher-spin impurity in the effective-coupling model simply shifts the quadratic curve to the left. This implies that for the higher-spin impurity case, the first transition point ({Delta }_{1}^{* }) is smaller, and the second transition ({Delta }_{2}^{* }) is more likely to happen when Δ 1 is still valid. This observation from the effective-coupling model, despite its over-simplicity, is consistent with the DMRG results as shown in Fig. 2a. Based on our findings, one may expect that higher-spin impurity cases such as spin-2, 5/2, and so on, also tend to bind the WI flux under the (Delta to sqrt{2{S}_{{rm{imp}}}}Delta) substitution. However, it is important to note that our effective model does not account for the effect of ({gamma }_{a}^{mu }), which may lead to unexpected results for even larger spins. Therefore, the precise nature of even higher-spin cases requires further numerical evidence to support the stability of the bound-flux sector.

Discussion

In this study, we investigated the behavior of Simp = 1 and Simp = 3/2 impurities in KSL, focusing on their effects on the flux-sector transitions and the stability of the bound-flux sector. Our analysis, using ED and DMRG methods, along with the phenomenological model in the Majorana representation, revealed several key findings: First, the local behavior of KSL around a magnetic impurity strongly depends on whether the impurity has a half-integer or integer spin. This dependence was demonstrated by considering the cases of Simp = 1 and Simp = 3/2, and analyzing their impact on flux-sector transitions. Second, magnetic impurities can bind Z2 fluxes in the lattice, similar to vacancies and quasivacancies. We observed a phase transition between bound-flux and zero-flux sectors by varying the impurity coupling strength, with the transition point dependent on the impurity’s spin magnitude. Third, for Simp = 3/2 impurities, a reentrant bound-flux sector was observed, remaining stable under finite magnetic fields. This stability increases with the coupling strength g. Fourth, our ED calculations for fermionic Hamiltonians in the triple-plaquette cluster revealed the energy differences between the bound-flux and zero-flux states, which are in good agreement with the flux-sector transition points identified in the DMRG study of the spin-basis Hamiltonian. Proper constraints by projection operators are essential for accurate flux-sector transition behavior.

The ability of spin-S impurities to bind Z2 fluxes and its stability against the external magnetic field has practical implications. When time-reversal symmetry is broken, flux binding to the impurity site results in the formation of localized Majorana zero modes, which are essential for realizing Ising anyons. These anyons exhibit non-Abelian statistics, making them valuable for topological quantum computation.

Compared to an Ising anyon bound at a vacancy13, the anyon found at the magnetic impurity site offers a more advantageous way to access low-energy Majorana-bound states through its magnetic channel. This unique feature of entangled Ising anyons can be observed in the dynamical correlation function of impurity spins, especially by focusing on the low-energy spectra.

Spectroscopy techniques with high spatial resolution, such as STM or nitrogen-vacancy (NV) center magnetometry, could further advance our understanding of the ground-state flux sector in Kitaev spin liquid phases. The low-energy spectrum carries crucial information about the ground-state flux sector, and both STM and NV center magnetometry are capable of probing these spectral features. STM might reveal the dynamical correlation function of impurity spins through spin-dependent tunneling21,23,24,30, while NV center magnetometry could detect magnetic noise linked to the same correlation function47 as well as emergent gauge field48. Observations of low-energy features, including Majorana zero modes, could potentially inform efforts toward realizing topological qubits composed of Ising anyons.

In addition, the even-odd effect on magnetization is another potential signature worth investigating. Local measurements of spin moments, such as those achievable via NV center magnetometry in the absence of a magnetic field, could provide clear evidence of this effect for a magnetic impurity embedded in Kitaev spin liquid phases. However, the spin Hamiltonian discussed in this work may be too simplified to accurately capture magnetic properties in real candidate materials11. Furthermore, introducing more than two magnetic impurities into the system could result in an RKKY-like interaction mediated by itinerant Majorana fermions in the bulk37,42,49. While this interaction could obscure or destabilize the even-odd effect, it may simultaneously promote magnetic ordering of impurity magnetic moments, even within an otherwise spin-liquid phase.

Methods

Implementation of DMRG

All DMRG calculations were performed using the NVIDIA Data Center GPU R470 Driver with the ITensorsGPU.jl package42. To ensure qualitative accuracy, determined by the expectation value of all plaquettes for both Wp in the bulk and WI (plus, even ({{mathbb{W}}}_{a}) for a half-integer impurity case) at impurity sites, we required an adequate bond dimension d depending on the impurity size and a good energy tolerance δE ≤ 1 × 10−7 while satisfying the cutoff at each sweep to be ≤ 1 × 10−9. For instance, in the single spin-3/2 impurity problem, we set the maximum bond dimension ({d}_{max }) to 3000 to ensure the cutoff condition. Additionally, we performed DMRG calculations five times independently for each parameter point and selected the result with the lowest ground state energy realizing reasonable flux expectation values.

Construction of projection operators for higher-spin

Here, we first review the Majorana representation for arbitrary spin size introduced by Ma38, and then construct projection operators applied in the ED calculation of the many-body Majorana Hamiltonian for the Simp = 3/2 case. The starting point is to consider a spin-(tilde{S}) operator as (2tilde{S}) of spin-1/2s:

$${tilde{S}}^{mu }=mathop{sum }limits_{a=1}^{2tilde{S}}{S}_{a}^{mu }=mathop{sum }limits_{a=1}^{2tilde{S}}frac{{sigma }_{a}^{mu }}{2},$$
(20)

where Sμ denotes the μ-component of spin-1/2 operator with σ being Pauli matrices, and a indicates the flavor degree of freedom. Then, the Majorana representation is applied for each spin-1/2 as ({sigma }_{a}^{mu }=i{gamma }_{a}^{mu }{gamma }_{a}^{0}) (The order of Majorana operators is slightly modified from the original to follow Kitaev’s Majorana representation12.). Here four Majorana fermions γx, γy, γz and γ0 are introduced for each flavor. As a result,

$${tilde{S}}^{mu }=frac{1}{2}mathop{sum }limits_{a=1}^{2tilde{S}}i{gamma }_{a}^{mu }{gamma }_{a}^{0}.$$
(21)

Note that for spin-1/2 operators, we use the usual symbol of Majorana fermions bx, by, bz and c instead of ({gamma }_{a}^{eta }) (η = x, y, z, 0). The local Z2 gauge field operators ({u}_{jk}^{mu }=i{b}_{j}^{mu }{b}_{k}^{mu }), which connect two spin-1/2s at sites j and k on the μ-bond, commute with the Hamiltonian H.

Since the Hilbert space is expanded in this representation- both for the bulk spin-1/2s and impurity sites-we need two kinds of projection operators to ensure the correct physical states are selected. The first one is required to ensure the commutation relations of Pauli matrices at each spin-1/2 or flavor in the Majorana representation, denoted by the condition ({C}_{gamma }^{a}:{D}_{a}={gamma }_{a}^{x}{gamma }_{a}^{y}{gamma }_{a}^{z}{gamma }_{a}^{0}=1) in ref. 38. Resulting operator is

$${P}_{a}=frac{1+{D}_{a}}{2}qquad (a=1,2,cdots ,,2tilde{S}).$$
(22)

Note that Pa satisfies the condition of a projection operator automatically as ({({P}_{a})}^{2}={P}_{a}).

The second projection operator, which is required only for higher-spin ((tilde{S} ,>, 1/2)) cases, mixes different flavors to ensure (| {tilde{S}}^{2}| =tilde{S}(tilde{S}+1)). This condition is represented as ({C}_{s}:{sum }_{mu }{(mathop{sum }nolimits_{a = 1}^{2tilde{S}}{gamma }_{a}^{mu }{gamma }_{a}^{0})}^{2}=-4tilde{S}(tilde{S}+1)) in ref. 38. For (tilde{S}=3/2), corresponding operator can be represented as

$${P}_{tilde{S} = 3/2}=-frac{1}{6}left[-3+sum _{mu }sum _{a > b}{gamma }_{a}^{mu }{gamma }_{a}^{0}{gamma }_{b}^{mu }{gamma }_{b}^{0}right].$$
(23)

It is worth mentioning that Pa is necessary for ({P}_{tilde{S} = 3/2}) to satisfy the condition of projection operator since

$${({P}_{tilde{S} = 3/2})}^{2}=-frac{1}{6}left[-3+sum _{mu }sum _{a > b},frac{4-{D}_{a}{D}_{b}}{3},{gamma }_{a}^{mu }{gamma }_{a}^{0}{gamma }_{b}^{mu }{gamma }_{b}^{0}right].$$
(24)

In the Majorana Hamiltonian with four-body interactions in the vicinity of the spin-3/2 impurity, these two kinds of projection operators ({P}_{tilde{S} = 3/2}) and Pa for a = 1, 2, 3 are essential for accurately calculating the energy difference between the two flux sectors.

Eigenenergies and eigenfunctions of the periodic ring

Here we provide details of the L = 12 periodic and antiperiodic rings used in the effective-coupling model. The energy-level spectrum calculated by Eq. (13) is shown in Fig. 4a. This spectrum can be verified by the diagonalization of the tight-binding model as well. The periodic ring, which corresponds to the zero-flux sector, contains two zero-energy modes. This is the crucial difference from the antiperiodic ring, which corresponds to the bound-flux sector. It is important to note that not all eigenmodes can couple to the impurity. First, the impurity fermion only directly hops to three nearest-neighbor sites, which are colored in red in Fig. 4b. Second, one can straightforwardly symmetrize the eigenmodes for each degenerate pair of modes, such that only some of the modes couple to the impurity fermion. Specifically, the 12-site model on a periodic ring has a C3v point-group symmetry, such that only three modes are subject to the transformation of the trivial irreducible representation and couple to the impurity site. Because the hopping integral is based on the overlap of one site (impurity) and the three neighboring sites simultaneously, one can simply assume ψimp = 1 as the impurity wavefunction and the sum of the amplitude of the three sites determines the effective coupling. Specifically, for a given eigenmode of the fermionic ring ψring, if

$$Delta sim sum _{j}{psi }_{{rm{imp}}},g,{psi }_{{rm{ring}},j}=sum _{j}g,{psi }_{{rm{ring}},j}=0$$
(25)

with j (3, 7, 11) in Fig. 4b, the eigenmode is simply decoupled to the impurity site. From Fig. 4b, we see that only the eigenmodes with ϵ = ± J/2 and one of the zero modes ϵ = 0 can couple to the impurity. This is numerically verified by the exact diagonalization of the 13-site tight-binding matrix. A similar analysis can be made on the antiperiodic ring for the bound-flux sector. However, since there is no zero-energy mode in the spectrum, the perturbative effect on the eigenmodes by the impurity fermion starts from the second-order correction. This leads to the qualitative argument based on our effective-coupling model, which is discussed in the main text.

Related Articles

Photonic topological insulators in femtosecond laser direct-written waveguides

Topological photonics attract significant interests due to their intriguing fundamental physics and potential applications. Researchers are actively exploring various artificial platforms to realize novel topological phenomena, which provides promising pathways for the development of robust photonic devices. Among these platforms, femtosecond laser direct-written photonic waveguides show unique ability to visualize intricate light dynamics in 2 + 1 dimensions, which rendering them ideal tools for investigating topological photonics. By integrating topological concepts into these waveguides, researchers not only deepen their understanding of topological physics but also provide potential methodology for developing advanced topological photonic integrated devices. In this review, we discuss recent experimental implementations of different topological phases within femtosecond laser direct-written photonic waveguides, as well as the fascinating physical phenomena induced by the interplay of topology with non-Hermiticity, nonlinearity and quantum physics are also introduced. The exploration of topological waveguide arrays shows great promise in advancing the field of topological photonics, providing a solid foundation for further research and innovation in this rapidly developing domain.

Quantum topological photonics with special focus on waveguide systems

In the burgeoning field of quantum topological photonics, waveguide systems play a crucial role. This perspective delves into the intricate interplay between photonic waveguides and topological phenomena, underscoring the theoretical underpinnings of topological insulators and their photonic manifestations. We highlight key milestones and breakthroughs in topological photonics using waveguide systems, alongside an in-depth analysis of their fabrication techniques and tunability. The discussion includes the technological advancements and challenges, limitations of current methods, and potential strategies for improvement. This perspective also examines the quantum states of light in topological waveguides, where the confluence of topology and quantum optics promises robust avenues for quantum communication and computing. Concluding with a forward-looking view, we aim to inspire new research and innovation in quantum topological photonics, highlighting its potential for the next generation of photonic technologies.

Topological excitations at time vortices in periodically driven systems

We consider two-dimensional periodically driven systems of fermions with particle-hole symmetry. Such systems support non-trivial topological phases, including ones that cannot be realized in equilibrium. We show that a space-time defect in the driving Hamiltonian, dubbed a “time vortex,” can bind π Majorana modes. A time vortex is a point in space around which the phase lag of the Hamiltonian changes by a multiple of 2π. We demonstrate this behavior on a periodically driven version of Kitaev’s honeycomb spin model, where ({{mathbb{Z}}}_{2}) fluxes and time vortices can realize any combination of 0 and π Majorana modes. We show that a time vortex can be created using Clifford gates, simplifying its realization in near-term quantum simulators.

2.5-dimensional topological superconductivity in twisted superconducting flakes

Multilayer flakes of two-dimensional materials were recently shown to be tunable by twisting monolayers on their surface. This raises the question whether qualitatively new phenomena can occur in such finite-thickness moiré systems. Here we demonstrate the emergence of distinct topological phases and transitions in N-layered flakes of nodal superconductors with a single monolayer twisted on top of it. We show that a c-axis current transforms the whole system into a chiral topological superconductor. Increasing the current drives a sequence of topological transitions between states characterized by a Chern number increasing from (sim {mathcal{O}}(N)) up to (sim {mathcal{O}}({N}^{2})), well beyond the additive effect of stacking N layers. We predict thickness-independent signatures of these states in the thermal Hall and tunneling microscopy measurements. Twisted superconductor flakes thus provide an example of a “2.5-dimensional” material where the synergy of two-dimensional layers extended in a third dimension realize states inaccessible in either monolayer or bulk materials.

Coulomb interactions and migrating Dirac cones imaged by local quantum oscillations in twisted graphene

Flat-band moiré graphene systems are a quintessential platform for investigating correlated phases of matter. Various interaction-driven ground states have been proposed, but despite extensive experimental effort, there has been little direct evidence that distinguishes between various phases, in particular near the charge neutrality point. Here we probe the fine details of the density of states and the effects of Coulomb interactions in alternating-twist trilayer graphene by imaging the local thermodynamic quantum oscillations with a nanoscale scanning superconducting quantum interference device. We find that the charging self-energy due to occupied electronic states is most important in explaining the high-carrier-density physics. At half-filling of the conduction flat band, we observe ferromagnetic-driven symmetry breaking, suggesting that it is the most robust mechanism in the hierarchy of phase transitions. Near charge neutrality, where exchange energy dominates over charging self-energy, we find a nematic semimetal ground state, which is theoretically favoured over gapped states in the presence of heterostrain. In this semimetallic phase, the flat-band Dirac cones migrate towards the mini-Brillouin zone centre, spontaneously breaking the threefold rotational symmetry. Our low-field local quantum oscillation technique can be used to explore the ground states of many strongly interacting van der Waals systems.

Responses

Your email address will not be published. Required fields are marked *