# Deconstructing Magnetization Noise: Degeneracies, Phases, and Mobile
Fractionalized Excitations in Tetris Artificial Spin Ice

M. Goryca1,2, X. Zhang3, J. Ramberger4, J. D. Watts4,5, C. Nisoli6, C.
Leighton4, P. Schiffer3,7, S. A. Crooker1 1National High Magnetic Field Lab,
Los Alamos National Laboratory, Los Alamos, NM 87545, USA  2Institute of
Experimental Physics, Faculty of Physics, University of Warsaw, Pasteura 5,
02-093 Warsaw, Poland  3Department of Applied Physics, Yale University, New
Haven, CT 06520, USA  4Department of Chemical Engineering and Materials
Science, University of Minnesota, Minneapolis, MN 55455, USA  5School of
Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA
6Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545,
USA  7Department of Physics, Yale University, New Haven, CT 06520, USA

###### Abstract

Direct detection of spontaneous spin fluctuations, or “magnetization noise”,
is emerging as a powerful means of revealing and studying magnetic excitations
in both natural and artificial frustrated magnets. Depending on the lattice
and nature of the frustration, these excitations can often be described as
fractionalized quasiparticles possessing an effective magnetic charge. Here,
by combining ultrasensitive optical detection of thermodynamic magnetization
noise with Monte Carlo simulations, we reveal emergent regimes of magnetic
excitations in artificial “tetris ice”. A marked increase of the intrinsic
noise at certain applied magnetic fields heralds the spontaneous proliferation
of fractionalized excitations, which can diffuse independently, without cost
in energy, along specific quasi-1D spin chains in the tetris ice lattice.

While the existence of elementary magnetic monopoles remains hypothetical,
spin (i.e., dipole) excitations in certain frustrated magnetic systems
Bramwell:2020 ; Skjaervo:2020 ; Rougemaille:2019 ; Nisoli:2021  can
fractionalize and separate into two delocalized “monopole-like” quasiparticles
that each carry an effective magnetic charge Castelnovo:2008 ; Mol:2009 ;
Morris:2009 ; Jaubert:2009 ; Ladak:2010 ; Castelnovo:2012 ; Perrin:2016 ;
Farhan:2019 . These fractionalized excitations are topologically protected,
can diffuse through the crystal lattice in thermal equilibrium, and can move
in response to applied magnetic fields, motivating studies of “magnetricity”
in analogy to electricity Bramwell:2009 ; Giblin:2011 ; Pollard:2012 ;
Vedmedenko:2016 ; Morley:2019 .

Mobile magnetic quasiparticles have been investigated in both natural spin ice
materials (such as the 3D pyrochlore Dy2Ti2O7 Bramwell:2009 ) and also in
engineered 2D arrays of nanomagnets known as artificial spin ice (ASI)
Ladak:2010 ; Mengotti:2011 ; Morgan:2011 ; Phatak:2011 . Recently, powerful
new detection modalities have emerged for revealing and studying these
fractionalized excitations, based on direct sensing of the very small – but
measurable – magnetization fluctuations that arise from their stochastic
creation, annihilation, and diffusion in thermal equilibrium. In Dy2Ti2O7,
spontaneous magnetization noise was detected at cryogenic temperatures via
SQUID magnetometry Dusad:2019 ; Samarakoon:2022 ; Hallen:2022 , from which the
presence of monopoles and their dynamics were inferred. Separately, in
archetypal square ASI lattices, magnetization noise from quasiparticle
kinetics was detected at room temperature via optical magnetometry Goryca:2021
; Goryca:2022 . Crucially, in these latter studies, the sudden appearance of
excess noise at certain applied magnetic fields revealed the presence of
phases rich in mobile magnetic charges.

Motivated by these studies, here we use optical magnetometry of spontaneous
noise to reveal new families of fractionalized excitations in the low-symmetry
frustrated ASI known as “tetris ice” Morrison:2013 ; Gilbert:2016 . By
applying small magnetic fields, tetris ice can be tuned through a variety of
complex spin configurations, whose degeneracies and boundaries are directly
revealed by noise. In particular, we find particularly intense and narrow
bands of noise for certain directions and ranges of applied field. Using Monte
Carlo simulations to deconstruct these noise signatures, these bands are shown
to herald novel regimes wherein magnetic quasiparticles delocalize and
proliferate, without cost in energy, along extended quasi-1D spin chains in
the lattice. These results demonstrate the power of noise-based studies to
probe microscopic details of complex magnetic phenomena, in this case
specifically the equilibrium kinetics driven by fractionalized magnetic
excitations.

Figure 1a shows a scanning electron microscope (SEM) image and illustration of
tetris ice. The individual islands have lateral dimensions 220 $\times$ 80 nm
and are made of ferromagnetic Ni0.8Fe0.2 using established lithographic
methods Gilbert:2016 . Each island behaves as a single Ising-like macrospin,
with magnetization orientation parallel or antiparallel to its long axis.
Crucially, the islands are very thin ($\approx$3.5 nm), so that at room
temperature they behave as thermally-active superparamagnets (i.e., their
magnetization direction thermally fluctuates in the absence of a strong
biasing magnetic field Kapaklis:2014 ; Chen:2019 ). This ensures that the ASI
can efficiently explore the huge manifold of possible magnetic configurations,
and remain at or near its lowest energy in thermal equilibrium.

Recent studies of tetris ice demonstrated that its low-energy configuration at
zero applied field ($B$=0) comprises ordered “backbones” (blue islands in Fig.
1a) that include the four-fold coordinated vertices ($z$=4), separated by
disordered “staircases” (red islands) containing $z$=3 and $z$=2 vertices
Gilbert:2016 . At $B$=0, the phenomenon of vertex frustration Morrison:2013
prevents all the staircase vertices from achieving their lowest-energy
configuration, leading to extensive degeneracy Gilbert:2016  and entropy-
driven ordering Saglam:2022 ; Miller:2022 .

Despite the intriguing physics observed in tetris ice at $B$=0, its properties
have never been explored in applied magnetic fields. However, given its lower
symmetry and the demonstrated emergence of complex collective behavior, a rich
set of phenomena in the presence of magnetic fields can be anticipated. We
investigate the field-dependent behavior of tetris ice by optically detecting
its intrinsic magnetization noise, using the approach shown in Fig. 1b. This
method is adapted from earlier studies of optically detected spin noise in
atomic, semiconductor, and ferromagnetic systems Crooker:2004 ; Crooker:2010 ;
Zapasskii:2013 ; Balk:2018 . Here, spontaneous magnetization fluctuations in
thermal equilibrium are passively detected via the Kerr rotation fluctuations
that they impart on a linearly-polarized laser reflected from the ASI surface.
These fluctuations are detected in real time by balanced photodiodes, and the
total noise power is computed via fast-Fourier transform methods. Small coils
apply magnetic fields $B_{x}$ and $B_{y}$ in the sample plane. By reorienting
the laser with respect to the sample, we measure fluctuations of either the
horizontal islands, $\langle(\delta m_{x})^{2}\rangle$, or the vertical
islands, $\langle(\delta m_{y})^{2}\rangle$.

Figures 1c and 1d show the main experimental results, which are field-
dependent maps of the spontaneous magnetization noise from the horizontal and
vertical islands in tetris ice. The maps reveal a complex structure, with both
sharp and diffuse bands of noise appearing along certain directions and
regions of $B_{x}$ and $B_{y}$. The dark areas, where noise is absent,
indicate regions of stable (non-fluctuating) magnetic order, which are
separated by noisy boundaries. As discussed below in Fig. 2, these boundaries
reveal where different magnetic configurations become energetically
degenerate. Note that these maps are not related by 90∘ rotation, reflecting
the absence of $C_{4}$ symmetry. Fig. 1d even shows an annular region of noise
surrounding $B\approx 0$ (but Fig. 1c does not), which is related to
fluctuations within the backbones and will be discussed later.

Of primary interest, and the main focus of this work, are the bright and
narrow bands in the noise map where spontaneous fluctuations are especially
strong (labeled A and B, indicated by dashed lines). These bands appear only
for sufficiently large applied fields, and follow specific diagonal
trajectories in the $(B_{x},B_{y})$ plane. Motivated by recent studies of
archetypal square ASI, wherein the sudden onset of spontaneous noise revealed
regimes of highly mobile monopoles, we therefore seek to identify the origin
of these intense fluctuations in tetris ice.

Figure 1: a) SEM image and diagram of tetris ice. The lattice comprises
‘backbones’ (blue islands) and ‘staircases’ (red islands). b) Experiment:
Magnetization fluctuations are passively detected via the polarization (Kerr
rotation) fluctuations that they impart on a linearly-polarized laser
reflected from the ASI. c) A map of the measured noise power from the
horizontal islands, $\langle(\delta m_{x})^{2}\rangle$, in thermal
equilibrium, vs. in-plane magnetic fields $B_{x}$ and $B_{y}$. d) Same, but
for the vertical islands, $\langle(\delta m_{y})^{2}\rangle$. Intense noise
along bands A and B herald the proliferation of fractionalized magnetic
excitations.

Crucial insight into these complex noise maps is provided by Monte Carlo (MC)
simulations of Glauber spin dynamics, following our earlier simulations of
noise in square and quadrupolar ASI Goryca:2021 ; Goryca:2022 . We simulated
32$\times$32 vertex tetris lattices with periodic boundary conditions, and
only nearest- and next-nearest neighbor interactions were considered ($J_{1}$
and $J_{2}$; see Fig. 1a). For each MC time step, $N$ single islands were
chosen randomly ($N$ is the number of islands in the lattice; cluster and loop
flips were not considered), and the spin flip acceptance probability was
$[1+\mathrm{exp}(\Delta/kT)]^{-1}$, where $\Delta$ is the energy difference
associated with the spin flip and $k$ is Boltzmann’s constant. We used
$J_{1}=1.7J_{2}$, and temperature $kT=0.7J_{2}$ that corresponds to room
temperature in our structures; both values are consistent with micromagnetic
MuMax3 simulations MuMax3  of these permalloy islands Goryca:2021 . These
values were also validated by directly comparing measured and simulated noise
maps. In particular, as discussed below, the positions of the various noise
bands is determined by $J_{1}$ and $J_{2}$, and their widths are determined in
part by $kT$. At each value of ($B_{x}$, $B_{y}$), $10^{6}$ MC annealing steps
were performed to ensure thermal equilibrium, and then the magnetization was
recorded for several million additional MC time steps. The average
magnetization $M$ and the thermodynamic fluctuations $\langle[\delta
M(t)]^{2}\rangle$ about this mean value were determined from the computed time
series.

Figure 2 shows a calculated map of the field-dependent magnetization,
revealing a rich tableau of magnetic configurations (color and brightness
indicate the direction and magnitude of M). Dotted black lines denote
boundaries between the different static configurations that are stabilized at
high field when $B_{x}$ and/or $B_{y}$ are large (labelled i-xii, and depicted
around the periphery of the map). Beginning in the upper-right corner of the
map where both $B_{x}$ and $B_{y}$ are large and positive (region i), all
horizontal and vertical islands are trivially magnetized along $+\hat{x}$ and
$+\hat{y}$, respectively, by the applied field. Proceeding clockwise around
the map (i.e., as $B_{y}\rightarrow 0$), a subset of the vertical islands
flips to $-\hat{y}$ when $B_{y}<J_{1}/\mu$, where $\mu$ is the moment of an
island; that is, when the Zeeman energy balances the interaction energy, and
the lattice transitions to configuration ii. The pink arrows in the diagram
next to region $ii$ show which islands have flipped upon transitioning from
$i\rightarrow ii$: here it is only the subset of vertical islands in the $z$=4
vertices; these have an ‘unbalanced’ number of horizontal neighbors.
Continuing clockwise around the perimeter, from region $ii\rightarrow
iii\rightarrow iv\rightarrow v~{}...$, different subsets of islands (pink
arrows) flip when crossing each boundary.

