Local spectroscopies across the superconductor-insulator transition

TL;DR

This study uses local spectroscopic measurements and Monte Carlo simulations to analyze quantum fluctuations across the superconductor-insulator transition, revealing spatially varying spectral and magnetic properties.

Contribution

It introduces a detailed local analysis of quantum fluctuations across the SIT using simulations and connects findings to experimental techniques.

Findings

Zero-energy spectral weight diminishes in insulating regions.

Correlations between spectral weight and compressibility are observed.

Quantum fluctuations influence local magnetic susceptibility across the transition.

Abstract

We explore the manifestation of quantum fluctuations across the superconductor-insulator transition (SIT) in three different local measurements: the two-particle density of states $P(\mathbf{r},\omega)$, compressibility $\kappa(\mathbf{r})$, and diamagnetic susceptibility $\chi(\mathbf{r})$. We perform Monte Carlo simulations on the 2D quantum Josephson junction array model and present local maps of these quantities as the system is tuned across the SIT. $P(\mathbf{r},\omega)$, obtained using Maximum Entropy analytic continuation techniques, shows strongly diminished zero-energy spectral weight in nearly-insulating islands, that also correlate with regions of suppressed $\kappa(\mathbf{r})$. We investigate the signatures of quantum fluctuations in the evolution of $\kappa(\mathbf{r})$ and $\chi(\mathbf{r})$ across the SIT, distinct from thermal fluctuations. We discuss the experimental implications of our results for scanning Josephson spectroscopy, compressibility, and scanning SQUID measurements.

Equations

$$\hat{H}=\frac{E_C}{2}\sum _i{\hat{n}}_i^2-E_J\sum _{⟨ij⟩}\cos ({\hat{\theta }}_i-{\hat{\theta }}_j)$$

$$G(\mathbf{r},{\mathbf{r}}^{\prime };\tau )=⟨e^{i\left(\hat{\theta }({\mathbf{r}}^{\prime },\tau )-\hat{\theta }(\mathbf{r},0)\right)}⟩$$

$$P(\mathbf{k},\omega )=-\frac{1}{\pi }\mathrm{Im}G(\mathbf{k},\omega ).$$

$$G(\mathbf{k},\tau )={\int }_{-\infty }^{\infty }\frac{d\omega }{\pi }\frac{e^{-\tau \omega }}{1-e^{-\beta \omega }}P(\mathbf{k},\omega ).$$

$$P(\omega )=\sum _{\mathbf{k}}P(\mathbf{k},\omega )=A(|\mathbf{r}-{\mathbf{r}}^{\prime }|=0,\omega )$$

$$⟨j_b⟩=-\frac{\partial \ln Z}{\partial A_b}=⟨K_b\sin ({\partial }_b\theta -A_b)⟩$$

$${{\rm Y}}_{b{b}^{\prime }}=\frac{\partial ⟨j_b⟩}{\partial A_{b^{\prime }}}$$

$${{\rm Y}}_{\tau \tau }(\mathbf{r},{\mathbf{r}}^{\prime };\tau ,{\tau }^{\prime })$$

$$=⟨-k_{\tau }(\mathbf{r},\tau )⟩\delta (\mathbf{r},\tau )-{\Lambda }_{\tau \tau }(\mathbf{r},\tau )$$

$$\kappa =\underset{\mathbf{k}\to 0}{\lim }{{\rm Y}}_{\tau \tau }(\mathbf{k},i{\omega }_n=0).$$

$$G(\mathbf{r},\tau )={\int }_{-\infty }^{\infty }\frac{d\omega }{\pi }\frac{e^{-\tau \omega }}{1-e^{-\beta \omega }}P(\mathbf{r},\omega ).$$

$$\kappa (\mathbf{r})$$

$$=\frac{1}{\beta }\sum _{{\mathbf{r}}^{\prime },\tau ,{\tau }^{\prime }}\left(⟨-k_{\tau }(\mathbf{r},\tau )⟩\delta ({\mathbf{r}}^{\prime },{\tau }^{\prime })-{\Lambda }_{\tau \tau }(\mathbf{r},{\mathbf{r}}^{\prime },\tau ,{\tau }^{\prime })\right).$$

$$⟨M⟩=\frac{\partial (T\log Z)}{\partial B}$$

$$\chi (\mathbf{r})=-\frac{\partial ⟨M(\mathbf{r})⟩}{\partial B}=⟨M(\mathbf{r})M⟩.$$

$$⟨M⟩=\frac{1}{2\beta }\sum _{⟨ij⟩,\tau }(x_i{y}_j-x_j{y}_i)j_{ij}^{\tau }$$

$$M(\mathbf{r})=\int d{\mathbf{r}}^{\prime }\cdot (\hat{z}\times \mathbf{j}({\mathbf{r}}^{\prime })).$$

$$Z=\mathrm{Tr}{e}^{-\beta \hat{H}}$$

$$Z$$

