Quantum critical behavior of the superfluid-Mott glass transition

# Quantum critical behavior of the superfluid-Mott glass transition

Thomas Vojta Department of Physics, Missouri University of Science and Technology, Rolla, Missouri 65409, USA    Jack Crewse Department of Physics, Missouri University of Science and Technology, Rolla, Missouri 65409, USA    Martin Puschmann Department of Physics, Missouri University of Science and Technology, Rolla, Missouri 65409, USA Institute of Physics, Technische Universität Chemnitz, 09107 Chemnitz, Germany    Daniel Arovas Department of Physics, University of California, San Diego, La Jolla, California 92093, USA    Yury Kiselev Department of Physics, University of California, San Diego, La Jolla, California 92093, USA
July 3, 2019
###### Abstract

We investigate the zero-temperature superfluid to insulator transitions in a diluted two-dimensional quantum rotor model with particle-hole symmetry. We map the Hamiltonian onto a classical -dimensional XY model with columnar disorder which we analyze by means of large-scale Monte Carlo simulations. For dilutions below the lattice percolation threshold, the system undergoes a generic superfluid-Mott glass transition. In contrast to other quantum phase transitions in disordered systems, its critical behavior is of conventional power-law type with universal (dilution-independent) critical exponents , , , , and . These values agree with and improve upon earlier Monte-Carlo results [Phys. Rev. Lett. 92, 015703 (2004)] while (partially) excluding other findings in the literature. As a further test of universality, we also consider a soft-spin version of the classical Hamiltonian. In addition, we study the percolation quantum phase transition across the lattice percolation threshold; its critical behavior is governed by the lattice percolation exponents in agreement with recent theoretical predictions. We relate our results to a general classification of phase transitions in disordered systems, and we briefly discuss experiments.

###### pacs:
05.30.Jp, 64.60.Cn, 74.81.-g, 67.85.Hj

## I Introduction

Zero-temperature phase transitions between superfluid and insulating ground states in systems of disordered interacting bosons are prototypical quantum phase transitions with experimental applications ranging from helium absorbed in vycor Crooker et al. (1983); Reppy (1984) to Josephson junction arrays van der Zant et al. (1992, 1996), superconducting films Haviland et al. (1989); Hebard and Paalanen (1990), doped quantum magnets in high fields Oosawa and Tanaka (2002); Hong et al. (2010); Yu et al. (2012), and to ultracold atoms in disordered optical lattices White et al. (2009); Krinner et al. (2013); D’Errico et al. (2014).

For generic disorder, the two bulk phases, viz. superfluid and Mott insulator, are separated by another phase, the Bose glass which is a compressible gapless insulator Fisher and Fisher (1988); Fisher et al. (1989); Pollet et al. (2009). It can be understood as the Griffiths phase Griffiths (1969); Thill and Huse (1995); Vojta (2006) of the superfluid-insulator transition in which rare large regions of local superfluid order coexist with the insulating bulk. The quantum phase transition between superfluid and Bose glass has been studied in great detail using various analytical and computational techniques. It has recently reattracted considerable attention because new analytical Weichman and Mukhopadhyay (2007) and numerical Priyadarshee et al. (2006); Meier and Wallin (2012); Ng and Sørensen (2015); Álvarez Zúñiga et al. (2015) findings have challenged the scaling relationFisher and Fisher (1988); Fisher et al. (1989) between the dynamical exponent and the space dimensionality (Refs. Weichman and Mukhopadhyay, 2007; Priyadarshee et al., 2006; Meier and Wallin, 2012; Ng and Sørensen, 2015; Álvarez Zúñiga et al., 2015 also contain long lists of references to earlier work.)

In the presence of particle-hole symmetry, the glassy Griffiths phase between superfluid and Mott insulator has a different character: it is the incompressible gapless Mott glass (also called the random-rod glass) Giamarchi et al. (2001); Weichman and Mukhopadhyay (2008). The quantum phase transition between superfluid and Mott glass has attracted less attention than the Bose glass transition. Moreover, the available quantitative results for two space dimensions do not agree with each other. Monte Carlo simulations of a link-current model Prokof’ev and Svistunov (2004) yielded a dynamical critical exponent and a correlation function exponent . 111Here, the numbers in parentheses are the errors of the last digits. A numerical strong-disorder renormalization group study of a particle-hole symmetric quantum rotor model gave , a correlation length exponent , and where is the order parameter susceptibility exponent Iyer et al. (2012). The Fisher relation then implies . Furthermore, a recent Monte Carlo study of a quantum rotor modelSwanson et al. (2014) reported good scaling by setting to its clean value which resulted in . All these models are expected to be in the same universality class. The critical behavior of the superfluid-Mott glass quantum phase transition in two dimensions must thus be considered an open question.

To address this question, we consider a site-diluted two-dimensional quantum rotor model with particle-hole symmetry. After mapping this Hamiltonian onto a classical -dimensional XY model with columnar defects, we perform large-scale Monte Carlo simulations for lattices with up to 11 million sites, averaging over to disorder configurations. The data are analyzed by a finite-size scaling techniqueGuo et al. (1994); Rieger and Young (1994); Sknepnek et al. (2004); *VojtaSknepnek06 that does not require prior knowledge of the dynamical exponent . We also include the leading corrections to scaling. Our results can be summarized as follows: The system features two distinct quantum phase transitions. For dilutions below the percolation threshold of the lattice, we find a superfluid-Mott glass transition characterized by universal (dilution-independent) critical behavior with exponent values , , , , and . The transition across the lattice percolation threshold falls into a different universality class. Its simulation data can be fitted well with the theory developed in Ref. Vojta and Schmalian, 2005a which yields critical exponents that can be expressed in terms of the classical percolation exponents and take the rational values , , , and .

