Randomness-induced quantum spin liquid behavior in the s=1/2 random J1-J2 Heisenberg antiferromagnet on the honeycomb lattice

TL;DR

This study shows that randomness can induce a gapless quantum spin liquid-like state in a frustrated honeycomb lattice Heisenberg model, suggesting such states may be common in real materials with inherent disorder.

Contribution

The paper demonstrates that randomness stabilizes a gapless QSL-like state over a wide parameter range in the honeycomb J1-J2 Heisenberg model, highlighting the role of disorder in quantum spin liquids.

Findings

Randomness induces a gapless QSL-like state in the model.

The random-singlet state is robust across various parameters.

Implications for experimental materials exhibiting QSL-like behaviors.

Abstract

We investigate the ground-state and finite-temperature properties of the bond-random $s=1/2$ Heisenberg model on a honeycomb lattice with frustrated nearest- and next-nearest-neighbor antiferromagnetic interactions, $J_1$ and $J_2$, by the exact diagonalization and the Hams--de Raedt methods. The ground-state phase diagram of the model is constructed in the randomness versus the frustration ($J_2/J_1$) plane, with the aim of clarifying the effects of randomness and frustration in stabilizing a variety of distinct phases. We find that the randomness induces the gapless quantum spin liquid (QSL)-like state, the random-singlet state, in a wide range of parameter space. The observed robustness of the random-singlet state suggests that the gapless QSL-like behaviors might be realized in a wide class of frustrated quantum magnets possessing a certain amount of randomness or inhomogeneity, without fine-tuning the interaction parameters. Possible implications to recent experiments on the honeycomb-lattice magnets Ba$_3$CuSb$_2$O$_9$ and 6HB-Ba$_3$NiSb$_2$O$_9$ exhibiting the gapless QSL-like behaviors are discussed.

Equations

$$\mathcal{H}=J_1\sum _{⟨i,j⟩}{j}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j+J_2\sum _{⟨⟨i,j⟩⟩}{j}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$S_{\mathbf{q}}$$

$$=\frac{1}{N}{\left[\sum _{i,j}⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩\cos \left(\mathbf{q}\cdot \left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\right)\right]}_J,$$

$$m_{AF}^2$$

$$=\frac{8}{N(N+4)}{\left[\sum _{\alpha }\sum _{i,j\in \alpha }⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩\right]}_J,$$

$$m^2=m_{\infty }^2+\frac{c_1}{\sqrt{N}},$$

$$m_{str}^2=\frac{8}{3N(N+4)}{\left[\sum _{\nu }\sum _{{\alpha }_{\nu }}\sum _{i,j\in {\alpha }_{\nu }}⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩\right]}_J,$$

$$\bar{q}=\frac{1}{N}\sqrt{{\left[\sum _{i,j}{⟨{\mathbf{S}}_i\cdot {\mathbf{S}}_j⟩}^2\right]}_J}.$$

$$S_{\mathbf{q}}(\omega )$$

$$=-\underset{\eta \to 0}{\lim }{\left[\frac{1}{\pi }\mathrm{Im}⟨(S_{\mathbf{q}}^z{)}^{†}\frac{1}{\omega +E_0+i\eta -\mathcal{H}}{S}_{\mathbf{q}}^z⟩\right]}_J,$$

Full text

I Introduction

The quantum spin liquid (QSL) state without any spontaneously broken Hamiltonian symmetry, which accompanies no magnetic long-range order (LRO) down to low temperatures, has long received much attention.QSL For the realization of such QSL state, geometrical frustration is considered to be essential, and frustrated magnets have been the main target of the quest for QSL materials. In particular, the $s=1/2$ organic triangular-lattice salts $\kappa$-(ET)2Cu2(CN)3,ETsalt ; ETsaltCv ; ET-Matsuda ; ET-ChargeGlass ; ET-Sasaki EtMe3Sb[Pd(dmit)2]2, dmitsalt ; dmitsalt-Matsuda ; dmitsaltCv ; dmit-ChargeGlass and more recently $\kappa$-H3(Cat-EDT-TTF)2 Isono1 ; Isono2 ; Ueda were reported to exhibit the QSL-like behaviors down to very low temperatures. The QSL states of these organic salts commonly exhibit gapless (or nearly gapless) behaviors characterized by, e.g., the low-temperature specific heat linear (or almost linear) in the absolute temperature $T$. ETsaltCv ; ET-Matsuda ; dmitsalt-Matsuda ; dmitsaltCv Another well-studied candidate of the QSL might be the $s=1/2$ kagome-lattice inorganic compound herbertsmithite ZnCu3(OH)6Cl2. Shores ; Helton ; Olariu ; Freedman ; Han ; Imai This kagome material was also reported to exhibit gapless QSL-like behaviors, Helton ; Olariu ; Han while a recent NMR study showed a nonzero spin gap. Imai

Despite such recent experimental progress, the true origin of the experimentally observed QSL-like behaviors still remains not fully understood and is under hot debate. In many theoretical studies, it has been assumed that the system is sufficiently clean so that the possible effect of randomness or inhomogeneity is negligible and unimportant. Meanwhile, one of the present authors (H.K.) and collaborators have claimed that the QSL-like behaviors recently observed in triangular-lattice organic salts and kagome-lattice herbertsmithite might be the randomness-induced one, the random-singlet state. Watanabe ; Kawamura ; Shimokawa The advocated random-singlet state is a gapless QSL-like state where spin singlets of varying strengths are formed in a hierarchical manner on the background of randomly distributed exchange interactions $J_{ij}$. The state may also be regarded as a sort of “Anderson-localized resonating valence bond (RVB) state”. Indeed, it was demonstrated that the random-singlet state exhibited the $T$-linear specific heat and the gapless susceptibility with an intrinsic Curie tail, accompanied by the gapless and broad features of the dynamical spin structure factor. Watanabe ; Kawamura ; Shimokawa

