Spin dynamics of the antiferromagnetic Heisenberg model on a kagome bilayer

TL;DR

This study explores the spin dynamics of a classical Heisenberg antiferromagnet on a kagome bilayer, revealing spin diffusion, magnon excitations, and effects of dilution, with implications for glassy behaviors in frustrated magnets.

Contribution

It provides the first detailed analysis of dynamical structure factors and dilution effects in kagome bilayer antiferromagnets using combined Monte Carlo and spin dynamics simulations.

Findings

Spin diffusion dominates low-frequency dynamics.

Persistent magnon excitations create half moon patterns.

Dilution leads to distinct dynamics of correlated and orphan spins.

Abstract

We study the spin dynamics of classical Heisenberg antiferromagnet with nearest neighbor interactions on a quasi-two-dimensional kagome bilayer. This geometrically frustrated lattice consists of two kagome layers connected by a triangular-lattice linking layer. By combining Monte Carlo with precessional spin dynamics simulations, we compute the dynamical structure factor of the classical spin liquid in kagome bilayer and investigate the thermal and dilution effects. While the low frequency and long wavelength dynamics of the cooperative paramagnetic phase is dominated by spin diffusion, weak magnon excitations persist at higher energies, giving rise the half moon pattern in the dynamical structure factor. In the presence of spin vacancies, the dynamical properties of the diluted system can be understood within the two population picture. The spin diffusion of the "correlated" spin clusters is mainly driven by the zero-energy weather-van modes, giving rise to an autocorrelation function that decays exponentially with time. On the other hand, the diffusive dynamics of the quasi-free "orphan" spins leads to a distinctive longer time power-law tail in the autocorrelation function. We discuss the implications of our work for the glassy behaviors observed in the archetypal frustrated magnet SrCr$_{9p}$Ga$_{12-9p}$O$_{19}$ (SCGO).

Equations

$$\mathcal{H}$$

$$\mathcal{H}=\frac{J}{2}\sum _{⊠}{\mathbf{L}}_{⊠}^2+\sum _{△}{\mathbf{L}}_{△}^2+\text{mbox}const.$$

$${\mathbf{L}}_{⊠}={\mathbf{L}}_{△}=0,$$

$$\frac{d{\mathbf{S}}_i}{dt}=-{\mathbf{S}}_i\times \frac{\partial \mathcal{H}}{\partial {\mathbf{S}}_i}=-J\sum _j{}^{\prime }{\mathbf{S}}_i\times {\mathbf{S}}_j,$$

$$\mathcal{S}(\mathbf{q},t)=⟨{\mathbf{S}}_{\mathbf{q}}(t)\cdot {\mathbf{S}}_{\mathbf{q}}^{\ast }(0)⟩,$$

$$\mathcal{S}(\mathbf{q},\omega )$$

$$A(t)=\frac{1}{N}\sum _i⟨{\mathbf{S}}_i(t)\cdot {\mathbf{S}}_i(0)⟩=\sum _{\mathbf{q}}\mathcal{S}(\mathbf{q},t),$$

$$\beta q^2\mathcal{S}(\mathbf{q},\omega )=\chi \frac{2D}{(\omega /q^2{)}^2+D^2},$$

Full text

I Introduction

The SCGO is one of the most intensely studied frustrated magnets obradors88 ; ramirez90 ; ramirez92 ; broholm90 ; uemura94 ; lee96a ; schiffer96 ; lee96b ; ramirez00 ; keren00 ; mendels00 ; bono05 ; iida12 ; yang15 . Thermodynamically, SCGO does not exhibit any signs of magnetic ordering down to temperatures $T_{g}=$ 3.5–7 K, depending weakly on the vacancy concentration $x=1-p$. Below $T_{g}$, the magnet enters an unconventional spin-glass phase. A cooperative paramagnetic regime, also known as a classical spin liquid, emerges at temperatures below the Curie-Weiss constant $\Theta_{\rm CW}\approx 500$ K. Geometrically, SCGO belongs to a class of frustrated Heisenberg antiferromagnets on the so-called bi-simplex lattices moessner98a ; moessner98b ; henley01 ; canals98 . These are networks of corner-sharing simplexes such as triangles and tetrahedra. Canonical examples include the pyrochlore moessner98a ; moessner98b and kagome chalker92 ; garanin99 ; henley09 ; chern13 antiferromagnets. In SCGO, the Cr3+ ions with spin $S=3/2$ reside on a two-dimensional lattice consisting of corner-sharing tetrahedra and triangles, known as the kagome bilayer or pyrochlore slab, as shown in Fig. 1. The strong short-range spin correlations in the low-temperature liquid phase result from the constraints of zero total spin in every simplex, a condition that minimizes the nearest-neighbor exchange interactions on such unit.

