Superfluid-Mott glass quantum multicritical point on a percolating lattice

We employ large-scale Monte Carlo simulations to study a particle-hole symmetric site-diluted quantum rotor model in two dimensions. The ground state phase diagram of this system features two distinct quantum phase transitions between the superfluid and the insulating (Mott glass) phases. They are separated by a multicritical point. The generic transition for dilutions below the lattice percolation threshold is driven by quantum fluctuations while the transition across the percolation threshold is due to the geometric fluctuations of the lattice. We determine the location of the multicritical point between these two transitions and find its critical behavior. The multicritical exponents read $z=1.72(2)$, $\beta/\nu=0.41(2)$, and $\gamma/\nu=2.90(5)$ . We compare our results to other quantum phase transitions in disordered systems, and we discuss experiments.


Introduction
Randomly diluted quantum many-particle systems feature two kinds of fluctuations at zero temperature.First, geometric fluctuations due to bond and/or site dilution of the lattice lead to a percolation problem with a geometric phase transition at the corresponding percolation threshold [1].Second, quantum fluctuations are caused by non-commuting operators in the many-particle Hamiltonian.Combining both often leads to rich ground-state phase diagrams with several different quantum phase transitions.
Prototypical examples of such behavior can be found in magnetic systems (for a review see, e.g., Ref. [2]).The diluted transverse-field Ising model displays long-range order for dilutions p up to the percolation threshold p c provided that the transverse-field h x is below the critical field h c (p).For dilutions above the percolation threshold, long-range order is not possible since the system consists of decoupled finite-size clusters.As a result, the p-h x ground state phase diagram displays two quantum phase transitions, separated by a multicritical point at (p c , h * ) [3,4,5,6,7].The generic transition across h c (p) for 0 < p < p c is driven by quantum fluctuations while the percolation transition across p c for h x < h * is driven by the geometric criticality of the lattice.Other quantum magnets behave similarly.The diluted square lattice Heisenberg quantum antiferromagnet is long-range ordered at the percolation threshold [8,9].If the quantum fluctuations are strengthened, for example via an interlayer coupling in a bilayer quantum antiferromagnet, a phase diagram with a multicritical point appears, as in the transverse-field Ising case [10,11,12,13].
In the present paper, we investigate the interplay of geometric and quantum fluctuations in a system of interacting bosons on a diluted lattice that undergoes a quantum phase transition arXiv:1608.06862v1[cond-mat.quant-gas]24 Aug 2016 between a superfluid and an insulating ground state.Experimental applications include Josephson junction arrays [14,15], granular superconducting films [16,17], helium absorbed in vycor [18,19], as well as doped quantum magnets in high fields [20,21,22].In the presence of disorder, the superfluid phase and the Mott insulator are separated by a "glassy" quantum Griffiths phase [23,24,25] in which rare regions of superfluid order exist in an insulating bulk system.For generic disorder (without special symmetries), this intermediate phase is the Bose glass, a compressible gapless insulator [26,27,28].If the Hamiltonian is particle-hole symmetric even in the presence of disorder, the intermediate phase between the superfluid and the Mott insulator is instead the so-called Mott glass (or random-rod glass) [29,30] which is an incompressible gapless insulator.
While the quantum phase transition between superfluid and Bose glass has been studied in great detail, the quantum phase transition between superfluid and Mott glass has attracted comparatively less attention.Early numerical work in two dimensions using quantum Monte Carlo simulations [31,32] and a strong-disorder renormalization group [33] produced inconclusive and partially contradicting results for the critical behavior of this transition.To resolve this issue, we recently performed large-scale Monte Carlo simulations of a site-diluted two-dimensional quantum rotor model with particle-hole symmetry [34].Similar to the quantum magnets discussed above, this system features two distinct quantum phase transitions.For dilutions p below the percolation threshold p c of the lattice, we find a superfluid-Mott glass transition characterized by a universal (dilution-independent) critical behavior.The transition across the lattice percolation threshold p c falls into a different universality class.
In the present paper, we briefly summarize these results and then focus on the multicritical point separating the generic quantum phase transition (p < p c ) from the percolation quantum phase transition across the lattice percolation threshold.The paper is organized as follows.We introduce in Sec. 2 the quantum rotor Hamiltonian and its mapping onto a classical XY model with columnar disorder.Section 3 is devoted to the Monte Carlo method and the data analysis techniques.Results for the generic transition, the percolation transition, and the multicritical point are discussed in Sec. 4. We conclude in Sec. 5.

