Collective spin dynamics of Z2 vortex crystals in triangular Kitaev-Heisenberg antiferromagnets

TL;DR

This paper identifies distinctive dynamical spin response signatures, including pseudo-Goldstone modes and multi-scattered intensity profiles, as key indicators of Z2 vortex crystals in layered triangular Kitaev-Heisenberg antiferromagnets.

Contribution

It reveals how to experimentally detect Z2 vortex crystals through their unique dynamical spin response signatures.

Findings

Presence of low-frequency pseudo-Goldstone phonon modes

Multi-scattered intensity profile at higher frequencies

Signatures of large vortex sizes and solitonic spin profiles

Abstract

We show that the mesoscopic incommensurate $\mathbb{Z}_2$ vortex crystals proposed for layered triangular anisotropic magnets can be most saliently identified by two distinctive signatures in dynamical spin response experiments: The presence of pseudo-Goldstone `phonon' modes at low frequencies $\omega$, associated with the collective vibrations of the vortex cores, and a characteristic multi-scattered intensity profile at higher $\omega$, arising from a large number of Bragg reflections and magnon bandgaps. These are direct fingerprints of the large vortex sizes and magnetic unit cells and the solitonic spin profile around the vortex cores.

Figures

Tables

Table 1. Table 1: The two types of magnetic unit cells depending on the period d 𝑑 d and the sign of K 𝐾 K . Spanning vector 𝐓 1 subscript 𝐓 1 {\bf T}_{1} and 𝐓 2 subscript 𝐓 2 {\bf T}_{2} , their mapping from one sublattice to another, and number 𝒩 m subscript 𝒩 𝑚 \mathcal{N}_{m} of spins in each magnetic unit cell.

sgn($K$)	$\begin{array}{c}\text{Condition}\\ \text{on}d\end{array}$	$\begin{array}{c}\text{Sublattice}\\ \text{mapping}\end{array}$	${𝐓}_1$	${𝐓}_2$	${𝒩}_m$
	$mod(d,3)=1$	(A,B,C) $\to$ (C,A,B)	$2d𝐚$	$2d𝐛$	$4d^2$
	otherwise	(A,B,C) $\to$ (A,B,C)	$2d\left(𝐚-𝐜\right)$	$2d\left(𝐛-𝐜\right)$	$12d^2$
$K>0$	$mod(d,3)=2$	(A,B,C) $\to$ (B,C,A)	$2d𝐚$	$2d𝐛$	$4d^2$
	otherwise	(A,B,C) $\to$ (A,B,C)	$2d\left(𝐚-𝐜\right)$	$2d\left(𝐛-𝐜\right)$	$12d^2$

Equations

$$\mathcal{H}={\sum }_{⟨ij⟩}\left(J{\mathbf{S}}_i\cdot {\mathbf{S}}_j+K{S}_i^{{\gamma }_{ij}}{S}_j^{{\gamma }_{ij}}\right).$$

$${\gamma }_{ij}=x,y,\text{or}z,$$

$${\mathcal{S}}^{\alpha \beta }(\mathbf{Q},\omega )=\int \text{dt}{e}^{-i\omega t}⟨S^{\alpha }(-\mathbf{Q},0)S^{\beta }(\mathbf{Q},t)⟩,$$

$$\mathcal{I}(\mathbf{Q},\omega )\propto \sum _{\alpha \beta }({\delta }_{\alpha \beta }-Q^{\alpha }{Q}^{\beta }/Q^2){\mathcal{S}}^{\alpha \beta }(\mathbf{Q},\omega ),$$

$${\mathbf{S}}_i\to {\mathbf{S}}_{\mathbf{R},\mu }\text{and}{\mathbf{r}}_i=\mathbf{R}+{\mathbf{\rho }}_{\mu },$$

$$\mathcal{H}=\frac{1}{2}\sum _{\mathbf{R}}\sum _{\mu =1}^{{\mathcal{N}}_m}\sum _{\nu =1}^{{\mathcal{N}}_m}{\mathbf{S}}_{\mathbf{R},\mu }^{\text{T}}\cdot {\mathcal{J}}_{\mu \nu }\cdot {\mathbf{S}}_{\mathbf{R}+{\mathbf{t}}_{\mu \nu },\nu },$$

$${\mathbf{S}}_{\mathbf{R},\mu }^{\text{T}}=(S_{\mathbf{R},\mu }^x,S_{\mathbf{R},\mu }^y,S_{\mathbf{R},\mu }^z),$$

$$\bm{\mathcal{J}}_{\mu\nu}=\left\{\!\!\begin{array}[]{rl}\bm{\mathcal{J}}_{\bm{\epsilon}},&\text{if}~{}\left(\mathbf{R}+\bm{\rho}_{\mu}\right)-\left(\mathbf{R}+{\bf t}_{\mu\nu}+\bm{\rho}_{\nu}\right)=\pm\bm{\epsilon},\\ 0,&{\text{otherwise}}\end{array}\right.\,,$$

$${\mathcal{J}}_{\mathbf{a}}=\left(\begin{array}{ccc}J& 0& 0\\ 0& J+K& 0\\ 0& 0& J\end{array}\right),{\mathcal{J}}_{\mathbf{b}}=\left(\begin{array}{ccc}J& 0& 0\\ 0& J& 0\\ 0& 0& J+K\end{array}\right),{\mathcal{J}}_{\mathbf{c}}=\left(\begin{array}{ccc}J+K& 0& 0\\ 0& J& 0\\ 0& 0& J\end{array}\right),$$

$$\{{\tilde{\mathbf{x}}}_i,{\tilde{\mathbf{y}}}_i,{\tilde{\mathbf{z}}}_i\}$$

$${\tilde{\mathbf{S}}}_{\mathbf{R},\mu }={\mathbf{U}}_{\mu }\cdot {\mathbf{S}}_{\mathbf{R},\mu }.$$

$${\mathbf{U}}_{\mu }=\left(\begin{array}{ccc}\cos {\theta }_{\mu }\cos {\varphi }_{\mu }& \cos {\theta }_{\mu }\sin {\varphi }_{\mu }& -\sin {\theta }_{\mu }\\ -\sin {\varphi }_{\mu }& \cos {\varphi }_{\mu }& 0\\ \sin {\theta }_{\mu }\cos {\varphi }_{\mu }& \sin {\theta }_{\mu }\sin {\varphi }_{\mu }& \cos {\theta }_{\mu }\end{array}\right).$$

$$\mathcal{H}=\sum _{\mathbf{R}}\sum _{\mu \nu }{\tilde{\mathbf{S}}}_{\mathbf{R},\mu }^{\text{T}}\cdot {\mathcal{T}}_{\mu \nu }\cdot {\tilde{\mathbf{S}}}_{\mathbf{R}+{\mathbf{t}}_{\mu \nu },\nu },$$

$$\begin{array}[]{c}\widetilde{S}_{{\bf R},\mu}^{x}\!\!\approx\!\!\sqrt{\frac{S}{2}}(a_{{\bf R},\mu}+a_{{\bf R},\mu}^{\dagger})\,,~{}\widetilde{S}_{{\bf R},\mu}^{y}\!\!\approx\!\!-i\sqrt{\frac{S}{2}}(a_{{\bf R},\mu}-a_{{\bf R},\mu}^{\dagger})\,,\\ \\ \widetilde{S}_{{\bf R},\mu}^{z}=S-a_{{\bf R},\mu}^{\dagger}a_{{\bf R},\mu}\,.\end{array}$$

$$\mathcal{H}={\mathcal{H}}_0+{\mathcal{H}}_1+{\mathcal{H}}_2+\mathcal{O}(S^{1/2}),$$

$${\mathcal{H}}_0=S^2\sum _{\mathbf{R}}\sum _{\mu \nu }{\mathcal{T}}_{\mu \nu }^{(3,3)}$$

$$\mathcal{H}_{1}=\sqrt{\frac{1}{2}}S^{3/2}\sum_{\mathbf{R}}\sum_{\mu\nu}\Big{[}\left(\mathcal{T}_{\mu\nu}^{(1,3)}-i\mathcal{T}_{\mu\nu}^{(2,3)}\right)a_{{\bf R},\mu}+\text{h.c.}\Big{]}$$

$${\mathcal{H}}_2$$

$$\begin{array}[]{c}f_{\mu\nu}=\mathcal{T}_{\mu\nu}^{(1,1)}-i\mathcal{T}_{\mu\nu}^{(1,2)}-i\mathcal{T}_{\mu\nu}^{(2,1)}-\mathcal{T}_{\mu\nu}^{(2,2)},\\ g_{\mu\nu}=\mathcal{T}_{\mu\nu}^{(1,1)}+i\mathcal{T}_{\mu\nu}^{(1,2)}-i\mathcal{T}_{\mu\nu}^{(2,1)}+\mathcal{T}_{\mu\nu}^{(2,2)}\,.\end{array}$$

$$a_{\mathbf{R},\mu }=\frac{1}{\sqrt{{\mathcal{N}}_m}}\sum _{\mathbf{q}}{e}^{i\mathbf{q}\cdot (\mathbf{R}+{\mathbf{\rho }}_{\mu })}{a}_{\mu ,\mathbf{q}},$$

$${\mathcal{H}}_2=E_{cl}/S+\sum _{\mathbf{q}}\sum _{\mu \nu }\frac{S}{4}{\mathcal{H}}_{2,\mathbf{q},\mu \nu }$$

$$\displaystyle\mathcal{H}_{2,\mathbf{q},\mu\nu}=f_{\mu\nu}\Big{[}e^{i\mathbf{q}\cdot\bm{\delta}_{\mu\nu}}a_{\mu,\mathbf{q}}a_{\nu,-\mathbf{q}}+e^{-i\mathbf{q}\cdot\bm{\delta}_{\mu\nu}}a_{\mu,-\mathbf{q}}a_{\nu,\mathbf{q}}\Big{]}+\text{h.c.}$$

$$\displaystyle+g_{\mu\nu}\Big{[}e^{i\mathbf{q}\cdot\bm{\delta}_{\mu\nu}}a_{\nu,\mathbf{q}}^{\dagger}a_{\mu,\mathbf{q}}+e^{-i\mathbf{q}\cdot\bm{\delta}_{\mu\nu}}a_{\mu,-\mathbf{q}}a_{\nu,-\mathbf{q}}^{\dagger}\Big{]}+\text{h.c.}$$

$$\displaystyle-2\mathcal{T}_{\mu\nu}^{(3,3)}\Big{[}a_{\mu,-\mathbf{q}}a_{\mu,-\mathbf{q}}^{\dagger}+a_{\mu,\mathbf{q}}^{\dagger}a_{\mu,\mathbf{q}}+a_{\nu,-\mathbf{q}}a_{\nu,-\mathbf{q}}^{\dagger}+a_{\nu,\mathbf{q}}^{\dagger}a_{\nu,\mathbf{q}}\Big{]},~{}~{}~{}~{}~{}~{}$$

$${\mathcal{H}}_2=E_{cl}/S+\sum _{\mathbf{q}}{\mathbf{x}}_{\mathbf{q}}^{†}\cdot {\mathbf{H}}_{\mathbf{q}}\cdot {\mathbf{x}}_{\mathbf{q}},$$

$${\mathbf{y}}_{\mathbf{q}}=(b_{1,\mathbf{q}},...,b_{{\mathcal{N}}_m,\mathbf{q}},b_{1,-\mathbf{q}}^{†},...,b_{{\mathcal{N}}_m,-\mathbf{q}}^{†}{)}^{\text{T}},$$

$$(\mathbf{g}\cdot {\mathbf{H}}_{\mathbf{q}})\cdot {\mathbf{T}}_{\mathbf{q}}={\mathbf{T}}_{\mathbf{q}}\cdot (\mathbf{g}\cdot {\mathbf{\Omega }}_{\mathbf{q}}),$$

$${\mathcal{S}}^{\alpha \beta }(\mathbf{Q},\omega )$$

$$S_{\mathbf{R},\mu }^{\alpha }=\sqrt{\frac{S}{2}}{\xi }_{\mu }^{\alpha }{a}_{\mathbf{R},\mu }+\sqrt{\frac{S}{2}}{\xi }_{\mu }^{{\alpha }^{\ast }}{a}_{\mathbf{R},\mu }^{†}+{\lambda }_{\mu }^{\alpha }(S-a_{\mathbf{R},\mu }^{†}{a}_{\mathbf{R},\mu }),$$

$$S_{\mu }^{\alpha }(-\mathbf{Q},0)=$$

$$=$$

$${\mathcal{S}}^{\alpha \beta }(\mathbf{Q},\omega )=$$

Full text

I Introduction

Recently a significant experimental and theoretical effort has been devoted to the understanding of correlated electron systems with 4d and 5d transition metal ions (like Ru3+ and Ir4+), characterized by effective $J_{\text{eff}}=1/2$ pseudospins, edge-sharing oxygen octahedra and tri-coordinated lattice geometries Jackeli and Khaliullin (2009); Rau et al. (2016); Hermanns et al. (2018); Winter et al. (2017a); Trebst (2017); Takagi et al. (2019). Owing to the strong interplay of spin-orbit coupling, crystal field and electronic correlations, these systems show a remarkable range of unconventional phases Rau et al. (2016); Trebst (2017); Hermanns et al. (2018); Winter et al. (2017a); Chaloupka et al. (2013); Rau et al. (2014); Sizyuk et al. (2014); Rousochatzakis et al. (2015); Rousochatzakis and Perkins (2017); Ducatman et al. (2018), including the renowned quantum spin liquids, possibly realized in $\alpha$-RuCl3 Plumb et al. (2014); Sears et al. (2015); Johnson et al. (2015); Majumder et al. (2015); Banerjee et al. (2016); Kasahara et al. (2018), the counter-rotating incommensurate spirals realized in the layered honeycomb $\alpha$-Li2IrO3 and its 3D analogues $\beta$- and $\gamma$-Li2IrO3 Biffin et al. (2014a, b); Modic et al. (2014); Takayama et al. (2015); Veiga et al. (2017); Williams et al. (2016); Breznay et al. (2017); Majumder et al. (2018), and a variety of complex multi-sublattice, single- and multi-${\bf Q}$ phases predicted under a magnetic field Janssen et al. (2016); Chern et al. (2017); Janssen and Vojta (2019).

The basic ingredient overarching the low-energy descriptions of such systems is the presence of bond-dependent anisotropic exchange, with the so-called Kitaev interactions Kitaev (2006); Jackeli and Khaliullin (2009) being the most prominent. As the bond dependence stems from spin-orbit coupling, such interactions are not limited to tri-coordinated lattices, but may also appear in other geometries, including the common frustrated geometries of the triangular, kagome, pyrochlore, and hyperkagome lattices Kimchi and Vishwanath (2014); Jackeli and Avella (2015); Rousochatzakis et al. (2016); Becker et al. (2015); Catuneanu et al. (2015); Li et al. (2015a); Shinjo et al. (2016); Kos and Punk (2017); Yao and Dong (2016, 2018); Kishimoto et al. (2018). In such lattices, the synergy of bond-dependent anisotropy and geometric frustration opens up the possibility for novel cooperative phases even when the anisotropy is not the dominant interaction, as in the above tri-coordinated systems.

Already the introduction of an infinitesimal Kitaev anisotropy $K$ in one of the simplest frustrated geometries, the triangular lattice [Fig. 1 (a)], highlights the prolificacy of this synergy Rousochatzakis et al. (2016): The three-sublattice 120∘ order of the triangular Heisenberg antiferromagnet (HAF) is immediately unstable under $K$, giving way to incommensurate crystals of $\mathbb{Z}_{2}$ vortices of mesoscopic size [Figs. 1 (c-f)], see also [Becker et al., 2015; Catuneanu et al., 2015; Shinjo et al., 2016; Trebst, 2017; Li et al., 2015a; Kos and Punk, 2017; Kishimoto et al., 2018]. Such vortices have been known Kawamura and Miyashita (1984) to be present in triangular HAFs as topological excitations, but here the bond-dependent anisotropy condenses them in the ground state via a commensurate-incommensurate (C-IC) nucleation mechanism Gennes (1975); Bak (1982); McMillan (1976); Schaub and Mukamel (1985). This is akin to the formation of magnetic domains Gennes (1975), Abrikosov vortices Abrikosov (1957); Essmann and Träuble (1967); Hess et al. (1989), blue phases in cholesteric liquid crystals Wright and Mermin (1989), skyrmions in chiral helimagnets  Bogdanov and Yablonskii (1989); Rößler et al. (2006); Mühlbauer et al. (2009); Tonomura et al. (2012); Yu et al. (2012, 2010); Seki et al. (2012); Janson et al. (2014), and other systems Suzuki (1983); Okubo et al. (2012); Leonov and Mostovoy (2015); Rosales et al. (2015); Hayami et al. (2016). Anisotropic antiferromagnets with hexagonal symmetry provide, therefore, a fertile ground for novel incommensurate phases with topological, particle-like properties.

While the prospect of realizing the $\mathbb{Z}_{2}$ vortex phase remains currently open (see Sec. VI below), here we explore the collective spin dynamics in this phase and demonstrate numerically how its presence can be most saliently observed in dynamical spectroscopic probes. To this end, we construct a large family of $\mathbb{Z}_{2}$ vortex crystals ($\mathbb{Z}_{2}$VC’s), for both positive and negative Kitaev anisotropy – with magnetic unit cell sizes extending up to $2028$ spin sites – and perform a semiclassical $1/S$ expansion to extract the magnon spectrum, the associated spin dynamical structure factors (DSF) $\mathcal{S}^{\alpha\beta}({\bf Q},\omega)$ (for all relative polarizations $\alpha,\beta=\{x,y,z\}$), and the corresponding inelastic neutron scattering (INS) intensity $\mathcal{I}({\bf Q},\omega)$.

The results close to the C-IC transition mirror two of the most distinctive features of the $\mathbb{Z}_{2}$VC phase, the large size of the vortices and their particle-like nature. Conceptually, both of these features stem from the C-IC nature of the transition from the ‘parent’ 120∘ state Rousochatzakis et al. (2016). The vortex size is large close to the transition because the vortices play the role of ‘discommensurations’ of the parent state, and their relative distance must diverge when we recover that state. This manifests in $\mathcal{I}({\bf Q},\omega)$ by a distinguished multi-fragmentation of the ‘parent’ magnon bands, arising from a high density of Bragg reflections.

The particle-like character of the vortices manifests at low frequencies $\omega$ via the presence of intense pseudo-Goldstone modes. These modes are associated with collective vibrations of the vortex cores around their equilibrium positions, and are thus analogous to phonons in crystals. Their appearance attests to the nonlinear character of the spin profile around the cores. As shown in [Rousochatzakis et al., 2016], the vortices arise by a special intertwining of three honeycomb superstructures of ferromagnetic (FM) domains [one for each sublattice of the ‘parent’ 120∘ phase, see Fig. 1 (d-f)], and this arrangement gives rise to abrupt, soliton-like modulations around the vortex cores. As demonstrated below [Fig. 2 (a)], the ground state energy landscape (as a function of the core positions) flattens significantly as we approach the C-IC transition, revealing a weak inter-particle potential at large distances. The pseudo-Goldstone modes (which are otherwise gapped out by the lattice cutoff) are thus a manifestation of the nonlinear spin profiles of the cores.

The paper is organized as follows. We begin with the definition of the model (Sec. II), a brief review of the main features of the $\mathbb{Z}_{2}$VC’s (Sec. III), and the iterative variational method used to obtain optimal crystals for given model parameters (Sec. IV). Our results for the collective spin dynamics and the corresponding predictions for the inelastic neutron scattering intensity are presented in Sec. V. A brief outlook is given in Sec. VI, while auxiliary information and technical details are relegated to Apps. A-D.

II Model

The Hamiltonian of the Heisenberg-Kitaev or $JK$-model Jackeli and Khaliullin (2009) on the triangular lattice reads

[TABLE]

Here $\langle ij\rangle$ denotes nearest neighbor lattice sites, ${\bf S}_{i}$ and ${\bf S}_{j}$ are pseudospins-1/2 degrees of freedom, and $J$ and $K$ denote the Heisenberg and Kitaev exchange parameter, respectively. The component $\gamma_{ij}$ is given by

[TABLE]

depending on whether $\langle ij\rangle$ belongs to the ‘xx’, ‘yy’ or ‘zz’ type of bonds, see Fig. 1 (a). The lattice plane is (111), and the vectors ${\bf a}$, ${\bf b}$ and ${\bf c}$ shown in Fig. 1 (a) point along ${\bf z}\!-\!{\bf x}$, ${\bf x}\!-\!{\bf y}$ and ${\bf y}\!-\!{\bf z}$, respectively. In what follows we use the parametrization $J\!=\!\cos\psi$ and $K\!=\!\sin\psi$ and restrict ourselves to the stability region $\psi\!\in\!(\tan^{-1}(-1/2)$, $\pi/2)$ of the vortex phase [shaded in Fig. 1 (b)] Rousochatzakis et al. (2016); Becker et al. (2015). We also set the lattice parameter $a\!=\!1$.

The ground state of the HAF point ($\psi\!=\!0$) of the $JK$ model (Eq. (1)) is the well-known 120∘, three-sublattice coplanar order Yafet and Kittel (1952), whose order parameter is that of a rigid rotator, i.e., SO(3). Classical analysis Rousochatzakis et al. (2016) shows that the 120∘ pattern is immediately unstable under an infinitesimal Kitaev interaction, giving way to a nontrivial, long-distance twisting of the SO(3) order parameter in both directions of the lattice plane, leading to localized $\mathbb{Z}_{2}$ vortices (see also Becker et al. (2015); Trebst (2017)). The cores of the vortices form a triangular superstructure whose period $d$ (the distance between the vortex cores) is determined by the competition between the Kitaev exchange $K$ and the Heisenberg exchange $J$. For small $|K|/J$, $d\propto J/|K|$, i.e., the distance between vortex cores goes to infinity at the HAF point, and the transition between the 120∘ order and the $\mathbb{Z}_{2}$VC phase is of the C-IC nucleation type Gennes (1975); Bak (1982); McMillan (1976); Schaub and Mukamel (1985).

III Main aspects of the $\mathbb{Z}_{2}$VC phase

Let us recall the main features of the $\mathbb{Z}_{2}$VC phase Rousochatzakis et al. (2016). First, the cores of the vortices are defects of the $120^{\circ}$ state, as they are associated with a finite FM canting and a reduced vector chirality. This means that the cores cost Heisenberg energy. However, the Kitaev energy around the cores is negative, which is why having cores with a given density is energetically favorable.

Second, the distance $d$ between the cores (Fig. 1 (c)) and the size of the vortices are infinite at $\psi\!=\!0$ (HAF point), and decreases monotonously as we depart from this point. The minimum values of $d$ are reached at the phase boundaries with the neighboring phases at $\tan\psi\!=\!-1/2$ ($d\!=\!1$) and $\psi\!=\!\pi/2$ ($d\!=\!2$). The magnetic unit cell contains $\mathcal{N}_{m}\!=\!4d^{2}$ or $\mathcal{N}_{m}\!=\!12d^{2}$ spins, depending, respectively, on whether the translation vectors of the state, $\mathbf{T}_{1}$ and $\mathbf{T}_{2}$ [see Fig. 1 (c)], map spins from one type of sublattice to another or not 111The second possibility has been overlooked in Ref. [Rousochatzakis et al., 2016]., see detailed discussion in Sec. IV and Table 1.

Third, the anatomy of the $\mathbb{Z}_{2}$VC can be best understood by visualizing separately the spins in the three sublattices of the 120∘ state, see Figs. 1 (d-f). In contrast to the 120∘ state, where all spins of a given sublattice are parallel to each other, forming a single FM domain of infinite size, here the spins of a given sublattice form a hexagonal superstructure of FM domains. The $\mathbb{Z}_{2}$ vortices then arise by the special way the three sublattice superstructures are intertwined with each other. In particular, the center of a FM domain in one sublattice (say $A$) coincides with vertices of the hexagonal superstructures in the other two sublattices ($B$ and $C$). Therefore, as we trace a closed loop around the center of a FM domain of $A$, the spins of $A$ remain roughly parallel along the loop, whereas the spins of $B$ and $C$ complete a $2\pi$-rotation, leading to a $\mathbb{Z}_{2}$ vortex, see bold arrows in Fig. 1 (e-f). The precise way the $2\pi$-rotation happens is related to the special role of the $[111]$ and $[100]$ directions, see color coding of Figs. 1 (d-f) and detailed discussion in Ref. [Rousochatzakis et al., 2016].

Finally, the $\mathbb{Z}_{2}$VC state preserves the discrete threefold rotation symmetry of the model. As we show below, this gives rise to three pairs of pseudo-Goldstone modes which are related to each other by threefold rotations. These modes track the first harmonic Bragg peaks in the static spin structure factor Rousochatzakis et al. (2016); Becker et al. (2015). Namely, they emanate from the corners of the BZ as we depart from the HAF point, and move towards the $\Gamma$ point for $K\!>\!0$, or the ${\bf M}$ points for $K\!<\!0$ [see Figs. 4-5 below].

IV Optimal crystals and variational minimization method

For each given $\psi$ inside the stability region of the $\mathbb{Z}_{2}$VC phase, the optimal value of $d$ can be obtained by the variational energy minimization scheme outlined in [Rousochatzakis et al., 2016]. In this approach, one exploits the fact that the $\mathbb{Z}_{2}$VC’s consist of three honeycomb superstructures of ferromagnetic domains, one for each of the three sublattices (A, B and C) of the HAF point. The majority of spins within each FM domain point along a specific direction in spin space, which happens to be one of the four $\langle 111\rangle$ axes. We therefore begin by constructing, for each given $\psi$, an initial state consisting of perfect FM domains (where all spins in the domain are parallel to each other and along the respective $\langle 111\rangle$ axis) with a size that corresponds to a fixed choice of $d$. Next, upon a random sampling, we sequentially rotate spins in the direction of their local mean fields. After a certain number of samplings, the system converges to a $\mathbb{Z}_{2}$VC and the corresponding energy per site $E_{0}(d)/N$ is extracted. This procedure is repeated for a series of different FM domain wall sizes, corresponding to different choices of $d$ (and always using appropriate clusters with periodic boundary conditions that accommodate the given superstructure). The resulting energies per site $E_{0}(d)/N$ are then plotted as a function of $d$ and one identifies the optimal crystal with the one associated with the minimum energy. Three examples were shown in Fig. 2 (a), for $\tan\psi\!=\!-0.34$, $-0.29$ and $-0.22$, for which the minimum energies per site are reached at $d\!=\!7$, $10$ and $19$, respectively. Following these steps we construct a large set of optimal crystals [see blue dots in Fig. 1 (b)], with $d$ extending from $1$ ($\tan\psi\!=\!-0.42$, $\mathcal{N}_{m}\!=\!4$ spins in the magnetic unit cell) to $d\!=\!19$ ($\tan\psi\!=\!-0.22$, $\mathcal{N}_{m}\!=\!1444$) for negative $K$, and from $d\!=\!2$ ($\tan\psi\!=\!4$, $\mathcal{N}_{m}\!=\!16$) to $d\!=\!13$ ($\tan\psi\!=\!0.46$, $\mathcal{N}_{m}\!=\!2028$) for positive $K$.

V Dynamical fingerprints of the $\mathbb{Z}_{2}$VC phase

The collective spin dynamics can now be studied, for each of these optimal crystals, using a numerical implementation of the Holstein-Primakoff transformation, followed by a generalized Bogoliubov transformation, and a numerical diagonalization that delivers the $\mathcal{N}_{m}$ magnon bands in the magnetic BZ. This is then used for the evaluation of

[TABLE]

where ${\bf S}(\mathbf{Q},t)$ is the Fourier transform of the total spin with ${\bf Q}$ in the first BZ of the lattice, and

[TABLE]

for further technical details see App. B.

V.1 Linear spin wave (LSW) expansion

Figure 3 shows the LSW dispersions for eight representative optimal $\mathbb{Z}_{2}$VC’s. The spectra are first obtained in the magnetic BZ and then plotted in the repeated scheme, along special symmetry directions in the lattice Brillouin zone (see hexagons in Fig. 6). Panel (d) shows the familiar result for the 120*∘*-order of the pure HAF ($\psi\!=0$), which can actually be considered as a $\mathbb{Z}_{2}$VC state with $d=\infty$ Mourigal et al. (2013). As we gradually move away from the HAF point, the size of the $\mathbb{Z}_{2}$-vortex becomes finite but still remains very large. Recall that the size of magnetic unit cell is $\mathcal{N}_{m}=4d^{2}$ or $12d^{2}$, depending on the orientation of the spanning vectors of the superlattice, see above. This explains the large number of $\mathcal{N}_{m}$ magnon bands that are visible in Fig. 3, except for panels (d) and (h). The figure also shows the band gaps between neighboring magnon bands, which result from Bragg reflections of the spin waves off the boundaries of the large magnetic unit cells. This high density of Bragg reflections and magnon band gaps is responsible for the multi-fragmented scattering profile announced above.

V.2 Spin dynamical structure factors (DSF)

Figures 4 (a-b) show the diagonal components ${\mathcal{S}}^{xx}(\mathbf{Q},\omega)$, ${\mathcal{S}}^{yy}(\mathbf{Q},\omega)$ and ${\mathcal{S}}^{zz}(\mathbf{Q},\omega)$ for two representative optimal crystals with large $d$, one at $\tan\psi\!=\!0.502$ ($d\!=\!11$, $\mathcal{N}_{m}\!=\!484$) and the other at $\tan\psi\!=\!-0.34$ ($d\!=\!7$, $\mathcal{N}_{m}\!=\!196$). First of all, it can be clearly seen that the three diagonal components are indeed related to each other by the threefold symmetry. Second, the overall shape of the DSF at intermediate and high $\omega$ (i.e., far enough from the corners of the lattice BZ) follows very roughly the shape of the three magnon bands of the DSF of the parent 120∘ order (see Fig. 3 (d) and top panels in Fig. 5, as well as Ref. Mourigal et al. (2013)). Equivalently, the unfolded (in the lattice BZ) $\mathcal{N}_{m}$ magnon bands of the $\mathbb{Z}_{2}$VC follow roughly the overall shape of the three ‘parent bands’. This is due to the fact that the magnon wavelengths in this part of the spectrum can be significantly smaller than the distance $d$ between vortices, and the short-distance fluctuations are still governed by the Heisenberg exchange. Despite this rough similarity, the huge number of Bragg reflections and associated band gaps (resulting from the large magnetic unit cell) give rise to a qualitative different DSF, with only a small portion of the bands standing out and an otherwise smeared out and multi-fragmented response.

The most intense modes in Figs. 4 (a-b) appear at low $\omega$, close to the corners of the BZ, where the magnon wavelengths become comparable to the distance $d$ between vortices. These intense modes are the collective, pseudo-Goldstone modes mentioned above, associated with the rigid vibrations of the vortex cores around their equilibrium positions. There are three ($\pm{\bf Q}$) pairs of such phonon-like modes [one for each diagonal component of $\mathcal{S}^{\alpha\alpha}({\bf Q},\omega)$], and their positions coincide with those of the first harmonics of the static structure factor Rousochatzakis et al. (2016); Becker et al. (2015), see yellow stars in the insets of Fig. 4. All in all, Fig. 4 therefore demonstrates the two most salient dynamical fingerprints of the $\mathbb{Z}_{2}$VC phase in the vicinity of the C-IC transition, the large vortex size and their particle-like character.

V.3 Evolution of the spectra with $K/J$

We now proceed to elucidate the way these features evolve as we move deeper into the $\mathbb{Z}_{2}$VC phase, and the vortices become smaller in size. To this end, we consider a series of eight representative $\mathbb{Z}_{2}$VC’s with decreasing $d$, four for $K\!>\!0$ [panels (a-b), $d\!=\!\{13,10,5,2\}$] and four for $K\!<\!0$ [panel (c), $d\!=\!\{19,9,2,1\}$]; For results on many more representative states see App. D. Figure 5 shows the associated $\mathcal{I}(\mathbf{Q},\omega)$, along with the intensity of the HAF point (top panels, $d\!=\!\infty$). The rough resemblance mentioned above, between the overall shape of the response with that of the parent state, persists down to $d\!=\!5$ and $\mathcal{N}_{m}\!=\!100$ for $K\!>\!0$, and down to $d\!=\!2$ and $\mathcal{N}_{m}\!=\!48$ for $K\!<\!0$. For smaller vortex sizes new features appear, such as the distinctively rich pattern for $\tan\psi\!=\!4$ (panel a, bottom) which is characteristic of strong Kitaev physics (see also below), and the two-band picture for $\tan\psi\!=\!-0.42$, which is characteristic of the neighboring $\widetilde{\text{FM}}$ state Rousochatzakis et al. (2016); Becker et al. (2015).

Turning to the evolution of the phonon-like modes, these must track the positions of the first harmonic Bragg peaks, as mentioned above. This is illustrated by the red dashed lines in panels (a) and (c). For $K\!>\!0$ (panel a), one of the phonon modes traces the path $\mathrm{K}\!\to\!\mathrm{O}_{1}\!\to\!\mathrm{O}_{2}\!\to\!\mathrm{O}_{3}\!\to\!\mathrm{O}_{4}$, while for $K\!<\!0$ (panel c) the phonon mode shown goes from $\mathrm{K}\to\mathrm{M}_{y}$, and similarly for the remaining phonon modes related by threefold rotations.

In addition to the phonon modes, we also find a second intense low-$\omega$ mode. For $K\!>\!0$, this is shown by the dashed yellow line in panel (a) and is elucidated further in panel (b). This mode traces the path $\mathrm{K}\!\to\!\mathrm{P}_{1}\!\to\!\mathrm{P}_{2}\!\to\!\mathrm{P}_{3}\!\to\!\mathrm{M}_{y}$, and is a precursor of an accidental, classical ground state degeneracy present at the Kitaev point ($\psi\!=\!\pi/2$) Rousochatzakis et al. (2016). This degeneracy is sub-extensive and manifests in the Fourier transform of the classical energy with lines of minima joining the ${\bf M}$ points of the BZ (e.g., the line $\mathrm{M}_{x}\to\mathrm{M}_{y}$). This is illustrated in panel (b) which shows the intensity along special horizontal cuts (panel b, top) parallel to $\mathrm{O}_{1}\mathrm{P}_{1}$ ($\tan\psi\!=\!0.46$), $\mathrm{O}_{2}\mathrm{P}_{2}$ ($\tan\psi\!=\!0.561$), $\mathrm{O}_{3}\mathrm{P}_{3}$ ($\tan\psi\!=\!0.93$) and $\mathrm{O}_{4}\to\mathrm{M}_{y}$ ($\tan\psi\!=\!4$). The intensities along these cuts show the development of an almost flat mode, which should ideally become completely flat at $\psi\!=\!\pi/2$ (Kitaev point). While quantum fluctuations eventually remove this degeneracy Jackeli and Avella (2015), the almost flat precursor of this physics away from the Kitaev point could still be observable.

V.4 Evolution of the spin gap with $K/J$

The presence of exchange anisotropy and the fact that there is no continuous translational symmetry implies that the crystallization of $\mathbb{Z}_{2}$ vortices into a superlattice comes with a finite spin gap $\Delta$. This is demonstrated in Fig. 2 (b), which shows the evolution of $\Delta$ vs $\psi$ for 17 optimal crystals. The gap is indeed finite everywhere inside the $\mathbb{Z}_{2}$VC phase. Its behaviour is non-monotonic and asymmetric with respect to the sign of $K$ (it is significantly larger for $K\!<\!0$ than for $K\!>\!0$). The softening of the gap in the vicinity of the C-IC transition ($\psi\!=\!0$) is in accord with the flattening of the ground state energy landscape [see scale in the vertical axis of Fig. 2 (a)] and the recovery of the true Goldstone mode at the HAF point.

Of particular interest is the region above $\psi\!=\!0.25\pi$, which shows not only a softening of the spin gap itself [Fig. 2 (b)], but also a significant accumulation of low-$\omega$ spectral weight [Fig. 5 (a, b)], reflecting the incipient frustrated Kitaev point. Strong quantum fluctuations may thus render this region susceptible to new collective physics that goes beyond our semiclassical analysis, see e.g., the recent study Maksimov et al. (2019) and Li et al. (2015a); Kos and Punk (2017).

VI Discussion

The prediction Rousochatzakis et al. (2016) that the coplanar 120∘ order of triangular Heisenberg antiferromagnets becomes immediately unstable under an infinitesimal Kitaev anisotropy, giving way to mesoscopic $\mathbb{Z}_{2}$ vortex crystals, has triggered a significant interest in the community Becker et al. (2015); Catuneanu et al. (2015); Shinjo et al. (2016); Trebst (2017); Li et al. (2015a); Kos and Punk (2017); Kishimoto et al. (2018); Maksimov et al. (2019), and remains to be explored and verified experimentally. At present, materials that have been discussed in this context, including the iridate Ba3IrTi2O9 Dey et al. (2012); Lee et al. (2017); Becker et al. (2015); Catuneanu et al. (2015), the mixed-valence iridate Ba3InIr2O9 Dey et al. (2017), and the rare-earth compound YbMgGaO4 (YMGO) Li et al. (2015b, c); Maksimov et al. (2019), suffer either from intrinsic disorder and impurities or additional complex anisotropic interactions Trebst (2017).

The $\mathbb{Z}_{2}$ vortex crystals can be detected by small-angle neutron or x-ray scattering methods, in analogy to 1D soliton lattices in modulated antiferromagnets (such as Ba2CuGe2O7 Zheludev et al. (1997)) or skyrmion lattices in chiral ferromagnetic helimagnets (such as MnSiMühlbauer et al. (2009) or Cu2OSeO3Langner et al. (2014)). Furthermore, the strongly inhomogeneous magnetization profile near the defected cores of the $\mathbb{Z}_{2}$ vortices would give rise to characteristic static hyperfine field distributions, which could be probed by NMR or $\mu$SR.

In this work, we have demonstrated that the $\mathbb{Z}_{2}$ vortex crystals can also be diagnosed in dynamical spectroscopic experiments in a more direct way. We have shown, in particular, that the collective spin dynamics of $\mathbb{Z}_{2}$ vortex crystals bears two of their most characteristic properties, the large vortex size and the nonlinear, particle-like nature of their defected cores. These show up with a characteristic multi-fragmented intensity profile at intermediate and high frequencies and a set of intense, fully fledged phonon-like modes at low frequencies. While certain aspects will be modified in higher orders of the $1/S$ expansion (for example, the characteristic high-frequency intensity profile will be further modified by the effect of the magnon decays which are known to be present for non-colinear magnetic orders Chernyshev and Zhitomirsky (2009); Winter et al. (2017b)), the main qualitative predictions can be used as ‘smoking guns’ for $\mathbb{Z}_{2}$ vortex crystals in appropriate materials.

Acknowledgement – We thank A. Chernyshev and P. Maksimov for helpful discussions. This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0018056. We also acknowledge the support of the Minnesota Supercomputing Institute (MSI) at the University of Minnesota.

Appendices

In these Appendices we provide auxiliary information and technical details on the magnetic unit cells of the $\mathbb{Z}_{2}$VC superstructures (App. A), the computation of the linear spin-wave spectra (App. B), the computation of the DSF and the INS intensities (App. C), as well as the INS profiles for a series of sixteen $\mathbb{Z}_{2}$VC’s (App. D).

Appendix A Magnetic unit cells

To perform the semiclassical $1/S$ expansion one needs to deduce the magnetic unit cell for each optimal $\mathbb{Z}_{2}$VC superstructure. It turns out that the spanning vectors, $\mathbf{T}_{1}$ and $\mathbf{T}_{2}$, of the magnetic unit cell are of two possible types, depending on the value of $d$ and the sign of the Kitaev interaction $K$. In the first type [Fig. 6(a)], the spanning vectors connect the centers of the domains belonging to the same sublattice, i.e., they connect (A,B,C)$\to$(A,B,C). This type of the magnetic unit cell encloses $\mathcal{N}_{m}\!=\!12\,d^{2}$ spins. In the second type [Fig. 6(b)], which was overlooked in Ref. Rousochatzakis et al. (2016), the spanning vectors connect one sublattice to another [(A,B,C)$\to$(B,C,A) for $K\!>\!0$ and (A,B,C)$\to$(C,A,B) for $K\!<\!0$]. This type of magnetic unit cell has $\mathcal{N}_{m}\!=\!4\,d^{2}$ spins. The conditions for $d$ and $K$ that give the two different types of magnetic unit cells, along with the associated spanning vectors and number of spins $\mathcal{N}_{m}$ in the magnetic unit cell are summarized in Table 1.

Appendix B Linear Spin Wave (LSW) analysis

In order to study the collective spin dynamics on top of a given optimal $\mathbb{Z}_{2}$VC, we must first relabel the spin sites $i\to({\bf R},\mu)$, where ${\bf R}=n_{1}{\bf T}_{1}+n_{2}{\bf T}_{2}$ is the position of the magnetic unit cell ($n_{1}$ and $n_{2}$ are integers), and $\mu\!=\!1-\mathcal{N}_{m}$ is the sublattice index inside the magnetic unit cell. Accordingly, we rewrite the spin ${\bf S}_{i}$ and its physical position ${\bf r}_{i}$ as

[TABLE]

respectively, where $\bm{\rho}_{\mu}$ is the sublattice vector associated to the $\mu$-th sublattice. The Hamiltonian is then written as

[TABLE]

where

[TABLE]

${\bf t}_{\mu\nu}$ is a primitive translation of the superlattice such that the spins at sites $i=({\bf R},\mu)$ and $j=({\bf R}+{\bf t}_{\mu\nu},\nu)$ interact with each other via $\bm{\mathcal{J}}_{\mu\nu}$, and

[TABLE]

where $\bm{\epsilon}\in\{{\bf a},{\bf b},{\bf c}\}$ and

[TABLE]

see Fig. 1. Next, for each site $i=({\bf R},\mu)$, we introduce local reference frames

[TABLE]

such that $\widetilde{\mathbf{z}}_{i}$ coincides with the direction of spin ${\bf S}_{i}$ in the classical ground state. The spin is then rotated into this local frame of reference by a unitary rotation matrix $\mathbf{U}_{\mu}$,

[TABLE]

The matrix $\mathbf{U}_{\mu}$ can be constructed using the polar and azimuthal angles $(\theta_{\mu},\phi_{\mu})$ associated with the direction of the spin in the classical ground state,

[TABLE]

Plugging into the Hamiltonian gives

[TABLE]

where $\bm{\mathcal{T}}_{\mu\nu}=\frac{1}{2}\bm{U}_{\mu}\bm{\mathcal{J}}_{\mu\nu}{\bf U}_{\nu}^{-\!1}$. Next, we perform a Holstein-Primakoff transformation Holstein and Primakoff (1940), and rewrite the spin operators $\widetilde{S}_{{\bf R},\mu}^{x}$, $\widetilde{S}_{{\bf R},\mu}^{y}$, and $\widetilde{S}_{{\bf R},\mu}^{z}$ in terms of bosonic creation and annihilation operators $a_{{\bf R},\mu}^{\dagger}$ and $a_{{\bf R},\mu}$ to lowest order as

[TABLE]

Then the Hamiltonian can be expanded in powers of $1/\sqrt{S}$,

[TABLE]

where the zeroth-order term

[TABLE]

represents the classical energy $E_{cl}$, the first-order term

[TABLE]

vanishes because we expand around the classical ground state, and the second-order term is

[TABLE]

where

[TABLE]

Using Fourier transform (where ${\bf q}$ belongs to the magnetic BZ)

[TABLE]

defining $\bm{\delta}_{\mu\nu}=({\bf R}+\bm{\rho}_{\mu})-({\bf R}+{\bf t}_{\mu\nu}+\bm{\rho}_{\nu})$, and symmetrizing with respect to $\mathbf{q}\rightarrow-\mathbf{q}$, we obtain

[TABLE]

where

[TABLE]

or in matrix form

[TABLE]

where ${\bf x}_{\mathbf{q}}=(a_{1,\mathbf{q}}\,,...\,,a_{\mathcal{N}_{m},\mathbf{q}}\,,a_{1,-\mathbf{q}}^{\dagger}\,,...\,,a_{\mathcal{N}_{m},-\mathbf{q}}^{\dagger})^{\text{T}}$, and ${\bf H}_{\mathbf{q}}$ is a $(2\mathcal{N}_{m})\times(2\mathcal{N}_{m})$ matrix. The diagonalization of ${\bf H}_{\mathbf{q}}$ involves introducing a new set of Bogoliubov quasiparticle operators Bogoliubov (1947); Blaizot and Ripka (1985),

[TABLE]

obtained from ${\bf x}_{\mathbf{q}}$ by a unitary canonical transformation ${\bf x}_{\mathbf{q}}=\mathbf{T_{q}}\cdot{\bf y}_{\mathbf{q}}$. The transformation must be such that the new bosons satisfy the bosonic commutation relation which, in terms of $\mathbf{T_{q}}$, gives the condition $\mathbf{T_{q}}^{\dagger}\cdot\mathbf{g}\cdot\mathbf{T_{q}}=\mathbf{g}$, where $\mathbf{g}=\text{diag}(\mathbf{I},-\mathbf{I})$ and $\mathbf{I}$ is a $\mathcal{N}_{m}\times\mathcal{N}_{m}$ unitary matrix. The matrix $\mathbf{T_{q}}$ can then be found by solving the eigenvalue equation (in matrix form)Blaizot and Ripka (1985)

[TABLE]

where $\bm{\Omega}_{\mathbf{q}}=\mathbf{T_{q}}^{\dagger}{\bf H}_{\mathbf{q}}\mathbf{T_{q}}=\text{diag}(\bm{\omega}_{\mathbf{q}},-\bm{\omega}_{\bf q})$, and $\bm{\omega}_{\mathbf{q}}$ is a diagonal matrix within elements $\{\omega_{1,\mathbf{q}},\omega_{2,\mathbf{q}},\ldots,\omega_{\mathcal{N}_{m},\mathbf{q}}\}$.

Appendix C Dynamical structure factor (DSF) and inelastic neutron scattering (INS) intensity

The DSF $S^{\alpha\beta}(\mathbf{Q},\omega)$ is given by the Fourier transform of the spin-spin correlations

[TABLE]

where the $\alpha$-th component of the spin on the sublattice $\mu$ is given by

[TABLE]

and $\xi_{\mu}^{\alpha}=[\mathbf{U}_{\mu}^{\!-\!1}]^{\alpha,1}-i[\mathbf{U}_{\mu}^{\!-\!1}]^{\alpha,2}$, $\lambda_{\mu}^{\alpha}=[\mathbf{U}_{\mu}^{\!-\!1}]^{\alpha,3}\;.$ Note that the third term in $S_{{\bf R},\mu}^{\alpha}$ can be dropped when calculating the DSF since this term only describes the reduction of the static ordered moment due to magnon population.

The Fourier transform of the spin component is given by

[TABLE]

where we used the relation $\mathbf{Q}=\mathbf{k}+{\bm{\tau}}$, where ${\bf Q}$ is the momentum transfer, $\mathbf{k}$ is a wavevector inside the first magnetic BZ, and ${\bm{\tau}}=n_{1}\mathbf{G}_{1}+n_{2}\mathbf{G}_{2}$ is a primitive vector of the reciprocal lattice of the superstructure, which satisfies $e^{i{\bm{\tau}}\cdot\mathbf{R}}=1$ for all ${\bf R}$. Then the DSF becomes

[TABLE]

where $\mathbf{V}_{\bm{\tau}}^{\alpha}$ is a vector array of coefficients given by

[TABLE]

Using the Bogoliubov transformation, we then obtain

[TABLE]

where the correlation functions of the bosonic quasiparticles are determined by

[TABLE]

where $n(\omega_{\gamma,\mathbf{k}})=[e^{\hbar\omega_{\gamma,\mathbf{k}}/(k_{B}T)}-1]^{-1}$ is the Bose factor at temperature T. At $T=0$, we therefore end up with

[TABLE]

Finally, the INS intensity $\mathcal{I}(\mathbf{Q},\omega)$ is given by the expressionSquires (2012)

[TABLE]

Appendix D Representative INS profiles

Figure 7 shows the evolution of the INS intensity $\mathcal{I}(\mathbf{Q},\omega)$ for sixteen representative $\mathbb{Z}_{2}$VC’s, as we depart away from the Heisenberg point ($\psi=0$) for both positive (panel c) and negative Kitaev interaction (panel d). The intensity profiles are shown along special symmetry directions in momentum space, see panels (a) and (b). The shift of the positions of the ‘phonon-like’ modes is highlighted by a red dashed curve. These modes follow the positions of the static structure factor. For $K>0$, the positions move from the corner of the BZ $\mathrm{K}$ towards the $\mathrm{\Gamma}$ point, whereas for $K<0$ they move along the directions $\mathrm{K}\rightarrow\mathrm{M}$. The yellow dashed line in panel (c) shows the accumulation of low-$\omega$ spectral weight as we approach the frustrated Kitaev point ($\psi=\pi/2$), see main text.