$$=\int \mathcal{D}\theta ⟨\theta ⟩e^{-\delta \tau \hat{H}}\dots e^{-\delta \tau \hat{H}}{e}^{-\delta \tau \hat{H}}\theta$$

$$\text{mathbbm}{1}_{\tau }=\int \left[\frac{1}{{\left(2\pi \right)}^{L^d}}\prod _{\mathbf{r}}d{\theta }_{\mathbf{r}}^{\tau }\right]|{\theta }^{\tau }⟩⟨{\theta }^{\tau }|$$

$$Z$$

$$=\int \mathcal{D}\theta \prod _{\tau =0}^{M-1}⟨{\theta }^{\tau +1}⟩e^{-\delta \tau \hat{H}}{\theta }^{\tau }\delta \left({\theta }^M-{\theta }^0\right)$$

$$Z\approx \int \mathcal{D}\theta \prod _{\tau =0}^{M-1}⟨{\theta }^{\tau +1}⟩e^{-\delta \tau \hat{T}}{e}^{-\delta \tau \hat{V}}{\theta }^{\tau }\delta \left({\theta }^M-{\theta }^0\right).$$

$$Z\approx \int \mathcal{D}\theta e^{\delta \tau E_J{\sum }_{⟨\mathbf{r}{\mathbf{r}}^{\prime }⟩,\tau }\cos \left({\theta }_{\mathbf{r}}^{\tau }-{\theta }_{{\mathbf{r}}^{\prime }}^{\tau }\right)}\prod _{\tau =0}^{M-1}{T}_{\tau }\delta \left({\theta }^M-{\theta }^0\right)$$

$$T_{\tau }$$

$$=\prod _{\mathbf{r}}\sum _{n_{\mathbf{r}}^{\tau }=-\infty }^{\infty }{e}^{-\delta \tau \frac{1}{2}{E}_C{n_{\mathbf{r}}^{\tau }}^2}{e}^{i{n}_{\mathbf{r}}^{\tau }\left({\theta }_{\mathbf{r}}^{\tau }-{\theta }_{\mathbf{r}}^{\tau +1}\right)}$$

$$S\left(\varphi \right)\stackrel{\alpha \to \infty }{\sim }{e}^{-\alpha }{e}^{\alpha \cos \varphi }$$

$$T_{\tau }\approx e^{(\delta \tau E_C{)}^{-1}{\sum }_{\mathbf{r}}\cos \left({\theta }_{\mathbf{r}}^{\tau }-{\theta }_{\mathbf{r}}^{\tau +1}\right)}.$$

$$Z\approx \int \mathcal{D}\theta e^{-S_{\mathrm{XY}}}$$

$$S_{\mathrm{XY}}=-K_{\tau }\sum _{\mathbf{r},\tau }\cos \left({\theta }_{\mathbf{r}}^{\tau }-{\theta }_{\mathbf{r}}^{\tau +1}\right)-K_{\mathbf{r}}\sum _{⟨\mathbf{r}{\mathbf{r}}^{\prime }⟩,\tau }\cos \left({\theta }_{\mathbf{r}}^{\tau }-{\theta }_{{\mathbf{r}}^{\prime }}^{\tau }\right)$$

$$K_{\tau }\equiv (\delta \tau E_C{)}^{-1}{K}_{\mathbf{r}}\equiv \delta \tau E_J.$$

$$e^{\alpha \cos \varphi }=\sum _{j=-\infty }^{\infty }{I}_j(\alpha )e^{ij\varphi }$$

Full text

I Introduction

Since its discovery, superconductivity has long been a rich environment for understanding the quantum nature of our universe. Superconductors themselves are macroscopic condensates of paired electrons that owe their coherence to quantum mechanics. It is no surprise then that superconductivity in disordered two-dimensional thin film materials has proven to be a rewarding playground for studying quantum phase transitions (QPTs).hebard1990 ; shahar1992 ; goldman1998 ; adams2004 ; sambandamurthy2004 ; steiner2005 ; stewart2007 ; gantmakher2010 These QPTs are driven by quantum fluctuations between two competing phases of matter at zero temperature tuned by a non-thermal parameter $g$.sachdev

For the superconductor-insulator transition (SIT) in thin films, theoryghosal1998 ; ghosal2001 ; bouadim2011 and experimentsacepe2008 ; sacepe2011 ; mondal2011 ; sherman2012 have both shown that the single-particle gap in the superconductor persists into the insulator. While the local amplitude remains finite and the single-particle density of states shows a robust gap, the coherence peaks become diminished and global superconductivity is lost. This suggests that Cooper pairs do not break up in the insulator, and in fact superconductivity is destroyed by the loss of global phase coherence between pairs on different islands. The SIT, then, is predominantly driven by quantum phase fluctuations.

We therefore use a bosonic description of a thin film superconductor in terms of a Josephson junction array (JJA) of superconducting islands, the Hamiltonian of which is related to the quantum XY model.wallin1994 ; sondhi1997 ; swanson2014 Here, amplitude fluctuations are ignored and we directly simulate the phase degrees of freedom using quantum Monte Carlo (QMC). This gives us access to quantum phase fluctuations, which have have been experimentally observed in global measurementscrane2007 ; sherman2015 ; poran2017 but have only recently been imaged locally using scanning SQUID techniques.kremen2018