The source of the randomness or inhomogeneity in real materials could be various. For triangular organic salts, in Ref. Watanabe , it was argued that the effective randomness or inhomogeneity might be self-generated for the spin degrees of freedom via the coupling between the spin and the charge, owing to the charge order inherent to these compounds and the associated slowing down of the electric-polarization degrees of freedom at each dimer molecule. Indeed, there reported some experimental evidence of such spatial inhomogeneity in the charge and spin distributions in triangular organic salts. ET-Sasaki ; Kawamoto ; Shimizu ; Nakajima For the kagome herbertsmithite, in Ref. Kawamura , it was suggested that the random Jahn–Teller (JT) distortion of the [Cu(OH)6]4- octahedra driven by the random substitution of magnetic Cu2+ for nonmagnetic ${\rm Zn^{2+}}$ on the adjacent triangular layer Freedman might give rise to the random modification of the exchange paths between various neighboring magnetic Cu2+ pairs on the kagome layer, leading to the random modulation of the exchange couplings between various neighboring Cu2+ pairs on the kagome layer.

The possible important role played by the randomness or inhomogeneity has also been reported in other triangular and kagome magnets as well. One example might be an inorganic triangular antiferromagnet Cs2Cu(Br1-xClx)4, a mixed crystal of Cs2CuBr4 and Cs2CuCl4. This random magnet was observed not to exhibit the magnetic LRO nor the spin-glass freezing down to very low temperatures in a certain range of $x$. Ono Another intriguing example might be the triangular-lattice organic salt $\kappa$-(ET)2Cu[N(CN)2]Cl, which exhibits the standard AF LRO, but exhibits the QSL-like behaviors when the randomness is artificially introduced by X-ray irradiation. Furukawa Another example might be the $s=1/2$ kagome-like Heisenberg AF, ZnCu3(OH)6SO4, where the spin-1/2 Cu2+ is located on the corrugated kagome plane. This kagome material contains a significant amount of randomness and was found to exhibit gapless QSL-like behaviors characterized by the $T$-linear specific heat. Li ; Zorko ; Zorko2 Thus, there now emerge growing experimental lines of evidence that the randomness-induced QSL-like state exists in nature.

From a theoretical viewpoint, one might naturally ask whether the combined effect of randomness and quantum fluctuations is sufficient to cause the QSL-like behaviors. A counter-example of this was reported for the case of the random-bond $s=1/2$ AF Heisenberg model on the square lattice. In contrast to the frustrated triangular- or kagome-lattice counterparts obeying the same form of randomness distribution, this unfrustrated model persists to exhibit the AF LRO up to the maximal randomness. RandomSquare1 ; Watanabe This observation certainly suggests that the frustration also plays an important role in realizing the QSL-like behaviors.

Under such circumstances, it is important to clarify the role of frustration along with that of randomness in realizing the random-singlet state in quantum magnets. In this paper, we wish to undertake such a study by investigating the properties of a model where the strengths of both randomness and frustration can independently be tuned. For this purpose, we choose the $s=1/2$ honeycomb-lattice Heisenberg model with the competing AF nearest- and next-nearest-neighbor interactions $J_{1}$ and $J_{2}$, as shown in Fig. 1. The honeycomb lattice is bipartite so that the $J_{1}$-only model is unfrustrated. Frustration is introduced via the competition between $J_{1}$ and $J_{2}$, the ratio $J_{2}/J_{1}$ representing the extent of frustration. The choice of the honeycomb lattice is motivated by the fact that the honeycomb lattice has only three nearest neighbors, a minimal number among various 2D lattices, and is subject to the enhanced effects of fluctuations possibly destroying the magnetic LRO.

In fact, the low-temperature properties of the regular $s=1/2$ $J_{1}$-$J_{2}$ Heisenberg model on the honeycomb lattice have attracted much interest. Fouet ; Mulder ; Mosadeq ; Lauchli ; White ; Ganesh ; Gong In particular, the ground-state phase diagram of the model was extensively studied as a function of $J_{2}$ by various methods, including the exact diagonalization (ED) Fouet ; Mosadeq ; Lauchli and the density matrix renormalization group (DMRG) White ; Ganesh ; Gong methods. When $J_{2}$ is small, the ground state is the standard two-sublattice AF state as illustrated in Fig. 2(a). When the strength of frustration exceeds the critical value $J_{2c}$, the AF state is destabilized and QSL-like states emerge. Fouet ; Mosadeq ; Lauchli ; White ; Ganesh ; Gong In fact, the critical $J_{2}$ value of the AF instability is around $0.2-0.25$, which is greater than the classical counterpart $J_{2c}=1/6$. Beyond $J_{2c}$, several types of nonmagnetic states have been discussed in the literature, some of which are illustrated in Fig. 3. They include the columnar-dimer state [Fig. 3(a)], which breaks both the lattice translational and the lattice rotational symmetries, the staggered-dimer state [Fig. 3(b)], which keeps the lattice translational symmetry but breaks the lattice rotational symmetry,Fouet ; Mosadeq ; White ; Ganesh and the plaquette state [Fig. 3(c)], which keeps the lattice rotational symmetry but breaks the lattice translational symmetry. Fouet ; Mosadeq ; Lauchli ; White ; Ganesh ; Gong These states are all variants of the valence-bond crystal (VBC) state with a nonzero spin gap. In addition, the possible occurrence of the $Z_{2}$ spin-liquid state was also discussed (see, e.g., Ref. Gong and references cited therein). The magnetically ordered state such as the stripe-ordered state illustrated in Fig. 2(b) was also invoked as a ground state of the larger $J_{2}$ region. Lauchli Our present study addresses the issue of what happens to these magnetic and nonmagnetic states of the regular model if one introduces the randomness.

Experimentally, the magnetic ordering of honeycomb-lattice AFs has attracted much recent interest. Owing to the bipartite nature of the lattice, some of them exhibit the standard AF order,Cu5SbO6 whereas some are reported not to exhibit the magnetic order.Azuma ; Matsuda ; Zhou ; Nakatsuji ; Mendels-Cu ; Ishiguro ; Do ; Hagiwara ; Katayama ; Cheng ; Darie ; Mendels-Ni Even if one sets aside the dimerized nonmagnetic state with a finite spin gap, occasionally induced by the uniform JT structural distortion of the lattice, interesting spin-liquid-like behaviors have been reported in several honeycomb magnets.

