Disorder effects in spiral spin liquids: Long-range spin textures, Friedel-like oscillations, and spiral spin glasses

Spiral spin liquids are correlated states of matter in which a frustrated magnetic system evades order by fluctuating between a set of (nearly) degenerate spin spirals. Here, we investigate the response of spiral spin liquids to quenched disorder in a $J_1-J_2$ honeycomb-lattice Heisenberg model. At the single-impurity level, we identify different order-by-quenched-disorder mechanisms and analyze the ensuing spin textures. In particular, we show that the latter generally display Friedel-like oscillations, which encode direct information about the spiral contour, i.e., the classical ground-state manifold. At finite defect concentrations, we perform extensive numerical simulations and characterize the resulting phases at zero temperature. As a result, we find that the competition between incompatible order-by-quenched-disorder mechanisms can lead to spiral spin glass states already at low to moderate disorder. Finally, we discuss extensions of our conclusions to nonzero temperatures and higher-dimensional systems, as well as their applications to experiments.

Figure 1

Schematic phase diagrams illustrating the response of spiral spin liquids to quenched bond disorder in the classical limit in (a) $D=2$ and (b) $D=3$ space dimensions. The axes represent temperature $T$ and disorder strength $\mathrm{\Delta }$; $T_c$ and ${\mathrm{\Delta }}_c$ correspond to the critical values above which the system enters a paramagnetic (PM) and a spin glass (SG) phase, respectively. The latter has a freezing temperature $T_f$ and, at intermediate $\mathrm{\Delta }$, exhibits a “spiral spin-glass” (SSG) regime. ${\mathrm{O}}_{1,2,3}$ represent possibly different phases with spiral LRO, selected either by quenched disorder (${\mathrm{O}}_1$ and ${\mathrm{O}}_3$) or thermal fluctuations (${\mathrm{O}}_2$). The red shading encodes the extent to which the static spin structure factors delineates the spiral contour or surface.

Figure 2

(a) Honeycomb lattice with unit cell (green box) and the primitive vectors ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$. For $J_1>0$, the angle between two spins inside the same unit cell is given by $\pi +{\theta }_{\mathbf{Q}}$. (b) Ground-state phase diagram of the clean classical $J_1-J_2$ Heisenberg model as function of $\alpha =J_2/\left|J_1\right|$, valid for both ferromagnetic and antiferromagnetic $J_1$. The label COL stands for “collinear” and can refer either to a ferromagnetic or Néel state, while ${\mathrm{SSL}}_1$ and ${\mathrm{SSL}}_2$ are different SSL regimes. (c) Spiral contours representing the ground-state manifold of the clean Hamiltonian for three different ratios $\alpha =J_2/\left|J_1\right|$. In each case, the spiral contour can be parameterized by a polar angle ${\varphi }_{\mathbf{Q}}$, which is measured with respect to the origin and $\mathbf{K}=(0,4\pi /\left(3\sqrt{3}a\right))$ in the ${\mathrm{SSL}}_1$ and ${\mathrm{SSL}}_2$ regimes, respectively.

Figure 3

Illustration of finite-size effects in the clean limit. Results are shown for $\alpha \equiv J_2/J_1=0.2$ and hexagonal samples with open zigzag boundary conditions. (a) Finite-size scaling of the ground-state energy per site (circles) fitted to a linear function (solid line). The dashed blue line indicates the energy of the infinite-system solution [Eqs. (2, 3, 4)] when placed into the cluster. In the limit $L\to \infty$, the results extrapolate to the value predicted by Eq. (5) (horizontal dashed line). (b) Static structure factor $\mathcal{S}\left(\mathbf{k}\right)$ of a single lowest-energy state for $L=40$. The red hexagon delimits a second Brillouin zone centered at $2\pi (1/3,1/\sqrt{3})$, while the white ring represents the spiral contour of the infinite system (shown as guide to the eye). The solid green and dashed cyan lines correspond to different high-symmetry directions. (c) Finite-size scaling of ${\max }_k\mathcal{S}\left(\mathbf{k}\right)$ normalized by the total number of sites $N$. (d) Full width at half maximum (FWHM) of the peak in $\mathcal{S}\left(\mathbf{k}\right)$ along a radial (blue circles) and tangential (red squares) cut through the second Brillouin zone.