In this paper we calculate local two-particle quantities such as compressibility, two-particle local density of states (LDOS), and local diamagnetic susceptibility in order to highlight their importance in observing local quantum fluctuations. Local spectroscopies have previously played an important role in identifying a variety of physical phenomena, including the use of scanning tunneling spectroscopy (STM) to map spatial inhomogeneities in high-$\mathrm{T_{c}}$ cuprates,howald2001 ; pan2001 ; lang2002 local compressibility measurements to observe electron-hole puddles in graphene,martin2008 and a combination of STM and spin-polarized STM to detect a possible signature of Majorana fermions in ferromagnetic chains on a superconductor.nadj-perge2014 Here we make predictions for experimental measurements that can be performed using scanning Josephson spectroscopy (SJS)randeria2016 and compressibility probes, which can be used to visualize quantum phase fluctuations.

The number-phase uncertainty principle that is at the heart of the SIT is also emerging as a paradigmatic model in quantum information. The “quantum phase slip” qubit describes a state with a well-defined flux in which number fluctuations introduce phase slips. The dual “Cooper pair box” is a qubit with well defined number of Cooper pairs in which Josephson tunneling creates number fluctuations.mooij2006 We discuss below how the qubit evolves across the SIT and connect these concepts to our results on spectroscopies.

Our main results are as follows:

(1) The global two-particle DOS $P(\omega)$ is peaked at $\omega=0$ in the superconducting phase due to the presence of a Cooper pair condensate. As the SIT is approached, this peak diminishes and the spectrum becomes gapped at the critical point, signaling a transition to a Cooper pair insulator. The global two-particle compressibility $\kappa$ is also finite in the superconducting phase where Cooper pairs are free to tunnel and vanishes at the critical point where they become localized.

(2) The local compressibility $\kappa(\mathbf{r})$ in a disordered system captures the onset of quantum fluctuations as the SIT is approached. In the superconducting phase, increasing phase fluctuations create pockets of localized Cooper pairs where the compressibility is small. The LDOS $P(\mathbf{r},\omega)$ in these regions exhibits an $\omega=0$ peak that is strongly suppressed whereas $P(\mathbf{r},\omega)$ outside of these pockets resembles the global $P(\omega)$ of a superconductor. As $g$ is increased toward the SIT, fluctuations of $\kappa(\mathbf{r})$ increase mirroring the presence of increasing phase fluctuations. Thermal fluctuations obfuscate the presence of these insulating pockets instead of increasing them in size, providing an easy way to separate quantum fluctuations from thermal fluctuations.

(3) The local diamagnetic susceptibility $\chi(\mathbf{r})$ also shows increasing fluctuations as the SIT is approached from the superconducting side. As a function of $T/T_{c}$, the standard deviation of $\chi(\mathbf{r})$ is peaked only in a narrow region around $T_{c}$ for $g$ deep in the superconducting phase. As $g$ is increased toward the critical point, this peak broadens around $T_{c}$, indicating that fluctuations of $\chi(\mathbf{r})$ are appearing well-below $T_{c}$. The fact that these extra fluctuations exist far below $T_{c}$ provide evidence that they are indeed of quantum origin.

II Model and Methods

A useful model for understanding quantum fluctuations across the SIT is the 2D JJA model with Hamiltonian

[TABLE]

where $\hat{n}_{i}$ and $\hat{\theta}_{i}$ are canonically conjugate Cooper pair number and phase operators, respectively, that satisfy the commutation relation $\left[\hat{\theta}_{i},\hat{n}_{j}\right]=i\delta_{ij}$. $E_{J}$ links phases on nearest neighbor sites via a Josephson coupling while $E_{C}$ represents the charging energy of Cooper pairs on each site. When $E_{J}$ is large, the phases align and the system is in a coherent superconducting phase. When the $E_{C}$ term dominates, the system favors a well-defined number eigenstate, leading to quantum phase fluctuations that destroy the superconducting order and transition the system to a bosonic insulating phase. Thus we use the ratio $g=E_{C}/E_{J}$ as a knob to tune the system across a QPT between a superconductor and an insulator. It is important to emphasize that loss of global phase coherence is responsible for destroying superconductivity in our model. We assume that fluctuations of the superconducting amplitude are small and that we are working at temperatures well below the pair-breaking scale $T^{\ast}$ of the superconducting island.

