Weak Versus Strong Disorder Superfluid-Bose Glass Transition in One Dimension

TL;DR

This study uses advanced simulations to analyze the superfluid to Bose glass transition in one-dimensional attractive hard-core bosons, revealing a critical line with two regimes distinguished by interaction strength and disorder effects.

Contribution

It identifies a Berezinskii-Kosterlitz-Thouless critical line with universal and non-universal regimes, highlighting the role of weak links and power-law distributions in the transition.

Findings

Universal critical Luttinger parameter $K_c=3/2$ at weak attraction.

Non-universal $K_c>3/2$ at stronger attraction with weak-link physics.

Power-law distributed weak links with a variable exponent $eta$ in strong disorder regime.

Abstract

Using large-scale simulations based on matrix product state and quantum Monte Carlo techniques, we study the superfluid to Bose glass-transition for one-dimensional attractive hard-core bosons at zero temperature, across the full regime from weak to strong disorder. As a function of interaction and disorder strength, we identify a Berezinskii-Kosterlitz-Thouless critical line with two different regimes. At small attraction where critical disorder is weak compared to the bandwidth, the critical Luttinger parameter $K_c$ takes its universal Giamarchi-Schulz value $K_{c}=3/2$. Conversely, a non-universal $K_c>3/2$ emerges for stronger attraction where weak-link physics is relevant. In this strong disorder regime, the transition is characterized by self-similar power-law distributed weak links with a continuously varying characteristic exponent $\alpha$.

Figures

Equations

$$P(J)\propto J^{\alpha },$$

$${\cal{H}}=\sum_{i}\Bigl{[}-\frac{1}{2}\left(b_{i}^{\dagger}b_{i+1}^{{\phantom{\dagger}}}+b_{i}^{\phantom{\dagger}}b_{i+1}^{\dagger}\right)+\Delta n_{i}n_{i+1}-\mu_{i}n_{i}\Bigr{]}.$$

$${\cal{H}}=\sum_{i}\Bigl{[}-\left(S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}\right)+\Delta S_{i}^{z}S^{z}_{i+1}-h_{i}S_{i}^{z}\Bigr{]}.$$

$$G_{ij}\equiv ⟨b_i^{†}{b}_j^{}⟩\propto |i-j{|}^{-\frac{1}{2K}},$$

$$K=\pi \sqrt{\kappa {\rho }_{\mathrm{SF}}}.$$

$$f_{\mathrm{BKT}}(L/{\ell }_0)=1+\frac{1}{K_c\ln (L/{\ell }_0)}.$$

$$\ln \left(P[\ln \mathcal{O}]\right)\sim (1+\alpha )\ln \mathcal{O},$$

$$\mathcal{H}=\sum_{i}\Bigl{[}-\frac{1}{2}b_{i}^{\dagger}b_{i+1}+\mathrm{H.c.}+\Delta n_{i}n_{i+1}+\mu_{i}n_{i}\Bigr{]},$$

$$⟨{\mathcal{H}}^2-⟨\mathcal{H}{⟩}^2⟩<ϵL,$$

$$G_{ij}=⟨b_i^{†}{b}_j⟩,$$

$$\displaystyle G_{ij}=\rho_{0}B_{0}\Big{[}\tilde{d}(i,j)\Big{]}^{-1/2K},$$

$$\tilde{d}(i,j)\equiv \frac{d(i+j|2L)d(i-j|2L)}{{\rho }_0^{-1}\sqrt{d(2i|2L)d(2j|2L)}},$$

$$d(i|L)=\frac{L}{\pi }|\sin (\pi i/L)|.$$

$$K=\frac{\pi }{2(\pi -\arccos \Delta )}.$$

Full text

Introduction — The understanding of quantum interacting disordered systems represents an extraordinary challenge in condensed matter and quantum statistical physics. This has been recently discussed in the context of the transition to many-body localization (MBL) Basko et al. (2006); Nandkishore and Huse (2015), for which a precise description is missing Luitz et al. (2015); Potter et al. (2015); Khemani et al. (2016). Besides all these efforts made to understand such a localization transition at high energy, there are still fundamental open issues regarding the low-energy part of the problem.

