# Critical Behavior and Collective Modes at the Superfluid Transition in
Amorphous Systems

Pulloor Kuttanikkad Vishnu  Department of Physics, Indian Institute of
Technology Madras, Chennai 600036, India.  Martin Puschmann  Department of
Physics, Missouri University of Science and Technology, Rolla, Missouri 65409,
USA.  Rajesh Narayanan  Department of Physics, Indian Institute of Technology
Madras, Chennai 600036, India<EMAIL_ADDRESS>Thomas Vojta
Department of Physics, Missouri University of Science and Technology, Rolla,
Missouri 65409, USA.

###### Abstract

We investigate the critical behavior and the dynamics of the amplitude (Higgs)
mode close to the superfluid-insulator quantum phase transition in an
amorphous system (i.e., a system subject to topological randomness). In
particular, we map the two-dimensional Bose-Hubbard Hamiltonian defined on a
random Voronoi-Delaunay lattice onto a (2+1)-dimensional layered classical XY
model with correlated topological disorder. We study the resulting model by
laying recourse to classical Monte Carlo simulations. We specifically focus on
the scalar susceptibility of the order parameter to study the dynamics of the
amplitude mode. To do so, we harness the maximum entropy method to perform the
analytic continuation of the scalar susceptibility to real frequencies. Our
analysis shows that the amplitude mode remains delocalized in the presence of
such topological disorder, quite at odds with its behavior in generic
disordered systems, where the randomness localizes the Higgs mode.
Furthermore, we show that the critical behavior of the topologically
disordered system is identical to that of its translationally invariant
counterpart, consistent with a modified Harris criterion. This suggests that
the localization of the collective excitations in the presence of disorder is
tied to the critical behavior of the quantum phase transition rather than a
simple Anderson-localization-type interference mechanism.

##  I Introduction

It is well established that quantum phase transitions in many-body systems are
deeply impacted by quenched disorder. Research over the past few decades has
shown that random disorder typified by vacancies, defects, and other types of
impurities can trigger a plethora of exciting novel phenomena such as quantum
Griffiths phases [1, 2, 3], infinite-randomness criticality [4, *Fisher95, 6,
*vojta_kotabage_prb_09], other unconventional critical points [8, *AKPR08,
10], as well as smeared transitions [11, *HoyosVojta08]. Classification
schemes [13, 14] for these phenomena have been put forward according to the
behavior of the disorder strength under coarse graining and on the importance
of rare disorder fluctuations (for reviews see, e.g., Refs. [15, *Vojta13,
*Vojta19]). Many of these exotic phenomena have now been realized
experimentally [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28].

The majority of the investigations in the literature on disorder effects at
quantum phase transitions have been restricted to the thermodynamic behavior.
Only recently, attention has focussed on the effects of impurities on the
dynamics of excitations hosted by many-body systems near quantum criticality.
Of particular interest is the role that impurities play for the properties of
the amplitude (Higgs) modes and the Goldstone modes that emerge when a
continuous symmetry is spontaneously broken.

Indeed, the question of disorder effects on the amplitude modes has gathered a
renewed impetus, as these modes can now be experimentally realized in various
magnetic, superconducting, cold atom, and charge-density-wave systems [29,
30]. In the case of clean systems, amplitude modes can be identified as sharp
excitations in the symmetry-broken phase, close to the critical point [31, 32,
33]. More recently, a number of studies investigating the effects of disorder
on collective excitation have come to the fore. A series of papers [34, 35,
36] has revealed that the amplitude mode obtained at a superfluid-Mott glass
transition is localized by the quenched impurities. Furthermore, it has been
shown that, in accordance with the Goldstone theorem, only the lowest-lying
Goldstone mode remains delocalized.

Even though there is ample evidence suggesting the localization of collective
modes, the mechanism underlying the localization is not fully resolved. Does
it hinge on a simple Anderson-localization type interference mechanism [37] or
is it tied to the critical behavior of the underlying quantum phase transition
via mode coupling effects?

To shed light on this question, we study the collective excitations,
particularly the amplitude mode near the superfluid-insulator transition on a
random Voronoi-Delaunay (VD) lattice [38], where the disorder is exemplified
by the random connectivity of the lattice sites. (This type of disorder is
often called topological disorder.) Random VD lattices are commonly found in
nature in the form of amorphous materials, foams, animal skins, muscle fibers,
and plant tissues and are used in different applications ranging from material
science to astronomy [39, 40, 38]. It is now known that the Harris criterion
[41] gets modified on such topological lattices because the disorder displays
strong anti-correlations. It thus decays qualitatively faster with increasing
length scale than generic quenched disorder [42]. The relevancy of the
disorder operator is now set by a modified criterion which states that in the
case of topological disorder, the clean critical behavior is destabilized by
impurity fluctuations if the inequality $(d+1)\nu>2$ is violated. Here $\nu$
is the correlation length critical exponent and $d$ is the space
dimensionality of the system. It is well known, for instance, that the random
connectivity disorder is irrelevant at the critical point of the 3D Ising
model [43, 44], in accordance with this modified Harris criterion. In
contrast, it was shown in Ref. [45] that Anderson localization on a two-
dimensional random VD lattice is not affected qualitatively by the disorder
anticorrelations and thus behaves analogously to the case of generic
randomness, i.e., all states are localized. This implies that analyzing the
amplitude mode on a random VD lattice offers a route to disentangle the
possible localization mechanisms.

In this paper, we therefore study a system of interacting bosons on a random
VD lattice across the superfluid-insulator quantum phase transition. The
purpose is twofold. First, we want to determine the quantum critical behavior
to verify that this system is indeed governed by the modified Harris criterion
[42] rather than the regular one. Second, we want to analyze the collective
excitations in the superfluid phase and their localization properties. To do
so, we employ a quantum-to-classical mapping to rewrite the Bose-Hubbard model
(in the particle-hole symmetric limit) onto a classical (2+1) dimensional XY
model [46]. The mapped classical model is studied using large-scale Monte
Carlo simulations. The scalar susceptibility of the order parameter is
obtained in terms of imaginary frequency and then converted to real frequency
using the maximum entropy method (MaxEnt) [47]. In accordance with the
modified Harris criterion, we find that the critical behavior of the
topologically disordered system is identical to that of its clean,
translationally invariant counterpart. In contrast, the superfluid-Mott glass
quantum phase transition in a generically disordered system was shown to be in
a novel universality class [48]. We further show that the amplitude mode
remains delocalized across the transition of the topologically disordered
system, while it was shown to be spatially localized for generic disorder [34,
35, 36]. This suggests that the localization of the amplitude mode is caused
by mode-coupling and renormalization effects and tied to the critical behavior
of the transition rather than a simple Anderson-localization mechanism
(because, in the latter case, one would expect localization on a random VD
lattice).

The rest of this paper is organized as follows. In Section II, we discuss the
model Hamiltonian and its mapping to the corresponding classical system. In
Section III, we discuss the methods used in our study, including the
construction of the random VD lattice, the classical Monte Carlo simulations,
and the maximum entropy method. The critical behavior of the system is
described in detail in Section IV. The same section also addresses the
localization properties of the amplitude modes. A semi-analytical
inhomogeneous mean-field analysis in support of our study is described in
Section V. Section VI is devoted to the effects of generic randomness on top
of the topological disorder of the random VD lattice. We conclude our findings
in Section VII.

##  II Model

Our starting point is the Bose-Hubbard model on a two-dimensional random VD
lattice. This lattice consists of a set of lattice sites at independent random
positions, connected by bonds that are defined by means of the VD construction
[38]. Details of the lattice construction will be given in Sec. III.1. The
Bose-Hubbard Hamiltonian is given by:

$H=\frac{1}{2}\sum_{i}U_{i}(\hat{n}_{i}-\tilde{n})^{2}-\sum_{<ij>}J_{ij}(a_{i}^{\dagger}a_{j}+\textrm{h.c.})~{}.$
(1)

Here $a_{i}$, $a_{i}^{\dagger}$, $\hat{n}_{i}$ are boson annihilation,
creation, and number operators, respectively. Further, $U_{i}$ is the Hubbard
interaction at site $i$, and $\tilde{n}$ is the average filling which we fix
at a large integer value to ensure particle-hole symmetry. Finally, $J_{ij}$
denotes the hopping amplitude between the nearest neighbor sites of the random
VD lattice. For pure connectivity disorder, we set $J_{ij}\equiv J$ for all
nearest neighbors. Additional generic disorder can be introduced by site or
bond dilutions or by drawing $J_{ij}$ from a random distribution. The
translationally invariant (clean) version of the model (1), defined, e.g., on
a square lattice, is known to host two phases, namely the Mott insulator and
superfluid phases, separated by a quantum phase transition in the 3D XY
universality class [49].