Considerable experimental efforts have been devoted to understanding the unusual spin glass phase in SCGO ramirez90 ; ramirez92 ; broholm90 ; uemura94 ; lee96b ; keren00 ; mendels00 ; bono05 ; yang15 . Despite the characteristic field-cooled and zero-field-cooled hysteresis in the bulk susceptibility, several dynamical properties of its glassy phase are distinctly different from those of conventional spin glasses. These include the quadratic $T^{2}$ behavior of the specific heat ramirez90 ; ramirez92 , the linear $\omega$-dependent dynamical susceptibility $\chi^{\prime\prime}$ lee96b , and a significantly weaker memory effect samarakoon16 . Taken together, these features suggest that SCGO belongs to a new state of glassy magnets, dubbed the spin jam yang15 ; samarakoon16 , that include several other magnetic compounds yang16 ; samarakoon17 . The source of this unusual dynamical phase in SCGO, however, remains to be clarified. One plausible scenario is that quantum fluctuations transform the macroscopic degeneracy associated with the classical spin liquid of the kagome bilayer into the rugged energy landscape of spin jam yang15 ; klich14 . It remains to be shown how the unusual glassy behaviors of the spin jam evolve from the spin dynamics of the cooperative paramagnet.

Toward this goal, we present in this paper the first systematic study of the dynamical properties of the bilayer-kagome classical spin liquid. By combining Monte Carlo simulations with energy-conserving Landau-Lifshitz dynamics, we compute the dynamical structure factor of the liquid regime. At the energy scales of the exchange interaction, we find signals of spin-wave excitations in the form of half moon pattern, replacing the pinch-point singularity of the static structure factor. On the other hand, the low-energy dynamics is dominated by spin-diffusion driven mostly by the zero-energy modes. The diffusion constant is found to depend weakly on temperature, but decrease significantly with increasing vacancy densities.

Our results will also serve as an important benchmark against which dynamical behaviors induced by other perturbations can be compared. Of particular interest are those perturbations, such as quantum order-by-disorder, that give rise to glassy dynamics characteristic of either the conventional spin-glass or the exotic spin-jam states. It is also worth noting that the dynamical properties of classical spin liquid has been extensively studied for Heisenberg antiferromagnets on both pyrochlore moessner98a ; moessner98b ; conlon09 and kagome lattices keren94 ; robert08 ; taillefumier14 ; bilitewski19 . Another aim of this paper is thus to compare the spin dynamics of bilayer kagome against these two well studied bi-simplex frustrated magnets.

The rest of the paper is organized as follows. In section II, we discuss the ground-state manifold of Heisenberg antiferromagnet on the kagome bilayer. We also outline the numerical framework that combines Monte Carlo simulation with energy-preserving Landau-Lifshitz dynamics method for computing the dynamical structure factor of a classical spin liquid. Magnetic excitations revealed from the dynamical structure factor are discussed in Sec. III. In particular, half-moon features, which are the dynamical manifestation of the famous pinch-point structure at finite energies, are highlighted. Systematic analysis of the low-energy spin dynamics, which is dominated by diffusive modes, is presented in Sec. IV. We present in Sec. V dynamical features due to quenched disorder introduced by vacancies. Of particular interest is the emergence of quasi-free orphan spins that interact with each other through a week effective interaction mediated by the background spin-liquid. We conclude in Sec. VI with a brief summary and outlook on future work.

II Model and method

We consider the classical Heisenberg model with nearest neighbor interactions on the kagome bilayer

[TABLE]

Here $J>0$ is the antiferromagnetic exchange, $\langle ij\rangle$ denotes nearest-neighbor pairs, and the classical spins $\mathbf{S}_{i}$ are unit vectors. The kagome bilayer has a Bravais triangular lattice with a unit cell consisting of two corner sharing tetrahedra of opposite orientation. The bases of the tetrahedra in the two kagome layers are connected by triangle units; see Fig. 1. The triangle and tetrahedron are the regular simplexes with $q=3$ and $q=4$ corners, respectively. Importantly, because of this corner-sharing simplex structure, the exchange interaction can also be expressed as a sum of the squared total spin of both types of simplexes