The most dramatic situation occurs in one dimension where quantum fluctuations plays a major role Giamarchi (2004). In the absence of disorder, interacting bosons may stabilize a quasi-ordered state at zero temperature, with a finite superfluid (SF) density and power-law decaying off-diagonal correlations $G_{ij}=\langle b^{\dagger}_{i}b^{\vphantom{\dagger}}_{j}\rangle\propto r^{-\frac{1}{2K}}$ at large $r\equiv|i-j|$, $K$ being the Luttinger liquid (LL) parameter Giamarchi (2004). In the presence of disorder, such a SF ground state is generally expected to be unstable towards localization, thus forming the so-called Bose glass (BG) phase Giamarchi and Schulz (1987, 1988); Fisher et al. (1989): an inhomogeneous gapless compressible fluid with exponentially suppressed correlations. Disordered interacting chains can be experimentally realized using spin ladders Klanjšek et al. (2008); Jeong et al. (2016), ultra-cold atoms Lucioni et al. (2011); D’Errico et al. (2014); Schreiber et al. (2015) and arrays of disordered Josephson junctions Haviland et al. (2001); Cedergren et al. (2017); Bard et al. (2017). Broadly speaking, Luttinger liquid physics is expected to occur in various contexts Cazalilla et al. (2011), e.g. quantum wires such as carbon nanotubes Bockrath et al. (1999); Ishii et al. (2003), fractional quantum Hall edge states Wen (1990), organic conductors Schwartz et al. (1998) and 4He confined in nanopores Del Maestro et al. (2011).

In their seminal work, based on a LL description and renormalization group (RG) calculations at weak disorder, Giamarchi and Schulz (GS) Giamarchi and Schulz (1987, 1988) have shown that a Berezinskii-Kosterlitz-Thouless (BKT) quantum phase transition should occur for a small but finite disorder strength, provided the clean LL parameter $K_{0}>3/2$. A recent two-loop calculation confirmed this scenario, with a universal jump at the BKT transition $K_{c}=3/2$ Ristivojevic et al. (2012, 2014). Below this value, SF is destroyed by quantum phase slips.

In order to go beyond this perturbatively weak disorder regime, an alternative picture has emerged Altman et al. (2004, 2008, 2010); Pielawa and Altman (2013), based on a strong disorder scenario where the localization transition is predicted to occur for $K_{c}>3/2$ with a different mechanism driven by rare weak links that effectively cut the disordered chain. In this approach, whose asymptotic validity relies on growing disorder under RG, a power-law tail for the distribution of Josephson couplings across weak links is expected

[TABLE]

with $\alpha>0$ ($<0$) in the SF (BG), and $\alpha=0$ at the transition. In a similar spirit, building on the idea of weak links, Pollet and co-workers Pollet et al. (2013, 2014); Yao et al. (2016) suggested a third scenario. Connecting GS to strong disorder physics, they included a renormalization of weak links à la Kane-Fisher Kane and Fisher (1992) and proposed a new “scratched-XY” universality class to describe the SF-BG transition for $K_{c}>3/2$. At weak-link criticality, their approach predicts that the critical LL parameter $K_{c}=\tilde{\alpha}/(\tilde{\alpha}-1)$ is governed by the non-universal bare exponent $\tilde{\alpha}$ of classical Josephson couplings Eq. (1). Despite several numerical studies Bouzerar and Poilblanc (1994); Schmitteckert et al. (1998); Hrahsheh and Vojta (2012); Pielawa and Altman (2013); Yao et al. (2016); Gerster et al. (2016), the extension of the weak disorder regime remains unclear, and there is still no real consensus concerning the scenario of the SF-BG transition at strong disorder.

In this Letter, we elucidate this controversy using high-precision, large-scale numerical simulations of one-dimensional hard-core bosons in a random potential with nearest-neighbor attraction. Such a model is relevant to recent experiments using spin ladders in solid state systems Klanjšek et al. (2008); Jeong et al. (2016).