In preparation for the Monte Carlo simulations, we rewrite the model (1) by
means of a quantum-to-classical mapping [46, 35, 34] onto a layered
(2+1)-dimensional classical XY model, see Fig. 1.

Figure 1: Layered (2+1)-dimensional random VD lattice. The same Voronoi layer
is repeated in the vertical (imaginary-time) direction.

Here, the third dimension represents imaginary time. The classical Hamiltonian
takes the form

$H_{\rm
cl}=-\sum_{\left<i,j\right>,\tau}J^{s}_{ij}\mathbf{S}_{i,\tau}\cdot\mathbf{S}_{j,\tau}-J^{t}\sum_{i,\tau}\mathbf{S}_{i,\tau}\cdot\mathbf{S}_{i,\tau+1}~{},$
(2)

where $\mathbf{S}_{j,\tau}$ is a classical $O(2)$ rotor, i.e., a two-component
unit vector, situated at spatial coordinate $i$ and imaginary time $\tau$. All
layers consist of the same two-dimensional VD tessellation. The spatial
interactions $J^{s}_{ij}$ couple Voronoi neighbors in the same layer (labeled
by the layer or imaginary time index $\tau$). $J^{t}$ elucidates the coupling
between equivalent sites in neighboring VD layers, i.e., in the imaginary time
direction. Due to the quantum-to-classical mapping, the disorder in the
classical Hamiltonian is perfectly correlated in the imaginary-time dimension.
This not only means that each layer contains the same VD tessellation, it also
implies that the interactions $J^{s}$ and $J^{t}$ do not depend on the
imaginary-time variable $\tau$.

The interaction constants appearing in the quantum and classical Hamiltonians
are related via $\beta_{c}J^{s}\sim J$ and $\beta_{c}J^{t}\sim 1/U$, where
$\beta_{c}=1/T$ is the inverse temperature of the classical model (with the
Boltzmann constant set to unity). This classical temperature does not
correspond to the actual temperature of the quantum system, which remains at
absolute zero. Thus, we can investigate the critical behavior of the
superfluid-insulator quantum phase transition occurring at zero temperature as
a function of the ratio $U/J$ by tuning the corresponding classical
temperature $T$ in the classical Hamiltonian (2). As we are interested in
universal properties, we set $J^{t}$ as well as all $J^{s}_{ij}$ between
Voronoi-Delaunay nearest neighbors to unity in the Monte Carlo simulations.

##  III Methods

###  III.1 Construction of the random Voronoi-Delaunay lattice

The random VD lattice consists of a set of lattice sites at random positions
uniformly distributed in a square box of linear size $L$ with periodic
boundary conditions. The point density is fixed at unity, i.e., a box of size
$L$ contains $L^{2}$ sites. The nearest neighbors of each site are determined
by means of the VD construction [38]. The Voronoi cell associated with a
lattice site is the area encompassing all points in the system area that are
closer to the given lattice site than to any other. Two lattice sites are
regarded as nearest neighbors if their Voronoi cells have a common edge. The
Delaunay triangulation is the dual of the Voronoi diagram and is defined as
the graph of all bonds linking neighboring pairs. The edges of the Voronoi
cells are the perpendicular bisectors of the edges in the Delaunay
triangulation, and the vertices of the Voronoi diagram are the circumcenters
of the triangles in the Delaunay triangulation.

In a random VD lattice, the number of nearest neighbors vary site to site.
These connectivity fluctuations introduce topological disorder into the
system. Barghathi and Vojta [42] showed that a topological constraint imposed
by the Euler equation in a two-dimensional triangulation graph introduces
strong anti-correlations between the connectivity fluctuations. This
suppresses the disorder effects and results in the modified Harris criterion,
$(d+1)\nu>2$, for the stability of a critical point against such topological
disorders. More details are discussed in Appendix A.

Our algorithm [45] for computing the Voronoi diagram and Delaunay
triangulation for a given set of random points follows a suggestion by
Tanemura et al. [50]. It relies on the striking “empty circumcircle property”
of a Delaunay triangulation which states that every triangular facet formed by
the nearest-neighbor bonds has an empty circumcircle, i.e., a circumcircle
that does not contain any other lattice sites. Our method finds the Voronoi
neighbors of a lattice site in two steps: (i) candidates for the neighbors are
identified based on their distance. (ii) Using these candidates, we then
construct all triangles with empty circumcircles for which the given site is
one of the vertices. The resulting algorithm is pretty efficient, finding the
Delaunay triangulation of $10^{6}$ sites takes about 30 seconds on an older
Intel core i5 CPU.

In order to study the effects of generic disorder on top of the topological
disorder of the random VD lattice, we also consider additional site dilution
which is explained in detail in Section VI.

###  III.2 Monte Carlo Simulations: Analysis of critical behavior

We use a combination of the Metropolis single-spin algorithm [51] and the
Wolff cluster algorithm [52] to perform our Monte Carlo simulations of the
classical Hamiltonian (2). A full Monte Carlo sweep is composed of a Wolff
cluster sweep and a Metropolis sweep over the entire lattice. The Wolff
cluster flips greatly reduce the critical slowing down, and the Metropolis
flips help to equilibrate isolated small clusters of lattice sites, which is
particularly useful in the presence of additional site dilution.

To reduce the computational cost in performing these Monte Carlo simulations,
we adopt the strategy outlined, e.g., in Refs. [53, 54, 35]. This approach is
based on averaging over a large number of disorder configurations while
keeping the number of measurement sweeps for each configuration relatively
small. This is possible because the Wolff algorithm leads to a short
equilibration time. Specifically, we perform simulations over about 1000 to
5000 disorder configurations (depending on system size and disorder strengths)
while keeping the number of equilibration sweeps to about 100 to 200 and the
number of measurement sweeps to 500. The strategy of simulating many disorder
configurations while keeping the number of measurements small reduces the
overall variance of various observables. As usual, the quality of the
equilibration is checked by comparing simulations with hot and cold starts.

Our primary observable is the order parameter

$\mathbf{m}=\frac{1}{N}\sum_{i,\tau}\mathbf{S}_{i,\tau},$ (3)

where $N$ is the total number of lattice sites. To identify the critical point
and to characterize the phase transition, we lay recourse to the Binder
cumulant of the order parameter, defined as

$U_{m}=\left[1-\frac{\langle|\mathbf{m}|^{4}\rangle}{3\langle|\mathbf{m}|^{2}\rangle^{2}}\right]_{dis}~{}.$
(4)

Here, $\langle\ldots\rangle$ denotes the thermodynamic (Monte Carlo) average
while the average over disorder realizations is represented by
$\left[\ldots\right]_{dis}$. Furthermore, we also compute the order parameter
susceptibility given by

$\chi=\frac{N}{T}\left[\langle|\mathbf{m}|^{2}\rangle-\langle|\mathbf{m}|\rangle^{2}\right]_{dis}~{},$
(5)

where $T$ is the classical temperature of the system.

Now, we discuss the finite-size scaling methods used to identify and analyze
the critical point. The introduction of quenched randomness breaks the
symmetry between the space and imaginary time directions. We thus cannot
expect the critical point to be Lorenz-invariant. Instead, the spatial
correlation length $\xi$ and the correlation length in imaginary time
direction $\xi_{\tau}$ must be treated as independent length scales.
Correspondingly, we need to treat the spatial system size $L$ and the size in
the imaginary time direction, $L_{\tau}$, as independent parameters. Note that
$L_{\tau}$ encodes the inverse physical temperature of the original quantum
model (1). Assuming conventional power-law dynamical scaling
$\xi_{\tau}\sim\xi^{z}$, the expected finite-size scaling form of the Binder
cumulant reads

$U_{m}=X_{U}(rL^{1/\nu},L_{\tau}/L^{z}),$ (6)

where $X_{U}$ is a dimensionless scaling function, $\nu$ is the correlation
length critical exponent, $z$ is the dynamical exponent, and
$r=(T-T_{c})/T_{c}$ is the reduced distance from criticality.

To find the critical temperature and to measure the value $z$ (which
determines the aspect ratio, i.e., the shape of the sample), we follow the
anisotropic scaling method described, e.g., in Refs. [55, 56]. We use the fact
that, at fixed $L$, the Binder cumulant as a function of $L_{\tau}$ has a
maximum at imaginary-time size $L_{\tau}^{max}$. The pair of system sizes
($L$, $L_{\tau}^{\rm max}$) gives the shape for which the correlations decay
about equally in the space and imaginary time directions; it is called the
optimal shape. According to the scaling form (6), the position of the maximum
at criticality, $r=0$, is expected to behave as

$L_{\tau}^{\rm max}\sim L^{z}.$ (7)

Moreover, because the Binder cumulant is a dimensionless quantity, its values
$U_{m}^{max}$ at criticality are independent of $L$. Hence, we can use this
method to estimate the transition temperature $T_{c}$.

The order parameter and order parameter susceptibility have analogous scaling
forms,

