# Anomalous superfluid density in a disordered charge density wave material:
Pd-intercalated ErTe3

Yusuke Iguchi1,2,3, Joshua A. Straquadine2,3, Chaitanya Murthy1,4, Steven A.
Kivelson1,2,4, Anisha G. Singh2,3, Ian R. Fisher1,2,3, and Kathryn A.
Moler1,2,3,4 1Stanford Institute for Materials and Energy Sciences, SLAC
National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, California
94025, USA
2Geballe Laboratory for Advanced Materials, Stanford University, Stanford,
California 94305, USA
3Department of Applied Physics, Stanford University, Stanford, California
94305, USA
4Department of Physics, Stanford University, Stanford, CA 94305, USA

###### Abstract

We image local superfluid density in single crystals of Pd-intercalated ErTe3
below the superconducting critical temperature, $T_{c}$, well below the onset
temperature, $T_{CDW}$, of (disordered) charge-density-wave order. We find no
detectable inhomogeneities. We observe a rapid increase of the superfluid
density below $T_{c}$, deviating from the behavior expected in conventional
Bardeen-Cooper-Schrieffer, and show that the temperature dependence is
qualitatively consistent with a combination of quantum and thermal phase
fluctuations.

PdxErTe3 is a model system for quasi-two-dimensional (2D) superconductivity
and for the competition between charge-density-wave (CDW) and superconducting
(SC) states. The superfluid density characterizes the phase stiffness of the
superconducting order parameter and determines the London penetration depth
$\lambda(T)$. In a conventional 3D Bardeen-Cooper-Schrieffer (BCS)
superconductor, the temperature dependence of the normalized superfluid
density, $n_{s}(T)=\lambda^{2}(0)/\lambda^{2}(T)$, is controlled by the
population of thermally excited Bogoliubov quasiparticles, and can be
calculated using the Bogoliubov–de Gennes equations [1] or the semi-classical
model [2]. At low temperatures, measurements of $n_{s}(T)$ provide information
about the superconducting gap structure $\Delta(T,\mathbf{k})$. At
temperatures close to $T_{c}$, however, the same theoretical considerations
imply that $dn_{s}(T)/dT|_{T\rightarrow T_{c}}$ is not very sensitive to the
gap structure, and changes somewhat but not dramatically in the strong-
coupling and/or dirty limits [3, 4].

$n_{s}(T)$ may have distinct features in quasi-2D conventional BCS
superconductors. When the superconducting coherence length $\xi$ is larger
than the film thickness, the Berezinskii-Kosterlitz-Thouless (BKT) theory
predicts an anomaly in the superfluid density at the BKT transition
temperature [5, 6, 7]. More generally, strong phase fluctuations may suppress
$T_{c}$ and increase $dn_{s}(T)/dT|_{T\rightarrow T_{c}}$ [8]. Such anomalies
have been observed in various ultra-thin film superconductors, including
Y1-xCaxBa2Cu3O7-δ [9], NbN [10], Pb [11], and $a$-MoGe [12].

We conducted measurements of the local diamagnetic susceptibility in PdxErTe3
($0<x<0.06$), a quasi-2D layered bulk superconductor, using scanning
superconducting quantum interference device (SQUID) microscopy (SSM) with
micron-scale spatial resolution. Our results show that the superfluid density
is homogeneous, with no detectable heterogeneity on micron scales.
Additionally, we find non-BCS-like temperature dependence of the superfluid
density with a steep slope $dn_{s}(T)/dT$ near $T_{c}$.

Figure 1:  Phase diagram of Pd-intercalated ErTe3. $T_{\mbox{\scriptsize
CDW1,2}}$ from Ref. [14]. $T_{c}$ determined by bulk ac susceptibility [14]
and STM [15]. $T_{c}$ obtained in this work (SSM) are plotted as red squares.

Recently, intertwined SC and CDW order has been observed in Pd-intercalated
ErTe3 [13, 14, 15]. The pristine ‘parent’ compound ErTe3 shows two, mutually
transverse, in-plane, unidirectional, incommensurate CDW states [16], with no
SC down to the measured lowest temperature, 100 mK [15]. Pd-intercalation
induces disorder in the crystal lattice, suppressing CDW formation and leading
to a SC ground state [Fig. 1] [14, 15]. In crystals with a Pd concentration
near $x=0.05$, long-range CDW is not observed [17]. Scanning tunneling
microscopy (STM) measurements of the tunneling conductance revealed a
homogeneous SC gap at length scales exceeding the SC coherence length, and
showed no direct correlation between the CDW and SC orders [15]. The
anisotropic in-plane coherence lengths were estimated as $\xi_{a}\sim 1500$ Å
and $\xi_{c}\sim 1000$ Å [15].

For this work, bulk single crystals of Pd-intercalated ErTe3 were grown using
the flux method [14]. We made images with a scanning SQUID susceptometer on
cleaved b-planes of PdxErTe3 at temperatures varying from 0.3 K to 3 K in a
Bluefors LD dilution refrigerator for samples with
$x=.003,.008,.018,.023,.029,.041,.054$. Our scanning SQUID susceptometer has a
pickup loop that measures the local magnetic flux $\Phi$ in units of the flux
quantum $\Phi_{0}$ [18] while scanning with a pickup loop-sample separation
$z$, which we call the height. The minimum $z$ can vary slightly between
cooldowns and is 800 nm in these measurements (supplemental material [19]).
The pickup loop is paired with a concentric field coil through which we apply
an ac current of $|I^{\text{ac}}|=$ 1 mA at a frequency of 1 kHz using an
SR830 Lock-in-Amplifier to produce a spatially varying localized ac magnetic
field [18]. We measure both quasi-static flux and the ac magnetic flux
$\Phi^{\text{ac}}$, and report the local ac susceptibility as
$\chi=\Phi^{\text{ac}}/|I^{\text{ac}}|$ in units of $\Phi_{0}$/A. SSM has been
employed to image inhomogeneous superfluid responses in unconventional
superconductors by detecting the local ac magnetic susceptibility [20, 21, 22,
23, 24, 25]. By measuring the dependence of the local susceptibility on the
scanning SQUID height, SSM enables estimation of the local London penetration
depth $\lambda$ [20, 26, 2, 28, 29, 30, 25].