Diluted quantum rotor model
The quantum rotor Hamiltonian is a prototypical model of interacting bosons.We consider a diluted square-lattice array of superfluid grains, coupled by Josephson junctions.The Hamiltonian reads where φi is the phase operator of site (grain) i, and ni is the number operator canonically conjugate to the phase.The sum over ij comprises pairs of nearest neighbors on the square lattice.The parameters U and J are the charging energy and the Josephson coupling, respectively.ni is the offset charge at site i.Site dilution is introduced via the quenched random variables i which are independent of each other.They can take the values 0 (vacancy) with probability p and 1 (occupied site) with probability 1 − p.
For generic values of the offset charges ni or generic non-integer filling, the Hamiltonian (1) is not particle-hole symmetric, and the intermediate phase between superfluid and Mott insulator is a Bose glass.To preserve the particle-hole symmetry, we set all offset charges ni = 0 and fix the particle number n at a large integer value.In this case, the disorder introduced by the site dilution is off-diagonal, and the intermediate phase is a Mott glass.
The qualitative behavior of the particle-hole symmetric site-diluted quantum rotor model can be easily understood [27,30,34].If the Josephson energy J dominates over the charging energy U and if the site dilution p is below the lattice percolation threshold p c , the ground state is a long-range ordered superfluid.Long-range order can be destroyed either by increasing the charging energy (which enhances the quantum fluctuations) or by raising the dilution beyond the percolation threshold (where the lattice decomposes into disconnected finite-size clusters).
Because we plan to perform Monte Carlo simulations using a highly efficient classical cluster algorithm, we now map the quantum rotor Hamiltonian (1) onto a classical model.In the presence of particle hole symmetry, the mapping leads to a classical XY model on a cubic lattice [35].The classical Hamiltonian reads Here i denotes a position in two-dimensional real space and t is the "imaginary time" coordinate.
The dynamical variable S i,t at each lattice site (i, t) is an O(2) unit vector.Note that the vacancy positions do not depend on the imaginary time coordinate t.The classical XY model ( 2) therefore features columnar defects, i.e., the disorder is perfectly correlated in the imaginary time direction.The values J s and J τ of the interaction energies as well as the "classical" temperature T are determined by the mapping from the original quantum rotor Hamiltonian (the classical temperature differs from the real physical temperature which is zero).As we are interested in the critical behavior which is expected to be universal, the precise values of J s and J τ are not important.Consequently, we fix J s = J τ = 1 and tune the classical XY model ( 2) through its phase transitions by varying either the classical temperature T or the dilution p.
In the absence of dilution, the Hamiltonian (2) represents the usual three-dimensional XY model.The clean two-dimensional superfluid-Mott insulator quantum critical point is therefore in the three-dimensional classical XY universality class which features a correlation length exponent ν ≈ 0.6717 [36].The stability of the clean critical behavior against disorder is governed by the Harris criterion [37] dν > 2. Here, d is the number of dimensions in which there is randomness; for columnar defects in a cubic lattice this implies d = 2.The Harris criterion is therefore violated, and the three-dimensional clean XY critical point is unstable against columnar defects.