$m=L^{-\beta/\nu}X_{m}(rL^{1/\nu},L_{\tau}/L^{z})$ (8)

and

$\chi=L^{\gamma/\nu}X_{\chi}(rL^{1/\nu},L_{\tau}/L^{z})~{},$ (9)

where $\beta$ and $\gamma$ are the order parameter and susceptibility critical
exponents. $X_{m}$ and $X_{\chi}$ are dimensionless scaling functions. Once
the optimal shapes are found (fixing the second argument of the scaling
functions), the finite-size scaling analysis to extract the critical exponent
proceeds as usual [57].

In our simulations, we use system sizes up to $L=128$ in the space direction
and up to $L_{\tau}=400$ in the imaginary-time direction.

###  III.3 Amplitude mode and Wick rotation

The amplitude mode is an oscillation of the order parameter magnitude in the
symmetry-broken (ordered) phase. Since the local degrees of freedom
$\mathbf{S}_{i,\tau}$ in the classical Hamiltonian (2) have fixed magnitude,
we define a fluctuating local order parameter magnitude $\rho$ via coarse
graining. To this end, we average $\mathbf{S}_{i,\tau}$ over a small cluster
consisting of site $i$ and all its (spatial) nearest neighbors,

$\rho(\mathbf{x}_{i},\tau)=\frac{1}{\kappa_{i}}\left|\mathbf{S}_{i,\tau}+\sum_{j}\mathbf{S}_{j,\tau}\right|~{}.$
(10)

Here, $\kappa_{i}$ is the number of sites in the local cluster. The
fluctuations of this order parameter amplitude can be gleaned from the
imaginary-time scalar susceptibility

$\chi_{\rho\rho}(\mathbf{x},\tau)=\langle\rho(\mathbf{x},\tau)\rho(0,0)\rangle-\langle\rho(\mathbf{x},\tau)\rangle\langle\rho(0,0)\rangle$
(11)

which becomes translationally invariant after disorder averaging. A Fourier
transformation then gives the scalar susceptibility
$\tilde{\chi}_{\rho\rho}(\mathbf{q},i\omega_{m})$ as a function of wave vector
$\mathbf{q}$ and Matsubara frequency $i\omega_{m}$.

The Monte Carlo simulations yield the imaginary time (or Matsubara frequency)
scalar susceptibility. To make contact with the experimentally relevant real-
time or real-frequency scalar susceptibility, we need to perform an analytic
continuation (Wick rotation),

$\chi_{\rho\rho}(\mathbf{q},\omega)=\tilde{\chi}_{\rho\rho}(\mathbf{q},i\omega_{m}\rightarrow\omega+i0^{+})~{}.$
(12)

Here, $\omega$ is the real frequency. The (critical part of the) real
frequency scalar susceptibility is expected to satisfy the generalized scaling
relation [32, 34],

$\chi_{\rho\rho}(\mathbf{q},\omega)=\omega^{[(d+z)\nu-2]/(\nu
z)}X(\mathbf{q}r^{-\nu},\omega r^{-\nu z})~{}.$ (13)

Here, $X$ is the scaling function. The Wick rotation of the scalar
susceptibility can be implemented as an integral transformation

$\tilde{\chi}_{\rho\rho}(\mathbf{q},i\omega_{m})=\frac{1}{\pi}\int_{0}^{\infty}d\omega\chi_{\rho\rho}^{{}^{\prime\prime}}(\mathbf{q},\omega)\frac{2\omega}{\omega_{m}^{2}+\omega^{2}}$
(14)

where $\chi_{\rho\rho}^{{}^{\prime\prime}}(\mathbf{q},\omega)={\rm
Im}\chi_{\rho\rho}(\mathbf{q},\omega)$ is the spectral function.
Unfortunately, this analytic continuation is an ill-posed problem and is
sensitive to errors from the numerical data. However, we can overcome this
problem by using the maximum entropy method (MaxEnt) [47]. More details of the
technique are given in Appendix B.

Note that we will focus on wave vector $\mathbf{q}=0$ in the following.
Nonzero wave vectors lead to extra numerical complications for our random VD
lattice because the coordinates of the lattice sites can take arbitrary real
values that vary from disorder realization to disorder realization. Thus, the
efficient fast-Fourier transformation algorithms cannot be applied directly.

##  IV Results

In this section, we present the results for the quantum critical behavior of
the superfluid-insulator transition in the Bose-Hubbard model (1) on a random
VD lattice. Furthermore, we study the localization properties of the amplitude
mode in the symmetry-broken (superfluid) phase in the presence of topological
disorder.

###  IV.1 Analysis of the critical behavior

As discussed in Sec. III.2, the first step towards analyzing the critical
behavior consists in finding the critical point by using the anisotropic
finite-size scaling of the Binder cumulant $U_{m}$. To do so, we compute
$U_{m}$ for a given spatial system size $L$ over a range of imaginary-time
system sizes $L_{\tau}$ at different temperatures around the phase transition.
We identify the critical temperature of the system as the temperature at which
the peak values of the $U_{m}$ vs. $L_{\tau}$ curves are independent of $L$
(see Fig. 2).

Figure 2: Binder cumulant $U_{m}$ as a function of the imaginary-time system
size $L_{\tau}$ for different spatial system sizes $L$ at $T=3.15480<T_{c}$,
$T=3.15505=T{c}$ and $T=3.15530>T_{c}$. The peak values are independent of $L$
at $T_{c}$ whereas they decrease with $L$ for $T>T_{c}$ and increase with $L$
for $T<T_{c}$. The statistical errors are of the order of the symbol size or
smaller.

By this method, the transition temperature of the classical Hamiltonian (2) is
found to be $T_{c}=3.15505(25)$.

To obtain the optimal shapes, we determine the position $L_{\tau}^{\rm max}$
for each $U_{m}$ vs. $L_{\tau}$ curve (by fitting a quadratic parabola in
terms of $\ln{L_{\tau}}$). The dynamical critical exponent $z$ is obtained by
fitting the optimal shapes to the relation $L_{\tau}^{\rm max}=aL^{z}$, as
shown in Fig. 3.

Figure 3: $L_{\tau}^{\rm max}$ (red circles) and $|\mathbf{m}|$ (blue
triangles) as functions of $L$ at $T_{c}$. The lines are fits to the
predictions of the scaling relations $L_{\tau}^{\rm max}=aL^{z}$ and
$m=aL^{-{\beta}/{\nu}}$. The statistical errors are of the order of the symbol
size or smaller.

Other critical exponents can be measured by analyzing systems of the optimal
shapes (fixing the argument $L_{\tau}/L^{z}$ in the scaling functions) at or
close to $T_{c}$. The order parameter critical exponent $\beta/\nu$ is
obtained from the decay of the order parameter with $L$ right at criticality,
$m=aL^{-{\beta}/{\nu}}$. Similarly, the critical exponent $\gamma/\nu$ of the
order parameter susceptibility is obtained from $\chi=aL^{{\gamma}/{\nu}}$ at
criticality, see Fig. 4.

Figure 4: $\chi$ (green circles) and $x_{L}$ (magenta triangles) as functions
of $L$ at $T_{c}$. The lines are fits to the power-law relations
$\chi=aL^{{\gamma}/{\nu}}$ and $x_{L}=aL^{{1}/{\nu}}$. The statistical errors
are of the order of the symbol size or smaller.

To get the correlation length critical exponent $\nu$, we first define the
observable, $x_{L}=|(d/dT)\ln{|\textbf{m}|}|$. The corresponding scaling form
can be obtained by differentiating the scaling form (8) of the order parameter
with respect to $T$. The resulting $L$ dependence at criticality is given by
$x_{L}=aL^{{1}/{\nu}}$.

The values of critical exponents resulting from this analysis are tabulated in
Table 1 along with the exponents obtained for the clean three-dimensional XY
universality class [58] and those for the superfluid-Mott glass transition
with generic disorder [48].

Exponent | Clean [58] | Generic disorder [48] | VD lattice
---|---|---|---
$\nu$ | 0.6717 | 1.16(5) | 0.672(8)
${\beta}/{\nu}$ | 0.519 | 0.48(2) | 0.520(4)
${\gamma}/{\nu}$ | 1.962 | 2.52(4) | 1.950(10)
$z$ | 1 | 1.52(3) | 1.008(9)
Table 1: Critical exponents of the superfluid-insulator quantum phase
transition on the random VD lattice, compared to the clean case [58] and the
case of generic disorder [48]. The quoted error bars include the error due to
the uncertainty of $T_{c}$ as well as the robustness of the fits against
removing data points at the ends of the $L$ range.

The table clearly shows that the critical exponents of the superfluid
transition on the random VD lattice agree with those of the equivalent
translationally invariant (clean) system rather than with the exponents of
superfluid transition in the presence of generic disorder. This is consistent
with the modified Harris criterion derived by Barghathi and Vojta [42], which
states that a clean critical point is stable against the topological disorder
of a random VD lattice if $(d+1)\nu>2$. The effects of additional generic
uncorrelated disorder on top of the topological disorder will be discussed in
Section VI.