Figure 2: Magnetic states of tetris ice, vs. $B_{x}$ and $B_{y}$. The color
and brightness indicate the direction and magnitude of the average
magnetization M, as computed by MC simulations. Dotted lines denote boundaries
between different magnetic configurations (labeled i-xii) that are stabilized
when $|B_{x,y}|$ is large. Moment configurations are shown; pink arrows
indicate the subset of islands that flip when transitioning from state
$i\rightarrow ii\rightarrow iii$ …, clockwise around the map’s perimeter.
Boundaries occur when Zeeman energies balance the local interaction energies
given by $J_{1,2}$. The dark central region is discussed in the Supplemental
Material.

From this diagram it is already possible to associate the experimentally-
measured bands of noise with boundaries between stable magnetic
configurations, and more importantly that the bands A and B of exceptionally
high noise occur at the internal diagonal $iii\leftrightarrow v$ boundary (or
equivalently, $ix\leftrightarrow xi$), and at the diagonal $xii\leftrightarrow
ii$ ($vi\leftrightarrow viii$) boundary.

Figure 3: Calculated maps of the equilibrium magnetization noise from a) the
horizontal and b) the vertical islands, showing excellent agreement with
experiment (cf. Fig. 1). c-e) Deconstructing noise: Noise maps from specific
subsets of horizontal islands, indicated by bright colors (blue and red denote
backbone and staircase islands). f-h) Same, but for subsets of the vertical
islands. i, j) Blue arrows indicate the 1D spin chains that flip when crossing
the diagonal $iii\leftrightarrow v$ and $vi\leftrightarrow viii$ boundaries,
respectively. These chains form the “spines” and “ribs” of the backbones.
Along these boundaries, the balance of Zeeman and interaction energies
restores the Z2 symmetry of the spin chains.

This association is confirmed in Figs. 3a and 3b, which show the calculated
noise power from the horizontal and vertical islands [$(\delta m_{x})^{2}$ and
$(\delta m_{y})^{2}$, respectively]. Crucially, the overall agreement with the
experimental data is remarkably good (cf. Fig. 1), including capturing both
narrow and diffuse bands of noise, as well as the small annulus of noise
surrounding $B\approx 0$ in $(\delta m_{y})^{2}$. Because the MC simulations
are validated in this way, we have confidence that they can be explored in
detail to reveal the correlated physics of the system. Most importantly, the
simulations allow us to determine which subset(s) of islands give rise to the
different noise signatures that are observed experimentally. In other words,
we can use the validated MC simulations to deconstruct the noise maps and
tease apart the fluctuations of individual moments in ways that are
experimentally intractable.

To this end, Figs. 3(c-e) show the calculated noise from specific subsets of
the horizontal islands (indicated by bright blue and red in the adjacent
diagrams). Similarly, Figs. 3(f-h) show noise from subsets of vertical
islands. From these noise maps it is clear that the horizontal and vertical
bands of noise appearing at large $B_{x,y}$ originate from subsets of
disconnected and uncorrelated islands that become thermally active when
crossing a horizontal or vertical configuration boundary in Fig. 2.

Figure 4: a) A thermal spin-flip of one island (green) in a ordered 1D “spine”
generates a dipole excitation and two adjacent vertices with magnetic charge
$q=\pm 2$. b) Subsequent flips of neighboring islands fractionalize the
excitation into two delocalized quasiparticles that move along the spine,
leaving behind a string of flipped islands. For $(B_{x},B_{y})$ along the
$iii\leftrightarrow v$ boundary this motion costs no net energy, and the
charges therefore diffuse independently and proliferate, leading to the excess
noise observed in band A. If $(B_{x},B_{y})$ is tuned off the boundary, the
$\pm$2 charges move preferentially in opposite directions and annihilate. c,d)
Similar phenomena along the “ribs”, at the $vi\leftrightarrow viii$ boundary
(band B). Note that ordered ribs end in $z$=3 vertices, with $q=\pm$1;
therefore excitations have $q=\mp$1 or $\pm 2$. e) Map of the density of
fractionalized (delocalized) quasiparticles.

Most importantly, these deconstructed maps reveal the origins of the intense
noise that emerges in the diagonal bands A and B of the experimental data.
Figures 3c and 3f, considered together, show that band A arises from
fluctuations of the horizontal and vertical islands comprising the extended 1D
spin chains that connect the $z=4$ vertices (i.e., the “spines” of the
backbones). As Fig. 3i shows, it is precisely and only these connected spin
chains that flip at the diagonal $iii\leftrightarrow v$ boundary (and
$ix\leftrightarrow xi$, by symmetry). Along this border, configurations iii
and v are energetically degenerate and the balance of Zeeman and interaction
energies restores the Z2 spin symmetry of the entire chain. Similarly, Figs.
3d and 3g show that noise in band B arises from the four-island chains that
form the “ribs” of the backbones; only these spin chains are degenerate at the
diagonal $ii\leftrightarrow xii$ (and $vi\leftrightarrow viii$) boundary (see
Fig. 3j).

Figures 4a,b show how fluctuations at the $iii\leftrightarrow v$ boundary
arise from spontaneous creation, fractionalization, and kinetics of magnetic
quasiparticles along the spin chains that form the spines of the backbones.
Vertices along these chains have coordination z=2 or z=4, and an effective
magnetic charge, $q$, given by the balance of inward- and outward-facing
moments. In region iii, all islands in the chain have stable orientation up or
right, and all vertices have $q=0$ (2-in, 2-out at the $z$=4 vertices, and
1-in, 1-out at the $z$=2 vertices). Just across the border, in region v, all
islands in the chain have orientation down or left, and again all vertices
have $q=0$. Along the border, however, the two configurations are degenerate
and Z2 symmetry of the spin chain is restored. Beginning with an ordered
chain, a thermal spin-flip of any island in the chain creates a dipole
excitation and a pair of adjacent vertices with $q=\pm 2$ (e.g., a $z$=4
vertex with 3-in, 1-out and a neighboring $z$=2 vertex with 0-in, 2-out).
Subsequent flips of neighboring islands cause the initial dipole excitation to
fractionalize into effective magnetic charges that separate and move along the
1D chain, leaving behind a string of flipped islands. Crucially, exactly along
the $iii\leftrightarrow v$ border this motion has no net cost in energy, and
spin excitations can therefore readily fractionalize into quasiparticles that
diffuse independently along the spin chain. The mean quasiparticle lifetime
becomes very long, resulting in their rapid proliferation and formation of a
regime rich in fractionalized excitations. Since the kinetics of these
quasiparticles is necessarily associated with flipping of islands, this
generates a telltale spike in the equilibrium noise, precisely as observed.

Similarly, Figs. 4c,d shows that noise in band B arises from the emergence and
kinetics of mobile charges along the short 4-island ribs of the backbones.
Figure 4e shows a map of the calculated density of fractionalized excitations
(i.e., quasiparticles separated by more than one lattice spacing). Bands A and
B of the measured noise maps therefore correspond to new regimes of mobile
spin excitations, stabilized by the interplay between Zeeman and interaction
energies.

Finally, note that the deconstructed noise maps in Figs. 3e,h confirm that
fluctuations at $B$=0 arise solely from the staircases, which at zero field
are vertex-frustrated, disordered, and extensively degenerate Morrison:2013 ;
Gilbert:2016 . In contrast, Figs. 3(c,d,f,g) confirm that the backbones do not
fluctuate at $B$=0, in line with the known stable type-I order of the $z$=4
vertices Morrison:2013 ; Gilbert:2016 . However, the small annular region of
noise that arises from the backbone spins reveals where the $z$=4 vertices
transition from type-I to stable type-II order as $|B_{x,y}|$ increases (see
Supplemental Material). As a final point of interest, we note that although
the horizontal backbone islands in Fig. 3c and 3d both generate annular
regions of noise, the total noise from all horizontal islands does not (see
Fig. 3a). This is because, as discussed in the Supplemental Material,
fluctuations of these island subsets are temporally anticorrelated, giving
minimal contribution to the net $(\delta m_{x})^{2}$, in agreement with
experiment. In contrast, fluctuations of the different vertical spins in Figs.
3f,g are temporally correlated, giving a large total $(\delta m_{y})^{2}$ in
Fig. 3b, again in agreement with experiment and further demonstrating the
utility of noise decomposition in ASIs.

We have shown that tetris ice can host a rich variety of field-induced
magnetic states. Comparing the configuration diagram (Fig. 2) and the
experimental noise maps (Fig. 1) demonstrates that noise studies are a potent
tool to directly detect and disentangle these states, where the critical
boundaries are characterized by high kinetic activity. Crucially, certain
noise signatures herald collective behaviors corresponding to the diffusion
and proliferation of fractionalized magnetic quasiparticles along specific
spin chains in the lattice. The fact that these phenomena are field-tunable
may open the door to future functionalities in sensing or beyond-Turing
computation Arava:2018 ; Caravelli:2020 ; Heyderman:2022 ; Gartside:2022 .
Future studies will concentrate on the detailed nature of the noise spectra,
and establish relations between the non-triviality of the disorder of certain
phases, tied to their collective behavior, and the ‘color’ of the noise
spectrum Goryca:2021 ; Hallen:2022 ; Klyuev:2017 ; Nisoli:EPL2021  at larger
frequencies.

We acknowledge support from the Los Alamos LDRD program, the US Department of
Energy (DOE) Quantum Science Center, the National Science Foundation (NSF)
DMR-2103711, and the US DOE Office of Basic Energy Sciences, Materials
Sciences and Engineering Division under Grant No. DE-SC0020162. The NHMFL is
supported by NSF DMR-1644779, the State of Florida, and the US DOE. Research
at the University of Warsaw leading to these results has received funding from
the Norwegian Financial Mechanism 2014-2021 under Grant No.
2020/37/K/ST3/03656 and from the Polish National Agency for Academic Exchange
within Polish Returns program under Grant No. PPN/PPO/2020/1/00030.