One example might be Bi3Mn4O12(NO3) (BMNO), an insulating $s=3/2$ system. This material remains paramagnetic down to low temperatures in zero field with a small spin correlation length, but exhibits the AF order upon applying even weak fields. Azuma ; Matsuda Since the Mn4+ spin is 3/2 here, the spin-liquid-like behaviors of BMNO might be expected to be essentially of classical origin. Indeed, the spin-liquid-like behavior and the field-induced AF of this material have successfully been explained within the classical (or semi-classical) picture as the “ring liquid” and the “pancake liquid” arising from the characteristic ring-like degeneracy in the wavevector space. Okumura The manner by which this degeneracy is lifted by fluctuations via the ‘order-from-disorder’ mechanism leading to the magnetic ordering was also studied. Okumura ; Mulder In particular, in Ref. Okumura , it was argued that, close to the AF instability, which was equal to $J_{2c}=1/6$ in the classical limit, the energy scale of this order from disorder could be very small such that the classical spin-liquid state (the ring liquid or the pancake liquid) might be stabilized down to very low temperatures.

A different type of spin-liquid-like behavior, likely to be of quantum origin, has recently been reported in the honeycomb-lattice-based compounds Ba3CuSb2O9, Zhou ; Nakatsuji ; Mendels-Cu ; Ishiguro ; Do ; Hagiwara ; Katayama and 6HB-Ba3NiSb2O9. Cheng ; Darie ; Mendels-Ni The former Cu2+ compound Ba3CuSb2O9 has $s=1/2$ forming a decorated honeycomb lattice, with a significant amount of randomness associated with the Cu-Sb ‘dumbbell’ orientation. Nakatsuji Stoichiometric samples tend to be hexagonal, whereas off-stoichiometric samples tend to be orthorhombic with a static JT transition occurring at a higher temperature $T\simeq 200$K. Nakatsuji ; Hagiwara Although the orthorhombic sample was reported to exhibit the spin freezing at a low temperature $T_{g}\simeq 110$ mK, the hexagonal sample does not order magnetically down to very low temperatures. Nakatsuji For such hexagonal samples, in Ref. Nakatsuji , two different scenarios were suggested. In one, the spin and orbital degrees of freedom are entangled and fluctuate together down to low temperatures forming a sort of spin-orbital liquid (the dynamical JT effect), Nakatsuji ; Ishiguro ; Do ; Hagiwara ; Katayama while, in the other, the static random JT distortion might lead to the ‘random-singlet state’. Mendels-Cu Many of the recent studies seem to point to the former scenario, while the NMR study of Ref. Mendels-Cu favored the second scenario, i.e., the ‘random-singlet state’ characterized by the two kinds of gaps, $\Delta_{1}$ and $\Delta_{2}$, $\Delta_{1}$ ($>\Delta_{2}$) being in the order of the main exchange interaction and $\Delta_{2}$ induced by applied fields. For the orthorhombic samples, the lattice symmetry is lowered from hexagonal to orthorhombic even macroscopically at temperatures below the static hexagonal-to-orthorhombic JT transition at 200 K, which might affect the underlying low-temperature physics. Indeed, the spin-glass freezing was reported at a low temperature of 110 mK in the orthorhombic sample. Nakatsuji

Hence, the situation for Ba3CuSb2O9 still remains not totally clear, including the question on how the random-singlet state identified for the bond-random triangular and kagome models, Watanabe ; Kawamura ; Shimokawa i.e., the randomness-induced gapless QSL-like state, is related or unrelated to the ‘random-singlet state’ discussed in the literature for the honeycomb-lattice-based AFs Ba3CuSb2O9 both in the hexagonal and orthorhombic samples.

By contrast, the latter Ni2+ compound 6HB-Ba3NiSb2O9 synthesized in high fields and at high temperatures has $s=1$. Cheng ; Darie ; Mendels-Ni Recent analysis has revealed that the structure of this compound is likely to be trigonal rather than hexagonal, forming a honeycomb lattice made up of triangular bilayers with a significant amount of randomness contained, where the nearest- and the next-nearest-neighbor interactions $J_{1}$ and $J_{2}$ arise from the alternating arrangement of the Ni-Sb dumbbell on adjacent triangular layers. Darie ; Mendels-Ni In this Ni2+ compound, the orbital degrees of freedom are absent unlike the Cu2+ counterpart, yet the system exhibits the QSL-like behavior characterized by the $T$-linear specific heat. One might wonder if the QSL-like behavior of this material might have any connection with the randomness-induced gapless QSL-like states of triangular and kagome AFs. Watanabe ; Kawamura ; Shimokawa This possible connection provides another motivation of our present study.

The structure of this article is as follows. In sect. II, we introduce our model on the honeycomb lattice, the random-bond $s=1/2$ $J_{1}$-$J_{2}$ Heisenberg model, and explain the computational method employed. In sect. III, we summarize the ground-state properties of the corresponding regular model on the basis of the previous works and some new data of our own. This sect. serves the basis of our study on the random model in the following sect. Sections IV and V are the main parts of the present paper. The ground-state properties of the random model are studied in sect. IV. The ground-state phase diagram is constructed in the frustration ($J_{2}/J_{1}$) versus the randomness (the parameter $\Delta$ is defined below) plane, and the properties of each phase are clarified. The finite-temperature properties of the model are studied in sect. V. Section VI is devoted to summary and discussion.

II Model and method

We consider the bond-random $s=1/2$ isotropic Heisenberg model on the honeycomb lattice with the AF nearest-neighbor and next-nearest-neighbor interactions $J_{1}$ and $J_{2}$. The Hamiltonian is given by

[TABLE]

