Effect of magnetoelastic coupling on spin-glass behavior in Heisenberg pyrochlore antiferromagnets with bond disorder

Motivated by puzzling aspects of spin-glass behavior reported in frustrated magnetic materials, we theoretically investigate effects of magnetoelastic coupling in geometrically frustrated classical spin models. In particular, we consider bond-disordered Heisenberg antiferromagnets on a pyrochlore lattice coupled to local lattice distortions. By integrating out the lattice degree of freedom, we derive an effective spin-only model, the bilinear-biquadratic model with bond disorder. The effective model is analyzed by classical Monte Carlo simulations using an extended loop algorithm. First, we discuss the phase diagrams in detail by showing the comprehensive Monte Carlo data for thermodynamic and magnetic properties. We show that the spin-glass transition temperature $T_{\mathrm{f}}$ is largely enhanced by the spin-lattice coupling $b$ in the weakly disordered regime. By considering the limit of strong spin-lattice coupling, this enhancement is ascribed to the suppression of thermal fluctuations in semidiscrete degenerate manifold formed in the presence of the spin-lattice coupling. We also find that, by increasing the strength of disorder $\Delta$, the system shows a concomitant transition of the nematic order and spin glass at a temperature determined by $b$, being almost independent of $\Delta$. This is due to the fact that the spin-glass transition is triggered by the spin collinearity developed by the nematic order. Although further-neighbor exchange interactions originating in the cooperative lattice distortions result in spin-lattice order in the weakly disordered regime, the concomitant transition remains robust with $T_{\mathrm{f}}$ almost independent of $\Delta$. We find that the magnetic susceptibility shows hysteresis between the field-cooled and zero-field-cooled data below $T_{\mathrm{f}}$, and that the nonlinear susceptibility shows a negative divergence at the transition. These features are common to conventional spin-glass systems. Meanwhile, we find that the specific heat exhibits a broad peak at $T_{\mathrm{f}}$, and that the Curie-Weiss temperature varies with $\Delta$, even in the region where $T_{\mathrm{f}}$ is insensitive to $\Delta$. In addition, we clarified that the concomitant transition remains robust against a substantial external magnetic field. These features are in clear contrast to the conventional spin-glass behavior. Furthermore, we show that the cubic susceptibility obeys a Curie-Weiss–type law and the estimated “Curie-Weiss” temperature gives a good measure of the spin-lattice coupling even in the presence of bond randomness. We also show, by studying single-spin-flip dynamics in the nematic phase, that the glassy spin dynamics may be observed at a rather high temperature in a realistic situation for weak disorder. All these results are discussed in comparison with experiments for typical pyrochlore magnets, such as ${\mathrm{Y}}_2{\mathrm{Mo}}_2{\mathrm{O}}_7$ and ${\mathrm{ZnCr}}_2{\mathrm{O}}_4$.

Figure 1

(a) 16-site cubic unit cell of the pyrochlore lattice. (b) Schematic illustration of the coupling of spins to a bond distortion ${\rho }_{ij}$. The antiferromagnetic exchange interactions are enhanced on shorter bonds. (c) Schematic illustration of a cooperative aspect of bond distortions; A shift of the B site while elongating (shortening) the AB (BC) bond enhances an antiferromagnetic spin correlation between the next-nearest-neighbor spins A and C.

Figure 2

Ground states of the model (5) for (a) $J_2^{\mathrm{coop}}=0$ and (b) $J_2^{\mathrm{coop}}>0$ ($b>0$ and $\Delta =0$). In (a), we show one of the macroscopically degenerate ground states (ice-rules configurations). The common axis of spins, which is denoted by a broken arrow $\vec{Q}$, is spontaneously selected below $T_{\mathrm{c}}$. In every tetrahedron, two of spins are aligned parallel to $\vec{Q}$ and the other two antiparallel to $\vec{Q}$. In (b), we present the $\vec{q}=0$ spin-lattice (Néel) order. The open and filled circles denote two nonequivalent sites with opposite spins.

Figure 3

(a) Schematic picture of a multiple valley structure in the spin-ice-type manifold. (b) Two different ice-rule states are shown. The hexagon with a bold dashed line denotes one of the shortest loops on which a flip of all spins transforms the ice-rule state to another ice-rule state. See the text for details.

Figure 4