We simulate this model using QMC two different ways as shown schematically in Fig. 1a. In both cases we use a quantum-classical mapping to map the quantum JJA Hamiltonian to a classical action.wallin1994 ; sondhi1997 First we map the 2D JJA to a (2+1)D XY model of classical phases, a language that is well suited for calculating the two-particle DOS $P(\omega)$ and compressibility $\kappa$. Simulations of the (2+1)D XY model are performed using a Wolff cluster algorithmwolff1989 on system sizes $64\times 64$ for global calculations and $24\times 24$ for local calculations. To calculate diamagnetic susceptibility $\chi$, we instead map the JJA to a (2+1)D integer current model (ICM) because the current basis is more natural for exploring the fluctuations of diamagnetic currents. We simulate the ICM using a worm algorithmprokofev2001 on a $64\times 64$ lattice for both global and local calculations. See Appendix A for more details on these mappings.

It is important to note that simulations calculating local quantities require a small amount of disorder to be introduced in order to create structure. Without disorder, local structure is washed out by Monte Carlo averaging. We introduce a small amount of disorder in the spatial bonds of each model by randomly removing a fraction $p=0.1$ of Josephson couplings $E_{J}$ throughout the lattice. This creates regions where insulating sites can nucleate and be detected by our local probes. There have previously been studies on disorder-tuned SITs,bouadim2011 ; swanson2014 however we emphasize that the disorder in our model is static and the fluctuations we see are ultimately caused by tuning $g$, not disorder.

Both the global $P(\omega)$ and local $P(\mathbf{r},\omega)$ are calculated by analytically continuing imaginary time data to real frequencies using the Maximum Entropy Method (MEM). gubernatis1991 ; sandvik1998 The procedure is delicate and we have performed extensive tests, including checking sum rules, to ensure the validity of our results (see Appendix B for more information).

III Results

III.1 Bosonic spectral function

The bosonic Green’s function in imaginary time is given by $G(\mathbf{r},\mathbf{r}^{\prime};\tau)=\left<\hat{a}^{\dagger}(\mathbf{r}^{\prime},\tau)\hat{a}(\mathbf{r},0)\right>$ where $\hat{a}$, $\hat{a}^{\dagger}$ are Cooper pair raising and lowering operators. In the language of the JJA we can rewritefisher1989 the raising operator in terms of its amplitude and phase as $\hat{a}^{\dagger}(\mathbf{r},\tau)=\sqrt{\hat{n}(\mathbf{r},\tau)}e^{i\hat{\theta}(\mathbf{r},\tau)}$, but since we are ignoring on-site amplitude fluctuations we can write the Green’s function purely in terms of the phase variable as

[TABLE]

which is a spin-spin correlation function in the classical XY representation.

The real frequency spectral function $P(\mathbf{k},\omega)$ is the imaginary part of the corresponding real frequency Green’s function $G(\mathbf{k},\omega)$

[TABLE]

However, since we are working with imaginary time in our QMC, we need a way to analytically continue $G(\mathbf{k},\tau)$ to a real frequency spectral function. This leads to the following relation between $G(\mathbf{k},\tau)$ and $P(\mathbf{k},\omega)$

[TABLE]

Solving this equation for $P(\mathbf{k},\omega)$ amounts to inverting a Laplace transform. However, performing this procedure on QMC data of $G(\mathbf{k},\tau)$ with error bars is non-trivial and requires the use of numerical analytic continuation techniques. In our work, we use MEM to obtain $P(\mathbf{k},\omega)$ from $G(\mathbf{k},\tau)$ and have validated our results by checking relevant sum rules.

We are particularly interested in the DOS which is the sum over momentum of the spectral function

[TABLE]

This amounts to performing MEM on $G(|\mathbf{r}-\mathbf{r}^{\prime}|=0,\tau)$ data. In Fig. 2a we plot $P(\omega)/\omega$ as a function of $g$ across the SIT. Since the form of (4) requires $A(-\omega)<0$, we plot $P(\omega)/\omega$ to obtain a quantity that more closely resembles an experimentally measured DOS. We see that $P(\omega)/\omega$ is strongly peaked at $\omega=0$ in the superconducting phase, corresponding to the existence of a Cooper pair condensate. As the SIT is approached, spectral weight shifts from the central peak to the finite energy modes on either side of $\omega$ until the system forms a gap for $g>g_{c}$. We plot the shift of this spectral weight from zero energy to finite energy modes in Fig. 2d. In Fig. 2c we plot the size of this gap $\omega_{\mathrm{gap}}$ as a function of $g$ along with other energy scales including the superfluid stiffness $\rho_{s}$ and the compressibility $\kappa$ (described in the next section). As we expect, $\rho_{s}$, $\kappa$, and $\omega_{\mathrm{gap}}$ go soft at $g_{c}$ from their respective sides of the transition.

III.2 Compressibility

The global compressibility $\kappa$ is a quantity that characterizes fluctuations of the Cooper pair number density operator $\hat{n}$. We can define the average current along a generic spacetime bond $b$ in the XY representation as

[TABLE]

where $Z=\operatorname{Tr}e^{-\beta\hat{H}}$ is the partition function, $A_{b}$ is the element of an externally applied vector potential along bond $b$, and $K_{b}$ is the coupling constant along that bond. We can then identify the average density $\braket{n}$ with the current along temporal bonds $\braket{j_{\tau}}$.wallin1994

