# Plasmonic gain in current biased tilted Dirac nodes

Sang Hyun Park  Department of Electrical & Computer Engineering, University of
Minnesota, Minneapolis, Minnesota, 55455, USA  Michael Sammon  Department of
Electrical & Computer Engineering, University of Minnesota, Minneapolis,
Minnesota, 55455, USA  Eugene Mele  Department of Physics and Astronomy,
University of Pennsylvania, Philadelphia, Pennsylvania, 19104, USA  Tony Low
<EMAIL_ADDRESS>Department of Electrical & Computer Engineering, University of
Minnesota, Minneapolis, Minnesota, 55455, USA

###### Abstract

Surface plasmons, which allow extreme confinement of light, suffer from high
intrinsic electronic losses. It has been shown that stimulated emission of
electrons can transfer energy to plasmons and compensate for the high
intrinsic losses. To-date, these realizations have relied on introducing an
external gain media coupled to the surface plasmon. Here, we propose that
plasmons in two-dimensional materials with closely located electron and hole
Fermi pockets can experience gain, when an electrical current bias is applied
along the displaced electron-hole pockets, without the need for an external
gain media. As a prototypical example, we consider WTe2 from the family of
1T′-MX2 materials, whose electronic structure can be described within a type-
II tilted massive Dirac model. We find that the nonlocal plasmonic response
experiences prominent gain for experimentally accessible currents on the order
of mA$\mu$m-1. Furthermore, the group velocity of the plasmon found from the
isofrequency curves imply that the amplified plasmons are highly collimated
along a direction perpendicular to the Dirac node tilt when the electrical
current is applied along it.

##  I Introduction

Surface plasmonsMaier2007  are collective excitations of electrons that result
in highly confined electromagnetic modes in metals or semiconductors. In two-
dimensional materialsLow2017 ; Basov2016 , such as grapheneJablan2009 ;
Low2014b , propagating surface plasmons have been observedFei2012 ; Chen2012 ;
Dai2015  and are accompanied by a variety of interesting phenomena such as
plasmonic waveguidingVakil2011 ; Liu2011a , topological plasmonsJin2017 ;
Jung2018 ; Jin2016b , chiral directional plasmonsSong2016 ; Kumar2016a , and
the plasmonic Fizeau drag effectDong2021a ; Zhao2021 ; Borgnia2015 ;
VanDuppen2016 , among many others. High quality factor surface plasmons in any
two-dimensional material is key to the observation of these near-field
plasmonic phenomena. However, the intrinsically high electronic losses of the
plasmons is a significant barrier to their performance and applicabilityNi2018
; Boltasseva2011 . Hence, it is highly desirable to overcome the fundamental
limitation on quality factor set by the electronic losses. One route to
overcoming this limitation is to introduce enough gain to compensate for the
intrinsic lossesXiao2010a .

Gain, or wave amplification, is a phenomenon that is ubiquitous in nature.
Rogue waves, Rijke tubes, and optical parametric amplification are all
examples in which a wave is amplified rather than damped as it propagates
through a mediumRijke1859 ; RAYLEIGH1878 ; Dysthe2008 ; Solli2007 ;
Baumgartner1979 . In most of these examples, the amplification is driven by
external pumping, such as a heater or laser that can inject energy into the
system. Early studies to impart gain to plasmons began with the realization
that stimulated emission can also amplify plasmons in a way similar to
amplification of lightBergman2003a . The first experimental demonstrations
were able to amplify plasmons by coupling to external media that were either
opticallyOulton2009 ; Noginov2009  or electricallyHill2009  pumped, realizing
plasmonic nanolasers with extremely compact field confinement. Optically
pumped setupsNezhad2010 ; Lu2012  require gain medium to be adjacent to the
surface plasmon mode. Electrically pumped systems encapsulate a typical
semiconductor laser system in conjunction with a plasmonic mode in a metal-
insulator-metal gap, accompanied by population inversion within the insulating
layer. However, integrating these approaches into on-chip nanophotonics has
been challenging due to difficulties in down-sizing the optical gain
components and issues related to contact resistance of small metallic
structuresAzzam2020 ; Ding2013 .

Here, we propose that plasmons in a two-dimensional material with separated
electron and hole pockets can experience internal plasmonic gain in the
presence of an applied electrical current. The plasmonic medium derives its
gain from nonlocal single particle transitions between the electron and hole
bands and does not require an external gain medium. The electrons(holes)
acquire a momentum anti-parallel(parallel) to the applied current direction.
The shifted electron and hole distributions open up a finite energy momentum
phase space window where population inversion can be achieved for single-
particle transitions. When this phase space overlaps with the plasmon
dispersion, plasmonic gain can occur. Such a condition for gain can be
satisfied in two-dimensional materials with low carrier density when the
electron and hole pockets have a separation on the order of the plasmon
momentum. Examples of materials that meet these prerequisites are found in the
1T′-MX2 family of transitional metal dichalcogenides such as WTe2 and
MoTe2Qian2014a ; Soluyanov2015 .

##  II Results

The inverted band structure found in many 1T′-MX2Qian2014a ; Soluyanov2015
materials follows the tilted 2D massive Dirac modelSodemann2015 ; Du2018

$H=t\eta k_{x}+v(k_{y}\sigma_{x}+k_{x}\sigma_{y})+(m/2-\alpha
k^{2})\sigma_{z}$ (1)

where $t$ is the tilt along the $k_{x}$ direction, $\eta=\pm 1$ is the node
index, $v$ is the Fermi velocity, $m$ is the gap size, and $\alpha$ is a term
that closes the Fermi surface at large $k$. This minimal model captures most
of the salient features of the unique inverted band structure. The energy
dispersion of eq. 1 given by $\varepsilon_{\pm}(\mathbf{k})=t\eta
k_{x}\pm[v^{2}k^{2}+(m/2-\alpha k^{2})^{2}]^{1/2}$ shows an interesting
dependence on the tilt angle parameterized by the ratio $t/v$ (see fig. 1a).
For $t/v<1$ the model describes a type-I Dirac node which has a single
electron or hole pocket at Fermi energies in the conduction or valence band.
For $t/v>1$ we find a type-II Dirac node which can support both an electron
and hole Fermi surface where $\hat{\mathbf{k}}_{e-h}$, the relative
displacement vector of the electron and hole pockets, is along the tilt
direction. Here, we will be focusing on the type-II Dirac node with
$E_{F}=0$eV which has electron and hole pockets of equal size. For
definiteness, the parameters of the Dirac model are set to fit the tight-
binding band structure of monolayer WTe2Hu2021a  and are given as
$v=2.427$eVÅ, $t=1.16v$, $m=0.138$eV, $\alpha=7.541$eVÅ2.

Gain in the plasmons can be seen from a nonlocal description of the electron
response. Hence we calculate the loss function
$L(\mathbf{q},\omega)=\textrm{Im}(1/\epsilon(\mathbf{q},\omega))$ which
reveals the plasmon excitations of the system. The dielectric response is
calculated from the random phase approximation to be
$\epsilon(\mathbf{q},\omega)=1-v_{c}\Pi(\mathbf{q},\omega)$ where $v_{c}$ is
the Coulomb potential and $\Pi(\mathbf{q},\omega)$ is the non-interacting
polarizability. The polarizability is given by a sum of the contributions from
two oppositely tilted nodes ($\Pi=\Pi^{+}+\Pi^{-}$) where the polarizability
of each node is given by

$\Pi^{\eta}(\mathbf{q},\omega)=\frac{g_{s}}{4\pi^{2}}\int
d^{2}\mathbf{k}\sum_{ss^{\prime}}\frac{f(E_{s\mathbf{k}}^{\eta})-f(E^{\eta}_{s^{\prime}\mathbf{k}^{\prime}})}{E^{\eta}_{s\mathbf{k}}-E^{\eta}_{s^{\prime}\mathbf{k}^{\prime}}+\hbar\omega+i\delta}|\langle\psi^{\eta}_{s\mathbf{k}}|\psi^{\eta}_{s^{\prime}\mathbf{k}^{\prime}}\rangle|^{2}.$
(2)

Here, $f(E_{\mathbf{k}})$ is the Fermi distribution,
$\mathbf{k}^{\prime}=\mathbf{k}+\mathbf{q}$, $E^{\eta}_{s\mathbf{k}}$ and
$|\psi^{\eta}_{s\mathbf{k}}\rangle$ are the eigenvalues and eigenvectors of
the tilted 2D massive Dirac Hamiltonian, and $g_{s}$ is the spin degeneracy. A
DC current bias $\mathbf{j}$ applied to the system may then be modeled by
modifying the Fermi distribution to
$f_{u}(E_{\mathbf{k}})=\left(e^{(E_{\mathbf{k}}-\mathbf{u}\cdot\mathbf{k}-E_{F})/k_{B}T}+1\right)^{-1}$Bistritzer2009
where $\mathbf{u}$ is a parameter induced by the current density. The modified
carrier occupation in valley $\eta=+1$ for a current applied along the
positive $x$-direction is shown in fig. 1a. The applied electrical current
requires that the electron and hole pockets in the $\eta=+1$($\eta=-1$) node
are shifted closer to(away from) each other (see inset of fig. 1a).

The loss function for a drift current $u/v=0.41$ and
$\mathbf{q}=q_{x}\hat{\mathbf{x}}$ is shown in fig. 1b. We first note that the
plasmon dispersion is skewed in the direction opposite to the drift current
and exhibits a very prominent non-reciprocal behavior. Moreover, the loss
function in the $-q_{x}$ plasmon branch becomes negative signifying gain in
the plasmon. The origin of gain can be understood by examining the imaginary
part of the polarizability from the $\eta=+1$ node shown in fig. 1c. The bias
creates two regions in which the imaginary part of the polarizability is
negative (blue regions labeled intra and inter in fig. 1c). Both these regions
originate from single particle transitions from a higher energy to lower
energy, i.e. an inverted population. The intraband transitions occur mainly at
large wavevectors of order $k_{F}$ and small frequencies that are distant from
the plasmonic spectrum. In contrast, the interband transitions between the
electron and hole pockets occur at much smaller wave vectors and much higher
frequencies. As the bias is increased, the interband gain region can overlap
with the plasmonic spectrum and impart gain. Indeed this can be seen by
examining the boundary for interband transitions (see dashed magenta line in
fig. 1b and c) connecting states above and below $E_{F}=0$ which matches well
with the onset of plasmonic gain as shown in fig. 1b. Therefore we conclude
that the observed plasmonic gain can be attributed to nonlocal interband
transitions.

Note that the plasmonic gain as shown by the loss function in fig. 1a does not
perfectly match the gain phase space suggested by the polarizability in fig.
1c. Specifically, the positive $q$ branch of the plasmon is not experiencing
any gain despite overlapping with the interband gain phase space. This can be
explained by taking into account the transitions from both the $\eta=1$ and
$\eta=-1$ valleys. fig. 1b has 4 regions that are separated by the interband
transition boundaries of each valley. Region 1 has no interband transitions
while regions 2 and 3 have transitions from $\eta=+1$ and $\eta=-1$
respectively. The plasmon branch that resides in region 2 therefore
experiences the gain from single particle interband transitions in the
$\eta=+1$ valley. Plasmons in region 4 are influenced by interband transitions
from both valleys, and since the $\eta=-1$ valley transitions only contribute
to loss (see SI), gain from the $\eta=+1$ valley is cancelled out.

Having shown that plasmonic gain can indeed emerge for the tilted Dirac
system, we now determine the threshold $u_{th}$ required for gain as a
function of the model parameters. Since $\alpha$ closes the Fermi surface and
therefore controls the size of the electron and hole pockets, we explore
$u_{th}$ as a function of $\alpha$ while leaving all other parameters fixed to
the values used for fig. 1. A larger $\alpha$ will in general lead to a
smaller Fermi pocket size. Let us first define the metric for quantifying the
gain in the system as

$g(\alpha,u)=\left|\int_{L<0}L(q,\omega)dqd\omega\right|$ (3)

where the integral is performed only over the region with gain ($L<0$).
Results for $g(\alpha,u)$ are shown in fig. 2a. A clear linear boundary is
observed between regions with and without plasmonic gain (a detailed analysis
of the linear phase boundary is given in the supplementary information). Since
the shift in the Fermi surface is proportional to both the parameter $u$ and
the unbiased $k_{F}$, we can expect that a smaller Fermi surface will require
a larger $u$ to reach the gain threshold. The threshold current density, which
is dependent on both $u$ and the carrier density, is found to be on the order
of a few mA$\mu$m-1 and also decrease as a function of $\alpha$ (see SI). To
take a closer look at what happens at the onset of gain, $g(\alpha,u)$ as a
function of $u$ at $\alpha=8$eVÅ2 is shown in fig. 2b. For $u<u_{th}$, the
slow increase in gain is due to intraband transitions as we can see from the
loss function (inset). For $u>u_{th}$, there is a sharp increase in gain which
originates from the onset of interband transitions imparting gain to the
plasmon. Note the difference in scale between the two loss function insets
which explains why such an abrupt transition is observed for $g(\alpha,u)$. We
now expand our description of the plasmon to a general
$\mathbf{q}=(q_{x},q_{y})$. At a given frequency, the complex plasmon
wavevector
$\mathbf{q}_{pl}=\mathbf{q}_{pl}^{\prime}+i\mathbf{q}_{pl}^{\prime\prime}$ can
be found by solving for the complex roots of $\epsilon(\mathbf{q},\omega)=0$.
For small $\mathbf{q}^{\prime\prime}$, the complex roots are found by solving