###  IV.2 Analysis of the amplitude mode

We now discuss the behavior of the amplitude (Higgs) mode close to the quantum
phase transition. It is known that, in the absence of disorder, the system
features a well-defined, soft-gapped amplitude mode characterized by a peak in
the spectral function
$\chi_{\rho\rho}^{{}^{\prime\prime}}(\mathbf{q}=0,\omega)$ at the Higgs energy
$\omega_{H}$ [33, 35, 34]. This peak survives as the quantum phase transition
is approached from the superfluid side. In agreement with the scaling form
(13), the Higgs energy vanishes as $\omega_{H}\sim r^{z\nu}$. In contrast, a
Higgs peak is not observable in the presence of generic disorder where the
scalar response is broad and non-critical [34, 35, 36], indicating that the
amplitude mode is localized.

In Fig. 5, we present the scalar spectral function
$\chi_{\rho\rho}^{{}^{\prime\prime}}(\mathbf{q}=0,\omega)$ for the Hamiltonian
(2) defined on a random VD lattice at different distances from criticality in
the superfluid phase.

Figure 5: Spectral function
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ as a function of
real frequency $\omega$ for different distances to the (classical) critical
temperature. The simulation is performed for a lattice with $L=L_{\tau}=128$
(reflecting the dynamical exponent $z=1$).

The data clearly demonstrate that the amplitude mode in this topologically
disordered system behaves analogously to the clean case, as
$\chi_{\rho\rho}^{{}^{\prime\prime}}(\mathbf{q}=0,\omega)$ features a well-
pronounced low-energy peak that softens as the critical point is approached.
No indications of spatial localization are observed. Additionally, the Higgs
peak in the scalar spectral function shows an excellent scaling collapse
according to Eq. (13) using the clean critical exponents, as seen in Fig. 6.

Figure 6: Scaling plot of the spectral density
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ on the superfluid
side of the quantum phase transition for the Voronoi lattice, as suggested by
the scaling form (13). The exponent values, $\nu=0.6717$ and $z=1$ with $d=2$,
belong to the clean XY universality class.

The results for this topologically disordered system are thus qualitatively
different from those obtained on a square lattice with site dilution [35, 34]
which represents the generic disorder case.

Our findings give us important information on the physical origins of the
localization mechanism of the amplitude mode. On the one hand, a tight-binding
model of non-interacting particles on a random VD lattice has been shown to
feature Anderson localization [45]. In other words, the connectivity disorder
anticorrelations of the random VD lattice do not qualitatively affect the
Anderson localization mechanism. On the other hand, we found here that the
amplitude mode close to the superfluid-insulator transition on a random VD
lattice does not become localized, in contrast to the generic disorder case.
This indicates that the localization of the amplitude mode observed in Refs.
[35, 34] may not be governed by an Anderson-type localization mechanism. The
modified Harris criterion [42] states that the topological disorder of the
random VD lattice is an irrelevant operator at the critical point of a (2+1)
dimensional XY model. Our results are therefore consistent with the
localization of the amplitude mode (close to the phase transition) being tied
to the critical behavior.

##  V Inhomogeneous mean-field theory

To gain further insight into the behavior of the collective excitations of the
bosonic Hubbard model (1) on a random VD lattice, we apply an inhomogeneous
mean-field theory [34, 36]. The excitations are obtained from expanding the
Hamiltonian about the spatially inhomogeneous mean-field solution. At the
Gaussian level, this results in noninteracting bosonic excitations. This
theory generalizes the approach of Refs. [59, 60] to the disordered case and
is related to the bond-operator method for Heisenberg magnets [61]. In the
presence of disorder, the approach captures Anderson localization physics but
not the mode-mode coupling and renormalizations near criticality, allowing us
to disentangle the mechanisms for the amplitude mode localization or lack
thereof.

To derive the mean-field theory, we truncate the local Hilbert space at
lattice site $j$ to three basis states, $|-_{j}\rangle,|0_{j}\rangle$, and
$|+_{j}\rangle$, representing the particle numbers $n_{j}=\bar{n}-1,\bar{n}$,
and $\bar{n}+1$, respectively. The variational ground state wave function is
written in product form, $|\Phi_{0}\rangle=\prod_{j}|\phi_{0j}\rangle$ with

$\displaystyle|\phi_{0j}\rangle$ $\displaystyle=$
$\displaystyle\cos(\theta_{j}/2)|0_{j}\rangle$ (15)
$\displaystyle+\sin(\theta_{j}/2)\left(e^{i\eta_{j}}|+_{j}\rangle+e^{-i\eta_{j}}|-_{j}\rangle\right)/\sqrt{2}~{}.$

The parameters (mixing angles) $\theta_{j}$ control the character of the
ground state, yielding a Mott insulator for $\theta_{j}=0$ and a superfluid
for $\theta_{j}>0$. The parameter $\eta_{j}$ represents the phase of the local
superfluid order parameter $\langle
a_{j}^{\dagger}\rangle\propto\psi_{j}=\sin\theta_{j}e^{-i\eta_{j}}$.

Minimizing the ground state energy $E_{0}=\langle\Phi_{0}|H|\Phi_{0}\rangle$
w.r.t. to the mixing angles $\theta_{j}$ and phases $\eta_{j}$ gives uniform
$\eta_{j}=\eta=\textrm{const}$ (which we set to zero in the following) while
the mixing angles fulfill the mean-field equations

$U_{i}\sin\theta_{i}=4\bar{n}\cos\theta_{i}\sum_{j}J_{ij}\sin\theta_{j}~{}.$
(16)

In the presence of randomness, this large set of coupled nonlinear equations
needs to be solved numerically [34, 36].

Figure 7 presents the resulting order parameter $\psi$ of the Bose-Hubbard
model (1) on a random VD lattice as a function of $U/J$.

Figure 7: Superfluid order parameter $\psi$ of the Bose-Hubbard model (1) on a
random VD lattice vs. $U/(\bar{n}J)$, computed from the mean-field equations
(16). Both the average local order parameter and the typical value (geometric
average) are shown. For comparison, the corresponding curve for a clean,
triangular lattice is also shown. The data are averages over 1000 disorder
configurations for a system size for $L=128$. The resulting statistical errors
are below the line thickness.

The data shows that the superfluid-insulator transition happens at
$U/(\bar{n}J)\approx 25.3$. This is close to the value for the regular (clean)
triangular lattice, $U/(\bar{n}J)=24$ which is included in the figure for
comparison. (This makes sense because the average coordination number of a
random VD lattice is exactly 6, identical to the coordination number of the
triangular lattice.) The figure also shows that the average order parameter
$(1/N)\sum_{j}\psi_{j}$ and the typical value $\exp[(1/N)\sum_{j}\ln\psi_{j}]$
are almost indistinguishable. This indicates that the order parameter
fluctuations are weak, much smaller than those of a site-diluted Bose-Hubbard
model at moderate dilutions (see Fig. 2 of Ref. [36]).

In order to investigate excitations of the mean-field ground state, we now
transform the basis in the Hilbert space associated with site $j$ from
$|-_{j}\rangle,|0_{j}\rangle$, $|+_{j}\rangle$ to a new orthonormal basis
consisting of $|\phi_{0j}\rangle,|\phi_{Hj}\rangle,|\phi_{Gj}\rangle$. The
states

$\displaystyle|\phi_{Hj}\rangle$ $\displaystyle=$
$\displaystyle\sin(\theta_{j}/2)|0_{j}\rangle-\cos(\theta_{j}/2)\left(|+_{j}\rangle+|-_{j}\rangle\right)/\sqrt{2}\quad$
$\displaystyle|\phi_{Gj}\rangle$ $\displaystyle=$ $\displaystyle
i\left(|+_{j}\rangle-|-_{j}\rangle\right)/\sqrt{2}$ (17)

correspond to fluctuations of the order parameter amplitude and phase,
respectively, w.r.t. to $|\phi_{0j}\rangle$. We now expand the Hamiltonian (1)
to quadratic order in the boson operators $b^{\dagger}_{Hj}$ and
$b^{\dagger}_{Gj}$ that create these states out of the fictitious vacuum
(ground state). The resulting fluctuation Hamiltonian decouples into amplitude
and phase parts, $H_{MF}=E_{0}+H_{H}+H_{G}$, which both take the form

$\displaystyle H_{\alpha}=\sum_{i}A_{\alpha i}b^{\dagger}_{\alpha i}b_{\alpha
i}+\sum_{\langle ij\rangle}B_{\alpha ij}(b^{\dagger}_{\alpha i}+b_{\alpha
i})(b^{\dagger}_{\alpha j}+b_{\alpha j})~{},~{}$ (18)