## References

  * (1) S. T. Bramwell and M. J. Harris, The history of spin ice, J. Phys.: Condens. Matter 32, 374010 (2020).
  * (2) S. H. Skjaervø, C. H. Marrows, R. L. Stamps, and L. J. Heyderman, Advances in artificial spin ice, Nat. Rev. Phys. 2, 13 (2020).
  * (3) N. Rougemaille and B. Canals, Cooperative magnetic phenomena in artificial spin systems: spin liquids, Coulomb phase and fragmentation of magnetism, Eur. Phys. J. B 92, 62 (2019).
  * (4) P. Schiffer and C. Nisoli, Artificial spin ice: Paths forward, Appl. Phys. Lett. 118, 110501 (2021).
  * (5) C. Castelnovo, R. Moessner, and S. L. Sondhi, Magnetic monopoles in spin ice, Nature 451, 42-45 (2008).
  * (6) L. A. Mól, R. L. Silva, R. C. Silva, A. R. Pereira, W. A. Moura-Melo, and B. V. Costa, Magnetic monopole and string excitations in two-dimensional spin ice, J. Appl. Phys. 106, 063913 (2009).
  * (7) D. J. P. Morris et al., Dirac Strings and Magnetic Monopoles in the Spin Ice Dy2Ti2O7, Science 326, 411 (2009).
  * (8) L. D. C. Jaubert and P. C. W. Holdsworth, Signature of magnetic monopole and Dirac string dynamics in spin ice, Nat. Phys. 5, 259 (2009).
  * (9) S. Ladak, D. E. Read, G. K. Perkins, L. F. Cohen, and W. R. Branford, Direct observation of magnetic monopole defects in an artificial spin-ice system, Nat. Phys. 6, 359 (2010).
  * (10) C. Castelnovo, R. Moessner, and S. L. Sondhi, Spin Ice, Fractionalization, and Topological Order, Annu. Rev. Condens. Matter Phys. 3, 35 (2012).
  * (11) Y. Perrin, B. Canals, and N. Rougemaille, Extensive degeneracy, Coulomb phase and magnetic monopoles in artificial square ice, Nature 540, 410-413 (2016).
  * (12) A. Farhan, M. Saccone, C. F. Petersen, S. Dhuey, R. V. Chopdekar, Y.-L. Huang, N. Kent, Z. Chen, M. J. Alava, T. Lippert, A. Scholl, and S. van Dijken, Emergent magnetic monopole dynamics in macroscopically degenerate artificial spin ice, Sci. Adv. 5, eaav6380 (2019).
  * (13) S. T. Bramwell, S. R. Giblin, S. Calder, R. Aldus, D. Prabhakaran, and T. Fennell, Measurement of the charge and current of magnetic monopoles in spin ice, Nature 461, 956-959 (2009).
  * (14) S. R. Giblin, S. T. Bramwell, P. C.W. Holdsworth, D. Prabhakaran, and I. Terry, Creation and measurement of long-lived magnetic monopole currents in spin ice, Nat. Phys. 7, 252 (2011).
  * (15) S. D. Pollard, V. Volkov, and Y. Zhu, Propagation of magnetic charge monopoles and Dirac flux strings in an artificial spin-ice lattice, Phys. Rev. B 85, 180402 (2012).
  * (16) E. Y. Vedmedenko, Dynamics of bound monopoles in artificial spin ice: how to store energy in Dirac strings, Phys. Rev. Lett. 116, 077202 (2016).
  * (17) S. A. Morley, J. M. Porro, A. Hrabec, M. C. Rosamond, D. A. Venero, E. H. Linfield, G. Burnell, M.-Y. Im, P. Fischer, S. Langridge, and C. H. Marrows, Thermally and field-driven mobility of emergent magnetic charges in square artificial spin ice, Sci. Rep. 9, 15989 (2019).
  * (18) E. Mengotti, L. J. Heyderman, A. F. Rodriguez, F. Nolting, R. V. Hugli, and H.-B. Braun, Real-space observation of emergent magnetic monopoles and associated Dirac strings in artificial kagome spin ice, Nat. Phys. 7, 68 (2011).
  * (19) J. P. Morgan, A. Stein, S. Langridge, and C. H. Marrows, Thermal ground state ordering and elementary excitations in artificial magnetic square ice, Nat. Phys. 7, 75 (2011).
  * (20) C. Phatak, A. K. Petford-Long, O. Heinonen, M. Tanase, and M. De Graef, Nanoscale structure of the magnetic induction at monopole defects in artificial spin-ice lattices, Phys. Rev. B 83, 174431 (2011).
  * (21) R. Dusad, F. K. K. Kirschner, J. C. Hoke, B. R. Roberts, A. Eyal, F. Flicker, G. M. Luke, S. J. Blundell, and J. C. S. Davis, Magnetic monopole noise, Nature 571, 234-239 (2019).
  * (22) A. M. Samarakoon et al., Anomalous magnetic noise in an imperfectly flat landscape in the topological magnet Dy2Ti2O7, PNAS 119, e2117453119 (2022).
  * (23) J. N. Hallén, S. A. Grigera, D. A. Tennant, C. Castelnovo, R. Moessner, Dynamical fractal and anomalous noise in a clean magnetic crystal, Science 378, 1218-1221 (2022).
  * (24) M. Goryca, X. Zhang, J. Li, A. L. Balk, J. D. Watts, C. Leighton, C. Nisoli, P. Schiffer, and S. A. Crooker, Field-Induced Magnetic Monopole Plasma in Artificial Spin Ice, Phys. Rev. X 11, 011042 (2021).
  * (25) M. Goryca, X. Zhang, J. D. Watts, C. Nisoli, C. Leighton, P. Schiffer, and S. A. Crooker, Magnetic field dependent thermodynamic properties of square and quadrupolar artificial spin ice, Phys. Rev. B 105, 094406 (2022).
  * (26) M. J. Morrison, T. R. Nelson, and C. Nisoli, Unhappy vertices in artificial spin ice: new degeneracies from vertex frustration, New J. Phys. 15, 045009 (2013).
  * (27) I. Gilbert, Y. Lao, I. Carrasquillo, L. O’Brien, J. D. Watts, M. Manno, C. Leighton, A. Scholl, C. Nisoli, and P. Schiffer, Emergent reduced dimensionality by vertex frustration in artificial spin ice, Nat. Phys. 12, 162 (2016).
  * (28) V. Kapaklis, U. B. Arnalds, A. Farhan, R. V. Chopdekar, A. Balan, A. Scholl, L. J. Heyderman, and B. Hjörvarsson, Thermal fluctuations in artificial spin ice, Nat. Nanotech. 9, 514 (2014).
  * (29) X. M. Chen, B. Farmer, J. S. Woods, S. Dhuey, W. Hu, C. Mazzoli, S. B. Wilkins, R. V. Chopdekar, A. Scholl, I. K. Robinson, L. E. De Long, S. Roy, and J. T. Hastings, Spontaneous Magnetic Superdomain Wall Fluctuations in an Artificial Antiferromagnet, Phys. Rev. Lett. 123, 197202 (2019).
  * (30) H. Saglam et al., Entropy-driven order in an array of nanomagnets, Nat. Phys. 18, 706-712 (2022).
  * (31) J. L. Miller, Entropy and order work together in an artificial spin ice, Phys. Today 75, 17 (2022).
  * (32) S. A. Crooker, D. G. Rickel, A. V. Balatsky, and D. L. Smith, Spectroscopy of spontaneous spin noise as a probe of spin dynamics and magnetic resonance, Nature 431, 49 (2004).
  * (33) S. A. Crooker, J. Brandt, C. Sandfort, A. Greilich, D. R. Yakovlev, D. Reuter, A. D. Wieck, and M. Bayer, Spin noise of electrons and holes in self-assembled (In,Ga)As quantum dots, Phys. Rev. Lett. 104, 036601 (2010).
  * (34) V. S. Zapasskii, A. Greilich, S. A. Crooker, Yan Li, G. G. Kozlov, D. R. Yakovlev, D. Reuter, A. D. Wieck, M. Bayer, Optical spectroscopy of spin noise, Phys. Rev. Lett. 110, 176601 (2013).
  * (35) A. L. Balk, F. Li, I. Gilbert, J. Unguris, N. A. Sinitsyn, and S. A. Crooker, Broadband Spectroscopy of Thermodynamic Magnetization Fluctuations through a Ferromagnetic Spin-Reorientation Transition, Phys. Rev. X 8, 031078 (2018).
  * (36) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, The Design and Verification of MuMax3, AIP Adv. 4, 107133 (2014).
  * (37) H. Arava, P. M. Derlet, J. Vijayakumar, J. Cui, N. S. Bingham, A. Kleibert, and L. J. Heyderman, Computational logic with square rings of nanomagnets, Nanotechnology 29, 265205 (2018).
  * (38) F. Caravelli and C. Nisoli, Logical gates embedding in artificial spin ice, New J. Phys. 22, 103052 (2020).
  * (39) L. J. Heyderman, Spin ice devices from nanomagnets, Nat. Nanotech. 17, 435 (2022).
  * (40) J. Gartside et al., Reconfigurable training and reservoir computing in an artificial spin-vortex ice via spin-wave fingerprinting, Nat. Nanotech. 17, 460 (2022).
  * (41) A. V. Klyuev, M. I. Ryzhkin, A. V. Yakimov, Statistics of Fluctuations of Magnetic Monopole Concentration in Spin Ice, Fluct. and Noise Lett. 16, 1750035 (2017).
  * (42) C. Nisoli, The color of magnetic monopole noise, Europhys. Lett. 135, 57002 (2021).

##  I Supplemental Material

###  I.1 Configurations of tetris ice in small magnetic fields near B
$\approx$ 0

Figure S1a shows a calculated zero-temperature magnetization map of tetris
ice, focusing on the region of small applied in-plane magnetic fields $B_{x}$
and $B_{y}$. Dotted lines denote the boundaries between different stable
magnetic configurations. For easier comprehension of the map, the net
magnetization can be decomposed into the magnetization of staircase and
backbone sublattices, shown in Figs. S1b and S1c, respectively, for which
moment configurations are also explicitly shown. For the staircase sublattice,
the four regions in the corners of the map correspond to fully polarized
configurations, where all islands are magnetized most closely along the
direction of the external magnetic field. The two central regions closest to
$B=0$ reflect the situation where the phenomenon of vertex frustration
prevents all the staircase vertices from achieving their lowest-energy state
and half of the $z=3$ vertices remain “unhappy”. In the case of the backbone
sublattice (Fig. S1c), the central grey region indicates where all $z$=4
vertices have type-I configuration (two spins in, two spins out, with no net
magnetization), which is surrounded by four stable phases where all $z$=4
vertices have type-II configuration (two spins in, two spins out, but with a
net diagonal magnetization). Noise along the annular perimeter of this central
region results from stochastic fluctuations within the backbones between
type-I and type-II order.

Figure S1: a) Calculated zero-temperature magnetization map of tetris ice. The
color and brightness indicate the direction and magnitude of the magnetization
M. Dotted black lines denote the boundaries between different stable magnetic
configurations. The net magnetization can be decomposed into the magnetization
of staircase and backbone sublattices, shown in panels b) and c),
respectively, for which moment configurations are explicitly shown. Inside the
central grey region of panel c), the $z=4$ vertices in the backbones exhibit
stable type-I order. Outside of this region they exhibit type-II order. The
red and blue dotted lines denote where increased noise is observed due to the
proliferation of fractionalized magnetic excitations, discussed in the main
text. The red lines follow $B_{y}=-2B_{x}$. The blue lines follow
$B_{y}=B_{x}-\left(J_{1}-J_{2}\right)/\mu$ for $B_{x}>0$, and
$B_{y}=B_{x}+\left(J_{1}-J_{2}\right)/\mu$ for $B_{x}<0$.

###  I.2 Mapping out temporal correlations between the noise from different
subsets of islands

Figure S2 shows an example of the detailed insight that Monte Carlo (MC)
simulations can provide into the temporal correlations between the noise from
different subsets of islands within a magnetic lattice. Figure S2a shows
calculated correlations between fluctuations of different subsets of
horizontal islands within the backbones (pink and green islands in the upper
panel) vs. magnetic field. The correlations in the map shown in the bottom
panel are calculated as $C=\left<\delta m_{\mathrm{P}}\delta
m_{\mathrm{G}}\right>$, where $\delta m_{\mathrm{P}}$ and $\delta
m_{\mathrm{G}}$ denote fluctuations of the pink and green island subsets,
respectively. Red colors indicate positive temporal correlation; i.e.,
stochastic spin-flips of (say) pink islands are accompanied by flips of the
green islands in the same direction within the same MC step. Blue colors
indicate temporal _anti_ correlation; i.e., flips of pink islands in a given
direction are accompanied by flips of green islands in the opposite direction,
within the same MC step.