Figure 2:  Homogeneous superfluid density on micron scales in PdxErTe3. (a)
Temperature dependence of local susceptibility images. (b) Normalized average
susceptibilities show sharp drops just below $T_{c}$. (c) The standard
deviation of the susceptibility shows only small peaks near $T_{c}$,
consistent with thermal drift. Inverted triangles indicate $T_{c}$ and solid
lines are numerical calculations, including thermal drifting of $\pm$5 mK
[19].

To investigate the inhomogeneity of superfluid response, we imaged the local
susceptibility at several temperatures. In all samples over the entire range
of Pd concentrations explored, we observed sharp and apparently homogeneous
transitions from the paramagnetic phase to the SC diamagnetic phase with
$T_{c}$’s in the range $T_{c}=0.8-2.8$ K [Fig. 2(a)]. The slight variation in
the observed paramagnetic (PM) susceptibility above $T_{c}$ among different Pd
concentrations could represent variation as a function of doping but could
also be due to differences in scan heights.

We analyze the susceptibility images by constructing a histogram of the number
of pixels with a given amplitude of $\chi$. The histograms show sharp peaks,
indicating a relatively homogeneous sample. The spacing between pixels is 300
nm, and each pixel samples a micron-scale area determined by the geometry of
the pickup loop and field coil. We choose a Gaussian function of the form
$\mathcal{N}\exp{(-(\chi-\beta)^{2}/2\gamma^{2})}$ to fit the peaks in the
histogram [Supplemental Fig. S1] [19]. The normalized susceptibility averaged
over the image is $<\bar{\chi}>\equiv\beta(T)/\beta$(0.5 K), and the upper
limit on the inhomogeneity of the superfluid response on micron scales is
characterized by the normalized standard deviation,
$\gamma(T)/\gamma(T>T_{c})$. Plotting $<\bar{\chi}>$ vs. $T$, we see that
$T_{c}$ as a function of Pd concentration [Fig. 2(b)] is consistent with
previous measurements based on bulk susceptibility and STM measurements [14,
15]. The upper limits on the inhomogeneity exhibit small peaks just below
$T_{c}$ [Fig. 2(c)] and are consistent with a slight thermal drift during the
scan [Supplemental Fig. S2] [19]. Thus, the superfluid response in PdxErTe3
($x$=0.003-0.054) is consistent with homogeneity on a micron scale.

Figure 3:  Local susceptibility vs. height provides $\lambda(T)$. (a) Height
dependence of local normalized susceptibility $\chi(z)/\phi_{s}$ in $x=0.041$
sample is well fitted by numerically calculated curves (solid lines) using
Supplemental Eq. (S4) with $\lambda(T)$ as a fitting parameter. The green-
color-filled area indicates the distance between the pickup loop’s center and
the sample surface when the SQUID tip touches the surface [19]. (b)
Temperature dependence of the penetration depth obtained from the fitting
results of Fig. 2(a) are plotted with an offset of 200 nm. (c) Estimated
penetration depth at $T=0$.

To determine the penetration depth, we measured susceptibility vs. height
[Fig. 3(a)]. The susceptibility is paramagnetic above $T_{c}$ and diamagnetic
below $T_{c}$. We fit the height dependence of the susceptibility [19] to a
model that assumes a circular pickup loop of radius $r^{\prime}$ and field
coil of radius $r$ at a height $z$ above the top of a thin film of thickness
$t$ on a substrate. The thin film is characterized by a London penetration
depth $\lambda$ and a paramagnetic permeability $\mu_{2}$. We estimate the
permeability $\mu_{2}=1.03\mu_{0}$, where $\mu_{0}$ is the permeability of
vacuum, by fitting the height dependence of the paramagnetic susceptibility
above $T>T_{c}$ to Supplemental Eq. (S4) with fixed parameters
$t$,$r^{\prime}$,$r$, and free parameter $\mu_{2}$. We then estimate
$\lambda(T)$ by fitting susceptibility vs. $z$ for each value of $T<T_{c}$ to
Supplemental Eq. (S4) with fixed parameters $t$,$r^{\prime}$,$r$, $\mu_{2}$, a
copper substrate permeability $\mu_{3}=\mu_{0}$, and free parameter
$\lambda(T)$.

The penetration depth does not depend strongly on temperature at low
temperatures [Fig. 3(b)]. We estimate $\lambda(T=0)$ across the doping series
to be in the range of 700-1000 nm, consistent with measurements of an isolated
vortex field [Supplemental Fig. S3] [19]. This penetration depth is a factor
of $3.5-5$ larger than the only other estimate of $\lambda$ in this material
of which we are aware, which was an indirect estimate from the lower critical
magnetic field at $T/T_{c}\sim 0.7$ for an $x=0.043$ sample [15]. The error
bars shown in the figure include all sources of error of which we are aware
(supplemental materials [19]). Interestingly, we did not observe a significant
dependence of $\lambda(T=0)$ on the Pd-intercalation concentration [Fig.
3(c)]. In BCS theory, $\lambda^{2}$ would be expected to decrease in
proportion to the mean free path [31], so the flat dependence of $\lambda(0)$
on $x$ suggests either that BCS theory does not apply or that $x$ is not the
main determining factor for the mean free path.

Using the obtained values of $\lambda$, we calculate the normalized superfluid
density $n_{s}(T)=\lambda^{2}(0)/\lambda^{2}(T)$. Our results reveal a rapid
increase of $n_{s}$ with decreasing temperature just below $T_{c}$ and a
slower increase at lower temperatures [Fig. 4(a)]. This temperature dependence
clearly deviates from the expectations of the conventional weak coupling
s-wave model (BCS model).

To investigate whether the anomalous temperature dependence of $n_{s}$ can be
simply attributed to details of the gap structure or strong coupling effects,
we consider an anisotropic s-wave model. In this model, the superconducting
gap is described as $\Delta(T,{\mathbf{k}})=\Delta_{0}(T)\times
g({\mathbf{k}})$, where $\Delta_{0}(T)$ represents the temperature dependence
of the gap, and $g({\mathbf{k}})$ its angular variation on the Fermi surface
[4]. The temperature dependence is approximated by the typical mean-field form
$\Delta_{0}(T)=\Delta_{0}(0)\tanh(\pi
T_{c}\sqrt{\alpha(T_{c}/T-1)}/\Delta_{0}(0))$, where $\Delta_{0}(0)$ is the
gap magnitude at $T=0$ and $\alpha$ is a parameter. For a gap with anisotropic
s-wave symmetry on a 2D cylindrical Fermi surface,
$g(\phi)=\sqrt{1-\varepsilon\sin^{2}\phi}$, where $\phi=0$ and $\pi/2$
correspond to the a and c axes, respectively, and
$\varepsilon=1-\left[\Delta_{c}(0)/\Delta_{a}(0)\right]^{2}$ (assuming that
$0<\Delta_{c}\leq\Delta_{a}$). We note that our model does not determine which
axis has a larger gap amplitude, as we take an angular average for the
normalized superfluid density. The fitting parameters in this model are
$\Delta_{0}(0)$, $\varepsilon$ and $\alpha$ [19], and the normalized
superfluid density is:

$n_{i}(T)=1-\frac{1}{2\pi T}\int_{0}^{2\pi}d\phi\,P_{i}(\phi)\\\
\int_{0}^{\infty}d\epsilon\,\cosh^{-2}\\!\left(\frac{\sqrt{\epsilon^{2}+\Delta^{2}(T,\phi)}}{2T}\right),$
(1)

where $i=a,c$, and $P_{a}=\cos^{2}\phi$, $P_{c}=\sin^{2}\phi$. We find that
our measured normalized superfluid density $n_{s}\simeq(n_{a}+n_{c})/2$ can
indeed be well fitted using Eq. (1) (for details of the fits, see the
supplemental material [19]) [Fig. 4]. However, the fitted parameter
$\alpha\sim 10$ is much larger than known models, such as $\alpha=1$
(isotropic s-wave) and $\alpha=2$ (s+g-wave) [4]. Moreover, the quality of the
fit strongly depends on the value of $\alpha$ rather than the anisotropy
$\varepsilon$ or the coupling constant $\Delta_{0}(0)/k_{B}T_{c}$. Thus our
fitting results suggest that the temperature-dependent superfluid density
cannot fit the BCS model.

Figure 4:  Comparison of the estimated normalized superfluid density from Fig.
3(b) to an anisotropic s-wave BCS model. (a) Superfluid density (dots) and
fits (solid lines), offset by 0.5. (b) Fitted values of $\alpha$ vs $x$.
Values $\alpha>>1$ are physically unrealistic for known BCS models. (c) Fitted
coupling constant $\Delta_{0}(0)/k_{B}T_{c}$ vs $x$. Figure 5:  Normalized
superfluid density in the quantum rotor model compared to experimental results
in PdxErTe3. (a) Variational solution of the quantum rotor model for several
values of the coupling constant $CJ$. (b) The experimental results from Fig.
4(a) compared with the results of the quantum rotor model from Fig. 5(a) and
Supplemental Fig. S5.

We next consider fluctuations, which can suppress $T_{c}$ and modify the
temperature dependence of the superfluid density. Quasi-2D electronic
structures can enhance fluctuations [15, 16]. A pure BKT scenario cannot be
applied here, as the sample thicknesses exceed the coherence length. Classical
phase fluctuations alone would destroy the SC order above $T_{\theta}=$ 7-14
K, estimated from formulas in Ref. [8] using $\xi$ = 100-150 nm and $\lambda$
= 700-1000 nm. Notably, this estimated $T_{\theta}$ is close to $T_{c}$,
suggesting that such phase fluctuations might significantly contribute to the
determination of $T_{c}$. (Note that Fang et al. estimated $T_{\theta}$ as 170
K, much larger than $T_{c}$, from $\lambda$ = 200 nm [15].) However,
superfluid density that is dominated by classical phase fluctuations would
exhibit a linear-$T$ dependence well below $T_{c}$ [33], not flattening until
quantum effects become important. Therefore, classical phase fluctuations
alone cannot explain our results.

Quantum phase fluctuations may modify this scenario. The small value of
$T_{\theta}$ and the quasi-2D character of the electronic structure likely
enhance the effectiveness of these fluctuations, which may be further enhanced
[34] by a degree of randomness of the inter-layer Josephson coupling produced
by the Pd intercalation. To determine whether a combination of quantum and
classical phase fluctuations might account for the observed anomalous
$T$-dependence of the superfluid density, we have studied a caricature of the
problem in terms of the quantum rotor model, governed by the Hamiltonian

$H=\sum_{j}{\frac{n_{j}^{2}}{2C}}-J\sum_{<i,j>}{\cos{(\theta_{i}-\theta_{j})}},$
(2)

where $n_{j}$ is the number of Cooper pairs on site $j$ and satisfies the
commutation relations
$\left[n_{i},n_{j}\right]=\left[e^{i\theta_{i}},e^{i\theta_{j}}\right]=0$ and
$\left[n_{i},e^{i\theta_{j}}\right]=\delta_{ij}\,e^{i\theta_{j}}$, $C$ is a
local capacitance which plays the role of an effective mass, and $J$ is a
measure of the phase stiffness within each plane. (This model omits many
possibly significant effects, including long-range Coulomb interactions and
dissipation stemming from the existence of quasiparticle excitations.) For
this model, we estimate the $T$-dependent superfluid density using the
variational method used in [35] (for details, see the supplemental material
[19]). The results for a range of coupling constants $C$ and $J$ capture some
of the salient features of our experimental findings, as shown in Fig. 5(b),
suggesting that strong quantum phase fluctuation are probably significant.

Finally, it is worth noting that $T_{c}$ displays a complex variation with $x$
as shown in Fig. 1, where $T_{c}$ initially rises rapidly with $x$ before
approximately ‘saturating’. The fact that $T_{c}$ does not decrease with the
disorder at $x>0.02$ might be attributed to Anderson’s theorem, but this
theorem does not explain the initial rise relative to zero Pd concentration
[36]. The $x$-dependence of $T_{c}$ likely reflects the complex interplay of a
variety of factors, including the competition between CDW formation and
superconductivity, the effects of disorder on the CDW state, and also the
influence of quantum phase fluctuations on the superconducting state.

In summary, we have used scanning SQUID susceptometry to examine, at the
microscopic level, the superfluid response on cleaved surfaces of Pd-
intercalated ErTe3. Our findings reveal that the superfluid response is
uniformly on a micron scale within the Pd-intercalation-induced
superconducting state, consistent with previous STM measurements. We also
observe an unexpectedly strong (relative to BCS) temperature dependence of the
superfluid density near $T_{c}$ for all Pd concentrations. To explain this
non-BCS-like temperature-dependent superfluid density in PdxErTe3, we employ
the quantum rotor model. Our results suggest that quantum phase fluctuations
suppress $T_{c}$ and determine the functional form of $\lambda(T)$ in
PdxErTe3. Moreover, our study highlights the potential of temperature-
dependent superfluid density as a valuable tool for investigating quantum
phase fluctuations in quasi-2D superconductors.

