Optical signatures of the superconducting Goldstone mode in granular aluminum: experiments and theory

TL;DR

This study combines THz experiments and theoretical modeling to reveal optically-active Goldstone modes in granular aluminum, showing how spatial inhomogeneities enable their optical detection and providing insights into collective superconducting excitations.

Contribution

The paper demonstrates the optical visibility of Goldstone modes in granular aluminum and models the inhomogeneities as a disordered Josephson array, linking experimental data with theory.

Findings

Observation of temperature-dependent sub-gap optical absorption.

Identification of Goldstone modes as the source of absorption.

Quantitative agreement between theory and experiment.

Abstract

Recent advances in the experimental growth and control of disordered thin films, heterostructures, and interfaces provide a fertile ground for the observation and characterisation of the collective superconducting excitations emerging below $T_c$ after breaking the $U(1)$ gauge symmetry. Here we combine THz experiments in a nano-structured granular Al thin film and theoretical calculations to demonstrate the existence of optically-active phase modes, which represent the Goldstone excitations of the broken gauge symmetry. By measuring the complex transmission trough the sample we identify a sizeable and temperature-dependent optical sub-gap absorption, which cannot be ascribed to quasiparticle excitations. A quantitative modelling of this material as a disordered Josephson array of nano-grains allows us to determine, with no free parameters, the structure of the spatial inhomogeneities induced by shell effects. Besides being responsible for the enhancement of the critical temperature with respect to bulk Al, already observed in the past, this spatial inhomogeneity provides a mechanism for the optical visibility of the Goldstone mode. By computing explicitly the optical spectrum of the superconducting phase fluctuations we obtain a good quantitative description of the experimental data. Our results demonstrate that nanograins arrays are a promising setting to study and control the collective superconducting excitations via optical means.

Figures

Equations

$$\frac{{\sigma }^{exc}}{{\sigma }_{dc}}\equiv \frac{{\sigma }_1^{exp}-{\sigma }_1^{BCS}}{{\sigma }_{dc}}$$

$${\mathcal{H}}_{XY}=-2\sum _{⟨i,j⟩}{J}_{ij}\left(S_i^{+}{S}_j^{-}+h.c.\right),$$

$$H_{HP}\simeq \frac{1}{2}\left[\sum _{i,\mu =x,y}{J}_{i,i+\mu }({\Delta }_{\mu }{\theta }_i{)}^2+\sum _{i,\mu =\pm x,\pm y}4J_{i,i+\mu }{L}_i^2\right]$$

$$S_{HP}$$

$$=$$

$$\omega =(4J{\xi }_0/\hslash )q=v_{s}q.$$

$$\hslash \bar{\omega }\sim 4J({\xi }_0\bar{q})\sim 4J$$

$$\sigma (\omega )$$

$$D$$

$$D_s$$

$${\sigma }_{reg}(\omega )$$

$$Z_{\alpha }=\frac{4}{E_{\alpha }}{\left(\sum _i{J}_{i,i+\hat{x}}{\Delta }_{\mu }{\varphi }_{\alpha ,i}\right)}^2,$$

$$J_{ij}=\frac{{\Delta }_i{\Delta }_j}{\beta }\frac{R_Q}{R_N}\sum _{l=-\infty }^{+\infty }\frac{1}{\sqrt{((\frac{(\pi (2l-1))}{\beta }{)}^2+{\Delta }_i^2)+((\frac{(\pi (2l-1))}{\beta }{)}^2+{\Delta }_j^2)}}$$

$$J_0=\frac{{\Delta }_0{R}_Q}{2R_N}\tanh \left(\frac{\beta {\Delta }_0}{2}\right),$$

$$R_Q=\frac{h}{4e^2}=6.45k\Omega$$

$$0<\omega /J<\gamma$$

$$J_s[K]=\frac{0.62\times \tilde{d}}{{\lambda }^2[\mu m^2]}$$

$$1=\frac{{\lambda }_{\mathrm{eff}}}{2}{\int }_{-{ϵ}_D}^{{ϵ}_D}\sum _n\frac{c_n{\nu }_{\gamma }({ϵ}^{\prime }-{ϵ}_n)F(T)}{\sqrt{{ϵ}^{\prime 2}+\Delta (R,T{)}^2}}d{ϵ}^{\prime }$$

$${\nu }_{\gamma }(ϵ)=\frac{\gamma }{\pi }\frac{1}{{ϵ}^2+{\gamma }^2}.$$

$$\gamma \sim \frac{4z{R}_Q}{R_N{\nu }_{TF}({ϵ}_{\mathrm{F}})},$$

$$⟨S_i^x⟩=\frac{1}{2}\tanh (\beta B_i),B_i=\sum _{j=i+z}{J}_{ij}\tanh (\beta B_j)$$

$$\begin{array}{c}{S}_i^x=1/2-a_i^{†}{a}_i\\ -S_i^z+i{S}_i^y={\left(1-a_i^{†}{a}_i\right)}^{1/2}{a}_i\simeq a_i\\ -S_i^z-i{S}_i^y=a_i^{†}{\left(1-a_i^{†}{a}_i\right)}^{1/2}\simeq a_i^{†}\end{array},$$

$${\mathcal{H}}_{XY}\simeq {\mathcal{H}}_{MF}+\frac{1}{2}\sum _{ij}\left(A_{ij}{a}_i^{†}{a}_j+B_{ij}{a}_i{a}_j\right)+h.c.,$$

$$A_{ij}$$

$$B_{ij}$$

$$S_i^x=\tanh (\beta B_i)/2-a_i^{†}{a}_i$$