Figure S2a shows correlations between two different subsets of horizontal
islands within tetris’ backbones. The pink and green islands are those shown
in Figs. 3c and 3d of the main text (respectively), where it was shown that
each subset, considered alone, generates an annular region of noise
surrounding $B\approx 0$. However, in combination they do not (see Fig. 3a).
This is because fluctuations of the pink and green islands are strongly
anticorrelated along the upper and lower horizontal borders of the annular
region. That is, when (say) a pink island flips, this rapidly drives a likely
flip of the nearby green island in the opposite direction (because they are
coupled via the connecting vertical island), leading to no net change in total
magnetization in the horizontal direction.

Conversely, Fig. S2b shows that the two indicated subsets of _vertical_
islands within the backbones generate noise that is strongly positively
correlated along the upper and lower borders of the annular region. As a
result, these horizontal sections of annular noise are clearly visible in the
total noise from the vertical islands shown in Fig. 3b, in agreement with the
experimental data. Interestingly, even though noise from these pink and green
subsets in Fig. S2b is _anti_ correlated along the approximately vertical
borders of the annular region, those sections are still visible in the total
noise shown in Fig. 3b. This is because the number of the pink islands in Fig.
S2b is twice the number of the green islands (which is not the case for the
horizontal islands in Fig. S2a). As a result, the noise generated by the two
island subsets does not cancel out completely, and the full annular region of
noise is still visible in Fig. 3b. This provides a fine example of how the
deconstruction of the noise into contributions from, and correlations between,
different subsets of moments is a powerful tool in the analysis of the
kinetics and coupling of complex magnets.

Figure S2: a) Calculated temporal correlations $C$ between fluctuations of the
pink and green subsets of horizontal islands (both subsets are components of
tetris’ backbones), vs. magnetic field. Red colors denote positive correlation
(flips of pink and green islands tend to occur in the same direction, within a
MC step), while blue colors denote _anti_ correlation (pink and green islands
tend to flip in opposite directions, within a MC step). Strong
anticorrelations along the upper and lower boundaries of the central annular
region indicates that noise from these two subsets largely cancels out, and
therefore does not appear in the total noise from all the horizontal islands
(see Fig. 3a of the main text). b) Same, but for temporal correlations between
fluctuations from different subsets of vertical islands within the backbones.


# Density of states and spectral function of a superconductor out of a
quantum-critical metal

Shang-Shun Zhang  School of Physics and Astronomy and William I. Fine
Theoretical Physics Institute, University of Minnesota, Minneapolis, MN 55455,
USA  Andrey V. Chubukov  School of Physics and Astronomy and William I. Fine
Theoretical Physics Institute, University of Minnesota, Minneapolis, MN 55455,
USA

###### Abstract

We analyze the validity of a quasiparticle description of a superconducting
state at a metallic quantum-critical point (QCP). A normal state at a QCP is a
non-Fermi liquid with no coherent quasiparticles. A superconducting order gaps
out low-energy excitations, except for a sliver of states for non-s-wave gap
symmetry, and at a first glance, should restore a coherent quasiparticle
behavior. We argue that this does not necessarily hold as in some cases the
fermionic self-energy remains singular slightly above the gap edge. This
singularity gives rise to markedly non-BCS behavior of the density of states
and to broadening and eventual vanishing of the quasiparticle peak in the
spectral function. We analyze the set of quantum-critical models with an
effective dynamical 4-fermion interaction, mediated by a gapless boson at a
QCP, $V(\Omega)\propto 1/\Omega^{\gamma}$. We show that coherent quasiparticle
behavior in a superconducting state holds for $\gamma<1/2$, but breaks down
for larger $\gamma$. We discuss signatures of quasiparticle breakdown and
compare our results with the data.

Introduction. Metals near a quantum critical point (QCP) display a number of
non-Fermi liquid properties like linear-in-$T$ resistivity, a broad peak in
the spectral function near $k_{F}$ with linear-in-$\omega$ width, singular
behavior of optical conductivity, etc [1, 2, 3, 4, 5, 6, 7, 8, 9, 10,
*jin2011link, 12, 13, 14, 15, 16, 17, 18, 19]. These properties are often
thought to be caused by the coupling of fermions to near-gapless fluctuations
of an order parameter, which condenses at a QCP [20, 21, 22, *finger_2001, 24,
25, 26, 27, 28, 29, 30, 31]. The same fermion-boson interaction gives rise to
superconductivity near a QCP [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, *son2,
43, 44, 45].

A superconducting order gaps out low-energy excitations, leaving at most a
tiny subset of gapless states for a non-$s-$wave order parameter. A general
belief has been that this restores fermionic coherence. A frequently cited
experimental evidence is the observed re-emergence of a quasiparticle peak
below $T_{c}$ in near-optimally doped cuprates (see e.g., Ref. [46]). From
theory side, the argument is that the fermionic self-energy in a
superconductor has a conventional Fermi-liquid form $\Sigma(\omega)\sim\omega$
at the lowest $\omega$, in distinction from a non-Fermi-liquid
$\Sigma(\omega)\propto\omega^{a}$ with $a<1$ in the normal state [47, 48,
*DellAnna2006, 50, 51, *metlitski2010quantum2, 53, 54, *Fitzpatrick_15, 56,
57, *sslee2, *lunts_2017, 60, 61]. In this paper, we analyze theoretically
whether fermions in a superconducting state at a QCP can be viewed as well-
defined coherent quasiparticles. We argue that this is not necessarily the
case as fermionic self-energy can still be singular on a real frequency axis
immediately above the gap edge. This singularity gives rise to markedly non-
BCS behavior of the density of states (DoS) and to broadening and eventual
vanishing of the quasiparticle peak.

Figure 1: Three possible forms of the electronic spectral function
$A({\bm{k}},\omega)$ at $T=0$ in a quantum critical superconductor at a small
but finite $k-k_{F}$ and in the absence of impurity broadening. (a):
$A({\bm{k}},\omega)$ vanishes at $\rvert\omega\rvert=\Delta$ and has a well-
defined peak at $\omega>\Delta$, (b): $A({\bm{k}},\omega)$ diverges at
$\rvert\omega\rvert=\Delta$, but it non-monotonic at larger $\omega$. The peak
in $A({\bm{k}},\omega)$ at $\rvert\omega\rvert>\Delta$ broadens, but still
exists. (c): $A({\bm{k}},\omega)$ diverges at $\rvert\omega\rvert=\Delta$, and
monotonically decreases at larger $\omega$. In case (a) fermions can be viewed
as well-defined quasiparticles, in case (c) the quasiparticle picture
completely breaks down. The case (b) is the intermediate one between (a) and
(c).  Figure 2:  (a) Exponents $\nu$ and $c$ for the leading and the
subleading terms in the expansion $D(\omega)\simeq
1+\alpha(\Delta-\omega)^{\nu}+\beta(\Delta-\omega)^{\nu+c}$, where
$D(\omega)=\Delta(\omega)/\omega$ and the gap edge $\Delta$ is the solution of
$D(\omega=\Delta)=1$. (b) Numerical result for $D(\omega)$ for $\gamma=0.8$.
Inset shows the power-law behavior near the gap edge with $\nu=1.18$,
consistent with (a). (c) Fermionic DoS at $T=0$ for $\gamma=0.35$ (thick green
line) and $\gamma=0.8$ (thin pink line). In both cases, the DoS vanishes below
the gap edge $\Delta$ and has a power-law singularity above it
$N(\omega)\propto 1/(\omega-\Delta)^{\nu/2}$, but the exponent $\nu$ is
different in the two cases, as we show in the right panel.

For superconductivity away from a QCP, mediated by a massive boson, numerous
earlier studies have found that the spectral function $A({\bm{k}},\omega)$ at
$T=0$ has a $\delta$-functional peak at
$\omega=(\Delta^{2}+(\xi_{\bm{k}}/Z)^{2})^{1/2}$, where
${\xi_{\bm{k}}}=v_{F}(k-k_{F})$ is a fermionic dispersion ($v_{F}$ is a Fermi
velocity), $\Delta$ is a superconducting gap, and $Z$ is an inverse
quasiparticle residue. A $\delta$-functional peak holds for momenta near the
Fermi surface, as long as $\omega<\Delta+\omega_{0}$, where $\omega_{0}$ is a
mass of a pairing boson in energy units. At larger $\omega$, fermionic damping
kicks in, and the peak broadens. The same physics leads to peak-dip-hump
behavior of $A({\bm{k}},\omega)$ as a function of $\omega$, observed most
spectacularly in near-optimally doped cuprate Bi2Sr2CaCu2O8+δ (see, e.g, Refs.
[62, 63]). At a QCP, the pairing boson becomes massless and $\omega_{0}$
vanishes. This creates a singular behavior near the gap edge at
$\omega=\Delta$, which holds even when ${\xi_{\bm{k}}}$ is finite.

A simple experimentation shows that there are three possible forms of
$A({\bm{k}},\omega)$, which we present in Fig. 1: it (i) either vanishes at
$\omega=\Delta$ and has a well-defined peak at $\omega>\Delta$ whose width at
small ${\xi_{\bm{k}}}$ is parametrically smaller than its energy; or (ii)
diverges at $\omega=\Delta$, but is non-monotonic at larger $\omega$ and
displays a broad maximum at some $\omega>\Delta$, or (iii) diverges at
$\omega=\Delta$ and monotonically decreases at larger $\omega$. In the first
case, fermions in a quantum-critical superconductor can be viewed as well-
defined quasiparticles; in the last case the quasiparticle picture completely
breaks down; the second case is the intermediate one between the other two.
Our goal is to understand under what circumstances $A({\bm{k}},\omega)$ of a
quantum-critical superconductor has one of these forms.

Model. For our study, we consider dispersion-full fermions, Yukawa-coupled to
a massless boson. We assume, like in earlier works (see, e.g., Refs. [64]),
that a boson is Landau overdamped, and its effective velocity is far smaller
than $v_{F}$. In this situation, the interaction that gives rise to non-Fermi
liquid in the normal state and to superconductivity, is a purely dynamical
$V(\Omega)$. The fermionic self-energy and the pairing gap, tuned into a
proper spatial pairing channel, are then determined by two coupled equations
in the frequency domain. At a QCP, $V(\Omega)$ is singular at vanishing
$\Omega$ in spatial dimension $D\leq 3$, and behaves as
$V(\Omega)\propto({\bar{g}}/\Omega)^{\gamma}$, where ${\bar{g}}$ is the
effective fermion-boson coupling, and the exponent $\gamma$ is determined by
the underlying microscopic model. The most studied models of this kind are of
fermions near an Ising-nematic or Ising/ferromagnetic QCP ($\gamma=1/3$) and
near an antiferromagnetic or charge density wave QCP $(\gamma=1/2)$. The same
effective interaction emerges for dispersion-less fermions in a quantum dot
coupled to Einstein bosons (the Yuakawa-SYK model) [65, *Schmalian_19a, 67,
68, 69]. For this last case, the exponent $\gamma$ is a continuous variable
$\gamma\in(0,1)$, depending on the ratio of fermion and boson flavors. An
extension of the Yukawa-SYK model to $\gamma\in(1,2)$ has recently been
proposed [70]. We follow these works and consider $\gamma$ as a continuous
variable. We note that the value of $\gamma$ is generally larger deep in a
superconducting state because of feedback from superconductivity on the
bosonic polarization. For simplicity, we neglect potential in-gap states
associated with non-$s$-wave pairing symmetry and focus on the spectral
function of fermions away from the nodal points and on features in the density
of states (DoS) above the gap edge. An extension to models with in-gap states
is straightforward.