($\alpha=H,G$). The coefficients $A_{\alpha i}$ and $B_{\alpha ij}$ are
functions of the mixing angles $\theta_{j}$. $H_{H}$ and $H_{G}$ can be
diagonalized independently by bosonic Bogoliubov transformations, $b_{\alpha
j}=\sum_{k}(u_{\alpha jk}d_{\alpha k}+v_{\alpha jk}^{\ast}d_{\alpha
k}^{\dagger})$. These Bogoliubov transformations have to be performed
numerically because of the randomness.

The localization properties of the eigenstates of $H_{H}$ and $H_{G}$ can be
characterized by their inverse participation numbers. Following Ref. [61], the
inverse participation number of eigenstate $k$ is given by
$P^{-1}(k)=\sum_{j}(|u_{\alpha j0}|^{2}-|v_{\alpha j0}|^{2})^{2}$, and the
corresponding generalized dimension reads $\tau_{2}(k)=\ln P(k)/\ln L$ (for
details see Ref. [36]). Figure 8 shows the dimension $\tau_{2}(0)$ of the
lowest-energy Goldstone and Higgs excitations of the Bose-Hubbard model (1) on
a random VD lattice as a function of $U/J$.

Figure 8: Generalized dimension $\tau_{2}(0)$ of the lowest-energy Goldstone
(solid lines) and Higgs modes (dashed lines), in Gaussian approximation, vs.
interaction $U/(\bar{n}J)$ for different system sizes $L$. The data are
averages over 1000 disorder configurations for a system size for $L=128$.

In the insulating phase, $U/(\bar{n}J)\gtrsim 25.3$, both excitations are
degenerate and strongly localized because $\tau_{2}$ rapidly decreases towards
zero with increasing system size. In the superfluid phase,
$U/(\bar{n}J)\lesssim 25.3$, in contrast, the two excitations behave
differently. The lowest Higgs mode remains strongly localized, whereas the
lowest Goldstone mode rapidly delocalizes; its $\tau_{2}$ increases with
increasing $L$ and approaches the embedding dimension $d=2$. The
delocalization of the Goldstone mode agrees with a general symmetry analysis
[62, 63] and with explicit results for Goldstone modes in a number of systems.

Interestingly, the behavior of the amplitude mode for the Bose-Hubbard model
on a random VD lattice agrees with that of the site-diluted Bose-Hubbard model
[34, 36]. In both cases, the mean-field theory predicts that the lowest Higgs
excitation is spatially localized. However, for the random VD lattice, this
does not agree with the results of the Monte Carlo simulations reported
earlier in Sec. IV.2 where it was found that the amplitude mode on the VD
lattice behaves just as in the clean case. As the inhomogeneous mean-field
approach captures the Anderson localization physics but not the mode-mode
coupling and renormalization of the full many-particle problem, the
discrepancy between the Monte Carlo results and the mean-field results for the
amplitude mode on the VD lattice provides further evidence for the amplitude
mode localization or lack thereof not being driven by an Anderson localization
mechanism.

##  VI Adding additional uncorrelated disorder

###  VI.1 Bose-Hubbard model on a site-diluted random VD lattice

The Monte Carlo results reported in Sec. IV demonstrated that the amplitude
mode remains delocalized in the presence of the topological disorder of a
random VD lattice. The results also suggested that this stems from the fact
that the topological disorder is an irrelevant perturbation at the clean
superfluid-insulator critical point, and the transition belongs to the same
universality class as the translationally invariant model. To provide further
evidence for this hypothesis, we now add generic uncorrelated disorder in the
form of site dilution to the topological disorder of the random VD lattice.
This uncorrelated disorder is a relevant perturbation at the clean superfluid-
insulator critical point. If the hypothesis is correct, we thus expect not
only a crossover to the disordered critical behavior of Ref. [48], we also
expect the amplitude mode to become localized.

After introducing random vacancies into the Bose-Hubbard Hamiltonian (1) on a
random VD lattice, we carry out the quantum-to-classical mapping and arrive at
a classical XY model on a layered VD lattice described by the Hamiltonian

$H_{\rm
dis}=-\sum_{\left<i,j\right>,\tau}J^{s}_{ij}\epsilon_{i}\epsilon_{j}\mathbf{S}_{i,\tau}\cdot\mathbf{S}_{j,\tau}-J^{\tau}\sum_{i,\tau}\epsilon_{i}\mathbf{S}_{i,\tau}\cdot\mathbf{S}_{i,\tau+1}~{}.$
(19)

Here, the independent quenched random variables $\epsilon_{i}$ take values 0
and 1 with probabilities $p$ and $1-p$, respectively. Since the positions of
the vacancies do not depend on the imaginary-time coordinate $\tau$, the
disorder is perfectly correlated in the imaginary-time direction (columnar
disorder perpendicular to the VD layers).

###  VI.2 Critical Behavior

We perform Monte Carlo simulations of the site-diluted XY model (19) for six
different vacancy concentrations, namely $p=0.125$, $0.200$, $0.270$ $0.300$,
$0.400$ and $0.415$. We analyze the data employing anisotropic finite-size
scaling as outlined in Sec. III.2. The resulting critical temperatures are
listed in Tab. 2.

Site dilution $p$ | Critical temperature $T_{c}$
---|---
0.000 | 3.15505(25)
0.125 | 2.850(2)
0.200 | 2.648(2)
0.270 | 2.448(2)
0.300 | 2.356(4)
0.400 | 2.020(4)
0.415 | 1.960(4)
Table 2: Transition temperatures of the XY model (19) with correlated site
dilution.

As explained in Sec. III.2, the dynamical exponent $z$ can be gleaned from the
dependence on the spatial system size $L$ of the peak position $L_{\tau}^{\rm
max}$ of the Binder cumulant curves right at the critical temperature. The
corresponding data for all dilution values are presented in Fig. 9.

Figure 9: Double logarithmic plot of $L_{\tau}^{\rm max}$ vs. $L$ at $T_{c}$
for various vacancy concentrations $p$. The solid lines represent a combined
fit to $L_{\tau}^{\rm max}=aL^{z}(1+bL^{-\omega})$ with universal $z$ and
$\omega$ but dilution-dependent coefficients $a$ and $b$. The fit gives
$z=1.50(3)$ and $\omega=1.2(2)$.

In contrast to the corresponding data for the undiluted random VD lattice
(Fig. 3), the data in Fig. 9 display significant deviations from pure power-
law behavior $L_{\tau}^{\rm max}\sim L^{z}$, in particular for lower dilutions
$p$. This deviation can be attributed to a cross-over from the clean critical
behavior to the disordered critical behavior in response to the (uncorrelated)
random-mass disorder generated by the vacancies. To account for this
crossover, we include a correction-to-scaling term via $L_{\tau}^{\rm
max}=aL^{z}(1+bL^{-\omega})$. The correction term is governed by the
irrelevant exponent $\omega$.

The exponents $z$ and $\omega$ are assumed to be universal (independent of the
dilution $p$) whereas the coefficients $a$ and $b$ depend on $p$.
Consequently, we perform a combined fit that includes the data for all six
dilutions, with universal $z$ and $\omega$ but nonuniversal $a$ and $b$; it
yields $z=1.50(3)$ and $\omega=1.2(2)$. The fit is of good quality, giving a
reduced $\tilde{\chi}^{2}\approx 0.53$. The leading corrections to scaling
seem to vanish around $p=0.400$, where the coefficient $b$ of the correction
term changes sign. This implies that for $p\approx 0.400$, a simple power-law
fit should suffice. Performing power-law fits for $p=0.400$ and $p=0.415$
yields $z=1.495$ and $z=1.529$, respectively. These values agree with the
estimate provided by the combined fit within the statistical errors.

Comparing our results for the site-diluted random VD lattice to those for the
site-diluted square lattice [48] (see also Tab. 1), we notice that values of
the dynamical exponent $z$ agree very well, whereas the values of the
irrelevant exponent $\omega$ differ from each other. This dichotomy can be
explained by the fact that the value of the exponent $z$ reflects the
asymptotic critical behavior governed by the finite-disorder renormalization
group fixed point. In contrast, the value of $\omega$ is related to the
leading irrelevant operator at this fixed point which can be different for
different models. Specifically, the Bose-Hubbard model on the diluted random
VD lattice contains additional (anti-correlated) topological disorder which
may affect the crossover scaling and thus the value of $\omega$.