###### Acknowledgements.

The authors thank A. Kapitulnik, A. Fang, and J. R. Kirtley for fruitful
discussion. This work was primarily supported by the Department of Energy,
Office of Science, Basic Energy Sciences, Materials Sciences and Engineering
Division, under Contract No. DE- AC02-76SF00515. Y.I. was partially supported
by a JSPS Overseas Research Fellowship. C.M. was supported in part by the
Gordon and Betty Moore Foundation’s EPiQS Initiative through GBMF8686.

## References

  * [1] T. Kita, Statistical Mechanics of Superconductivity (Springer, Tokyo, 2015).
  * [2] B. S. Chandrasekhar and D. Einzel, Ann. Phys. 2(new volume 505), 535 (1993).
  * [3] D. Xu, S.K. Yip, and J.A. Sauls, Phys. Rev. B 51, 16233 (1995).
  * [4] A. Maisuradze, R. Gumeniuk, W. Schnelle, M. Nicklas, C. Baines, R. Khasanov, A. Amato, and A. Leithe-Jasper, Phys. Rev. B 86, 174513 (2012).
  * [5] L. Benfatto, C. Castellani, and T. Giamarchi, Phys. Rev. B 77, 100506(R) (2008).
  * [6] V. L. Berezinskii, J. Exp. Theor. Phys. 34 610 (1972).
  * [7] J. M. Kosterlitz and D. J. Thouless, J. Phys. C: Solid State Phys. 5 L124 (1972).
  * [8] V. J. Emery and S. A. Kivelson, Nature 374, 434 (1995).
  * [9] I. Hetel, T. R. Lemberger, R. Mohit, Nat. Phys. 3, 700(2007).
  * [10] A. Kamlapure, M. Mondal1, M. Chand, A. Mishra, J. Jesudasan, V. Bagwe, L. Benfatto, V. Tripathi, and P. Raychaudhuri, Appl. Phys. Lett. 96, 072509 (2010).
  * [11] H. Nam, H. Chen, T. Liu, J. Kim, C. Zhang, J. Yong, T. R. Lemberger, P. A. Kratz, J. R. Kirtley, K. A. Moler, P. W. Adams, A. H. MacDonald, and C.-K. Shih, Proc. Natl Acad. Soc. 113, 10513(2016).
  * [12] S. Mandal, S. Dutta, S. Basistha, I. Roy, J. Jesudasan, V. Bagwe, L. Benfatto, A. Thamizhavel, and P. Raychaudhuri, Phys. Rev. B, 102, 060501(R) (2020).
  * [13] J. B. He, P. P. Wang, H. X. Yang, Y. J. Long, L. X. Zhao, C. Ma, M. Yang, D. M. Wang, X. C. Shangguan, M. Q. Xue, P. Zhang, Z. A. Ren, J. Q. Li, W. M. Liu and G. F. Chen, Supercond. Sci. Technol. 29, 065018(2016).
  * [14] J. A. W. Straquadine, F. Weber, S. Rosenkranz, A. H. Said, and I. R. Fisher, Phys. Rev. B 99, 235138(2019).
  * [15] A. Fang, A. G. Singh, J. A. W. Straquadine, I. R. Fisher, S. A. Kivelson, and A. Kapitulnik, Phys. Rev. Res. 2, 043221 (2020).
  * [16] N. Ru, C. L. Condron, G. Y. Margulis, K. Y. Shin, J. Laverock, S. B. Dugdale, M. F. Toney, and I. R. Fisher, Phys. Rev. B 77, 035114 (2008).
  * [17] A. Fang, J. A. W. Straquadine, I. R. Fisher, S. A. Kivelson, and A. Kapitulnik, Phys. Rev. B 100, 235446(2019).
  * [18] J. R. Kirtley, L. Paulius, A. J. Rosenberg, J. C. Palmstrom, C. M. Holland, E. M. Spanton, D. Schiessl, C. L. Jermain, J. Gibbons, Y.-K.-K. Fung, M. E. Huber, D. C. Ralph, M. B. Ketchen, G. W. Gibson Jr., and K. A. Moler, Rev. Sci. Instrum. 87, 093702 (2016).
  * [19] See Supplemental Material at [URL will be inserted by publisher] for (Sec.1) the details of calculations of inhomogeneity of superfluid response, (Sec.2) the details of isolated vortex field measurements,(Sec. 3) the estimate of London penetration depth from the susceptibility measurements, (Sec.4) the details of fitting the temperature-dependent superfluid density, (Sec.5) and the detailed calculations of superfluid density in the quantum rotor model.
  * [20] C. W. Hicks, T. M. Lippman, M. E. Huber, J. G. Analytis, J.-H. Chu, A. S. Erickson, I. R. Fisher, and K. A. Moler, Phys. Rev. Lett. 103, 127003 (2009).
  * [21] B. Kalisky, J. R. Kirtley, J. G. Analytis, J.-H. Chu, A. Vailionis, I. R. Fisher, and K. A. Moler, Phys. Rev. B 81, 184513 (2010).
  * [22] C. A. Watson, A. S. Gibbs, A. P. Mackenzie, C. W. Hicks, and K. A. Moler, Phys. Rev. B 98, 094521 (2018).
  * [23] I. P. Zhang, J. C. Palmstrom, H. Noad, L. Bishop-Van Horn, Y. Iguchi, Z. Cui, E. Mueller, J. R. Kirtley, I. R. Fisher, and K. A. Moler, Phys. Rev. B 100, 024514 (2019).
  * [24] L. B.-V. Horn, Z. Cui, J. R. Kirtley, and K. A. Moler, Rev. Sci. Instrum. 90, 063705 (2019).
  * [25] Y. Iguchi, H. Man, S. M. Thomas, F. Ronning, P. F. S. Rosa, and K. A. Moler, Phys. Rev. Lett. 130, 196003 (2023).
  * [26] L. Luan, T. M. Lippman, C. W. Hicks, J. A. Bert, O. M. Auslaender, J.-H. Chu, J. G. Analytis, I. R. Fisher, and K. A. Moler, Phys. Rev. Lett. 106, 067001 (2011).
  * [27] J. R. Kirtley, B. Kalisky, J. A. Bert, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, J. H. Ngai, Y. Segal, F. J. Walker, C. H. Ahn, and K. A. Moler, Phys. Rev. B 85, 224518 (2012).
  * [28] T. M. Lippman, B. Kalisky, H. Kim, M. A. Tanatar, S. L. Bud’ko, P. C. Canfield, R. Prozorov, and K. A. Moler, Physica C 483, 14 (2012).
  * [29] J. A. Bert, K. C. Nowack, B. Kalisky, H. Noad, J. R. Kirtley, C. Bell, H. K. Sato, M. Hosoda, Y. Hikita, H. Y. Hwang, and K. A. Moler, Phys. Rev. B 86, 060503(R) (2012).
  * [30] Y. Iguchi, I. P. Zhang, E. D. Bauer, F. Ronning, J. R. Kirtley, and K. A. Moler, Phys. Rev. B 103, L220503 (2021).
  * [31] M. Tinkham, Introduction to superconductivity, 2nd Ed. (Dover Publications, 2004).
  * [32] R. Prozorov and R. W. Giannetta, Supercond. Sci. Technol. 19, R41 (2006).
  * [33] E. W. Carlson, S.A. Kivelson, V. J. Emery, and E. Manousakis, Phys. Rev. Lett. 83, 612 (1999).
  * [34] E. Nakhmedov, O. Alekperov, and R. Oppermann, Phys. Rev. B 86, 214513 (2012).
  * [35] S. Chakravarty, G.-L. Ingold, S. Kivelson, and A. Luther, Phys. Rev. Lett. 56, 2303 (1986).
  * [36] P.W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).