Monte Carlo method
To determine the critical behavior of the classical XY Hamiltonian (2), we perform Monte Carlo simulations using the highly efficient Wolff cluster algorithm [38] which greatly reduces the critical slowing down close to the phase transition.In addition to the Wolff updates, we also perform conventional Metropolis updates [39] to equilibrate small disconnected clusters of lattice site that can occur at higher dilutions p and may be missed by the Wolff algorithm.
Our simulations cover several dilutions p = 0, 1/8, 1/5, 2/7, 1/3, 9/25 and the lattice percolation threshold p c = 0.407253.Linear sizes range from L = 10 to 150 in the two space directions and L τ = 6 to 1792 in the imaginary time direction.All observables are averaged over 10,000 to 50,000 disorder configurations, depending on system size.Each sample is equilibrated for 100 full Monte Carlo sweeps, followed by a measurement period of 500 sweeps with a measurement taken after every sweep.(A full Monte Carlo sweep consists of a Metropolis sweep and a Wolff sweep, the latter is defined as a number of cluster flips such that the total number of flipped spins equals the number of lattice sites.)The choice of rather short runs for a large number of disorder configurations helps reducing the overall statistical error of the data, as discussed in Refs.[13,40,41,42,43].To eliminate systematic biases in observables stemming from the short simulation runs, we use improved estimators (see, e.g., the appendix of Ref. [43]).
We analyze the Monte Carlo data by finite-size scaling.As the columnar disorder in the classical XY model breaks the symmetry between the space and imaginary time directions, observables generally do not scale in the same way with the system sizes L and L τ .We consider, for example, the average Binder cumulant Here, [. ..] dis refers to the disorder average, . . .denotes the Monte Carlo average for each sample, and m = (1/N ) i,τ S i,τ is the order parameter (with N being the number of lattice sites).The finite-size scaling form of the Binder cumulant, depends on two independent arguments.Here, r = (T − T c )/T c is the distance from criticality, z is the dynamical critical exponent, and X gav is the dimensionless scaling function. 1n the absence of a value for the dynamical exponent z, the usual analysis of the Binder cumulant (which identifies the critical point with the crossing of the curves for different L) breaks down because the correct sample shapes are not known.We therefore follow the method devised in Refs.[13,42,45,46].The Binder cumulant has a maximum as function of L τ at fixed L and T .At the critical temperature, the peak position L max τ scales as L z .This allows us to use an iterative approach that finds the value of z together with the critical point.The method begins with a guess of the dynamical exponent and the corresponding sample shapes.The approximate crossing of the g av vs. T curves for these shapes results in an estimate for the critical temperature T c .We then consider g av as a function of L τ for fixed L at this temperature.The peak positions L max τ provide improved estimates for the optimal sample shapes.We iterate these steps three or four times to converge the values of T c and z with reasonable accuracy.
Once the optimal sample shapes are known (which fixes the scaling combination L τ /L z ), the finite-size scaling analysis of other thermodynamic observables such as the order parameter m and its susceptibility χ proceeds normally.Their finite-size scaling forms read where X m and X χ are dimensionless scaling functions, and β and γ are the order parameter and susceptibility critical exponents, respectively.Correlation lengths in space and imaginary time directions are calculated, as usual, from the second moment of the spin correlation function [47,48,49].The finite-size scaling forms of the reduced correlation lengths are given by 4. Results

Phase diagram
The phase diagram of the classical Hamiltonian (2) in the dilution-(classical) temperature plane that arises from our Monte Carlo simulations is presented in Fig. 1.The generic phase boundary (for p < p c ) was determined in Ref. [34] by performing simulations at dilutions p = 0, 1/8, 1/5, 2/7, 1/3, and 9/25.In that work, we also confirmed the location of the percolation phase transition by carrying out simulations at p = p c and T = 1.0.The position of the multicritical point follows from our present simulations (see below); it is located at p = p c = 0.407253 and T = T * = 1.337.Qualitatively, this phase diagram is very similar to the phase diagrams of the diluted quantum magnets discussed in the introductory section.Long-range (superfluid) order is restricted to dilutions p ≤ p c and sufficiently weak quantum fluctuations (which are represented by the classical temperature T after the mapping of the rotor Hamiltonian onto the classical XY model).Importantly, as for the quantum magnets, long-range order survives on the critical percolation cluster (at p = p c ) up to a nonzero classical temperature T * .This causes the existence of two distinct phase transitions (generic and percolation) separated by the multicritical point that is the focus of the present paper.Note that the disordered (insulating) phase can be decomposed into two regions.If the classical temperature is below the transition temperature T c (0) of the undiluted system, large spatial regions without vacancies can show local superfluid order even if the bulk is insulating.This is the Griffiths phase of the transition, i.e., the Mott glass.For classical temperatures above T c (0), local order is impossible; the system is thus a conventional Mott insulator.