Having found the optimal sample shapes, we can now determine the exponents
$\beta/\nu$, $\gamma/\nu$, and $\nu$ by analyzing the system size dependence
of the order parameter, the order parameter susceptibility, and
$x_{L}=|(d/dT)\ln{|\textbf{m}|}|$ at the critical temperature, just as in Sec.
IV.1. To capture deviations from pure power-law behavior, we include
corrections to scaling terms, and fit the data with the functional forms
$m=aL^{-\beta/\nu}(1+bL^{-\omega})$, $\chi=aL^{\gamma/\nu}(1+bL^{-\omega})$,
and $x_{L}=aL^{1/\nu}(1+bL^{-\omega})$. For each quantity, we perform a
combined fit of the curves for all dilutions with universal exponent values
but dilution-dependent coefficients $a$ and $b$. (In some cases, we need to
exclude smaller system sizes to improve the quality of fit and achieve an
acceptable $\bar{\chi}^{2}$.) As an example, the analysis of $x_{L}$ leading
to a value for the exponent $\nu$ is shown in Fig. 10.

Figure 10: Double logarithmic plot of $x_{L}=|d\ln{|\textbf{m}|}/dT|$ as
function of $L$ near $T_{c}$ for various dilutions $p$. The solid lines
represent a combined fit of the data for all dilutions with
$x_{L}=aL^{1/\nu}(1+bL^{-\omega})$ with universal $\nu$ and $\omega$ but
nonuniversal coefficients $a$ and $b$. It gives $\nu=1.08(2)$ and
$\omega=0.52(5)$. The dotted lines mark smaller system sizes that are not
included in the fit.

These fits give the exponent values $\beta/\nu=0.50(3)$, $\gamma/\nu=2.47(4)$,
and $\nu=1.08(2)$. They all agree within their error bars with the exponents
reported for the site-diluted square lattice system [48] (see also Tab. 1).
This confirms the hypothesis that adding uncorrelated randomness on top of the
topological disorder destabilizes the clean critical behavior and causes the
system to flow to the generic finite-disorder fixed point of Ref. [48]. As in
the case of $L_{\tau}^{\rm max}$, the values of the irrelevant exponent
$\omega$ emerging from the combined fits of $m$, $\chi$, and $x_{L}$ do not
agree well with their counterparts in the diluted square-lattice system. This
can again be explained by the existence of additional irrelevant operators
111We also emphasize that extracting subleading corrections from fits is more
difficult numerically than determining the leading terms, and additional
uncertainties arise if more than one correction term contributes..

###  VI.3 Amplitude Mode

We now turn our attention to the behavior of the amplitude mode under the
combined influence of the topological disorder of the random VD lattice and
the uncorrelated disorder due to dilution. To do so, we measure the imaginary
time scalar susceptibility (11) within our Monte Carlo simulations and
analytically continue it to real time and frequency by means of the maximum
entropy method applied to the transformation (14), in analogy to Sec. IV.2.
The resulting scalar spectral functions
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ in the superfluid
phase close to criticality are presented in Fig. 11 for two different
dilutions, $p=0.125$ and $p=0.400$.

Figure 11: Spectral density
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ in the superfluid
phase as a function of real frequency $\omega$ for dilutions (a) $p=0.125$ and
(b) $p=0.400$. The simulations were performed at the optimal shapes (a)
$L=120$ and $L_{\tau}=256$ (b) $L=88$ and $L_{\tau}=256$.

For the weaker dilution of $p=0.125$ (Fig. 11(a)), a peak in the spectral
density is still visible for the larger distances to criticality, but its
magnitude is already suppressed compared to the case of pure topological
disorder shown in Fig. 5. With decreasing distance from criticality, this peak
rapidly shrinks and broadens. Moreover, the peak position does _not_ soften as
expected from Eq. (13) as the critical point is approached. Thus, the scalar
response of the diluted system violates naive scaling. This becomes even more
obvious for the higher dilution of $p=0.400$ (Fig. 11(b)) for which the scalar
spectral function does not show any low-energy peak. Instead,
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ is dominated by a
broad hump centered at a microscopic energy scale. This hump is non-critical,
i.e., it is completely insensitive to the distance from criticality.

In Fig. 12, we compare the spectral functions for several dilutions at a fixed
distance from criticality, $r=-0.02$, along with the non-diluted lattice,
i.e., the $p=0$ case.

Figure 12: Spectral density
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\mathbf{q}=0,\omega)$ as a function of
real frequency at a fixed distance from criticality $r=|T-T_{c}|/T_{c}=-0.02$,
for the dilutions $p=0$, $0.270$, and $0.415$, for $L=128$ and the optimal
$L_{\tau}$ for each $p$.

For the case of the lattice with only topological disorder, we see (as in Sec.
IV.2) a sharp Higgs peak. Introducing site dilution rapidly suppresses this
peak and replaces it by a broad non-critical hump.

The broadening and suppression of the Higgs peak demonstrated in Figs. 11 and
12 is completely analogous to the behavior observed at the superfluid-Mott
glass transition on a diluted square lattice [34, 35]. To explain this
behavior, we note that every observable close to a critical point has a
singular part and a non-singular part. The singular part is governed by the
long-wavelength low-energy critical fluctuations and generically follows
scaling whereas the non-singular part is due to non-critical microscopic
degrees of freedom. The scalar susceptibility can thus be decomposed as

$\chi_{\rho\rho}(\mathbf{q},\omega)=\chi_{\rho\rho}^{\rm
reg}(\mathbf{q},\omega)+\omega^{[(d+z)\nu-2]/(\nu
z)}X(\mathbf{q}r^{-\nu},\omega r^{-\nu z})~{}$ (20)

where $\chi_{\rho\rho}^{\rm reg}(\mathbf{q},\omega)$ represents the non-
singular part, and the second term is the scaling form of the singular part.
At the clean critical point (which also controls the transition in the system
with a purely topological disorder, see Sec. IV.1), $z=1$ and $\nu=0.672$.
Consequently, the exponent $[(d+z)\nu-2]/\nu z$ which governs the magnitude of
the singular part of $\chi_{\rho\rho}$ takes a value very close to zero,
meaning that this magnitude does not change as the critical point is
approached. On the other hand, at the finite-disorder fixed point controlling
the transitions in the diluted systems, $z=1.52$ and $\nu=1.16$. The exponent
$[(d+z)\nu-2]/\nu z$ therefore takes a strongly positive value of about 1.18.
This implies that the singular part is rapidly suppressed as the distance from
criticality decreases, and the scalar response becomes dominated by the
noncritical microscopic excitations contained in the non-singular part. For
the weakest dilution we studied here, $p=0.125$, the suppression of the Higgs
peak with decreasing distance from criticality can be seen directly in Fig.
11(a). For the higher dilution of $p=0.400$ shown in Fig. 11(b), the
suppression is much stronger, and only the non-singular part of
$\chi_{\rho\rho}$ can be seen for all temperatures.

##  VII Conclusion

In summary, we investigated the quantum critical behavior and the collective
excitations near the two-dimensional superfluid-insulator quantum phase
transition of interacting bosons on a random VD lattice. This was accomplished
by mapping the Bose-Hubbard model onto a classical XY Hamiltonian on a
(2+1)-dimensional layered VD lattice. This classical model was then simulated
by means of Monte Carlo simulations to identify the universality class of the
superfluid-insulator transition. We also studied the behavior of the scalar
spectral function to analyze the localization properties of the amplitude mode
in the presence of topological disorder. To disentangle possible mechanisms
for the amplitude mode localization, the Monte Carlo simulations were
complemented by an inhomogeneous mean-field theory.

Our results can be summarized as follows. In consonance with the modified
Harris criterion [42] (but in disagreement with the regular one [41]), we
found that the superfluid-insulator transition in the presence of topological
disorder imposed by the random VD lattice belongs to the clean 3D XY
universality class. On the superfluid side of the quantum phase transition,
the scalar spectral function behaves in the same way as in the clean
(translationally invariant) case. Specifically, the spectral function shows a
sharp peak at a Higgs energy $\omega_{H}$, which softens as the transition is
approached and conforms to the naive scaling expectations [32, 34]. This
implies that the amplitude mode is not localized by the topological disorder.
As the random VD lattice is known to support Anderson localization for a
system of non-interacting particles [45], our results indicate that amplitude
mode localization observed in the presence of generic (uncorrelated) disorder
[34, 35] does not stem from an Anderson-localization-type mechanism. Instead,
it appears to be tied to the critical behavior of the quantum phase transition
via the scale dimension of the scalar susceptibility.

This observation also agrees with our results for the Bose-Hubbard model on a
random VD lattice with additional site dilution. The uncorrelated disorder
introduced by the vacancies destabilizes the clean critical behavior and
causes the transition to belong to the disordered universality class
identified in Ref. [48]. At the same time, the sharp Higgs peak in the scalar
spectral function is destroyed, and the scalar response becomes broad and non-
critical.