Supplemental Material for
“Anomalous superfluid density in a disordered charge density wave material:
Pd-intercalated ErTe3 ”
by Iguchi $et$ $al.$

##  I Inhomogeneity of superfluid response

Figure S1:  Number of pixels showing the same susceptibility at $T$ = 1.0-2.6
K in x = 0.029. Solid lines are fitting curves of
$\mathcal{N}\exp(-(x-\beta)^{2}/2\gamma^{2})$.

We estimated the noise in the inhomogeneity
$\gamma(T)/\gamma(>T_{\mbox{\scriptsize c}})$ due to the experimental noise of
susceptibility measurement and the thermal drift effect. The experimental
noise of susceptibility, $\chi_{N}=\pm$ 0.025 $\Phi_{0}$/A, does not depend on
the temperature, thus this noise is considered as $C_{\mbox{\scriptsize
ex}}=1$. If we assumed the drifting temperature range was $\Delta T_{N}=\pm$5
mK, the noise of the thermal drift effect was estimated as,

$C_{\mbox{\scriptsize th}}=(\Delta\chi(T_{m+1})+\Delta\chi(T_{m}))\times\Delta
T_{N},$ (S1)

where

$\Delta\chi(T_{m})=\frac{\chi(T_{m+1})-\chi(T_{m})}{T_{m+1}-T_{m}}\frac{1}{\bar{\chi}_{N}},$
(S2) $\bar{\chi}_{N}=\frac{\chi_{N}}{\beta(0.5\mathrm{K})}.$ (S3)

The estimated total noises of $C_{\mbox{\scriptsize tot}}=C_{\mbox{\scriptsize
ex}}+C_{\mbox{\scriptsize th}}$ were quantitatively consistent with the
observed data in Fig. 2(c).

Figure S2:  (a) Normalized averaged susceptibility with experimental noise of
$\pm 0.025~{}\Phi_{0}$/A. (b) Normalized averaged susceptibility with
experimental noise of $\pm 0.025~{}\Phi_{0}$/A and thermal drift of $\pm 5$
mK.

##  II Isolated vortex field

We observed an isolated vortex at $x=0.29$ at 0.5 K, which was shown in Fig.
A3. The cross section of the magnetic flux was consistent with the numerical
simulation of a point source magnetic monopole field with total flux
$\Phi=\Phi_{0}$, which includes the SQUID structure[1]. In this model, the
magnetic monopole is set at $z=0$, where the sample surface is at $z=\lambda$,
and and the center of pickup loop is at $z=z_{0}+\lambda$, where $z_{0}=800$
nm and the magnetic monopole field is defined as
$H(\vec{r})=\Phi_{0}z/\mu_{0}r^{3}$. The best fit curve used $\lambda$ = 1.0
$\mu$m ($\pm$0.2 $\mu$m). This is consistent with the result from the
susceptibility height dependence(Fig. 3(c)).

Figure S3:  Cross section of an isolated vortex field observed was fitted by
the point source model. In the inset, the dashed line shows the cross section
and scale bar is 5 $\mu$m.

##  III Estimate the London penetration from Susceptibility measurements

To estimate the local $\lambda$, we fit the observed height dependence of
$\chi$ to an expression [Eq. (S4)] developed by Kirtley et al. for the case of
homogeneous and isotropic $\lambda$ [2]. This approach is valid if $\lambda$
varies slowly on the relevant length scales and is approximately isotropic in-
plane (the out-of-plane $\lambda$ does not appear in the result for $\chi$ in
this case [3]). In the model, the sample surface is the $z=0$ plane. The
pickup loop and field coil are at $z>0$ in vacuum, where the permeability is
$\mu_{0}$. In the sample ($0\geq z\geq-t$, where $t$ is the sample thickness),
the London penetration depth is $\lambda$ and the permeability is $\mu_{2}$.
Below the sample ($-t>z$), there is a non-superconducting substrate with a
permeability $\mu_{3}$. The radius of the field coil and the pickup loop are
$r$ and $r^{\prime}$, respectively. By solving Maxwell’s equations and the
London equation for the three regions in the limit of $r^{\prime}\ll r$
($r=0.79~{}\mu$m and $r^{\prime}\sim 0.1~{}\mu$m), one obtains the SQUID
height dependence of the susceptibility $\chi(z)$ as

