Quadrupolar quantum criticality on a fractal
Jonathan D'Emidio Simon Lovell Ribhu K. Kaul




Department of Physics & Astronomy, University of Kentucky, Lexington, KY-40506-0055

Department of Physics & Astronomy, University of
  Kentucky, Lexington, KY-40506-0055


Department of Physics & Astronomy, University of Kentucky, Lexington, KY-40506-0055

We study the ground state  ordering of quadrupolar ordered S=1
magnets as a function of spin dilution probability p on the
triangular lattice. In sharp contrast to the
ordering of S=1/2 dipolar Néel
magnets on percolating clusters, we find that the quadrupolar magnets are
quantum disordered at the percolation threshold, p=p^*. Further we
find that long-range quadrupolar order is present for all p<p^* and
vanishes first exactly at p^*. Strong evidence
for scaling behavior close to p^* points to an
unusual quantum criticality without fine tuning that arises from an interplay of quantum fluctuations and randomness.

Quadrupolar quantum criticality on a fractal
    Ribhu K. Kaul
    December 30, 2023
============================================







§ INTRODUCTION


Quantum spin models with random vacancies provide a rich
playground for the study of the combined effects of strong interactions and
quenched disorder that are relevant to experimental measurements on a
number of doped magnetic
alloys <cit.>. Magnetic order in higher dimensions is generally stable to a small concentration of vacancies. On the other
 hand, if its moments are diluted beyond the percolation
 threshold, the system breaks up into zero-dimensional clusters and
hence magnetic order must be lost. What is the nature of the
 magnetic quantum phase transition as the dilution is varied? 

The answer to this question is most thoroughly understood for
the randomly diluted transverse field Ising
model. The expected ground state phase diagram <cit.>
as a function of the transverse field g and
site dilution probability p, has three kinds of
quantum phase transitions, Fig. <ref>. For g>g_ C the
transition takes place at a value of p smaller than the percolation threshold p_c, the critical
point is then described by the random transverse field Ising model
which possesses an infinite randomness fixed point <cit.>; the physics of percolation
however plays no role (“R” in Fig. <ref>). For
g<g_ C on the other hand percolating clusters are magnetically ordered and the
singularities at the transition are determined by those of classical
percolation (“P” in
Fig. <ref>) <cit.>. 
Finally exactly at g=g_ C a critical point is obtained: Crucial to our discussion, quantum criticality
at the percolation threshold  (“C” in Fig. <ref>) requires fine tuning of
the quantum fluctuations and hence “C” is not the generic transition for
Ising magnets on dilution.

Another family of quantum spin models where this question has been
addressed in detail are bipartite Néel ordered Heisenberg models. For random
depletion of the square lattice it was found that percolating clusters
have long range order <cit.>, the destruction is thus
akin to Fig. <ref> “P”. The surprising stability of
the Néel order has been traced back to the effect of uncompensated Berry phases
which result in “orphan” moments <cit.>.  In the bilayer geometry,
with bond dilution that circumvents random Berry phases, a
phase diagram similar to Fig. <ref> was found <cit.>, with quantum
phase transitions corresponding to “P” <cit.> in some regions and
corresponding to “R” in others, in contrast to the Ising case the critical phenomena here
is controlled by a finite disorder fixed point with conventional scaling <cit.>. The critical point “C” which
requires fine tuning has also been studied <cit.>. A
variety of phase diagrams and critical phenomena can be accessed by
dilution of Néel magnets in geometries different from the single and
bilayer <cit.>.





The two examples of dipolar ordered magnets (Ising and Heisenberg)
make it clear that the details such as symmetry of the order parameter
and Berry phases play a crucial role in
determining how quantum magnetism is destroyed by random dilution. In this
work we address the dilution transition in a different kind of
system: quadrupolar
ordered (also referred to as spin nematic) magnets in their most common S=1 realization <cit.>. The study
of the quadrupolar phase of S=1 magnets has become increasingly
popular motivated by their possible sighting in some triangular lattice
Ni based magnets (see e.g. <cit.>). Experimental studies of Zn replacement of Ni
in NiGa_2S_4 provide a direct experimental motivation for the site
diluted S=1 magnets we study here <cit.>. As we
shall describe below, in contrast to what has been
observed for Néel order in S=1/2 magnets and transverse field Ising models, we
find here that quadrupolar order vanishes exactly at the percolation threshold without any
fine tuning, resulting in new quantum critical behavior.





§ MODEL