In previous studies of the $\gamma$-model, we focused on the novel
superconducting behavior at $\gamma>1$, when the pairing interaction is
attractive on the Matsubara axis, while on the real axis Re$V(\Omega)$ is
repulsive [71, 72]. We argued that this dichotomy gives rise to phase slips of
the gap function on the real axis. Here, we restrict ourselves to $\gamma\leq
1$, when this physics is not present and, hence, does not interfere with the
analysis of the validity of a quasiparticle description in a superconducting
state.

Pairing gap and quasiparticle residue. For superconductivity mediated by a
dynamical interaction, the paring gap $\Delta(\omega)$ and the inverse
quasiparticle residue $Z(\omega)$ are functions of the running real fermionic
frequency $\omega$. We define the gap edge $\Delta$ (often called the gap)
from the condition $\Delta(\omega)=\omega$ at $\omega=\Delta$.

For our purposes, it is convenient to introduce
$D(\omega)=\Delta(\omega)/\omega$. The gap edge is at $|D|=1$. The equation
for $D(\omega)$ that we need to solve is

$\displaystyle\omega B(\omega)D(\omega)=A(\omega)+C(\omega),$ (1)

where $B(\omega)$ and $A(\omega)$ are regular functions of $\omega$ (see [73,
74]). The $C(\omega)$ term depends on the running $D(\omega)$,

$C(\omega)=\bar{g}^{\gamma}\sin\frac{\pi\gamma}{2}\int_{0}^{\omega}\frac{d\Omega}{\Omega^{\gamma}}\frac{D(\omega-\Omega)-D(\omega)}{\sqrt{D^{2}(\omega-\Omega)-1}}.$
(2)

Its presence makes Eq. (S3) an integral equation. The inverse residue
$Z(\omega)$ is expressed via $D(\omega^{\prime})$ as

$\displaystyle
Z(\omega)=B(\omega)+\frac{\bar{g}^{\gamma}\sin\frac{\pi\gamma}{2}}{\omega}\int_{0}^{\omega}\frac{d\Omega}{\Omega^{\gamma}}\frac{1}{\sqrt{D^{2}(\omega-\Omega)-1}}$
(3)

and is readily obtained once $D(\omega)$ is known.

Figure 3:  Spectral function $A({\bm{k}},\omega)$ at $T=0$ for four
representative $\gamma$. The broadening in the plots is intrinsic. (a-d):
color-coded plot at negative $\omega$, as measured by the ARPES intensity at T
= 0. (e-f): constant-${\bm{k}}$ cuts of $A({\bm{k}},\omega)$ at
$\xi_{\bm{k}}=0$ and at $\xi_{\bm{k}}=\pm 4\bar{g}$. For $\gamma<1/2$, the
spectral function has a sharp quasiparticle peak at
$\omega+\Delta\propto\xi^{2}_{{\bm{k}}}$. For $\gamma>1/2$, the peak moves to
$\omega+\Delta\propto|\xi_{{\bm{k}}}|^{1/(1-\gamma)}$ and broadens up, which
eventually disappears (see text).

At $\gamma=0$, which models a BCS superconductor, $C(\omega)=0$ and
$D(\omega)=A(\omega)/(\omega B(\omega))$ is a regular function of frequency.
Near the gap edge at $\omega>0$, $D(\omega)-1\sim\omega-\Delta$ and
$Z(\omega)\approx Z(\Delta)\equiv Z$. We assume and then verify that
$D(\omega)$ remains regular in some range of $\gamma>0$. Substituting
$D(\omega)-1\sim\omega-\Delta$ into (S7) for $\gamma>0$, we obtain
$C(\omega)-C(\Delta)\sim(\omega-\Delta)^{3/2-\gamma}$. We see that $C(\omega)$
is non-analytic near the gap edge, but for $\gamma<1/2$, the exponent
$3/2-\gamma$ is larger than one. In this situation, the non-analytic term in
$C(\omega)$ generates a non-analytic term in $D(\omega)$ of order
$(\omega-\Delta)^{3/2-\gamma}$, which is smaller than the regular
$\omega-\Delta$ term. Evaluating the prefactors, we obtain slightly above the
gap edge, at $\omega=\Delta+\delta$

$\displaystyle
D^{\prime}(\Delta+\delta)=1+\alpha\delta+A\cos[\pi(3/2-\gamma)]\delta^{3/2-\gamma},$
$\displaystyle
D^{\prime\prime}(\Delta+\delta)=-A\sin[\pi(3/2-\gamma)]\delta^{3/2-\gamma},$
(4)

where $\alpha\sim 1/{\bar{g}}$, $A=\sqrt{{\alpha\over
2}}{{\bar{g}}^{\gamma}\sin(\pi\gamma/2)\over\Delta B(\Delta)}J(\gamma,1)$ and
$J(\gamma,\nu)$ is expressed via Beta functions:

$J(\gamma,\nu)=B(1-\gamma,\gamma-1-{\nu\over 2})-B(1-\gamma,\gamma-1+{\nu\over
2}).$ (5)

For $\gamma>1/2$, $3/2-\gamma>1$, and the calculation of $D(\omega)$ has to be
done differently. We find after straightforward analysis that the leading
$\delta$-dependent term in $D(\Delta+\delta)$ is non-analytic and of order
$\delta^{\nu}$, where $\nu$ is the solution of $J(\gamma,\nu)=0$. The exponent
$\nu\approx 1+0.67(\gamma-1/2)$ for $\gamma\approx 1/2$ and $\nu\approx 1.3$
for $\gamma=1$. The subleading term in $D(\Delta+\delta)$ scales as
$\delta^{\nu+c}$, where $c>0$ is approximately linear in $\gamma-1/2$. In Fig.
2, we plot $\nu(\gamma)$ and $c(\gamma)$ along with the numerical results of
$D(\omega)$ for a representative $\gamma=0.8$. The exponent $\nu$ extracted
from this numerical $D(\omega)$ is $1.18$, which matches perfectly with the
analytical result. The behavior at $\gamma=1/2$ is special, and we discuss it
in Ref. [73].

Substituting $D(\Delta+\delta)$ into the formula for $Z(\omega)$, Eq. (S8), we
obtain

$\displaystyle Z^{\prime}(\Delta+\delta)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!Z(\Delta)\\!+\\!B\cos(\pi(\gamma+\nu/2-1))\delta^{1-\gamma-\nu/2},$
(6) $\displaystyle Z^{\prime\prime}(\Delta+\delta)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!B\sin(\pi(\gamma+\nu/2-1))\delta^{1-\gamma-\nu/2}.$ (7)

where
$B=\frac{\bar{g}^{\gamma}\sin\frac{\pi\gamma}{2}}{\Delta\sqrt{2\alpha}}B(1-\gamma,\frac{\nu}{2}+\gamma-1)$.
For $\gamma<1/2$, $Z(\omega)=Z(\Delta)+O(\delta^{1/2-\gamma})$ is
approximately a constant near the gap edge. For $\gamma>1/2$, the inverse
residue diverges at the gap edge, indicating a qualitative change in the
system behavior.

Spectral function and DoS. The spectral function and the DoS per unit volume
are given by

$\displaystyle
A({\bm{k}},\omega)=-\frac{1}{\pi}{\text{Im}}G_{R}({\bm{k}},\omega),$
$\displaystyle
N(\omega)=\frac{1}{V}\sum_{\bm{k}}A({\bm{k}},\omega)=N_{F}\omega{\text{Im}}\sqrt{\frac{1}{\Delta^{2}(\omega)-\omega^{2}}},$
(8)

where the retarded Green’s function $G_{R}({\bm{k}},\omega)=-(\omega
Z(\omega)+{\xi_{\bm{k}}})/({\xi_{\bm{k}}^{2}}+(\Delta^{2}(\omega)-\omega^{2})Z^{2}(\omega))$.
ARPES intensity is proportional to $A({\bm{k}},\omega)n_{F}(\omega)$, which at
$T=0$ selects negative $\omega$. At $\gamma=0$ (BCS limit), $N(\omega)\sim
1/(\omega-\Delta)^{1/2}$, and the spectral function has a $\delta$-functional
peak at $\omega=(\Delta^{2}+({\xi_{\bm{k}}}/Z)^{2})^{1/2}$. In Fig. 2 (c,d),
we show the DoS $N(\omega)$, obtained from the numerical solution of the full
gap equation (S3) for representative $\gamma=0.35$ and $0.8$. We see that in
both cases the DoS describes a gapped continuum, but there is a qualitative
difference in the behavior near the gap edge: for $\gamma=0.35$, $N(\omega)$
has the same $1/\delta^{1/2}$ singularity as for $\gamma=0$, and for
$\gamma=0.8$ the DOS behaves as $1/\delta^{0.59}$, which perfectly matches the
analytical form $\delta^{-\nu/2}$, given that $\nu=1.18$ for $\gamma=0.8$.

The spectral function $A({\bm{k}},\omega)$ is shown in Fig. (3). For
comparison with ARPES, we set $\omega$ to be negative:
$\omega=-(\Delta+\delta)$. For any $\gamma$, there is no frequency range,
where $A({\bm{k}},\omega)$ is a $\delta$-function, simply because the bosonic
mass vanishes at a QCP. Still, for $\gamma<1/2$,
$D(-(\Delta+\delta))-1\propto\delta$ and $Z(-(\Delta+\delta))\approx
Z(-\Delta)=Z(\Delta)$. In this situation, the spectral weight on the Fermi
surface, integrated over an infinitesimally small range around
$\omega=-\Delta$ immediately above the real axis, is finite, like in BCS case.
Away from the Fermi surface, the spectral function vanishes as
$|\omega+\Delta|^{1/2-\gamma}$ at the gap edge and displays a quasiparticle
peak at $\omega\approx-(\Delta^{2}+({\xi_{\bm{k}}}/Z(\Delta))^{2})^{1/2}$. The
peak is well defined at small $\delta$ as its width $O(\delta^{1/2-\gamma})$
is parametrically smaller than its frequency. This is the same behavior as in
Fig. 1 (a). For $\gamma>1/2$, the situation is qualitatively different. Now
$Z(-\Delta-\delta)$ diverges at $\delta\to 0$ and
$D(-\Delta-\delta)-1\sim|\delta|^{\nu}\ll|\delta|$. In this case, the integral
of $A({\bm{k}}_{F},\omega)$ over an infinitesimally small range around
$\omega=-\Delta$ vanishes, which can be interpreted as a vanishing of a
quasiparticle peak. At finite $\xi_{\bm{k}}$, the spectral function diverges
at the gap edge as $1/|\omega+\Delta|^{\gamma/2+\gamma-1}$. For $\gamma$
slightly above $1/2$, $A({\bm{k}},\omega)$ is non-monotonic and possess a
broad maximum at
$|\omega+\Delta|\sim\left({{\xi_{\bm{k}}}/{\bar{g}}^{\gamma}}\right)^{1\over
1-\gamma}$. This is the same behavior as in Fig. 1 (b). For larger $\gamma$,
the maximum disappears, and $A({\bm{k}},\omega)$ monotonically decreases at
$|\omega|>\Delta$. This is the same behavior as in Fig. 1 (c). For small
$\xi_{\bm{k}}$, the maximum disappears at $\gamma\sim 0.9$. For larger
$\xi_{\bm{k}}$, it disappears at smaller $\gamma$, first for positive
$\xi_{{\bm{k}}}$ (see Fig. 4).