In the regime where the clean LL parameter $K_{0}>3/2$, we map out the entire phase diagram of model Eq. (2) as a function of interaction and disorder strengths (see Fig. 1). We clearly identify a BKT transition line separating SF and BG regimes by a jump $K_{c}$ of the LL parameter. While $K_{c}=3/2$ for weak attraction and small disorder, the situation is qualitatively different when disorder exceeds the bandwidth, and weak links start to proliferate. At criticality we clearly find $K_{c}>3/2$ and scale-invariant power-law distributions for the SF density and the off-diagonal correlations. We strikingly observe a similar power-law tail exponent $\alpha$ for both quantities, which validates a weak-link scenario. Our results for $\alpha$ and $K_{c}$ along the critical line are different from the aforementioned predictions Pielawa and Altman (2013); Yao et al. (2016). Fig. 1 summarizes our results.

Model and numerical approaches — We study hard-core bosons on a chain, governed by the Hamiltonian

[TABLE]

Hard-core bosons can be exactly mapped onto spin-1/2 using a Matsubara-Matsuda transformation, yielding the corresponding XXZ chain in a magnetic field

[TABLE]

The random fields $h_{i}={\mu_{i}-\Delta}$ are i.i.d. random variables box-distributed $\in[-h,h]$. The half-bandwidth ${\rm{BW}}=1+\Delta$ is the relevant energy scale: in the absence of disorder, the XXZ chain is a polarized ferromagnet (an insulator in the bosonic language) for fields above ${\rm{BW}}$ and below $-{\rm{BW}}$. In the presence of disorder, if $h>\rm{BW}$ (above the dashed line in Fig. 1 top) there is a non-zero probability $P(r)\sim\left(1-{\rm{BW}}/h\right)^{r}$ to find a locally insulating region of length $r$, across which the SF response is exponentially suppressed, thus defining a weak link. In this strong disorder regime, it is already known that the classical local order parameter is power law-distributed Lemarié et al. (2013). There, rare disorder configurations may qualitatively change the physics. Conversely, $h<\rm{BW}$ defines the weak disorder regime.

We have explored the ground state properties of model (2) with two state-of-the-art many-body numerical techniques. The matrix product state (MPS) formulation of the density-matrix renormalization group Schollwöck (2011) has been used to capture the off-diagonal correlator $G_{ij}=\langle b^{\dagger}_{i}b^{\vphantom{\dagger}}_{j}\rangle$ for open chains up to $L=300$ sites. We have further used the quantum Monte Carlo (QMC) algorithm in its stochastic series expansion version Syljuåsen and Sandvik (2002) to estimate the SF density $\rho_{\rm SF}$ and the compressibility $\kappa$ for closed rings up to $L=512$. In the latter case, ground state convergence is achieved at very low temperature using the $\beta$-doubling scheme Sandvik (2002). Disorder averaging is performed over a large number of independent samples, from several hundred up to a few thousand. Numerical details are discussed in the Supplementary Material sup .

Luttinger parameter and BKT transition — We first focus on the LL parameter $K$ which can be computed using MPS from the one-body density matrix Cazalilla (2004)

[TABLE]

and with QMC using the compressibility $\kappa$ and the SF density $\rho_{\rm SF}$ through the hydrodynamic relation Giamarchi (2004)

[TABLE]

In the clean case $K_{0}(\Delta,h=0)=(\pi/2)(\pi-\arccos\Delta)^{-1}$ is exactly known from the Bethe ansatz solution of the XXZ chain Yang and Yang (1966) (dashed blue line in the lower panel of Fig. 1). We therefore focus on the attractive part of the phase diagram $\Delta\leq-0.5$ where $K_{0}\geq 3/2$ such that a disorder-induced SF-BG transition is expected.