Schematic phase diagrams of model (5) for three typical cases: (a) in the absence of the spin-lattice coupling $b=0$, (b) $b>0$ and $J_2^{\mathrm{coop}}=0$, and (c) $b>0$ and $J_2^{\mathrm{coop}}>0$. Our interest here is how the SG transition is induced by the bond disorder $\Delta$ in the competition with the nematic and Néel orderings, as shown in (b) and (c).

Figure 5

Flowchart of the MC simulation in this study. Each MC step consists of a lattice sweep by single-spin updates, loop flips, and a replica exchange between neighboring temperatures.

Figure 6

$\Delta -T$ phase diagrams obtained by MC simulation at $b=0.2$: (a) $J_2^{\mathrm{coop}}=0.0$ and (b) $J_2^{\mathrm{coop}}=0.075$. The nematic ($T_{\mathrm{c}}$), antiferromagnetic ($T_{\mathrm{N}}$), and SG transition temperatures ($T_{\mathrm{f}}$) are denoted by squares, triangles, and circles, respectively. In (a), $T_{\mathrm{f}}$ coincides with $T_{\mathrm{c}}$ for $\Delta \gtrsim 0.3$, suggesting a multicritical point at $\Delta \simeq 0.3(\simeq b)$. The cross in (a) denotes $T_{\mathrm{f}}$ for $b=0$ and $\Delta =0.1$ [20, 21]. See the text for details.

Figure 7

Schematic $\Delta -T$ phase diagram for $b>0$ and $J_2^{\mathrm{coop}}=0$ [see Fig. 6]. The nematic transition temperature $T_{\mathrm{c}}$ is almost independent of $\Delta$ as $T_{\mathrm{c}}\simeq b$. On the other hand, the SG transition temperature $T_{\mathrm{f}}$ increases linearly with $\Delta$ in the small-$\Delta$ region. $T_{\mathrm{f}}$ is enhanced by $b$ compared to that in the bilinear limit ($b=0$) denoted by the dotted line. The nematic and SG transitions merge into a concomitant transition for $\Delta \gtrsim b$.

Figure 8

The specific heat $C$, spin collinearity $Q^2$, and SG susceptibility ${\chi }_{\mathrm{SG}}$ calculated at $b=0.2$ for (a) $\Delta =0.0$, (b) $0.1$, (c) $0.2$, (d) $0.3$, (e) $0.5$, and (f) $0.8$. The data are calculated for the system sizes ranging from $L=2$ (128 spins) to $L=8$ (8192 spins).

Figure 9

(a) $L$ dependence of the peak values of the specific heat $C$, $C_{\mathrm{peak}}$. The lines are power-law fitting by $C_{\mathrm{peak}}\propto L^p$. The obtained values of $p$ are shown in Table 1. (b) System-size dependencies of the peak temperature of $C$. The lines represent the extrapolations of $C_{\mathrm{peak}}$ to the bulk limit with $T_{\mathrm{c}}\left(L\right)-T_{\mathrm{c}}\propto 1/L^3$. The obtained values of $T_{\mathrm{c}}$ are summarized in Table 1.

Figure 10

Scaling collapse of ${\chi }_{\mathrm{SG}}$ calculated with $J_2^{\mathrm{coop}}=0$ at (a) $\Delta =0.1$, (b) 0.2, (c) 0.3, (d) 0.5, and (e) 0.8. The estimated $T_{\mathrm{f}}$ and the critical exponents are shown in Table 1.

Figure 11

Scaling collapse of ${\chi }_{\mathrm{SG}}$ calculated for the Ising model given in Eq. (21). We obtained $T_{\mathrm{f}}=0.151\left(2\right)$, $\gamma =0.79\left(2\right)$, and $\nu =0.86\left(2\right)$ for $L=2,3$, and 4.

Figure 12

Finite-size scaling analysis of the nematic transition at $\Delta =0.8$ and $J_2^{\mathrm{coop}}=0$. (a) Scaling collapse of ${\chi }_Q$ for the data in the range of $0.225\le T\le 0.275$. We obtained $T_{\mathrm{f}}=0.248\left(2\right)$, ${\gamma }_Q=1.5\left(1\right)$, and ${\nu }_Q=0.80\left(2\right)$. (b) Semilogarithmic plot of the $L$ dependence of $C_{\mathrm{peak}}$. (c) $1/C_{\mathrm{peak}}$ as a function of $1/L$.

Figure 13