Figure 4

(a) Schematic showing the angles ${\theta }_{\mathbf{Q}}$ between NN spins in the same unit cell and ${\mathrm{\Theta }}_3$ between NNN spins separated by a vector ${\mathbf{a}}_3$. The green and red lines highlight the directions of the NN and NNN bond defects we considered for the ObQD arguments. (b) Momentum-dependence of the function ${\theta }_{\mathbf{k}}$ (3). Solid hexagons delimit Brillouin zones, and the inner dashed hexagon corresponds to the spiral contour at the Lifshitz point $\alpha =1/2$, which separates the ${\mathrm{SSL}}_1$ and ${\mathrm{SSL}}_2$ regimes. (c) Variation of $\mathbf{k}·{\mathbf{a}}_3$ with $k_y$. For wave vectors $\mathbf{k}$ in the spiral contour, this sets the angle ${\mathrm{\Theta }}_3$ between NNN spins, which is why we limited the results to the compact interval $]-\pi ,\pi ]$. We also took $k_y$ as the $y$ axis to match the visualization in (b). The horizontal dashed lines indicate the locations of the upper and lower edges of the $\alpha =1/2$ spiral contour. (b) and (c) can be used to determine under which circumstances bond defects do not produce a texture, for details see text.

Figure 5

ObQD results for a single vacancy. [(a) and (b)] Variation along the spiral contour of the energy correction ${\mathcal{H}}_{\mathrm{tex}}/J_1$ due to the spin texture for $\alpha =0.2$ and $0.8$, respectively. (c) Difference $\mathrm{\Delta }{\mathcal{H}}_{\mathrm{tex}}$ between the energy corrections at different high-symmetry points as a function of $\alpha$. (d) Same as (c), but zoomed into the vicinity of the first Lifshitz point $\alpha =1/6$. The sign change at $\alpha \approx 0.2148$ signals a shift in the wave vectors $\mathbf{Q}$ selected by the vacancy.

Figure 6

Different types of two-vacancy clusters on the honeycomb lattice with (a) NN and (b) NNN vacancies.

Figure 7

Friedel-like oscillations in the out-of-plane texture $\langle S_{i\mu }^y\rangle /\left(S\right|\delta J_1\left|\right)$ due to a weak NN bond defect in a system with [(a) and (b)] $\alpha =0.2$ and [(c) and (d)] $0.3$, obtained from the calculation outlined in Appendix pp4. The left (right) column displays results for the $\mu =0(\mu =1)$ sublattice, respectively. In all cases, the impurity is located at the origin and the data were obtained by performing the integral in Eq. (D25) numerically.

Figure 8

Numerical results for the four different parameter sets listed in Table 1. Top: Disorder-averaged static structure factors for systems of linear size $L=40$. The elements in the plots are the same as those in Fig. 3, except for the dotted white lines. These represent the spiral contour for the average ratio $\overline{\alpha }$ in each parameter set. (Middle and bottom) Finite-size scaling of the normalized height and FWHM of the Bragg peaks. The linear fits shown discarded the data points corresponding to the two smallest system sizes. Magnetic LRO is seen for sets A and C, but not for set B. For set D see text.

Figure 9

Same as Fig. 8, but for the parameter sets listed in Table 4. The disorder-averaged structure factors shown here are symmetrized by additionally averaging over $2\pi /3$ rotations. There is no magnetic LRO in any of the cases.

Figure 10

Histograms of the energy per site, $E/N$, of converged solutions of the classical minimization algorithm. The horizontal axis represents deviations relative to the lowest energy, $\mathrm{\Delta }E/N=max[{10}^{-8},(E-E_{\min })/N]$, while the vertical axis corresponds to the frequency with which solutions in a given energy window were obtained after ${10}^3$ different initializations. The arrows highlight the positions of the peaks with $\mathrm{\Delta }E/N<{10}^{-5}$. These histograms map out the local minima of the energy landscape for finite clusters of [(a)–(d)] the clean system and [(e)–(f)] single disorder realizations of set F. In both cases, the clusters had open zigzag boundary conditions.