The rest of the paper is organized as follows. Section II introduces the quantum rotor Hamiltonian, the mapping to the classical XY model, and the finite-size scaling technique. Monte Carlo simulations for both the generic () transition and the percolation transition are discussed in Sec. III. We conclude in Sec. IV.

## Ii Theory

### ii.1 Diluted rotor model

The starting point is a site-diluted quantum rotor model on a square lattice given by the Hamiltonian

 H=U2∑iϵi(^ni−¯ni)2−J∑⟨ij⟩ϵiϵjcos(^ϕi−^ϕj) . (1)

Here, is the number operator at site , is the phase operator, and and represent the charging energy and the Josephson coupling, respectively. is the offset charge at site . In the Josephson term, refers to pairs of nearest neighbors. The quenched random variables implement the site dilution. They are independent of each other and take the values 0 (vacancy) with probability and 1 (occupied site) with probability .

As we are interested in the superfluid-Mott glass transition, we set all offset charges to zero and consider commensurate (integer) filling . In this case, the disorder is purely off-diagonal, and the model is particle-hole symmetric. The qualitative features of its phase diagram are well understood Fisher et al. (1989); Weichman and Mukhopadhyay (2008). If the charging energy dominates, , the ground state is a Mott insulator. In the opposite limit, , the ground state is a superfluid as long as the dilution is below the lattice percolation threshold . For , the lattice consists of disconnected clusters and a long-range ordered superfluid state is impossible.

In the case of particle-hole symmetry, the quantum rotor model (1) can be mapped Wallin et al. (1994) onto a classical -dimensional XY model on a cubic lattice having the Hamiltonian

 Hcl=−Js∑⟨i,j⟩,tϵiϵjSi,t⋅Sj,t−Jτ∑i,tϵiSi,t⋅Si,t+1 (2)

where is an O(2) unit vector at the lattice site with spatial coordinate and “imaginary time” coordinate . The coupling constants and are determined by the original quantum rotor Hamiltonian (1) with being an effective “classical” temperature, not equal to the real physical temperature. (The physical temperature of the quantum system (1) maps onto the inverse system size in imaginary time direction of the classical model.) Due to universality, the exact values of and are not important for the critical behavior. We therefore set and drive the XY model (2) through the transition by varying the classical temperature . Because the vacancy positions do not depend on the imaginary time coordinate , the defects in the classical model (2) are columnar, i.e., the disorder is perfectly correlated in the imaginary time direction (see Fig. 1).

In the clean undiluted limit , the Hamiltonian (2) simplifies to the usual three-dimensional XY model. The correlation length critical exponent of the three-dimensional XY universality class takes the value (see, e.g., Ref. Campostrini et al., 2006). This value violates the Harris criterion Harris (1974) where is the number of dimensions in which there is randomness. Consequently, the three-dimensional clean XY critical point is unstable against columnar defects, and we expect the diluted system to feature a different critical behavior.

### ii.2 Anisotropic finite-size scaling

Finite-size scaling Barber (1983); Cardy (1988) is a powerful tool for analyzing Monte Carlo data. Particularly useful are quantities of scale dimension zero such as the (average) Binder cumulant

 gav=[1−⟨|m|4⟩3⟨|m|2⟩2]dis, (3)

where is the order parameter ( denotes the number of lattice sites). refers to the disorder average and denotes the Monte Carlo average for each sample. In an isotropic system with a single relevant length scale, it takes the scaling form . Here is the linear system size, is the distance from criticality, and is a scaling function. This scaling form implies that vs.  curves for systems of different sizes all cross at criticality, , having the value . This can be used to find the critical point with high accuracy. Moreover, the slopes of the vs.  curves at vary as which can be used to measure .

As the quenched disorder in our Hamiltonian (2) breaks the symmetry between the space and imaginary time directions, we need to distinguish the linear system size in the two space directions from the size in the imaginary time direction. ( corresponds to the inverse physical temperature of the original quantum model (1).) If the putative disordered critical point fulfills conventional power-law dynamical scaling, the finite-size scaling form of the average Binder cumulant then reads

 gav(r,L,Lτ)=Xgav(rL1/ν,Lτ/Lz) (4)

where is the dynamical critical exponent, and is the dimensionless scaling function which now depends on two arguments. Note that some quantum phase transitions in disordered systems feature exotic activated dynamical scaling instead of power-law scaling, for example the ferromagnetic transition in the random transverse-field Ising model Fisher (1992); *Fisher95, the pairbreaking superconductor-metal quantum phase transition Hoyos et al. (2007); *VojtaKotabageHoyos09; Del Maestro et al. (2008); *DRHV10; Xing et al. (2015), and magnetic transitions in itinerant systems Ubaid-Kassis et al. (2010); Vojta (2010). For activated dynamical scaling, the scaling combination in Eq. (4) needs to be replaced by where is the tunneling exponent. Based on the classification of disordered quantum phase transitions developed in Refs. Vojta and Schmalian, 2005b; Vojta, 2006, we do not expect the superfluid-Mott glass transition to show activated scaling. We will return to this point in the concluding section.