$$-S_i^z+i{S}_i^y=\sqrt{\tanh (\beta B_i)}{\left(1-a_i^{†}{a}_i\right)}^{1/2}{a}_i\simeq \sqrt{\tanh (\beta B_i)}{a}_i$$

$$-S_i^z-i{S}_i^y=a_i^{†}\sqrt{\tanh (\beta B_i)}{\left(1-a_i^{†}{a}_i\right)}^{1/2}\simeq \sqrt{\tanh (\beta B_i)}{a}_i^{†}$$

$$A_{ij}$$

$$B_{ij}$$

$$a_i=\sum _{\alpha }\left(u_{\alpha i}{\gamma }_{\alpha }+v_{\alpha i}{\gamma }_{\alpha }^{†}\right)$$

$$a_i^{+}=\sum _{\alpha }\left(v_{\alpha i}{\gamma }_{\alpha }+u_{\alpha i}{\gamma }_{\alpha }^{†}\right)$$

$${\mathcal{H}}_{PS}=\sum _{\alpha }{E}_{\alpha }{\gamma }_{\alpha }^{†}{\gamma }_{\alpha }\text{ ,}$$

$${\theta }_i$$

$${\varphi }_{\alpha i}$$

Full text

I Introduction

Zero resistance at finite temperature in systems where momentum is not conserved, one of the defining features of a superconductor, is strictly related schrieffer ; nagaosa to the phase rigidity of the complex order parameter. Phase rigidity is a typical consequence of a spontaneous symmetry breaking: the Hamiltonian that describes the superconductor is invariant under a phase rotation of the electronic degrees of freedom but in the ground state the macroscopic order parameter chooses a particular value of the electronic phase, and therefore breaks spontaneously this rotational $U(1)$ symmetry. According to the Goldstone theorem, the collective excitation connecting the possible degenerate ground states must be massless at long wavelengthschrieffer ; nagaosa .

In principle, this collective mode, usually termed the Goldstone mode nambu ; goldstone , should manifest as a low-energy excitation. In the context of superconductivity it should appear as a sub-gap excitation. However, despite intensive research, it has not yet conclusively been observed experimentally. The standard explanation is based on the idea of AndersonAnderson1958 , proposed shortly after the BCS theory: Coulomb interactions, present in any material, boost the typical frequency of the Goldstone mode to the plasma energy scale, well above the energy gap, so that it cannot be thermally excited at temperatures below the superconducting (SC) critical temperature $T_{c}$. A second issue is that typical spectroscopic measurements probe the system in the long-wavelength regime where the phase mode is decoupled from the transverse electromagnetic field, so that one cannot observe it in the ac conductivity. However, these conclusions only hold for a spatially homogeneous SC state, since disorder and inhomogeneity can affect both the spectrum of the phase mode and its optical visibility.

In recent years, due to experimental advances in the growth and control of SC thin films, the electromagnetic response of conventional superconductors have been studied with unprecedented precision and in a broad range of disorder strengths. The improved experimental resolution revealed that in some conventional $s$-wave superconductors, as NbN, InOx and granular Al, a finite absorption can be found even below the threshold $2\Delta$ for the Cooper-pair breaking, with $\Delta$ being the SC gap armitage_prb07 ; driessen12 ; driessen13 ; bachar_jltp14 ; frydman_natphys15 ; samuely15 ; armitage_prb16 ; scheffler16 ; pracht_prb16 . One possible interpretationstroud_prb00 ; cea_prb14 ; trivedi_prx14 ; seibold_cm17 of these experiments points out to the relevance of the Goldstone mode, made optically active by the spontaneous inhomogeneity of the SC ground state that has been observed in strongly-disordered filmssacepe08 ; sacepe10 ; mondal10 ; sacepe11 ; chand12 ; kam13 ; roditchev13 ; roditchev_natphys14 ; samuely16 ; leridon_prb16 ; brun_review17 . This interpretation is still under debate for two main reasons. From one side, the explicit calculation of the phase-mode absorption has been performed so farstroud_prb00 ; cea_prb14 ; trivedi_prx14 ; seibold_cm17 under the assumption that long-range Coulomb forces should not be included. This relies on the expectation that the optical conductivity is the response to the local field, so it is irreducible with respect to the Coulomb interactionbelitz_prb89 ; cea_prb14 , i.e. it only depends on the undressed sound-like phase spectrum. Nonetheless, disorder could still mix the reducible and irreducible response, and an explicit estimate of this effect is still lacking. From the other side, other mechanisms could potentially induce sub-gap excitations in the proximity of the superconductor-to-insulator transition (SIT). One interesting proposalfrydman_natphys15 is that large enough disorder can push the SC system towards an effective strong-coupling regime, where the energy scale of the Higgs mode would lie below the one for Cooper-pair breaking, making it visible as a sub-gap excitation. However, explicit calculationscea_prl15 within disordered fermionic models did not find yet evidence of a sharp, subgap Higgs mode in the strong-coupling regime emerging near the SIT.