where ${\bm{S}}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ is an $s=1/2$ spin operator at the $i$-th site on the honeycomb lattice, the sums $\Braket{i,j}$ and $\Braket{\Braket{i,j}}$ are taken over all nearest-neighbor and next-nearest-neighbor pairs on the lattice where periodic boundary conditions are assumed, while $j_{ij}\geq 0$ is the random variable obeying the bond-independent uniform distribution between $[1-\Delta,1+\Delta]$ with $0\leq\Delta\leq 1$. Hereafter, we put $J_{1}=1$ and $J_{2}/J_{1}=J_{2}>0$. Then, the parameter $J_{2}$ represents the degree of frustration, because the frustration in the present model is exclusively borne by the competition between the nearest- and next-nearest-neighbor interactions. Our present choice of the bond-independent uniform distribution for $j_{ij}$ is just for simplicity, whereas, in real materials, the distribution could be more complex and correlated. The parameter $\Delta$ represents the extent of the randomness: $\Delta=0$ corresponds to the regular case and $\Delta=1$ to the maximally random case. Again, the extent of the randmness $\Delta$ is taken to be common between $J_{1}$ and $J_{2}$ just for simplicity. By tuning the parameters $\Delta$ and $J_{2}$, we can control the degrees of both the randomness and the frustration independently.

The ground-state properties of the model are computed by the ED Lanczos method. We treat finite-size clusters with the total number of spins $N$ up to $N\leq 32$ (all even-$N$ samples with $8\leq N\leq 32$), with periodic boundary conditions being applied. All clusters studied are commensurate with the two-sublattice AF order illustrated in Fig. 2(a), and with the staggered-dimer order illustrated in Fig. 3(b). The clusters of $N=8$, 12, 16, 20, 24, 28, and 32 are commensurate with the stripe order of Fig. 2(b), among which $N=8$, 24, and 32 possess the threefold rotational symmetry of the bulk honeycomb lattice. (More generally, the clusters of $N=8$, 14, 18, 24, 26, and 32 possess the threefold rotational symmetry of the bulk honeycomb lattice.) The clusters of $N=12$, 18, 24, and 30 are commensurate with the columnar-dimer and the plaquette orders shown respectively in [Figs. 3(a) and 3(c)].

The numbers of independent bond realizations, $N_{s}$, used in the configurational or sample average are $N_{s}=100$, 50, 25, 16, and 10 for $N=8-24$, 26, 28, 30, and 32 for the order parameter, the spin gap, and the static spin structure factor, whereas $N_{s}=100$, 100, and 25 for $N=18,24$, and 32 for the dynamical spin structure factor, respectively. Error bars are estimated from sample-to-sample fluctuations.

The finite-temperature properties are computed by the Hams–de Raedt method, HamsRaedt where the thermal average is replaced by the average over a few pure states produced via the imaginary time-evolution of initial vectors. The method enables us to calculate various finite-temperature properties at nearly the same computational cost as the Lanczos method. Our finite-temperature computation is performed for the size $N=24$, where the averaging is made over 30 initial vectors and 50 independent bond realizations. Error bars of physical quantities are estimated from the scattering over both samples and initial states.

III Ground-state properties of regular model

In this section, we present the ground-state properties of the regular $(\Delta=0)$ $J_{1}$-$J_{2}$ model on the honeycomb lattice. There already exist some ED works on this model, Fouet ; Mosadeq ; Lauchli even to sizes larger than our largest one $N=32$. Even so, we feel that presenting some of our data in the form appropriate for later comparison with the random model might be useful.

III.1 Possible phases and phase diagram

First, we summarize in Fig. 4 the ground-state phase diagram of the regular $J_{1}-J_{2}$ honeycomb model in the region of $0\leq J_{2}\leq 0.5$, obtained on the basis of previous calculations and our present ED calculations. Three distinct phases appear, i.e., the two-sublattice AF phase, and the nonmagnetic gapped I and II phases. For the smaller $J_{2}$ region of $J_{2}\lesssim 0.25$, the standard two-sublattice AF state is stabilized owing to the bipartite nature of the honeycomb lattice. When $J_{2}$ is increased beyond a critical value of $J_{2c}\simeq 0.25$, the nonmagnetic gapped state, the gapped I state, appears. The state is magnetically disordered with a finite spin gap. As a candidate of this gapped I phase, many previous studies have pointed to the plaquette phase, although no complete consensus has emerged. Fouet ; Mosadeq ; Lauchli The plaquette state maintains the threefold lattice rotational symmetry but breaks the translational symmetry as can be seen from Fig. 3(c). Another candidate of the gapped I phase might be the $Z_{2}$ spin liquid phase. Gong Most studies reported the critical $J_{2}$ value of around $0.2-0.25$. Mosadeq ; Lauchli ; White ; Ganesh When $J_{2}$ is further increased beyond $J_{2}=0.35$-0.4, there occurs a transition from the gapped I phase to another gapped phase (gapped II phase), where the state is still magnetically disordered with a finite spin gap. As a candidate of this gapped II phase, many previous studies have pointed to the staggered-dimer state, Fouet ; Mosadeq ; White ; Ganesh also called the lattice nematic state.Mulder The staggered-dimer state maintains the translational symmetry, but breaks the threefold rotational symmetry as can be seen from Fig. 3(b). Another candidate of the gapped II phase might be the magnetic stripe phase shown in Fig. 2(b). Lauchli

We perform some ED calculations to obtain additional information on the ground-state properties of the regular model, and obtain the results basically supporting the phase diagram of Fig. 4. In Fig. 5, we show the computed ground-state spin structure factor $S_{\bm{q}}$ defined by

[TABLE]

where ${\bm{S}}_{\bm{q}}=\sum_{j}{\bm{S}}_{j}e^{i{\bm{q}}\cdot{\bm{r}}_{j}}$ is the Fourier transform of the spin operator, ${\bm{r}}_{j}$ is the position vector at the site $j$, ${\bm{q}}$ is the wavevector, while $\langle\cdots\rangle$ and $[\cdots]_{J}$ represent the ground-state expectation value (or the thermal average at finite temperatures) and the configurational average over $J_{ij}$ realizations. The $J_{2}$ values are (a) $J_{2}=0.15$, (b) 0.3, and (c) 0.5.