QMC results are displayed in Fig. 2 panel (a) for both $\rho_{\rm SF}$ and $\kappa$ as a function of disorder strength $h$ for interaction $\Delta=-3/4$. In this case, a transition occurs (discussed below) at $h_{c}=0.375(15)$ between a SF regime at weak disorder and a BG above $h_{c}$. While the SF response vanishes in the BG at the thermodynamic limit $L\to\infty$, the compressibility remains smooth across the transition, with almost no size dependence. On the other hand $\rho_{\rm SF}(L)$, and therefore also $K(L)$ from Eq. (5), displays quite strong finite size effects which are naturally expected for a BKT transition with logarithmic corrections at criticality of the form $K(L)=K_{c}\times f_{\rm BKT}(L/\ell_{0})$ Weber and Minnhagen (1988); Laflorencie et al. (2001); Hsieh et al. (2013), with $\ell_{0}$ a non-universal length scale:

[TABLE]

Assuming such a scaling form Eq. (6) sup , the SF-BG transition is clearly detected sup . In Fig. 2 (b), the transitions for various values of $\Delta$ are characterized by different crossings obtained with a single free parameter, namely the $O(1)$ length scale $\ell_{0}$, as indicated on the plot. From such an analysis, we can extract both $h_{c}$ and $K_{c}$ which are reported in Fig. 1. The critical boundary $h_{c}$ vs. $\Delta$ shows two distinct behaviors: (i) in the strong disorder regime $h_{c}>\,$BW, $K_{c}>3/2$ is not universal ; (ii) in the weak disorder regime $h_{c}<\,$BW, the LL parameter takes the universal value $K_{c}=3/2$ in agreement with GS scenario. While it is more difficult to locate the SF-BG critical point from MPS simulations, the one-body density matrix $G$ (4) has a decay controlled by $K$ 111The LL parameter $K$ from QMC is a bulk quantity which displays logarithmic corrections of the form Eq. (6) at criticality whereas $K$ extracted from the decay of $G$ is an exponent whose finite size corrections are not known. (see also sup ), whose value is in very good agreement with QMC estimates, as shown in Fig. 1 (green symbols, bottom).

Critical distributions: emergence of weak-link physics — In order to further explore the qualitative differences between the two critical regimes, it is very instructive to go beyond disorder averages and investigate the critical fluctuations. Two quantities of significant interest are the distributions of the one-body density matrix $G$, see Eq. (4), (computed from MPS) and of the superfluid density $\rho_{\mathrm{SF}}$ (obtained from QMC). Using MPS, we have considered the distribution $P_{r}[\ln\tilde{G}]$ at fixed $r\equiv|i-j|$ in the bulk of the system 222$r\gg 1,|i-L/2|\ll L/2,|j-L/2|\ll L/2$. Here $\ln{\tilde{G}}=\ln G-{\overline{\ln G}}$, where ${\overline{\cdots}}$ denotes disorder averaging. The subtraction of the typical value $\overline{\ln G}$ means that we are considering the fluctuations of $G$ on top of the characteristic decay $\propto r^{-1/2K}$. These distributions are self-similar for different $r$ sup , which allows us to consider the averaged distribution $Q\equiv\overline{P_{r}}$ over $20$ different values of $r\in[20,40]$. Using QMC, we study the distribution $P[\ln{\tilde{\rho}}_{\rm SF}]$ where $\ln{\tilde{\rho}}_{\rm SF}=\ln{{\rho}}_{\rm SF}-{\overline{\ln{{\rho}}_{\rm SF}}}$, which again explores the fluctuations around the typical SF response.

Representative results are shown in Fig. 3 for four different critical points. In the case of GS criticality (rightmost panel: $\Delta=-0.625,h_{c}=0.32(3),K_{c}=1.48(4)$), the two distributions are clearly distinct and their narrowing with increasing system size signals self-averaging. However, this is no longer the case in the strong disorder regime. Instead of self-averaging, we observe self-similar distributions with exponential tails

[TABLE]