It is then clear that an experimental confirmation of the role of Goldstone modes would require a weakly-coupled and inhomogeneous superconductor, where the inhomogeneity is not driven by too large disorder, so that the sub-gap features cannot be attributed to other unrelated effects. It would be also desirable to establish a quantitative relation between the experimental source of inhomogeneity and the parameters of the theoretical model, that presently is not known for inhomogeneities induced by impurities. A promising settings for this analysis is provided by films made of well-coupled nanograins. Indeed, recent research Bose2010 ; Garcia-Garcia2008 ; Garcia-Garcia2011 ; Brihuega2011 ; Shanenko2006 ; Shanenko2007 in single isolated superconducting nano-grains of Sn has found that shell effects, induced by fluctuations in the density of states, lead to large changes in the SC gap by very small changes in the grain size, provided that the grain shape is sufficiently symmetric. As a result, a Josephson array of SC nano-grains can be highly inhomogeneous even when the SC transition occurs very far from the SIT. Unlike inhomogeneities induced by impurities it is possible to know experimentally Lerer2014 ; Deutscher1973 ; Deutscher1973a ; pracht_prb16 , with relative accuracy, the grain size, its distribution, and the density of grains, that control the resistivity of the material and the inhomogeneity of its SC properties. An additional advantage of the nanograins, as compared to strongly-disordered superconductors and arrays made of artificial Josephson junctions, is that in this case the capacitive charging energy of each grain dominates over the junction capacitance. In this situation the long-range part of the Coulomb interaction is automatically screened, and the phase mode remains an acoustical mode at long wavelengthfazio_review01 ; halperin_prb89 ; ribeiro2014 ; garcia_prb14 .

A prototype system belonging to this category is granular Al. The typical grain size, a few nanometers, and its shape, spherical, is optimalgarcia_prb14 to observe strong shell effects and consequently strong inhomogeneities in the array even in the low- resistivity region, where the average coupling among grains is strong enough to give a dc conductivity consistent with metallic behavior. An indirect proof of the emerging inhomogeneity is provided by the experimentally observed enhancement of the critical temperature Abeles1966 ; Deutscher1973 ; Lerer2014 ; pracht_prb16 with respect to bulk Al. Indeed, it has been arguedgarcia_prb14 that the local enhancement of the superconducting gap due to shell effects can increase, within a percolative scheme, the critical temperature of the array.

A natural question to ask is whether inhomogeneities induced by shell effects can also lead to the observation of collective Goldstone excitations. Here we answer this question positively. More specifically, we first provide direct experimental evidence of a broad sub-gap resonance in the ac conductivity of granular Al film deep in the metallic region. In this sample corrections to the bulk mean-field due to shell effects are still substantial to induce inhomogeneity in the system, but the film is metallic enough to suppress the corrections to the critical temperature due to quantum phase fluctuationsgarcia_prb14 . We model granular Al as an array of nano-grains where the inhomogeneity of the Josephson coupling between grains is computed microscopically starting from the grain-size distribution. The obtained distribution of the local coupling is then used as an input to compute both the $T_{c}$ of the array and the optical response of the disordered phase modes. The same model can then explain two striking experimental observations in the system: (i) the enhancement of $T_{c}$ with respect to bulk Al and (ii) the emergence of an optically-active Goldstone mode, visible as an extra sub-gap absorption. The good quantitative agreement between the theoretical calculations and the experiments provides a strong support to the relevance of Goldstone modes in inhomogeneous superconductors, opening interesting perspectives for applications in artificially-designed inhomogeneous systems.

The paper is organized as follows: In Sec. II we describe the experimental setup and the measurements in a granular Al sample in the metallic regime. Sec. III describes the theoretical model used to compute the SC properties of the arrays. In Sec. IV we compare explicitly the theoretical results with the experimental data for the THz conductivity. The concluding remarks are presented in Sec. V. Additional technical details on the theoretical calculations are provided in Appendix A and B.

II Experiments

To access the frequency-dependent ($\nu$) complex conductivity $\sigma(\nu)=\sigma_{1}(\nu)+i\sigma_{2}(\nu)$, we performed phase-sensitive measurements Pracht2013 of the complex transmission coefficient of tunable and coherent THz radiation passing through a granular Al film of $d=40$ nm thickness. The film was grown on a MgO2 dielectric $10\times 10$ mm2 substrate held at 77 K via thermal evaporation in partial pressure of O2. The O2 concentration while deposition was adjusted such to produce a granular morphology characterized by a dc-resistivity of $\rho_{dc}=0.26\times 10^{-3}\Omega$cm at 5 K and a strongly enhanced $T_{c}$ of 2.74 K. We simultaneously fit the amplitude and phase of the transmitted THz radiation with Fresnel equations via $\sigma_{1}$ and $\sigma_{2}$ (detailed information on the analysis procedure is found in Pracht2013 ). Note that this approach is not based on any microscopic model for the charge carrier dynamics. The low-temperature data for several samples at different disorder levels have been shown in a previous publication pracht_prb16 .

As $\sigma_{1}$ equally probes the coupling of light to single-particle excitations as well as to optically-active collective modes of the order parameter, it can be difficult to clearly disentangle different absorption channels. To outline the effect of the Goldstone mode, we focus on the real (dissipative) part $\sigma_{1}$ and define the excess conductivity with respect to the Mattis-Bardeen BCS prediction

[TABLE]

where $\sigma_{dc}$ is the dc conductivity measured right above $T_{c}$. We infer $\sigma_{1}^{BCS}$ by fitting $\sigma_{1}^{exp}$ at frequencies above the energy gap, i.e. the excitations into the quasiparticle continuum, to the standard Mattis-Bardeen functional mb ; pracht_prb16 , see Fig.  1. Interestingly, this yields a surprisingly good description of $\sigma_{1}^{exp}(\nu)$ for $\nu\geq 2\Delta/(hc)$ and, at the same time, reveals (and quantifies) the notable excessive absorption $\sigma_{1}^{exc}(\nu)$ at frequencies below twice the SC gap. As one can see in Fig. 1, as the temperature rises from $T=1.6$ K towards $T_{c}=2.74$ K, the extra absorption gets continuously suppressed, but nonetheless it still represents a significant contribution to $\sigma_{1}^{exp}$ aside the thermally-excited quasiparticles. On the other hand, the remaining features of the spectra are rather conventional. For example, the inset of Fig. 1 shows the temperature dependence of the gap, as extracted from the Mattis-Bardeen fit. It can be very well fitted to a BCS-ilke behavior, giving an extrapolated value at $T=0$ equal to $\Delta=0.4$ meV, so that $\Delta/T_{c}=1.78$ has the conventional weak-coupling value, despite the enhancement of the $T_{c}$ as compared to bulk Al.

