Site dilution of quantum spins in the honeycomb lattice

We consider the effect of random site dilution on a honeycomb lattice of quantum spins described by the antiferromagnetic Heisenberg spin-$S$ model. Using linear spin-wave theory, we compute the zero-temperature magnetization and density of states as a function of dilution up to the classical percolation threshold. Several subtle issues regarding the treatment of quasidivergent zero-energy modes, which appear in the real-space formulation of the spin-wave problem, are clarified. For $S>1∕2$, the spin-wave theory is well defined in the sense that results at all dilution concentrations are consistent with the underlying assumptions of the theory. For $S=1∕2$, however, the approximation breaks down. In this case, we have studied the effect of dilution on the staggered magnetization using the stochastic series expansion Monte Carlo method. Two main results are to be stressed from the Monte Carlo calculation: (i) an improved estimate for the staggered magnetization of the undiluted system $m_{\mathrm{av}}(L\to \infty )=0.2677\left(6\right)$ and (ii) a finite value of the staggered magnetization of the percolating cluster at the classical percolation threshold, showing that there is no quantum critical transition driven by dilution in the Heisenberg model. We have used the computed staggered magnetization and density of states to analyze neutron scattering experiments and Néel temperature measurements of two quasi-two-dimensional diluted honeycomb systems: (i) ${\mathrm{Mn}}_p{\mathrm{Zn}}_{1-p}\mathrm{P}{\mathrm{S}}_3$ (a diluted $S=5∕2$ system) and (ii) $\mathrm{Ba}{\left({\mathrm{Ni}}_p{\mathrm{Mg}}_{1-p}\right)}_2{\mathrm{V}}_2{\mathrm{O}}_8$ (a diluted $S=1$ system). We have found that our calculations are in good agreement with the experimental data.

Figure 1

(Color online) A finite size honeycomb lattice showing the periodic boundary conditions used in the numerical calculations.

Figure 2

(Color online) An example of the cluster formation for a particular disorder realization in a lattice with $L=10$. The original lattice is shown in (a); in (b) each site is chosen to be diluted with probability $1-p_c$, with $p_c=0.697043$; in (c) the existent clusters (12 in this case) have been determined with rigid boundary conditions; in (d) the largest cluster found in (c) (blue) is augmented by the periodic boundary conditions.

Figure 3

(Color online) The solid (blue) lines show the distribution of the number of sites in the largest cluster of a randomly site-diluted honeycomb lattice. From top to bottom, the three panels correspond to dilution levels $x=0.1,0.5,0.8$. From left to right, the peaks correspond to linear system sizes $L=4,5,\dots ,18$. The vertical dotted (red) lines indicate the average cluster sizes ${\overline{N}}_c$, computed as per Eq. (44).

Figure 4

(Color online) The effective scaling dimension of the largest cluster takes one of two values in the $L\to \infty$ limit: $D_{\mathrm{eff}}=2$ $(0⩽x<1)$ or $D_{\mathrm{eff}}=91∕48$ $(x=1)$. For $x\lesssim 0.5$, $D_{\mathrm{eff}}$ is close to its asymptotic value at all system sizes. When $x$ is close to 1, very large system sizes are necessary to reach the asymptotic regime. The figure inset shows the largest-cluster size distribution at percolation plotted in reduced coordinates. Each curve is computed as a histogram over ${10}^5$ disorder realizations for system sizes $L=5,6,\dots ,48$. As $L\to \infty$, the finite-size results converge to a smooth scaling function (one not dissimilar from that of the square-lattice case; see Fig. 2 of Ref. 5).

Figure 5

(Color online) Finite-size scaling of the average quantum correction to the staggered magnetization $⟨\delta m_z(p,L)⟩$ for different values of dilution $x=(1-p)∕(1-p_c)$. Each point was obtained after ${10}^5$ disorder realizations of lattices with equal number of sites in each sublattice. Also shown for $x=1$ is the result obtained when the realized lattices are not constrained to have $N_a=N_b$: $*$ zero modes were subtracted and the highest amplitude (nonzero) mode (see text) was subtracted if $N_a\ne N_b$; $**$ only zero modes were subtracted. For $x=0$ the RSS result was obtained by a reciprocal space sum using the analytical result (Ref. 24).

Figure 6