at low values of ${\cal O}={\tilde{G}}$ or ${\tilde{\rho}}_{\rm SF}$, which corresponds to the aforementioned power-law distribution Eq. (1). Moreover, we strikingly observe a matching between the distribution of ${\tilde{G}}$ and ${\tilde{\rho}}_{\rm SF}$, which are independently obtained. The agreement improves when the critical disorder $h_{c}$ grows compared to the BW (see Fig. 3). This result can be interpreted as a smoking gun for weak-link physics in the strong disorder regime, characterized by a well-defined scale-invariant exponent $\alpha$.

The weak link exponent $\alpha$ can be estimated from both $G$ (MPS) and $\rho_{\mathrm{SF}}$ (QMC) data using Eq. (1). Fig. 1 (red symbols, bottom) shows quantitative agreement, in particular at large $h_{c}/{\rm{BW}}$. When approaching the GS regime, weak-link physics becomes irrelevant, and $\alpha$ diverges. The color map in the phase diagram (Fig. 1 top) shows the evolution of $\alpha$ extracted from $G$ for $L=100$ chains. Note that in the BG regime this map has to be understood as a finite-size snapshot. Indeed, above the critical line self-averaging and self-similarity are lost and the distributions get broader with increasing system sizes sup . This is true throughout the entire BG regime, thus supporting the fact that despite the two different critical mechanisms, there is only a single BG phase in which distributions get flatter, with an asymptotic behavior $\alpha\to-1$, in agreement with an infinite disorder scenario Altman et al. (2010).

Conclusion. — Our numerical results convincingly show that the SF-BG transition for attractive hard-core bosons in a random potential follows two different routes, while BKT universality is maintained. As summarized in Fig. 1, the GS scenario with a universal jump of the LL parameter $K_{c}=3/2$ occurs over a finite region of the phase diagram which corresponds to a critical disorder strength $h_{c}$ smaller than the bandwidth BW. In this weak disorder scenario, self-averaging is observed at the transition. On the other hand, for stronger attraction the bandwidth is smaller, and when $h>{\rm{BW}}$ rare insulating regions start to appear, yielding weak links across the system. In this strong disorder regime, the SF-BG transition occurs with a non-universal (interaction-dependent) jump of $K_{c}>3/2$. Moreover, the critical fluctuations have a power-law distribution Eq. (1) with a single exponent $\alpha$ governing both the SF density $\rho_{\mathrm{SF}}$ and the one-body density matrix $G$. Contrary to the strong disorder RG approach by Altman et al. Altman et al. (2010), the weak-link exponent $\alpha$ continuously varies along the critical line, and diverges when approaching the GS regime. Our numerics also contrast with the “scratched-XY” universality of Pollet et al. Yao et al. (2016) where the LL parameter is predicted to obey $K_{c}={\rm{max}}(\frac{3}{2},{\frac{\tilde{\alpha}}{\tilde{\alpha}-1}})$, which disagrees with our results, assuming the bare exponent $\tilde{\alpha}=\alpha$ 333We justify this through the observed scale invariance at criticality sup .

Our conclusions clearly call for new analytical approaches to describe the strong disorder regime, as well as for the boundary between weak link and GS criticalities. As a final comment, one expects that the richness of the critical properties at play for the zero-temperature localization transition may influence the MBL physics at higher energy in the attractive regime we have studied. Moreover, weak link physics being a central concept for understanding rare regions at play for the MBL problem Agarwal et al. (2017); Nahum et al. (2017), we have provided a quantitative proof for such bottlenecks in the ground state. This should pave the way towards a better understanding of rare regions effects in quantum disordered systems. Experimental verification of our results should be within reach using for instance solid state spin ladder systems or disordered cold atom setups.

Appendix A Supplementary material

Appendix B Technical details (MPS)

In the following, we will outline the numerical approach for the matrix product state (MPS) simulations. As a reminder, we seek to obtain the ground state at half filling of the following Hamiltonian:

[TABLE]