$\displaystyle\textrm{Re}\left[\epsilon(\mathbf{q}^{\prime}_{pl})\right]=0$
(4a)
$\displaystyle\mathbf{q}^{\prime\prime}_{pl}\cdot\left[\nabla_{\mathbf{q^{\prime}}}\textrm{Re}(\epsilon)\right]_{\mathbf{q}^{\prime}_{pl}}=-\textrm{Im}\left[\epsilon(\mathbf{q}_{pl}^{\prime})\right].$
(4b)

Note, however, that eq. 4b does not uniquely define the direction of
$\mathbf{q}_{pl}^{\prime\prime}$. $\mathbf{q}_{pl}^{\prime\prime}$ is taken to
be parallel to the group velocity as the gain or loss is expected to be in the
direction of the energy flow. The isofrequency curves for $u/v=0$ and
$u/v=0.41$ are shown in fig. 3. As observed for the loss function in fig. 1, a
finite drift velocity breaks reciprocity and induces gain in the plasmon
spectrum. The isofrequency curves shown in fig. 3c,d further reveal that both
the gain experienced by the plasmon for non-zero drift velocity and the
topology of the isofrequency contour evolve with frequency. We find that at
low frequency the isofrequency contour is separated into two distinct pockets,
and that the amplified plasmons have a group velocity primarily along
$\hat{\mathbf{k}}_{e-h}$. As the frequency is increased above
$\hbar\omega=0.15$eV, these pockets merge into a single closed isofrequency
surface which is relatively ‘flat’ in the region of gain and has a group
velocity perpendicular to $\hat{\mathbf{k}}_{e-h}$. Hence, highly directional
propagation of the amplified plasmons perpendicular to
$\hat{\mathbf{k}}_{e-h}$ is expected.

With the complex plasmon wavevectors found, we now examine the real space
behavior of the plasmons launched by a point source. The plasmon field can be
written as a superposition of unidirectional plane-waves

$\mathbf{E(r)}\propto\sum_{\mathbf{q}_{pl}}e^{i\mathbf{q}^{\prime}_{pl}\cdot\mathbf{r}}e^{-\mathbf{q}^{\prime\prime}_{pl}\cdot\mathbf{r}}\Theta(\mathbf{v}_{g}\cdot\mathbf{r}).$
(5)

$\Theta$ is a step function that describes a one-sided plane wave propagating
in the direction of the group velocity. Results for the plane wave
superposition for $\hbar\omega=0.16$eV with drift velocities $u<u_{th}$ and
$u>u_{th}$ are shown in fig. 4. For both cases, the plasmon propagation
direction is roughly along the $y$-axis as expected from the isofrequency
plots. It is important to emphasize that this is different from what one would
expect from the phase velocity, which lies primarily in the
$-q_{x}$-direction. For $u>u_{th}$, amplification of the plasmon in the
direction of $\mathbf{v}_{g}$ is clearly observed.

##  III Discussion & Conclusion

In summary, we have demonstrated that plasmonic gain may be achieved by
applying an electric current bias to materials with a type-II Dirac node band
stucture. While we have used parameters that were fitted to the band structure
of WTe2 for most results presented in this work, the gain phenomena observed
is generic and can be anticipated for semimetals with an electron hole pocket
separation on the order of the plasmon wavevector. General requirements for
observing gain are as follows. The material must have electron and hole
pockets with a momentum space separation on the order of the plasmon wave
vector and the current must be applied along the direction of separation. As a
result, gain can be imparted to plasmons when its dispersion overlaps with the
current driven population inverted single particle transition phase space. We
have shown in fig. 1 that for the type-II Dirac node, while both intra- and
interband transitions allow for plasmon emission, the intraband transitions
are well separated from the plasmon dispersion in the $(\mathbf{q},\omega)$
space. The interband transitions on the other hand overlap with the plasmon
dispersion. Also, recall that contributions from both the $\eta=+1$ and
$\eta=-1$ valleys must be included. A given plasmon emission process
$(\mathbf{q},\omega)$ in valley $\eta=+1$ may correspond to plasmon absorption
in valley $\eta=-1$ which overall cancels out the effect of gain to the
plasmon. Hence the emission process that imparts gain to the plasmon must not
have an absorbing counterpart in the opposite valley. When these conditions
are met, we have shown that internal plasmonic gain may be imparted, thus
providing a convenient route towards overcoming the intrinsic loss in the
plasmon. The current density required to achieve gain is found to be on the
order of a few mA$\mu m^{-1}$ which is achievable in monolayer
grapheneSon2018a  and thin film WTe2Mleczko2016 . Furthermore, we find that
the directionality of the plasmon propagation becomes highly dependent on
frequency and current magnitude. Thus our proposed setup presents a simple
pathway for overcoming plasmonic loss for tightly confined two-dimensional
plasmons while also providing control over their propagation direction.
Propagation of these plasmons should be readily observable using a typical
near-field scanning optical microscope setupDong2021a .

###### Acknowledgements.

All authors acknowledge support by the National Science Foundation, NSF/EFRI
Grant No. EFRI-1741660.

Figure 1: Tilted 2D massive Dirac model and loss function calculation. a, A
representative energy band structure of eq. 1 for the $\eta=+1$ node with a
schematic representation of intraband and interband single particle
transitions. Parameters fitted to the WTe2 tight-binding band structure are
$v=2.427$eVÅ, $t=1.16v$, $m=0.138$eV, $\alpha=7.541$eVÅ2. We have further set
$E_{F}=0$eV and $u=1$eVÅ. Inset shows the tight-binding band structure of the
XTe2 family where the red boxed region is the region we are modeling by the
Dirac model. b, Loss function $L(q_{x},\omega)$ calculation for parameters
given above. The magenta(black) dashed lines show the interband transition
boundaries for the $\eta=+1$($\eta=-1$) valley. c, Polarizability of the
$\eta=+1$ valley. The two regions with negative $\textrm{Im}\Pi^{+}$ originate
from intra and interband transitions. Figure 2: Threshold drift velocity for
gain. a, The gain metric $g(\alpha,u)$ is shown as a 2D color plot. A clear
linear boundary between regions with and without plasmonic gain is observed.
The red dashed line is a linear fitting to the transition boundary that gives
$\alpha=50.99u/v-6.508$. b, A slice of the 2D color plot at $\alpha=8$eVÅ2
(white dashed line in fig. 2a) is shown. The inset band structures
schematically show the transitions responsible for $u/v=0.165$ and $u/v=0.33$
which are below and above $u_{th}$ respectively. The inset loss functions are
scaled to only represent the regions over which $g(\alpha,u)$ is integrated as
defined by $g(\alpha,u)=\left|\int_{L<0}L(q,\omega)dqd\omega\right|$. A black
dashed line is drawn at $u_{th}/v=0.285$. Figure 3: Isofrequency curves for
the plasmon dispersion. a, Isofrequency curves for $u=0$. All other parameters
for the Dirac model are identical to parameters used in fig. 1. Color scale
quantifies the quality factor $q_{1}/q_{2}$ which gives the number of cycles
the plasmon goes through before decaying or amplifying by $e$. c, Isofrequency
curves for $u=1$eVÅ($u/v=0.41$). Negative values of the quality factor denote
regions of the plasmon dispersion that experience gain. b,d, Top view of the
isofrequency curves shown in fig. 3a,c. Arrow direction shows direction of
increasing frequency.  Figure 4: Fields excited by a point dipole. a,b,
Isofrequency curves at $\hbar\omega=0.16$eV for drift velocity $u/v=0.37$ and
$u/v=0.41$ respectively. c,d, Fields generated by eq. 5 for isofrequency
curves fig. 4a and b respectively. The solid white inset plots show the field
amplitude along the white dashed lines and is normalized to the field at the
origin.

## References

  * (1) Maier, S. A.  _Plasmonics: Fundamentals and Applications_ (Springer US, New York, NY, 2007).
  * (2) Low, T. _et al._ Polaritons in layered two-dimensional materials.  _Nat. Mater._ 16, 182–194 (2017).
  * (3) Basov, D. N., Fogler, M. M. & García De Abajo, F. J.  Polaritons in van der Waals materials.  _Science_ 354 (2016).
  * (4) Jablan, M., Buljan, H. & Soljačić, M.  Plasmonics in graphene at infrared frequencies.  _Phys. Rev. B_ 80, 245435 (2009).
  * (5) Low, T. & Avouris, P.  Graphene plasmonics for terahertz to mid-infrared applications.  _ACS Nano_ 8, 1086–1101 (2014).
  * (6) Fei, Z. _et al._ Gate-tuning of graphene plasmons revealed by infrared nano-imaging.  _Nature_ 486, 82–85 (2012).
  * (7) Chen, J. _et al._ Optical nano-imaging of gate-tunable graphene plasmons.  _Nature_ 487, 77–81 (2012).
  * (8) Dai, S. _et al._ Graphene on hexagonal boron nitride as a tunable hyperbolic metamaterial.  _Nat. Nanotechnol._ 10, 682–686 (2015).
  * (9) Vakil, A. & Engheta, N.  Transformation Optics Using Graphene.  _Science_ 332, 1291–1294 (2011).
  * (10) Liu, M. _et al._ A graphene-based broadband optical modulator.  _Nature_ 474, 64–67 (2011).
  * (11) Jin, D. _et al._ Infrared Topological Plasmons in Graphene.  _Phys. Rev. Lett._ 118, 245301 (2017).
  * (12) Jung, M., Fan, Z. & Shvets, G.  Midinfrared Plasmonic Valleytronics in Metagate-Tuned Graphene.  _Phys. Rev. Lett._ 121, 86807 (2018).
  * (13) Jin, D. _et al._ Topological magnetoplasmon.  _Nat. Commun._ 7, 13486 (2016).
  * (14) Song, J. C. & Rudner, M. S.  Chiral plasmons without magnetic field.  _Proc. Natl. Acad. Sci. U. S. A._ 113, 4658–4663 (2016).
  * (15) Kumar, A. _et al._ Chiral plasmon in gapped Dirac systems.  _Phys. Rev. B_ 93, 041413 (2016).
  * (16) Dong, Y. _et al._ Fizeau drag in graphene plasmonics.  _Nature_ 594, 513–516 (2021).
  * (17) Zhao, W. _et al._ Efficient Fizeau drag from Dirac electrons in monolayer graphene.  _Nature_ 594, 517–521 (2021).
  * (18) Borgnia, D. S., Phan, T. V. & Levitov, L. S.  Quasi-Relativistic Doppler Effect and Non-Reciprocal Plasmons in Graphene.  _arXiv_ 1512.09044 (2015).
  * (19) Duppen, B. V., Tomadin, A., Grigorenko, A. N. & Polini, M.  Current-induced birefringent absorption and non-reciprocal plasmons in graphene.  _2D Mater._ 3, 015011 (2016).
  * (20) Ni, G. X. _et al._ Fundamental limits to graphene plasmonics.  _Nature_ 557, 530–533 (2018).
  * (21) Boltasseva, A. & Atwater, H. A.  Low-Loss Plasmonic Metamaterials.  _Science_ 331, 290–291 (2011).
  * (22) Xiao, S. _et al._ Loss-free and active optical negative-index metamaterials.  _Nature_ 466, 735–738 (2010).
  * (23) Rijke, P. L.  Notiz über eine neue Art, die in einer an beiden Enden offenen Röhre enthaltene Luft in Schwingungen zu versetzen.  _Ann. der Phys. und Chemie_ 183, 339–343 (1859).
  * (24) Rayleigh.  The Explanation of Certain Acoustical Phenomena.  _Nature_ 18, 319–321 (1878).
  * (25) Dysthe, K., Krogstad, H. E. & Müller, P.  Oceanic rogue waves.  _Annu. Rev. Fluid Mech._ 40, 287–310 (2008).
  * (26) Solli, D. R., Ropers, C., Koonath, P. & Jalali, B.  Optical rogue waves.  _Nature_ 450, 1054–1057 (2007).
  * (27) Baumgartner, R. & Byer, R.  Optical parametric amplification.  _IEEE J. Quantum Electron._ 15, 432–444 (1979).
  * (28) Bergman, D. J. & Stockman, M. I.  Surface Plasmon Amplification by Stimulated Emission of Radiation: Quantum Generation of Coherent Surface Plasmons in Nanosystems.  _Phys. Rev. Lett._ 90, 4 (2003).
  * (29) Oulton, R. F. _et al._ Plasmon lasers at deep subwavelength scale.  _Nature_ 461, 629–632 (2009).
  * (30) Noginov, M. A. _et al._ Demonstration of a spaser-based nanolaser.  _Nature_ 460, 1110–1112 (2009).
  * (31) Hill, M. T. _et al._ Lasing in metal-insulator-metal sub-wavelength plasmonic waveguides.  _Opt. Express_ 17, 11107 (2009).
  * (32) Nezhad, M. P. _et al._ Room-temperature subwavelength metallo-dielectric lasers.  _Nat. Photonics_ 4, 395–399 (2010).
  * (33) Lu, Y.-J. _et al._ Plasmonic Nanolaser Using Epitaxially Grown Silver Film.  _Science_ 337, 450–453 (2012).
  * (34) Azzam, S. I. _et al._ Ten years of spasers and plasmonic nanolasers.  _Light Sci. Appl._ 9, 90 (2020).
  * (35) Ding, K. & Ning, C. Z.  Fabrication challenges of electrical injection metallic cavity semiconductor nanolasers.  _Semicond. Sci. Technol._ 28, 124002 (2013).
  * (36) Qian, X., Liu, J., Fu, L. & Li, J.  Quantum spin Hall effect in two-dimensional transition metal dichalcogenides.  _Science_ 346, 1344–1347 (2014).
  * (37) Soluyanov, A. A. _et al._ Type-II Weyl semimetals.  _Nature_ 527, 495–498 (2015).
  * (38) Sodemann, I. & Fu, L.  Quantum Nonlinear Hall Effect Induced by Berry Curvature Dipole in Time-Reversal Invariant Materials.  _Phys. Rev. Lett._ 115, 216806 (2015).
  * (39) Du, Z. Z., Wang, C. M., Lu, H.-Z. & Xie, X. C.  Band Signatures for Strong Nonlinear Hall Effect in Bilayer WTe2.  _Phys. Rev. Lett._ 121, 266601 (2018).
  * (40) Hu, M., Ma, G., Wan, C. Y. & Liu, J.  Realistic tight-binding model for monolayer transition metal dichalcogenides of 1T’ structure.  _Phys. Rev. B_ 104, 035156 (2021).
  * (41) Bistritzer, R. & MacDonald, A. H.  Hydrodynamic theory of transport in doped graphene.  _Phys. Rev. B_ 80, 085109 (2009).
  * (42) Son, S.-K. _et al._ Graphene hot-electron light bulb: incandescence from hBN-encapsulated graphene in air.  _2D Mater._ 5, 011006 (2017).
  * (43) Mleczko, M. J. _et al._ High Current Density and Low Thermal Conductivity of Atomically Thin Semimetallic WTe2.  _ACS Nano_ 10, 7507–7514 (2016).