As observed beforepracht_prb16 the presence of extra sub-gap absorption is a general feature of granular Al samples regardless the level of inter-grains coupling, measured by the value of the normal-state resistivity. This has to be contrasted with homogeneously disordered films of conventional superconductorfrydman_natphys15 ; driessen12 ; armitage_prb16 ; practh2012 , where significant deviations from the Mattis-Bardeen prediction only occurs at relatively large disorder levels, where the dc conductivity displays an insulating behavior. To make a quantitative comparison let us consider for example the NbN films of Ref. armitage_prb16 . Here a consistent sub-gap absorption, comparable to the one reported in Fig. 1, appears only for film resistance $\rho_{dc}=1.2\times 10^{-3}\Omega$cm, with film thickness $60-120$ nm. For our granular Al sample the sub-gap absorption can be clearly distinguished already for a resistivity value about one order of magnitude smaller, $\rho_{dc}=0.26\times 10^{-3}\Omega$cm (and $d=40$ nm).

III Theoretical model

III.1 Optical response of disordered phase modes

As a starting point for the theoretical analysis of the optical conductivity we need a proper model for the phase degrees of freedom in a inhomogeneous superconductor. Our main aim is to study low-energy properties for which it is not necessary to consider amplitude fluctuations. We will also consider the quasi-two dimensional limit for computational simplicity. A reasonable description for the phase fluctuations in a granular system is then provided by the two-dimensional (2D) quantum XY model matsubara ; anderson59 for $S=1/2$ pseudospins with bonds on a square lattice:

[TABLE]

where the sum is extended to all the pairs of nearest neighbour sites and $J_{ij}$ is the hopping that depends on the resistance of the model and the value of the SC gap, i.e. the amplitude of order parameter, in neighboring grains. The local gap can be computed in the mean-field limit as a function of the grain sizegarcia_prb14 . We first discuss the outcomes of the model (2) for a generic $J_{ij}$. In the second part of the section we describe the details of the calculation of $J_{ij}$.

The idea underlying the present approach is the well-known Anderson mappinganderson59 between a $s$-wave superconductor and the pseudospin Hamiltonian (2). Here the $S_{i}$ are spin $1/2$ operators, such that the in-plane spin component describes the pairing operator, $S^{-}_{i}=c_{i\downarrow}c_{i\uparrow}$, while the $z$ component describes the local density, $S^{z}_{i}=\frac{1}{2}\left(\sum_{\sigma}c^{\dagger}_{i\sigma}c_{i\sigma}-1\right)$, with $c_{i\sigma}$ fermionic annihilation operators. Thus within such a spin-like picture of superconductivity the SC order appears as the spontaneous magnetization of the $XY$ component, and $J_{ij}>0$ is the hopping amplitude, which represents the energetic gain to move a Cooper pair from a given site $j$ to a nearest neighbor site $i$.

An exact treatment of Eq. (2) is beyond current analytical techniques. In order to proceed we compute the mean-field ground state, and fluctuations with respect to the mean-field ground state, by mapping the spins into bosonic operators by means of the usual Holstein-Primakov approximation (see Appendix B). After the mapping, the evaluation of the superconducting phase-mode spectrum is equivalent to the evaluation of the spin-wave spectrum in an ordinary ferromagnetic spin model. In practice, this is equivalent to study the following quantum phase-only model:

[TABLE]

where $L_{i}$ are the quantum operators canonically conjugated to the phase operators $\theta_{i}$ and $\Delta_{\mu=x,y}$ is the discrete phase gradient in the $\mu$ direction. In the homogeneous case $J_{i,i+\mu}=J$ one can easily see that Eq.  (3) describes a sound-like phase mode. Indeed, deriving the action corresponding to the Hamiltonian via the usual identification $16JL_{i}^{2}\rightarrow(\hbar\partial_{t}\theta)^{2}/16J$ one immediately obtains

[TABLE]

so that the phase-modes spectrum is sound-like with typical velocity $v_{s}={4J\xi_{0}}/\hbar$:

[TABLE]

Here the scale $\xi_{0}$ is the coherence length, which represents the typical length scale over which the coarse-grained model (2) is valid. For Al nanograins the effective quantum phase model is usually takenribeiro2014 ; garcia_prb14 as a Josephson-junction array in the limit where the self-capacitance dominates over the capacitance of the junctions, so that non-local $L_{i}L_{j}$ quantum terms in Eq. (3) are absent. In this regime the long-range part of the Coulomb force is screened, and the phase mode is not converted into a plasmonfazio_review01 ; halperin_prb89 . For samples in the metallic limit, as the one we are considering, one should also account for the additional dissipation mechanism provided by tunnelling of electrons between neighbouring grains.halperin_prb89 ; ribeiro2014 This effect renormalizes the local capacitance to the effective scale $1/J$,ribeiro2014 ; garcia_prb14 exactly as found in Eq. (3) starting directly from the quantum pseudospin model (2).