Recently, Beattie-Hauser and Vojta [65] studied the scalar susceptibility in a
site-diluted three-dimensional classical XY model. The dilution was
uncorrelated in all three directions, which means this model cannot be
interpreted as arising from the quantum-to-classical mapping of a diluted
quantum Hamiltonian (which would lead to disorder correlated in the imaginary
time direction). Nonetheless, its scalar susceptibility can be compared with
the scalar susceptibility of the mapped Hamiltonian (2) in imaginary time
(i.e., before the Wick rotation). The three-dimensional classical model with
uncorrelated vacancies fulfills the conventional Harris criterion $d\nu>2$, if
just barely. Its critical behavior, therefore belongs to the clean 3D XY
universality class. Beattie-Hauser and Vojta [65] determined that the scalar
susceptibility of this model in the long-range ordered phase fulfills naive
scaling. This provides additional evidence that the presence or absence of
unconventional behavior of the scalar susceptibility in a disordered system is
controlled by its scale dimension at the transition in question.

Our results have broader implications for disordered quantum phase
transitions. Based on the scaling form (20), one can make a general prediction
about the fate of the scalar susceptibility at any quantum phase transition
controlled by a finite-disorder fixed point. At such a fixed point, the
correlation length exponent $\nu$ must fulfill the inequality $d\nu>2$ [66].
Consequently, the exponent $[(d+z)\nu-2]/\nu z$ which controls the amplitude
of the scalar susceptibility will be positive and larger than unity. This
implies that, at any finite-disorder critical point, the singular part of
scalar susceptibility is suppressed as the transition is approached, and the
scalar response is dominated by non-critical microscopic fluctuations.

In many experimental systems, the (bare) disorder is expected to be only
moderately strong. To interpret experimental data, it will therefore be
important to understand the crossover from clean to disordered behavior as the
transition is approached. This crossover is explicitly visible in Fig. 11(a),
but a quantitative analysis requires the study of larger, more weakly
disordered systems close to criticality. This is numerically very expensive
and thus remains a task for the future.

###### Acknowledgements.

The simulations were performed on the Pegasus and Foundry clusters at Missouri
S&T and the Aqua cluster at IIT Madras. T.V. acknowledges support from the
National Science Foundation under Grants No. DMR-1828489 and OAC-1919789 as
well as Grants No. PHY-1748958 and PHY-2309135 to the Kavli Institute for
Theoretical Physics. P.K.V. acknowledges support via the IIE Travel award by
the Global Engagement Office, IIT Madras. R.N. and P.K.V. also acknowledge
funding from the Center for Quantum Information Theory in Matter and
Spacetime, IIT Madras and from the Department of Science and Technology, Govt.
of India, under Grant No. DST/ICPS/QuST/Theme-3/2019/Q69, as well as support
from the Mphasis F1 Foundation via the Centre for Quantum Information,
Communication, and Computing (CQuICC).

## Appendix

##  Appendix A Connectivity disorder fluctuations in the random VD lattice

The quenched disorder in the Bose-Hubbard model (1) stems from the random
connectivity of the underlying random VD lattice. At the mean-field level, the
local coupling strength is proportional to the number of nearest neighbors
(coordination number) that a given lattice site has. The spatial fluctuations
of the coordination number of a random VD lattice were studied in detail in
Ref. [42]. To this end, a large random VD lattice was divided into blocks of
linear size $L_{b}$. (This can be done either via the real space positions of
the lattice sites or via the link distance, the number of bonds separating two
sites. Both partitions lead to the same results.) A block-averaged
coordination number for block no. $k$ can be defined by

$Q_{k}=\frac{1}{N_{b,k}}\sum_{i\in k}q_{i},$ (21)

where $q_{i}$ is the coordination number of site $i$, and $N_{b,k}$ is the
number of sites in block $k$.

The spatial disorder fluctuations can be characterized by the variance of the
block averaged coordination number $\sigma^{2}_{Q}$ and its dependence on the
block size $L_{b}$. As the number of lattice sites in a block behaves as
$L_{b}^{d}$, one might naively expect the central limit theorem result
$\sigma^{2}_{Q}\sim L^{-d}$. However, the analysis in Ref. [42] revealed a
more rapid decrease, $\sigma^{2}_{Q}\sim L^{-(d+1)}$. This rapid decay of the
disorder fluctuations under coarse graining stems from a topological
constraint imposed by the Euler equation of a two-dimensional graph,
$N-E+F=\chi$. Here, $N$ is the number of lattice sites (vertices), $E$ is the
number of edges and $F$ is the number of facets in the graph. $\chi$ is the
Euler characteristic which is equal to $0$ for periodic boundary conditions
(torus topology). If every facet is a triangle, as is the case for the VD
triangulation, $3F=2E$ because each triangle has three edges, and each edge is
shared by two triangles. Consequently, $E=3N$ which implies that the total
coordination does not fluctuate, and the average coordination number is
exactly 6 for any disorder realization. Thus, if there is a spatial region in
the random VD with above average connectivity, there must be another region
with below average connectivity to compensate. The topological constraint thus
introduces anti-correlations in the connectivity disorder which lead to the
more rapid decay of $\sigma^{2}_{Q}$.

Repeating the derivation of the Harris criterion with the relation
$\sigma^{2}_{Q}\sim L^{-(d+1)}$ rather than the central limit theorem result
$\sigma^{2}_{Q}\sim L^{-d}$ leads to the inequality $(d+1)\nu>2$ for the
stability of a clean critical point [42].

##  Appendix B Maximum Entropy Method

The analytic continuation of the scalar susceptibility from imaginary to real
frequencies is given by Eq. (14) which yields the Matsubara susceptibility
$\tilde{\chi}_{\rho\rho}(i\omega_{n})$ as an integral over the real-frequency
spectral function $\chi^{{}^{\prime\prime}}_{\rho\rho}(\omega)$. The inversion
of this relation is an ill-posed problem and very sensitive to Monte Carlo
noise. We can overcome this difficulty by using the maximum entropy method
[47]. The idea is to use Bayesian inference to determine the most probable
spectral density for a given set of Monte Carlo data for the imaginary-time
susceptibility. Mathematically, the problem reduces to a minimization problem
for the cost function

$Q=\Delta-\alpha S~{}$ (22)

The first term in the cost function,

$\Delta=\left(\tilde{\chi}_{\rho\rho}-K\chi^{{}^{\prime\prime}}_{\rho\rho}\right)\Sigma^{-1}\left(\tilde{\chi}_{\rho\rho}-K\chi^{{}^{\prime\prime}}_{\rho\rho}\right)~{},$
(23)

measures how well the Matsubara susceptibility corresponding to
$\chi^{{}^{\prime\prime}}_{\rho\rho}(\omega)$ reproduces the input data. Here,
$K$ is a discretized version of the integration kernel
$2\omega/(\omega_{m}^{2}+\omega^{2})$ in Eq. (14). It is given by

$K(\omega,\omega_{m})=\frac{\Delta\tau\sinh(\Delta\tau\omega)}{\cosh(\Delta\tau\omega)-\cos(\Delta\tau\omega_{m})}$
(24)

where $\Delta\tau$ is the imaginary time step [67, 35]. The matrix
$\Sigma_{mn}=\left<\tilde{\chi}_{\rho\rho}(i\omega_{m})\tilde{\chi}_{\rho\rho}(i\omega_{n})\right>-\left<\tilde{\chi}_{\rho\rho}(i\omega_{m})\right>\left<\tilde{\chi}_{\rho\rho}(i\omega_{n})\right>$
contains the covariance matrix elements of the scalar susceptibility Monte
Carlo data.

The second term in the cost function, the entropy $S$, smoothes the resulting
spectral function and prevents the overfitting of the Monte Carlo noise. It is
defined as

$S=-\sum_{\omega}\chi^{{}^{\prime\prime}}_{\rho\rho}(\omega)\ln{\chi^{{}^{\prime\prime}}_{\rho\rho}(\omega)}.$
(25)

The free parameter $\alpha$ controls the relative weights of the two cost
function terms. The optimum value of this parameter can be obtained using a
particular version of the L-curve method [68, 69]. It involves maximizing the
curvature $k=d^{2}\Delta/d(\ln{\alpha})^{2}$. This maximum marks the crossover
from fitting the actual information in the Monte Carlo data to fitting the
noise. More details of our implementation of the maximum entropy method are
described in Refs. [34, 35].