# Investigating Power laws in Deep Representation Learning

Arna Ghosh  Arnab Kumar Mondal  Kumar Krishna Agrawal  Blake Richards

###### Abstract

Representation learning that leverages large-scale labelled datasets, is
central to recent progress in machine learning. Access to task relevant labels
at scale is often scarce or expensive, motivating the need to learn from
unlabelled datasets with self-supervised learning (SSL). Such large unlabelled
datasets (with data augmentations) often provide a good coverage of the
underlying input distribution. However evaluating the representations learned
by SSL algorithms still requires task-specific labelled samples in the
training pipeline. Additionally, the generalization of task-specific encoding
is often sensitive to potential distribution shift. Inspired by recent
advances in theoretical machine learning and vision neuroscience, we observe
that the eigenspectrum of the empirical feature covariance matrix often
follows a power law. For visual representations, we estimate the coefficient
of the power law, $\alpha$, across three key attributes which influence
representation learning: learning objective (supervised, SimCLR, Barlow Twins
and BYOL), network architecture (VGG, ResNet and Vision Transformer), and
tasks (object and scene recognition). We observe that under mild conditions,
proximity of $\alpha$ to 1, is strongly correlated to the downstream
generalization performance. Furthermore, $\alpha\approx 1$ is a strong
indicator of robustness to label noise during fine-tuning. Notably, $\alpha$
is computable from the representations without knowledge of any labels,
thereby offering a framework to evaluate the quality of representations in
unlabelled datasets.

Machine Learning, ICML

##  1 Introduction

Representation learning using deep neural networks (DNN) is central to recent
progress in machine learning (ML). The most common approach used for learning
effective representations from visual data is supervised learning, wherein
parameters of a model are trained with labels to optimize performance of a
particular task. However, supervised training requires large-scale human-
annotated dataset, which are often expensive and difficult to collect, thereby
restricting it’s application to a narrow domain. Recent advances in self-
supervised learning (SSL), i.e. learning useful representations without
relying on labels, provide early evidence of task-agnostic representations
learned by imposing structural constraints on the features (Chen et al., 2020;
Tian et al., 2020; Zbontar et al., 2021). As such, evaluating representations
in SSL systems entails computing performance on downstream tasks while fixing
the feature map (usually with linear readouts on intermediate layers). While
SSL algorithms learn these representations without labelled data, assessing
the quality is dependent on access to downstream labels. This motivates a
natural question, _is it possible to assess the quality of visual
representations in deep neural networks without labels?_

Figure 1: Analysing large populations of neural activation in response to
visual stimuli (Stringer et al., 2019), suggests an eigenspectrum with
$\lambda_{i}\propto i^{-\alpha}$, with coefficent $\alpha\approx 1$. For
intermediate representations learned by DNNs, we evaluate the sample
covariance matrix
$\Sigma_{N}(\mathbf{f})=\frac{1}{N}\sum_{i=1}^{N}\mathbf{f}_{\theta}(x_{i})\mathbf{f}_{\theta}(x_{i})^{\top}$
and investigate the eigenspectrum of $\Sigma_{N}{\mathbf{(f)}}$.

To answer this question, we must first define desirable characteristics of
“good” representations. Arguably, the most important attribute of a good
representation is its usability for a wide range of downstream tasks. In other
words, it would be ideal if we could somehow determine whether a given mapping
for representations will allow us to achieve high levels of performance even
when we expose the network to different data distributions.

Notably, biological brains possess representations that are useful for a
multitude of downstream tasks. What are the properties of biological neural
representations that may permit this general usability? Recent findings in
vision neuroscience have revealed an interesting property of learned
representations in the brain (Stringer et al., 2019) that may point to the
answer. Specifically, representations in primary visual cortex (V1 in mice)
were found to obey a power law, i.e. the variance explained by the $n^{th}$
principal component of the covariance matrix scaled as $1/n$. In recent work
(Kong et al., 2022) observe similar structure in representations from V1 in
macaque monkeys and relate robustness to eigenspectrums of neural activations.
Put another way, the coefficient of the power law (often denoted $\alpha$),
has been observed close to 1 in biological brains. Concurrent work in DNNs has
shown that enforcing $\alpha\approx 1$ in representations makes them more
robust to adversarial attacks (Nassar et al., 2020). Together, these results
suggest that an $\alpha$ close to 1 may indicate a model that has the
potential to exhibit robustness to noise and good generalization.

We explore this idea and examine the extent to which an $\alpha$ value close
to 1 correlates with downstream performance on new data distributions and
finetuning under loss functions. We show that across a variety of DNN
architectures, learning objectives, and tasks, models that exhibit good
generalization performance exhibit $\alpha$ values close to 1. Furthermore, we
show that $\alpha$ values close to 1 are also correlated with good transfer
learning performance. These results suggest that proximity of $\alpha$ to 1 is
a potential measure of representation quality. Importantly, this measure can
be calculated without any labels. As such, it provides a means to assess SSL
even when labelled datasets are unavailable.

###  1.1 Related Work

Evaluating representations and model quality We note a substantial body of
work aiming to empirically characterize structure of emergent representations
in DNN, without requiring labels (Nguyen et al., 2020; Raghu et al., 2021).
One such index that quantifies the similarity of representations across layers
(of same or different models) is Centered Kernel Alignment (CKA) (Kornblith et
al., 2019). While CKA does not provide explicit guidance for downstream
performance, (Martin et al., 2021) show that in-distribution generalization
gap can be predicted using a different index based on the model’s parameters.
In particular, they show that Empirical Spectral Density (ESD) of weight
matrices for many DNNs obey a power-law, with the coefficient of decay being
predictive of in-distribution performance. In the present work, we explore
similar indices that potentially correlate with out-of-distribution
generalization by examining the eigenspectrum of _activations_.

Generalization in Overparameterized Models Modern neural networks often have
significantly more parameters than number of training samples, challenging
classical understanding of the bias-variance tradeoff. Overparameterization
permits neural networks to overfit to noise in training data, without
impairing their generalization to unseen data. In recent work, Bubeck et al.
proved that for the interpolator (function implemented by the neural network)
to be smooth, overparameterization is a necessary condition (Bubeck & Sellke,
2021).

Furthermore, this benign overfitting phenomenon in an overparameterized linear
regression problem has been linked to the power law coefficient of the input
covariance matrix (Bartlett et al., 2020). Specifically, Bartlett et al.
showed that for an infinite-dimensional linear regression problem, benign
overfitting is possible iff the eigenspectrum satisfies a power law (upto
polylog factors). More recently, Lee et al. found that the tail eigenvalues of
infinite-width network kernels exhibit a power law decay (Lee et al., 2020).
Following this, Tripuraneni et al. explored high-dimensional random feature
regression settings and analytically showed a dependence between eigenspectrum
decay rate of the feature covariance matrix and generalization error
(Tripuraneni et al., 2021). While these characterizations provide a
theoretical understanding of generalization error in the asymptotic or random
feature settings, corresponding questions in the finite dimensional DNN
trained with gradient descent are open problems.

For deep linear networks trained using gradient descent, the eigenvalues of
input covariance determine the generalization error dynamics (Advani et al.,
2020). Advani et al. demonstrated that small eigenvalues determine the
convergence of training dynamics as well as the overfitting error at
convergence. For overparameterized 2-layer neural networks, Arora et al.
provided a fine grained analysis of generalization bounds (Arora et al.,
2019). In contrast, we study modern DNN architectures and explore the
covariance structure of their learned features on visual recognition tasks.

Main Contributions Our core contributions include:

  1. 1.

establishing that across architectures (VGG, ResNet, ViT), intermediate layers
exhibit representations where the eigenspectrum follows a power law.

  2. 2.

empirical verification that the proximity of $\alpha$ to 1 is strongly
correlated with downstream performance across multiple key attributes
(backbone architecture, pretraining objective & downstream task) which
influence representation learning.

  3. 3.

demonstrating that self-supervised learning algorithms learn representations
that are robust to finetuning under noisy labels. Additionally, improvement in
performance across finetuning correlates strongly with the proximity of
$\alpha$ to 1.

##  2 Preliminaries

We are interested in evaluating quality of features learned by modern neural
network architectures, especially in the high-dimensional space of visual
representations. Formally, we consider inputs
$\mathbf{x}_{i}\in\mathcal{X}\subset\mathbb{R}^{d}$, and learn mappings
$\mathbf{f}:\mathbb{R}^{d}\to\mathbb{R}^{D}$ such that
$\mathbf{f}(\mathbf{x})$ is a vector of $D$-dimensional features. For
instance, sample images from the ImageNet dataset have dimensions
$(3,224,224)$, i.e $d\approx 0.15M$.

We consider DNNs as our function approximators, where each architecture
implicitly defines a function class
$\mathcal{F}=\\{\mathbf{f}_{\theta}:\theta\in\Theta\\}$ where
$\Theta\subset\mathbb{R}^{p}$ is the feasible set of model parameters (e.g
bounded $\Theta,{[-B,B]}^{p}$ for some $B\in\mathbb{R})$. The search for
_good_ representations usually poses the following optimization problem:

_With a dataset $\mathcal{D}_{\mathrm{pretrain}}$ from some data distribution
$\mathbb{P}_{\mathrm{pretrain}}$ (potentially with labels), search for
_optimal_ parameters $\theta^{*}$ such that_

$\displaystyle\theta^{*}=\arg\min_{\theta\in\Theta}\mathcal{L}(\mathbf{f}_{\theta},D_{\mathrm{pretrain}})$
(1)

_where $\mathcal{L}(\mathbf{f},D)$ is the learning objective._

The above optimization problem is usually non-convex, and often use gradient
based optimizers to find an approximate solution $\hat{\theta}$. With
$\hat{\theta}$, we evaluate the quality of representations on a downstream
task with $\mathcal{D}_{\mathcal{T}}$ of $N$ samples from
$\mathbb{P}_{\mathcal{T}}$.