In weakly-disordered and homogeneous BCS superconductors the phase modes are optically inert at Gaussian level. Indeed, despite the fact that they are crucial to restore the gauge invariance of the longitudinal responseschrieffer ; cea_prb14 ; seibold_cm17 , they do not contribute to the transverse physical response, which reduces to the superfluid contribution at $\omega=0$ . As a consequence the real part of the optical conductivity is well described by the usual Mattis-Bardeen expressionmb , with a superfluid delta peak at $\omega=0$, followed by a finite absorption above the threshold $2\Delta$ where Cooper-pair breaking by the electro-magnetic field begins. However, this picture is modified at strong disorder, and in general when the system displays a spatial inhomogeneity, since the disordered phase modes give rise to finite-frequency absorptionstroud_prb00 ; cea_prb14 ; trivedi_prx14 already in the Gaussian approximation. A simple argument to understand this effect is sketched in Fig. 2. In granular Al, shell effects induce a spatial inhomogeneity in the local stiffness of the model (2). In this situation the phase fluctuations of the local SC order parameter induced by the incoming radiation are larger when the phase rigidity (i.e. $J_{ij}$) is smaller. This leads to an inelastic response which results in a partial absorption of the e.m. field, i.e. in a finite value of the real part of the ac optical conductivity. In other words, the phase modes at the characteristic finite momentum $\xi_{0}\bar{q}\sim 1$ set by the inhomogeneity become mixed to the transverse physical response at ${\bf q}=0$. The typical energy scale of the absorption is then

[TABLE]

so that it is given by the overall scale $J$ of the local Josephson couplings in Eq. (2), which will be in the following our main fitting parameter.

The above qualitative arguments can be made quantitative by computing explicitly the optical conductivity of the model (2) within the Holstein-Primakoff approximation (3). Following the general procedure outlined in Ref. cea_prb14 the optical conductivity e.g. in the $x$ direction is given by (details are given in Appendix B)

[TABLE]

with $N$ number of lattice sites and $D$, $D_{s}$ denote the diamagnetic and superfluid weight, respectively. Here $E_{\alpha}$ are the eigenvalues of the disordered phase spectrum of the model (3), and $Z_{\alpha}$ represents their electrical dipoles. These are given explicitly by

[TABLE]

where the $\phi_{\alpha,i}$ are connected to the eigenvectors of each $E_{\alpha}$ mode (see Eq. (33) below). When the system is homogeneous, i.e. $J_{i,i+\mu}=J$, each $Z_{\alpha}$ is proportional to the total gradient of the phase over the system, which vanishes for periodic boundary conditions. In this case, the optical conductivity Eq. (8) consists only of the superfluid peak at $\omega=0$, with $D_{s}=D=4J$. However disorder in the local Josephson couplings between grains induces a finite electrical dipole for the phase modes, leading to a regular part $\sigma_{reg}$ of the optical conductivity (8), responsible for the finite-frequency absorption. The exact form of this absorption depends in general on the disorder of the local stiffnesses $J_{i,j}$ used in the model (2). To make contact with the structure of granular Al, we will compute in the next section the distribution of the $J_{ij}$ due to shell effects in an array of nanograins.

III.2 Calculation of $J_{ij}$ and P($J_{ij}$)

As was mentioned previously, we model granular aluminum as an array of nano-grains with different sizedeutscher_prb80 ; garcia_prb14 . This implies that, because of shell-effects, each grain has a different superconducting gap and critical temperature, leading in turn to an inhomogeneous distribution of local stiffnesses. The strength of the inhomogeneities is controlled by the effective electron-phonon coupling, the grain size distribution and the strength of the coupling among grains. Typically the smaller and more isolated the grain is, the stronger are the inhomogeneities of the sample. In situations in which BCS theory applies, the distribution of the local stiffness $J_{ij}$ of the model (2) is given by the Ambegaokar-Baratoff expression Ambegaokar1963 :

[TABLE]

where $\Delta_{i}(L)$ is the value of the superconducting order parameter of a grain of size $L$ and $\beta=(k_{B}T)^{-1}$ with $k_{B}$ the Boltzmann’s constant. We compute $\Delta_{i}(L)$ by solving numerically the BCS gap equation for a spherical grain of size $L$ where we also took into account that the grain is open, namely, electrons can hop from one grain to another. The hopping is more likely as the normal state resistivity becomes smaller. Effectively this hopping smears out the spectral density that weakens shell effects, especially in the deep metallic regime where the electron hops to another grain very quickly. Details of the calculation are given in the Appendix A.

The distribution of $J_{ij}$ is ultimately determined by the experimental distribution of grain sizes. From experiments Deutscher1973 ; Lerer2014 we know that for the range of resistances of interest the distribution of sizes is close to a log normal distribution with an average diameter $2$ nm and variance of $0.5$ nm. This is the distribution that we will use in the rest of the paper. Despite the broad distribution of sizes, the fact that the typical size is only $2$nm cast doubts on the applicability of BCS in most grains as for $L\sim 2$ nm the mean level spacing is much larger than the superconducting gap, violating thus Anderson’s criterion. However we note that this criterion is only applicable for isolated grains where the spectrum is discrete with no imaginary part. The smearing of the spectral density mentioned above has also the effect of bringing back enough spectral weight close to the Fermi energy so that a mean field approach is applicable. This is especially true in the good metallic limit we are interested in.

We are now ready for the computation of the distribution $P(J_{ij})$. We will normalise it with respect to the hopping $J_{0}$ in the limit of a homogeneous array:

[TABLE]