where $b_{i}^{(\dagger)}$ destroys (creates) a hard-core boson at site $i$ on a chain of length $L$, $n_{i}\equiv b_{i}^{\dagger}b_{i}$ and $\Delta$ is the strength of nearest-neighbour interactions. The on-site disorder is given by $\mu_{i}\in[-h,h]$ with a uniform distribution of the disorder. MPS calculations are performed in the (conserved) sector $\sum_{i}\langle n_{i}\rangle=L/2$, such that exact half-filling is imposed. We use Open Source Matrix Product States (OSMPS) Wall and Carr (2012); Jaschke et al. (2017). Since the ground state algorithm uses a variational approach starting from an initial educated guess, convergence to the ground state is not formally guaranteed. In particular, there is the risk of the algorithm getting “stuck” in a metastable state, which is known to affect disordered systems that may have several local minima of the energy. To alleviate this risk, we perform as many DMRG sweeps as necessary up to a maximum of $10$ (for a clean system, usually only a few sweeps are necessary), using a maximum bond dimension of $\chi=300$ until the following criterion is met:

[TABLE]

where the tolerance per site $\epsilon=10^{-10}$. Numerically, our results are obtained by finding the ground state of (S1) using these parameters for a total of $R$ realizations using an embarrassingly parallel approach and storing the full $G_{ij}\equiv\langle b_{i}^{\dagger}b_{j}\rangle$ for each realization.

Determining the Luttinger liquid parameter $K$.— We describe in detail the procedure to determine the decay of the single-particle density matrix or correlator $G$:

[TABLE]

which in the superfluid (SF) regime obeys an algebraic decay as $\propto r^{-1/2K}$ with $r\equiv|i-j|$. In a finite system with open (infinite wall) boundary conditions, $G$ will vanish at the boundary of the system. To correct for this, we use a method based on conformal field theory Cazalilla (2004). According to this theory, the value of $G$ in a finite system is described to leading order by:

[TABLE]

where $d$ is the “chord function:”

[TABLE]

$B_{0}$ is a non-universal constant and $\rho_{0}$ is the average density (these constants are not important for our purposes). The numerical procedure to determine $K$ now consists of computing $G_{ij}^{(\nu)}$ for each realization $\nu=1\ldots R$ and all $i,j\in[1,L]$. Next, we compute both the typical value $G_{ij,\mathrm{typ}}\equiv\exp(\frac{1}{R}\sum_{\nu}\ln G_{ij}^{(\nu)})$ and the average $G_{ij,\mathrm{avg}}\equiv\frac{1}{R}\sum_{\nu}G_{ij}^{(\nu)}$. $K$ can now be determined using the Levenberg-Marquardt least-squares fitting algorithm using either the typical or average $G$.

Although the approach (S4) captures the most important finite-size corrections, it is still more accurate in the bulk than near the edges of the system or for very short-range corrections. Hence, we consider only $r>10$ and discard values of $G$ with $i,j$ close to the edge of the system, specifically we consider only $|i-L/2|<L/3$ and $|j-L/2|<L/3$. In the clean case, the exact value for $K$ is known from the Bethe Ansatz solution Yang and Yang (1966):

[TABLE]

Benchmarking to the non-interacting case at $\Delta=-0.75$ gives a numerical result that is approximately $0.1\%$ off from this exact result. An example of this procedure to determine $K$ is shown in Figure S1 for $\Delta=-3/4$ and $h=0.375$, at the transition from SF to BG. For this value and system sizes $L=100$, $L=200$ and $L=300$ we obtain $K=1.76(5)$, $K=1.74(3)$ and $K=1.74(3)$ using $761$, $368$ and $280$ realizations respectively. Although no reliable finite-size scaling procedure is known to extrapolate to $L\rightarrow\infty$, these numerical results (consistent with results for other values of $\Delta,h$) suggest that the finite-size corrections converge to the thermodynamic result more quickly than QMC results using the superfluid density and compressibility, so that the finite-size result for $L=300$ is fairly close to the QMC estimate in the thermodynamic limit (see Fig. 1 of the main manuscript). This also opens up avenues for accurate experimental measurements of $K$ using finite systems yielding only a modest overestimate of $K$, provided an accurate method to measure the correlations $\langle b_{i}^{\dagger}b_{j}\rangle$ is used.