How can one perform a finite-size scaling analysis of Monte Carlo data based on the scaling form (4) of the average Binder cumulant? If the value of is known, the analysis is as simple as in the isotropic case: One chooses system sizes and such that were is a constant. Then the vs.  curves for systems of different sizes cross at criticality [with the value ] which can be used to locate the critical point. However, in the absence of a value for , this approach breaks down because the correct shapes (aspect ratios) of the samples are not known.

A method for finding the correct sample shape within the simulations Guo et al. (1994); Rieger and Young (1994); Sknepnek et al. (2004); *VojtaSknepnek06 can be based on the following property of the Binder cumulant: For fixed , as a function of has a peak at position and value . The peak position marks the optimal sample shape, where the ratio behaves like the corresponding ratio of the correlation lengths in time and space directions, . (If the aspect ratio deviates from the optimal one, the system can be decomposed into independent units either in space or in time direction, and thus decreases.) At criticality, must be proportional to , fixing the second argument of the scaling function . This implies that the peak value at criticality is independent of and that the vs.  curves of samples of the optimal shape () cross at .

In our simulations, we use an iterative approach. We start from a guess for and the corresponding sample shapes. The approximate crossing of the vs.  curves for these samples gives an estimate for . At this temperature, we next analyze as a function of for fixed . The values of give improved estimates for the optimal sample shapes and thus for . After iterating this procedure three or four times, the values of and will have converged with reasonable accuracy.

Once and are determined, the finite-size scaling analysis proceeds as usual, based on the scaling forms

 m = L−β/νXm(rL1/ν,Lτ/Lz) , (5) χ = Lγ/νXχ(rL1/ν,Lτ/Lz) (6)

for the order parameter and its susceptibility . Here, and are dimensionless scaling functions, and and are the order parameter and susceptibility critical exponents, respectively.

In addition to these thermodynamic quantities, we also calculate the correlation lengths and is the space and imaginary time directions, respectively. They are obtained, as usual, from the second moment of the spin-spin correlation function Cooper et al. (1982); Kim (1993); Caracciolo et al. (2001) and can be expressed in terms of the Fourier transform of the correlation function,

 ξs = ⎡⎢⎣(~G(0,0)−~G(qs0,0)q2s0~G(qs0,0))1/2⎤⎥⎦dis , (7) ξτ = ⎡⎢⎣(~G(0,0)−~G(0,qτ0)q2τ0~G(0,qτ0))1/2⎤⎥⎦dis . (8)

Here, and are the minimum values of the wave numbers and that fit into a system of linear size and in space and imaginary time direction, respectively. The reduced correlation lengths and have scale dimension zero, their scaling forms therefore read

 ξs/L = Xξs(rL1/ν,Lτ/Lz) , (9) ξτ/Lτ = Xξτ(rL1/ν,Lτ/Lz) . (10)

## Iii Monte Carlo simulations

### iii.1 Overview

Our Monte Carlo simulations of the classical XY model (2) combine the Wolff cluster algorithm Wolff (1989) with conventional Metropolis updates Metropolis et al. (1953). Specifically, a full Monte Carlo sweep consists of a Metropolis sweep over the lattice followed by a Wolff sweep. (A Wolff sweep is defined as a number of cluster flips such that the total number of flipped spins equals the number of lattice sites.) The Wolff algorithm greatly reduces the critical slowing down, and the Metropolis updates equilibrate small disconnected clusters of sites that are missed in the construction of the Wolff clusters (this becomes important at higher dilutions ).

We simulate systems with linear sizes up to in space direction and up to in the imaginary time direction at dilutions , 1/8, 1/5, 2/7, 1/3, 9/25 and the percolation threshold .

The simulation of disordered systems requires a high numerical effort because many samples with different disorder configurations need to be studied to compute averages, variances, and distributions of observables. For good performance, one must thus carefully optimize the number of samples (i.e., disorder configurations) and the number of measurements during the simulation of each sample. Based on the consideration in Refs. Ballesteros et al., 1998a, b; Sknepnek et al., 2004; *VojtaSknepnek06; Zhu et al., 2015, we have chosen rather short runs of full sweeps per sample (with a measurement after each sweep) but large numbers of disorder configurations ranging from to 50 000 depending on the system size. The equilibration period is taken to be 100 full sweeps, significantly longer than the actual equilibration times that reach 30 to 40 sweeps at maximum. Short Monte Carlo runs can lead to biases in some of the observables. To eliminate these, we have implemented improved estimators along the lines discussed in the appendix of Ref. Zhu et al., 2015.

The phase diagram resulting from these simulations is shown in Fig. 2.

The critical temperature decreases with increasing dilution from its clean value , as expected. For dilutions above the percolation threshold , the lattice consists of disconnected finite-size clusters. Therefore, long-range superfluid order is impossible. Right at , there is an infinite cluster of dimension where is the dimensionality of the critical percolation cluster in two dimensions, and the extra 1 stems from the imaginary time direction. As is larger than the lower critical dimension of the XY model, the XY model on the critical percolation cluster orders below a multicritical temperature . This implies that the phase boundary coincides with the classical percolation threshold for (see also Ref. Vojta and Hoyos, 2008). We thus identify two different phase transitions, (i) the generic superfluid-Mott glass transition for and (ii) a percolation transition across the lattice percolation threshold.