Figure 11

Same as Fig. 8 but for systems with $p=5\%$ of vacancies and a maximum system size of $L=40$. The disorder-averaged structure factors shown here are symmetrized by additionally averaging over $2\pi /3$ rotations.

Figure 12

Illustration of scalar chiralities ${\chi }_{ijk}^1$ in a ground state of a clean system with $\alpha =0.2$, open zigzag boundary conditions, and linear system size $L=28$. The left and right figures correspond to two different choices for triplets of sites such that $i1, j1$, and $k1$ are next-nearest neighbors. In both cases, the values of ${\chi }_{ijk}^1$ are represented by triangles whose vertices point toward these three sites. The size of each triangle is proportional to the absolute values of ${\chi }_{ijk}^1$, while the binary color code represents its sign.

Figure 13

Plots of the lower branch ${\epsilon }_{\mathbf{k}0}$ of the linear spin-wave dispersion for $\alpha =0.2$ with [(a) and (b)] ${\varphi }_{\mathbf{Q}}=0$ and [(c) and (d)] $\pi /2$. In upper row, the dashed gray hexagons indicate the boundary of the first Brillouin zone, while the yellow vector corresponds to the ordering wave vector (a) $\mathbf{Q}=(2/3){cos}^{-1}[1/{\left(4\alpha \right)}^2-5/4]\widehat{\mathbf{x}}\approx 0.8353\widehat{\mathbf{x}}$ and (c) $\mathbf{Q}=(2/\sqrt{3}){cos}^{-1}[1/\left(4\alpha \right)-1/2]\widehat{\mathbf{y}}\approx 0.8345\widehat{\mathbf{y}}$. Plots (b) and (d) depict the radial dependence of the dispersion along rays with polar angles $\phi$ ranging from 0 (blue) to $\pi /2$ (gray).

Figure 14

Linear-response results for a horizontal (${\mathbf{R}}_{nm}=\mathit{0}$) NN bond defect with $\delta J_1<0$, obtained by using the ${\varphi }_{\mathbf{Q}}=0$ spin spiral as a reference state for the spin-wave expansion. [(a) and (b)] Real and imaginary parts, respectively, of ${\tilde{\chi }}_{\mathbf{k}0}^{11}(\omega =0)/\left(\right|\delta J_1\left|S\right)$. The dashed hexagons indicate the boundaries of the Brillouin zone. [(c) and (d)] Real-space configuration of $\langle S_{i0}^1\rangle /\left|\delta J_1\right|$ on a linear and logarithmic scale $\mathrm{sgn}\left(z\right)ln\left(\right|z/z_{\min }\left|\right)$, respectively. These in-plane textures correspond to a periodic system with $M=256$ and are depicted such that an impurity lies at the center of the unit cell outlined in black. (e) Finite-size scaling of the position $y_n$ of the $n$-th zero along the yellow line $x=48a$ in (d). The dashed black lines indicate the slope of power laws proportional to $\sqrt{M}$ and $M$. (f) Decay of the texture along the green high-symmetry line in (d). The color code was selected such that the brighteness of the shade of red increases with system size. In (b) and (c), white colors occur when the susceptibility exceeds the limits of the plot legends.

Figure 15

Similar to Fig. 14, but for a vertical (${\mathbf{R}}_{nm}={\mathbf{a}}_3$) NNN bond defect on the $\mu =0$ sublattice. To keep consistency with the ObQD result, the spin-wave calculations were performed by using the ${\varphi }_{\mathbf{Q}}=\pi /2$ coplanar spiral as a reference state.

Figure 16

Same as Fig. 14, but for a single vacancy on the $\mu =0$ sublattice.

Figure 17

Same as Fig. 16, but for $\alpha =0.3$.

Figure 18

Spatial decay of the texture generated by a weak NNN bond defect for $\alpha =30$. The two dashed lines correspond to power laws $1/y$ and $1/\sqrt{y}$, indicating that the behavior expected for a triangular-lattice antiferromagnet ($\alpha \to \infty$) emerges below a characteristic distance $R_{\mathrm{cross}}\left(\alpha \right)$ from the impurity.