Contributions to $⟨\delta m_z(p_c,L)⟩$ from the lower nonzero energy mode (upper panel) and the lower energy mode higher than the lower nonzero energy mode (lower panel). The average was taken over ${10}^4$ disordered honeycomb lattices, with $N_a-N_b=\pm 1$, at $p=p_c$.

Figure 7

(Color online) Average quantum-mechanical factor $m_{\mathrm{QM}}\left(x\right)$ vs dilution $x=(1-p)∕(1-p_c)$ for different values of the spin $S$. Panel (a) shows the results for the honeycomb lattice and panel (b) for the square lattice.

Figure 8

(Color online) Average magnetic moment per magnetic site as a function of dilution $p$. The linear spin-wave result for the $S=5∕2$ Heisenberg antiferromagnet in the honeycomb lattice (red squares) is compared with neutron scattering data on ${\mathrm{Mn}}_p{\mathrm{Zn}}_{1-p}\mathrm{P}{\mathrm{S}}_3$ from Ref. 2 (gray circles).

Figure 9

(Color online) $T_N\left(p\right)∕T_N\left(0\right)$ vs $p$ for $S=5∕2$. The results obtained by numerically solving Eq. (52) (squares) and the mean-field result of Eq. (56) (diamonds) are compared with experimental results on ${\mathrm{Mn}}_p{\mathrm{Zn}}_{1-p}\mathrm{P}{\mathrm{S}}_3$ from Ref. 2 (circles).

Figure 10

(Color online) $T_N\left(p\right)∕T_N\left(0\right)$ vs $p$ for $S=1$. The results obtained by numerically solving Eq. (52) (squares) and the mean-field result of Eq. (56) (diamonds) are compared with experimental results on $\mathrm{Ba}{\left({\mathrm{Ni}}_p{\mathrm{Mg}}_{1-p}\right)}_2{\mathrm{V}}_2{\mathrm{O}}_8$ from Ref. 19 (circles).

Figure 11

(Color online) DOS of the Heisenberg antiferromagnetic model in the linear spin-wave approximation for the square lattice. The four panels correspond to different dilution levels $x$. An energy mesh with spacing 0.01 in units of $JS$ was used. These results were obtained applying the recursion method to systems with $L=128$, and averaging over 200 to 400 disorder realizations. The dotted line is the clean limit DOS.

Figure 12

(Color online) DOS of the Heisenberg antiferromagnetic model in the linear spin-wave approximation for the honeycomb lattice. The four panels correspond to different dilution levels $x$. An energy mesh with spacing 0.01 in units of $JS$ was used. These results were obtained applying the recursion method to systems with $L=128$, and averaging over 200 to 400 disorder realizations. The dotted line is the clean limit DOS.

Figure 13

(Color online) Low-energy behavior of the DOS of the Heisenberg antiferromagnetic model in the linear spin-wave approximation for the honeycomb (left) and square (right) lattices at percolation. An energy mesh with spacing $5\times {10}^{-4}$ in units of $JS$ was used. These results were obtained applying the recursion method to systems with $L=128$, and averaging over 800 disorder realizations.

Figure 14

(Color online) The staggered magnetization of the undiluted honeycomb lattice ($x=0$, $N_c=N=2L^2$) is extrapolated to the thermodynamic limit following Eqs. (58a, 58b). A simultaneous fit of the two data sets yields the value $m_{\mathrm{av}}(L\to \infty )=0.2677\left(6\right)$.

Figure 15

(Color online) The main plot shows an extrapolation to the thermodynamic limit of twice the $z$-projected staggered magnetization for various dilution levels $x$ (as indicated by the symbols in the upper-left legend). The lines drawn through the data points represent a global fit to Eqs. (58) in which $m_{\mathrm{QM}}\left(x\right),a_1\left(x\right),b_1\left(x\right),\dots$, are treated as powerseries in $x$ and varied. The resulting function $m_{\mathrm{QM}}\left(x\right)$ appears as the solid (pink) line in the figure inset alongside Monte Carlo results for the square lattice (from Ref. 5) and spin-wave results for the honeycomb lattice. The (red) errorbars indicate the values of $m_{\mathrm{QM}}$ extrapolated from each fixed-$x$ dataset taken individually.

Figure 16

(Color online) The disorder-averaged coordination number $\overline{z}\left(x\right)$ is plotted as a function of dilution level for infinite square and honeycomb lattices. The difference between the two lattice types narrows as $x\to 1$. The undiluted square lattice is 33% more coordinated than the honeycomb lattice. At percolation, it is only 12% more so.