The generalized electromagnetic response tensor

[TABLE]

describes the response of a current $j_{b}$ along a spacetime bond $b$ to an externally applied vector potential $A_{b^{\prime}}$ along a bond $b^{\prime}$. While the superfluid stiffness $\rho_{s}$ can be obtained from the static, transverse long wavelength limit of the spatial response function $\Upsilon_{xx}$, the compressibility is given by the static long wavelength limit of the temporal response function $\Upsilon_{\tau\tau}$ kubo2003

[TABLE]

where $\left<-k_{\tau}(\mathbf{r},\tau)\right>=\left<K_{\tau}\cos(\partial_{\tau}\theta(\mathbf{r},\tau))\right>$ and the temporal current-current correlator $\Lambda_{\tau\tau}(\mathbf{r},\tau)=\left<K_{\tau}^{2}\sin(\partial_{\tau}\theta(\mathbf{r},\tau))\sin(\partial_{\tau}\theta(0,0))\right>$. Note that since we are performing linear response we take the limit $A_{\tau}\rightarrow 0$, and we make use of translational symmetry in the second line.

The compressibility $\kappa$ is then the static, long wavelength limit of $\Upsilon_{\tau\tau}(\mathbf{r},\tau)$

[TABLE]

Note that this amounts to calculating the response of the pair number density $\braket{n}\sim\braket{j_{\tau}}$ with respect to an externally applied electric potential $\phi\sim A_{\tau}$, which is the usual definition of compressibility.

In Fig. 2c we show $\kappa$ as a function of $g$. $\kappa$ is finite in the superconducting phase where Cooper pairs are phase coherent and are able to tunnel across the system. $\kappa$ decreases as $g_{c}$ is approached and vanishes in the Mott insulating phase where Cooper pairs become localized by phase fluctuations.

III.3 Local quantities

We next turn our attention to calculations of the LDOS $P(\mathbf{r},\omega)$ and local compressibility $\kappa(\mathbf{r})$. $P(\mathbf{r},\omega)$ is related to the local Green’s function (2) at $\mathbf{r}=\mathbf{r}^{\prime}$ by inverting (4) once again

[TABLE]

The local compressibility $\kappa(\mathbf{r})$ is given by the local response function $\Upsilon_{\tau\tau}(\mathbf{r},\tau)$ in (9)

[TABLE]

In order to extract local structure from our QMC simulations, we break translational symmetry by introducing a small fraction of bond disorder. In Fig. 3a we show a local map of $\kappa(\mathbf{r})$ at $g=3.6$, near the SIT. We see that the system forms puddles where $\kappa(\mathbf{r})$ is significantly suppressed. These incompressible regions are locations where a large density of bonds have been cut, resulting in the formation of insulating islands. In Fig. 3b we plot the corresponding $P(\mathbf{r},\omega)$ in two representative regions. We see a strong spatial correlation between the strength of $\kappa(\mathbf{r})$ and the distribution of low-energy spectral weight. The superconducting region highlighted in blue is highly compressible and has most of its spectral weight peaked strongly around $\omega=0$, reflecting the strength of the superconducting condensate. The behavior of $P(\mathbf{r},\omega)$ in this region matches that of the global $P(\omega)$ shown in purple, which reflects that of a globally phase-coherent superconductor. On the other hand, in the region highlighted in red with small compressibility, we see that the $\omega=0$ peak in $P(\mathbf{r},\omega)$ is highly suppressed . Here we can see evidence of a two-particle gap beginning to form as spectral weight shifts from $\omega=0$ to finite energy modes, indicative of an emerging Cooper pair insulator.

It is important to emphasize that the emergence of insulating islands shown in Fig. 3a is caused by quantum phase fluctuations due to proximity to a quantum critical point. To illustrate this, in Fig. 4a we also plot $\kappa(\mathbf{r})$ as a function of $g$ and temperature $T$. As $g$ increases toward the SIT, fluctuations in $\kappa(\mathbf{r})$ increase, leading to an increase in the size and prevalence of incompressible islands. Interestingly, as $T$ increases the fluctuations in $\kappa(\mathbf{r})$ become smoothed out and the insulating islands become smaller. This is due to the fact that $\kappa$ is sensitive specifically to number fluctuations. While quantum number fluctuations are expected to decrease as $g$ increases, thermal number fluctuations increase as $T$ increases due to higher energy number states becoming available. This is clearly seen in Fig. 4b, where the distribution of $\kappa(\mathbf{r})$ broadens significantly with increasing $g$, but narrows slightly with increasing $T$. This difference in behavior as $\kappa(\mathbf{r})$ evolves with $g$ and $T$ provides a way to distinguish quantum fluctuations from thermal fluctuations. We propose that an experiment that measures $\kappa(\mathbf{r})$ across the SIT will be able to directly image the presence of quantum fluctuations. The fact that these fluctuations are in fact quantum phase fluctuations is further confirmed by our results on the diamagnetic susceptibility.

III.4 Diamagnetic susceptibility