[TABLE]

Here $\mathbf{L}_{\boxtimes}=\sum_{i\in\boxtimes}\mathbf{S}_{i}$ denotes total spin in a tetrahedron, $\mathbf{L}_{\triangle}=\sum_{i\in\triangle}\mathbf{S}_{i}$ denotes total spins of a triangle, and $\sum_{\boxtimes}$ and $\sum_{\triangle}$ indicate summation over tetrahedra and triangles, respectively, in the kagome-bilayer lattice. One can immediately see that the exchange energy is minimized by the condition that total spin of every simplexes is zero:

[TABLE]

The ground-state condition is confirmed by our Monte Carlo simulations. The fact that a macroscopic number of spin configurations satisfy the minimum energy condition leads to a classical spin liquid regime at temperatures $T\lesssim J$. Indeed, our Monte Carlo simulations show no signs of phase transition down to temperatures $T\approx 0.001J$, consistent with previous studies arimori01 ; kawamura02 ; sen11 ; sen12 . Instead, a spin-disordered phase with strong short-range correlation is obtained at low temperatures.

In general, there are two types of spin dynamics in the liquid regime. At short time scales, or high frequencies ($\omega\sim J$), there are spin-wave excitations corresponding to small amplitude deviations from the ground-state manifold. These excitations are similar to the magnons in unfrustrated magnets. On the other hand, the macroscopic number of zero modes, or the weather-vane modes, that connect different ground states dominate the long-time dynamical behaviors of the frustrated bi-simplex antiferromagnet. The resultant drifting of the system in the ground-state manifold gives rise to spin-diffusion behaviors and an exponential decaying spin autocorrelation. In the following, we discuss our simulation results within this general picture.

The equation of motion for classical spins is given by the Landau-Lifshitz equation

[TABLE]

where the prime indicates summation is restricted to the nearest neighbors of the $i$-th spin. Here we numerically integrate the Landau-Lifshitz equation to compute the dynamical structure factor of the classical spin liquid. Low temperature Monte Carlo simulations are first used to obtain spin configurations in equilibrium of a specified temperature. These are then used as the initial states for the energy-conserving precession dynamics simulations. An efficient semi-implicit integration algorithm mentink10 is employed to integrate the above Landau-Lifshitz equation. The high efficiency of the algorithm comes from fact that it preserves the spin length at every time step and the energy values are well conserved with time irrespective of the step size or the time span of the simulation. From the numerically obtained spin trajectories $\mathbf{S}_{i}(t)$, we compute the dynamical correlation function $\mathcal{S}(\mathbf{q},t)$

[TABLE]

where $\mathbf{S}_{\mathbf{q}}(t)\equiv\sum_{i}\mathbf{S}_{i}(t)\exp(i\mathbf{q}\cdot\mathbf{r}_{i})/\sqrt{N}$ is the spatial Fourier transform of the instantaneous spin configuration, and $\langle\cdots\rangle$ denotes the ensemble average over independent initial states of a given temperature. The dynamical structure factor is then given by

[TABLE]

which is essentially the space-time Fourier transform of the spin-spin correlator $C_{ij}(t)\equiv\langle\mathbf{S}_{i}(t)\cdot\mathbf{S}_{j}(0)\rangle$.

III Magnons and half moon patterns

The intensity plot of the scaled dynamical structure factor $\beta\,\mathcal{S}(\mathbf{q},\omega)$, where $\beta=1/T$, is shown in Fig. 2 for four different temperatures. The spin excitations are clearly dominated by the low-energy quasi-static fluctuations that extend over most of the Brillouin zone; also see Fig. 3 for the density plots of $\mathcal{S}(\mathbf{q},\omega)$ in the reciprocal space at constant energies. At low temperatures, the similar patterns of the quasi-static excitations indicating nontrivial scaling behaviors to be discussed below. Moreover, the relatively weak excitations at higher energies $\omega\gtrsim J$ result from the magnon fluctuations in the vicinity of an instantaneous ground state. Contrary to the kagome antiferromagnets robert08 ; taillefumier14 , no sharp propagating modes can be seen in the dynamical structure factor of the bilayer kagome.