We briefly remark on the difference between using the typical and average $G$ for determining $K$. Numerically, we find that in the SF regime they yield almost identical results. However, in the BG regime they start to deviate and moreover the estimate for $K$ develops a strong system size dependence. The point at which the typical $K_{\mathrm{typ}}$ and average $K_{\mathrm{avg}}$ estimates for $K$ start to deviate from each other yields a qualitative criterion for the SF-BG transition which is in good agreement with estimates of the critical disorder from QMC. However, an accurate determination of the critical disorder would require detailed knowledge of the scaling of $|K_{\mathrm{typ}}-K_{\mathrm{avg}}|$ close to the transition, which is a priori unknown.

As a final remark, we stress that the decay of the correlator as shown in Fig. S1 should not be confused with the distribution of $P_{r}[\ln G]$ as discussed in the main text. The distributions in the latter case are obtained by looking at a fixed $r=|i-j|$ (or fixed $\tilde{d}$, which yields similar distributions for sufficiently small variations of $r$ or $\tilde{d}$) and then collecting all values of $G$ (not just the typical value) within a narrow window around $r$ or $\tilde{d}$. Thus, one takes a small “slice” of values in Fig. S1 and then computes the probability density function $P_{r}[\ln G]$ using $R$ values per each $G_{ij}$ inside the slice. The total distribution $Q=\sum_{r}\overline{P_{r}}$ is then obtained by summing and averaging several such contributions in consecutive slices, while subtracting the typical value to correct for the influence of the decay $\propto r^{-1/2K}$ so that consecutive $P_{r}$’s are self-similar. We illustrate this procedure in Fig. S1 (bottom panel). After all these steps, only a modest number of realizations is required for very accurate distributions (see Fig. 3 of the main paper).

Appendix C Technical details (QMC)

The Quantum Monte Carlo (QMC) simulations have been performed using the Stochastic Series Expansion (SSE) method Syljuåsen and Sandvik (2002). While this is a rather standard technique for non-frustrated quantum spin models in the absence of disorder, in the presence of randomness one has to pay some attention regarding several issues, as discussed in previous works Sandvik (2002); Laflorencie et al. (2004):

Equilibration and thermalization.— In order to achieve a fast equilibration, the $\beta$-doubling scheme Sandvik (2002) is used, cooling down the system following $\beta=3,\,6,\,12,\cdots,3*2^{p}$, with $p_{\rm max}=10-11$, depending on the system size. At each temperature step, $N_{\rm eq}=10^{3}$ MC steps are used.

• 2)

Ground state convergence.— The convergence to $T=0$ properties is also achieved using the $\beta$-doubling procedure, as exemplified below in Fig. S2 for $\Delta=-3/4$ and various disorder strengths.

• 3)

Monte Carlo vs. sample-to-sample fluctuations.— In the strong disorder regime, using $N_{\rm mes}=10^{3}$ MC steps of measurement is enough since the distributions displays fat tails. On the other hand, in order to correctly sample such power-law tails, several thousands of independent disordered samples have been used. In the weak disorder regime, where distributions are much narrower, we have used $10^{4}$ and up to $10^{5}$ MC steps in order to correctly capture self-averaging.

We also note that simulations have been performed in the grand-canonical ensemble, where particle number can fluctuate. Regarding the random field configurations, we did not impose any constraint such as $\sum_{i}h_{i}=0$ over a sample. Instead, each sample has an independent configuration taken in a uniform box distribution in the range $[-h,h]$. We do not expect the distinction between canonical vs. microcanonical ensemble for the disorder Monthus (2004) to change the critical properties, as demonstrated for a related study in two dimensions Álvarez Zúñiga et al. (2015); Ng and Sørensen (2015).

Appendix D Broadening of the superfluid stiffness distribution in the Bose glass regime

As discussed in the main text, the weak-link distributions in the Bose glass regime (at $h>h_{c}$) get broader with increasing system size. This is shown in Figure S3 where QMC results for $P(\ln\rho_{\rm sf})$ are displayed for two representative values of interaction $\Delta=-0.5$ and $-0.875$ at the same disorder strength $h=0.4$. Here we see that both cases are qualitatively similar, despite the fact that