The specific heat $C$, spin collinearity $Q^2$, SG susceptibility ${\chi }_{\mathrm{SG}}$, and the square of sublattice magnetization $m_{\mathrm{s}}^2$, calculated at $b=0.2$ and $J_2^{\mathrm{coop}}=0.075$: (a) $\Delta =0.0$, (b) 0.2, (c) 0.5, and (d) 0.8.

Figure 14

Extrapolation of the peak temperatures of the specific heat $C$ to the bulk limit. The data are calculated with $b=0.2$ and $J_2^{\mathrm{coop}}=0.075$.

Figure 15

Squared sublattice magnetization $m_{\mathrm{s}}^2$ as a function of $\Delta$ at (a) $T=0.25$, (b) 0.2, and (c) 0.1. We take $b=0.2$ and $J_2^{\mathrm{coop}}=0.075$.

Figure 16

Scaling collapses of the SG susceptibility ${\chi }_{\mathrm{SG}}$ and nematic susceptibility ${\chi }_Q$ at (a) $\Delta =0.5$ and (b) 0.8 for $J_2^{\mathrm{coop}}=0.075$ and $b=0.2$ in the plateau regime [see Fig. 6]. We take $b=0.2$ and $J_2^{\mathrm{coop}}=0.075$. The estimated values for the transition temperatures and the critical exponents are presented in Table 3.

Figure 17

$T$ dependence of the inverse of the magnetic susceptibility $\chi$ calculated at $b=0.2$ and $J_2^{\mathrm{coop}}=0$ for $L=3$. The data are plotted for $0.0\le T\le 1.5$ and $0.2\le T\le 1.0$ in (a) and (b), respectively. The lines in (a) and (b) denote the fits by the Curie-Weiss law in Eq. (23) in the range of $1.2\le T\le 1.5$ and $0.6\le T\le 0.9$, respectively.

Figure 18

Field-cooled (FC, filled squares) and zero-field-cooled (ZFC, open circles) susceptibilities calculated at $b=0.2$ and $J_2^{\mathrm{coop}}=0$. The data are calculated at $H=0.1$ for $L=3$.

Figure 19

$T$ dependence of the inverse of the cubic susceptibility ${\overline{\chi }}_3$ calculated at $J_2^{\mathrm{coop}}=0$. The system size is $L=1$ (128 spins). (a) The data for different $b$ at $\Delta =0$. The straight lines denote the fits by high-$T$ asymptotic behavior in Eq. (25). The inset shows the estimated ${\theta }_3$ is shown as a function of $b$. (b) The data for different $\Delta$ at $b=0.2$. The straight lines denote the fits by Eq. (25). The inset shows the $\Delta$ dependence of the estimated ${\theta }_3$.

Figure 20

$T$ dependence of the nonlinear susceptibility ${\chi }_3$ calculated for $b=0.2$ and $J_2^{\mathrm{coop}}=0$ at (a) $\Delta =0.0$, (b) $\Delta =0.1$, and (c) $\Delta =0.8$. The vertical dashed, dotted, solid lines denote the nematic transition temperature, the SG transition temperature, and the transition temperature of the concomitant transition, respectively (see Table 1).

Figure 21

Scaling collapse of the nonlinear susceptibility ${\chi }_3$ at $\Delta =0.8$. The data are taken from Fig. 20. We obtained the best fit with $\gamma =2.8\left(5\right)$ and $\nu =1.7\left(3\right)$.

Figure 22

(a) $H$ dependence of the specific heat $C$ at $\Delta =0.5$. (b) $H$ dependence of the ZFC (open circles) and FC (filled squares) susceptibilities at $\Delta =0.8$. See also Fig. 18. The data are calculated at $J_2^{\mathrm{coop}}=0$ for the system size $L=3$.

Figure 23

Autocorrelation functions calculated with $\Delta =0.1$, $b=0.2$, and $J_2^{\mathrm{coop}}=0.0$ (linear regime). The MC dynamics exhibits a dynamical freezing below $T_{\mathrm{c}}\simeq b$.

Figure 24

Comparison of thermalization processes of $q_{\mathrm{EA}}^2$ from a disordered configuration with and without using the extended loop update. The data are taken at $\Delta =0.1$ and $T=0.08$ for the model (5) with $b=0.2$ and $J_2^{\mathrm{coop}}=0$ in the system size $L=2$. The relaxation is remarkably accelerated by the loop update.