The static structure factor, corresponding to Fig. 3(a) with $\omega=0$, exhibits sharp pinch points which are a hallmark of highly correlated spin liquid in bi-simplex frustrated magnets. The source of these singularities can be attributed to the ground-state constraints Eq. (3), which translates into a solenoid condition $\nabla\cdot\mathbf{B}=0$ for an emergent “magnetic” or flux field that is a coarse-grained representation of the spin configuration. This in turn gives rise to an anisotropic dipolar-like correlation of the flux field, which manifests itself as the pinch-point singularity in the reciprocal space garanin99 ; isakov04 ; henley05 .

At finite temperatures, the width of the pinch-point is roughly proportional to $\sqrt{T}$ conlon10 . Interestingly, the pinch point is also smeared with increasing $\omega$, and is replaced by the so-called half moon pattern at $\omega\gtrsim J$, as shown in Fig. 3 (b) and (c). Similar features, called the “excitation rings” have been observed in the finite-energy dynamical structure factor of the coplanar spin liquid phase of kagome robert08 . It has been pointed out that the half-moon can be viewed as the pinch-point with a dispersive dynamical flux field yan18 . These crescent patterns at high energies are the remnants of the propagating magnons mentioned above. Compared with the coplanar phase in kagome, the half-moon feature is much weaker in the liquid phase of bilayer kagome, indicating less rigid local structures in the instantaneous ground state.

IV Spin diffusion

The relatively weak half moon excitations also indicate a dominating spin diffusive dynamics in bilayer kagome. In general, spin diffusion dominates the excitation spectrum of disordered Heisenberg systems in the hydrodynamic limit halperin69 ; muller88 . In frustrated magnets, this diffusion results from the macroscopic number of zero-energy modes in the instantaneous ground state, causing the system to wander around the degenerate manifold. One particular manifestation of this diffusion is the decay of the spin autocorrelation function

[TABLE]

where again $\langle\cdots\rangle$ is the thermal average, which is achieved through averaging over independent initial states from Monte Carlo simulations. Fig. 4(a) shows $A(t)$ as a function of time for various temperatures obtained from a $L=30$ system. The decay of the autocorrelation function is found to be exponential $A(t)\sim\exp(-t/\tau)$ in the low temperature regime, and the numerically extracted time constant $\tau$ is shown in Fig. 4(b) as a function of temperature.

The nearly linear segment in the log-log plot suggests a power-law dependence $\tau\sim T^{-\zeta}$, where the numerically obtained exponent $\zeta=0.924\pm 0.015$, which is close to 1 as predicted by a soft-spin Langevin dynamics model for frustrated magnets with macroscopic ground-state degeneracy conlon09 . The exponential decay with $\tau\sim 1/T$ is consistent with the zero-mode driven spin-diffusion scenario moessner98a ; conlon09 , since the zero modes have no intrinsic energy scales, and the only relevant one is set by the inverse temperature. This result is also in stark contrast to the high-$T$ conventional paramagnet in which the spin-diffusion is shown to produce a power-law tail in the autocorrelation function muller88 ; gerling90 ; gerling90b ; constantoudis97 ; bagchi13 .

While the microscopic mechanisms of spin-diffusion could be thermal or quantum fluctuations, or the large number of zero modes in frustrated systems, fundamentally the diffusive spin dynamics is related to the fact that the total spin density $\mathbf{m}=\sum_{i}\mathbf{S}_{i}/N$ is a constant of the equation of motion. By combining the continuity equation $\partial\mathbf{m}/\partial t+\nabla\cdot\mathbf{j}=0$ with a phenomenological Fick’s law for local spin current $\mathbf{j}=-D\nabla\mathbf{m}$, one arrives at the familiar diffusion equation for the magnetization density. In the hydrodynamic regime, this introduces a diffusion timescale $\tau_{d}=1/Dq^{2}$ for perturbations characterized by wave vector $\mathbf{q}$. This is indeed confirmed by our dynamical simulations. Fig. 5(a) shows the time dependence of the dynamical correlation function $\mathcal{S}(\mathbf{q},t)$ for various wave vectors. Each curve is obtained after averaging over 500 independent initial states from Monte Carlo simulations. The correlation function is found to decay exponentially with time: $\mathcal{S}(\mathbf{q},t)\sim\exp(-t/\tau_{d})$, where the numerically extracted relaxation time, shown in Fig. 5(b), is isotropic in the reciprocal space and follows nicely the expected behavior $\tau_{d}^{-1}=Dq^{2}$ for wave vectors close to the Brillouin zone center.