$\displaystyle
R(\hat{\theta},\mathcal{T})=\mathcal{L}(\mathbf{f}_{\hat{\theta}},D_{\mathcal{T}})$
(2)

A concrete example is the pretraining on ImageNet dataset, where we use a
pretrained VGG16 model to extract features from images. These representations
are finetuned with a linear readout on MIT-67 dataset.

###  2.1 Covariance estimation and eigenspectrum

For parameterized functions $\mathbf{f}_{\theta}:\mathcal{X}\to\mathbb{R}^{D}$
, the empirical feature covariance matrix, $\Sigma_{N}(\mathbf{f_{\theta}})$
is defined as

$\displaystyle\Sigma_{N}(\mathbf{f_{\theta}})$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbf{f}_{\theta}(x_{i})\mathbf{f}_{\theta}(x_{i})^{T}$

The eigenspectrum of $\Sigma_{N}(\mathbf{f_{\theta}})$ informs us about the
variance explained by each principal component of the space spanned by
representations $\mathbf{f_{\theta}(x)}$. Using the spectral decomposition
theorem on symmetric matrices, $\Sigma=U\Lambda U^{T}$, where $\Lambda$ is a
diagonal matrix with nonnegative entries, and $U$ is a matrix whose columns
are the eigenvectors of $\Sigma$. Without loss of generality we assume that
$\lambda_{1}\geq\lambda_{2}...\geq\lambda_{m}$, where $m=\min(N,D)$ is the
rank of $\Sigma_{N}(\mathbf{f_{\theta}})$.

The eigenspectrum of a covariance matrix is said to follow a _power law_
$PL(\alpha)$, or _zeta distribution_ if for
$\lambda_{j}\in[\lambda_{min},\lambda_{max}]$, the eigenvalues $\lambda_{j}$
are all nonnegative, and

$\displaystyle\lambda_{j}\propto j^{-\alpha}$

for some $\alpha>0$. Here, $\alpha$ is the _slope_ of the power law, and is
referred to as the _coefficient of decay_ of the eigenspectrum. Intuitively,
small $\alpha$ (typically $\alpha\leq 1$) suggests a dense encoding, while a
high $\alpha$ (rapid decay) corresponds to a sparse encoding.

###  2.2 Network Architectures

We examine representation learning primarily in the context of visual
downstream tasks. As such, we look at neural network architectures which
follow different design principles, with the shared goal of pushing the
envelope on generalization and robustness.

Classic Convolutional Neural Networks Deep Convolutional Neural Networks
(CNNs) have been instrumental in early successes in image classification. In
particular, (Simonyan & Zisserman, 2014) proposes a class of VGG-X
architectures, characterized by a sequence of convolutional layers, two fully
connected layers and a softmax readout layer.

Residual Networks with Skip Connections Learning with residual connections (He
et al., 2016) propose a breakthrough in the ability to train very deep neural
networks ($\sim$1000 layers) by adding a residual connection to _block of
convolutions_. Recent findings of (Raghu et al., 2021) suggest that early and
later residual blocks learn qualitatively different representations of the
same input. In particular, the early layers are shown to learn more localized
information about the input, while the downstream layers with a larger
receptive field are able to capture more global information.

Vision Transformers Recent advances with self-attention modules in
transformers (Dosovitskiy et al., 2020) propose a different paradigm for
representation learning. With patch-based image embedding and self-attention,
even the early layers have a global receptive field, allowing the model to
store representations at different scales via multiple heads. While the later
layers of the network learn specialization to have globally relevant features,
(Raghu et al., 2021) suggests that the representations learned by the early
layers are indeed across different scales.

###  2.3 Learning objectives

As outlined in Equation 1, models in the same function class are often
optimized under different learning objectives, to incorporate inductive
biases. We broadly differentiate between supervised (with labels) and self-
supervised (without labels) pretraining objectives.

Supervised Learning Since we consider primarily classification tasks, models
trained with this objective learn to predict the labels corresponding to
inputs in $\mathcal{X}$. To recover the maximum likelihood estimator, the
cross entropy loss is often used to measure the error of the model, as
$\mathcal{L}_{\mathrm{sup}}=-\frac{1}{N}\sum_{i=1}^{N}y_{i}\log
p_{\theta}(x_{i})$ where
$p_{\theta}(x_{i})=\mathrm{softmax}(\mathbf{g}_{\phi}(\mathbf{f_{\theta}}(x_{i})))$
denotes the logits and $\mathbf{g}_{\phi}:\mathbb{R}^{D}\to\mathbb{R}^{k}$
denote the readout layer. Here, $k$ is the number of classes.

Self-Supervised Learning A common theme in several recently proposed SSL
algorithms is to impose invariance constraints on the representations, and
rely on different views of the data to learn these embeddings. In our
experiments, we evaluate features pretrained under two classes:

  1. 1.

Dual-Networks Motivated by the invariance of representations under benign
augmentations, the _dual-networks_ learning objective can be summarised as

$\mathcal{L}_{\mathrm{dual}}=\frac{1}{2}\|\mathbf{f}(\mathbf{x}^{A};\theta)-\mathbf{f}(\mathbf{x}^{B};\theta_{0})\|^{2}$
(3)

In particular, we consider SimCLR, BYOL (Chen et al., 2020; Grill et al.,
2020) as pretraining objectives, and evaluate the representations on
downstream task performance.

  2. 2.

Efficient Encoding Inspired by the efficient coding hypothesis (Barlow et al.,
1961), the _Barlow Twins_ learning objective (Zbontar et al., 2021) proposes
imposing a _soft-whitening_ constraint as the learning objective.

$\displaystyle\mathcal{L}_{\mathrm{BT}}$
$\displaystyle=\underset{\text{invariance}}{\underbrace{\sum_{i}(1-\mathcal{C}(\mathbf{f_{\theta}})_{ii})^{2}}}+\lambda\underset{\text{redundancy-
reduction}}{\underbrace{\sum_{i}\sum_{j\neq
i}\mathcal{C}(\mathbf{f_{\theta}})_{ij}^{2}}}$ (4) where
$\displaystyle\mathcal{C}(\mathbf{f_{\theta}})_{ij}=\frac{\sum_{N}\mathbf{f}(x^{A};\theta)_{i}\mathbf{f}(x^{B};\theta)_{j}}{\sqrt{\sum_{N}\mathbf{f}(x^{A};\theta)_{i}^{2}\sum_{N}\mathbf{f}(x^{B};\theta)_{j}^{2}}}$

Notably with sufficient large $\lambda$, the model would impose an $\alpha=0$
constraint on the representations.

##  3 Power-Laws in Deep Representation Learning

Algorithm 1 Computing Eigenspectrum Decay

⬇

# model: encoder network

# layer: layer index

# B: batch size

# N: size of evaluation split

# D: dimensionality of the representations

cov = torch.zeros(D, D)

for batch in eval_dataloader:

# extract features

feats = model.feats(batch,layer=layer) # BxD

# aggregate covariates per batch

cov += torch.mm(feats.T, feats) / N # DxD

# compute eigenspectrum

eigenspectrum = torch.linalg.eigvals(cov)

# find decay coefficient (where $\lambda_{i}\propto i^{-\alpha}$)

alpha = fit_powerlaw(eigenspectrum)

Incorporating structure in deep representation learning (e.g inductive biases
via model architectures) has been extremely effective in reducing
generalization error of learned models across domains. To explain these
empirical findings, recent work in theoretical machine learning attempts to
analyse bounds on generalization error in restricted settings. While we
primarily focus on finite-width neural networks on real-world datasets, we
begin by highlighting generalization bounds in the infinite feature
dimensional regime (Bartlett et al., 2020) for linear regression with Gaussian
features. We extend such properties of the min-norm interpolant to finite-
width linear regression with power-law in covariates, optimized using gradient
descent to motivate our experiments in neural networks. For the rest of the
section, we denote feature maps as $\mathbf{f}_{\theta}(x)\in\mathbb{R}^{D}$
where $\mathbf{f}_{\theta}:\mathcal{X}\to\mathbb{R}^{D}$, and the readout
network as $\mathbf{g}_{\phi}:\mathbb{R}^{D}\to\mathbb{R}^{k}$, where $k$ is
target dimensionality. For simplicity of analysis, we consider linear readouts
unless explicitly mentioned, i.e. $\mathbf{g}_{\phi}(x)=x^{T}\phi$.

###  3.1 $D\to\infty$ regime

The eigenspectrum of the feature covariance matrix sheds light on the relation
between the smoothness of mapping $\mathbf{f}_{\theta}$ and good
generalization in the asymptotic regime when $D\to\infty$. In this setting,
Stringer et al. show that the kernel function associated with
$\mathbf{f}_{\theta}$ is continuous iff $\alpha>1$ (see Thm 3. in
Supplementary of (Stringer et al., 2019) for proof).

######  Lemma 3.1.

(Stringer et al., 2019) Let $\lambda_{1}\geq\lambda_{2}\geq...\geq 0$ be the
eigenspectrum of $\Sigma_{\mathbf{f}}$ and $K(x,x^{\prime})$ be the kernel
function corresponding to the mapping
$\mathbf{f}_{\theta}:\mathcal{X}\to\mathbb{R}^{D}$. If
$\frac{1}{D}\mathbb{E}_{x}[\mathbf{f}_{\theta}(x)^{T}\mathbf{f}_{\theta}(x)]$
is finite, ie. finite variance in the representation space, and $K$ is
continuous, then $\lambda_{n}=o(n^{-1})$.

The proof for Lemma 3.1 entails showing that trace of the continuous kernel
$K$ is the total variance in representation space, which in turn equals trace
of $\Sigma({\mathbf{f}})$, i.e. sum of all eigenvalues, $\lambda_{i}$’s.
Following this equality, it can be shown that for finite variance,
$\lambda_{n}\leq\epsilon\frac{1}{n}$ for all $n\geq N$ and some $\epsilon>0$.
Interestingly, Lemma 3.1 shows a direct relation between $\alpha$ and
smoothness of the representation kernel and offering an insight into why
having slow eigenspectrum decay (or low decay coefficients) could be
pathological.

In more recent work, (Bartlett et al., 2020) studied the linear regression
setting with Gaussian features, and proved that the min-norm solution provides
good generalization performance iff the eigenspectrum of
$\Sigma(\mathbf{f_{\theta}})$ follows a power law (upto polylogarithm factors)
with $\alpha=1$. More formally,

######  Lemma 3.2.

(Bartlett et al., 2020) In a linear regression problem parameterized as
$\hat{y}=\mathbf{f}_{\theta}(x)^{T}\psi$, let the risk of a min-norm solution
be defined as

$R[\psi]:=\mathbb{E}_{x,y}\left[(y-\mathbf{f}_{\theta}(x)^{T}\psi)^{2}-(y-\mathbf{f}_{\theta}(x)^{T}\psi^{*})^{2}\right]$
(5)

where $\psi^{*}$ denotes the optimal parameter vector. If the $n^{th}$
eigenvalue of $\Sigma({\mathbf{f}})$ follows a power law, i.e.
$\lambda_{n}=n^{-\alpha}\ln(n+1)^{-\beta}$, then $R[\psi]$ is small iff
$\alpha=1$ and $\beta>1$

We refer the reader to (Bartlett et al., 2020) for the proof of Lemma 3.2.
Taken together, Lemma 3.1 and Lemma 3.2 offers a normative explanation for
$\alpha$ being close to 1, i.e. for representations to be smooth as well as
offering a basis for achieving competitive performance on downstream tasks.

Figure 2: Overparameterized linear regression, with inputs drawn from Gaussian
distribution with power-law in the covariance matrix. Models are trained with
gradient descent with initialization $\psi_{0}=\mathbf{0}$. Note that
$\alpha>1$ particularly suffers from high train, test MSE loss, with low
generalization error with $\alpha\approx 1$.

###  3.2 Finite dimensional models

The asymptotic regime of infinite width is an excellent framework to study
theoretical properties of DNN representations. Despite having a large number
of parameters, practical DNNs always possess finite dimensional
representations, making it important to investigate the implications of such
results in finite width models. As such, Lemma 3.1 is hinged on the fact that
the total variance in the representations is bounded and equals the sum of
eigenvalues, under continuous kernels. In finite dimensions, the sum of the
eigenvalues will be bounded as long as $\alpha\geq 0$, so $\alpha$ might not
necessarily be indicative of the geometry in the representations. Similarly,
Lemma 3.2 allows for a wide range of eigenvalue sequences alongside powerlaws
with $\alpha>1$(refer to Thm 6 in (Bartlett et al., 2020) for details).

These finite dimensional analogies to Lemma 3.1 and Lemma 3.2 raise the
question: _Does $\alpha$ sufficiently larger than 1, still allow efficient
learning and strong generalizability?_