where $\Delta_{0}$ is the bulk value of the gap in the homogeneous case and we defined the quantum of resistance as

[TABLE]

while $R_{N}$ is the normal-state resistance per square of the array. For the sample shown in Fig. 1 $\rho_{dc}=263\mu\Omega$cm, that for $d=40$ nm gives $R_{N}=65.75\Omega\simeq 0.01R_{Q}$, which is deep in the metallic region. In Fig. 3 we depict the corresponding distribution $P(J_{ij})$ for different temperatures. It is strongly bimodal with a peak at zero hopping corresponding to bonds where the gap vanishes. The other peak has relatively fat tails centred at a $J$ larger than the bulk one at that temperature. Once the distribution of local stiffnesses is known, we can also estimate the global critical temperature as induced by percolation of spheres, which is known to be a good approximation in the limit of small resistance, i.e $R_{N}/R_{Q}\ll 1$ garcia_prb14 . The existence of large tails in the distribution of the stiffness implies that the global $T_{c}$ given by percolation is substantially enhanced. Indeed, by mapping the local $J_{ij}$ to local $T_{c}$ values one can estimate the $T_{c}$ of the array as the temperature where the percolation threshold is reached, as discussed in Ref. garcia_prb14 . By using this approach we get here $T_{c}\sim 1.75T_{c}^{\rm bulk}$. Even if this is still smaller than the experimental value $T_{c}/T_{c}^{0}=2.3$ ($T_{c}=2.74$ K) it goes in the right direction.

IV Comparison between theory and experiments

After computing the distribution of the local stiffness, we use it to obtain the optical conductivity according to Eq. (8) above. As we mentioned previously, we need to fix the overall scale $J$ of the stiffness, since in general $\sigma(\omega)$ will be different from zero in a range of values

[TABLE]

where $\gamma$ depends on the model for the disorder, i.e. on the form of $P(J_{ij})$. The experiments observe a substantial sub-gap absorption in a range of energies such that $\gamma J\lesssim 2\Delta$, i.e. below the threshold for quasiparticle creation. As a consequence we will fix $J$ in order to have our simulation to overlap with the experimental data. We note that our model only describes the absorption by Goldstone modes. Thus, to make contact with the experimentally-determined excess conductivity (1) we must use a proper normalization. Since our model is taken for computational simplification as purely two-dimensional, the conductivity is a multiple of the 2D quantum of absorption $\sigma_{0}=e^{2}/\hbar=0.25\times 10^{-3}\Omega$. To translates it in a 3D conductivity we should then divide our 2D conductivity by a transverse length scale $\tilde{d}$ of the order of the distance between layers. Notice that, for the same reason, our overall scale $J$ is expected to be quantitatively lower than the stiffness $J_{0}$ of the whole array, as given by Eq. (14) above.

To test the effects of the SC inhomogeneity induced by the granular nature of the film we compare in Fig. 4 the optical conductivity obtained from the $P(J_{ij}/J)$ distribution of Fig. 3 with the experimental data, by using the overall strength of the equivalent 2D homogeneous coupling $J$ as the only free parameter. Even though the distribution of the $P(J_{ij})$ obtained by considering shell effects is itself temperature dependent due to the temperature variation of the SC order parameter $\Delta_{i}$, we fixed here its shape at $T=0$ and let the system evolve thermally according to the mean-field solution of the model (2). This approach is equivalent to account for the effects due to phase fluctuations alone, without including also the ones due to thermal suppression of the pairing, which is not present in the phase-only model (2). Details on the $T$ dependence are given in Appendix B. As one can see in Fig. 4 the calculations account reasonably well for the experimental data, considering that the distribution of the $P(J_{ij}/J)$ has been determined microscopically without any free parameter, except for the overall scale $J$ setting the local 2D stiffness, that we determine as $J=0.4$cm*-1*. The transverse length scale is relatively small, $\tilde{d}=0.6$ Å, but still consistent with our 2D limit.

By closer inspection of Fig. 3 one sees that roughly speaking the $P(J_{ij})$ resembles a bimodal distribution, with about $p=0.2$ spectral weight in the bonds with $J_{ij}=0$. We then tested the possibility to reproduce the experiments with a similar but simpler diluted distribution, where $J_{ij}$ has a probability $p=0.2$ of having value $J_{ij}=0$ and probability $1-p$ of having a value $J_{ij}=J$. The result for the bimodal distribution is shown in Fig. 5. For the diluted model the support of the $P(J_{ij}/J_{0})$ is more compact, so the absorption occurs for a smaller value of $\gamma$ in Eq. (16). This implies that both $J=0.9$ cm*-1* and $\tilde{d}=1.24$ Å  turn out to be about a factor of two larger. At present, the available experimental data do not allow to seriously discriminate between the two distributions, even though the diluted model displays a slightly better agreement for the temperature evolution.

Finally, we can use the extracted information on the values of $J$ and $\tilde{d}$ to make contact with the measured value of the magnetic penetration depth $\lambda$, that can be obtained by extrapolating the measured $\sigma_{2}(\omega)$ at zero frequency. More specifically, the inverse penetration depth is connected to the superfluid stiffness $J_{s}=\hbar n_{s}^{2d}/4m$ of the sample bypracht_prb16

[TABLE]