More generally, here we try to understand our results using the hydrodynamic theory of the paramagnetic state, which suggests a generalized dynamical susceptibility: $\chi(\mathbf{q},\omega)=-\chi(\mathbf{q})\,{Dq^{2}}/(Dq^{2}-i\omega)$ foster75 ; halperin69 , where $\chi(\mathbf{q})$ is the static susceptibility at wave vector $\mathbf{q}$ and $D$ is the spin diffusion coefficient. The dynamical structural factor is obtained through the fluctuation-dissipation theorem: $\mathcal{S}(\mathbf{q},\omega)\approx 2[n_{B}(\omega)+1]\,{\rm Im}\chi(\mathbf{q},\omega)$, where $n_{B}(\omega)=1/(e^{\beta\omega}-1)$. In the $\omega\ll T$ regime, assuming $\chi(\mathbf{q})\approx\chi$ is a constant for small $q$, the dynamical structure factor can be expressed in a scaling form

[TABLE]

A similar result can be obtained from the Langevin soft-spin model conlon09 . By plotting $\beta q^{2}\mathcal{S}$ versus $\omega/q^{2}$, we find nice data collapsing from curves of different wave vectors, as shown in Fig. 6 (a) and (b), indicating a static susceptibility that indeed weakly depends on $\mathbf{q}$ for wave vectors close to zone center. On the other hand, we find that the collapsing of data points from different temperatures is not very satisfactory. Instead, we fit the collapsed data points from each temperature with the Lorentzian scaling function in Eq. (8) and extract both the spin diffusion coefficient $D$ and static susceptibility $\chi$. The temperature dependence of these two quantities are shown in Fig. 6(c) and (d). The spin-diffusion coefficient decreases quite appreciably with temperature, while the susceptibility remains roughly the same within the error bars.

V Dilution effects

We next investigate the effect of dilution on the spin dynamics of the liquid phase. Previous studies have indicated that dilution with non-magnetic vacancies does not induce the spin-glass behavior of SCGO henley01 ; shender93 . In fact, the condition Eq. (3) is satisfied for every simplex, for both tetrahedron and trinagle, in the ground states even for strong dilution henley01 . Consequently, a macroscopic degeneracy remains and the low-$T$ phase seems well approximated by a Coulombic classical spin liquid. To demonstrate this, we compute the dynamical structure factor of the diluted kagome bilayer using a combination of Monte Carlo simulations with the energy-conserving Landau-Lifshitz dynamics simulations. Fig. 7 shows the computed $\mathcal{S}(\mathbf{q},\omega)$ at $T=0.01J$ for four different vacancy concentrations. In addition to the thermal average over independent initial states, the $\mathcal{S}(\mathbf{q},\omega)$ of the diluted system is computed with a further average over the disorder, or different vacancy configurations. Interestingly, we find no dramatic change to the calculated $\mathcal{S}(\mathbf{q},\omega)$ even for vacancy density as high as $x=0.5$. The quasi-static excitations show similar patterns for all concentrations, although both the energy of spin-wave-like excitations at large $\omega$ and the bandwidth of the quasi-static excitations are slightly reduced with increasing vacancy concentrations.

Focusing on the small-$\omega$ and $\mathbf{q}$ regime, we found that the dynamical structure factor is still well approximated by the scaling function of Eq. (8), as shown in Fig. 8(a) and (b), indicating a dominating spin diffusion behavior. The diffusion coefficient $D$ extracted from the data-point collapsing is plotted in Fig. 8(c) as a function of $x$. The reduced diffusivity with increasing vacancy concentration indicates a longer relaxation time $\tau_{d}=1/Dq^{2}$, or a slower dynamics, caused by the disorder, although the system remains liquid-like. Fig. 8(d) shows the extracted static susceptibility $\chi$ in the $q\to 0$ limit versus vacancy concentration $x$. This trend is consistent with the two population picture since the quasi-free orphan spins dominate the low-$T$ static susceptibility, hence $\chi$ increases with the vacancy concentration.