$\displaystyle\chi(z)/\phi_{s}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}dx\,e^{-2x\bar{z}}xJ_{1}(x)$ (S4)
$\displaystyle\left[\frac{-(\bar{q}+\bar{\mu}_{2}x)(\bar{\mu}_{3}\bar{q}-\bar{\mu}_{2}x)+e^{2\bar{q}\bar{t}}(\bar{q}-\bar{\mu}_{2}x)(\bar{\mu}_{3}\bar{q}+\bar{\mu}_{2}x)}{-(\bar{q}-\bar{\mu}_{2}x)(\bar{\mu}_{3}\bar{q}-\bar{\mu}_{2}x)+e^{2\bar{q}\bar{t}}(\bar{q}+\bar{\mu}_{2}x)(\bar{\mu}_{3}\bar{q}+\bar{\mu}_{2}x)}\right],$

where $\phi_{s}=A\mu_{0}/2\Phi_{0}r$ is the self inductance between the field
coil and pickup loop, $A$ is the effective area of the pickup loop,
$\bar{z}=z/r$, $J_{1}$ is the Bessel function of first order, $\bar{t}=t/r$,
$\bar{q}=\sqrt{x^{2}+\bar{\mu}_{2}(r/\lambda)^{2}}$, and
$\bar{\mu}_{2}=\mu_{2}/\mu_{0}$. For the bulk sample on a copper substrate
($\bar{t}\gg 1,\bar{\mu}_{3}=1$), the observed height dependence only depends
on $\lambda$, $\mu_{2}$ and the SQUID structure. The observed susceptibility,
which is normalized by the mutual inductance 55 $\Phi_{0}$/A at far from the
sample surface, at different heights allowed us to estimate the local
$\lambda$ by fitting to Eq. (S4) with $\bar{t}=10$, $\bar{\mu}_{2}=1.03$, and
$\bar{\mu}_{3}=1$. The permeability $\bar{\mu}_{2}$ is estimated by fitting
the height dependence of the paramagnetic susceptibility at $T>T_{c}$. The
scan height $z$ is determined as $z=z_{cal}(V-V_{0})$, where
$z_{cal}=0.5~{}\mu$m/V is estimated by fitting the height dependence of the
paramagnetic susceptibility, $V_{0}=400$ nm is the distance between the pickup
loop’s center and the sample surface when the SQUID tip touches the surface
($V=0$), determined by optical measurements at room temperature, $V=0$ is
determined by detecting a kink in capacitance measurements. The error value of
$z$ causes the biggest error for determining $\lambda(0)$. The errors in
determining $V=0$ and $V_{c}$ are $\sim\pm$50 nm and $\sim\pm$100 nm,
respectively. Error-bars in Figs. 3(c) and 4(a) are calculated by using
$\pm$100 nm uncertainty in $z$. Note that we used $V_{0}=500$ nm for
$x=$0.018, 0.023, and 0.029 because the SQUID tip inadvertently picked up a
small particle while measuring the $x=0.018$ sample, prior to measuring the
other two. An estimate of 100 nm was determined by comparing the
susceptibility amplitude of samples with particles to those without (data not
shown).

##  IV Fitting the temperature-dependent superfluid density

We consider the 2D ellipsoid form of the superconducting gap as the
anisotropic $s$-wave pairing symmetry. 2D ellipsoid form is expressed by
$(x^{2}/a^{2})+(y^{2}/b^{2})=1$ where $x=a\cos{\phi}$ and $y=b\sin{\phi}$. In
the polar coordinate, this form is expressed by
$r^{2}=x^{2}+y^{2}=a^{2}\left[1-\left(1-(b/a)^{2}\right)\sin^{2}{\phi}\right]$.
Thus we introduce the anisotropic $s$-wave symmetry model as
$\Delta(T,{\mathbb{k}})=\Delta_{a}(T)\times\sqrt{1-[1-(\Delta_{c}(0)/\Delta_{a}(0))^{2}]\sin^{2}{\phi}}=\Delta_{0}(T)\times
g(\phi,\varepsilon)$, where $\varepsilon=1-(\Delta_{c}(0)/\Delta_{a}(0))^{2}$.

To fit our experimental results of the temperature-dependent superfluid
density, for simple estimate, we use the approximate formula of

$\Delta^{\prime}_{0}(T)=\Delta_{0}(T)/k_{\rm B}T_{\rm
c}=\Delta^{\prime}_{0}(0)\tanh(\pi\sqrt{a(T_{\rm
c}-T)/T}/\Delta^{\prime}_{0}(0))$ (S5)

with fitting parameters $a$ and $\Delta^{\prime}_{0}$ [5, 4].

The temperature-dependent superfluid density is expressed by the semi-
classical model as Eq. (2) in the main paper. To fit the experimentally
obtained superfluid density $n^{exp}_{s}(T)$ by using the simulation model
$n^{sim}_{s}(T,a,\varepsilon,\Delta_{a}(0),\lambda(0))$, we calculate the
error square sum
$\Xi=\sum_{i}[n_{s}(T_{i})-n^{sim}_{s}(T_{i},a,\varepsilon,\Delta_{a}(0),\lambda(0))]^{2}$.
Figures S4 estimate for the uncertainties in the parameters
$a,\Delta_{c}(0)/\Delta_{a}(0),\Delta_{a}(0)/k_{\rm B}T_{\rm c},$ and
$\lambda(0)$ [2]. The errors in parameters are determined by the double amount
of the global minimum $\Xi_{\rm min}$. The errors in $\lambda(0)$ are
estimated as less than $\pm 40$ nm, which is much shorter than the systematic
errors due to uncertainty in determining the scan height, $\pm 100$ nm, that
is used in Fig. 3(c).

Figure S4:  The error square sum $\Xi$ normalized by the global minimum
$\Xi_{\rm min}$ are plotted for (a-g) fitting parameter $a$ and the coupling
amplitude $\Delta_{a}(0)/k_{B}T_{c}$ or (h-n) anisotropy
$\Delta_{c}(0)/\Delta_{a}(0)$ and $\lambda(0)$ in Eq. (S5) for each
experimental results. Red dots are the global minimum points. Error bars in
Figs. 4(b,c) are determined by the range of $\Xi_{\rm min}<2$. Figure S5:  The
best fitting results for all measured superfluid density. (a-g) Temperature-
dependent superconducting gap with best fitting parameters of
$a,\Delta_{a}(0)/k_{B}T_{c},\Delta_{c}(0)/\Delta_{a}(0)$.

##  V Computing the superfluid density from quantum phase fluctuations

###  V.1 The problem