Phase fluctuations increase both as a function of temperature and a function of $g$. The question becomes how to separate thermal phase fluctuations from quantum phase fluctuations. A well-known property of a superconductor is the fact that it generates diamagnetic supercurrents in the presence of a magnetic field. In general the magnetization generated by an external field is related to the free energy by

[TABLE]

where $-T\log Z$ is the free energy in terms of the partition function $Z$ and $B$ is an external applied magnetic field. We can calculate the corresponding local diamagnetic susceptibility $\chi(\mathbf{r})$, which is sensitive to phase fluctuations, from linear response using a Kubo formula

[TABLE]

The local diagmagnetic susceptibility is the response of a local induced magnetization $M(\mathbf{r})$ to a global applied magnetic field $B$. This amounts to calculating the correlator between the local magnetization $M(\mathbf{r})$ and the global magnetization $M$. To obtain $\chi(\mathbf{r})$, we perform QMC simulations in the dual ICM representation of the quantum JJA model. This representation is more suited to calculating quantities involving the supercurrent since the QMC configurations themselves are given in terms of integer currents (see Appendix A2).

In the language of the ICM, we can write the total magnetization, assuming a uniform $B$-field in the $\hat{z}$-direction, as

[TABLE]

where $j_{ij}^{\tau}$ is an integer current on a spatial bond connecting sites $i$ and $j$ on timeslice $\tau$. The spatial pattern that $j_{ij}^{\tau}$ is summed along comes from the corresponding vector potential of the uniform $B$-field, $\mathbf{A}=\frac{B}{2}(-x,y,0)$. $M(\mathbf{r})$ is calculated from the discrete analog of the magnetization current $\mathbf{j}(\mathbf{r})=\Delta\times\mathbf{M}(\mathbf{r})$ where $\mathbf{j}$ are the integer currents. Inverting this for $M_{z}(\mathbf{r})=M(\mathbf{r})$ gives us

[TABLE]

In Fig. 5a we show local maps of $\chi(\mathbf{r})$ as a function of $g$ and $T/T_{c}$. Here, $T_{c}$ is the temperature at which the superconductor transitions to a normal state but still contains pairs; it should not be confused with the pair-breaking transition, that occurs at a much higher temperature. As expected, we observe that fluctuations of $\chi(\mathbf{r})$ increase with $T/T_{c}$ due to thermal fluctuations. Similarly, fluctuations of $\chi(\mathbf{r})$ also increase with increasing $g$. In both cases we see pockets with near-zero susceptibility beginning to form as phase fluctuations destroy the ability for coherent supercurrents to exist. However, interestingly, as $g$ increases the fluctuations in $\chi(\mathbf{r})$ appear at lower values of $T/T_{c}$. The fact that fluctuations appear well below $T_{c}$ when the system is near a quantum critical point suggests that the fluctuations with increasing $g$ are predominantly caused by quantum phase fluctuations.

This shows up as a broadening of the standard deviation of $\chi(\mathbf{r})$ across $T/T_{c}$ as shown in Fig. 5b. For small $g$, the standard deviation is peaked mainly around $T_{c}$, suggesting that only thermal fluctuations are relevant in this regime. However as $g$ increases the temperature range where standard deviation is large broadens to include temperatures well below $T_{c}$ reflecting the importance of quantum fluctuations. This broadening of the temperature range of $\chi(\mathbf{r})$ fluctuations was recently observed in scanning SQUID experiments of the thin film superconductor NbTiN.kremen2018 In the experiment, a scanning SQUID was used to directly image the local $\chi(\mathbf{r})$ and the corresponding standard deviation as a function of temperature and film thickness was found to qualitatively match that of theory. Importantly, this was the first time quantum fluctuations were directly imaged in an experiment.

IV Conclusion and Outlook

We have presented three different calculations of local two-particle response functions across the SIT: density of states $P(\mathbf{r},\omega)$, compressibility $\kappa(\mathbf{r})$, and diamagnetic susceptibility $\chi(\mathbf{r})$. We have shown how these quantities can be used to probe the local structure of fluctuations across the SIT. Particularly, for both $\chi(\mathbf{r})$ and $\kappa(\mathbf{r})$ we see an increase in local quantum fluctuations as the system approaches the critical point from the superconducting side, independent of thermal fluctuations

Our aim is to connect these calculations to spectroscopic measurements that can be performed in an experiment. We have already discussed how $\chi(\mathbf{r})$ can be measured using scanning SQUID techniques. We next turn our attention to measurements of $P(\omega)$ and $\kappa(\mathbf{r})$. It has previously been shown that it is possible capture the local structure of the superconducting order parameter using a combination of scanning tunneling spectroscopy (STM) and scanning Josephson spectroscopy (SJS) using a superconducting Pb tip.randeria2016 There, a suppression of the zero-energy peak measured in the SJS conductance was found on impurity sites, similar to our results on disorder sites of the JJA. We propose an experiment that uses SJS in conjunction with local compressibility measurementsmartin2008 to map out the evolution of quantum fluctuations across the SIT as we have done in Fig. 3 and 4.