In the following sections, we discuss the critical behaviors of these transitions in detail. To test our codes, we have also studied the clean limit using system sizes up to sites. By analyzing the crossings of the Binder cumulant and the reduced correlation length, we find a critical temperature . Finite-size scaling then gives the critical exponents , , and . Within their errors, they agree well with high-precision results for the three-dimensional XY universality class Campostrini et al. (2006).

As a further test for the universality of the (generic) critical behavior, we also perform exploratory simulations of a soft-spin version of the classical Hamiltonian. They are discussed in Sec. III.4.

### iii.2 Generic superfluid-Mott glass transition

To analyze the critical behavior of the generic transition occurring for , we consider five different dilutions, , 1/5, 2/7, 1/3, and 9/25. As described in Sec. II.2, we use an iterative procedure that consists of two types of simulation runs. The first are runs right at for systems with several different for each . Finite-size scaling of the Binder cumulant at as a function of and gives the optimal sample shapes and the dynamical exponent . In the second set of simulations, we vary the temperature over a range in the vicinity of , but we consider only the optimal shapes found in the first part. Finite-size scaling of the order parameter, susceptibility, Binder cumulant, and correlation length as functions of and then yields the critical exponents , , and .

The inset of Fig. 3 shows the Binder cumulant as a function of for several to 100 at the estimated critical temperature for dilution .

As expected at the critical point, the maximum Binder cumulant for each of the curves does not depend on . (The remaining weak variation visible in the figure can be attributed to corrections to scaling, see below.) To generate a scaling plot that tests the scaling form (4), we now fit each vs.  curve with an inverted parabola in . The vertex of this parabola yields the position of the maximum and its value . When plotting vs.  the data scale very well, as can be seen in the resulting scaling plot in the main panel of Fig. 3. This demonstrates that the Binder cumulant fulfills Eq. (4) with high accuracy. We have created the corresponding scaling plots for all the other dilutions, , 1/5, 2/7, and 9/25, with analogous results. 222For low dilutions , the parabola fits of vs.  are affected by corrections to scaling for small and . We thus slightly adjust and to further improve the quality of the data collapse onto a common master curve. This applies to the four smallest system sizes for and 1/5 and the three smallest sizes for . The resulting change of the value of is about 0.01, well below the error due to the uncertainty in .

To determine the dynamical critical exponent , we now analyze the dependence of the positions of the maximum on . According to Eq. (4), we expect the power-law dependence . In Fig. 4, we plot vs.  for all dilutions .

The curves show significant deviations from pure power-law behavior, in particular for the smaller dilutions, indicating that the crossover from clean to disordered critical behavior is slow. The resulting corrections to scaling are strong and cannot be neglected. Pure power-law fits of the data would therefore only yield effective, scale-dependent exponents. To determine the true asymptotic exponents, we include the leading corrections to scaling via the ansatz with universal (dilution-independent) critical exponents and but dilution-dependent prefactors and . The exponent values resulting from a combined fit of the data for all five dilutions are and . The fit is of good quality giving . [We denote the reduced sum of squared errors of the fit (per degree of freedom) by to distinguish it from the susceptibility .] The fit is also robust against removing complete data sets or removing points from the upper or lower end of each set. Interestingly, the leading corrections to scaling appear to vanish somewhere between and 9/25, as the prefactor of the correction term changes sign. Correspondingly, pure power-law fits of the and 9/25 data yield and 1.546, respectively. These values are close to the estimate from the combined fit and nicely bracket it on both sides. An additional significant source of errors is the uncertainty of the critical temperature. To assess its effect on the dynamical exponent, we repeat the vs.  analysis (for dilutions and 9/25) at temperatures slightly above and below our estimated (, roughly at the boundaries of our confidence intervals). This leads to shifts in of about 0.01 to 0.02. Our final estimate for the dynamical critical exponent therefore reads .

To find the remaining critical exponents, we now turn to the Monte Carlo runs that use the optimal sample shapes . According to Eqs. (5) and (6), and can be obtained from the dependence of the order parameter and susceptibility at of the optimally shaped samples. As we expect corrections to scaling to be important, we again include subleading terms in our fit functions, for the order parameter and for the susceptibility. Here , , and are the universal, dilution-independent critical exponents while the coefficients and are again non-universal. (Note that and generally differ from quantity to quantity; we use the same symbols to avoid cluttering up the notation too much.) When performing fits of our data to these expressions, we noticed, however, that the quality of the fits is extremely sensitive to small changes of the estimates for (much more so than in the analysis of the dynamical exponent above). To determine higher accuracy estimates of , we use the criterion that the value of at criticality should approach a dilution-independent constant with at a universal critical point. Varying until this criterion is fulfilled yields improved estimates for the critical temperatures, viz.  for , for , for , for , and for . We estimate the error of these values to be about . Figure 5 shows the resulting dependence on .

In the large- limit, approaches the value 0.599(2). Note that the non-monotonic behavior of for weak dilutions suggests that at least two corrections to scaling terms contribute at small .