Generic and percolation transitions
In this subsection, we briefly summarize the critical behaviors of the generic and percolation transitions.To find the critical exponents of the generic transition, we analyzed [34] Monte Carlo data for dilutions p = 1/8, 1/5, 2/7, 1/3, and 9/25.By including correction-to-scaling terms that account for deviations from pure power-law behavior, we could fit all these data using universal, dilution-independent exponents.The values of the dynamical exponent z, the scale dimension β/ν of the order parameter, the corresponding susceptibility exponent γ/ν, and the correlation length exponent ν are listed in Table 1.Other exponents can be found using various scaling relations.
The critical behavior of the transition across the lattice percolation threshold was analyzed earlier by means of a scaling theory [12] that relates its exponents to the classical percolation exponents (which are known exactly in two space dimensions).The resulting exponent values Table 1.Critical behavior of the superfluid-Mott glass quantum phase transitions.The exponents of the generic (p < p c ) transition were found numerically in Ref. [34].The exponents of the percolation transition were determined from a scaling theory [12] and numerically confirmed in Ref. [34].Their values are exact.The multicritical exponents are the result of the present paper.  1 as well.Our Monte Carlo data for the percolation transition (taken at p = p c and T = 1.0) agree with these predictions [34].

Multicritical point
We now turn to the main topic of this paper, viz., the multicritical point that separates the generic and percolation transitions.To find the multicritical temperature T * and the dynamical exponent z, we employ the iterative procedure described in Sec., m, and χ vs. L for samples of the optimal shape at the multicritical temperature T * = 1.337.The statistical errors are well below the symbol size.The lines represent power-law fits (see text).They are dotted in the regions not included in the fit.temperature of T * = 1.337.(The error ∆T * ≈ 0.003 to 0.005 is estimated heuristically from the quality of the crossings of the g av vs. T curves for optimally shaped samples.)The figure demonstrates that the maximum Binder cumulant g max av for each of the curves does not depend on L (at least for the larger L), it takes a universal value of about 0.6171(3).This is the behavior expected at a (multi)critical point.The deviations for smaller L can be attributed to finite-size effects, i.e., to corrections to scaling.The scaling plot presented in the main panel of Fig. 2 confirms that g av fulfills the scaling form (4) with high precision (slight deviations for small L and L τ stem, again, from corrections to scaling).
Once the optimal shapes are known, the dynamical critical exponent z can be found by analyzing the dependence of L max τ on L, as shown in Fig. 3.The data can be fitted well using a pure power-law dependence L max τ ∼ L z giving a value z = 1.725 (3).The statistical errors of L max τ are determined from 1000 synthetic data sets.Each set is obtained from the original data by adding Gaussian random noise corresponding the uncertainties of the original data.The power-law fit is of good quality, it yields χ2 ≈ 0.74.[The reduced sum of squared errors of the fit (per degree of freedom) is denoted by χ2 to distinguish it from the order parameter susceptibility χ.]In addition to the statistical error, z also has an error stemming from the uncertainty in T * .To estimate it, we repeat the L max τ vs. L analysis at the higher temperature T * + 0.01, significantly further away from T * than our estimated T * error.This leads to a shift of z of about 0.025.Our final estimate for the dynamical critical exponent thus reads z = 1.72 (2).
The critical exponents β/ν and γ/ν can be determined from the L-dependence of the order parameter and the susceptibility, respectively, for samples of optimal shape at the multicritical temperature.The corresponding plots are also shown in Fig. 3.The order parameter data for L > 20 can be fitted well by a pure power law m ∼ L −β/ν giving β/ν = 0.411(1) with χ2 ≈ 0.99.If the leading correction-to-scaling term is included via m = aL −β/ν (1 + bL −ω ), the fit range can be extended to all L > 10, yielding an almost unchanged exponent value β/ν = 0.406(4) with χ2 ≈ 1.26.The largest source of error is again the uncertainty in T * .To estimate it, we repeat the analysis at T * ± ∆T * with ∆T * = 0.003.This leads to shifts of β/ν of about 0.02.Our final estimate therefore is β/ν = 0.41 (2).
The susceptibility data show stronger deviations from pure power-law behavior.A fit to χ ∼ L γ/ν works for sizes L > 36 and yields γ/ν = 2.905(4) with χ2 ≈ 0.66.Extending the fit range to sizes L > 28 gives γ/ν = 2.880(3) but leads to unacceptably large χ2 values.The fit becomes unstable when a correction-to-scaling term is included via χ = aL γ/ν (1 + bL −ω ).Due to these instabilities γ/ν has a larger error.Our final value is γ/ν = 2.90 (5) where the error bar is estimated heuristically from the robustness of the fit; it includes the error due to the uncertainty in T * obtained from repeating the analysis at T * ± 0.003.