Since the presence of vacancies does not change the liquid nature or the frustrated spin-interactions in the kagome bilayer, it is unclear whether the non-magnetic impurities introduce any new dynamical effect. On the other hand, the so-called orphan spins due to the dilution are known to induce nontrivial effects on the equilibrium properties of the kagome bilayer schiffer97 . The orphan spin corresponds to defect triangular simplex with only one surviving spin and two non-magnetic sites henley01 ; moessner99 . An example of the orphan spin is shown in Fig. 9. One can also think of orphan spin as connecting a $q=3$ triangular simplex and a $q=1$ point simplex, which is the spin itself. The orphan spin behaves as a quasi-free spin with a fractionalized length $S/2$ when perturbed by a magnetic field sen11 ; sen12 . Experimentally, these seemingly isolated free spins in diluted SCGO produce a Curie-like component in the static susceptibility even at temperatures well below $\Theta_{\rm CW}$ schiffer97 ; moessner99 . Detailed Monte Carlo simulations uncover a complex spin texture surrounding the defect simplex whose total spins indeed sum to $S/2$ sen11 ; sen12 . The fractionalized spin-texture also provides a natural explanation for the short-range oscillating signal observed in nuclear magnetic resonance limot02 .

An intuitive argument for the fractionalized $S/2$ orphan spins was originally given by Henley from the viewpoint of bi-simplex structure henley01 . Because each spin is shared by two simplexes in the bi-simplex lattices such as the kagome-bilayer, the total magnetization can be written as $\mathbf{M}_{\rm tot}=\frac{1}{2}\sum_{\alpha}\mathbf{L}_{\alpha}$, where $\alpha$ now runs over tetrahedral, triangular, and $q=1$ simplexes in the presence of orphan spins. In the ground states, total spin of each tetrahedron and triangle simplex remains zero, as evidenced by Monte Carlo simulations henley01 . As a result, the total magnetization of the system becomes $\mathbf{M}_{\rm tot}=\frac{1}{2}\sum_{\alpha}^{q=1}\mathbf{L}_{\alpha}$, where now the summation is restricted to $q=1$ single-point simplexes. As shown above, such $q=1$ simplex is just the orphan spin itself, so we have $\mathbf{M}_{\rm tot}=\frac{1}{2}\sum_{i\in{\rm orphan}}\mathbf{S}_{i}$, which also means that each orphan spin can be viewed as a quasi-free spin with a fractionalized length $S/2$ when in a magnetic field henley01 .

A natural question then is what is the dynamical manifestation of these orphan spins. To this end, we examine the spin-spin autocorrelation function $A(t)$ defined in Eq. (7). Fig. 10 shows the semi-log plot of autocorrelation functions with and without vacancies obtained from our dynamical simulations. In both cases, the initial decay of the autocorrelation can be well described by an exponential function, i.e. $A(t)\sim e^{-t/\tau}$ for small $t$. However, while the exponential decay persists to longer time-scales in the non-diluted system, the autocorrelation function of the diluted magnet exhibits a long-time tail, indicating a significantly reduced decline rate of the spin-autocorrelation.

This two-stage relaxation of the autocorrelation function can be understood in the framework of the two population picture schiffer97 discussed previously, namely, the classical spin liquid of kagome bilayer can be viewed of consisting of the “correlated” population which forms momentless clusters ($\mathbf{L}_{\boxtimes}=\mathbf{L}_{\triangle}=0$) and the population of quasi-free “orphan” spins that weakly interact with each other sen12 . Of course, at very strong dilution, the set of free spins also include those completely isolated magnetic ions henley01 . Dynamically, these two populations of spins are expected to behave differently. As discussed in Sec. IV, spin diffusion in the classical spin liquid, which is mainly driven by the zero energy modes, results in an autocorrelation function $A(t)$ which decays exponentially with time. On the other hand, since the vacancy-induced orphan spins can be viewed as nearly free spins, one expects their dynamical behavior to be similar to that of uncorrelated paramagnet. Earlier works have shown that diffusion of Heisenberg spins in an uncorrelated paramagnet leads to a power-law tail in the autocorrelation function, i.e. $A(t)\sim 1/t^{\alpha}$ muller88 ; gerling90 ; gerling90b ; constantoudis97 ; bagchi13 , where the exponent $\alpha$ depends strongly on the dimensionality. For 2D Heisenberg magnet, it is estimated to be $\alpha\approx 1.05\pm 0.025$ muller88 .