Just as the bilinear Heisenberg model is the archetype for
realizing Néel order, the biquadratic Heisenberg interaction is the
archetypical model system that realizes quadrupolar order in S=1
magnets. We shall consider the following Hamiltonian,

    H = - ∑_⟨ ij⟩ J_ij (S⃗_i ·S⃗_j )^2,

 where ⟨ ij ⟩ denote nearest neighbors in the triangular
 lattice, and S⃗_i are the S=1 Pauli matrices on the i^
   th site. The model possesses an explicit physical SO(3) internal spin
 symmetry.  We use standard uncorrelated random site and bond
 dilution: For site dilution J_ij=|J| r_i r_j where r_i=0 with
 probability p and r_i=1 with probability 1-p. For bond dilution
 J_ij=|J| χ_ij where χ_ij=0 with probability p and
 χ_ij=1 with probability 1-p. While site dilution is closer
 to experimental realizations such as Zn doped NiGa_2S_4,
 bond dilution provides an alternate access to the percolation
 threshold. 

Our numerical results were obtained
 using stochastic series expansion <cit.> with an efficient sign-problem free algorithm
 that has been described and tested
 previously <cit.>. We calculate averages on finite size
 clusters at finite-temperature, T. 
To study the quadrupolar order in the magnet we define the quadrupolar
order parameter through the equal time two point
correlation function of the quadrupolar order parameter.


    O_Q^2= ⟨1/N_c^2∑_a,i,j⟨Q̂^aa_i Q̂^aa_j⟩⟩_p,

where the quadrupolar order parameter is defined in terms of the local S=1
Pauli matrices as
Q̂^ab= Ŝ^a Ŝ^b+ Ŝ^bŜ^a/2-
2/3δ^ab, N_c is the size of the cluster, ⟨…⟩ is the quantum
mechanical average over the thermal density matrix and ⟨…⟩_p is a disorder average over realization of the ensemble of clusters.





§  NUMERICAL RESULTS


It is now well established that the model
Eq. (<ref>) without any depletion (p=0) has quadrupolar
order with ferromagnetically aligned directors <cit.> (i.e. O_Q^2 in Eq. <ref>
is finite in the ground state in the thermodynamic limit). We begin our numerical study by asking the
following question: As p is increased does the quadrupolar order parameter
eventually vanish like “R”, “C” or “P” (referring to Fig. <ref>).  As we have discussed in
the introduction,
generic expectations would be
either “R” or “P”.


 


In Fig. <ref> we show extrapolations of the quadrupolar
order parameter to
the thermodynamic limit at the percolation threshold. We define two different ensembles for
finite-size scaling which allow us to approach the thermodynamic limit in two
different ways, with fixed cluster size N_c (on an infinite
underlying lattice) and the largest cluster on a finite L× L
lattices, following <cit.>. We have also studied
both bond and site dilution at the percolation threshold. Consistently across
all finite-size scaling schemes we find that the order parameter vanishes
at the percolation threshold. For a similar scaling analysis for the
Néel order, see Fig. 10 in Ref. <cit.>
which clearly shows the ordering of S=1/2 Heisenberg model on the
percolating cluster in two dimensions. In contrast here we find the
quadrupolar order vanishes in the thermodynamic limit. This data clearly eliminates the
“P” possibility since this requires the percolating clusters to be
magnetically ordered.   We have taken great care to make sure our
data is equilibrated and in the limit of T=0. The ground state limit
requires extraordinarily low temperatures becasue of the weak links
that connect percolating clusters. Details are provided
in the Appendix. 





Now that we have shown that the percolating cluster are magnetically
disordered, the generic expectation is that the quadrupolar order
vanishes before the percolation threshold is reached at some p<p^*
as illustrated for “R”  in Fig. <ref>. 
To address this question quantitively we study a quadrupolar
“Binder” ratio R_Q≡⟨ O_Q^4
  ⟩/⟨ O_Q^2 ⟩^2. R_Q is expected
to be monotonically decreasing
with L in a quadrupolar ordered phase. When graphed as a function of the tuning
parameter p, R_Q data for different L are expected to cross at the phase
transition to a non-magnetic phase. Our zero temperature data for the
case of site dilution in Fig. <ref> clearly shows that the R_Q crosses and as the thermodynamic limit is reached  the crossing point approaches the percolation
threshold p^*=0.5 with high accuracy. This establishes that the
quadrupolar order vanishes continuously as the percolation threshold
is approached. This allows us to eliminate the
possibility “R”, in which the order parameter vanishes before the
percolation threshold is reached. We hence conclude that the order
parameter vanishes first precisely at the percolation threshold p^*, as
illustrated in the cartoon “C” shown in Fig. <ref>.





The vanishing of quadrupolar order right at the
percolation threshold raises the interesting possibility that at p=p^* the
system is quantum critical and the quadrupolar correlation possess
scale invariance. 


Random systems are well known to display a range of
novel scaling behavior. We present evidence however that 
our model possesses conventional power law scaling. In order to test the scaling hypothesis, we first study the dynamic scaling in
imaginary time of the order parameter at p=p^*. In
Fig. <ref> we show scaling collapse of the order
parameter data at the percolation threshold p=p^* for the maximum sized cluster
on an L× L lattice as a function of
the putative scaling variable L^zT. For a quantum critical system,
we expect the scaling form O_Q^2 = 1/L^1+η F ( L ^z
  T). Note the absence of a tuning parameter (we only require the
  system to be at the percolation threshold) that would normally be
  expected for instance
  from previously studied phase diagrams of quantum rotor models in
  Fig. <ref> (the tuning parameter there is the choice of
  g=g_ C). We find an excellent collapse with small finite size
  corrections to the exponents.