with $\tilde{d}$ given in Å. According to Eq. (8) above, the superfluid stiffness $J_{s}$ is given by the values of $D_{s}$ obtained numerically, as $J_{s}=D_{s}/4$. For the case of the microscopic $P(J_{ij})$ shown in Fig. 4 we obtain that $J_{s}=0.625J=0.35$ K, so that $\lambda^{-2}=0.93$ ($\mu$m)-2, while for the diluted model we have that $J_{s}=0.57J=0.72$ K but since $\tilde{d}$ is larger we get analogously $\lambda^{-2}=0.9$ ($\mu$m)-2. In both cases we then obtain a very good agreement with the measured value $\lambda^{-2}_{exp}=1.06$ ($\mu$m)-2, reinforcing the consistency of the overall theoretical analysis.

V Conclusions

Our work presents a comprehensive experimental and theoretical analysis of the unconventional superconducting THz response of granular Al. By measuring the complex conductivity of the system across the critical temperature we have clearly established the existence of a sizeable optical absorption that cannot be ascribed neither to thermally excited quasiparticles nor to the Cooper pairs broken by the electromagnetic field. Such extra absorption represents a striking violation of the usual Mattis-Bardeen paradigm for disordered BCS superconductors. At the same time, the critical temperature of the film greatly exceeds the one of bulk Al. We have shown that both features can be quantitatively described by modelling granular aluminum as an array of Al nano-grains, where the local Josephson couplings $J_{ij}$ between neighboring grains present a wide distribution with fat tails. Such inhomogeneity, induced by shell effects in each nano-grain and by the distribution of grain sizes, can be computed microscopically by using as input parameters the well-known BCS values of bulk Al. The inhomogeneity of the array has two crucial consequences. On one side, the non-vanishing probability of large local values of the SC gap and then of the local Josephson coupling $J_{ij}$ explains, within a percolative scheme, the enhancement of $T_{c}$ with respect to the homogeneous case. On the other side the inhomogeneity makes the SC phase mode optically active, explaining the anomalous sub-gap absorption. While the SC inhomogeneity has been also proven to emerge in homogeneously disordered films in proximity of the insulating state, in granular Al film it is a natural consequence of the confinement of superconductivity at the nanoscale. As such, it also manifests in arrays of well-coupled grains, where the overall resistivity still preserves a metallic behavior. The nanostructure has also another advantage, that has been emphasized in the above discussion: it makes the local charging effects stronger, screening the long-range Coulomb forces and leaving intact the sound-like dispersion of the Goldstone mode. The inherent inhomogeneity can then make the phase mode optically active, explaining the anomalous sub-gap absorption observed experimentally.

Our results provide thus strong evidence that the Goldstone mode can be observed in metallic superconducting nano-grains, provided that the superconducting state is sufficiently inhomogeneous. The optical signatures of the phase modes are not universal, but depend on the probability distribution of the local, inhomogeneous Josephson couplings. On a wider perspective, one can imagine to design an artificial array of Josephson junctions to explore the evolution of the phase response as a function of the Coulomb screening, controlled by the relative strength of the local charging energy with respect to the junction capacitance. This approach would ultimately address the long-standing issue of the interplay between inhomogeneity and Coulomb interactions in homogeneously disordered systems, and their effect on the optical visibility of the Goldstone mode.

Acknowledgements.

This project has been supported by the Deutsche Forschungsgemeinschaft (DFG), the German-Israeli-Foundation (GIF), the Italian MIUR (PRINRIDEIRON-2012X3YFZ2), the Italian MAECI under the Italian-India collaborative project SUPERTOP-PGR04879, and the Graphene Flagship. U.S.P. acknowledges a fellowship of the Studienstiftung des deutschen Volkes. AMG acknowledges partial financial support from a QuantEmX grant from ICAM and the Gordon and Betty Moore Foundation through Grant GBMF5305

Appendix A Numerical calculation of the superconducting gap in isolated and open spherical nano-grains

We have carried out a fully numerical calculation of the BCS energy gap with an effectively complex spectrum that accounts, in a spherical grain of radius $R$, for the possibility of tunnelling to other grains. The resulting gap equation is

[TABLE]

with $\lambda_{\rm eff}=\lambda(1+\frac{3\pi}{4k_{F}R})$. This additional factor accounts for the leading correction to the bulk coupling constant as a consequence of non-trivial matrix elements Garcia-Garcia2008 ; Garcia-Garcia2011 ,

[TABLE]

The spectrum $\epsilon_{n}$ is obtained from the zeros of the Bessel function of order $n+1/2$. In the Fermi golden rule approximation,

[TABLE]

and $F(T)=\tanh(\sqrt{\epsilon^{\prime 2}+\Delta(R,T)^{2}}/2T)$.

For the case of granular Al considered in the present work we use the well-known parameters for bulk Al, i.e. $\lambda\approx 0.18$, $\epsilon_{D}=36.4$meV, $k_{F}\approx 17.5{\rm nm}^{-1}$ and $\epsilon_{F}=11.63$ eV, $\xi\sim 1600$nm. We can then compute numerically the superconducting gap $\Delta(R,T)$ as a function of the grain size distribution and of the normal resistance $R_{N}$, that is controlled by the tunnelling rate. Once determined the local $\Delta_{ij}$ the Josephson coupling between neighboring grains are given by Eq. (13) in the main text.

Appendix B Calculation of the optical conductivity in the disordered $XY$ model

Let us first of all show the derivation of Eq. (3) from the quantum pseudospin model (2). As a starting point we compute the mean-field ground state, obtained by assuming that the spins align in the plane along a given, say $x$, direction $\langle S_{i}^{x}\rangle\neq 0$. By introducing the Weiss field $B_{i}=\frac{1}{2}\sum_{j}J_{ij}\langle S_{j}^{x}\rangle$ the Hamiltonian (2) can be approximated at mean-field level with $\mathcal{H}_{MF}=-\frac{1}{2}\sum_{i}B_{i}S_{i}^{x}$. One can then easily derive the self-consistent equations

