Quantum-critical scaling properties of the two-dimensional random-singlet state

Using quantum Monte Carlo simulations, we study effects of disorder on the $S=1/2$ Heisenberg model with exchange constant $J$ on the square lattice supplemented by multispin interactions $Q$. It was found recently [L. Liu et al., Phys. Rev. X 8, 041040 (2018)] that the ground state of this $J-Q$ model with random couplings undergoes a quantum phase transition from the Néel antiferromagnetic state into a randomness-induced spin-liquid-like state that is a close analog to the well known random-singlet (RS) state of the Heisenberg chain with random couplings. This 2D RS state arises from a spontaneously symmetry-broken fourfold degenerate columnar valence-bond solid that is broken up by the disorder into finite domains, with spinons localized at topological defects. The interacting spinons form a critical collective many-body state without magnetic long range order but with the mean spin-spin correlations decaying with distance $r$ as $r^{-2}$, as in the one-dimensional RS state. The dynamic exponent $z\ge 2$, varying continuously with the model parameters. In this work, we further investigate the properties of the RS state in the $J-Q$ model with random $Q$ couplings. We study the temperature dependence of the specific heat and various susceptibilities for large enough systems to reach the thermodynamic limit. We also analyze the size dependence of the critical magnetic order parameter and its susceptibility in the ground state. For all these quantities, we find consistency with the conventional quantum-critical scaling laws when the condition implied by the $r^{-2}$ form of the spin correlations is imposed. In particular, all the different quantities can be explained by the same value of the dynamic exponent $z$ at fixed model parameters. We argue that the RS state identified in the $J-Q$ model corresponds to a generic renormalization group fixed point that can be reached in many quantum magnets with random couplings and that it has already been observed experimentally.

Figure 1

Illustration of the terms in the 2D $J-Q$ model studied in this paper. Sites on the square lattice are shown as blue dots. The light blue bars represent the $J$ terms while the groups of three connected red bars represent the products of three singlet projectors constituting each $Q$ term.

Figure 2

Internal energy density vs temperature for system sizes $L=32$, 64, and 128, averaged over ${10}^3$ disorder realizations. The raw data are shown in (a) and the results normalized by $E\left(T_0\right)$, with $T_0=0.01$, are shown in (b). The insets show magnified graphs of the low-$T$ data. In (a) we do not show the error bars for $L=64$ in order to avoid clutter (the size of these errors is between those for $L=32$ and $L=124$). In the inset of (b), we have marked the lowest and highest points on the vertical axis by $y_1=0.99994$ and $y_2=1.00004$.

Figure 3

Normalized energy $E\left(T\right)/E\left(T_0\right)$ obtained by averaging over more than ${10}^4$ disorder realizations with $L=256$ and $T_0=1/150$. The blue curves show a fit of the data up to $T_{\max }=0.1$ to Eq. (6). The exponent $\alpha$ is graphed versus $T_{\max }$ in Fig. 4. Panel (a) shows the whole range of date and (b) focuses on the results for $T$ up to 0.05.

Figure 4

Exponent $\alpha =2/z+1$ in the power law fit, Eq. (6), of the normalized energy as a function of the highest temperature $T_{\max }$ included. The data are from Fig. 3. Based on the value of the reduced ${\chi }^2$ values, which are shown in the inset, the fits are statistically acceptable for $T_{\max }\le 0.1$. The result $\alpha =1.503\left(6\right)$ from the susceptibility fits in Sec. 4 is shown as the blue line.

Figure 5

Temperature dependence of the disorder-averaged uniform (a) and local (b) susceptibility for systems of size $L=256$. The red and blue curves are fits to the form $f\left(T\right)=a+b{T}^{-c}+d{T}^e$, where $c=1-2/z>0$ controls the leading (divergent) behavior. The values of the leading exponents are $c=0.494\left(9\right)$ for ${\chi }_{\mathrm{u}}$ and $c=0.500\left(8\right)$ for ${\chi }_{\mathrm{loc}}$. In the correction term the exponent $e$ is positive, with $e=0.97\left(6\right)$ for ${\chi }_{\mathrm{u}}$ and $e=0.60\left(2\right)$ for ${\chi }_{\mathrm{loc}}$. The black solid curves show the fitted functions without the second power laws, and the dashed curves show only the leading term $b{T}^{-c}$. The error bars on the QMC data are smaller than the graph symbols.

Figure 6

Dependence of the staggered susceptibility on the inverse temperature for system sizes $L=128$ and 256. The predicted linear form, Eq. (16), is fully consistent with the data starting from $\beta \approx 50$, as shown with the fitted line.

Figure 7

Log-log plot showing the size dependence of $S\left(\pi \right), {\chi }_s$, and ${\chi }_s/S\left(\pi \right)$ for system sizes up to $L=48$. The blue and black curves show fits to the form $Y=a{L}^z+b{L}^c$. The inset shows a log-linear plot of $S\left(\pi \right)$, making clear the logarithmic divergence of the form Eq. (18), shown by the fitted line.

Figure 8

Effective dynamic exponent defined using two system sizes, $L$ and $2L$, according to Eq. (20) for the quantities ${\chi }_{\mathrm{s}}$ (red symbols) and ${\chi }_{\mathrm{s}}/S\left(\pi \right)$ (blue symbols). The effective system size is defined as $L_{\mathrm{eff}}={\left(2L^2\right)}^{1/2}$. The solid lines show linear fits in $1/L_{\mathrm{eff}}$ giving the estimates of the extrapolated exponents $z=z_{\mathrm{eff}}(L\to \infty )=3.96\left(17\right)$ from ${\chi }_s$ and $z=3.94\left(14\right)$ from ${\chi }_s/S\left(\pi \right)$ the text. The dashed curves are obtained from the direct fits of the size dependent quantities in Fig. 7 to a form with two power laws, which were subsequently used in Eq. (20). The resulting extrapolated dynamic exponent with error bars are indicated close to $1/L_{\mathrm{eff}}=0$.