Local spectroscopies across the superconductor-insulator transition

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

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)$

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)$

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

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.

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

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

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

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)$

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.

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

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

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.

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

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

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

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

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.

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

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

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

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

In this form the partition function becomes

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)$

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

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

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}$

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

This leads to the fully mapped classical partition function

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

where the couplings are given by

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$.

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}$

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

and inserting into (30) $Z$ becomes

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

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

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

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$.