An alternate test of scaling can be made by varying p away from the
percolation threshold. In Fig. <ref> we study the scaling of
the order parameter with the deviation from the percolation theshold
p-p^*. We find excellent scaling behavior assuming the simple
scaling form O_Q^2 = 1/L^1+η G (
(p-p^*)L^1/ν) . The value of η so obtained is in excellent
agreement with the estimate obtained from the previous scaling analysis 
at the percolation threshold Fig. <ref>.










§ DISCUSSION


We have presented extensive evidence that in two dimensions
quadrupolar order vanishes continuously as a function of dilution
right at the percolation threshold. At the percolation threshold the
systems shows conventional scaling behavior symptomatic of a finite
disorder fixed point. Our study is an unusual example of quantum
criticality at the percolation threshold without any fine tuning. We
contrast this with the generic phase diagram Fig. <ref>
where at the percolation threshold, tuning of quantum fluctuations to
a special value if required to observe quantum criticality.

An interesting open question is whether the critical exponents we have
found are universal or can vary continuously as non-universal
properties of the magnet are varied at the percolation threshold. Such
varying exponents were reported in a study of specific diluted
dimerized S=1/2 magnet <cit.>. Interestingly, the dynamic critical exponent z we find from our scaling analysis is
numerically somewhat close to
d_f=91/48≈ 1.8958… (the fractal dimension of a percolating
cluster). A number of previous works have found or predicted such scaling at the
percolation threshold and our finding could be consistent with such
behavior <cit.>,
though numerically z=2 would also be consistent with our value
within errors. 
A complete theory of the unusual behavior and scaling we have found here is
an interesting direction for future work.

We acknowledge financial support fron NSF DMR-1611161. The numerical
simulations reported here were carried out on the DLX cluster at
University of Kentucky and by resources allocated by XSEDE.











§ EQUILIBRATION OF DISORDERED CLUSTERS



In order to ensure the proper equilibration of our disordered clusters, we have adopted the same equilibration and measurement protocol as in Ref. sandvik2002:s12perc.  For each disorder realization we begin by performing N_e equilibration sweeps at some large initial temperature (β≈ 1, setting |J|=1), followed by N_m more measurements sweeps, then again N_e equilibration sweeps, and finally N_m measurement sweeps.  Once this process is complete we perform β doubling on our configurations (see appendix <ref>) and start the process again.  Thus for each disorder realization we have two separate measurement cycles at many different values of β.

Separating out two distinct measurement segments allows us to check the equilibration of our disordered clusters.  In order to achieve the best QMC averages with minimal computational time, we have set N_m=2 N_e.  We can then study the percent difference of our two measurement segments (each averaged over disorder realizations) as a function of N_m to determine its optimal value, which is given in Fig. <ref>.  We observe that at the percolation threshold (where the clusters are most fragmented), the difference between our first and second measurement cycles becomes statistically insignificant when N_m ≈ 64.  We therefore cautiously set N_m=200 throughout the course of our numerical studies.


 





§ ZERO TEMPERATURE CONVERGENCE



The most computationally expensive component of our numerical studies
is to converge our configurations to the ground state on disordered
clusters.  In order to achieve this, extremely large values of the
inverse temperature are needed (relative to that of a clean lattice).
Furthermore, the free energy landscape of configurations is rugged,
and thus quenching a randomly initialized starting configuration to
low temperature abruptly may leave it stuck in a local minimum.  In
order to circumvent this issue, and to efficiently reach the low
temperatures required with minimal equilibration, we implement the
β doubling procedure <cit.>.  The procedure works as follows: given a disorder realization, we begin the equilibration and measurement cycles at some initial high temperature (β_0 on the order of unity).  Configurations at high temperatures are relatively easy to equilibrate, since there are few operators acting in the operator string.  After the equilibration and measurement cycles, the value of β and the QMC configuration are both doubled.  In terms of the configuration, this corresponds to repeating the operator sequence twice.  This procedure gives another valid partition function configuration that is close to being equilibrated at an inverse temperature 2β_0.  From here the equilibration and measurement sequence is again carried out.  The process is continued until the final target value β_max is reached.  

In order to illustrate the convergence of our zero temperature data, in Fig. <ref> we show the binder ratio for max clusters on an L=32 lattice near the site percolation threshold for all of our β doubled values.  We see that very large values of β are required to converge to the ground state.  This is a reflection of the fact that z is close to 2 at this disordered quantum critical point, meaning that doubling the size of the lattice would require quadrupling β to remain near the ground state.

 