To answer this question, we narrow our focus to gradient-based optimization
techniques usually used to train DNNs. In particular, from the optimization
perspective, Advani et al. showed that for deep linear regression in high
dimensions, the time required for training and the steady-state generalization
error are both $\mathcal{O}(\frac{1}{\lambda_{min}})$ (Advani et al., 2020). A
key difference from our work is that they assume the inputs are drawn from an
isotropic Gaussian distribution. Instead, we investigate the generalization of
linear regression on $\mathbf{f}_{\theta}(x)$, which has a power law structure
in its covariates. First, in the linear regression setting, Advani et al.
emphasize

######  Lemma 3.3.

Let $\hat{y}=\mathbf{f}_{\theta}(x)^{T}\psi$ be a finite dimensional linear
regression problem where $\psi$ is learned using gradient descent in order to
optimize the training error,
$\mathbb{E}_{x,y}[(y-\mathbf{f}_{\theta}(x)^{T}\psi)^{2}]$, where
$(x,y)\sim\mathcal{D}_{train}$, i.e. the training dataset. If we assume
isotropic Gaussian noise in targets $y$, then the time required by gradient
descent to minimize the training error,
$T_{convergence}=\mathcal{O}(\frac{1}{\lambda_{min}})$

Lemma 3.3 effectively states that small eigenvalues in $\Sigma_{\mathbf{f}}$
impedes training using gradient descent. We extend Lemma 3.3 to explicitly
state training convergence time as the following theorem:

######  Theorem 3.4.

Let $\hat{y}=\mathbf{f}_{\theta}(x)^{T}\psi$ be an overparameterized linear
regression problem where $\psi$ is learned using gradient descent in order to
optimize the training error,
$\mathbb{E}_{x,y}[(y-\mathbf{f}_{\theta}(x)^{T}\psi)^{2}]$, where
$(x,y)\sim\mathcal{D}_{train}$. If we assume power law distribution in
eigenspectrum of representations at $\mathbf{f}_{\theta}$, i.e.
$\lambda_{n}=\frac{c}{n^{\alpha}}\quad\forall n\geq n^{*}$, where
$n^{*}\in\\{1,2...N\\}$, then the time required by gradient descent to
minimize the training error, $T_{convergence}=\mathcal{O}(N^{\alpha})$ where N
is number of training samples.

The outline of the proof builds on key results relating to gradient descent
dynamics from (Shah et al., 2018). We show that gradient descent updates, when
$\psi$ is initialized to 0, yield a recursive relation for $\psi(k)$, i.e.
$\psi$ after $k$ update steps. Plugging this relation in the gradient
formulation, we show that the update step length along the $n^{th}$ principal
direction of $\mathbf{f}_{\theta}(x)$ shrinks exponentially with a decay rate
proportional to $\lambda_{n}$. Therefore, the time to convergence in training
is controlled by the smallest eigenvalue which, by design, follows the power
law. In sum, Theorem 3.4 provides an explanation against arbitrarily large
values of $\alpha$. Taken together, Theorem 3.4 suggests that $\alpha=1$ might
be beneficial, where these representations form basis for gradient-based
optimization on downstream task performance.

An Additional Motivating Example Before exploring the link between
$\alpha\approx 1$ and generalization performance of deep networks, we first
consider its relationship to the finite dimensional regression setting as in
(Bartlett et al., 2020). Specifically, Theorem 6 of their paper states that
the conditions for “benign overfitting”, wherein a model can perfectly fit
noisy training data without any subsequent loss of performance on testing
data, may be looser in finite dimensions as opposed to the necessary and
sufficient condition of $\alpha=1$ in infinite dimensions. To empirically test
this in the finite high dimensional setting, we examine linear least squares
regression using different covariate structures for the input data. Formally,
we consider covariates $\\{x_{i}\\}_{i=1}^{N}$, such that
$x_{i}\in\mathbb{R}^{d}$ is sampled from a Gaussian distribution with
covariance structure $\Sigma=\mathrm{diag}\\{\lambda_{1},...\lambda_{d}\\}$
where $\lambda_{j}\sim PL(\alpha)$, i.e $\lambda_{j}\propto cj^{-\alpha}$. We
assume access to the corresponding labels $\\{y_{i}\\}_{i=1}^{N}$ generated
under a teacher function $\theta^{*}$, such that
$y_{i}=x_{i}^{T}\theta^{*}+\epsilon_{i}$. We find that in this scenario there
is a clear relationship between the proximity of $\alpha$ to 1 and the
presence of benign overfitting. As shown in Figure 2, when $\alpha$ is close
to 1, the training loss is low, but the validation loss is also low. Thus,
when $\alpha$ is close to 1, the generalization properties are at their best.
This example thus provides another hint that $\alpha\approx 1$ is a potential
measure for how well a model will be able to generalize.

In our experiments, we investigate the following questions on vision
classification tasks, when $\mathbf{f}_{\theta}$ is a neural network:

  1. 1.

Representation Quality: How does $\alpha$ vary across backbone architectures
and pretraining learning objectives? Are representations with $\alpha\approx
1$ likely to enjoy better out-of-distribution performance?

  2. 2.

Generalization across tasks: Is $\alpha$ informative of generalization when
evaluated on different downstream tasks?

  3. 3.

Robustness: On finetuning with noisy labels, does task-performance correlate
to $\alpha$?

##  4 Experimental Setup

Driven by this observation in linear regression settings, we investigate
whether $\alpha\approx 1$ is a good indicator of generalization performance in
DNNs. We investigate this relationship between $\alpha$ and downstream task
performance across different network architectures and pretraining loss
functions. In this work, we restrict our focus to visual learning tasks,
specifically object recognition and scene recognition tasks.

In each experiment, we evaluate the $\alpha$ of emergent representations,
$\mathbf{f}(x)$, at intermediate layers of a DNN pretrained on ImageNet (Deng
et al., 2009) by computing the covariance matrix, $\Sigma_{\mathbf{f}}$, and
fitting a power law on the eigenspectrum (refer to Algorithm 1). Our
pretrained models are taken from PyTorch Hub (Paszke et al., 2019) and timm
(Wightman, 2019). Given that we observed no significant difference between the
observed $\alpha$ values in the train and test sets, we refer to this
empirical estimate as the $\alpha$ for the dataset. To estimate the capacity
of intermediate representations in solving the downstream task, we train a
linear readout layer, $\mathbf{g}(.)$, from representations to target logits.
Intuitively, this comes down to establishing a relationship between the
manifold geometry and linear separability of representations (Chung et al.,
2018). Thereafter, we observe the correlation between estimated $\alpha$ and
the linear readout performance. Notably, we also tried non-linear
$\mathbf{g}(.)$ and observed a similar trend in results (see Appendix).

###  4.1 Feature Backbones

In this section, we investigate the relationship between $\alpha$ of the
representation covariance matrix and object recognition performance in DNNs
with different backbones and the role of depth. In order to do so, we examine
varying depth configurations within network architectures across three
generations of models on the STL-10 dataset (Quattoni & Torralba, 2009).
First, we observe deep Convolutional Neural Networks (CNNs) without any
residual connection as our first family of models. Specifically, we choose
three different configurations of VGG-Net (Simonyan & Zisserman, 2014), namely
VGG-13, VGG-16 and VGG-19. We inspect representations that are input to the
dropout and MaxPool layers during the forward pass of the network. Second, we
consider Deep Residual Networks (He et al., 2016) which are widely used in
computer vision. We inspect the representations that are input to each of the
residual blocks as well as the Adaptive average pool in ResNet-13, ResNet-50
and ResNet-101 during their respective forward passes. Finally, owing to the
recent success of transformers in object recognition tasks we consider Vision
Transformers (ViT) (Dosovitskiy et al., 2020) as the third family of models,
namely ViT-Base/8 , ViT-Large/16 and ViT-Huge/14. Unlike VGG and ResNet, we
only look at features corresponding to the [CLS] token in the intermediate
layers because it summarizes the entire input image and is used in practice
for class prediction. For all model architectures, we use the weights obtained
from pretraining on ImageNet.

Figure 3: Performance on STL-10 improves for representations with $\alpha$
approaching 1 across different model architectures. In each model, solid
points correspond to the later layers. Correlation coefficient for the greater
than 1 and less than 1 regimes are $\rho_{>1}=-0.922$, ${}^{*}p<0.05$ and
$\rho_{<1}=0.493$, ${}^{*}p<0.05$ respectively.

Figure 3 illustrates the relation between performance and $\alpha$ for all the
nine architectures across intermediate layer representations, as described
above. We found that, while most intermediate representations in CNNs (with or
without residual connections) exhibit $\alpha<1$, representations in ViTs
mostly exhibit $\alpha>1$. Nevertheless, representations extracted from the
deepest layers of all the models exhibit $\alpha$ value in the proximity of 1,
irrespective of the total depth of each model (see Appendix Figure 9).
Furthermore, the performance on downstream task increases with depth. This is
unsurprising because all networks were trained to perform object recognition
on ImageNet (Deng et al., 2009) and thereby would have leveraged hierarchical
processing to learn features that are tuned towards object recognition.
Surprisingly enough, we observe a strong significant correlation between
$\alpha$ and performance on the STL-10 dataset, i.e. a different data
distribution than the training dataset, across layers and model architectures
($\rho=-0.922$, ${}^{*}p<0.05$ for representations exhibiting $\alpha>1$ and
$\rho=-0.922$, ${}^{*}p<0.05$ for representations exhibiting $\alpha<1$). It
is worth noting here that the correlation was weaker for the earliest layers
of each model. We believe that early layers learn more task invariant features
that reflect the statistics of natural images (Kornblith et al., 2019; Zeiler
& Fergus, 2014) and therefore lack task relevant information in their
representations. Taken together, this observation confirms our hypothesis that
$\alpha$ is a good indicator of out-of-distribution generalization performance
when representations possess task relevant information.

###  4.2 Learning objective

In this section, we first aim to understand how the $\alpha$ value changes
across the layers of a fixed architecture DNN when trained with different
learning objectives. We take a ResNet-50 model (He et al., 2016) pre-trained
using three different SSL algorithms, namely SimCLR (Chen et al., 2020), BYOL
(Grill et al., 2020) and Barlow Twins (Zbontar et al., 2021), and the
supervised learning loss objectives on ImageNet-1k(Deng et al., 2009) dataset.
We use a similar procedure as before to extract representations from the
network and estimate $\alpha$.

Figure 4: The performance using linear readouts on STL-10 is strongly
correlated to proximity of $\alpha$ to 1. ${}^{*}p<0.05$. (All models have
ResNet-50 backbone.)

Similar to results in the previous section, all networks irrespective of the
pretraining loss function, exhibit $\alpha$ closer to 1 in the deeper layers
in contrast to intermediate layers (see Appendix Figure 10). This surprising
result indicates that although the pretraining loss function was different,
representations extracted from deepest layers are reflective of the object
semantics in natural images. Furthermore, Figure 4 illustrates the strong
correlation between $\alpha$ and generalization performance on STL-10 across
all pretraining loss functions. Together with results from the previous
section, we validate our hypothesis that the representations that demonstrate
good out-of-distribution generalization performance are characterized by
$\alpha$ close to 1.

Figure 5: Performance on MIT67 is strongly correlated to proximity of $\alpha$
to 1 across different model architectures, except ResNet models. Correlation
coefficient for the greater than 1 and less than 1 regimes are
$\rho_{>1}=-0.667$, ${}^{*}p<0.05$ and $\rho_{<1}=0.708$, ${}^{*}p<0.05$
Figure 6: For fixed backbone architecture (ResNet-50), models trained with
different learning objectives (supervised, self-supervised). We evaluate
intermediate features with linear probe on MIT67, observing a strong
performance when $\alpha\propto 1$. ${}^{*}p<0.05$.

###  4.3 Scene vs Object Recognition

So far, we have observed that $\alpha$ serves as a good indicator for out-of-
distribution generalization performance, under mild considerations111When
representations contain task relevant information, when we change two key
components of Section 2: $\mathbf{f}_{\theta}$ and $\mathcal{L}$. In this
section, we change the evaluation task from object recognition to scene
recognition. It is known that DNNs perceive objects and scenes differently
based on their architecture (Nguyen et al., 2020). We investigate whether
$\alpha$ is a good measure to observe when the downstream task is different
from the training task. Scene recognition is fundamentally different from the
object classification task as the model needs to focus on global features in
the entire image in contrast to a local features in the image which contains
the object (Oliva & Torralba, 2006). We follow the same procedure as the above
sections, but on the MIT67 indoor scene recognition dataset (Quattoni &
Torralba, 2009).

Figure 5 and Figure 6 demonstrate that the relation between $\alpha$ and out-
of-distribution generalization performance holds for most models on the scene
recognition task. An outlier is the supervised ResNet family of models. This
behavior can be attributed to the extent of local information possessed by
Residual networks, owing to their fully convolutional architecture, which does
not enable the resultant representations to capture information relevant to
scene recognition (Raghu et al., 2021). In other words, representations of
ResNet models lack the task relevant information. Taken together, our results
demonstrate that $\alpha$ can be used as a measure to estimate generalization
performance across all three key elements of DNN models when the
representations possess information relevant to downstream task.