For the $J_{2}$-values smaller than $J_{2c}\simeq 0.25$, a rather sharp peak corresponding to the standard two-sublattice AF order appears at the $\Gamma$ point located at ${\bm{q}}=\pi(\frac{2}{\sqrt{3}},\frac{4}{3})$, where the length unit is taken here to be the nearest-neighbor distance of the honeycomb lattice. An example for the case of $J_{2}=0.15$ is shown in Fig. 5(a). As $J_{2}$ is increased to $J_{2}=0.3$, the AF peak at the $\Gamma$ point becomess broadened, but without any other peaks appearing. The $J_{2}=0.3$ point lies in the gapped I phase in the phase diagram of Fig. 4 so that the broad peak observed here is likely to be associated with the magnetic short-range order (SRO), which will be confirmed by our data of the AF order parameter shown below. Interestingly, at a still larger value of $J_{2}=0.5$ shown in Fig. 5(c) lying in the gapped II phase in the phase diagram of Fig. 4, the peak remains broad, but its position moves from the $\Gamma$ point to the M points located at ${\bm{q}}=\pi(\frac{2}{\sqrt{3}},0)$, $\pi(\frac{1}{\sqrt{3}},3)$, and $\pi(\frac{1}{\sqrt{3}},-3)$. From the lattice symmetry, there exist three equivalent but independent $M$ points, which, in terms of the ordering pattern, corresponds to the three distinct types of the stripe order Lauchli illustrated in Fig. 2(b). The observed broadness of the M-point peaks indicates that the stripe order remains to be a SRO here, not a true LRO, which will be confirmed also by our data of the stripe order parameter shown in Fig. 6(b).

In any case, the observation that the $S_{\bm{q}}$ peaks appear at mutually distinct positions at $J_{2}=0.3$ and 0.5 strongly suggests that the phases at $J_{2}=0.3$ and at $0.5$ are indeed different, the gapped I and II phases, and there is a phase transition between them.

In Fig. 6(a), we show the AF order parameter, the squared sublattice magnetization $m_{AF}^{2}$ associated with the standard two-sublattice order (the $\Gamma$-point order), plotted versus $1/\sqrt{N}$ for various values of $J_{2}$. It is defined by

[TABLE]

where $\alpha=A,B$ denotes the two triangular sublattices of the original honeycomb lattice, and the sum over $i,j\in\alpha$ is taken over all sites $i,j$ belonging to the sublattice $\alpha$. When the system retains a relevant magnetic long-range order $m_{\infty}$, the size dependence proportional to $1/\sqrt{N}$ is expected from the spin-wave analysis, i.e.,

[TABLE]

where $c_{1}$ is a constant. As can be seen from the figure, a linear extrapolation of our finite-$N$ data indicates that the standard AF LRO vanishes at around $J_{2}=0.25$, consistently with the phase diagram shown in Fig. 4. Lauchli

We also compute the magnetic order parameter associated with the stripe order (the M-point order) $m_{str}^{2}$, defined by

[TABLE]

where $\nu=1,2,$ and 3 refer to the three distinct types of stripe order associated with the threefold rotational symmetry of the lattice. The result is shown in Fig. 6(b) as a function of $1/\sqrt{N}$ for $J_{2}=0.3,0.4$, and 0.5, where we restrict the sizes only to $N=8$, 24, and 32, which retain the threefold lattice rotational symmetry being commensurate with the stripe order. As can be seen from the figure, $m_{str}^{2}$ is always extrapolated to negative values for $N\rightarrow\infty$, which indicates that the stripe order remains a SRO in both the gapped I and II phases.

We also compute the spin-gap energy, i.e., the energy difference $\Delta E$ between the ground state and the triplet first-excited state. It turns out that $\Delta E$ is extrapolated to zero for the AF state of $J_{2}\lesssim 0.25$, while to a positive nonzero value for the gapped I and II states of $J_{2}\gtrsim 0.25$. Some of the data will be shown in sect. 4 below (see Fig. 12).

IV Random model: Phase diagram and ground-state properties

IV.1 Phase diagram

In this section, we present our numerical results on the random $(\Delta\not=0)$ model in the region of $0\leq J_{2}\leq 0.5$. We first show our main result in Fig. 7, the ground-state phase diagram of the bond-random $s=1/2$ $J_{1}-J_{2}$ Heisenberg model on the honeycomb lattice in the frustration ($J_{2}$) versus the randomness ($\Delta$) plane. The $\Delta=0$ line corresponds to the phase diagram of the regular model shown in Fig. 4 in the previous section. In fact, when the randomness $\Delta$ is sufficiently weak, the phase diagram turns out to be qualitatively unchanged from that of the regular model.

In Fig. 7, four distinct phases are identified. Three of them, i.e., the AF phase, and the gapped phases I and II, have already been identified in the regular model, while, when the strength of the randomness exceeds a critical value $\Delta_{c}(J_{2})$, the fourth phase, the random-singlet phase, appears. In fact, the random-singlet phase turns out to be stabilized for a wide range of parameter space for moderate or strong randomness. As will be shown shortly, this phase is basically of the same nature as the one recently identified in the $s=1/2$ random triangular and kagome Heisenberg models. Watanabe ; Kawamura ; Shimokawa

Below, we show our numerical data for various observables including the order parameters, the spin-gap energy, and the static and dynamical structure factors for the $J_{2}$-regions of (i) $0\leq J_{2}\leq 0.25$ and (ii) $0.25\leq J_{2}\leq 0.5$, separately.

IV.2 Region $0\leq J_{2}\leq 0.25$

First, we investigate the $0\leq J_{2}\leq 0.25$ region where the standard AF order is stable in the regular model of $\Delta=0$. In Fig. 8, we show the squared sublattice magnetization $m_{AF}^{2}$ plotted versus $1/\sqrt{N}$ for the cases of (a) $J_{2}=0$ and (b) $J_{2}=0.15$, for various values of the randomness $\Delta$ spanning from the regular case of $\Delta=0$ to the maximal randomness case of $\Delta=1$.

For $J_{2}=0$, i.e., for the nearest-neighbor model without frustration, Fig. 8(a) shows that $m_{\infty}^{2}$ is extrapolated to positive values for any $\Delta$, indicating the AF LRO stabilized up to the maximal randomness. This behavior of the random honeycomb-lattice nearest-neighbor model differs from that of the random triangular-lattice nearest-neighbor model where the AF order becomes unstable at a finite amount of randomness of $\Delta_{c}\simeq 0.5$, giving way to the random-singlet state for stronger randomness $\Delta>\Delta_{c}$.Watanabe