(i) at $\Delta=-0.5$, $h/{\rm{BW}}=0.8$, corresponding to “weak” disorder;

(ii) at $\Delta=-0.875$, $h/{\rm{BW}}=3.2$, which is far above the BW line, i.e. in the “strong” disorder regime.

From the behavior of these distributions, qualitatively similar to what was observed for the Bose glass regime in two dimensions Álvarez Zúñiga et al. (2015), there is no distinction between “weak” and “strong” disorder regime, and we expect only a single Bose glass phase with $\alpha\to-1$ when $L\to\infty$.

Appendix E Critical distributions and power law $\alpha$

As discussed in the main manuscript, the distributions of the superfluid density $\tilde{\rho}_{\mathrm{SF}}$ and correlator $Q[\ln\tilde{G}]$ show two qualitatively different regimes: a self-similar regime if $h_{c}>\rm BW$, and a self-averaging regime for $h_{c}<\rm BW$. In Figures S4-S5 we show two characteristic examples of critical distributions for QMC and MPS data separately. Since $\rho_{\mathrm{SF}}$ and $G_{ij}$ are entirely different quantities, there is no a priori reason to expect the distributions to be similar. That they are, nonetheless, very similar at the transition in the weak-link regime shows the crucial importance of the power-law tails associated with weak links in this regime. Indeed, in the weak disorder regime $h_{c}<\rm BW$ the distributions are very different although still self-averaging for both QMC and MPS data.

It is also instructive to consider the behavior of the power law $\alpha$ near the superfluid-Bose glass transition, as shown in Figure S6. A qualitative difference between the distributions of $\rho_{\mathrm{SF}}$ and $G$ is that the distributions for $\rho_{\mathrm{SF}}$ are self-similar only exactly at the transition, while for MPS data they remain self-similar in the superfluid regime. This means that the value of $\alpha$ shown in the color contour plot of Fig. 1 of the main manuscript should be a good approximation of the value of the power law in the superfluid regime in the thermodynamic limit. While we do not observe a change in the power law tail as a function of $L$ for $Q[\ln G]$ in the BG regime, we do observe a broadening of the peak of the distribution (not shown), which should eventually remove the power law in the limit $L\rightarrow\infty$. On the other hand, the data for $\rho_{\mathrm{SF}}$ shows this trend already for modest $L$ (see Figure S6).

The prediction for the critical $\alpha$ can be compared to the prediction of Yao et al. Yao et al. (2016), assuming a direct correspondence between the bare power law coefficient $\tilde{\alpha}$ of the weak-link theory of Yao *et al. *and our numerical results for $\alpha$. Such a direct correspondence is justified on the basis of the self-similarity of the distributions at criticality that we observe. In Figure S7 we show the comparison between both predictions for the power law $\alpha$ at the critical LL parameter $K_{c}$.

Appendix F Finite size scaling of the superfluid density at criticality

Here we want to discuss the logarithmic corrections at the SF-BG transition. While the standard BKT transition leads to the following finite size corrections Weber and Minnhagen (1988); Hsieh et al. (2013) $f_{\rm BKT}(L/\ell_{0})=1+\frac{1}{2\ln(L/\ell_{0})+\ln\left[\ln(L/\ell_{0})\right]}$, the integration of the Giamarchi-Schulz renormalization equations leads at the lowest order to a slightly different form: $f_{\rm BKT}(L/\ell_{0})=1+\frac{1}{K_{c}\ln(L/\ell_{0})}$.

In the main text (Fig. 2), we have therefore used the latter finite size expression to analyse our QMC data for the Luttinger parameter in the critical region. In Fig. S8 we show two panels for our QMC data at $\Delta=-0.55$ using the aforementioned logarithmic corrections. Given our finite size samples, we have observed that both forms give a similarly good description, and the critical parameters ($h_{c},K_{c}$) remain unchanged within error bars.