Using the improved critical temperatures, we now proceed to determine and . Figure 6 shows the order parameter at as a function of for all dilutions .

The combined fit of all data to is of good quality () if the smallest system sizes are excluded (see figure). Interestingly, the sizes that need to be excluded are exactly those for which in Fig. 5 appears to be dominated by the second subleading correction to scaling term.) The exponents resulting from the fit read and . To assess the error arising from the uncertainty in , we repeat the analysis for temperatures with . This leads to shifts of of about 0.01. Our final estimate therefore reads .

The system-size dependence of the order parameter susceptibility at criticality is presented in Fig. 7 for all dilutions .

After excluding the smallest system sizes (see figure), the combined fit of all data to is again of good quality () and yields the exponents and . After including potential errors from the uncertainty in and the fit range, the final exponent estimate is .

So far, the analysis has focused on the behavior right at . To find a complete set of critical exponents, we now determine the correlation length exponent which requires off-critical data. Figure 8 shows the temperature dependence of the Binder cumulant and the reduced correlation length for systems of optimal shape but different sizes at dilution .

Both quantities have scale dimension zero, therefore, the curves for different system sizes are expected to cross at the critical temperature . The figure demonstrates that the crossings for both quantities shift with increasing , reflecting significant corrections to scaling. According to Eqs. (4) and (8), the slopes and at the critical temperature vary as with system size. To extract the slopes, we fit straight lines (for ) or quadratic parabolas (for ) to the data close to . The resulting slopes are shown as a function of system size in Figs. 9 and 10, respectively.