[TABLE]

where the sum over $j$ in the above equation is extended to all the $z=4$ nearest-neighbours. Once the ground state has been determined, fluctuations above it can be described in the Holstein-Primakoff approximation. At $T=0$ this corresponds to introducing the bosonic annihilation (creation) operators $a_{i}$ ($a^{\dagger}_{i}$), related to the spins by:

[TABLE]

where we oriented explicitly the quantization axis along $x$. Substituting (22) into (2), at Gaussian level we get the following Hamiltonian:

[TABLE]

with

[TABLE]

At finite temperature the form of the Holstein-Primakov relations (22) has to be modified, to remove the contribution of the thermal excitations leading to the suppression of the order parameter (21) from the definition of the bosonic operators $a$. This implies that Eq.s (22) have to be replaced by

[TABLE]

At the same time we have to rescale the matrix elements (24)-(25) appearing in (23) as:

[TABLE]

The Hamiltonian (23) can be diagonalized via a standard Bogolubov transformation for bosons:

[TABLE]

so that

[TABLE]

with the energies $E_{\alpha}\geq 0$ and the coefficients $u$, $v$ determined by solving the secular equations: $\left[\gamma_{\alpha},\mathcal{H}_{PS}\right]=E_{\alpha}\gamma_{\alpha}$.

To describe the bosonic excitations in terms of collective modes and complete the mapping into the quantum Hamiltonian (3) we need to make one step further by defining the ”phase operators” as:

[TABLE]

The $\theta_{i}$’s are the quantum operators associated with the phase fluctuations of the SC order parameter: this identification, that will be formally justified below by the coupling of the Gauge field in the original Hamiltonian (2), can be understood from a semi-classical argument. Since, as we explained below Eq. (2), the operator associated with to the local order parameter is $S^{-}=S^{x}-iS^{y}$ if we put $S^{-}\simeq|\Delta_{i}|\left(1+i\theta_{i}\right)$, this obviously implies (32), since in the pseudospin mapping $|\Delta_{i}|=\langle S^{x}_{i}\rangle$. Once the phase operators $\theta_{i}$ have been identified, we can equally define their conjugate operators

[TABLE]

with $\ell_{\alpha}\equiv(u_{\alpha i}+v_{\alpha i})/\sqrt{2}$, which satisfy the usual commutation relations $[\theta_{i},L_{j}]=i\delta_{ij}$. By means of the definitions (32)-(34) one can easily derive Eq. (3), where the coefficient of the $(\nabla_{\mu}\theta_{i})^{2}$ term is replaced in general at finite temperature by $J_{i,i+\mu}\rightarrow J_{i,i+\mu}\tanh(\beta B_{i})/2\tanh(\beta B_{i+\mu})/2$.

To derive the current-current response function we need to couple the superconductor to the gauge field A. By exploiting once more the mapping between the fermions and the pseudospin operators, and by using the well-known Peierls substitution $c_{i}\rightarrow c_{i}e^{-ie\int^{{\bf r}_{i}}{\bf A}\cdot d{\bf l}}$ for the fermionic operators, we immediately get that the gauge field enters the pseudospin Hamiltonian (2) as

[TABLE]

with the factor of two taking into account the double charge of a Cooper-pair and $A_{i}^{\mu}\equiv{\bf A}({\bf r}_{i})\cdot\hat{\bf r}_{\mu}$. Here $e$ is the electron charge, and we set the sound velocity $c=1$. By means of the relation (35) one can derive the current operator in terms of the pseudospin operators, and then use again the Holstein-Primakoff transformations (26) to express it in terms of the spin-wave bosonic operators $a_{i},a_{i}^{\dagger}$, or equivalently the phase-momentum operators (32)-(34). In particular, it can be easily shown that (35) is equivalent to substitute $\Delta_{\mu}\theta_{i}\rightarrow\Delta_{\mu}\theta_{i}-2eA^{\mu}_{i}$ inside the Hamiltonian (3), thus justifying a posteriori the role of the $\theta_{i}$ in representing the phase of the SC order parameter.

The operator associated with the current flowing across the link $(i,i+\mu)$ then reads:

[TABLE]

where, as usual, the first term is the diamagnetic part, linearly proportional to the gauge field, while the second one defines the paramagnetic current operator. By using the Kubo formula, the current in linear-response theory can then be written as

[TABLE]

where the electromagnetic kernel $K^{\mu\nu}_{ij}$ is computed explicitly as:

[TABLE]

In the case of a uniform field along the $x$ axis $\mathbf{A}=A(t)\hat{x}$, we can define the optical conductivity by averaging over the space coordinates: $\sigma(\omega)=-\frac{i/N}{(\omega+i0^{+})}\sum_{i}\frac{\langle I^{x}_{i}(\omega)\rangle}{A(\omega)}$, with $N$ the number of lattice sites. The real part of the conductivity reduces then to Eq. (8) in the main text. To compute it we determine numerically at each temperature the eigenvalues $E_{\alpha}$ and the eigenvectors (29)-(30) of the Hamiltonian (31) for a given disorder configuration of the $J_{ij}$ couplings (using $N=30$), and we then average the result over 100 disorder configurations. The self-consistent equation (21) gives a vanishing order parameter for a given mean-field temperature $T_{MF}$. To make a direct comparison with the experimental data we then plot our numerical results as a function of $T/T_{MF}$ for the corresponding value of $T/T_{c}$ of the experimental measurements.