By contrast, as can be seen from Fig. 8(b), for $J_{2}=0.15$, there exists a finite critical randomness of $\Delta_{c}\simeq 0.6$ beyond which the AF LRO vanishes. This observation, together with the previous data on the random triangular model, Shimokawa demonstrates that a certain amount of frustration is necessary to destabilize the AF LRO by introducing the randomness. We have made a similar analysis for other values of $J_{2}\leq 0.25$, and draw the AF-random-singlet phase boundary (marked by red points) in the phase diagram of Fig. 7.

In order to investigate the possible appearance of other types of magnetic order, we compute the spin freezing parameter $\bar{q}$ defined by

[TABLE]

This quantity can detect any type of static spin order, even including the random one such as the spin-glass order. The computed $1/\sqrt{N}$-dependence of $\bar{q}$ is shown in Fig. 9 for the cases of (a) $J_{2}=0$ and (b) $J_{2}=0.15$. The inset exhibits a magnification of the larger $N$ region. Note that, when the system possesses the AF LRO, $\bar{q}$ also becomes nonzero. The interest here is whether $\bar{q}$ could be nonzero in the parameter region without the AF LRO. As can be seen from these figures, whether the extrapolated $\bar{q}$ is positive or negative (zero) well correlates with the behavior of $m_{AF}^{2}$ shown in Fig. 8, indicating that no magnetically ordered state other than the standard AF order appears in the parameter range studied. In particular, no indication of the spin-glass order stabilized by the introduced randomness is obtained. This means that the state observed at $J_{2}=0.15$ for a stronger randomness of $\Delta\gtrsim 0.6$ is a nonmagnetic state without any static spin order.

We also compute the mean spin-gap energy $\Delta E=[\Delta E]_{J}$ of the random model (data not shown here), to find that $\Delta E$ is always extrapolated to vanishing values within the error bar, indicating that the nonmagnetic state found in the $\Delta\gtrsim 0.6$ region for $J_{2}=0.15$ is gapless, distinct from the gapped I and II phases.

We also compute the static and dynamical spin structure factors. In Fig. 10, we show the static spin structure factors for the maximally random case of $\Delta=1$ for various $J_{2}$-values. For the case of $J_{2}=0.15$, as can be seen from Fig. 10(a), the peaks appear at the same $q$-points as the corresponding regular one, i.e., at the $\Gamma$ point, although they are much broadened as compared with the regular case, reflecting the SRO nature of the associated AF order.

The dynamical spin structure factor $S_{\bm{q}}(\omega)$ is defined by

[TABLE]

where $E_{0}$ is the ground-state energy, and $\eta$ is a phenomenological damping factor taking a sufficiently small positive value. We employ the continued fraction method to compute $S_{\bm{q}}(\omega)$, Gagliano putting $\eta=0.02$.

The $\omega$-dependence of $S_{\bm{q}}(\omega)$ computed at the $\Gamma$ point and at the maximal randomness of $\Delta=1$ is shown in Fig. 11 for the cases of (a) $J_{2}=0$ and (b) $J_{2}=0.15$. In both Figs. 11(a) and 11(b), a sharp peak is observed in the small-$\omega$ region. In the regular case of $\Delta=0$, this peak is a finite-size counterpart of the delta-function peak expected at $\omega=0$ in the thermodynamic limit. Close inspection reveals that the small-$\omega$ peak in the maximally random case of $\Delta=1$ exhibits mutually different behaviors between (a) $J_{2}=0$ and (b) $J_{2}=0.15$. In Fig. 11(a), the peak tends to grow in height, while its position tends to move toward $\omega=0$ as the size $N$ is increased. By contrast, in Fig. 11(b), the peak height tends to saturate even when the system size $N$ is increased. This difference corresponds to the observation that in the maximally random case, the AF LRO is present for $J_{2}=0$, but is absent for $J_{2}=0.15$. Another interesting observation is that $S_{\bm{q}}(\omega)$ at $J_{2}=0.15$ and $\Delta=1$, corresponding to the nonmagnetic state, exhibits a very broad background component and a tail extending to larger values of $\omega$, in addition to a relatively sharp peak at a smaller $\omega$. Such a broad feature is a characteristic of the random-singlet state as previously studied in detail for the cases of the bond-random triangular and kagome models. Shimokawa In fact, our data of $S_{\bm{q}}(\omega)$ resembles the corresponding one of the random triangular model. This is the second indication that the nonmagnetic state observed for the $\Delta\gtrsim 0.6$ region at $J_{2}=0.15$ is indeed a random-singlet state.

IV.3 Region $0.25\leq J_{2}\leq 0.5$

Next, we move to the larger-$J_{2}$ region of $0.25\leq J_{2}\leq 0.5$, which corresponds to the gapped I and II phases of the regular model of $\Delta=0$. In order to obtain information of the possible magnetic LRO, we compute the size dependence of the freezing parameter $\bar{q}$, and the result is shown in Fig. 14 for (a) $J_{2}=0.3$ and (b) $J_{2}=0.5$, each corresponding to the gapped I and gapped II phases, respectively. As can be seen from the figures, $\bar{q}$ is extrapolated to zero within the error bar for any value of $\Delta$, indicating that the ground state in this region is always nonmagnetic. (Precisely speaking, for the case of $\Delta=1$, the fit using all the data points yields a slightly positive $\bar{q}$, while the fit using only larger-$N$ data points of $N\geq 24$ yields a vanishing $\bar{q}$ within one $\sigma$.)

In Fig. 12, we show the size dependence of the spin-gap energy $\Delta E$ for (a) $J_{2}=0.3$ and (b) $J_{2}=0.5$. Interestingly, the extrapolated $\Delta E$ is zero, i.e., the system is gapless for larger $\Delta>\Delta_{c}$, while it becomes nonzero, i.e., gapped for smaller $\Delta<\Delta_{c}$. (In our data fit of Fig. 12(b), the data of $N=18$, 26, and 30 for $J_{2}=0.3$ are excluded since they largely deviate from other data.) The values of $\Delta_{c}$ are estimated to be $\Delta_{c}\simeq 0.2$ for $J_{2}=0.3$ and $\Delta_{c}\simeq 0.4$ for $J_{2}=0.5$. The gapped states for smaller $\Delta$ correspond to the gapped I and II states discussed in the previous section. The change observed between the gapless and gapped behaviors on increasing $\Delta$ suggests the occurrence of a randomness-induced phase transition. The gapless nonmagnetic phase stabilized at $\Delta>\Delta_{c}$ is likely to be the random-singlet state.