Conclusions
In conclusion, we have reported the results of large-scale Monte Carlo simulations of a sitediluted particle-hole symmetric quantum rotor model in two dimensions.They were aimed at investigating the quantum phase transitions between the superfluid and insulating (Mott glass) ground states of disordered interacting bosons.The ground state phase diagram in the dilution-quantum fluctuation plane features two distinct quantum phase transitions, (i) the generic transition that occurs for dilutions below the lattice percolation threshold and is driven by quantum fluctuations, and (ii) the transition across the lattice percolation threshold which is driven by geometric fluctuations.The two phase transition lines meet at a multicritical point at which both quantum and geometric fluctuations become critical.
In the absence of disorder, the superfluid-Mott insulator transition belongs to the threedimensional classical XY universality class.Its correlation length exponent ν ≈ 0.6717 violates the Harris criterion dν > 2. The critical exponents of the diluted system are therefore expected to differ from the clean ones.This expectation is fulfilled by our Monte Carlo results, summarized in Table 1, for the generic transition, the percolation transition, and the multicritical point.In all cases, the transition is of conventional "finite-disorder" type (featuring power-law dynamical scaling) rather than an infinite randomness critical point [50,51,52,53] or a smeared transition [54,55].
We have considered here the two-component quantum rotor Hamiltonian that is a model of disordered interacting bosons.The analogous three-component quantum rotor Hamiltonian which arises as an effective low-energy theory of a bilayer Heisenberg antiferromagnet was studied in Refs.[13,42].Both systems have very similar phase diagrams, but their exponent values differ.For example, the dynamical exponent z of the two-component model considered here takes the values 1.52 and 1.72 for the generic transition and the multicritical point, respectively.The corresponding values for the three-component case are 1.31 and 1.54.Other exponents show similar differences.Note, however, that the exponents of the percolation quantum phase transition do not depend on the number of rotor components, as they are determined by the geometric criticality of the lattice only.
It is also interesting to compare the influence of the space dimensionality.The onedimensional equivalent of the rotor Hamiltonian (1) undergoes a Kosterlitz-Thouless quantum phase transition [56] in the two-dimensional classical XY universality class in the absence of disorder.This transition fulfills the Harris criterion (as ν = ∞), weak disorder therefore does not change the critical behavior.At larger disorder, the character of the transition does change even though the details of its critical behavior are still controversially discussed (see, e.g., Refs.[57,58,59,60] and references therein).In contrast, the clean rotor model in three space dimensions violates the Harris criterion because its correlation length exponent takes the meanfield value 1/2.We therefore expect a scenario similar to the two-dimensional case considered in the present paper.
Experimentally, there are several ways to realize diluted systems of interacting bosons.These include granular superconductors (whose superconductor-insulator quantum phase transition is believed to be bosonic in nature, at least in some systems), ultracold atoms (where the dilution could perhaps be realized via mixtures of heavy and light atoms), or certain magnetic systems such as diluted anisotropic spin-1 antiferromagnets [61].The last example is particularly promising because the required particle-hole symmetry appears naturally as a consequence of the up-down symmetry of the spin Hamiltonian in the absence of an external magnetic field.

Figure 1 .
Figure 1.(Color online) Phase diagram of the (2+1)-dimensional XY model (2) as a function of dilution p and classical temperature T .The black dots mark the numerically determined critical points on the generic and percolation phase boundaries.The red star denotes the multicritical point (MCP) at p = p c = 0.407253 and T = T * = 1.337 that separates the generic and percolation transitions.The lines are guides for the eye only.

3 .Figure 2 .Figure 3 .
Figure 2. (Color online) Binder cumulant g av as function of L τ for several L at the multicritical temperature T * = 1.337.The relative statistical error of g av is between 0.05% and 0.08%.Inset: Raw data g av vs. L τ .Main panel: Scaling plot g av /g max av vs. L τ /L max τ .