Table 1: Classification accuracies on STL-10, on finetuning with noisy-labels
(15% label-noise) in ResNet-50.

Finetune Acc. Finetune Acc. Objective (noiseless) $\uparrow$ (noise=15%)
$\uparrow$ % Drop $\downarrow$ BarlowTwins 91.81 $\pm$ 1.4 88.45$\pm$ 0.34
$\sim$ 3.36 SimCLR 88.83 $\pm$ 0.18 85.6$\pm$ 0.2 $\sim$ 3.2 Supervised 88.39
$\pm$ 0.41 81.19 $\pm$ 1.01 $\sim$ 8.20

###  4.4 Robustness to finetuning with noisy labels

Ideally, “good” representations should exhibit robustness under finetuning
with noisy labels. To test the robustness of visual representations learned
under different objectives, and its relation to $\alpha$, we extract features
from pretrained models and finetune them by training on a dataset with noisy
labels. We evaluate representations that emerged when training with different
learning objectives, but with a fixed backbone (ResNet-50) and the same
downstream dataset (STL-10). Concretely, we extract features from layer-75 of
the ResNet-50 (ResNet-50/L75) by freezing gradient propagation upto layer-75,
and finetune the later layers for image classification.

In Figure 7, we plot the accuracy against $\alpha$ across multiple epochs of
finetuning, where $\alpha$ is evaluated for the layer just before linear
readouts (post-adaptive pooling). We find that as performance on the
validation set improves across epochs, $\alpha$ approaches 1. When compared to
finetuning with noiseless targets, Table 1 highlights the relative robustness
of self-supervised learning objectives when compared with supervised learning.
This finding supports our hypothesis on variability in $\alpha$ during
finetuning correlates with performance improvement.

Figure 7: Representations from ResNet-50/L75 finetuned for 30 epochs, on
STL-10 dataset with label_noise=15%. Note that with increasing performance,
$\alpha$ approaches 1 from either side of the spectrum (solid points are later
epochs). Consequently, SSL objectives with representations s.t $\alpha$ closer
to 1 compared to supervised training, are more robust (see Table 1).
${}^{*}p<0.05$. Check Figure 13 for details.

##  5 Discussion

Summary Our experiments suggest a strong correlation between decay coefficient
for sample eigenspectrum of representations, $\alpha$, and generalization
performance in tasks central to visual perception. In particular, from Figures
3 and 4 we note that across network architectures and pretraining objectives,
classification accuracy on STL-10 improves as $\alpha$ approaches 1.
Additionally, Table 1 suggests that pretraining with self-supervised learning
objectives provides representations which are relatively robust under noisy
finetuning.

Redundancy in ViT representations Unlike all other models that were
considered, representations from early layers of ViT had a rapid eigenspectrum
decay with $\alpha>1$ (see Figure 3). The transformer architecture has a
notable difference by design, i.e. early layers possess global receptive field
context via self-attention on patch embeddings. Raghu et al. found that early
ViT layers incorporate both local and global information (Raghu et al., 2021).
Based on these insights, one intuitive interpretation of our results is that
the representations have a low effective rank, and encode redundant
information relevant across multiple scales.

Limitations While the role of $\alpha$ and its relationship to generalization
performance is better understood in the asymptotic setting for linear
regression, similar questions in finite dimensional nonlinear models are
unanswered. Moreoever, the complexity of computing eigenvalues scales
$\mathcal{O}(D^{3})$ where $D$ is the dimensionality of the representations.
It is also worth noting that the empirical correlation was weaker for the
earliest layers in each model. We believe that early layers learn more task
invariant features such as corners and edges (Kornblith et al., 2019; Zeiler &
Fergus, 2014), and therefore lack relevant information for downstream task.
Therefore, $\alpha$ might not necessarily be reflective of generalization
performance in poorly-trained models.

Future Directions Learning efficiently at scale from unlabelled datasets poses
an exciting open problem in deep representation learning. We hope this work
opens new perspectives on design of learning objectives and model
architectures to learn task-agnostic features. Like smooth interpolation
(Bubeck & Sellke, 2021), we hope that overparameterization and geometry in
high-dimensions provide hints to a more principled understanding of
generalization in deep neural networks.

## Acknowledgement

The authors would like to thank Zahraa Chorghay and Colleen Gillon for their
aesthetic contribution to the manuscript and figures. This research was
enabled in part by support provided by Mila (mila.quebec/en/) and Compute
Canada (www.computecanada.ca).

## References

  * Advani et al. (2020) Advani, M. S., Saxe, A. M., and Sompolinsky, H.  High-dimensional dynamics of generalization error in neural networks.  _Neural Networks_ , 132:428–446, 2020.
  * Arora et al. (2019) Arora, S., Du, S., Hu, W., Li, Z., and Wang, R.  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  In _International Conference on Machine Learning_ , pp. 322–332. PMLR, 2019.
  * Barlow et al. (1961) Barlow, H. B. et al.  Possible principles underlying the transformation of sensory messages.  _Sensory communication_ , 1(01), 1961.
  * Bartlett et al. (2020) Bartlett, P. L., Long, P. M., Lugosi, G., and Tsigler, A.  Benign overfitting in linear regression.  _Proceedings of the National Academy of Sciences_ , 117(48):30063–30070, 2020.
  * Bubeck & Sellke (2021) Bubeck, S. and Sellke, M.  A universal law of robustness via isoperimetry.  _arXiv preprint arXiv:2105.12806_ , 2021.
  * Chen et al. (2020) Chen, T., Kornblith, S., Norouzi, M., and Hinton, G.  A simple framework for contrastive learning of visual representations.  In _International conference on machine learning_ , pp. 1597–1607. PMLR, 2020.
  * Chung et al. (2018) Chung, S., Lee, D. D., and Sompolinsky, H.  Classification and geometry of general perceptual manifolds.  _Physical Review X_ , 8(3):031003, 2018.
  * Deng et al. (2009) Deng, J., Dong, W., Socher, R., Li, L.-J., Li, K., and Fei-Fei, L.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 248–255, 2009.  doi: 10.1109/CVPR.2009.5206848.
  * Dosovitskiy et al. (2020) Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * Grill et al. (2020) Grill, J.-B., Strub, F., Altché, F., Tallec, C., Richemond, P. H., Buchatskaya, E., Doersch, C., Pires, B. A., Guo, Z. D., Azar, M. G., et al.  Bootstrap your own latent: A new approach to self-supervised learning.  _arXiv preprint arXiv:2006.07733_ , 2020.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Kong et al. (2022) Kong, N. C., Margalit, E., Gardner, J. L., and Norcia, A. M.  Increasing neural network robustness improves match to macaque v1 eigenspectrum, spatial frequency preference and predictivity.  _PLOS Computational Biology_ , 18(1):e1009739, 2022.
  * Kornblith et al. (2019) Kornblith, S., Norouzi, M., Lee, H., and Hinton, G.  Similarity of neural network representations revisited.  In _International Conference on Machine Learning_ , pp. 3519–3529. PMLR, 2019.
  * Lee et al. (2020) Lee, J., Schoenholz, S., Pennington, J., Adlam, B., Xiao, L., Novak, R., and Sohl-Dickstein, J.  Finite versus infinite neural networks: an empirical study.  _Advances in Neural Information Processing Systems_ , 33:15156–15172, 2020.
  * Martin et al. (2021) Martin, C. H., Peng, T. S., and Mahoney, M. W.  Predicting trends in the quality of state-of-the-art neural networks without access to training or testing data.  _Nature Communications_ , 12(1):1–13, 2021.
  * Nassar et al. (2020) Nassar, J., Sokol, P., Chang, S., and Harris, K.  On 1/n neural representation and robustness.  _Advances in Neural Information Processing Systems (NeurIPS)_ , 2020\.
  * Nguyen et al. (2020) Nguyen, T., Raghu, M., and Kornblith, S.  Do wide and deep networks learn the same things? uncovering how neural network representations vary with width and depth.  _arXiv preprint arXiv:2010.15327_ , 2020.
  * Oliva & Torralba (2006) Oliva, A. and Torralba, A.  Chapter 2 building the gist of a scene: the role of global image features in recognition.  _Progress in Brain Research_ , pp. 23–36, 2006.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Quattoni & Torralba (2009) Quattoni, A. and Torralba, A.  Recognizing indoor scenes.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , pp. 413–420, 2009.  doi: 10.1109/CVPR.2009.5206537.
  * Raghu et al. (2021) Raghu, M., Unterthiner, T., Kornblith, S., Zhang, C., and Dosovitskiy, A.  Do vision transformers see like convolutional neural networks?  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Shah et al. (2018) Shah, V., Kyrillidis, A., and Sanghavi, S.  Minimum norm solutions do not always generalize well for over-parameterized problems.  _stat_ , 1050:16, 2018.
  * Simonyan & Zisserman (2014) Simonyan, K. and Zisserman, A.  Very deep convolutional networks for large-scale image recognition.  _arXiv preprint arXiv:1409.1556_ , 2014.
  * Stringer et al. (2019) Stringer, C., Pachitariu, M., Steinmetz, N., Carandini, M., and Harris, K. D.  High-dimensional geometry of population responses in visual cortex.  _Nature_ , 571(7765):361–365, 2019.
  * Tian et al. (2020) Tian, Y., Yu, L., Chen, X., and Ganguli, S.  Understanding self-supervised learning with dual deep networks.  _arXiv preprint arXiv:2010.00578_ , 2020.
  * Tripuraneni et al. (2021) Tripuraneni, N., Adlam, B., and Pennington, J.  Covariate shift in high-dimensional random feature regression.  _arXiv preprint arXiv:2111.08234_ , 2021.
  * Wightman (2019) Wightman, R.  Pytorch image models.  https://github.com/rwightman/pytorch-image-models, 2019.
  * Zbontar et al. (2021) Zbontar, J., Jing, L., Misra, I., LeCun, Y., and Deny, S.  Barlow twins: Self-supervised learning via redundancy reduction.  _arXiv preprint arXiv:2103.03230_ , 2021.
  * Zeiler & Fergus (2014) Zeiler, M. D. and Fergus, R.  Visualizing and understanding convolutional networks.  In _European conference on computer vision_ , pp. 818–833. Springer, 2014.

##  Appendix A Proofs

In this section, we present a formal proof of Theorem 3.4. In order to do so,
we will use a lemma pertaining to iterative expression of the linear
regression parameters over training epochs. This lemma is inspired by the
results presented in (Shah et al., 2018)

######  Lemma A.1.

Let $\hat{y}=x^{T}w$ be a finite dimensional linear regression problem where
$w$ is learned using gradient descent in order to optimize the training error,

$\mathcal{L}=\mathbb{E}_{x,y}[(y-\hat{y})^{2}]=\mathbb{E}_{x,y}[(y-x^{T}w)^{2}]$
(6)

where $(x,y)\sim\mathcal{D}_{train}$, i.e. the training dataset. Then $w_{k}$,
i.e. $w$ after training for $k$ epochs can be written as

$w_{k}=X^{T}(XX^{T})^{-1}\left[I-(I-\eta XX^{T})^{k}\right]Y$ (7)

where $X,Y$ indicate the entire training dataset, i.e.
$X\in\mathbb{R}^{N\times d}$ and $Y\in\mathbb{R}^{N}$

###### Proof.

We start with the gradient of the regression loss function for a defined
training set denoted by $(X,Y)$ in a vectorized notation:

$\displaystyle F(w_{k})$ $\displaystyle=|Y-Xw_{k}|^{2}$
$\displaystyle\implies\nabla F(w_{k})$ $\displaystyle=X^{T}(Xw_{k}-Y)$ (8)

We assume that the weights are initialized at 0, i.e. $w_{0}=0$. Using the
gradient descent update:

$\displaystyle w_{k+1}$ $\displaystyle=w_{k}-\eta\nabla F=w_{k}-\eta
X^{T}(Xw_{k}-Y)=(I-\eta X^{T}X)w_{k}+\eta X^{T}Y$ $\displaystyle w_{1}$
$\displaystyle=(I-\eta X^{T}X)w_{o}+\eta X^{T}Y=\eta X^{T}Y$
$\displaystyle\text{Let }w_{k}$ $\displaystyle=\eta X^{T}u_{k}Y\implies
u_{1}=I$ $\displaystyle w_{2}$ $\displaystyle=\eta X^{T}u_{2}Y$
$\displaystyle=(I-\eta X^{T}X)\eta X^{T}Y+\eta X^{T}Y=\eta X^{T}[(I-\eta
XX^{T})+I]Y$ $\displaystyle\implies u_{2}$ $\displaystyle=(I-\eta XX^{T})+I$
$\displaystyle u_{k}$ $\displaystyle=(I-\eta
XX^{T})u_{k-1}+I=\sum_{i=0}^{k}(I-\eta XX^{T})^{i-1}$
$\displaystyle=(I-(I-\eta XX^{T}))^{-1}(I-(I-\eta
XX^{T}))\sum_{i=0}^{k}(I-\eta XX^{T})^{i-1}$ $\displaystyle=(\eta
XX^{T}))^{-1}(I-(I-\eta XX^{T}))\sum_{i=0}^{k}(I-\eta XX^{T})^{i-1}$
$\displaystyle=\frac{1}{\eta}(XX^{T})^{-1}\sum_{i=0}^{k}[(I-\eta
XX^{T})^{i-1}-(I-\eta XX^{T})^{i}]$
$\displaystyle=\frac{1}{\eta}(XX^{T})^{-1}[I-(I-\eta XX^{T})^{k}]$
$\displaystyle\implies w_{k}$ $\displaystyle=\eta
X^{T}u_{k}Y=X^{T}(XX^{T})^{-1}[I-(I-\eta XX^{T})^{k}]Y$ (10)

∎

This proves Lemma A.1. We will now use this lemma to prove Theorem 3.4. From
this result, we can also write:

$\Delta w_{k}=\eta X^{T}(I-\eta XX^{T})^{k}Y$ (11)

We restate the theorem from the main text. Note that the notations are
simplified here from the theorem statement to improve readability.

######  Theorem A.2.

Let $\hat{Y}=X^{T}w$ be a finite dimensional linear regression problem where
$w$ is learned using gradient descent. If we assume power law distribution in
eigenspectrum of $X$, i.e. $\lambda_{n}=\frac{c}{n^{\alpha}}\forall
n\in\\{1,2...N\\}$, then the time required by gradient descent to minimize the
training error, $T_{convergence}=\mathcal{O}(N^{\alpha})$

###### Proof.

Using result from Lemma A.1, it is clear that the gradient converges to 0 if
$\lambda_{1}<\frac{1}{\eta}$ where $\lambda_{1}$ is the leading eigenvalue of
$XX^{T}$.

Thus, $\eta<\frac{1}{\lambda_{1}}$, i.e. small learning rate setting. So, we
set $\eta=\frac{\hat{\eta}}{\lambda_{1}}$ where $\hat{\eta}<1$. Plugging this
in Equation 11

$\Delta
w_{k}=\frac{\hat{\eta}}{\lambda_{1}}X^{T}(I-\frac{\hat{\eta}}{\lambda_{1}}XX^{T})^{k}Y$
(12)

Let $X=U\wedge^{\frac{1}{2}}V^{T}$ denote the singular value decomposition
(SVD), which implies $XX^{T}=U\wedge U^{T}$. Using the SVD, we get
$(I-\frac{\hat{\eta}}{\lambda_{1}}XX^{T})^{k}=(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge
U^{T})^{k}$. It is worth noting that eigenvalues and eigenvectors of
$(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge U^{T})$ are related to that of
$XX^{T}$ as shown below:

$\displaystyle(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge U^{T})u_{i}$
$\displaystyle=u_{i}-\frac{\hat{\eta}}{\lambda_{1}}(U\wedge U^{T}u_{i})$
$\displaystyle=u_{i}-\frac{\hat{\eta}}{\lambda_{1}}\lambda_{i}u_{i}\quad\quad\text{[Using
$U^{T}U=I$]}$ $\displaystyle\implies(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge
U^{T})u_{i}$
$\displaystyle=(1-\frac{\hat{\eta}}{\lambda_{1}}\lambda_{i})u_{i}$ (13)

Using Equation 13, we can write $(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge
U^{T})$ in the eigendecomposition form as $U\tilde{\wedge}U^{T}$ where
$\tilde{\lambda_{i}}=1-\hat{\eta}\frac{\lambda_{i}}{\lambda_{1}}$. Thus,
$\left(I-\frac{\hat{\eta}}{\lambda_{1}}U\wedge
U^{T}\right)^{k}=U\tilde{\wedge}^{k}U^{T}$. Plugging this in Equation 12, we
get:

$\Delta
w_{k}=\frac{\hat{\eta}}{\lambda_{1}}V\wedge^{\frac{1}{2}}U^{T}U\tilde{\wedge}^{k}U^{T}Y=\frac{\hat{\eta}}{\lambda_{1}}V\wedge^{\frac{1}{2}}\tilde{\wedge}^{k}S$
(14)

where $S=U^{T}Y\in\mathbb{R}^{N}$. For the $i^{th}$ element, we get $\Delta
w_{k}^{(i)}=\frac{\hat{\eta}}{\lambda_{1}}\sum_{j}v_{i,j}\sqrt{\lambda_{j}}\tilde{\lambda_{j}}^{k}S_{j}=\frac{\hat{\eta}}{\lambda_{1}}\sum_{j}v_{i,j}\sqrt{\lambda_{j}}(1-\hat{\eta}\frac{\lambda_{j}}{\lambda_{1}})^{k}S_{j}$.
Since all other factors remain constant across training, i.e. do not change
with $k$, the convergence of gradient descent depends on the factors
$\left(1-\hat{\eta}\frac{\lambda_{j}}{\lambda_{1}}\right)^{k}$. Note that we
define gradient descent to converge when $\Delta w_{k}^{(i)}\approx 0$
$\forall$ $i$. Therefore, the limiting factor that determines rate of
convergence is $(1-\hat{\eta}\frac{\lambda_{j}}{\lambda_{1}})^{k}$, which in
turn is limited by the smallest eigenvalue factor:
$\frac{\lambda_{N}}{\lambda_{1}}$.
Assuming $\frac{\lambda_{N}}{\lambda_{1}}\ll
1\implies\hat{\eta}\frac{\lambda_{N}}{\lambda_{1}}\ll 1$ as
$\hat{\eta}<1\implies(1-\hat{\eta}\frac{\lambda_{N}}{\lambda_{1}})^{k}\approx
1-k\hat{\eta}\frac{\lambda_{N}}{\lambda_{1}}$.
Hence the convergence time,
$k^{*}=\mathcal{O}(\hat{\eta}\frac{\lambda_{1}}{\lambda_{N}})=\mathcal{O}(\frac{\lambda_{1}}{\lambda_{N}})$
If $\lambda_{i}$ follows power law, i.e, $\lambda_{i}=ci^{-\alpha}$ and
$\frac{\lambda_{N}}{\lambda_{1}}=N^{-\alpha}$ then
$k^{*}=\mathcal{O}(N^{\alpha})$ i.e. $k^{*}$ grows exponentially with
$\alpha$. ∎

##  Appendix B Experimental results

Figure 8: Performance of non-linear readout on STL-10 is strongly correlated
to proximity of $\alpha$ to 1 across different pretraining loss functions. The
trends are similar to the performance of a linear readout, as shown in Figure
4. ${}^{*}p<0.05$. Figure 9: $\alpha$ for intermediate layer representations
from different backbone architectures demonstrates the contrasting
representations learned by CNNs and ViT. Figure 10: $\alpha$ for intermediate
layer representations from networks trained using different loss functions
show similar trends. Representations from deeper layers exhibit $\alpha$
closer to 1 as compared to middle layer representations. Figure 11: $\alpha$
for intermediate layer representations from different backbone architectures
in MIT67. Representations learned by ViT is qualitatively different from those
is CNNs both in object and scene recognition datasets. Figure 12: $\alpha$ for
intermediate layer representations from networks trained using different loss
functions show similar trends in MIT67. Representations from deeper layers
exhibit $\alpha$ closer to 1 as compared to intermediate layer
representations. Figure 13: Representations from ResNet-50/L75 finetuned for
30 epochs, on STL-10 dataset with label_noise=15%. Note that with increasing
performance, $\alpha$ approaches 1 from either side of the spectrum (solid
points are later epochs). Consequently, SSL objectives with representations
s.t $\alpha$ closer to 1 compared to supervised training, are more robust (see
Table 1). ${}^{*}p<0.05$


Faculty of Mathematics and Physics, University of Ljubljana, Slovenia Jozef
Stefan Institute, Ljubljana, Slovenia

# Spectroscopy of hadrons with heavy quarks from lattice QCD

Sasa Prelovsek

###### Abstract

Lattice QCD results on hadrons with heavy quarks are briefly reviewed. The
focus is on the spectrum of conventional and exotic hadrons. Structure of
certain conventional hadrons is addressed as well.

##  1 Introduction

Experiments revealed that hadrons with the following minimal contents exist:
mesons $\bar{q}q$, baryons $qqq$, tetraquarks $\bar{q}q\bar{q}q$, pentquarks
$\bar{q}qqqq$ and hybrid mesons $\bar{q}Gq$. The first two sectors correspond
to conventional hadrons, while the last three are referred to as exotic
hadrons. We will briefly review lattice results on some of these states
111Only few references are cited due to the page limit. Other references are
listed in the slides. .

##  2 Hadron spectroscopy with Lattice QCD

The spectrum of hadrons (below, near, or above threshold) is extracted from
the energies $E_{n}$ of QCD eigenstates $|n\rangle$ on a finite and
discretized lattice in Eucledian space-time. The eigen-energies $E_{n}$ are
determined from the time-dependence of two-point correlation functions
$\langle O_{i}(t_{E})O_{j}^{\dagger}(0)\rangle$ $=\sum_{n}\langle
O_{i}|n\rangle~{}e^{-E_{n}t_{E}}\langle n|O^{\dagger}_{j}\rangle~{},$ where
operators $O$ create/annihilate the hadron system with a given quantum number
of interest.

The masses of strongly stable hadrons well below threshold are obtained as
$m\\!=\\!E_{n}|_{\vec{p}=0}$. These have already been determined and agree
well with the experiment, e.g. [1, 2].

Figure 1:  Extracting resonances and near-threshold bound states from one-
channel scattering.

The masses of hadrons near threshold and hadronic resonances have to be
inferred from the scattering of two hadrons $H_{1}H_{2}$, which is encoded in
the scattering amplitude $T(E)$. The simplest example is a one-channel
scattering in partial wave $l$ sketched in Fig. 1. Lüscher has shown that the
energy $E$ of a two-hadron eigenstate in finite volume renders $T(E)$ at that
energy in infinite volume [3]. This relation leads to $T(E)$ for real $E$,
which is then analytically continued to complex energies. A pole in $T(E)$
indicates the presence of a state, while its position renders its mass
$m=\mathrm{Re}(E)$ and the width $\Gamma=-2~{}\mathrm{Im}(E)$. A resonance
corresponds to a pole away from the real axes. A bound state corresponds to a
pole below the threshold: the state is referred to the bound state if the pole
occurs for positive imaginary momenta $p=i|p|$ and a virtual bound state if it
occurs for $p=-i|p|$, where $p$ denotes the magnitude of the 3-momentum in the
center-of-mass frame. The majority of the scattering studies are still
performed at $m_{\pi}>m_{\pi}^{phy}$ and at a single lattice spacing.

Figure 2:  A resonance that decays via two channels $R\to
H_{1}H_{2},~{}H_{1}^{\prime}H_{2}^{\prime}$ has to be inferred from the
scattering of two coupled channels.

The resonances that decay via several strong decay channels $R\to
H_{1}H_{2},~{}H_{1}^{\prime}H_{2}^{\prime},..$ have to be extracted from the
coupled-channel scattering (Fig. 2). The scattering matrix for two coupled
channels contains three unknown functions of energy. It is customary to
parametrize their energy dependence in order to extract them from the eigen-
energies using the Lüscher’s formalism.

The formalism and progress in addressing decays $R\to H_{1}H_{2}H_{2}$ was
reviewed at the last edition of this conference [4].

##  3 Spectroscopy of various hadron sectors

The majority of the discovered exotic hadrons contain heavy quarks as those
are more likely to form quasi-bound states due to small kinetic energies. Most
of them decay strongly and the theoretical challenge to study them increases
with the number of decay channels.

Figure 3:  Poles and masses for a charmonium-like system with $I=0$ from
lattice [5]. LHCb discovery of a $X(3960)$ [6] composed of $\bar{c}s\bar{s}c$.
Possible binding mechanism.

Charmonium-like states $\mathbf{\bar{c}c}$,
$\mathbf{\bar{c}c\bar{q}q^{\prime}}$, $\mathbf{\bar{c}cuud}$

The spectrum of charmonium-like states with $I=0$ was extracted from the
coupled channels $\bar{D}D-\bar{D}_{s}D_{s}$ at $m_{\pi}\simeq 280~{}$MeV [5]
(Fig. 3). All the states except for two (indicated by magenta arrows) appear
to be conventional charmonia $\bar{c}c$ . In addition, two exotic scalar
states are predicted near both thresholds. The heavier one has a large
coupling to $\bar{D}_{s}D_{s}$ and a small coupling to $\bar{D}D$ \- it likely
corresponds to a state $X(3960)$ composed of $\bar{c}s\bar{s}c$ recently
discovered by LHCb [6]. The two additional scalars were not found at
$m_{\pi}\simeq 390~{}$MeV in [7].