We imagine that it is only quantum and thermal phase fluctuations that
determine the evolution of the superfluid density, i.e. local pairing is
assumed to exist over a broader range of $T$. The simplest such model is the
quantum rotor model, which we will adopt for simplicity—although more
realistic models can be treated in similar fashion.

We thus consider the Hamiltonian

$H=\sum_{j}\frac{n_{j}^{2}}{2C}-J\sum_{\langle
i,j\rangle}\cos(\theta_{i}-\theta_{j}),$ (S6)

where $n_{j}$ is the number of Cooper pairs on site $j$ and satisfies the
commutation relations:

$\left[n_{i},n_{j}\right]=\left[e^{i\theta_{i}},e^{i\theta_{j}}\right]=0$ (S7)

and

$\left[n_{i},e^{i\theta_{j}}\right]=\delta_{ij}\,e^{i\theta_{j}}.$ (S8)

This can be expressed in terms of the imaginary time effective action

$S=\int_{0}^{\beta}d\tau\left\\{\frac{C}{2}\sum_{j}\dot{\theta}_{j}^{2}-J\sum_{\langle
i,j\rangle}\cos(\theta_{i}-\theta_{j})\right\\}.$ (S9)

###  V.2 Variational solution

We treat the problem using the same variational method as used (for a related
problem) in Ref. 6. We introduce a variational effective action,

$S_{\text{tr}}=\int_{0}^{\beta}d\tau\left\\{\frac{C}{2}\sum_{j}\dot{\theta}_{j}^{2}+\frac{\mu}{2}\sum_{\langle
i,j\rangle}(\theta_{i}-\theta_{j})^{2}\right\\},$ (S10)

and choose $\mu$ to minimize the variational free energy. The trial action is
quadratic and can be diagonalized by inserting Fourier expansions of the
fields (here $\omega_{n}=2\pi n/\beta$ is a bosonic Matsubara frequency):

$\theta_{j}(\tau)=\frac{1}{\beta}\sum_{n}\int\frac{d^{d}k}{(2\pi)^{d}}\
e^{i\omega_{n}\tau-i\vec{k}\cdot\vec{r}_{j}}\,\theta_{\vec{k},n}.$ (S11)

Since $\theta_{j}(\tau)$ is real,
$\theta_{\vec{k},n}=(\theta_{-\vec{k},-n})^{*}$.

(Note that here we have completely ignored the fact that $\theta$ is an
angular variable, i.e. we have neglected vortices. We should in principle
account for this fact by including linear terms of the form $2\pi m\tau/\beta$
in the Fourier expansion of $\theta(\tau)$, where $m$ is an integer winding
number. However, this would complicate the analysis, and presumably only alter
the results significantly when $\mu$ is very small—perhaps small enough to be
pre-empted by a first-order transition to the disordered state. This justifies
the present “spin-wave” approximation in which we ignore the periodicity of
$\theta$.)

Inserting the Fourier expansion of $\theta_{j}(\tau)$, the variational
effective action becomes

$S_{\text{tr}}=\frac{C}{2\beta}\sum_{n}\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left(\omega_{n}^{2}+\Omega_{\vec{k}}^{2}\right)\left|\theta_{\vec{k},n}\right|^{2},$
(S12)

with

$\Omega_{\vec{k}}^{2}=\frac{\mu}{2C}\sum_{j=\text{n.n.}\,i}\left|1-e^{i\vec{k}\cdot\vec{R}_{ij}}\right|^{2}=\frac{\mu}{C}\sum_{a=1}^{d}\sin^{2}(k_{a}/2),$
(S13)

where the second equality is for the $d$-dimensional hypercubic lattice.

Since the trial action is quadratic and diagonal in $\theta_{\vec{k},n}$, one
has

$\langle\theta_{\vec{k},n}\theta^{*}_{\vec{k}^{\prime},m}\rangle_{\text{tr}}=(2\pi)^{d}\delta(\vec{k}-\vec{k}^{\prime})\,\delta_{n,m}\
\frac{\beta}{C(\omega_{n}^{2}+\Omega_{\vec{k}}^{2})}.$ (S14)

It follows that

$\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}=\frac{1}{\beta
C}\sum_{n}\int\\!\frac{d^{d}k}{(2\pi)^{d}}\
\frac{1}{\omega_{n}^{2}+\Omega_{\vec{k}}^{2}}\left|1-e^{i\vec{k}\cdot\vec{R}_{ij}}\right|^{2}.$
(S15)

Performing the Matsubara sum over $n$ yields

$\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}=\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left[\frac{2f(\Omega_{\vec{k}})+1}{2C\Omega_{\vec{k}}}\right]\left|1-e^{i\vec{k}\cdot\vec{R}_{ij}}\right|^{2},$
(S16)

where $f(\omega)=(e^{\beta\omega}-1)^{-1}$ is the Bose occupation factor.
Moreover, because the trial action is quadratic,

$\langle\exp[i(\theta_{i}-\theta_{j})]\rangle_{\text{tr}}=\exp\\!\left[-\frac{1}{2}\langle[\theta_{i}-\theta_{j}]^{2}\rangle_{\text{tr}}\right].$
(S17)

We now use the variational inequality

$F\leq F_{\text{tr}}+\frac{1}{\beta}\langle
S-S_{\text{tr}}\rangle_{\text{tr}}.$ (S18)

The right-hand side is (using the formulae above)

$F_{\text{tr}}-\sum_{\langle
i,j\rangle}\left\\{Je^{-\frac{1}{2}\langle[\theta_{i}-\theta_{j}]^{2}\rangle_{\text{tr}}}+\frac{\mu}{2}\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}\right\\}.$
(S19)

The derivative of this expression with respect to $\mu$ should vanish. Using
the fact that $\partial F_{\text{tr}}/\partial\mu=\frac{1}{2}\sum_{\langle
i,j\rangle}\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}$, we obtain
the self-consistency condition for $\mu$:

$\mu=J\exp\left[-\frac{1}{2}\langle[\theta_{i}-\theta_{j}]^{2}\rangle_{\text{tr}}\right].$
(S20)

The trial action $S_{\text{tr}}$ describes a collection of harmonic
oscillators with mass $C$ and frequencies $\Omega_{\vec{k}}$ at temperature
$1/\beta$. Thus, the trial free energy per site is

$\frac{F_{\text{tr}}}{V}=\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left\\{\frac{1}{2}\Omega_{\vec{k}}+\frac{1}{\beta}\log\\!\left(1-e^{-\beta\Omega_{\vec{k}}}\right)\right\\}.$
(S21)