We are also interested in connecting with recent developments in quantum information and quantum computing. As emphasized earlier, the essence of the SIT is the number-phase uncertainty principle. This can be used to define two types of dual qubits based on whether the phase dominates with quantum phase slips (QPS) disrupting that order (“QPS” qubit) on the SC site or whether the number of Cooper pairs is well-defined with Cooper pair tunneling disrupting the order (“Cooper pair box” qubit) on the insulating side.mooij2006 The behavior of the Cooper pair spectral function in Fig. 6 tracks the evolution of the qubit across the SIT .

QPS qubit: The QPS qubit is dominated by the inductive energy $E_{J}=\Phi_{0}^{2}/{(2L)}$, where $L$ is the inductance of the loop and $\Phi_{0}=h/{(2e)}$ is the SC flux quantum. The charging term mixes states with different fluxoid number $f=\Phi/\Phi_{0}$ where $\Phi$ is the flux through the loop and lifts the degeneracy at half-integer values of $f$. A current biased Josephson junction can be modeled as the dynamics of the phase in a slanted washboard potential. The phase is trapped in one of the minima yielding a zero voltage state, a superconductor, until the current exceeds a critical value.

Cooper pair box qubit: In the dual regime, the insulator has a fixed number of Cooper pairs on each island and is dominated by the charging scale $E_{C}=(2e)^{2}/{(2C)}$ where $C$ is the capacitance of the island. Josephson tunneling mixes states with $n$ and $n+1$ Cooper pairs on an island and lifts the degeneracy at half integer values. In the voltage-biased configuration, charge is trapped in a potential minimum resulting in a zero current state, an insulator, for voltages below a critical value.

Once we understand how a qubit behaves in different regimes, it in fact becomes a device to measure and quantify the degree of fluctuations, both thermal and quantum, across QPTs. We expect these ideas will motivate developments in quantum measurement.

Acknowledgements

H.K. would like to acknowledge support from the Israel US bi-national foundation grant no. 2014325. N.T. acknowledges support from the DOE-BES grant DE-FG02- 07ER46423. We would like to thank Yen Lee Loh for helpful discussions.

Appendix A: Quantum-Classical Mapping

Here we summarize the mapping from the quantum JJA Hamiltonian (1) to two different classical actions for use in Monte Carlo simulations, the XY model $S_{\mathrm{XY}}$ and the integer current model $S_{\mathrm{ICM}}$. This procedure is carried out as detailed in Wallin et. al.wallin1994

Appendix A1: Mapping to Classical XY Model

The imaginary time path integral formalism allows us to map the JJA Hamiltonian to a corresponding classical one which can be easily simulated in Monte Carlo. The process begins with taking the partition function of the quantum system

[TABLE]

and breaking the inverse temperature $\beta$ into $M$ imaginary timeslices of width $\delta\tau$, $\beta=M\delta\tau$. The partition function, written in the phase basis, then becomes

[TABLE]

where the product in $\mathcal{D}\theta$ runs over all the sites of a lattice of size $L^{d}$ and the states $\ket{\theta}=\prod_{\mathbf{r}}\ket{\theta_{\mathbf{r}}}$. We can now insert $M-1$ complete sets of states

[TABLE]

between each exponential in (20), using the imaginary timeslice index $\tau$ to keep track of each identity.

In this form the partition function becomes

[TABLE]

where we have absorbed the differential and prefactors of the identity into $\mathcal{D}\theta$.

We can now make the interpretation that the partition function (23) is a sum over paths of phase configurations $\ket{\theta}$ propagating across imaginary time, which is periodic in $\beta$. Thus far everything we have done has been exact. The form of (23) allows us to make an important simplifying approximation. If we choose a suitably large number of timeslices $M$ such that $\delta\tau$ is small, we can perform a Suzuki-Trotter decompositionsuzuki1993 on the exponential $e^{-\delta\tau\hat{H}}=e^{-\delta\tau\left(\hat{T}+\hat{V}\right)}=e^{-\delta\tau\hat{T}}e^{-\delta\tau\hat{V}}+\mathcal{O}\left(\delta\tau^{2}\right)$ and drop terms of $\mathcal{O}\left(\delta\tau^{2}\right)$

[TABLE]

Here, $\hat{T}$ and $\hat{V}$ represent the $E_{C}$ and $E_{J}$ terms in $\hat{H}$ respectively. Since $\ket{\theta}$ are eigenstates of $\hat{V}$, the second exponential pulls out the familiar classical XY cosine term

[TABLE]

where what remains is to evaluate the $\hat{T}$ term $T_{\tau}=\braket{\theta^{\tau+1}}{e^{-\delta\tau\hat{T}}}{\theta^{\tau}}$. Using the fact that the eigenstates of $\hat{T}$ are number states $\ket{n^{\tau}}=\prod_{\mathbf{r}}\ket{n_{\mathbf{r}}^{\tau}}$, we can insert a complete set of number states for each timeslice and obtain