Further information can be obtained from the static and dynamical spin structure factors. Figure 10 has shown the static spin structure factor for the maximally random case of $\Delta=1$ for various $J_{2}$-values. For the case of $J_{2}=0.5$, as can be seen from Fig. 10(c), the peak appears at the same $q$-points as the corresponding regular case, i.e., at the M points. Meanwhile, for the case of $J_{2}=0.3$ shown in Fig. 10(b), the peak structure itself is hardly discernible.

In Fig. 13, we show the $\omega$-dependence of the corresponding dynamical spin structure factor $S_{\bm{q}}(\omega)$ for $\Delta=1$, i.e., at (a) $J_{2}=0.3$ computed at the $\Gamma$ point, and at (b) $J_{2}=0.5$ computed at the M points. As can be seen from the figures, the computed $S_{\bm{q}}(\omega)$ exhibits much less peaky behavior as compared with the regular case, possessing a very broad component with a long tail extending to a larger $\omega$. This feature is a characteristic of the random-singlet state of the random triangular and kagome models studied in Refs. Kawamura and Shimokawa . (The data of Fig. 13 appear to resemble the random kagome model more than the random triangular model.) Anyway, such resemblance also justifies our present identification of the randomness-induced gapless nonmagnetic state stabilized for larger $\Delta$ as the random-singlet state, as observed in the $0\leq\Delta\leq 0.25$ region.

Collecting all the results above, we construct the ground-state phase diagram in the $J_{2}$-$\Delta$ plane as shown in Fig. 7. The phase boundary between the AF state and the random-singlet state (red points in Fig. 7) is determined from the behavior of $m_{AF}^{2}$, while that between the gapped I and II phases and the random-singlet phase (blue points in Fig. 7) is determined from the spin gap $\Delta E$. The phase boundary between the gapped I and II phases (green points in Fig. 7) is determined from the peak location of $S_{\bm{q}}$.

V Random model: Finite-temperature properties

In this section, we investigate the finite-temperature properties of the random $J_{1}$-$J_{2}$ Heisenberg model on the honeycomb lattice, focusing on the specific heat and the uniform susceptibility. To compute the thermal average of physical quantities, we employ the Hams–de Raedt method. HamsRaedt Following the previous sections, we present the results for the regions of (i) $0\leq J_{2}\leq 0.25$ and (ii) $0.25\leq J_{2}\leq 0.5$, separately.

V.1 Region $0\leq J_{2}\leq 0.25$

In Fig. 15, we show the temperature dependence of the specific heat per spin for (a) $J_{2}=0$ and (b) $J_{2}=0.15$, and of the susceptibility per spin for (c) $J_{2}=0$ and (d) 0.15.

As can be seen from Fig. 15(b), in the region of the random-singlet state, e.g., at $J_{2}=0.15$ and $\Delta=1$, the computed low-temperature specific heat exhibits a $T$-linear behavior as generically expected for the random-singlet state,Watanabe ; Kawamura while it exhibits a stronger curvature in the AF state. Likewise, as can be seen from Fig. 15(d), the susceptibility tends to exhibit a gapless behavior with a Curie tail in the region of the random-singlet state, just as expected for that state.Watanabe ; Kawamura Hence, our finite-temperature data also supports the identification that the ground state at $J_{2}=0.15$ and $\Delta=1$ is the random-singlet state.

V.2 Region $0.25\leq J_{2}\leq 0.5$

Next, we turn to the region $0.25\leq J_{2}\leq 0.5$. In Fig. 16, we show the temperature dependence of the specific heat for (a) $J_{2}=0.3$ and (b) $J_{2}=0.5$, and of the susceptibility for (c) $J_{2}=0.3$ and (d) $J_{2}=0.5$.

As can be seen from Fig. 16(a) and (b), in the region of the random-singlet state, e.g., at $\Delta=0.6$ and 1 for both $J_{2}=0.3$ and 0.5, the computed low-$T$ specific heat tends to exhibit a $T$-linear behavior as expected for the random-singlet state.Watanabe ; Kawamura

As mentioned in sect. III, the gapped I and II phases stabilized in this regime for weaker randomness are likely to be the plaquette and the staggered-dimer states, each of which spontaneously breaks the $Z_{3}$ symmetry, the lattice translation in the former and the lattice rotational in the latter. In 2D, a spontaneously breaking $Z_{3}$ symmetry might accompany a finite-$T$ transition with a divergent specific heat, possibly lying in the universality class of the three-state Potts model, as was pointed out in Ref. Mulder for the latter case. As can be seen from Figs. 16(a) and 16(b), the computed specific heat exhibits a double-peak structure in this region, where the lower-temperature peak might correspond to such a finite-temperature $Z_{3}$-symmetry-breaking transition. In that case, the low-temperature peak would diverge in the thermodynamic limit. Unfortunately, our lattice size $N=24$ presently available is too small to prove or disprove such a theoretical expectation.

As can be seen from Figs. 16(c) and 16(d), the susceptibility in the region of the random-singlet state tends to exhibit a gapless behavior with a Curie tail. Watanabe ; Kawamura This also lends support to our identification of the random-singlet state in the phase diagram. In Fig. 16(d), the Curie tail for $\Delta=1$ is very weak and is hardly visible in the temperature range studied. However, a finite fraction of samples turns out to have triplet ground states, suggesting that a weak Curie tail would eventually show up at sufficiently low temperatures even in this case.

VI Summary and discussion