It follows that the variational inequality is $F\leq\tilde{F}_{\text{tr}}$,
where

$\frac{\tilde{F}_{\text{tr}}}{V}=\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left\\{\frac{1}{2}\Omega_{\vec{k}}+\frac{1}{\beta}\log\\!\left(1-e^{-\beta\Omega_{\vec{k}}}\right)\right\\}-\frac{1}{2}\sum_{j=\text{n.n.}\,i}\left\\{Je^{-\frac{1}{2}\langle[\theta_{i}-\theta_{j}]^{2}\rangle_{\text{tr}}}+\frac{\mu}{2}\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}\right\\}.$
(S22)

Inserting the self-consistency condition for $\mu$, this reduces to

$\frac{\tilde{F}_{\text{tr}}}{V}=\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left\\{\frac{1}{2}\Omega_{\vec{k}}+\frac{1}{\beta}\log\\!\left(1-e^{-\beta\Omega_{\vec{k}}}\right)\right\\}-\frac{\mu}{2}\sum_{j=\text{n.n.}\,i}\left[1+\frac{1}{2}\langle(\theta_{i}-\theta_{j})^{2}\rangle_{\text{tr}}\right].$
(S23)

The right hand side should be compared with the same quantity computed in the
disordered state. For the disordered state we take the free action, i.e.
$S_{\text{free}}$ is the same as $S_{\text{tr}}$ but with $\mu=0$. The
corresponding free energy per site is

$\frac{F_{\text{free}}}{V}=-\frac{1}{\beta}\log[\Theta(e^{-\beta/2C})],$ (S24)

where

$\Theta(x)\equiv\sum_{n=-\infty}^{\infty}x^{n^{2}}.$ (S25)

In addition, since the sites are decoupled in the free limit, $\langle
S-S_{\text{free}}\rangle_{\text{free}}=0$.

###  V.3 Summary

The self-consistency condition for $\mu$ can be written as

$\mu=J\exp\left[-\ \frac{\beta}{C}\ W\big{(}\beta\sqrt{\mu/C}\big{)}\right],$
(S26)

where

$W(z)=\frac{1}{4d}\int\\!\frac{d^{d}k}{(2\pi)^{d}}\
\frac{1}{z}\,[g(\vec{k})]^{1/2}\coth\\!\left(\frac{z}{2}\,[g(\vec{k})]^{1/2}\right)$
(S27)

and $g(\vec{k})$ is the structure factor

$g(\vec{k})=\sum_{a=1}^{d}\sin^{2}(k_{a}/2).$ (S28)

Eq. (S26) must be solved numerically for $\mu$.

The variational free energy per site can be written as

$\frac{\tilde{F}_{\text{tr}}}{V}=\frac{1}{\beta}\
Y\big{(}\beta\sqrt{\mu/C}\big{)}-\mu d,$ (S29)

where

$Y(z)=\int\\!\frac{d^{d}k}{(2\pi)^{d}}\left\\{\frac{z}{2}\,[g(\vec{k})]^{1/2}\left[1-\frac{1}{2}\coth\\!\left(\frac{z}{2}\,[g(\vec{k})]^{1/2}\right)\right]+\log\\!\left(1-e^{-z[g(\vec{k})]^{1/2}}\right)\right\\}.$
(S30)

This should be compared to the free energy per site of the disordered state:

$\frac{F_{\text{free}}}{V}=-\frac{1}{\beta}\log[\Theta(e^{-\beta/2C})],\qquad\Theta(x)\equiv\sum_{n=-\infty}^{\infty}x^{n^{2}}.$
(S31)

The calculation results are plotted as a function of $C/\beta$ with coupling
constants $CJ$ in Fig S6, where the transition temperature $C\beta_{c}\propto
T_{c}$ was determined by the self-consistent solution of Eq. (S26) with weak
coupling constants $CJ<0.55$ or suppressed due to
$F_{\text{free}}<\tilde{F}_{\text{tr}}$ with strong coupling constants
$CJ>0.55$. From these results, the normalized temperature
($T/T_{c}=\beta_{c}/\beta$) dependence of the normalized superfluid density
($n_{s}=\mu(CJ,T/T_{c})/\mu(CJ,0)$) are obtained as shown in Fig. S7.

Figure S6:  Temperature dependence of the calculated normalized superfluid
density $\mu(CJ,C/\beta)/J$ in the toy model of Eq. (S6). The grey data are
the calculated results where $F_{\text{free}}<\tilde{F}_{\text{tr}}$. Figure
S7:  Normalized temperature dependence of the superfluid density
$\mu(CJ,T/T_{c})/\mu(CJ,0)$ in the toy model of Eq. (S6). (a) With weak
coupling constants $CJ<0.55$, the ordered state is stable
($F_{free}>\tilde{F}_{tr}$). (b) With strong coupling constants $CJ>0.55$, the
superfluid density has a linear temperature dependence near $T_{c}$.

## References

  * [1] J. R. Kirtley, L. Paulius, A. J. Rosenberg, J. C. Palmstrom, D. Schiessl, C. L. Jermain, J. Gibbons, C. M. Holland, Y.-K.-K. Fung, M. E. Huber, M. B. Ketchen, D. C. Ralph, G. W. Gibson Jr. and K. A Moler, Supercond. Sci. Technol. 29, 124001(2016).
  * [2] J. R. Kirtley, B. Kalisky, J. A. Bert, C. Bell, M. Kim, Y. Hikita, H. Y. Hwang, J. H. Ngai, Y. Segal, F. J. Walker, C. H. Ahn, and K. A. Moler, Phys. Rev. B 85, 224518 (2012).
  * [3] Kogan V. G., Phys. Rev. B 68, 104511 (2003).
  * [4] R. Prozorov and R. W. Giannetta, Supercond. Sci. Technol. 19, R41 (2006).
  * [5] F. Gross, B. S. Chandrasekhar, D. Einzel, K. Andres, P. J. Hirschfeld, H. R. Ott, J. Beuers, Z. Fisk, and J. L. Smith, Z. Physik B - Condensed Matter 64, 175–188 (1986).
  * [6] S. Chakravarty, G.-L. Ingold, S. Kivelson, and A. Luther, Phys. Rev. Lett. 56, 2303 (1986).