Comparison with ARPES The behavior shown in Fig. 4 is our result in some range
of $\gamma>1/2$. For positive $\xi_{\bm{k}}$ (i.e., outside the Fermi
surface), the spectral function has a single non-dispersing maximum at the gap
edge, except for the smallest $\xi_{{\bm{k}}}$, while for negative
$\xi_{\bm{k}}$, $A({\bm{k}},\omega)$ has a kink at the gap edge
$\omega=-\Delta$ and a dispersing maximum at
$\omega=-\Delta-O\left(|{\xi_{k}}|^{1/(1-\gamma)}\right)$. This behavior is
consistent with the ARPES data for Bi2201, Ref. [75]. The data shows that the
spectral function near the antinode, where our analysis is valid, displays an
almost non-dispersing maximum at positive $\xi_{\bm{k}}$, while for negative
$\xi_{\bm{k}}$ it displays a non-dispersing kink at the same energy and a
dispersing maximum at larger $|\omega|$. We associate the non-dispersing
feature at both positive and negative $\xi_{\bm{k}}$ with the gap edge
$\Delta$, and associate the dispersing maximum, observed in [75] at
$\xi_{\bm{k}}<0$, with the dispersing maximum in Fig. 4.

Figure 4: (a) Spectral function $A({\bm{k}},\omega)$ at positive and negative
$\xi_{\bm{k}}=\pm 4{\bar{g}}$ at $\gamma=0.6$. To account for impurity
scattering, we convoluted the spectral function with a Lorentzian of width
$\sim 0.03{\bar{g}}$. (b) Spectral function at a set of discrete momenta. It
displays a non-dispersing gap edge singularity (green dots) and a dispersing
maximum (blue circles). This theoretical $A({\bm{k}},\omega)$ is consistent
with the ARPES data for Bi2201, Ref. [75] (see text).

Discussion and summary. In this work, we analyzed the applicability of
quasiparticle description of a superconducting state which emerges out of a
non-Fermi liquid at a metallic QCP. We considered the model with an effective
dynamical 4-fermion interaction $V(\Omega)\propto 1/\Omega^{\gamma}$, mediated
by a gapless boson at a QCP and analyzed the spectral function and the DoS for
$\gamma\in(0,1)$. Interaction $V(\Omega)$ gives rise to a non-Fermi liquid in
the normal state with self-energy $\Sigma(\omega)\propto\omega^{1-\gamma}$ and
to pairing below some finite $T_{c}$. A superconducting order gaps out low-
energy excitations and, at a first glance, should restore fermionic coherence.
We found, however, that this holds only for $\gamma<1/2$. For larger $\gamma$
the spectral function and the DoS exhibit qualitatively different behavior
than that in a superconductor with coherent quasiparticles. (different power-
laws). We argued that the quasiparticle peak broadens up and completely
disappears for $\gamma$ close to one.

Away from a QCP, a pairing boson is massive and at the lowest energies a
Fermi-liquid description holds already in the normal state and continue to
hold in a superconductor. In particular, in the immediate vicinity of the gap
edge, the system displays a BCS-like behavior for all $\gamma$. Still, the
system behavior over a broad frequency range is governed by the physics at a
QCP, as numerous experiments on the cuprates and other correlated systems
indicate. We argued that our results are quite consistent with the ARPES data
for Bi2201 [75, 12].

###### Acknowledgements.

We acknowledge with thanks useful conversations with a number of our
colleagues. This work was supported by the U.S. Department of Energy, Office
of Science, Basic Energy Sciences, under Award No. DE-SC0014402.

## References

  * Martin  _et al._ [1990] S. Martin, A. T. Fiory, R. Fleming, L. Schneemeyer, and J. V. Waszczak, Physical Review B 41, 846 (1990).
  * Löhneysen  _et al._ [1994] H. v. Löhneysen, T. Pietrus, G. Portisch, H. Schlager, A. Schröder, M. Sieck, and T. Trappmann, Physical review letters 72, 3262 (1994).
  * Norman  _et al._ [1998] M. R. Norman, H. Ding, M. Randeria, J. C. Campuzano, T. Yokoya, T. Takeuchi, T. Takahashi, T. Mochiku, K. Kadowaki, P. Guptasarma, _et al._ , Nature 392, 157 (1998).
  * Shen and Sawatzky [1999] Z.-X. Shen and G. Sawatzky, physica status solidi (b) 215, 523 (1999).
  * Grigera  _et al._ [2001] S. Grigera, R. Perry, A. Schofield, M. Chiao, S. Julian, G. Lonzarich, S. Ikeda, Y. Maeno, A. Millis, and A. Mackenzie, Science 294, 329 (2001).
  * Damascelli  _et al._ [2003] A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod. Phys. 75, 473 (2003).
  * Inosov  _et al._ [2007] D. S. Inosov, J. Fink, A. A. Kordyuk, S. V. Borisenko, V. B. Zabolotnyy, R. Schuster, M. Knupfer, B. Büchner, R. Follath, H. A. Dürr, W. Eberhardt, V. Hinkov, B. Keimer, and H. Berger, Phys. Rev. Lett. 99, 237002 (2007).
  * Daou  _et al._ [2009] R. Daou, N. Doiron-Leyraud, D. LeBoeuf, S. Li, F. Laliberté, O. Cyr-Choiniere, Y. Jo, L. Balicas, J.-Q. Yan, J.-S. Zhou, _et al._ , Nature Physics 5, 31 (2009).
  * Cooper  _et al._ [2009] R. Cooper, Y. Wang, B. Vignolle, O. Lipscombe, S. Hayden, Y. Tanabe, T. Adachi, Y. Koike, M. Nohara, H. Takagi, _et al._ , Science 323, 603 (2009).
  * Sarkar  _et al._ [2017] T. Sarkar, P. Mandal, J. Higgins, Y. Zhao, H. Yu, K. Jin, and R. L. Greene, Physical Review B 96, 155449 (2017).
  * Jin  _et al._ [2011] K. Jin, N. Butch, K. Kirshenbaum, J. Paglione, and R. Greene, Nature 476, 73 (2011).
  * Hashimoto  _et al._ [2014] M. Hashimoto, I. M. Vishik, R.-H. He, T. P. Devereaux, and Z.-X. Shen, Nature Physics 10, 483 (2014).
  * Kaminski  _et al._ [2015] A. Kaminski, T. Kondo, T. Takeuchi, and G. Gu, Philosophical Magazine 95, 453 (2015), https://doi.org/10.1080/14786435.2014.906758 .
  * Keimer  _et al._ [2015] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, Nature 518, 179 (2015).
  * Taillefer [2010] L. Taillefer, Annual Review of Condensed Matter Physics 1, 51 (2010), https://doi.org/10.1146/annurev-conmatphys-070909-104117 .
  * Legros  _et al._ [2019] A. Legros, S. Benhabib, W. Tabis, F. Laliberté, M. Dion, M. Lizaire, B. Vignolle, D. Vignolles, H. Raffy, Z. Li, _et al._ , Nature Physics 15, 142 (2019).
  * Varma [2020] C. M. Varma, Rev. Mod. Phys. 92, 031001 (2020).
  * Hartnoll and Mackenzie [2022] S. A. Hartnoll and A. P. Mackenzie, Rev. Mod. Phys. 94, 041002 (2022).
  * Harada  _et al._ [2022] K. Harada, Y. Teramoto, T. Usui, K. Itaka, T. Fujii, T. Noji, H. Taniguchi, M. Matsukawa, H. Ishikawa, K. Kindo, D. S. Dessau, and T. Watanabe, Phys. Rev. B 105, 085131 (2022).
  * Löhneysen  _et al._ [2007] H. v. Löhneysen, A. Rosch, M. Vojta, and P. Wölfle, Rev. Mod. Phys. 79, 1015 (2007).
  * Valla  _et al._ [1999] T. Valla, A. Fedorov, P. Johnson, B. Wells, S. Hulbert, Q. Li, G. Gu, and N. Koshizuka, Science 285, 2110 (1999).
  * Abanov  _et al._ [2003] A. Abanov, A. V. Chubukov, and J. Schmalian, Advances in Physics 52, 119 (2003).
  * Abanov  _et al._ [2001a] A. Abanov, A. V. Chubukov, and J. Schmalian, Journal of Electron spectroscopy and related phenomena 117, 129 (2001a).
  * Fink  _et al._ [2006] J. Fink, A. Koitzsch, J. Geck, V. Zabolotnyy, M. Knupfer, B. Büchner, A. Chubukov, and H. Berger, Phys. Rev. B 74, 165102 (2006).
  * Shibauchi  _et al._ [2014] T. Shibauchi, A. Carrington, and Y. Matsuda, Annu. Rev. Condens. Matter Phys. 5, 113 (2014).
  * Restrepo  _et al._ [2022] F. Restrepo, U. Chatterjee, G. Gu, H. Xu, D. K. Morr, and J. C. Campuzano, Communications Physics 5, 45 (2022).
  * Scalapino [2012] D. J. Scalapino, Rev. Mod. Phys. 84, 1383 (2012).
  * Efetov  _et al._ [2013] K. B. Efetov, H. Meier, and C. Pepin, Nature Physics 9, 442 (2013).
  * Lederer  _et al._ [2017] S. Lederer, Y. Schattner, E. Berg, and S. A. Kivelson, Proceedings of the National Academy of Sciences 114, 4905 (2017).
  * Hartnoll  _et al._ [2011] S. A. Hartnoll, D. M. Hofman, M. A. Metlitski, and S. Sachdev, Phys. Rev. B 84, 125115 (2011).
  * Wang and Torroba [2017] H. Wang and G. Torroba, Phys. Rev. B 96, 144508 (2017).
  * Bonesteel  _et al._ [1996] N. E. Bonesteel, I. A. McDonald, and C. Nayak, Phys. Rev. Lett. 77, 3009 (1996).
  * Abanov  _et al._ [2001b] A. Abanov, A. V. Chubukov, and A. M. Finkel’stein, EPL (Europhysics Letters) 54, 488 (2001b).
  * Chubukov  _et al._ [2020] A. V. Chubukov, A. Abanov, Y. Wang, and Y.-M. Wu, Annals of Physics 417, 168142 (2020).
  * Metlitski  _et al._ [2015] M. A. Metlitski, D. F. Mross, S. Sachdev, and T. Senthil, Phys. Rev. B 91, 115111 (2015).
  * Perali  _et al._ [1996] A. Perali, C. Castellani, C. Di Castro, and M. Grilli, Phys. Rev. B 54, 16216 (1996).
  * Moon and Sachdev [2009] E. G. Moon and S. Sachdev, Phys. Rev. B 80, 035117 (2009).
  * Wang and Chubukov [2013] Y. Wang and A. V. Chubukov, Phys. Rev. Lett. 110, 127001 (2013).
  * Wang  _et al._ [2016] Y. Wang, A. Abanov, B. L. Altshuler, E. A. Yuzbashyan, and A. V. Chubukov, Phys. Rev. Lett. 117, 157001 (2016).
  * Wang  _et al._ [2017] H. Wang, S. Raghu, and G. Torroba, Phys. Rev. B 95, 165137 (2017).
  * Son [1999] D. T. Son, Phys. Rev. D 59, 094019 (1999).
  * Chubukov and Schmalian [2005] A. V. Chubukov and J. Schmalian, Phys. Rev. B 72, 174520 (2005).
  * Abanov and Chubukov [2020] A. Abanov and A. V. Chubukov, Phys. Rev. B 102, 024524 (2020).
  * Yuzbashyan  _et al._ [2022] E. A. Yuzbashyan, M. K.-H. Kiessling, and B. L. Altshuler, Phys. Rev. B 106, 064502 (2022).
  * Zhang  _et al._ [2022] S.-S. Zhang, Y.-M. Wu, A. Abanov, and A. V. Chubukov, Physical Review B 106, 144513 (2022).
  * Kaminski  _et al._ [2000] A. Kaminski, J. Mesot, H. Fretwell, J. C. Campuzano, M. R. Norman, M. Randeria, H. Ding, T. Sato, T. Takahashi, T. Mochiku, K. Kadowaki, and H. Hoechst, Phys. Rev. Lett. 84, 1788 (2000).
  * Oganesyan  _et al._ [2001] V. Oganesyan, S. A. Kivelson, and E. Fradkin, Phys. Rev. B 64, 195109 (2001).
  * Metzner  _et al._ [2003] W. Metzner, D. Rohe, and S. Andergassen, Phys. Rev. Lett. 91, 066402 (2003).
  * Dell’Anna and Metzner [2006] L. Dell’Anna and W. Metzner, Phys. Rev. B 73, 045127 (2006).
  * Rech  _et al._ [2006] J. Rech, C. Pépin, and A. V. Chubukov, Phys. Rev. B 74, 195126 (2006).
  * Metlitski and Sachdev [2010a] M. A. Metlitski and S. Sachdev, Physical Review B 82, 075127 (2010a).
  * Metlitski and Sachdev [2010b] M. A. Metlitski and S. Sachdev, Physical Review B 82, 075128 (2010b).
  * Mross  _et al._ [2010] D. F. Mross, J. McGreevy, H. Liu, and T. Senthil, Phys. Rev. B 82, 045121 (2010).
  * Raghu  _et al._ [2015] S. Raghu, G. Torroba, and H. Wang, Phys. Rev. B 92, 205104 (2015).
  * Fitzpatrick  _et al._ [2015] A. L. Fitzpatrick, S. Kachru, J. Kaplan, S. Raghu, G. Torroba, and H. Wang, Phys. Rev. B 92, 045118 (2015).
  * Klein  _et al._ [2020] A. Klein, A. V. Chubukov, Y. Schattner, and E. Berg, Phys. Rev. X 10, 031053 (2020).
  * Lee [2009] S.-S. Lee, Phys. Rev. B 80, 165102 (2009).
  * Dalidovich and Lee [2013] D. Dalidovich and S.-S. Lee, Phys. Rev. B 88, 245106 (2013).
  * Schlief  _et al._ [2017] A. Schlief, P. Lunts, and S.-S. Lee, Phys. Rev. X 7, 021010 (2017).
  * Punk [2016] M. Punk, Phys. Rev. B 94, 195113 (2016).
  * Maslov and Chubukov [2010] D. L. Maslov and A. V. Chubukov, Phys. Rev. B 81, 045110 (2010).
  * Dessau  _et al._ [1991] D. S. Dessau, B. O. Wells, Z.-X. Shen, W. E. Spicer, A. J. Arko, R. S. List, D. B. Mitzi, and A. Kapitulnik, Phys. Rev. Lett. 66, 2160 (1991).
  * Kordyuk  _et al._ [2002] A. A. Kordyuk, S. V. Borisenko, T. K. Kim, K. A. Nenkov, M. Knupfer, J. Fink, M. S. Golden, H. Berger, and R. Follath, Phys. Rev. Lett. 89, 077003 (2002).
  * [64] S.-S. Zhang, E. Berg, and A. V. Chubukov, arXiv:2301.01873, and references therein .
  * Esterlis and Schmalian [2019] I. Esterlis and J. Schmalian, Phys. Rev. B 100, 115132 (2019).
  * Hauck  _et al._ [2020] D. Hauck, M. J. Klug, I. Esterlis, and J. Schmalian, Annals of Physics 417, 168120 (2020).
  * Wang [2020] Y. Wang, Phys. Rev. Lett. 124, 017002 (2020).
  * Chowdhury and Berg [2020] D. Chowdhury and E. Berg, Phys. Rev. Research 2, 013301 (2020).
  * Classen and Chubukov [2021] L. Classen and A. Chubukov, Phys. Rev. B 104, 125120 (2021).
  * [70] J. Schmalian, private communication.
  * Wu  _et al._ [2021a] Y.-M. Wu, S.-S. Zhang, A. Abanov, and A. V. Chubukov, Phys. Rev. B 103, 024522 (2021a).
  * Wu  _et al._ [2021b] Y.-M. Wu, S.-S. Zhang, A. Abanov, and A. V. Chubukov, Phys. Rev. B 103, 184508 (2021b).
  * [73] See the Supplementary Information for more details.
  * Marsiglio  _et al._ [1988] F. Marsiglio, M. Schossmann, and J. P. Carbotte, Phys. Rev. B 37, 4965 (1988).
  * He  _et al._ [2011] R.-H. He, M. Hashimoto, H. Karapetyan, J. Koralek, J. Hinton, J. Testaud, V. Nathan, Y. Yoshida, H. Yao, K. Tanaka, _et al._ , Science 331, 1579 (2011).