To verify the above picture, we present detailed examination of both the short-time and long-time behaviors of spin autocorrelation function for a diluted kagome bilayer with a vacancy density $x=0.3$. We performed extensive Monte Carlo and Landau-Lifshitz dynamics simulations on a $L=48$ lattice (with total number of spins $N=7\times L^{2}=16128$). Over 50 independent realizations of the disorder were constructed, and for each vacancy configuration, 100 independent initial spin states are prepared at the simulation temperatures. Fig. 11(a) shows semi-log plot of spin-autocorrelation at various temperatures. At short time scale, the decrease of $A(t)$ can be reasonably approximated by an exponential decay $A(t)\sim e^{-t/\tau(T)}$, similar to the undiluted case, with a temperature-dependent decay time constant $\tau$. The numerically extracted relaxation time $\tau$, shown in Fig. 12(a), again exhibits a power-law dependence on temperature $\tau\sim 1/T^{\zeta}$, with an exponent $\zeta\sim 0.952\pm 0.017$, which is similar to the undiluted case. As discussed in the previous section, the spherical approximation for the classical spin liquid predicts an exponent $\zeta=1$. It is unclear whether the deviation here is due to finite-size effect or the soft-spin approximation.

At longer time scales, the decay of the autocorrelation function slows down and turns into a power-law tail, $A(t)\sim A_{\infty}+\mathcal{C}(T)/t^{\alpha}$, with the same exponent $\alpha$ for different temperatures. Interestingly, as shown in Fig. 12(b), the amplitude of this power-law tail also exhibits a power-law dependence $\mathcal{C}\sim 1/T^{\eta}$, with an exponent $\eta=2.123\pm 0.021$. It is also worth noting that the decay of spin-autocorrelation saturates to a small but non-zero value at large times, as shown in Fig. 12(a). Similar results, which can be attributed to finite-size effect, have been reported in the spin-dynamics of uncorrelated classical Heisenberg chains constantoudis97 ; bagchi13 .

To better understand the power-law decay and the origin of the nonzero asymptotic $A_{\infty}$ in the diluted systems, we consider the dynamics of an orphan spin. At any finite temperatures, the total spin of individual simplex does not vanish identically, hence the ground-state condition Eq. (3) is not strictly satisfied. Indeed, the fluctuation of simplex magnetization is given by $\langle\mathbf{L}_{\alpha}^{2}\rangle\sim T/J$ moessner98a ; sen12 . This also indicates a non-zero coupling between orphan spins and the background correlated spin liquid. This residual coupling leads to incoherent precession of orphan spins and the exponential decay of the orphan-spin autocorrelation function. On the other hand, as shown in Ref. sen12 , there is an emergent effective interaction between the orphan spins. At low enough temperatures, their collective dynamics induced by this residual interaction thus slows down the exponential decay of autocorrelation function that is caused by coupling to the background spin liquid, and turns it into a power-law decay, similar to the anomalous spin diffusion in classical Heisenberg magnets at high temperatures muller88 .

Interestingly, the exponent $\alpha\approx 2.1$ obtained from our numerical fitting is significantly different from that of the 2D paramagnet. This unusual result could be attributed to the complex interaction between the orphan spins. As demonstrated in Ref. sen12, , there is an emergent Heisenberg exchange interaction between the orphan spins, that is determined by the charge-charge correlator of the underlying Coulomb spin liquid. Moreover, the sign (FM vs AFM) depends on whether the two orphan spins belong to the same kagome layer or not. It has been speculated whether this complex and potentially frustrated interaction might lead to glassy dynamics at low temperatures. Indeed, although it is believed that the autocorrelation function of spin-glass exhibits an stretched exponential decay at temperatures above the glass transition $T_{g}$, general scaling rules near the glass transition point imply a cutoff power-law palmer84 ; sompolinsky82 ; ogielski85 ; pinkvos90 ; keren96 , such as the Ogielski form $A(t)\sim t^{-\alpha}\,\exp[-(\lambda t)^{\beta}]$, where the parameter $\lambda\to 0$ as $T\to T_{g}$. If the kagome-bilayer can be viewed as exhibiting a glass transition at $T=0$, as conventional 2D spin-glasses, the autocorrelation function might be dominated by a power-law behavior at intermediate time scale before it is cut off by the stretched exponential. Further larger scale simulations are required to investigate this scenario.