## References

  * Guo  _et al._ [1996] M. Guo, R. N. Bhatt, and D. A. Huse, Phys. Rev. B 54, 3336 (1996).
  * Rieger and Young [1996] H. Rieger and A. P. Young, Phys. Rev. B 54, 3328 (1996).
  * Vojta and Schmalian [2005] T. Vojta and J. Schmalian, Phys. Rev. B 72, 045438 (2005).
  * 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).
  * Altman  _et al._ [2004] E. Altman, Y. Kafri, A. Polkovnikov, and G. Refael, Phys. Rev. Lett. 93, 150402 (2004).
  * Altman  _et al._ [2008] E. Altman, Y. Kafri, A. Polkovnikov, and G. Refael, Phys. Rev. Lett. 100, 170402 (2008).
  * Vojta  _et al._ [2011] T. Vojta, J. A. Hoyos, P. Mohan, and R. Narayanan, J. Phys.: Condens. Matter 23, 094206 (2011).
  * Vojta [2003] T. Vojta, Phys. Rev. Lett. 90, 107202 (2003).
  * Hoyos and Vojta [2008] J. A. Hoyos and T. Vojta, Phys. Rev. Lett. 100, 240601 (2008).
  * Motrunich  _et al._ [2000] O. Motrunich, S.-C. Mau, D. A. Huse, and D. S. Fisher, Phys. Rev. B 61, 1160 (2000).
  * Vojta and Hoyos [2014] T. Vojta and J. A. Hoyos, Phys. Rev. Lett. 112, 075702 (2014).
  * Vojta [2006] T. Vojta, J. Phys. A: Math. Gen. 39, R143 (2006).
  * Vojta [2013a] T. Vojta, AIP Conf. Proc. 1550, 188 (2013a).
  * Vojta [2019] T. Vojta, Ann. Rev. Condens. Mat. Phys. 10, 233 (2019).
  * Ubaid-Kassis  _et al._ [2010] S. Ubaid-Kassis, T. Vojta, and A. Schroeder, Phys. Rev. Lett. 104, 066402 (2010).
  * Shi  _et al._ [2014] X. Shi, P. V. Lin, T. Sasagawa, V. Dobrosavljević, and D. Popović, Nat. Phys. 10, 437 (2014).
  * 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).
  * Shen  _et al._ [2016] S. Shen, Y. Xing, P. Wang, H. Liu, H. Fu, Y. Zhang, L. He, X. C. Xie, X. Lin, J. Nie, and J. Wang, Phys. Rev. B 94, 144517 (2016).
  * Xing  _et al._ [2017] Y. Xing, K. Zhao, P. Shan, F. Zheng, Y. Zhang, H. Fu, Y. Liu, M. Tian, C. Xi, H. Liu, J. Feng, X. Lin, S. Ji, X. Chen, Q.-K. Xue, and J. Wang, Nano Lett. 17, 6802 (2017).
  * Wang  _et al._ [2017] R. Wang, A. Gebretsadik, S. Ubaid-Kassis, A. Schroeder, T. Vojta, P. J. Baker, F. L. Pratt, S. J. Blundell, T. Lancaster, I. Franke, J. S. Möller, and K. Page, Phys. Rev. Lett. 118, 267202 (2017).
  * Saito  _et al._ [2018] Y. Saito, T. Nojima, and Y. Iwasa, Nat. Commun. 9, 778 (2018).
  * Zhang  _et al._ [2019] C. Zhang, Y. Fan, Q. Chen, T. Wang, X. Liu, Q. Li, Y. Yin, and X. Li, NPG Asia Mater. 11, 76 (2019).
  * Lewellyn  _et al._ [2019] N. A. Lewellyn, I. M. Percher, J. Nelson, J. Garcia-Barriocanal, I. Volotsenko, A. Frydman, T. Vojta, and A. M. Goldman, Phys. Rev. B 99, 054515 (2019).
  * Reiss  _et al._ [2021] P. Reiss, D. Graf, A. A. Haghighirad, T. Vojta, and A. I. Coldea, Phys. Rev. Lett. 127, 246402 (2021).
  * Kaur  _et al._ [2022] S. Kaur, H. K. Kundu, S. Kumar, A. Dogra, R. Narayanan, T. Vojta, and A. Bid, Temperature dependent cloaking of the quantum griffiths singularity in LaScO3/SrTiO3 heterostructures (2022), arXiv:2206.10215 [cond-mat.supr-con] .
  * Shimano and Tsuji [2020] R. Shimano and N. Tsuji, Ann. Rev. Condens. Mat. Phys. 11, 103 (2020).
  * Pekker and Varma [2015] D. Pekker and C. Varma, Ann. Rev. Condens. Mat. Phys. 6, 269 (2015).
  * Podolsky  _et al._ [2011] D. Podolsky, A. Auerbach, and D. P. Arovas, Phys. Rev. B 84, 174522 (2011).
  * Podolsky and Sachdev [2012] D. Podolsky and S. Sachdev, Phys. Rev. B 86, 054508 (2012).
  * Gazit  _et al._ [2013a] S. Gazit, D. Podolsky, and A. Auerbach, Phys. Rev. Lett. 110, 140401 (2013a).
  * Puschmann  _et al._ [2020] M. Puschmann, J. Crewse, J. A. Hoyos, and T. Vojta, Phys. Rev. Lett. 125, 027002 (2020).
  * Crewse and Vojta [2021] J. Crewse and T. Vojta, Phys. Rev. B 104, 014511 (2021).
  * Puschmann  _et al._ [2021] M. Puschmann, J. C. Getelina, J. A. Hoyos, and T. Vojta, Ann. Phys. 435, 168526 (2021), special issue on Localisation 2020.
  * Anderson [1958] P. W. Anderson, Phys. Rev. 109, 1492 (1958).
  * Okabe  _et al._ [2000] A. Okabe, B. Boots, K. Sugihara, and S. N. Chiu,  _Spatial Tessellations. Concepts and Applications of Voronoi Diagrams_ (Wiley, Chichester, 2000).
  * Weaire and Rivier [1984] D. Weaire and N. Rivier, Contemp. Phys. 25, 59 (1984).
  * Aurenhammer [1991] F. Aurenhammer, ACM Comput. Surv. 23, 345 (1991).
  * Harris [1974] A. B. Harris, J. Phys. C 7, 1671 (1974).
  * Barghathi and Vojta [2014] H. Barghathi and T. Vojta, Phys. Rev. Lett. 113, 120602 (2014).
  * Janke and Villanova [2002] W. Janke and R. Villanova, Phys. Rev. B 66, 134208 (2002).
  * Lima  _et al._ [2008] F. W. S. Lima, U. M. S. Costa, and R. N. Costa Filho, Phys. A: Stat. Mech 387, 1545 (2008).
  * Puschmann  _et al._ [2015] M. Puschmann, P. Cain, M. Schreiber, and T. Vojta, Eur. Phys. J. B 88, 314 (2015).
  * Wallin  _et al._ [1994] M. Wallin, E. S. Sørensen, S. M. Girvin, and A. P. Young, Phys. Rev. B 49, 12115 (1994).
  * Jarrell and Gubernatis [1996] M. Jarrell and J. Gubernatis, Phys. Rep. 269, 133 (1996).
  * Vojta  _et al._ [2016] T. Vojta, J. Crewse, M. Puschmann, D. Arovas, and Y. Kiselev, Phys. Rev. B 94, 134501 (2016).
  * Weichman and Mukhopadhyay [2008] P. B. Weichman and R. Mukhopadhyay, Phys. Rev. B 77, 214516 (2008).
  * Tanemura  _et al._ [1983] M. Tanemura, T. Ogawa, and N. Ogita, J. Comput. Phys. 51, 191 (1983).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
  * Wolff [1989] U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
  * Vojta and Sknepnek [2006] T. Vojta and R. Sknepnek, Phys. Rev. B 74, 094415 (2006).
  * Ballesteros  _et al._ [1998] 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 (1998).
  * 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).
  * Cardy [1988] J. L. Cardy, ed.,  _Finite-size scaling_ (North-Holland, 1988).
  * Campostrini  _et al._ [2006] M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B 74, 144506 (2006).
  * Altman and Auerbach [2002] E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 (2002).
  * Pekker  _et al._ [2012] D. Pekker, B. Wunsch, T. Kitagawa, E. Manousakis, A. S. Sørensen, and E. Demler, Phys. Rev. B 86, 144527 (2012).
  * Vojta [2013b] M. Vojta, Phys. Rev. Lett. 111, 097202 (2013b).
  * Gurarie and Chalker [2002] V. Gurarie and J. T. Chalker, Phys. Rev. Lett. 89, 136801 (2002).
  * Gurarie and Chalker [2003] V. Gurarie and J. T. Chalker, Phys. Rev. B 68, 134207 (2003).
  * Note [1] We also emphasize that extracting subleading corrections from fits is more difficult numerically than determining the leading terms, and additional uncertainties arise if more than one correction term contributes.
  * Beattie-Hauser and Vojta [2023] R. Beattie-Hauser and T. Vojta, arXiv:2311.07457 (2023), unpublished.
  * Chayes  _et al._ [1986] J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, Phys. Rev. Lett. 57, 2999 (1986).
  * Gazit  _et al._ [2013b] S. Gazit, D. Podolsky, A. Auerbach, and D. P. Arovas, Phys. Rev. B 88, 235108 (2013b).
  * Hansen and OĹeary [1993] P. C. Hansen and D. P. OĹeary, SIAM J. Sci. Comput. 14, 1487 (1993).
  * Bergeron and Tremblay [2016] D. Bergeron and A.-M. S. Tremblay, Phys. Rev. E 94, 023303 (2016).