Supplementary information for “Density of states and spectral function of a
superconductor out of a quantum-critical metal”
by Shang-Shun Zhang and Andrey V Chubukov

##  I Gap equation along the real-frequency axis and its solution

We will use the approach pioneered by Marsiglio, Shossmann, and Carbotte [1].
In this approach, one first solves non-linear gap equation along the Matsubara
axis, which can be done rather straightforwardly as the gap function
$\Delta(\omega_{m})$ can be chosen to be real for all frequencies and is a
regular function of $\omega_{m}$ even when the pairing boson is massless. One
then uses this $\Delta(\omega_{m})$ as an input for the equation for complex
$\Delta(\omega)$ along the real frequency axis.

The non-linear integral equation for
$D(\omega_{m})=\omega_{m}\Delta(\omega_{m})$ on the Matsubara axis, and the
equation for the inverse quasiparticle residue
$Z(\omega_{m})=1+\Sigma(\omega_{m})/\omega_{m}$ ($\Sigma(\omega_{m})$ is the
fermionic self-energy), have the form

$\omega_{m}D(\omega_{m})=\pi
T\sum_{\omega_{m}^{\prime}}\frac{\left(D(\omega_{m}^{\prime})-D(\omega_{m})\right)\text{sgn}(\omega_{m}^{\prime})}{\sqrt{1+D^{2}(\omega_{m}^{\prime})}}V(\omega-\omega_{m}^{\prime}),$
(S1) $Z(\omega_{m})=1+\frac{1}{\omega_{m}}\pi
T\sum_{\omega_{m}^{\prime}}\frac{\text{sgn}(\omega_{m}^{\prime})}{\sqrt{1+D^{2}(\omega_{m}^{\prime})}}V(\omega-\omega_{m}^{\prime}),$
(S2)

where $V(\Omega_{m})=\left({\bar{g}}/{\rvert\Omega_{m}\rvert}\right)^{\gamma}$
is the same as in the main text ($\bar{g}$ is an effective fermion-boson
coupling, and $\gamma$ depends on the underlying microscopic model).This set
of equations has a non-zero solution $D(\omega_{m})$ below a finite pairing
temperature $T_{p}\sim{\bar{g}}$. Fig. S1 shows the numerical solution for
$\Delta(\omega_{m})$ at $T=10^{-6}{\bar{g}}\ll T_{p}$, for different $\gamma$.
We see from the figure that $\Delta(\omega_{m})$ approaches a constant value
at small frequencies and decays as $\omega_{m}^{-\gamma}$ at high frequencies.
This behavior holds for all $\gamma$ and can be easily verified analytically.

Figure S1: Numerical resulls for the gap function $\Delta(\omega_{m})$ along
the Matsubara axis. The calculation is performed at temperature
$T=10^{-6}{\bar{g}}$ using the hybrid-frequency method [2].

The gap equation along the real frequency axis is

$\omega B(\omega)D(\omega)=A(\omega)+C(\omega),$ (S3)

(Eqn. (1) in the main text). This equation is obtained by using the spectral
representation of an analytic function on the upper frequency half-plane

$f(i\omega_{m})=\frac{1}{\pi}\int dx\frac{\text{Im}f(x)}{x-i\omega_{m}}$ (S4)

and, where possible, keeping $D(\omega_{m})$ as an input function. This
approach was pioneered for electron-phonon interaction in Refs. [3, 4, 5, 6]
for the electron-phonon problem.

For our case, the functions $A(\omega)$ and $B(\omega)$ are directly expressed
via $D(\omega_{m})$ along the Matsubara axis as

$\displaystyle A(\omega)$
$\displaystyle=\frac{1}{2}\int_{0}^{\infty}d\omega_{m}\frac{D(\omega_{m})}{\sqrt{1+D^{2}(\omega_{m})}}$
$\displaystyle\times\left(\frac{\bar{g}^{\gamma}}{(\omega_{m}+i\omega)^{\gamma}}+\frac{\bar{g}^{\gamma}}{(\omega_{m}-i\omega)^{\gamma}}\right),$
(S5) $\displaystyle B(\omega)$
$\displaystyle=1+\frac{i}{2\omega}\int_{0}^{\infty}d\omega_{m}\frac{1}{\sqrt{1+D^{2}(\omega_{m})}}$
$\displaystyle\times\left(\frac{\bar{g}^{\gamma}}{(\omega_{m}+i\omega)^{\gamma}}-\frac{\bar{g}^{\gamma}}{(\omega_{m}-i\omega)^{\gamma}}\right).$
(S6)

and $C(\omega)$ is given by

$C(\omega)=\bar{g}^{\gamma}\sin\frac{\pi\gamma}{2}\int_{0}^{\omega}\frac{d\Omega}{\Omega^{\gamma}}\frac{D(\omega-\Omega)-D(\omega)}{\sqrt{D^{2}(\omega-\Omega)-1}},$
(S7)

(Eqn (3) in the main text). This function depends on the running
$D(\omega-\Omega)$, which makes Eq. (S3) an integral equation. The inverse
residue $Z(\omega)$ is expressed via $D(\omega^{\prime})$ as

$\displaystyle
Z(\omega)=B(\omega)+\frac{\bar{g}^{\gamma}\sin\frac{\pi\gamma}{2}}{\omega}\int_{0}^{\omega}\frac{d\Omega}{\Omega^{\gamma}}\frac{1}{\sqrt{D^{2}(\omega-\Omega)-1}}$
(S8)

(Eqn (4) in the main text) and is readily obtained once $D(\omega)$ is known.

The gap equation along the real-frequency axis has an iterative structure in
the sense that $D(\omega)$ depends on $D(\omega^{\prime})$ at
$\omega^{\prime}<\omega$. This allows us to solve this equation iteratively,
using the low-frequency form $D(\omega)\simeq\Delta(0)/\omega$ as an input,
with $\Delta(0)\equiv\Delta(\omega_{m}=\pi T)$. In Fig. S2 we show the results
for $D(\omega)$ and $Z(\omega)$ for three representative values of $\gamma$.
In all cases, $D(\omega)$ and $Z(\omega)$ are real below the gap edge
$\omega=\Delta$ and are complex above the gap edge, where $\Delta$ is defined
as $\Delta(\omega)=1$ at $\omega=\Delta$. We see that for $\gamma>1/2$,
$Z(\omega)$ diverges at the gap edge. We use this fact in the main text in the
analysis of the spectral function.

Figure S2: Numerical results for $D(\omega)$ and $Z(\omega)$ for
$\gamma=0.25$, $\gamma=0.5$ and $\gamma=0.8$.

##  II The case of $\gamma=1/2$

In the main text we argued that for $\gamma<1/2$, the function
$D(\Delta-\delta)-1\propto\delta$, where $\delta=\Delta-\omega$, and the
correction scales as $\delta^{3/2-\gamma}$. More specifically, we found
iteratively that

$D(\Delta-\delta)=1+\delta\sum_{n=0}^{\infty}\alpha_{n}\delta^{n\epsilon}$
(S9)