The power-law tail and the associated collective behaviors also depend strongly on the density $x$ of orphan spins. The effective interaction between two such defects separated by a distance $r$ is given by $J_{\rm eff}(r)\sim T\mathcal{J}(r/\xi(T))$, where $\xi\sim 1/\sqrt{T}$ is the temperature-dependent correlation length of the background spin liquid, and the function $\mathcal{J}(r)\sim\exp(-r)$ decays exponentially at large distances sen12 . The $T$-linear prefactor here indicates the entropic origin of the effective interaction, namely $J_{\rm eff}$ arises from conformational entropy of the fluctuating background correlated spins. Since the average distance between orphan spins scales as $\ell\sim 1/\sqrt{x}$, one thus obtains an average interaction $\overline{J}_{\rm eff}\sim J_{\rm eff}(\ell)\sim T\,\exp(-\sqrt{T/x})$, which becomes exponentially weak at small vacancy percentages. Despite this weakened interaction, the collective behavior of orphan spins would set in at a temperature that is of the order of the effective interaction. This gives the condition: $T^{*}\sim\,\overline{J}_{\rm eff}(T^{*})$. Using the expression for $\overline{J}_{\rm eff}$ above, one thus obtain a characteristic temperature $T^{*}\sim x$ that decreases linearly with the reduced defect density. Physically, the thermal correlation length at this $T^{*}$ is comparable to inter-orphan-spin distance.

VI Discussion and outlook

To summarize, we have extensively characterized the spin dynamics in the liquid phase of Heisenberg antiferromagnet on the kagome bilayer, which is relevant for the frustrated magnet SCGO. By computing the dynamical structure factor at different temperatures and dilutions, we show that the spin excitations are dominated by spin diffusion in the low energy, long time regime. The spin diffusion constant depends weakly on temperature, but decreases with dilution. Another interesting result is the half moon pattern of the dynamical structure factor with energy $\omega\gtrsim J$. Similar features have recently been observed in some pyrochlore compounds, it remains to be seen whether these remnants of the propagating spin waves can be observed in SCGO. Our simulations on diluted bilayer kagome shows that spin diffusion remains the dominant process in the presence of site disorder. This result further confirms, from the dynamical viewpoint, that site-disorder itself does not immediately cause glassy behaviors in the classical spin liquid, although the diffusion relaxation time becomes longer with increasing disorder. However, for disorder due to non-magnetic vacancies, the presence of so-called orphan spins results in an intriguing power-law tail in the spin-autocorrelation function. This power-law slow dynamics indicates that the system might be on the verge of a glass transition, which could be induced by other perturbations.

As discussed above, our work offers important benchmark for future dynamics studies of kagome bilayer that include other perturbations. Of particular interest are those perturbations that might transform the classical spin liquid into either the conventional spin glass or the more exotic spin jam. Indeed, since the diffusive spin dynamics in highly frustrated magnets is mainly driven by the zero-energy modes, one expects a diminishing diffusivity when the number of such zero modes is significantly reduced. For example, the entropic barrier in the coplanar phase of kagome reduces the continuous weather-van modes to discrete zero modes defined on closed loops. It has been proposed that the much slower relaxation of these discrete loops might give rise to glassiness without intrinsic disorder in kagome cepas12 ; cepas14 . However, the coplanar phase induced by thermal order-by-disorder seems to remain a classical spin liquid taillefumier14 . A transition into the glassy regime might still occur at a lower temperature when the dynamics is dominated by quantum tunneling of loops cepas12 .

Contrary to kagome Heisenberg antiferromagnets, there is no thermal induced coplanar or collinear phase in kagome bilayer. On the other hand, it has been proposed by one of us and co-authors in Ref. klich14 that a coplanar regime, in which spins in each tetrahedron are collinear, can be induced by quantum fluctuations. Moreover, different coplanar ground states can be mapped to discrete hexagonal tilings. Importantly, there is no continuous weather-van modes in this coplanar regime, and the only zero-energy modes are system-wide extended strings klich14 . As jamming transition often occurs in such constrained discrete models, the resultant coplanar phase is dubbed the spin jam klich14 ; yang15 . It is argued that quantum fluctuations transform the degenerate classical ground-state manifold into a rugged landscape that is different from that of conventional spin glass. While this spin-jam picture seem to explain some properties of SCGO and other similar glassy magnets, such as the much weaker memory effect samarakoon16 ; samarakoon17 , an important open question is to see how dynamical behaviors characteristic to spin-jam evolve from the classical spin liquid, which will be left for future study.

Acknowledgements.

We thank Anjana Samarakoon and Israel Klich for useful discussions and for collaborations on related projects. DZ and SHL were supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0016144. The author also acknowledge the support of Advanced Research Computing Services at the University of Virginia.