Both the ground state and the finite-temperature properties of the random-bond $s=1/2$ $J_{1}-J_{2}$ Heisenberg model on the honeycomb lattice are investigated by the ED and the Hams–de Raedt methods. The ground-state phase diagram is constructed in the randomness ($\Delta$) versus the frustration ($J_{1}/J_{2}$) plane in order to obtain insight into the role of randomness and frustration in stabilizing various phases. Without frustration, i.e., for $J_{2}=0$, the AF LRO is kept stable up to the maximal randomness of $\Delta=1$. Hence, frustration plays a role in destabilizing the AF order. In other words, just the randomness is insufficient to induce the random-singlet state. In the phase diagram, we found three types of nonmagnetic states stabilized. For the regular and the weakly random cases, we found, on increasing $J_{2}$, the gapped I and gapped II phases, each likely to be the plaquette and the staggered-dimer phases, respectively. For the case of stronger randomness, we generically found a randomness-induced gapless QSL-like state, a random-singlet state, essentially of the same type as previously identified for the random triangular and kagome models. The observed robustness of the random-singlet state suggests that the gapless QSL-like behaviors might be realized in a wide class of frustrated quantum magnets possessing a certain amount of randomness or inhomogeneity, without fine-tuning the interaction parameters.

Note that we discussed the random-singlet state, of the type where the high symmetry of the underlying lattice is observed at the macroscopic level, even though each $J_{ij}$ realization microscopically breaks the high symmetry of the lattice. In some other situations, the lattice symmetry might be lowered even at the macroscopic level via, e.g., the possible uniform JT distortion as in the case of orthorhombic samples of Ba3CuSb2O9. Such a symmetry lowering would often enhance the spatially aligned singlet formation such as the VBC, and tends to induce a finite spin gap, i.e., it tends to induce the gapped nonmagnetic state rather than the gapless nonmagnetic state.

On the basis of these findings, we now wish to discuss the possible experimental implications of our present results, particularly regarding the honeycomb-lattice-based magnets 6HB-Ba3NiSb2O9 and Ba3CuSb2O9. First, we wish to discuss the QSL-like behavior recently observed in the $s=1$ honeycomb AF 6HB-Ba3NiSb2O9. Cheng ; Mendels-Ni ; Darie This compound is found to exhibit gapless QSL-like behaviors, accompanied by the $T$-linear low-$T$ specific heat and the gapless susceptibility with a Curie tail. As mentioned in Ref. Mendels-Ni , the exchange interaction in this material might be of the $J_{1}$-$J_{2}$ honeycomb-type, where a considerable amount of structural disorder exists. The upper limit of the exchange disorder was estimated to be $\Delta\lesssim 0.25$. Mendels-Ni Although the ratio $J_{2}/J_{1}$ is not known precisely, the structure of this compound suggests a moderately large $J_{2}$. Of course, the $s=1$ nature of this compound might somewhat modify the phase diagram of the $s=1/2$ model obtained here. For example, in the case of the regular model, the $J_{2}$-value at the AF-plaquette phase boundary was estimated to be $J_{2c}\simeq 0.3$, Gong-s=1 ; Li-s=1 which should be compared with the corresponding estimate $J_{2c}\simeq 0.25$ for the $s=1/2$ model. Even with such an uncertainty, the existence of a significant amount of exchange randomness of $\Delta\simeq 0.25$ is likely to locate 6HB-Ba3NiSb2O9 lying in the random-singlet state or close to its phase boundary: see Fig. 7. Hence, a good possibility seems to exist that the experimentally observed gapless QSL-like behavior of 6HB-Ba3NiSb2O9 might indeed be that of the random-singlet state. In fact, in Ref. Mendels-Ni , within the modeling of 6HB-Ba3NiSb2O9 by the random triangular model, the possibility of the random-singlet state as proposed in Ref. Watanabe for the random triangular model was examined, but led to a negative result, arguing that the extent of the randomness deduced for 6HB-Ba3NiSb2O9 is smaller than the critical randomness of the random triangular model, $\Delta_{c}\simeq 0.5$. Yet, as mentioned above, this compound is likely to be better modeled as the random $J_{1}-J_{2}$ honeycomb model with moderately large $J_{2}$, to which our present study suggests $\Delta_{c}\simeq 0.25$ much smaller than the corresponding value of the random triangular model. Hence, the modeling of 6HB-Ba3NiSb2O9 as the $J_{1}$-$J_{2}$ honeycomb model leaves a good possibility of the random-singlet state. Of course, care has to be taken in a truly quantitative comparison, because the type of randomness assumed here is simplified in that it obeys a simple uniform distribution and is taken to be completely bond-independent.

For the hexagonal sample of Ba3CuSb2O9, in Ref. Nakatsuji , two scenarios were proposed, i.e., a random static JT-distortion-driven ‘random-singlet state’ versus a dynamical JT-distortion-driven ‘spin-orbital liquid state’. The former would essentially be of the same character as the random-singlet state discussed in the present paper. Meanwhile, many recent experimental studies point to the second scenario, although the NMR study Mendels-Cu and the density-functional calculation Shanavas suggested the first scenario. If the second scenario applies, to properly understand the QSL-like behaviors observed in the hexagonal sample of Ba3CuSb2O9, one would need to consider the effect of fluctuating orbital degrees of freedom. Ishiguro ; Hagiwara ; Katayama ; Nasu ; Mila ; Nasu2

For the orthorhombic sample, on the other hand, a symmetry-lowering static JT distortion occurring at a higher temperature might somewhat modify the nature of the low-temperature spin state (also called ‘random-singlet state’ in Ref. Do ) from the one discussed in the present paper, with an enhanced character of the VBC and a finite spin gap. Quantitative details of the gapless/gapful issue would depend on the competition between the extent of the randomness and that of the uniform distortion in $J_{ij}$, and needs further clarification.

In any case, the randomness-induced gapless QSL-like state, the random-singlet state, prevails in quantum magnets on a variety of frustrated lattices, including not only the $J_{1}-J_{2}$ honeycomb lattice as studied here, but also the triangular and the kagome lattices, if a certain amount of randomness or inhomogeneity is introduced in some way or another. This would further extend our concept of the QSL state as a novel state of matter.

Acknowledgements

The authors wish to thank T. Shimokawa and T. Okubo for valuable discussion. This study was supported by JSPS KAKENHI Grant Number JP25247064. Our code was based on TITPACK Ver.2 coded by H. Nishimori. We are thankful to ISSP, the University of Tokyo, and to YITP, Kyoto University, for providing us with CPU time.