where $\epsilon=1/2-\gamma$ and $\alpha_{0}=O(1/{\bar{g}})$. The expression
for $\alpha_{1}$ is presented in the main text, after Eq. (4). It is
proportional to $J(\gamma,1)=B(1-\gamma,\gamma-3/2)-B(1-\gamma,\gamma-1/2)$,
where $B(a,b)$ is a Beta function ($B(a,b)=\Gamma(a)\Gamma(b)/\Gamma(a+b)$).
For small $\epsilon$ (i.e., for $\gamma\leq 1/2$), $\alpha_{1}\sim
J(\gamma,1)\sim 1/\epsilon$. For the next term in (S9) we find $\alpha_{2}\sim
1/\epsilon^{2}$, and so on.

We see that the perturbative expansion in $\delta^{\epsilon}$ in (S9) holds
for $(\delta/{\bar{g}})^{\epsilon}/\epsilon\leq 1$. Outside this range, all
terms in Eq. (S9) are relevant. As $\gamma$ approaches $1/2$ from below and
$\epsilon$ decreases, the perturbative regime shrinks to exponentially small
$\delta<{\bar{g}}\exp(-|\log{\epsilon}|/\epsilon)$.

To understand the form of $D(\omega)$ outside the perturbative regime, we
express $(\delta/{\bar{g}})^{\epsilon}$ as
$e^{\epsilon\log{(\delta/{\bar{g}})}}$ and expand (S10) in powers of
$\log{(\delta/{\bar{g}})}$. We obtain

$\displaystyle
D(\Delta-\delta)=1+\delta\sum_{n=0}^{\infty}{\tilde{\alpha}}_{n}(\log{\frac{\delta}{{\bar{g}}}})^{n}$
(S10)

where ${\tilde{\alpha}}_{0}=\alpha_{0}+\alpha_{1}+\alpha_{2}+..$,
${\tilde{\alpha}}_{1}=\epsilon\alpha_{1}+2\epsilon\alpha_{2}+..$,
${\tilde{\alpha}}_{2}=2\epsilon^{2}\alpha_{2}+...$. We see that each
${\tilde{\alpha}}_{n}$ is a series, in which the first term is independent on
$\epsilon$, and the others diverge as powers of $1/\epsilon$, because
$\alpha_{1}\sim 1/\epsilon$, $\alpha_{2}\sim 1/\epsilon^{2}$, and so on.

We now argue that singular parts of ${\tilde{\alpha}}_{n}$ can be neglected.
The argument is two-fold. First, in the calculations, the $1/\epsilon$
divergencies originate from the divergence of $J(1/2-\epsilon,1)\approx
1/(2\epsilon)$. This divergence is regularized by a finite boson mass, such
that strictly at $\epsilon=0$, one has $J(1/2,1)=0$ instead of infinity.
Second, if we assume that ${\tilde{\alpha}}_{0}$ in (S10) remains finite at
$\epsilon=0$ and substitute the trial
$D(\Delta-\delta)=1+\delta{\tilde{\alpha}}_{0}$ in the gap equation at
$\gamma=1/2$ and compute iteratively the next term in $D(\Delta-\delta)$, we
find it in the form ${\tilde{\alpha}}_{1}\delta\log{(\delta/{\bar{g}})}$ with
a finite $\tilde{\alpha}_{1}=\sqrt{{\bar{g}}\tilde{\alpha}_{0}}/(4\Delta
B(\Delta))$. Extending the iterative analysis, we find that all
${\tilde{\alpha}}_{n}$ are finite at $\gamma=1/2$ (i.e., $\epsilon=0$), as we
anticipated.

Figure S3: The real and imaginary parts of the gap function near the gap edge
$\omega=\Delta$ for $\gamma=1/2$, obtained by solving the non-linear gap
equation numerically. The leading term in $D^{\prime}(\Delta-\delta)-1$, shown
in the inset of (a), is linear in $\delta\Delta-\omega$, and the subleading
scales as $\delta/|\log{|\delta|}|$, as is confirmed by the linear relation in
panel (a), The imaginary part $D^{\prime\prime}$ appears at negative $\delta$
above the gap edge. The numerical result in panel (b) clearly shows the
scaling relation $D^{\prime\prime}\sim\delta/\log^{2}{(|\delta|/{\bar{g}})}$,
expected from the Kramers-Kronig relation with $D^{\prime}$.

We didn’t manage to sum up analytically the logarithmic series in (S10). The
numerical solution for $D(\Delta-\delta)$ for $\gamma=1/2$ shows that
$D(\Delta-\delta)-1$ remains linear in $\delta$ (see Fig. S2 (c)), and the
corrections scale as $1/|\log{(|\delta|/{\bar{g}})}|$ (see Fig. S3 (a)). By
Kramers-Kronig relation, this implies that at negative $\delta$, when
$\omega>\Delta$ is above the gap edge, the imaginary part of $D(\omega)$
scales as
$D^{{}^{\prime\prime}}(\Delta+|\delta|)\propto\delta/\log^{2}{(|\delta|/{\bar{g}})}$
(the same form is obtained by just noticing that
$\log({-|\delta|})=\log({-(\omega-\Delta+i0)})=\log{|\delta|}-i\pi$). This
form of $D^{{}^{\prime\prime}}(\Delta+|\delta|)$ is consistent with our
numerical solution above the gap edge, Fig. S3 (b). The solution clearly shows
that the ratio $D^{{}^{\prime\prime}}(\Delta+|\delta|)/\delta$ decreases at
the smallest $\delta$.

Figure S4: (a) The spectral function $A({\bm{k}},\omega)$ for $\gamma=1/2$.
(b) Constant $\xi_{\bm{k}}$ cuts along the blue lines in panel (a).

We next use the result for $D(\omega)$ to obtain the inverse quasi-particle
residue near the gap edge. Substituting $D(\Delta-\delta)\approx
1+{\tilde{\alpha}_{0}}\delta$ into (S8), we obtain at $\gamma=1/2$

$\displaystyle Z(\Delta-\delta)={1\over
2\Delta}\sqrt{{{\bar{g}}\over\tilde{\alpha}_{0}}}|\log{\delta}|.$ (S11)

Analytically continuing this function to negative $\delta$, i.e., to $\omega$
above the threshold, we obtain

$\displaystyle Z(\Delta+|\delta|)={1\over
2\Delta}\sqrt{{{\bar{g}}\over\tilde{\alpha}_{0}}}\left(|\log{|\delta|}|+i\pi\right),$
(S12)

Note that the imaginary part of $Z(\omega)$ jumps to a finite value at
$\omega$ infinitesimally above the threshold. This behavior is consistent with
the numerical solution for $Z(\omega)$, see Fig. S2 (d).

Finally, we use the results for $D(\omega)$ and $Z(\omega)$ and compute the
spectral function near the gap edge. On the Fermi surface, the spectral
function at negative $\omega$ and $|\omega|>\Delta$ takes the form

$\displaystyle
A({\bm{k}}_{F},\omega)\propto{1\over\rvert\omega+\Delta\rvert\log({\bar{g}}/\rvert\omega+\Delta\rvert)}.$
(S13)

Slightly away from the Fermi surface, the spectral function has a peak at
$|\omega|=\Delta+\delta_{\bm{k}}$ where
$\delta_{\bm{k}}\sim(\xi^{2}_{\bm{k}}/{\bar{g}})/\log^{2}(|{\bar{g}}/\xi_{\bm{k}}|)$.
The peak width scales as $\delta_{\bm{k}}/\log(|{\bar{g}}/\xi_{\bm{k}}|)$ and
is logarithmically smaller than the energy variation $|\omega|-\Delta$. Also,
for any non-zero $\xi_{\bm{k}}$, the spectral function jumps at the gap edge
to a finite value of order $1/\xi_{\bm{k}}^{2}$. In Fig. S4 we show the
numerical result for the spectral function. It is consistent with the behavior
we just described.

##  III Universal form of the spectral function at $1/2<\gamma<1$

For frequencies $\omega$ near the gap edge and for momenta near the Fermi
surface, when $\xi_{\bf k}$ is much smaller than $|\omega Z(\omega)|$, a
straightforward calculation shows that for $\gamma>1/2$, the spectral function
can be expressed as a scaling function of $\xi_{\bf
k}/|\omega+\Delta|^{1-\gamma}{\bar{g}}^{\gamma}$ (we set $\omega<0$). Namely,

$A({\bm{k}},\omega)\propto\frac{1}{\rvert\omega+\Delta\rvert^{\frac{\nu}{2}+1-\gamma}}\Phi\left(\frac{\xi_{\bm{k}}}{\rvert\omega+\Delta\rvert^{1-\gamma}{\bar{g}}^{\gamma}}\right),$
(S14)

where

$\Phi(x)\equiv\frac{x^{2}+Q_{\gamma}^{2}\sin[\pi(\nu-c)]/\sin(\pi
c)}{\left(x^{2}+Q_{\gamma}^{2}\cos(2\pi\gamma)\right)^{2}+Q_{\gamma}^{4}\sin^{2}(2\pi\gamma)}$
(S15)

with $Q_{\gamma}=\sin(\pi\gamma/2)B(1-\gamma,\nu/2+\gamma-1)$. In Fig. S5, we
plot the dimensionless function $\Phi(x)$ for different values of $\gamma$. At
$x\gg 1$, $\Phi(x)\sim 1/x^{2}$; at $x\ll 1$, $\Phi(x)\sim$ const. For
$\gamma<\gamma_{c}\simeq 0.9$, function $\Phi(x)$ contains a local maximum at

$\displaystyle x_{*}^{2}\sim\sqrt{(u-v)^{2}+w^{2}}-u,$ (S16)

where $u=Q_{\gamma}^{2}\sin[\pi(\nu-c)]/\sin(\pi c)$,
$v=Q_{\gamma}^{2}\cos(2\pi\gamma)$, and $w=Q_{\gamma}^{2}\sin(2\pi\gamma)$.
This maximum can be interpreted as an over-damped, but still existing quasi-
particle peak. At $\gamma>\gamma_{c}$, the function $\Phi(x)$ monotonically
decreases with $x$. In this case, the quasiparticle description breaks down
completely.

Figure S5: The function $\Phi(x)$, Eqn (S15), for different values of
$\gamma$.

We emphasize that this behavior holds only for small enough $\xi_{\bf k}$. For
larger $\xi_{\bf k}$, the spectral function does depend on the sign of
$\xi_{\bf k}$ and as $\gamma$ increases, the quasiparticle behavior gets
completely destroyed first for positive $\xi_{\bf k}$ and then, at larger
$\gamma$, for negative $\xi_{\bf k}$.

## References

  * Marsiglio  _et al._ [1988] F. Marsiglio, M. Schossmann, and J. P. Carbotte, Phys. Rev. B 37, 4965 (1988).
  * Wu  _et al._ [2020] Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Phys. Rev. B 102, 024525 (2020).
  * Marsiglio and Carbotte [1991] F. Marsiglio and J. P. Carbotte, Phys. Rev. B 43, 5355 (1991), for more recent results see F. Marsiglio and J.P. Carbotte, “Electron-Phonon Superconductivity”, in “The Physics of Conventional and Unconventional Superconductors”, Bennemann and Ketterson eds., Springer-Verlag, (2006) and references therein; F. Marsiglio, Annals of Physics 417, 168102-1-23 (2020).
  * Karakozov  _et al._ [1991] A. Karakozov, E. Maksimov, and A. Mikhailovsky, Solid State Communications 79, 329 (1991).
  * Combescot [1995] R. Combescot, Phys. Rev. B 51, 11625 (1995).
  * Wu  _et al._ [2019] Y.-M. Wu, A. Abanov, Y. Wang, and A. V. Chubukov, Phys. Rev. B 99, 144512 (2019).