The exponent is now obtained from fits of the slopes to the form . In the case of the reduced correlation length (Fig. 9) a combined fit of all dilutions is of good quality after the smallest system sizes have been excluded ( and yields as well as . The corresponding fit of the slopes of the Binder cumulant has a somewhat poorer quality ( and is not very stable with respect to adding and removing data points at the ends of the interval. The resulting exponents and have therefore larger errors. In addition to the slopes of the Binder cumulant and the reduced correlation length at , we have also studied the slopes of and (not shown). After we account for the differences between all these estimates and include potential errors from the uncertainty in (by repeating the analysis at temperatures ) we arrive at the final estimate . This value fulfills the inequality Chayes et al. (1986) .

The critical exponents , , and are not independent of each other as they must fulfill the hyperscaling relation where is the space dimensionality. Our values, , , and fulfill this relation within their error bars. We also note that all our estimates for the leading irrelevant exponent are roughly consistent with each other, giving us confidence that our results represent true asymptotic rather than effective critical exponents.

### iii.3 Percolation transition

We now turn to the percolation transition that occurs when the system is tuned through the percolation threshold at low (classical) temperatures (see Fig. 2). The critical behavior of this transition stems from the critical geometry of the percolating lattice while the dynamical fluctuations of the rotor variables are uncritical and “just go along for the ride” (the rotor model on each of the percolation clusters is locally ordered). Vojta and Schmalian Vojta and Schmalian (2005a) developed a theory of this percolation quantum phase transition. It predicts critical behavior governed by the lattice percolation exponents. For two space dimensions it yields the exact exponent values , , , and .

To test these predictions, we perform simulations at dilution and temperature . These calculations require a particularly high numerical effort, because the large value of leads to a rapid increase with of the optimal system size in imaginary time direction. We have thus restricted the simulations to sizes up to and using between 10 000 and 50 000 disorder configurations.

The data analysis proceeds in analogy to Sec. III.2. We obtain from the maxima of the Binder cumulant as a function of at fixed . In Fig. 11, we present a plot of vs. .

The data can be fitted with high quality () to the predicted power law . After having found , we calculate the order parameter and susceptibility right at criticality for optimally shaped samples of different sizes. The resulting data are also presented in Fig. 11. The susceptibility data can be fitted well to the predicted power law giving . The exponent is very small, corresponding to a slow decay of the order parameter with . Subleading corrections are thus much more visible as indicated by the curvature of the vs.  curve in Fig. 11. We have therefore fitted the order parameter to . This fit is again of high quality, with .

Our simulation data thus agree nearly perfectly with the critical behavior predicted in Ref. Vojta and Schmalian, 2005a.

### iii.4 Soft-spin model

We also consider a soft-spin version of the classical Hamiltonian to test whether or not its critical exponents agree with those of the hard-spin model analyzed above, as is expected from universality. The soft-spin Hamiltonian reads

 Hsoft = −∑⟨i,j⟩,tϵiϵjSi,t⋅Sj,t−∑i,tϵiSi,t⋅Si,t+1 (11) −12∑i,tϵi|Si,t|2+∑i,tϵi(|Si,t|2)2

where now represents an unrestricted two-component vector. We perform Monte-Carlo simulations of this soft-spin model using the efficient Worm algorithm Prokof’ev and Svistunov (2001), studying dilutions and 0.337. The system sizes range from to 24 with fixed at using the dynamical exponent value found in Sec. III.2 333We actually use which is close to the effective dynamical exponent found for the system size range and dilution considered..

We now analyze the correlation length in imaginary time direction (equivalent to the inverse energy gap of the corresponding quantum model) on the disordered side of the phase transition. According to Eq. (10), its scaling form for samples of shape can be written as . Thus, if we plot vs. , the data for different sizes and temperatures should all fall onto a single master curve. Figure 12 presents such a plot for two site dilutions , with the critical exponents and fixed at the values found in Sec. III.2.

Within their statistical errors, the data scale well. Consequently, even though we have not independently determined the critical exponents of the soft-spin model (11), the Monte Carlo data are compatible with the critical behavior found earlier.

## Iv Conclusions

In summary, we have carried out large-scale computer simulations to determine the critical behavior of the superfluid-Mott glass quantum phase transition in two space dimensions. To this end, we have mapped a quantum rotor model with commensurate filling and off-diagonal disorder onto a (2+1)-dimensional classical XY model with columnar defects. We have then analyzed this classical system by means of Monte Carlo methods.

The corresponding clean superfluid-Mott insulator transition is in the three-dimensional XY universality class; its correlation length exponent violates the Harris criterion with . The clean critical behavior is therefore expected to be unstable against the columnar disorder. Accordingly, we have found that the critical behavior of the superfluid-Mott glass transition differs from that of the clean superfluid-Mott insulator transition.

In contrast to other quantum phase transitions in disordered systems Fisher (1992); *Fisher95; Hoyos et al. (2007); *VojtaKotabageHoyos09; Del Maestro et al. (2008); *DRHV10; Xing et al. (2015); Ubaid-Kassis et al. (2010); Vojta (2010), the superfluid-Mott glass transition features a conventional finite-disorder critical point whose dynamical scaling is characterized by a power-law relation between the correlation lengths in the space and time directions (rather than an infinite-randomness critical point with activated dynamical scaling for which would grow exponentially with ). This result agrees with the general classification of phase transitions in disordered systems based on the rare region (or defect) dimensionality Vojta and Schmalian (2005b); Vojta (2006). In terms of the mapped, classical Hamiltonian (2), the rare regions in our problem are one-dimensional rods with XY order-parameter symmetry. As the lower-critical dimension of the classical XY model is , the rare region dimensionality fulfills , putting the system into the conventional class A of the classification.

For the generic transition occurring at dilutions below the lattice percolation threshold , our Monte Carlo data are described well by a universal critical behavior with dilution-independent critical exponents. The numerical estimates of the exponent values are summarized in Table 1 and compared to earlier results in the literature.

Our results are in reasonable agreement with (but more accurate than) Monte Carlo simulations of a link-current model Prokof’ev and Svistunov (2004) that is expected to be in the same universality class as our Hamiltonian. The results in Ref. Iyer et al., 2012 were obtained using a numerical implementation of the strong-disorder renormalization group. This method is expected to give approximate rather than exact results at a conventional finite-disorder critical point such as the one under consideration here. In view of this, the agreement of and can be considered satisfactory. However, the values of , , and (that all involve the scale dimension of the order parameter) are far away from the Monte Carlo results in this work and in Ref. Prokof’ev and Svistunov, 2004. Our findings are also incompatible with the clean value that was assumed in Ref. Swanson et al., 2014.

It is interesting to consider the evolution of the dynamical exponent with the order parameter dimensionality. The deviation of from the clean value, which is for any number of components, can be understood as a measure of the strength of the disorder effects. In the (2+1)-dimensional Heisenberg model (three order parameter components) with columnar defects, the exponent takes the value Sknepnek et al. (2004); Vojta and Sknepnek (2006) . The (2+1)-dimensional XY model (two components) studied in the present paper has , while the corresponding Ising model Motrunich et al. (2000) (one component) features activated scaling that corresponds to . The value of thus increases monotonically with decreasing order parameter dimensionality.

In addition to the generic superfluid-Mott glass transition that occurs for dilutions , we have also investigated the percolation quantum phase transition across . Here, our Monte Carlo data agree very well with the predictions of the scaling theory by Vojta and Schmalian Vojta and Schmalian (2005a).

Potential routes to study the superfluid-Mott glass transition in experiment include disordered bosonic systems in ultracold atoms as well as dirty and granualar superconductors (for some superconductor-insulator transitions, there is experimental and numerical evidence for the bosonic nature of the transition). In these systems, it may be hard, though, to fulfill the condition of exact particle-hole symmetry in the presence of disorder. Statistical particle hole symmetry may be easier to achieve, but it is not fully resolved whether or not it would destabilize the Mott glass and turn it into a Bose glass Weichman and Mukhopadhyay (2008); Altman et al. (2010); Wang et al. (2015).

Another type of experimental systems that contain Mott-glass physics are diluted anisotropic spin-1 antiferromagnets Roscilde and Haas (2007). In this case, the particle-hole symmetry appears naturally as it is a consequence of the up-down symmetry of the spin Hamiltonian in the absence of an external magnetic field. Such a magnetic realization of a Mott glass (albeit in three dimensions) was recently observed in bromine-doped dichloro-tetrakis-thiourea-nickel (DTN) Yu et al. (2012).

## Acknowledgements

This work was supported in part by the NSF under Grant Nos. DMR-1205803 and DMR-1506152 as well as by funds from the UCSD Academic Senate. M.P. acknowledges support by an InProTUC scholarship of the German Academic Exchange Service. We thank Snir Gazit, Gil Refael, and Nandini Trivedi for helpful discussions.

## References

• Crooker et al. (1983) B. C. Crooker, B. Hebral, E. N. Smith, Y. Takano,  and J. D. Reppy, Phys. Rev. Lett. 51, 666 (1983).
• Reppy (1984) J. D. Reppy, Physica B+C 126, 335 (1984).
• van der Zant et al. (1992) H. S. J. van der Zant, F. C. Fritschy, W. J. Elion, L. J. Geerligs,  and J. E. Mooij, Phys. Rev. Lett. 69, 2971 (1992).
• van der Zant et al. (1996) H. S. J. van der Zant, W. J. Elion, L. J. Geerligs,  and J. E. Mooij, Phys. Rev. B 54, 10081 (1996).
• Haviland et al. (1989) D. B. Haviland, Y. Liu,  and A. M. Goldman, Phys. Rev. Lett. 62, 2180 (1989).
• Hebard and Paalanen (1990) A. F. Hebard and M. A. Paalanen, Phys. Rev. Lett. 65, 927 (1990).
• Oosawa and Tanaka (2002) A. Oosawa and H. Tanaka, Phys. Rev. B 65, 184437 (2002).
• Hong et al. (2010) T. Hong, A. Zheludev, H. Manaka,  and L.-P. Regnault, Phys. Rev. B 81, 060410 (2010).
• Yu et al. (2012) R. Yu, L. Yin, N. S. Sullivan, J. S. Xia, C. Huan, A. Paduan-Filho, N. F. O. Jr, S. Haas, A. Steppke, C. F. Miclea, F. Weickert, R. Movshovich, E.-D. Mun, B. L. Scott, V. S. Zapf,  and T. Roscilde, Nature 489, 379 (2012).
• White et al. (2009) M. White, M. Pasienski, D. McKay, S. Q. Zhou, D. Ceperley,  and B. DeMarco, Phys. Rev. Lett. 102, 055301 (2009).
• Krinner et al. (2013) S. Krinner, D. Stadler, J. Meineke, J.-P. Brantut,  and T. Esslinger, Phys. Rev. Lett. 110, 100601 (2013).
• D’Errico et al. (2014) C. D’Errico, E. Lucioni, L. Tanzi, L. Gori, G. Roux, I. P. McCulloch, T. Giamarchi, M. Inguscio,  and G. Modugno, Phys. Rev. Lett. 113, 095301 (2014).
• Fisher and Fisher (1988) D. S. Fisher and M. P. A. Fisher, Phys. Rev. Lett. 61, 1847 (1988).
• Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein,  and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
• Pollet et al. (2009) L. Pollet, N. V. Prokof’ev, B. V. Svistunov,  and M. Troyer, Phys. Rev. Lett. 103, 140402 (2009).
• Griffiths (1969) R. B. Griffiths, Phys. Rev. Lett. 23, 17 (1969).
• Thill and Huse (1995) M. Thill and D. A. Huse, Physica A 214, 321 (1995).
• Vojta (2006) T. Vojta, J. Phys. A 39, R143 (2006).
• Weichman and Mukhopadhyay (2007) P. B. Weichman and R. Mukhopadhyay, Phys. Rev. Lett. 98, 245701 (2007).
• Priyadarshee et al. (2006) A. Priyadarshee, S. Chandrasekharan, J.-W. Lee,  and H. U. Baranger, Phys. Rev. Lett. 97, 115703 (2006).
• Meier and Wallin (2012) H. Meier and M. Wallin, Phys. Rev. Lett. 108, 055701 (2012).
• Ng and Sørensen (2015) R. Ng and E. S. Sørensen, Phys. Rev. Lett. 114, 255701 (2015).
• Álvarez Zúñiga et al. (2015) J. P. Álvarez Zúñiga, D. J. Luitz, G. Lemarié,  and N. Laflorencie, Phys. Rev. Lett. 114, 155301 (2015).
• Giamarchi et al. (2001) T. Giamarchi, P. Le Doussal,  and E. Orignac, Phys. Rev. B 64, 245119 (2001).
• Weichman and Mukhopadhyay (2008) P. B. Weichman and R. Mukhopadhyay, Phys. Rev. B 77, 214516 (2008).
• Prokof’ev and Svistunov (2004) N. Prokof’ev and B. Svistunov, Phys. Rev. Lett. 92, 015703 (2004).
• (27) Here, the numbers in parentheses are the errors of the last digits.
• Iyer et al. (2012) S. Iyer, D. Pekker,  and G. Refael, Phys. Rev. B 85, 094202 (2012).
• Swanson et al. (2014) M. Swanson, Y. L. Loh, M. Randeria,  and N. Trivedi, Phys. Rev. X 4, 021007 (2014).
• Guo et al. (1994) M. Guo, R. N. Bhatt,  and D. A. Huse, Phys. Rev. Lett. 72, 4137 (1994).
• Rieger and Young (1994) H. Rieger and A. P. Young, Phys. Rev. Lett. 72, 4141 (1994).
• Sknepnek et al. (2004) R. Sknepnek, T. Vojta,  and M. Vojta, Phys. Rev. Lett. 93, 097201 (2004).
• Vojta and Sknepnek (2006) T. Vojta and R. Sknepnek, Phys. Rev. B. 74, 094415 (2006).
• Vojta and Schmalian (2005a) T. Vojta and J. Schmalian, Phys. Rev. Lett. 95, 237206 (2005a).
• Wallin et al. (1994) M. Wallin, E. S. Sorensen, S. M. Girvin,  and A. P. Young, Phys. Rev. B 49, 12115 (1994).
• Campostrini et al. (2006) M. Campostrini, M. Hasenbusch, A. Pelissetto,  and E. Vicari, Phys. Rev. B 74, 144506 (2006).
• Harris (1974) A. B. Harris, J. Phys. C 7, 1671 (1974).
• Barber (1983) M. N. Barber, in Phase Transitions and Critical Phenomena, Vol. 8, edited by C. Domb and J. L. Lebowitz (Academic, New York, 1983) pp. 145–266.
• Cardy (1988) J. Cardy, ed., Finite-size scaling (North Holland, Amsterdam, 1988).
• Fisher (1992) D. S. Fisher, Phys. Rev. Lett. 69, 534 (1992).
• Fisher (1995) D. S. Fisher, Phys. Rev. B 51, 6411 (1995).
• Hoyos et al. (2007) J. A. Hoyos, C. Kotabage,  and T. Vojta, Phys. Rev. Lett. 99, 230601 (2007).
• Vojta et al. (2009) T. Vojta, C. Kotabage,  and J. A. Hoyos, Phys. Rev. B 79, 024401 (2009).
• Del Maestro et al. (2008) A. Del Maestro, B. Rosenow, M. Müller,  and S. Sachdev, Phys. Rev. Lett. 101, 035701 (2008).
• Del Maestro et al. (2010) A. Del Maestro, B. Rosenow, J. A. Hoyos,  and T. Vojta, Phys. Rev. Lett. 105, 145702 (2010).
• Xing et al. (2015) Y. Xing, H.-M. Zhang, H.-L. Fu, H. Liu, Y. Sun, J.-P. Peng, F. Wang, X. Lin, X.-C. Ma, Q.-K. Xue, J. Wang,  and X. C. Xie, Science 350, 542 (2015).
• Ubaid-Kassis et al. (2010) S. Ubaid-Kassis, T. Vojta,  and A. Schroeder, Phys. Rev. Lett. 104, 066402 (2010).
• Vojta (2010) T. Vojta, J. Low Temp. Phys. 161, 299 (2010).
• Vojta and Schmalian (2005b) T. Vojta and J. Schmalian, Phys. Rev. B 72, 045438 (2005b).
• Cooper et al. (1982) F. Cooper, B. Freedman,  and D. Preston, Nucl. Phys. B 210, 210 (1982).
• Kim (1993) J. K. Kim, Phys. Rev. Lett. 70, 1735 (1993).
• Caracciolo et al. (2001) S. Caracciolo, A. Gambassi, M. Gubinelli,  and A. Pelissetto, Eur. Phys. J. B 20, 255 (2001).
• Wolff (1989) U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
• Metropolis et al. (1953) N. Metropolis, A. Rosenbluth, M. Rosenbluth,  and A. Teller, J. Chem. Phys. 21, 1087 (1953).
• Ballesteros et al. (1998a) H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Muñoz Sudupe, G. Parisi,  and J. J. Ruiz-Lorenzo, Phys. Rev. B 58, 2740 (1998a).
• Ballesteros et al. (1998b) H. G. Ballesteros, L. A. Fernandez, V. Martin-Mayor, A. Munoz Sudupe, G. Parisi,  and J. J. Ruiz-Lorenzo, Nucl. Phys. B 512, 681 (1998b).
• Zhu et al. (2015) Q. Zhu, X. Wan, R. Narayanan, J. A. Hoyos,  and T. Vojta, Phys. Rev. B 91, 224201 (2015).
• Vojta and Hoyos (2008) T. Vojta and J. A. Hoyos, in Recent Progress in Many-Body Theories, edited by J. Boronat, G. Astrakharchik,  and F. Mazzanti (World Scientific, Singapore, 2008) p. 235.
• (59) For low dilutions , the parabola fits of vs.  are affected by corrections to scaling for small and . We thus slightly adjust and to further improve the quality of the data collapse onto a common master curve. This applies to the four smallest system sizes for and 1/5 and the three smallest sizes for . The resulting change of the value of is about 0.01, well below the error due to the uncertainty in .
• Chayes et al. (1986) J. T. Chayes, L. Chayes, D. S. Fisher,  and T. Spencer, Phys. Rev. Lett. 57, 2999 (1986).
• Prokof’ev and Svistunov (2001) N. Prokof’ev and B. Svistunov, Phys. Rev. Lett. 87, 160601 (2001).
• (62) We actually use which is close to the effective dynamical exponent found for the system size range and dilution considered.
• Motrunich et al. (2000) O. Motrunich, S. C. Mau, D. A. Huse,  and D. S. Fisher, Phys. Rev. B 61, 1160 (2000).
• Altman et al. (2010) E. Altman, Y. Kafri, A. Polkovnikov,  and G. Refael, Phys. Rev. B 81, 174528 (2010).
• Wang et al. (2015) Y. Wang, W. Guo,  and A. W. Sandvik, Phys. Rev. Lett. 114, 105303 (2015).
• Roscilde and Haas (2007) T. Roscilde and S. Haas, Phys. Rev. Lett. 99, 047205 (2007).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters