Dynamical Conductivity Across The Disorder-Tuned Superconductor-Insulator Transition

A quantum phase transition is a dramatic event marked by large spatial and temporal fluctuations, where one phase of matter with its ground state and tower of excitations reorganizes into a completely different phase. We provide new insight into the disorder-driven superconductor-insulator transition (SIT) in two dimensions, a problem of great theoretical and experimental interest, with the dynamical conductivity \sigma(\omega) and the bosonic (pair) spectral function P(\omega) calculated from quantum Monte Carlo simulations. We identify characteristic energy scales in the superconducting and insulating phases that vanish at the transition due to enhanced quantum fluctuations, despite the persistence of a robust fermionic gap across the SIT. Disorder leads to enhanced absorption in \sigma(\omega) at low frequencies compared to the SIT in a clean system. Disorder also expands the quantum critical region, due to a change in the universality class, with an underlying T=0 critical point with a finite low-frequency conductivity.

The interplay of superconductivity and localization has proven to be a rich and intriguing problem, especially in two dimensions [1][2][3][4][5][6][7] . Both paradigms stand on the shoulders of giants -the BCS theory of superconductivity and the Anderson theory of localization. Yet, when the combined effects of superconductivity and disorder are considered, both paradigms break down, even for s-wave superconductors.
It has been shown [8][9][10] in model fermionic Hamiltonians with attraction between electrons and disorder arising from random potentials, that the single-particle density of states continues to show a hard gap across the disorder-driven quantum phase transition and that pairs continue to survive into the insulating state. The superconducting transition temperature T c , however, does decrease with increasing disorder and vanishes at a critical disorder signaling a superconductor-insulator transition (SIT). These theoretical predictions are supported by scanning tunneling spectroscopy experiments [11][12][13][14] and by magnetoresistance oscillations 6 in disordered thin films.
Recent conductivity measurements at frequencies well within the superconducting gap (0-20 GHz) [15][16][17][18][19][20] have observed low-frequency features that cannot be accounted for by pair-breaking mechanisms. A theoretical understanding of the low-frequency dynamical conductivity is vital for understanding the role of fluctuations and for guiding future experiments that probe the SIT.
The robustness of the single-particle gap across the SIT suggests that the low-energy physics near the SIT can be described by an effective "bosonic" Hamiltonian, the disordered quantum XY model, where the relevant degrees of freedom are the phases of the local superconducting order parameter. This model is also relevant for ultracold atomic gases in optical lattices where the transition is tuned by changing the tunneling of bosons compared to their on-site repulsion [21][22][23][24] . More recently, it has also become possible to include disorder in optical lattices using speckle patterns. By increasing the strength of the disorder potential it could be possible to drive quantum phase transitions from a superfluid to a Bose glass [25][26][27][28] ; our results are also relevant for such experiments.
We 1 op (R) = hc R# c R" i generated in the presence of large disorder, as we now explain.
We show in Fig. 5 that even for 'homogeneous' disorder, that is, an uncorrelated random potential V (R) (Fig. 5a), the pairing amplitude 1 op (R) exhibits an emergent 'granular' structure (shown in Fig. 5b). The system self-organizes into superconducting islands, on the scale of the coherence length, with finite 1 op (R), interspersed with insulating regions where 1 op (R) is negligible. The spatial variations of spectral features (asymmetry and coherence peaks) in this inhomogeneous state were already discussed above in connection with Fig. 4.
The close connection between inhomogeneity and energy gaps is made clear in Fig. 5b,c, which demonstrates two striking facts. We see that (1) there is an energy gap in the LDOS at every site, and (2) small gaps ! dos (R) in the LDOS are spatially correlated with large 1 op (R) SC islands.
A simple way to understand these results is to use the pairingof-exact-eigenstates approach generalized to highly disordered systems 15 . In the limit of weak attraction, pairing leads to a gap in the low-energy DOS in the underlying Anderson insulator and leads to the islands with non-zero 1 op and a small energy gap. On the other hand, the insulating sea corresponds to the higher-energy strongly localized states in the system.
From this perspective one can see that the gap ! dos , observed in the spatially average DOS, initially decreases with increasing disorder owing to a reduction in the DOS near the chemical potential in our model. (In a real material, the coupling will also decrease 29 with disorder.) However, at high disorder, the gap grows (consistent with Fig. 1) like ! dos ⇡ |U |/(2⇠ 2 loc ), where ⇠ loc is the single-particle localization length 15 . This is due to the enhanced effective attraction between fermions confined to a smaller localization volume ⇠ 2 loc . The phase stiffness (or superfluid density) ⇢ s (T = 0), on the other hand, decreases monotonically with disorder as the SC islands become smaller and the Josephson coupling between islands becomes weaker. Thus, even if one starts with a weak-coupling BCS superconductor with ! dos ⌧ ⇢ s , disorder will necessarily Comparison with experiments. We describe the c between our predictions and experiments on the disor SIT in systems such as indium oxide, titanium nitride, an nitride films, for which our theory seems to be the most ap First, let us discuss the insulating side of the SIT. The e a gap in the insulator implies activated transport, cons early measurements on amorphous InO x films 5 . Furtherm is evidence for pairs on the insulating side of the tran specially patterned amorphous bismuth films.
Recent scanning tunnelling microscpy (STM) exper directly relevant to our predictions on the supercondu of the SIT. Experiments on homogeneously disordered T have shown that, whereas T c goes to zero at the SIT, gap ! dos remains finite, in agreement with Fig. 1 33 have all found a persisting up to many times T c . In particular, they marked suppression of the low-energy DOS togeth destruction of coherence peaks above T c , in complete with our predictions.
We hope that future STM experiments will study the anticorrelation that we predict between the heig coherence peaks (associated with large pairing amplitud small energy gaps in the local DOS. The obvious quantu scaling between T c and ⇢ s (0) at the SIT, well studied different systems 34 , also remains to be tested experim s-wave superconducting films.

Conclusion
In conclusion, we have obtained detailed insights and p for observable properties of the highly disordered superc and insulating states in 2D films, and of the transitio these states. Although we focused on s-wave SC film not escaped our attention that aspects of our resu 5t, T = 0, and n = 0.875. Note the emergent 'granular' structure where the pairing amplitude 'self-organizes' into superconducting islands on the scale of the coherence length, even though the 'homogeneous' disorder potential in a varies on the scale of a lattice spacing. c, Local energy gap ! dos (R) from BdG, defined as the smallest energy at which the local DOS is non-zero (N(R,!) > 0.004). Note that this gap is finite everywhere and that the smallest gaps occur on the SC islands defined by the largest pairing amplitude.
1 op (R) = hc R# c R" i generated in the presence of large disorder, as we now explain.
We show in Fig. 5 that even for 'homogeneous' disorder, that is, an uncorrelated random potential V (R) (Fig. 5a), the pairing amplitude 1 op (R) exhibits an emergent 'granular' structure (shown in Fig. 5b). The system self-organizes into superconducting islands, on the scale of the coherence length, with finite 1 op (R), interspersed with insulating regions where 1 op (R) is negligible. The spatial variations of spectral features (asymmetry and coherence peaks) in this inhomogeneous state were already discussed above in connection with Fig. 4.
The close connection between inhomogeneity and energy gaps is made clear in Fig. 5b,c, which demonstrates two striking facts. We see that (1) there is an energy gap in the LDOS at every site, and (2) small gaps ! dos (R) in the LDOS are spatially correlated with large 1 op (R) SC islands.
A simple way to understand these results is to use the pairingof-exact-eigenstates approach generalized to highly disordered systems 15 . In the limit of weak attraction, pairing leads to a gap in the low-energy DOS in the underlying Anderson insulator and leads to the islands with non-zero 1 op and a small energy gap. On the other hand, the insulating sea corresponds to the higher-energy strongly localized states in the system.
From this perspective one can see that the gap ! dos , observed in the spatially average DOS, initially decreases with increasing disorder owing to a reduction in the DOS near the chemical potential in our model. (In a real material, the coupling will also decrease 29 with disorder.) However, at high disorder, the gap grows (consistent with Fig. 1) like ! dos ⇡ |U |/(2⇠ 2 loc ), where ⇠ loc is the single-particle localization length 15 . This is due to the enhanced effective attraction between fermions confined to a smaller localization volume ⇠ 2 loc . The phase stiffness (or superfluid density) ⇢ s (T = 0), on the other hand, decreases monotonically with disorder as the SC islands become smaller and the Josephson coupling between islands becomes weaker. Thus, even if one starts with a weak-coupling BCS superconductor with ! dos ⌧ ⇢ s , disorder will necessarily Comparison with experiments. We describe the connection between our predictions and experiments on the disorder-tuned SIT in systems such as indium oxide, titanium nitride, and niobium nitride films, for which our theory seems to be the most appropriate. First, let us discuss the insulating side of the SIT. The existence of a gap in the insulator implies activated transport, consistent with early measurements on amorphous InO x films 5 . Furthermore, there is evidence for pairs on the insulating side of the transition 8 in specially patterned amorphous bismuth films.
Recent scanning tunnelling microscpy (STM) experiments are directly relevant to our predictions on the superconducting side of the SIT. Experiments on homogeneously disordered TiN films 18 have shown that, whereas T c goes to zero at the SIT, the STM gap ! dos remains finite, in agreement with Fig. 1 33 have all found a pseudogap persisting up to many times T c . In particular, they observe a marked suppression of the low-energy DOS together with a destruction of coherence peaks above T c , in complete agreement with our predictions.
We hope that future STM experiments will study in detail the anticorrelation that we predict between the height of the coherence peaks (associated with large pairing amplitude) and the small energy gaps in the local DOS. The obvious quantum critical scaling between T c and ⇢ s (0) at the SIT, well studied in rather different systems 34 , also remains to be tested experimentally in s-wave superconducting films.

Conclusion
In conclusion, we have obtained detailed insights and predictions for observable properties of the highly disordered superconducting and insulating states in 2D films, and of the transition between these states. Although we focused on s-wave SC films, it has not escaped our attention that aspects of our results bear a 1 op (R) = hc R# c R" i generated in the presence of large disorder, as we now explain.

NATURE PHYSICS
We show in Fig. 5 that even for 'homogeneous' disorder, that is, an uncorrelated random potential V (R) (Fig. 5a), the pairing amplitude 1 op (R) exhibits an emergent 'granular' structure (shown in Fig. 5b). The system self-organizes into superconducting islands, on the scale of the coherence length, with finite 1 op (R), interspersed with insulating regions where 1 op (R) is negligible. The spatial variations of spectral features (asymmetry and coherence peaks) in this inhomogeneous state were already discussed above in connection with Fig. 4.
The close connection between inhomogeneity and energy gaps is made clear in Fig. 5b,c, which demonstrates two striking facts. We see that (1) there is an energy gap in the LDOS at every site, and (2) small gaps ! dos (R) in the LDOS are spatially correlated with large 1 op (R) SC islands.
A simple way to understand these results is to use the pairingof-exact-eigenstates approach generalized to highly disordered systems 15 . In the limit of weak attraction, pairing leads to a gap in the low-energy DOS in the underlying Anderson insulator and leads to the islands with non-zero 1 op and a small energy gap. On the other hand, the insulating sea corresponds to the higher-energy strongly localized states in the system.
From this perspective one can see that the gap ! dos , observed in the spatially average DOS, initially decreases with increasing disorder owing to a reduction in the DOS near the chemical potential in our model. (In a real material, the coupling will also decrease 29 with disorder.) However, at high disorder, the gap grows (consistent with Fig. 1) like ! dos ⇡ |U |/(2⇠ 2 loc ), where ⇠ loc is the single-particle localization length 15 . This is due to the enhanced effective attraction between fermions confined to a smaller localization volume ⇠ 2 loc . The phase stiffness (or superfluid density) ⇢ s (T = 0), on the other hand, decreases monotonically with disorder as the SC islands become smaller and the Josephson coupling between islands becomes weaker. Thus, even if one starts with a weak-coupling BCS superconductor with ! dos ⌧ ⇢ s , disorder will necessarily Comparison with experiments. We describe the connection between our predictions and experiments on the disorder-tuned SIT in systems such as indium oxide, titanium nitride, and niobium nitride films, for which our theory seems to be the most appropriate. First, let us discuss the insulating side of the SIT. The existence of a gap in the insulator implies activated transport, consistent with early measurements on amorphous InO x films 5 . Furthermore, there is evidence for pairs on the insulating side of the transition 8 in specially patterned amorphous bismuth films.
Recent scanning tunnelling microscpy (STM) experiments are directly relevant to our predictions on the superconducting side of the SIT. Experiments on homogeneously disordered TiN films 18 have shown that, whereas T c goes to zero at the SIT, the STM gap ! dos remains finite, in agreement with Fig. 1 33 have all found a pseudogap persisting up to many times T c . In particular, they observe a marked suppression of the low-energy DOS together with a destruction of coherence peaks above T c , in complete agreement with our predictions.
We hope that future STM experiments will study in detail the anticorrelation that we predict between the height of the coherence peaks (associated with large pairing amplitude) and the small energy gaps in the local DOS. The obvious quantum critical scaling between T c and ⇢ s (0) at the SIT, well studied in rather different systems 34 , also remains to be tested experimentally in s-wave superconducting films.

Conclusion
In conclusion, we have obtained detailed insights and predictions for observable properties of the highly disordered superconducting and insulating states in 2D films, and of the transition between these states. Although we focused on s-wave SC films, it has not escaped our attention that aspects of our results bear a FIG. 1. The emergent inhomogeneity of the local pairing amplitude ∆(r) in a disordered superconductor in the left panel and the robustness of the single particle gap 8-10 across the SIT suggests an effective low-energy description in terms of a disordered quantum XY model shown on the right. The quantum phase transition occurs when long range phase coherence is lost between weakly connected "superconducting islands" tuned by the ratio Ec/EJ of charging energy to Josephson coupling as well as by disorder, modeled by removing a fraction p of the Josephson bonds.
anisotropic classical 3D XY model [29][30][31] and simulate the model using Monte Carlo methods. We focus on the behavior of two dynamical quantities of fundamental significance, the conductivity σ(ω) and the boson ("pair") spectral function P (ω) obtained by analytic continuation from imaginary time using the maximum entropy method supplemented by sum rules. Disorder is introduced into the quantum model by breaking bonds ("Josephson couplings") on a 2D square lattice with a probability p. We compare the results of the disorder-driven SIT with the clean system 29,32 , where the SIT is tuned by E c /E J , the charging energy relative to the Josephson coupling.
Our main results are as follows. (1) The conductivity Re σ(ω) in the clean superconductor shows absorption above a threshold ω Higgs that can be associated with the scale of the Higgs (amplitude) mode. As we approach the SIT from the superconducting (SC) side, both the superfluid stiffness ρ s and the Higgs scale fermionic energy gap remains finite across the transition.
(2) In the insulating state of the clean system, we find a threshold ω σ for absorption in Re σ(ω) and show that it is twice the gap ω B in the bosonic spectral function Im P (ω)/ω. We show that both these scales go soft on approaching the SIT from the insulating side. Furthermore, in the insulator, Im σ(ω) becomes negative at low frequencies, indicating "capacitive" response.
(3) The low-frequency spectral weight in Re σ(ω) for the disordered system is greatly enhanced relative to its clean counterpart, so that there is no clear optical gap in the vicinity of the SIT, despite the existence of a non-zero fermonic energy gap. We find that enhanced quantum phase fluctuations and rare regions generate low-frequency spectral weight for ω well below the clean ω Higgs scale in the SC state, and well below the clean ω σ = 2ω B scale on the insulating side. (4) The spectral function Im P (ω)/ω has a characteristic peak in the insulator, whose energy ω B is a measure of the inverse coherence time scale for bosonic (pair) excitations. The vanishing of the superfluid stiffness ρ s on the SC side and the vanishing of ω B from the insulating side are shown to demarcate the quantum critical regime at the SIT for both the clean and the disordered system. (5) The low-frequency conductivity σ * in the quantum critical regime between the SC and the insulator can be estimated meaningfully from the integrated spectral weight over a frequency range of the order of the temperature (see Eq. 18). We find σ * 0.5(4e 2 /h) at the disorder-driven SIT in comparison to σ * 0.4(4e 2 /h) at the SIT in the pure system, in good agreement with recent studies 32-35 of the disorder-free problem. Model: The quantum XY model is equivalent to a Josephson-junction array, with the Hamiltonian where the number operatorn i at site i is canonically conjugate to the phase operatorθ i . Here E c is the charging energy. The Josephson couplings are J ij = E J with probability (1 − p) and J ij = 0 with probability p. The clean system (p = 0) is a coherent superconductor when E J dominates over E c , with phases aligned across all the junctions. However, large E c /E J favors a well-defined number eigenstate, leads to strong phase fluctuations, and drives the system into an insulating state. Thus E c /E J can be used to tune across the SIT in the clean system. A quantum phase transition can also be induced by increasing disorder p (bond dilution) for fixed E c /E J . (Fig. 3(h)). Thus Eq. 1 is a simple yet non-trivial model that describes a disorder-tuned SIT with a dynamical exponent z = 1.
Our results are obtained from calculations of the superfluid stiffness ρ s , the complex conductivity σ(ω), and shows a crossover from "inductive" (ω Im σ(ω) = ρs > 0) to "capacitative" (ω Im σ(ω) < 0) behavior at small ω across the transition. (c) The boson spectral function Im P (ω)/ω, which has a peak centered about zero frequency in the superconductor, develops a characteristic scale ωB in the insulator that grows with disorder.
We use the Kubo formula for the complex conductivity σ(ω) expressed in terms of Λ xx (q = 0, τ ) and transform the imaginary-time QMC results to real frequency using the maximum entropy method (MEM); see Appendix B. We have checked our results extensively using sum rules and compared the MEM results with direct estimates in imaginary time, as described in detail below. Similarly, we use QMC methods to calculate the imaginary time correlation function P (r, τ ) = a † (r, τ )a(0, 0) , where the bosonic creation operator is a † = exp iθ(r, τ ), and we obtain the spectral function Im P (ω) using the MEM. Superconductor: We first discuss the SC and insulating state in both the clean and disordered systems, before turning to the quantum critical point. The SC state is characterized by a non-zero superfluid stiffness ρ s (see Fig. 2). We use our calculated ρ s to test the sum rule for the MEM-derived optical conductivity. The total spectral weight is given by ∞ 0 dω Re σ(ω) = π −k x /2, where −k x is the kinetic energy. We find that ∞ 0 + dω Re σ(ω) (note the lower limit of 0 + ) calculated from the MEM result differs from −k x by an amount that is exactly accounted for by the delta function ρ s δ(ω). We have checked this sum rule both in the clean and the disordered systems (see Appendix B).
In the clean superconductor ( Fig. 3(a)), Re σ(ω) shows finite spectral weight above a threshold. Note that in the bosonic model, the cost of making electron-hole excitations is essentially infinite (i.e., much larger than all scales of interest). Phase fluctuations of the order parameter, Ψ = A exp(iθ), lead to a current j ∼ Im Ψ * ∇Ψ ∼ |A| 2 ∇θ. This then leads to the absorption threshold 36,37 for creating a massive amplitude excitation (Higgs mode) and a massless phase excitation (phonon). Hence, we identify the threshold in Re σ(ω) with the Higgs scale ω Higgs . We emphasize that even though the microscopic model (1) has only phase degrees of freedom, its longwavelength behavior upon coarse-graining contains both amplitude (Higgs) and phase fluctuations (phonons and vortices). In addition, one can show that Re σ(ω) has a ω 5 tail at low energies arising from three-phonon absorption in a clean SC. The large power-law suppression, together with a very small numerical prefactor 38 , however, makes this spectral weight too small to be visible in our numerical results for Re σ(ω).
As E c /E J is tuned to reach the SIT in the clean system, ρ s decreases and vanishes at the transition; see Fig. 2. We also find that the Higgs scale goes soft upon approaching the quantum critical point, as expected.
The disordered SC results differ in several ways from those of the clean system. First, the superfluid stiffness ρ s is reduced by disorder, vanishing at the SIT upon tuning the transition by disorder p. An important difference is the absence of a discernible Higgs threshold in Re σ(ω) for the disordered SC; see Fig. 3(b). Qualitatively we can understand this by the fact that once disorder breaks momentum conservation even single-phonon absorption is permitted and one no longer needs a multi-phonon process for absorption. The effect of long-range Coulomb interactions, which change the phonon dispersion (∼ q) to that of a 2D plasmon (∼  transform in the reactive response Im σ(ω) = ρ s /ω can indeed be measured. In the SC, the finite low-frequency absorption in Re σ(ω) (due to the single-phonon processes discussed above) causes ω Im σ(ω) to deviate from a constant, as is evident in Fig. 4. Our results are qualitatively similar to what has been seen in recent experiments, which, however, have focused on finitetemperature transitions in weakly disordered samples 20 .
Insulator: The clean insulator shows a hard gap in Re σ(ω) with an absorption threshold that we denote by ω σ ; see Fig. 3(e). To gain insight into this gap, we look at the boson spectral function Im P (ω)/ω in Fig. 3(g), which too shows a hard gap ω B , the analog of what was dubbed ω pair in Ref. 10. The simplest process contributing to the conductivity is described diagrammatically as the convolution of two boson Greens functions leading to ω σ = 2ω B as seen in Fig. 2. We also see that both of these energy scales go soft as the SIT is approached from the insulating side. In addition, there is a welldefined peak in Im P (ω)/ω at a characteristic scale ω B (Fig. 3(g)), which also goes soft at the SIT (Fig. 2). In contrast to the hard gap of the clean system, the dirty insulator exhibits absorption down to arbitrarily low frequencies (see Fig. 3(f)), which is, at least in part, due to rare regions. This then raises the question: what is the characteristic energy scale that goes soft as one approaches the SIT from the insulating side? We find that this scale is the location of the low-energy peak at ω B in the boson spectral function Im P (ω)/ω, whose evolution with disorder is most readily seen in the "slingshot-like" plot in Fig. 4(c). The corresponding changes in Re σ(ω) are shown in Fig. 4(a). We also note that there is a marked change in Im σ(ω) across the SIT. We see from Fig. 4(b) that it changes sign at low frequencies from an inductive (Im σ(ω) > 0) to a capacitive (Im σ(ω) < 0) response going through the disorder-tuned SIT. Quantum criticality: We have already discussed the various scales that go soft on approaching the SIT from either side. The results for the clean system, with SIT tuned by E c /E J , are summarized in Fig. 2. We now analyze the results for the disordered system. The finite temperature QMC data, taken at face value, suggest a finite separation between the disorder values at which ρ s goes to zero from the SC side and the characteristic boson scale ω B vanishes from the insulating side; see Fig. 5  (a,b). We emphasize that this intermediate region is not a Bose metal separating the SC and insulator, but rather the quantum critical region. As shown in Fig. 5 (c), the SC transition temperature T c , at which ρ s vanishes, and the crossover scale T * , at which ω B vanishes, define this fan-shaped critical region. (We have used z = 1 in scaling the system size as we go down in temperature in Fig. 5 (c); see Appendix A.) Both T c and T * extrapolate to zero at the same critical disorder p c ≈ 0.337 (for the chosen value of E c /E J ) with the scaling |p − p c | zν where z = 1 and ν = 0.96 ± 0.06.
Finally, we turn to the important question of the universal conductivity at the SIT [29][30][31][32]39 . The d.c. limit requires ω → 0 first and then T → 0, which is not possible when analytically continuing Matsubara data 40 . What we can meaningfully do is to exploit quantum critical scaling and sum rules. The MEM results (i) satisfy the conductivity sum rule, which integrates over all frequencies (see Appendix B), and (ii) are reliable for high frequencies ω > 2πT . Taking the difference of integrated spectral weights, we can reliably estimate σ * = (2πT ) −1 2πT 0 + dω Re σ(ω, T ; p). We may now use the universal scaling form 40 Re σ(ω, Thus σ * is a T -independent universal constant at the quantum critical point p = p c and closely related to the low frequency conductivity measured in experiments. Another estimate of the low-frequency conductivity comes directly from the current correlator σ * Λ = β 2 Λ xx (q = 0, τ = β/2)/π at the largest available value of imaginary time 41 . The σ * estimates obtained by the two methods show good agreement ( Fig. 6(a)) and provide a non-trivial check on the analytic continuation.
In Fig. 6 (b) we plot σ * (T ; p) as a function of p for various temperatures. In the superconductor (p < p c ) the conductivity increases with decreasing T , while the opposite trend is observed in the insulator (p > p c ). Precisely at the SIT p = p c , we find a T -independent crossing point which also allows us to estimate the critical σ * . Another way to scale the data is to plot σ * (T ; p) as a function of the scaling variable |p − p c |T −1/zν . We find data collapse for p c = 0.337 and zν = 0.96 (consistent with Fig. 5) with a critical value of σ * ≈ 0.5σ Q . For a detailed comparison of the critical exponents and σ * with previous results 42 , see Appendix C. Conclusions: We have presented calculations of the complex dynamical conductivity σ(ω) and the boson spectral function P (ω) across the SIT driven by increasing the charging energy E c /E J as well as by increasing disorder p. By comparison of the clean and disordered problems, we see the effect of disorder on the Higgs scale ω Higgs in the superconductor and on the threshold ω σ in the insulator, in generating low frequency weight in absorption in both superconducting and insulating phases, and in expanding the region over which critical fluctuations are observable. In the literature, an insulating phase of bosons with disorder has been referred to as a compressible Bose glass phase away from particle-hole symmetry 21 or an incompressible Mott glass phase with particle-hole symmetry 42 . We work with a particle-hole symmetric system, and while we see evidence of a gap-like scale in the insulator, we also find a low-frequency tail in the absorption, presumably arising from rare regions. In this respect our insulator seems more akin to a Bose glass. It is important to emphasize that the effects we have calculated have required going beyond mean field theories, even those that included emergent granularity due to the microscopic disorder, by focussing on the role of fluctuations of the order parameter. We have calculated the effect of these fluctuations, both amplitude and phase, on experimentally accessible observables using QMC methods coupled with maximum entropy methods, constrained by sum rules. Recently the AdS-CFT holographic mapping has been used to obtain the dynamical conductivity at the disorder-free bosonic quantum critical point 34,35 . Our focus here has been on the evolution of the dynamical quantities in both the phases, superconducting and insulating, and across the disorderdriven SIT, for which the holographic formalism has not yet been developed.
Our calculations have laid the foundation for key signatures in dynamical response functions across quantum phase transitions. Though we have focused on the disorder-driven s-wave SIT in thin films, the ideas are equally relevant for a diverse set of problems, including: (i) unconventional superconductors like the high Tc cuprates that have a quantum critical point tuned by doping, (ii) SIT at oxide interfaces like LaAlO 3 /SrTiO 3 , (iii) SIT in the next generation of weakly coupled layered materials like dichalcogenide monolayers, and (iv) bosons in optical lattices with speckle disorder. since n i = −id/dθ i . The partition function can be expressed as the coherent-state path integral Z = D[θ]e −S with action 29 For a slowly varying phase, this becomes is an integer, then the middle term does not contribute to the free energy because β 0 ∂ τ θ i = 2π×integer. In this special case of particle-hole symmetry, the dynamical exponent is z = 1. Away from this particle-hole symmetry point, the first derivative term remains in the action, and for the pure system z = 2. Upon including disorder V i in the diagonal potential, recent Monte Carlo simulations 43 have obtained z = 1.83 ± 0.05.
The model we have studied has bond disorder J ij that respects particle-hole symmetry. In this case the dynamical exponent is expected to remain z = 1 as argued in Refs. 44 and 45. We also note that a recent Monte Carlo study of the (1+1)D JJA also concluded that z = 1 in the presence of bond disorder 46 . In order to definitively establish the value of z, two-parameter finite-size scaling with varying aspect ratios of L τ and L is necessary. Within the scope of the analysis presented here the good scaling collapse of our data for the bond-disordered model, shown in Figs. 5 and 6, is indeed consistent with z = 1.
Our Monte Carlo simulations are performed by mapping it Eq. 1 onto an anisotropic 3D classical XY model with Hamiltonian 31 by performing a Trotter decomposition of imaginary time into L τ slices of width ∆τ such that the inverse temperature β = L τ ∆τ ; r and r are points in the 2D plane and τ j denotes the j th imaginary time slice; and the dimensionless coupling constants are K τ = 1/∆τ E c and K 0 = ∆τ E J .
We perform Monte Carlo simulations using the efficient Wolff cluster update method 47 . In all of our simulations, we have set K 0 = 0.1, which we have checked to be sufficiently small to remove the error from the Trotter decomposition. For the clean system, we performed simulations on lattices of size 256×256 with L τ = 64. For the disorder tuned transition, we worked at fixed E c /E J = 3.0. Simulations at different temperatures have been performed by changing the number of imaginary time slices L τ from 32 to 128; for each L τ , we fix L = L τ since the dynamical exponent is z = 1. All disorder results have been averaged over 100 disorder realizations.

APPENDIX B: DYNAMICAL OBSERVABLES AND ANALYTIC CONTINUATION
We calculate the imaginary time, or equivalently the Matsubara frequency (ω n = 2nπ/β), current-current correlation function where the paramagnetic current in our model is given by j x (r, τ ) ≡ K 0 sin [θ(r +x, τ ) − θ(r, τ )]. The conductivity is related to the analytic continuation of Λ xx at q = 0 where −k x is the average kinetic energy along bonds in the x-direction. Re σ(ω) is then given by The superfluid stiffness ρ s is obtained from the difference between the transverse and longitudinal limits of the current-current correlation function ρ s /π = Λ xx (q x → 0, q y = 0, iω n = 0) − Λ xx (q x = 0, q y → 0, iω n = 0) and the sum rule −k x = Λ xx (q x → 0, q y = 0, iω n = 0). Finally, Re σ(ω) obeys the optical conductivity sum rule 2 ∞ 0 + dω Re σ(ω) = π −k x − ρ s , which serves as a nontrivial check on our analytic continuation results.
(1) Analytic continuation of Λ(τ ): The imaginary time correlation function Λ xx (τ ) calculated in our Monte Carlo simulations is related to its real-frequency counterpart (and subsequently to σ(ω)) through To extract the real frequency data, we have employed the maximum entropy method (MEM) 48 to invert this Laplace transform. We have performed extensive tests on our Maximum Entropy routine; further details can be found in the supplemental material of Ref. 10.
In addition to these tests on the MEM routine itself, whenever possible, we have checked those characteristics of the spectra obtained via the MEM against features which can be directly calculated for the Monte Carlo correlation functions. Gapped functions, either in σ(ω) and P (ω), have recognizable exponential decays in the imaginary time correlation functions corresponding to the gap scale, whereas spectra without a gap correspond to correlation functions with no discernible gap scale in the τ data, see Fig. 7. The extracted gap scales from Λ(τ ) or P (τ ) track consistently with those scales coming from the real frequency functions obtained after performing the analytic continuation. This is true both in the clean system and in the disordered system, where the presence of even small disorder makes the reading off a Higgs scale in the superconductor unreliable, consistent with the analytically continued results.
We have also carefully checked that the sum rule on σ(ω) is verified. For lattice systems, the optical conductivity sum rule is This includes the spectral weight contained in the deltafunction response proportional to the superfluid stiffness ρ s . The regular part of the spectrum (which we obtain from analytic continuation) satisfies We emphasize that this sum rule is not built into the MEM routine, and provides an independent verification of the procedure. The sum rule is shown in Fig. 8 for both the clean and the disorder driven transitions. which we invert using the MEM in exactly the same way as for σ(ω).
The boson spectral functions obeys the following two sum rules This second sum rule can be seen by integrating both sides Eq. 13 over τ from 0 to β. Note that, in the limit of T → 0, Eq. 14 reduces to a sum rule on Im P (ω) itself ∞ 0 dω Im P (ω) = 1. Results for the sum rules are shown in Fig. 8.

APPENDIX C: UNIVERSAL CONDUCTIVITY AND CRITICAL EXPONENTS
There have been many attempts to calculate the value of the so-called universal conductivity at the superconductor-insulator quantum phase transition. We will only focus on those models expected to be in the same universality class as our model (z = 1); a more complete history can be found in Ref. 49 and the references therein.
Since the conductivity is a universal function of ω/T at the critical point 40 , there are different and possibly distinct limiting values of σ(ω/T ) σ(0) = σ(ω → 0, T = 0) (16) σ(∞)= σ(ω = 0, T → 0) (17) In our paper, we have proposed another universal quantity that can be reliably extracted from the numerics as explained in the text. We will express all σ values in units of σ Q = 4e 2 /h. For disorder-free models in the (2+1)D XY universality class, our value of σ * ≈ 0.4 is consistent with the recent result of Ref. 32, where they found σ(0) = 0.45 ± 0.05 using Padé approximates to analytically continue MC data at the critical point, modified from the previous estimate 29 of σ(0) = 0.285 ± 0.02 obtained by extrapolation of the current-current correlation function for ω n → 0. More recently, groups 34,35 have used holographic continuation to perform analytic continuation at the critical point. They find σ(∞) = 0.32 and σ(∞) = 0.359(4) respectively.
For the disorder-tuned transition, we have obtained σ * ≈ 0.50 (19) ν = 0.96 ± 0.06 (20) where z = 1 by definition for the model. There have been only a few results on disordered transitions that can meaningfully be compared to our work. A Monte Carlo study of the (2+1)D XY model with onsite charging energy disorder 31 found z = 1.07 ± 0.03, ν ≈ 1, and σ(0) = 0.27 ± 0.04 obtained by extrapolation of Λ xx for ω n → 0. We expect that using analytic continuation could modify this estimate. Studies of the disordered quantum rotor model using strong disorder renormalization group theory 42 have found ν = 1.09 ± 0.04, although they have not looked at the universal conductivity.