[TABLE]

where in the last line we have used the fact that the overlap of phase and number states takes the form $\braket{\theta_{\mathbf{r}}^{\tau^{\prime}}}{n_{\mathbf{r}}^{\tau}}=e^{-in_{\mathbf{r}}^{\tau}\theta_{\mathbf{r}}^{\tau^{\prime}}}$ up to a normalization factor that can be absorbed into $\mathcal{D}\theta$.

The sum in $T_{\tau}$ now has the form $S\left(\phi\right)=\sum_{n=-\infty}^{\infty}e^{-\frac{1}{2\alpha}n^{2}}e^{in\phi}$ which can be rewritten using the Poisson summation formula in terms of the summand’s Fourier transform as $S\left(\phi\right)=\sum_{k=-\infty}^{\infty}\sqrt{2\pi\alpha}e^{-\frac{\alpha}{2}\left(\phi-2\pi k\right)^{2}}$.

Note that in this relation $\alpha=(\delta\tau E_{C})^{-1}$. In the limit of large $\alpha$ (small $\delta\tau$) this periodic sum of Gaussians serves as the Villain approximationvillain1976 of the function $e^{\alpha\cos\phi}$

[TABLE]

which allows $T_{\tau}$ to take the desired form (the constant prefactor can be absorbed into $\mathcal{D}\theta$)

[TABLE]

This leads to the fully mapped classical partition function

[TABLE]

This is the partition function for a classical anisotropic XY model in ($d+1$)-dimensions with action

[TABLE]

where the couplings are given by

[TABLE]

We simulate this model with Monte Carlo using an efficient Wolff algorithm where we tune $g=E_{C}/E_{J}$ by altering the couplings $K_{\tau}$ and $K_{\mathbf{r}}$ and tune temperature by changing the number of timeslices $M$. All global simulations are run on a $64\times 64$ lattice at temperature $T=0.15625E_{J}$ corresponding to $M=64$. All local simulations are run on a $24\times 24$ lattice at the same temperature and at bond disorder $p=0.1$.

Appendix A2: Mapping to Integer Current Model

To map to the ICM, we start with the partition function (30) using the form of $S_{\mathrm{XY}}$ given by (31). The goal is to sum over all $\theta$ configurations and rewrite $Z$ as a sum over classical integer current configurations $j$. To do this, we use use the Fourier series expansion of $e^{\alpha\cos\phi}$

[TABLE]

where $I_{j}$ is the modified Bessel function of order $j$. Using the shorthand $\partial_{b}\theta_{\mathbf{r}}^{\tau}$ to denote a phase difference along a generic spacetime bond $b\in(\hat{\tau},\hat{x},\hat{y})$, we get

[TABLE]

and inserting into (30) $Z$ becomes

[TABLE]

To obtain the ICM, we simply have to evaluate the remaining integrals over $\theta$. Note that for each $\theta$ these integrals take the form

[TABLE]

where $\mathbf{\Delta}$ is the discrete gradient applied across spacetime. This effectively puts a Kirchhoff’s law-like constraint on the currents $j_{\mathbf{r},b}^{\tau}$; the sum of the spacetime currents $j_{\mathbf{r},b}^{\tau}$ flowing into a site $(\tau,\mathbf{r})$ must equal the sum of currents flowing out.

We finally obtain a partition function written in terms of only the $j_{\mathbf{r},b}^{\tau}$ integer currents

[TABLE]

where $\sum_{[j_{\mathbf{r},b}^{\tau}]}\>^{\prime}\>$ is a set of sums over all integer currents $j_{\mathbf{r},b}^{\tau}$ in spacetime subject to the constraint $\mathbf{\Delta}\cdot\mathbf{j}_{\mathbf{r}}^{\tau}=0$ and $S_{\mathrm{ICM}}$ takes the form

[TABLE]

with corresponding couplings $K_{\mathbf{r}}$ and $K_{\mathbf{\tau}}$ given by (32). $S_{\mathrm{ICM}}$ is a model of integer currents flowing on a spacetime lattice subject to a divergenceless constraint. We perform Monte Carlo simulations on this model using a worm algorithm tuning $g$ and $T$ the same way we do in the XY model. All simulations are performed on a $64\times 64$ spatial lattice at varying temperatures. Local simulations are performed on the same lattice size with a spatial bond disorder $p=0.1$.

Appendix B: Sum Rules for Maximum Entropy Method

As shown in (4), the spectral function $P(\omega)$ is related to the Green’s function $G(\omega)$ by the inversion of a Laplace transform

[TABLE]

To check the validity of our analytic continuation procedure, we can derive several sum rules from this relation. First, note that $G(\tau=0)=1$ must be true. This leads to the first sum rule

[TABLE]

The second sum rule can be obtained by integrating both sides of (40) with respect to $\tau$ from [math] to $\beta$

[TABLE]

We plot the results of the sum rules as a function of $g$ in Fig. 7 and find our analytically continued data to be in good agreement with them.