A candidate for pentaquark $P_{c}=\bar{c}cud$ with $J^{P}=1/2^{-}$ was found
$6\pm 3~{}$MeV below threshold $\bar{D}\Sigma_{c}$ threshold [8]. It appears
as a bound state in one-channel scattering amplitude $\bar{D}\Sigma_{c}$ (Fig.
4), where the lower-lying channel $J/\psi p$ was omitted.

Figure 4:  $\bar{D}\Sigma_{c}$ scattering in $P_{c}$ channel [8], LHCb
discovery and possible binding mechanism.

Hadrons with two heavy quarks: $\mathbf{QQ\bar{q}\bar{q}^{\prime}}$

The $bb\bar{u}\bar{d}$ and $bb\bar{s}\bar{d}$ tetraquarks with
$J^{P}\\!=\\!1^{+}$ are expected to reside significantly below strong decay
thresholds (Fig. 5). This is a reliable conclusion based on a number of
lattice simulations and model-based calculations (listed in the slides). The
hadron $bb\bar{u}\bar{d}$ with such a deep binding is expected to have a small
size, which indicates the dominance of
$[bb]_{\bar{3}_{c}}^{S=1}[\bar{u}\bar{d}]_{3_{c}}^{S=0}$. So far, this the
only tetraquark where lattice finds a strong support for the dominance of the
diquark antidiquark Fock component.

Figure 5:  Left: The binding energies of $bb\bar{u}\bar{d}$ and
$bb\bar{u}\bar{s}$ with $J^{P}=1^{+}$ from lattice (see references in the
slides). Right: The dependence of the binding energy on $m_{b}$ and $m_{u,d}$
[9, 10].  Figure 6:  The $DD^{*}$ scattering in $T_{cc}$ channel and the pole
locations from lattice [11, 12, 13].

Tetraquarks $QQ\bar{q}\bar{q}^{\prime}$ ($Q\\!=\\!c,b~{},q\\!=\\!u,d,s$) with
other flavors and spin-parities are expected near or above strong decay
thresholds $H_{1}H_{2}$, as suggested also by Fig. 5. In order to prove the
existence of such a state, one needs to extract the scattering matrix and
establish a pole in it as shown in Fig. 1.

The doubly charm tetraquark $T_{cc}=cc\bar{u}\bar{d}$ was discovered by LHCb
just $0.4$ MeV below $DD^{*}$ threshold [14], it has $I\\!=\\!0$ and most
likely $J^{P}\\!=\\!1^{+}$. The lattice results on $DD^{*}$ scattering from
three simulations at different $m_{\pi}$ are shown in Fig. 6. All simulations
find significant attraction; the attraction decreases with increasing
$m_{\pi}$. This implies that $T_{cc}$ would-be bound state from experiment
converts to virtual bound state or a resonance at heavier $m_{\pi}$. Virtual
bound state poles are indeed found in simulations [11, 12, 13] when assuming
effective range approximation. Relaxing this approximation and taking into
account the effect of left-hand cut from one-pion exchange leads to a pair of
virtual bound states or a resonance [15]. The long-range potential is
dominated by $\pi\pi$ exchange in [11]. The dominant attraction in this
channel is attributed to $\rho$ exchange in [13].

Hadrons with a single heavy quark

The charmed scalar mesons would form a $SU(3)$ flavor triplet in Fig. 7a
according to the quark model. However, a new paradigm is supported by the
effective field theories based on HQET and ChPT, combined with the lattice
results as well as the experimental data (see slides for many references).
According to this paradigm, the spectrum features $c\bar{q}$ as well as
$c\bar{q}~{}\bar{q}q$ Fock components ($q=u,d,s$). The latter decomposes to
the multiplets $\bar{3}\oplus 6\oplus 5$ in the $SU(3)$ flavor limit. The
attractive interactions within the anti-triplet and the sextet suggest the
existence of hadrons with flavors indicated by cyan circles in Fig. 7b, with
two pair of poles for $I=1/2$ charmed mesons. The lower one resides at
$2.1\\!-\\!2.2~{}$GeV in agreement with the lattice simulations and it is a
natural partner of $D_{s0}^{*}(2317)$. The heavier pole at
$2.4\\!-\\!2.5~{}$GeV is suggested by the EFT re-analysis of the lattice and
experimental data.

The simulation of $D^{*}\pi$ scattering finds one axial charmed meson
dominated by s-wave and the other by d-wave coupling to $D^{*}\pi$ [16] (Fig.
8). This is in line with HQET prediction and significantly different
experimental decay widths.

Figure 7:  Scalar charmed mesons according to the quark model (a) and the new
paradigm (b).  Figure 8:  Poles related to charmed mesons with
$J^{P}=1^{+},2^{+}$ from [16].

Bottomonium-like states $\mathbf{\bar{b}b}$, $\mathbf{\bar{b}Gb}$,
$\mathbf{\bar{b}b\bar{q}q^{\prime}}$

Figure 9:  (a,b) The static potential of the quarkonium system that accounts
also for the coupling to a pair of heavy mesons [17]. (c,d) Poles and
composition of the bottomonium-like states [18] from analogous static
potential [19].  Figure 10:  The attractive static potential between $B$ and
$\bar{B}^{*}$ [20] is likely responsible for the existence of
$Z_{b}\simeq\bar{b}b\bar{d}u$ [21].  Figure 11:  The potentials related to
hybrids from [22]. Masses of $\bar{b}Gb$ hybrids (update of [23, 24]) from
older potentials [25] (blue) and from relativistic $b$ quarks [26] (violet).

These systems can be studied using relativistic, non-relativistic or static
b-quarks. Here I focus on the last option, where two heavy quarks and
additional light degrees of freedom are investigated via the Born-Oppenheimer
approximation. The eigen-energies at fixed distance between heavy quarks
render the potential $V(r)$. Motion of heavy quarks within this potential is
then studied with Schrödinger-type equation. The aim is to determine whether
bound states or resonances exist.

The static potential for $\bar{b}b$ system with $I=0$ that accounts also for
the coupling to a pair of heavy mesons [17] is shown in Fig. 9b. Coupled-
channel Schrödinger equation with analogous potential (from an earlier
calculation [19]) renders the poles related to bottomonium-like states and
their composition in terms of $\bar{b}b$, $\bar{B}B$ and $\bar{B}_{s}B_{s}$
[18] (Figs. 9c,d). These are dominated by the conventional Fock component
$\bar{b}b$, except for state $n\\!=\\!5$, which is likely related to
unconventional $\Upsilon(10750)$ discovered by Belle.

The observed $Z_{b}$ resonances with flavor content $\bar{b}b\bar{d}u$ are
challenging for rigorous treatment since the lowest decay channel is
$\Upsilon_{b}\pi$, while they reside at the higher threshold $B\bar{B}^{*}$.
This was taken into account in the extraction of the potential between $B$ and
$\bar{B}^{*}$ in Fig. 10, which is attractive at small distances [27, 20, 28].
This attraction is likely responsible for the existence of the exotic $Z_{b}$.

The excited potentials for $\bar{b}b$ with certain spin-parities in Fig. 11
are relevant for hybrid mesons $\bar{b}Gb$. The masses from these potentials
within the Born-Oppenheimer approach [23, 24] agree with those obtained using
relativistic $b$-quarks [26].

Figure 12: Probing structure of hadrons with charm quarks in simulations [29,
2, 30].

##  4 Structure of hadrons with heavy quarks

The charge distribution of charmonia was probed via EM and other currents in
[29], while various decay constants that probe the wave function were
extracted in [2]. The mass decomposition with respect to various terms of the
QCD Hamiltonian was determined for charmed baryons in [30] (Fig. 12).

##  5 Looking ahead

The drawback of spectroscopy on Eucledian lattice is suppressed contribution
$e^{-E_{n}t_{E}}$ of the excited states. The evolution of the QCD systems in
real Minkovsky time $e^{-iE_{n}t_{M}}$ is one of the long-term goals for
quantum computers. Such time evolution for tetraquark and pentquark systems
has already been studied in one-dimensional QCD with one quark flavor on a
quantum computer [31] (Fig. 13).

Figure 13:  Tetra- and penta-quarks in one-dimensional QCD from quantum
computer [31].

##  6 Conclusions

Experiments have provided great discoveries of new conventional as well as
around thirty exotic hadrons. I have reviewed the theoretical challenge to
understand the spectroscopic properties of various hadron sectors from lattice
QCD. This approach renders masses of hadrons that are strongly stable, as well
as most of the hadrons that are slightly below the strong decay threshold or
decay strongly via one decay channel. The theoretical challenge increases with
the number of open decay channels. It seems impossible to address the high-
lying states like $Z_{c}(4430)$ with current lattice methods, while many
interesting physics conclusions are already available for certain lower-lying
states.

Acknowledgments

I acknowledge the support from ARRS research core funding No. P1-0035.

## References

  * [1] D. Hatton et al., HPQCD coll., Phys. Rev. D 103, 054512 (2021), 2101.08103
  * [2] J. Koponen et al., HPQCD coll., Mex. Fis. Suppl. 3, 0308018 (2022), 2204.02137
  * [3] M. Lüscher, Nucl. Phys. B 364, 237 (1991)
  * [4] F. Romero-López, PoS LATTICE2022, 235 (2023), 2212.13793
  * [5] S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, S. Piemonte, JHEP (2021), 2011.02542
  * [6] Aaij, Roel et al, LHCb, PRL131 (2022) 071901, 2210.15153
  * [7] D.J. Wilson, C.E. Thomas, J.J. Dudek, R.G. Edwards (2023), 2309.14071
  * [8] H. Xing, J. Liang, L. Liu, P. Sun, Y.B. Yang (2022), 2210.08555
  * [9] A. Francis et al., Phys. Rev. D 99, 054505 (2019), 1810.10550
  * [10] B. Colquhoun et al., PoS LATTICE2021, 144 (2022)
  * [11] Y. Lyu, S. Aoki, T. Doi, T. Hatsuda, Y. Ikeda, J. Meng (2023), 2302.04505
  * [12] M. Padmanath, S. Prelovsek, Phys. Rev. Lett. 129, 032002 (2022), 2202.10110
  * [13] S. Chen et al., Phys. Lett. B 833, 137391 (2022), 2206.06185
  * [14] R. Aaij et al. (LHCb), Nature Phys. 18, 751 (2022), 2109.01038
  * [15] M.L. Du, A. Filin, V. Baru, X.K. Dong, E. Epelbaum, F.K. Guo, C. Hanhart, A. Nefediev, J. Nieves, Q. Wang, Phys. Rev. Lett. 131, 131903 (2023), 2303.09441
  * [16] N. Lang, D.J. Wilson (Hadron Spectrum), PRL 129, 252001 (2022), 2205.05026
  * [17] J. Bulava et al., Phys. Lett. B 793, 493 (2019), 1902.04006
  * [18] P. Bicudo, N. Cardoso, L. Mueller, M. Wagner (2022), 2205.11475
  * [19] G.S. Bali et al., Phys. Rev. D 71, 114513 (2005), hep-lat/0505012
  * [20] S. Prelovsek, H. Bahtiyar, J. Petkovic, Phys. Lett. B 805, 135467 (2020), 1912.02656
  * [21] A. Garmash et al. (Belle), Phys. Rev. Lett. 116, 212001 (2016), 1512.07419
  * [22] C. Schlosser, M. Wagner, Phys. Rev. D 105, 054503 (2022), 2111.00741
  * [23] N. Brambilla et al., Phys. Rev. D99, 014017 (2019), 1805.07713
  * [24] N. Brambilla, W.K. Lai, J. Segovia, J. Castellà, PRD 101, 054040 (2020), 1908.11699
  * [25] K.J. Juge, J. Kuti, C.J. Morningstar, NPB Proc. Suppl. 63, 326 (1998), hep-lat/9709131
  * [26] S.M. Ryan, D.J. Wilson (Hadron Spectrum), JHEP 02, 214 (2021), 2008.02656
  * [27] A. Peters, P. Bicudo, K. Cichy, M. Wagner, 1602.07621
  * [28] M. Sadl, S. Prelovsek, Phys. Rev. D 104, 114503 (2021), 2109.08560
  * [29] J. Delaney, C.E. Thomas, S.M. Ryan (2023), 2301.08213
  * [30] J.B. Li, L.C. Gui, W. Sun, J. Liang, W. Qin (2022), 2211.04713
  * [31] Y.Y. Atas et al., Phys. Rev. Res. 5, 033184 (2023), 2207.03473