Dark Matter Bound States:
A Window into the Early Universe

Juri Smirnov1,

1 Department of Mathematical Sciences, University of Liverpool, Liverpool, L69
7ZL, United Kingdom

<EMAIL_ADDRESS>
14th International Conference on Identification of Dark Matter
Vienna, Austria, 18-22 July 2022 10.21468/SciPostPhysProc.?

## Abstract

Despite significant experimental sensitivity to point-like, weakly interacting
particles at the electroweak mass scale, dark matter has not been found yet.
This could hint at a more complex dark sector with multiple states or
composite dark particles, much like the complexity of the standard model.
Currently, our experimental sensitivity to such scenarios is limited by a lack
of theoretical understanding. The investigation of Dark Matter systems with
unstable or stable bound states provides a playground that leads to the
development of important tools needed for the understanding of dark sectors
with complex phenomena and allows to experimentally test well-motivated DM
scenarios.

###### Contents

  1. 1 Introduction
  2. 2 Unstable Dark Matter Bound States
    1. 2.1 Effects on the Dark Matter Abundance
    2. 2.2 New Signals from late-time Bound-State formation
  3. 3 Stable Dark Matter Bound States
    1. 3.1 Freezeout in Confining Dark Sectors
    2. 3.2 New Experimental Signatures
  4. 4 Conclusion

##  1 Introduction

The particle nature of dark matter is one of the foremost, central problems of
high-energy physics. But despite overwhelming astrophysical evidence for its
existence, key properties, such as its interaction strength, and its precise
impact on structure formation remain unknown. Similarly, the production
mechanism of dark matter is a longstanding, central problem of cosmology.
Importantly, in specific particle-physics scenarios, these two problems are
tightly linked. Thus, experimental discovery and mass determination would
provide a new building block for particle physics, pointing toward new
structures, symmetries, and energy scales that lie beyond the standard model,
while also providing deep new insights into how to model the evolution of the
universe.

Why, despite unprecedented experimental sensitivity — for example with LZ [1,
2], and FermiLAT [3] —, have we failed to detect this new particle? The
assumption that dark matter (DM) is a single, elementary particle with an
electroweak scale mass (e.g., a vanilla WIMP, or weakly interacting massive
particle) may be too simple. The dark sector could be complex, just like the
standard model (SM). It could feature new strong dynamics, number-changing
interactions, multiple stable dark states, and composite objects. In such
scenarios, the stable relics of the dark sector would be produced very
differently from the simplest thermal relics, and feature very different
properties at present times, changing how we should search for them. Even for
TeV-scale WIMPs, self-interaction effects can lead to multiple bound states
with severe effects on production and phenomenology.

##  2 Unstable Dark Matter Bound States

The WIMP remains one of the most compelling DM scenarios, and its parameter
space above the GeV scale is wide open [4]. It would be negligent to not fully
explore it experimentally. As discussed in Ref. [5], this classic scenario can
be subject to complex phenomena, such as bound-state formation. This
possibility is still under-explored, understanding those effects will strongly
improve experimental sensitivity, and help attack the WIMP window from the
high-mass end [6, 7].

Studies of a $U(1)$ toy mode provide an interesting example to understand in
which situations bound-state formation during freezeout will become relevant
[8, 9, 10]. However, in a concrete model of weakly interacting massive
particles (WIMPs), the situation can become more complex, especially in the
case of non-Abelian force mediators. An important group of publications
investigated bound-state formation effects in the simplest models of WIMPs.
For example, the electroweak triplet was studied in Refs. [11, 12]. Further
studies of thermal effects in this framework have been performed in [13, 14,
15, 16].

In terms of searches, experiments such as the Fermi-LAT, the Cherenkov
telescope H.E.S.S [17] as well as the planned CTA [18], and LHASSO [19]
experiments, will soon be dominating the search for thermally produced WIMP-
like dark matter. Note that also neutrino telescopes, such as Antares and
IceCube [20] will gain in sensitivity and begin to place severe bounds on dark
matter annihilation into neutrino final states.

A particular advantage of all those experiments is, that they rely on exactly
the same process which is responsible for the thermal relic production in the
early universe and thus provide the best way to address the WIMP-like thermal
production hypothesis. As will be discussed below the study of dark matter
bound-state formation will supply the needed theoretical framework to allow
the indirect detection experiments to progress faster into the open parameter
space.

Figure 1: Electron positron annihilation exemplifies non-perturbative effects
at low velocities. Top to bottom: tree-level annihilation, annihilation with
Sommerfeld enhancement, and annihilation via positronium formation. Circle
size schematically depicts the relative size of the cross sections.

###  2.1 Effects on the Dark Matter Abundance

An important input for the experimental searches is the target annihilation
reaction rate, provided by the production process. In a detailed study,
important methodologies have been developed [5], that can be used to describe
dark-matter bound-state formation for arbitrary gauge groups. The developed
formalism relies on the non-relativistic expansion and the dipole
approximation, which is well applicable to all weakly coupled models. Further
studies have been performed by taking into account finite temperature effects
within the framework of thermal quantum field theory [21, 22].

In Fig. 1, I sketch the general idea of the bound-state effect on the dark
matter freezeout, using electron-positron annihilation as an example. Here, it
is well known, that at low velocities the non-relativistic potential between
the particles will enhance the annihilation rate. Furthermore, the total rate
is dominated by the on-shell formation of positronium, an $e^{+}e^{-}$ bound
state, that consequently decays to photons. The situation is similar in DM
systems. Here, however, early universe conditions have to be taken into
account. At large temperatures, the forming bound states are quickly broken by
thermal collisions, but when temperatures fall, the decay of forming bound
states dominates and provides an additional DM depletion channel. Thus on
general grounds the total DM annihilation rate increases, which also increases
the DM mass at which the correct relic abundance is reproduced.

This effect on the freezeout process can be important in many scenarios. Two
cases should be distinguished, one is the scenario where the mediating
particle is massless, as is the case in color co-annihilation [23, 24]. The
second case becomes relevant, where despite the fact that the potential has a
finite range, bound states can be formed. This is the case when the Bohr
radius of the forming bound state is smaller than the Yukawa radius, set by
the mediator mass, a condition that can be formulated as $M_{V}<\alpha M_{\rm
DM}$, where $M_{\rm DM}$ is the DM mass, $\alpha$ the effective coupling
strength and $M_{V}$ is the mass of the mediating particle. This condition is
satisfied in the case of DM which carries an electroweak ($SU(2)_{L}$) charge.
In particular, EW multiplets that are larger than the triplet are increasingly
strongly affected [25, 7]. A concrete example is the $5$, or quintuplet,
representation of electroweak DM. This representation is particularly
interesting since it is automatically stable without the need to postulate an
additional symmetry [26]. In this case, as we showed in Ref. [5], the expected
DM mass is increased by almost a factor of two and is expected to be around
$14$ TeV. I will discuss interesting phenomenological implications, using this
example candidate in the next section.

Figure 2: The low-energy gamma-ray spectrum of an electroweak $5-$plet DM
candidate, which shows the gamma-ray lines from capture photons on top of the
continuum fragmentation spectrum based on Refs. [25].

###  2.2 New Signals from late-time Bound-State formation

As known from SM observations, the formation of bound states is accompanied by
particle emission to compensate for the released binding energy. Given that DM
is an element of the electroweak multiplet photon emission is possible in the
bound-state formation process, which provides an excellent, new signature for
indirect DM searches. Typical binding energies scale as
$E_{B}\sim\kappa\alpha^{2}M_{\rm DM}$, where $\kappa<1$ is a factor that
depends on the bound state quantum numbers. Thus the typically emitted energy
of the capture-photon will be orders of magnitude lower than the DM mass. This
allows probing very heavy DM candidates by studying photon emission at lower
energies.

In Fig. 2, I show the line spectra of bound-state formation processes in a
model with a Majorana fermion charged under the $SU(2)_{L}$, which lives in
the 5-element multiplet, the quintuplet [5, 25]. The relative line height and
positions allow for a unique reconstruction of DM mass and gauge charge, or in
other words its multiplicity under the electroweak gauge group. This is a
unique situation, where an indirect DM search can provide exact particle
physics properties of the DM candidate. This spectroscopic method will become
particularly relevant for sensitivity improvement of gamma-ray experiments to
heavy dark matter candidates. This method will greatly boost the sensitivity
of H.E.S.S [17], LHASSO [19], as well as the planned CTA [18], to high DM
masses.

##  3 Stable Dark Matter Bound States

Building blocks of SM matter are known to be composite. In particular, the
structure of the strong interaction that leads to composite nucleons has a
conserved global symmetry, which is responsible for the long proton lifetime,
the baryon number. It is plausible that the stability of DM particles
originates from a similar principle, and DM is a stable composite object, see
for example Refs. [27, 28]. Interestingly, the simple assumption that DM is a
composite object opens up the possibility that it is made of constituents
charged under the strong nuclear force [29], an option that has been
previously deemed impossible.

Changing to the simple assumption that DM is a composite particle, can
significantly affect the production in the early universe and the detection
prospects at late times. In particular, since DM which is composite
necessarily features multiple states, including excitations, the understanding
of its production process can provide much more insight into the conditions of
the early universe. Given the high temperatures at which the product is
expected to proceed, this provides a new window into the early universe, at
unprecedentedly early times.

###  3.1 Freezeout in Confining Dark Sectors

As a concrete example of a composite DM model, I will discuss a dark QCD type
model, with the main difference being that the constituent dark quarks have
masses well above the confinement sale. In the confined state, this sector has
a stable particle, which is the lightest baryon, stabilized by global, dark
baryon-number conservation. The dark mesons, carry no baryon number and are
typically unstable. The lightest particles are glueballs, composite gluonic
states, that are also unstable but can feature intermediate lifetimes, due to
a suppressed decay operator. This scenario has been investigated in Refs. [28,
30], and perturbative bound-state formation was included in the computation.
However, the effects of the phase transition (PT), which in this case is a
first-order PT, have not been incorporated. Recently, we demonstrated that
indeed the PT effect changes the expected relic abundance by orders of
magnitude.

Figure 3: Schematic representation of the first order PT and its effect on the
dark baryon abundance. Top row: bubble nucleation, bubble growth, percolation,
and the formation of pockets with trapped free quarks. Bottom row: shrinking
of the pockets, recoupling of the bound-state formation interactions and
baryon leakage from the pocket, confinement with a reduced population of
color-neutral dark baryons.

In Fig. 3, I show how a first-order phase transition can significantly affect
the number of relic particles in the universe, by compressing them in
shrinking false vacuum pockets. In a recent study [31, 32], my collaborators
and I demonstrated that in a large class of dark matter models with a
confining phase transition this effect indeed strongly suppresses the dark-
sector relic abundance. In light of our findings, many concrete studies of
confining dark sectors need to be re-evaluated.

Furthermore, if light degrees of freedom with a relativity long life-time are
present in the theory their decay can lead to further dilution of the DM
states [33], as discussed in Ref. [34] such degrees of freedom, the glueballs,
are unavoidable in this type of theories, and in large fractions of the
parameter space have lifetimes up to a second. This leads to an additional DM
dilution and moves the expected mass scale to the multi PeV scale. This mass
is orders of magnitude above the value expected for a thermal DM candidate
that annihilates with a cross section that respects the unitarity bound [35,
25].

###  3.2 New Experimental Signatures

In previous work, we have proposed that metastable, excited DM bound states
can lead to phenomena similar to radioactive decay in the dark sector [28].
The important conclusion of this finding would be, that ultra-heavy thermal
dark matter can be searched for in indirect detection experiments. This was
previously deemed hopeless, as the low number density suppresses the
annihilation signal. A signal that is due to decay of the dark state, which is
metastable, scales linearly with the number density is thus much stronger, and
allows probing even very heavy DM scenarios.

Another possibility to search for dark sector relics with large masses arises
if the particles have large elastic cross sections, of the order of a barn,
with standard model particles. Such strong interactions might have been missed
by direct detection and other searches. In particular, the low background
detectors, such as XENON1T would be blind to such strong interactions [36].
Since the dark sector particles would scatter multiple times when passing
through the detector, the signals would be rejected as background.
Implementing coincidence measurement techniques can open up the window on the
strong interacting relics with little additional technical effort. In
particular as liquid noble gas, and neutrino detectors gain in volume one
would expect coincident signals to play a role, and dedicated analyses are
being developed [37].

Intriguingly, the compression mechanism, discussed above, opens up the
possibility of forming compact macroscopic objects, see for example [38, 39].
Despite their low number density, their existence can lead to striking
signatures, such as white dwarf star ignition citeGraham:2018efk. Recent
observations of a new class of supernovae events, the so called Calcium-rich
gap transients, could result from such white dwarf-DM collisions [40]. This
would explain their low progenitor masses, which seems to be well below the
Chandrasekhar limit.

Composite DM models can provide, on the one hand, large SM-DM cross sections,
and on the other, feature final states of DM annihilation with lifetimes up to
a second. This opens up two new avenues to search for DM using extrasolar
planets and brown dwarfs. One search is based on energy deposition and heating
of the objects [41]. The other possibility is to search for annihilation
products of the secondary decays of dark sector particles with intermediate
lifetimes [42]. Both searches are significantly affected by DM distribution in
the celestial objects, and in particular populations of captured DM close to
the surface can lead to striking signatures [43]. The research area of
extrasolar planets is rapidly growing. The new infrared telescope JWST has
already started data collection [44], and with other infrared observatories
that will be launched in the near future, we expect excellent data in the
coming few years.

##  4 Conclusion

I have discussed how a simple change in our assumptions about dark matter
substrucutre, can have profound implications on the properties of the dark
sector. In particular, leaving behind the single particle paradigm can lead to
a very different freezeout phenomenology, and thermally produced dark matter
states with masses well above the unitarity bound. Finally, the dark sector
bound-state phenomena allow for investigating those scenarios in a variety of
experimental searches, many of which have a great short-term optimization
potential.

## Acknowledgements

I thank the organizers of the iDM 2022 conference and in particular Dr.
Valentyna Mokina, for organizing this community effort for preserving the
discussions of iDM 2022 in form of this proceeding series.

#### Funding information

I acknowledge funding by the ERC under grant Number 742104.

## References

  * [1] J. Aalbers _et al._ ,  _First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment_ (2022),  2207.03764.
  * [2] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo and F. S. Queiroz,  _The waning of the WIMP? A review of models, searches, and constraints_ ,  Eur. Phys. J. C 78(3), 203 (2018),  10.1140/epjc/s10052-018-5662-y,  1703.07364.
  * [3] M. Ackermann _et al._ ,  _Searching for Dark Matter Annihilation from Milky Way Dwarf Spheroidal Galaxies with Six Years of Fermi Large Area Telescope Data_ ,  Phys. Rev. Lett. 115(23), 231301 (2015),  10.1103/PhysRevLett.115.231301,  1503.02641.
  * [4] R. K. Leane, T. R. Slatyer, J. F. Beacom and K. C. Y. Ng,  _GeV-scale thermal WIMPs: Not even slightly ruled out_ ,  Phys. Rev. D 98(2), 023016 (2018),  10.1103/PhysRevD.98.023016,  1805.10305.
  * [5] A. Mitridate, M. Redi, J. Smirnov and A. Strumia,  _Cosmological Implications of Dark Matter Bound States_ ,  JCAP 05, 006 (2017),  10.1088/1475-7516/2017/05/006,  1702.01141.
  * [6] S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo and L. Vittorio,  _The last complex WIMPs standing_ ,  Eur. Phys. J. C 82(11), 992 (2022),  10.1140/epjc/s10052-022-10918-5,  2205.04486.
  * [7] S. Bottaro, D. Buttazzo, M. Costa, R. Franceschini, P. Panci, D. Redigolo and L. Vittorio,  _Closing the window on WIMP Dark Matter_ ,  Eur. Phys. J. C 82(1), 31 (2022),  10.1140/epjc/s10052-021-09917-9,  2107.09688.
  * [8] B. von Harling and K. Petraki,  _Bound-state formation for thermal relic dark matter and unitarity_ ,  JCAP 12, 033 (2014),  10.1088/1475-7516/2014/12/033,  1407.7874.
  * [9] M. Cirelli, P. Panci, K. Petraki, F. Sala and M. Taoso,  _Dark Matter’s secret liaisons: phenomenology of a dark U(1) sector with bound states_ ,  JCAP 05, 036 (2017),  10.1088/1475-7516/2017/05/036,  1612.07295.
  * [10] F. Kahlhoefer and E. Urdshals,  _On dark atoms, massive dark photons and millicharged sub-components_ ,  Phys. Lett. B 807, 135601 (2020),  10.1016/j.physletb.2020.135601,  2001.04492.
  * [11] E. Braaten, E. Johnson and H. Zhang,  _Zero-range effective field theory for resonant wino dark matter. Part I. Framework_ ,  JHEP 11, 108 (2017),  10.1007/JHEP11(2017)108,  1706.02253.
  * [12] P. Asadi, M. Baumgart, P. J. Fitzpatrick, E. Krupczak and T. R. Slatyer,  _Capture and Decay of Electroweak WIMPonium_ ,  JCAP 02, 005 (2017),  10.1088/1475-7516/2017/02/005,  1610.07617.
  * [13] T. Binder, K. Mukaida and K. Petraki,  _Rapid bound-state formation of Dark Matter in the Early Universe_ ,  Phys. Rev. Lett. 124(16), 161102 (2020),  10.1103/PhysRevLett.124.161102,  1910.11288.
  * [14] J. Harz and K. Petraki,  _Higgs-mediated bound states in dark-matter models_ ,  JHEP 04, 130 (2019),  10.1007/JHEP04(2019)130,  1901.10030.
  * [15] J. Harz and K. Petraki,  _Radiative bound-state formation in unbroken perturbative non-Abelian theories and implications for dark matter_ ,  JHEP 07, 096 (2018),  10.1007/JHEP07(2018)096,  1805.01200.
  * [16] I. Baldes, F. Calore, K. Petraki, V. Poireau and N. L. Rodd,  _Indirect searches for dark matter bound state formation and level transitions_ ,  SciPost Phys. 9(5), 068 (2020),  10.21468/SciPostPhys.9.5.068,  2007.13787.
  * [17] H. Abdalla _et al._ ,  _Searches for gamma-ray lines and ’pure WIMP’ spectra from Dark Matter annihilations in dwarf galaxies with H.E.S.S_ ,  JCAP 11, 037 (2018),  10.1088/1475-7516/2018/11/037,  1810.00995.
  * [18] A. Morselli,  _The Dark Matter Programme of the Cherenkov Telescope Array_ ,  PoS ICRC2017, 921 (2018),  10.22323/1.301.0921,  1709.01483.
  * [19] G. Di Sciascio,  _The LHAASO experiment: from Gamma-Ray Astronomy to Cosmic Rays_ ,  Nucl. Part. Phys. Proc. 279-281, 166 (2016),  10.1016/j.nuclphysbps.2016.10.024,  1602.07600.
  * [20] A. Albert _et al._ ,  _Combined search for neutrinos from dark matter self-annihilation in the Galactic Center with ANTARES and IceCube_ ,  Phys. Rev. D 102(8), 082002 (2020),  10.1103/PhysRevD.102.082002,  2003.06614.
  * [21] T. Binder, A. Filimonova, K. Petraki and G. White,  _Saha equilibrium for metastable bound states and dark matter freeze-out_ ,  Phys. Lett. B 833, 137323 (2022),  10.1016/j.physletb.2022.137323,  2112.00042.
  * [22] T. Binder, B. Blobel, J. Harz and K. Mukaida,  _Dark matter bound-state formation at higher order: a non-equilibrium quantum field theory approach_ ,  JHEP 09, 086 (2020),  10.1007/JHEP09(2020)086,  2002.07145.
  * [23] J. Ellis, J. L. Evans, F. Luo and K. A. Olive,  _Scenarios for Gluino Coannihilation_ ,  JHEP 02, 071 (2016),  10.1007/JHEP02(2016)071,  1510.03498.
  * [24] C. Gross, A. Mitridate, M. Redi, J. Smirnov and A. Strumia,  _Cosmological Abundance of Colored Relics_ ,  Phys. Rev. D 99(1), 016024 (2019),  10.1103/PhysRevD.99.016024,  1811.08418.
  * [25] J. Smirnov and J. F. Beacom,  _TeV-Scale Thermal WIMPs: Unitarity and its Consequences_ ,  Phys. Rev. D 100(4), 043029 (2019),  10.1103/PhysRevD.100.043029,  1904.11503.
  * [26] M. Cirelli, N. Fornengo and A. Strumia,  _Minimal dark matter_ ,  Nucl. Phys. B 753, 178 (2006),  10.1016/j.nuclphysb.2006.07.012,  hep-ph/0512090.
  * [27] O. Antipin, M. Redi, A. Strumia and E. Vigiani,  _Accidental Composite Dark Matter_ ,  JHEP 07, 039 (2015),  10.1007/JHEP07(2015)039,  1503.08749.
  * [28] A. Mitridate, M. Redi, J. Smirnov and A. Strumia,  _Dark Matter as a weakly coupled Dark Baryon_ ,  JHEP 10, 210 (2017),  10.1007/JHEP10(2017)210,  1707.05380.
  * [29] V. De Luca, A. Mitridate, M. Redi, J. Smirnov and A. Strumia,  _Colored Dark Matter_ ,  Phys. Rev. D 97(11), 115024 (2018),  10.1103/PhysRevD.97.115024,  1801.01135.
  * [30] N. A. Dondi, F. Sannino and J. Smirnov,  _Thermal history of composite dark matter_ ,  Phys. Rev. D 101(10), 103010 (2020),  10.1103/PhysRevD.101.103010,  1905.08810.
  * [31] P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer and J. Smirnov,  _Thermal squeezeout of dark matter_ ,  Phys. Rev. D 104(9), 095013 (2021),  10.1103/PhysRevD.104.095013,  2103.09827.
  * [32] P. Asadi, E. D. Kramer, E. Kuflik, G. W. Ridgway, T. R. Slatyer and J. Smirnov,  _Accidentally Asymmetric Dark Matter_ ,  Phys. Rev. Lett. 127(21), 211101 (2021),  10.1103/PhysRevLett.127.211101,  2103.09822.
  * [33] P. Asadi, T. R. Slatyer and J. Smirnov,  _WIMPs without weakness: Generalized mass window with entropy injection_ ,  Phys. Rev. D 106(1), 015012 (2022),  10.1103/PhysRevD.106.015012,  2111.11444.
  * [34] P. Asadi, E. D. Kramer, E. Kuflik, T. R. Slatyer and J. Smirnov,  _Glueballs in a thermal squeezeout model_ ,  JHEP 07, 006 (2022),  10.1007/JHEP07(2022)006,  2203.15813.
  * [35] K. Griest and M. Kamionkowski,  _Unitarity Limits on the Mass and Radius of Dark Matter Particles_ ,  Phys. Rev. Lett. 64, 615 (1990),  10.1103/PhysRevLett.64.615.
  * [36] C. V. Cappiello, J. I. Collar and J. F. Beacom,  _New experimental constraints in a new landscape for composite dark matter_ ,  Phys. Rev. D 103(2), 023019 (2021),  10.1103/PhysRevD.103.023019,  2008.10646.
  * [37] J. Aalbers _et al._ ,  _A Next-Generation Liquid Xenon Observatory for Dark Matter and Neutrino Physics_ (2022),  2203.02309.
  * [38] S. Bottaro, M. Costa and O. Popov,  _Asymmetric accidental composite dark matter_ ,  JHEP 11, 055 (2021),  10.1007/JHEP11(2021)055,  2104.14244.
  * [39] K. Kawana and K.-P. Xie,  _Primordial black holes from a cosmic phase transition: The collapse of Fermi-balls_ ,  Phys. Lett. B 824, 136791 (2022),  10.1016/j.physletb.2021.136791,  2106.00111.
  * [40] J. Smirnov, A. Goobar, T. Linden and E. Mörtsell,  _White Dwarfs in Dwarf Spheroidal Galaxies: A New Class of Compact-Dark-Matter Detectors_ (2022),  2211.00013.
  * [41] R. K. Leane and J. Smirnov,  _Exoplanets as Sub-GeV Dark Matter Detectors_ ,  Phys. Rev. Lett. 126(16), 161101 (2021),  10.1103/PhysRevLett.126.161101,  2010.00015.
  * [42] R. K. Leane, T. Linden, P. Mukhopadhyay and N. Toro,  _Celestial-Body Focused Dark Matter Annihilation Throughout the Galaxy_ ,  Phys. Rev. D 103(7), 075030 (2021),  10.1103/PhysRevD.103.075030,  2101.12213.
  * [43] R. K. Leane and J. Smirnov,  _Floating Dark Matter in Celestial Bodies_ (2022),  2209.09834.
  * [44] B. E. Miles, B. A. Biller, P. Patapis, K. Worthen, E. Rickman, K. K. W. Hoch, A. Skemer, M. D. Perrin, C. H. Chen, S. Mukherjee, C. V. Morley, S. E. Moran _et al._ ,  _The JWST Early Release Science Program for Direct Observations of Exoplanetary Systems II: A 1 to 20 Micron Spectrum of the Planetary-Mass Companion VHS 1256-1257 b_ ,  arXiv e-prints arXiv:2209.00620 (2022),  2209.00620.


# Biharmonic distance of graphs

Yulong Wei1 Rong-hua Li2***Corresponding authors.
E-mail address<EMAIL_ADDRESS>(Y. Wei<EMAIL_ADDRESS>(R-H.
Li), yangweihua@tyut.edu.cn(W. Yang). Weihua Yang1

1Department of Mathematics, Taiyuan University of Technology, Taiyuan, 030024,
China
2School of Computer Science & Technology, Beijing Institute of Technology,
Beijing, 100081, China

Abstract Lipman et al. [ACM Transactions on Graphics 29 (3) (2010), 1–11]
introduced the concept of biharmonic distance to measure the distances between
pairs of points on a 3D surface. Biharmonic distance has some advantages over
resistance distance and geodesic distance in some realistic contexts.
Nevertheless, limited work has been done on the biharmonic distance in the
discrete case. In this paper, we give some characterizations of the biharmonic
distance of a graph. Some basic mathematical properties of biharmonic distance
and biharmonic index are established.

Keywords Biharmonic distance; Laplacian matrix; Pseudoinverse Laplacian;
Kirchhoff index; Hypercube

##  1 Introduction

Measuring the distances between pairs of vertices of a graph is a classical
problem and thus several distance measures have previously been proposed and
are commonly used in many fields. For example, the geodesic distance is using
the length of the shortest path to measure distance which is intuitive and
useful. The geodesic distance of any two vertices of a graph can be computed
by the breadth-first search algorithm. The diameter of a graph is the maximal
geodesic distance between vertices. Network topology is always represented by
a graph, where vertices represent processors and edges represent links between
processors. The diameter is an important parameter for interconnection
networks. As another example, resistance distance, which originated in
electrical circuit theory, has played a prominent role in circuit theory [11],
chemistry [15, 18], combinatorial matrix theory [3, 22] and spectral graph
theory [6]. Moreover, resistance distance has extensive applications ranging
from quantifying biological structures [18], distributed control systems [4]
and power grid systems [20].

The concept of biharmonic distance was first proposed in [16] as a measure of
distance between two points on a curve surface. It has some advantages over
geodesic distance and resistance distance as a metric that incorporates both
local and global graph structure in some realistic contexts, e.g., computer
graphics, geometric processing, and shape analysis. Moreover, it can be used
to measure the robustness of the second order noisy consensus problem without
leaders [2]. In 2018, Yi et al. [23] defined the biharmonic distance index
(i.e., biharmonic Kirchhoff index [24], or biharmonic index [21]) of a graph
to further describe the behavior of second-order consensus dynamics and
established a connection between the biharmonic distance of a graph and its
second-order network coherence. Since direct computation of biharmonic
distances is computationally infeasible for huge networks with millions of
vertices, Zhang et al. [25] developed a nearly linear-time algorithm to
approximate all diagonal entries of pseudoinverse of the square of graph
Laplacian matrix by the Johnson-Lindenstrauss lemma and Laplacian solvers.
However, limited work has been done on the biharmonic distance of a graph.

In this paper, we investigate some properties on the biharmonic distance and
biharmonic index of an undirected connected simple graph. First, we give some
characterizations of the biharmonic distance of a graph and obtain some bounds
for it. Second, we reveal a relationship between biharmonic index and
Kirchhoff index and determine the unique graph having the minimum biharmonic
index among the connected graphs with given order. As applications, we study
the biharmonic distance of some special graphs, such as the complement of a
graph, the Cartesian product of two graphs and the Cayley graph of a finite
abelian group. In particular, the combinatorial expression of biharmonic
distance between any two vertices of hypercubes are determined.

The rest of this paper is organized as follows. Section 2 introduces some
terminology and preliminaries. Some characterizations of the biharmonic
distance of a graph are given in Section 3. The biharmonic index of a graph is
investigated in Section 4. Some applications are given in Section 5. Section 6
concludes the paper and proposes some open problems.

##  2 Terminology and preliminaries

Let $G$ be a connected graph with $V(G)=\\{v_{1},\ldots,v_{n}\\}$ and
$E(G)=\\{e_{1},\ldots,e_{m}\\}$ and let $G^{c}$ be the complement of a graph
$G$. The neighborhood $N_{G}(v)$ of a vertex $v$ in $G$ is the set of vertices
adjacent to $v$. The degree of a vertex $v$ in $G$ is $|N_{G}(v)|$. A graph
$G$ is $k$-regular if $|N_{G}(v)|=k$ for any $v\in V(G)$. We refer readers to
[5] for terminology and notation unless stated otherwise.

The adjacency matrix of the graph is denoted by $A$ and $D$ represents the
diagonal matrix of vertex degrees. The Laplacian matrix of $G$ is the matrix
$D-A$, denoted by $L$ (or $L(G)$). Throughout, $O$ denotes the all-$0$ matrix,
$\mathbf{0}$ the all-$0$ vector and $\mathbf{1}$ the all-$1$ vector. Suppose
that $\delta_{u}$ is the elementary unit vector with a $1$ in coordinate $u$.
In the whole paper, we write $\delta_{uv}$ instead of $\delta_{u}-\delta_{v}$.

The following result shows that the spectrum of $L$ determines the number of
connected components in a graph.

######  Lemma 2.1 ([8])

The multiplicity of $0$ as an eigenvalue of $L$ is equal to the number of
connected components in $G$.

Thus, if $G$ is a connected graph, then the multiplicity of $0$ as an
eigenvalue of $L$ is $1$. In fact, the matrix $L$ is symmetric and positive
semidefinite, which means $L$ has a spectral decomposition as

$L=\sum_{k=1}^{n}\lambda_{k}z_{k}z_{k}^{\top},$

where $0=\lambda_{1}<\lambda_{2}\leq\cdots\leq\lambda_{n}$ are its $n$
positive eigenvalues, and $z_{1},\ldots,z_{n}$ are the corresponding mutually
orthogonal unit eigenvectors.

The following two lemmas are important to obtain our main results.

######  Lemma 2.2 ([9])

Let $G$ be a graph with $n$ vertices and $E(G)\neq\emptyset$. Then
$\lambda_{2}(G)=\cdots=\lambda_{n}(G)$ if and only if $G$ is isomorphic to
$K_{n}$.

######  Lemma 2.3 (Matrix-Tree Theorem [3])

Let $G$ be a graph. Then the cofactor of any element of $L$ equals the number
of spanning trees of $G$.

Let $B$ be an $m\times n$ matrix. A matrix $H$ of order $n\times m$ is said to
be a _$g$ -inverse matrix_ of $B$ if $BHB=B$. The _Moore Penrose inverse_ of
an $n\times m$ real matrix $M$, denoted by $M^{+}$, is the $m\times n$ real
matrix that satisfies the following equations:

$MM^{+}M=M,M^{+}MM^{+}=M^{+},(MM^{+})^{\top}=MM^{+},(M^{+}M)^{\top}=M^{+}M.$

We call the Moore Penrose inverse of $L$ pseudoinverse Laplacian and write
$L^{2+}$ instead of $(L^{+})^{2}$. Similarly, the matrix $L^{+}$ also has a
spectral decomposition as

$L^{+}=\sum_{k=2}^{n}\lambda_{k}^{-1}z_{k}z_{k}^{\top},$ (1)

where $\lambda_{k}$ is the $k$th smallest eigenvalue of the Laplacian matrix
$L$ of a connected graph $G$ for $k=2,\ldots,n$ and $z_{2},\ldots,z_{n}$ are
the corresponding mutually orthogonal unit eigenvectors. Actually, the
eigenvalues of $L^{+}$ are
$0,\lambda_{n}^{-1},\lambda_{n-1}^{-1},\ldots,\lambda_{2}^{-1}$ [3].

The notion of biharmonic distance was first proposed in [16] as a measure of
distance between two points $u$ and $v$ on a curve surface:

$d_{B}^{2}(u,v)=g_{d}(u,u)+g_{d}(v,v)-g_{d}(u,v),$

where $g_{d}$ is the discrete Green’s function of the discretized, normalized
bilaplacian $\widetilde{L}^{2}$, equivalent to the pseudoinverse of
$\widetilde{L}^{2}$, and $\widetilde{L}$ is the normalization of Laplacian
$L$. Fitch and Leonard [13] defined the biharmonic distance between two
vertices $i$ and $j$ in the graph $G$, which they denoted $\gamma_{i,j}$,
analogously without normalizing $L$:

$\gamma_{i,j}=L_{ii}^{2+}+L_{jj}^{2+}-2L_{ij}^{2+}.$

In this paper we adopt the definition in [23], which is a fusion of the above
two conceptions, defined as follows.

######  Definition 2.4 ([16, 13, 23])

For a graph $G$, the biharmonic distance $d_{B}(u,v)$ between two vertices $u$
and $v$ is defined by

$d_{B}^{2}(u,v)=L_{uu}^{2+}+L_{vv}^{2+}-2L_{uv}^{2+},$ (2)

where $L^{2+}_{ij}$ means the $(i,j)$-entry of $L^{2+}$.

Biharmonic distance is a metric, which satisfies the following properties:
non-negativity, nullity, symmetry, and triangle inequality [16]. The
biharmonic index [23], as a distance-based graph invariant, is used for
describing the behavior of second-order consensus dynamics. The concept of
biharmonic index of graphs is defined as follows.

######  Definition 2.5 ([23, 21])

For a connected graph $G$, the biharmonic index of $G$ is equal to

$\frac{1}{2}\sum_{u\in V(G)}\sum_{v\in V(G)}d_{B}^{2}(u,v),$

denoted by $B(G)$.

##  3 The biharmonic distance of a graph

In this section, we will give some characterizations of the biharmonic
distance of a graph. Given a connected graph $G$, suppose that $\lambda_{k}$
is the $k$th smallest eigenvalue of the Laplacian matrix of $G$ for
$k=1,2,\ldots,n$ and $z_{1},z_{2},\ldots,z_{n}$ are the corresponding mutually
orthogonal real unit eigenvectors. Since $L\mathbf{1}=\mathbf{0}$, set
$z_{1}=\frac{1}{\sqrt{n}}\mathbf{1}$. Let $B$ be an $m\times n$ matrix. If
$S\subset\\{1,\ldots,m\\}$, $T\subset\\{1,\ldots,n\\}$, then $B[S|T]$ will
denote the submatrix of $B$ determined by the rows corresponding to $S$ and
the columns corresponding to $T$. The submatrix obtained by deleting the rows
indexed in $S$ and the columns indexed in $T$ will be denoted by $B(S|T)$.
When $S=\\{i\\}$, $T=\\{j\\}$ are singletons, then $B(S|T)$ is denoted by
$B(i|j)$. The _rank_ of a matrix $B$ is denoted by $\textrm{r}(B)$.

###  3.1 The spectral characterization

###### Theorem 3.1

Let $G$ be a connected graph with $n$ vertices. Then the biharmonic distance
between two vertices $u$ and $v$ is

$d_{B}(u,v)=(\sum_{k=2}^{n}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2})^{\frac{1}{2}}.$
(3)

Proof. Note that $z_{1},z_{2},\ldots,z_{n}$ is a normal orthogonal basic of
$\mathbf{R}^{n}$. Then by formulas (1) and (2),

$\displaystyle d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\delta_{uv}^{\top}L^{2+}\delta_{uv}$ $\displaystyle=$
$\displaystyle\delta_{uv}^{\top}(L^{+}L^{+})\delta_{uv}$ $\displaystyle=$
$\displaystyle\delta_{uv}^{\top}(\sum_{k=2}^{n}\lambda_{k}^{-1}z_{k}z_{k}^{\top})(\sum_{k=2}^{n}\lambda_{k}^{-1}z_{k}z_{k}^{\top})\delta_{uv}$
$\displaystyle=$
$\displaystyle\delta_{uv}^{\top}(\sum_{k=2}^{n}\lambda_{k}^{-2}z_{k}z_{k}^{\top})\delta_{uv}$
$\displaystyle=$
$\displaystyle\sum_{k=2}^{n}\lambda_{k}^{-2}(\delta_{uv}^{\top}z_{k}z_{k}^{\top}\delta_{uv})$
$\displaystyle=$
$\displaystyle\sum_{k=2}^{n}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}.$

Thus, we obtain the desired result. $\hfill\Box\vskip 6.0pt plus 2.0pt minus
2.0pt$

Let $S=\\{k\mid 2\leq k\leq n\\}$, $\sigma_{2}=\\{k\mid 3\leq k\leq
n~{}\textmd{and}~{}\lambda_{k}>\lambda_{2}\\}$ and $\sigma_{n}=\\{k\mid 2\leq
k\leq n-1~{}\textmd{and}~{}\lambda_{k}<\lambda_{n}\\}$. The eigenspace of $L$
corresponding to $\lambda_{k}$ is denoted by $V_{\lambda_{k}}$ and the
orthogonal complement of $V_{\lambda_{k}}$ is denoted by
$V_{\lambda_{k}}^{\bot}$ for $1\leq k\leq n$.

###### Theorem 3.2

Let $G$ be a connected graph with $n$ vertices. Then the biharmonic distance
between two distinct vertices $u$ and $v$ satisfies that

$\sqrt{2}\lambda_{n}^{-1}\leq d_{B}(u,v)\leq\sqrt{2}\lambda_{2}^{-1},$

and $d_{B}(u,v)=\sqrt{2}\lambda_{j}^{-1}$ if and only if $G\cong K_{n}$ or
$\delta_{uv}\in V_{\lambda_{k}}^{\bot}$ for any $k\in\sigma_{j}$, where
$j\in\\{2,n\\}$.

Proof. Since

$\displaystyle\sum_{k=2}^{n}(z_{k}(u)-z_{k}(v))^{2}$ $\displaystyle=$
$\displaystyle\sum_{k=2}^{n}z_{k}^{2}(u)+\sum_{k=2}^{n}z_{k}^{2}(v)-2\sum_{k=2}^{n}z_{k}(u)z_{k}(v)$
$\displaystyle=$
$\displaystyle(\sum_{k=1}^{n}z_{k}^{2}(u)-z_{1}^{2}(u))+(\sum_{k=1}^{n}z_{k}^{2}(v)-z_{1}^{2}(v))-2(\sum_{k=1}^{n}z_{k}(u)z_{k}(v)-z_{1}(u)z_{1}(v))$
$\displaystyle=$
$\displaystyle(1-\frac{1}{n})+(1-\frac{1}{n})-2(0-\frac{1}{\sqrt{n}}\times\frac{1}{\sqrt{n}})$
$\displaystyle=$ $\displaystyle 2,$

by formula (3), we have

$d_{B}(u,v)\leq\lambda_{2}^{-1}(\sum_{k=2}^{n}(z_{k}(u)-z_{k}(v))^{2})^{\frac{1}{2}}=\sqrt{2}\lambda_{2}^{-1}.$

Note that the Laplacian eigenvalues of $K_{n}$ are $n$ with multiplicities
$n-1$ and $0$ with multiplicity $1$. Then
$d_{B}(u,v)=(\sum_{k=2}^{n}n^{-2}(z_{k}(u)-z_{k}(v))^{2})^{\frac{1}{2}}=\dfrac{\sqrt{2}}{n}$
for any $u$, $v\in V(K_{n})$. Therefore, if $G\cong K_{n}$, then the equality
holds.

If $\delta_{uv}\in V_{\lambda_{k}}^{\bot}$ for any $k\in\sigma_{2}$, then
$z_{k}(u)-z_{k}(v)=\delta_{uv}^{\top}z_{k}=0$ for any $z_{k}\in
V_{\lambda_{k}}$ and $k\in\sigma_{2}$. Thus,

$\displaystyle d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\sum_{k=2}^{n}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle=$
$\displaystyle\sum_{k\in\sigma_{2}}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}+\sum_{k\in
S\setminus\sigma_{2}}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}$ $\displaystyle=$
$\displaystyle\sum_{k\in\sigma_{2}}(\lambda_{k}^{-2}\times 0)+\sum_{k\in
S\setminus\sigma_{2}}\lambda_{2}^{-2}(z_{k}(u)-z_{k}(v))^{2}$ $\displaystyle=$
$\displaystyle\lambda_{2}^{-2}\sum_{k=2}^{n}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle=$ $\displaystyle 2\lambda_{2}^{-2}.$

Thus, the equality holds.

On the other hand, if $d_{B}(u,v)=\sqrt{2}\lambda_{2}^{-1}$ and $G$ is not
isomorphic to $K_{n}$, then by Lemma 2.2, $\sigma_{2}\neq\emptyset$. Assume to
the contrary that there exists an integer $k_{1}\in\sigma_{2}$ such that
$\delta_{uv}\notin V_{\lambda_{k_{1}}}^{\bot}$. Then
$z_{k_{1}}(u)-z_{k_{1}}(v)=\delta_{uv}^{\top}z_{k_{1}}\neq 0$. Thus,

$\displaystyle d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\sum_{k=2}^{n}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle=$
$\displaystyle\lambda_{k_{1}}^{-2}(z_{k_{1}}(u)-z_{k_{1}}(v))^{2}+\sum_{k\in
S\setminus\\{k_{1}\\}}\lambda_{k}^{-2}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle<$
$\displaystyle\lambda_{2}^{-2}(z_{k_{1}}(u)-z_{k_{1}}(v))^{2}+\sum_{k\in
S\setminus\\{k_{1}\\}}\lambda_{2}^{-2}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle=$
$\displaystyle\lambda_{2}^{-2}\sum_{k=2}^{n}(z_{k}(u)-z_{k}(v))^{2}$
$\displaystyle=$ $\displaystyle 2\lambda_{2}^{-2},$

a contradiction.

By the similar arguments, we have $d_{B}(u,v)\geq\sqrt{2}\lambda_{n}^{-1}$ and
equality holds if and only if $G\cong K_{n}$ or $\delta_{uv}\in
V_{\lambda_{k}}^{\bot}$ for any $k\in\sigma_{n}$.

As mentioned above, we complete this proof. $\hfill\Box\vskip 6.0pt plus 2.0pt
minus 2.0pt$

Figure 1: The structures of $K_{4}^{-}$ and $W_{5}$.

###### Example 3.1

Let $K_{4}^{-}$ be the graph obtained from $K_{4}$ by deleting an edge and let
$W_{5}$ be a wheel graph with $5$ vertices (see Figure 1). The eigenvalues of
the Laplacian matrix of $K_{4}^{-}$ are $0$, $2$, $4$, $4$. Then
$\sigma_{4}=\\{2\\}$. Suppose
$x=(x_{u_{1}},\ldots,x_{u_{4}})^{\top}\in\mathbf{R}^{4}$. Solve the system of
linear equations $(2I-L(K_{4}^{-}))x=\mathbf{0}$ and obtain that
$x_{u_{1}}=x_{u_{3}}=0$. Then $d_{B}(u_{1},u_{3})=\dfrac{\sqrt{2}}{4}$ in
$K_{4}^{-}$. The eigenvalues of the Laplacian matrix of $W_{5}$ are $0$, $3$,
$3$, $5$, $5$. Then $\sigma_{2}=\\{4,5\\}$. Suppose
$y=(y_{v_{1}},\ldots,y_{v_{5}})^{\top}\in\mathbf{R}^{5}$. Solve the system of
linear equations $(5I-L(W_{5}))x=\mathbf{0}$ and obtain that
$y_{v_{2}}=y_{v_{4}}$. Then $d_{B}(v_{2},v_{4})=\dfrac{\sqrt{2}}{3}$ in
$W_{5}$.

###  3.2 The determinant characterization

Let $L$ be the Laplacian matrix of a connected graph $G$ with $n$ vertices. In
this section, we first give a formula for biharmonic distance of $G$ in terms
of the $g$-inverse matrix of $L^{2}$ and then by this formula, an expression
for biharmonic distance between any two distinct vertices of $G$ is
established in terms of an $(n-2)\times(n-2)$ principal minor of $L^{2}$ and
the number of spanning trees of $G$.

###### Lemma 3.3

Let $G$ be a connected graph with $n$ vertices. If $u,v\in V(G)$, then

$d_{B}^{2}(u,v)=\delta_{uv}^{\top}H\delta_{uv},$

where $L$ is the Laplacian matrix of $G$ and $H$ is a $g$-inverse matrix of
$L^{2}$.

Proof. $\mathbf{Claim~{}1.}$ $L^{2+}=(L^{2})^{+}$.

By formula (1),
$L^{2+}=L^{+}L^{+}=(\sum_{k=2}^{n}\lambda_{k}^{-1}z_{k}z_{k}^{\top})^{2}=\sum_{i=2}^{n}\lambda_{k}^{-2}z_{k}z_{k}^{\top}=(L^{2})^{+}$.

$\mathbf{Claim~{}2.}$
$\delta_{uv}^{\top}H_{1}\delta_{uv}=\delta_{uv}^{\top}H_{2}\delta_{uv}$ for
any two $g$-inverse matrices $H_{1}$ and $H_{2}$ of $L^{2}$.

We use Col$(B)$ to denote the column space of a matrix $B$. Since $L$ is a
symmetric matrix, we have
$\textrm{r}(L^{2})=\textrm{r}(L^{\top}L)=\textrm{r}(L)$. Thus,
$\dim$(Col($L^{2}))=\dim$(Col$(L))$. Note that
Col$(L^{2})$$\subseteq$Col$(L)$. Then Col$(L^{2})$$=$Col$(L)$. Since $G$ is a
connected graph, by Lemma 2.1, $\textrm{r}(L)=n-1$. Since
$L\mathbf{1}=\mathbf{0}$ and $\delta_{uv}^{\top}\mathbf{1}=0$,
$\delta_{uv}\in$Col$(L)$. Then $\delta_{uv}$$\in$Col$(L^{2})$. Therefore,
there exists a vector $z$ such that $L^{2}z=\delta_{uv}$. Thus, considering
$H_{1}$, $H_{2}$ any two $g$-inverse matrices of $L^{2}$,

$\displaystyle\delta_{uv}^{\top}H_{1}\delta_{uv}-\delta_{uv}^{\top}H_{2}\delta_{uv}$
$\displaystyle=$ $\displaystyle\delta_{uv}^{\top}(H_{1}-H_{2})\delta_{uv}$
$\displaystyle=$ $\displaystyle(L^{2}z)^{\top}(H_{1}-H_{2})(L^{2}z)$
$\displaystyle=$ $\displaystyle z^{\top}((L^{2})^{\top}(H_{1}-H_{2})L^{2})z$
$\displaystyle=$ $\displaystyle z^{\top}(L^{2}H_{1}L^{2}-L^{2}H_{2}L^{2})z$
$\displaystyle=$ $\displaystyle z^{\top}(L^{2}-L^{2})z$ $\displaystyle=$
$\displaystyle 0.$

By Claims 1 and 2, we have
$d_{B}^{2}(u,v)=\delta_{uv}^{\top}L^{2+}\delta_{uv}=\delta_{uv}^{\top}(L^{2})^{+}\delta_{uv}=\delta_{uv}^{\top}H\delta_{uv}$,
where $H$ is a $g$-inverse matrix of $L^{2}$. $\hfill\Box\vskip 6.0pt plus
2.0pt minus 2.0pt$

###### Theorem 3.4

Let $G$ be a connected graph with $n$ vertices. Then the biharmonic distance
between two distinct vertices $u$ and $v$ is

$d_{B}(u,v)=\dfrac{\sqrt{\det(L^{2}(u,v|u,v))}}{\sqrt{n}\tau(G)},$

where $L$ is the Laplacian matrix of $G$ and $\tau(G)$ is the number of
spanning trees of $G$.

Proof. Let $C=L^{2}(v|v)$. Then by Cauchy-Binet Formula and Lemma 2.3, we have

$\displaystyle\det C$ $\displaystyle=$ $\displaystyle\det(L^{2}(v|v))$
$\displaystyle=$ $\displaystyle\sum_{S\subseteq
V(G),|S|=n-1}(\det(L[V(G)\setminus\\{v\\}|S]))^{2}$ $\displaystyle=$
$\displaystyle n\tau^{2}(G)$ $\displaystyle>$ $\displaystyle 0.$

We construct the following matrix $H$: In $L^{2}$, replace entries in the
$v$th row and column by zeros and replace $L^{2}(v|v)$ by $C^{-1}$.

$\mathbf{Claim.}$ $L^{2}HL^{2}=L^{2}$.

By the construction of $H$, there exists a permutation matrix $P$ such that
$P^{\top}HP=\begin{bmatrix}0&O\\\ O&Q^{\top}C^{-1}Q\end{bmatrix}$ and
$P^{\top}L^{2}P=\begin{bmatrix}b_{11}&b^{\top}\\\ b&Q^{\top}CQ\end{bmatrix}$,
where $Q$ is a real $(n-1)\times(n-1)$ matrix such that $Q^{\top}Q=I$. Then

$\displaystyle P^{\top}L^{2}HL^{2}P$ $\displaystyle=$
$\displaystyle(P^{\top}L^{2}P)(P^{\top}HP)(P^{\top}L^{2}P)$ $\displaystyle=$
$\displaystyle\begin{bmatrix}b_{11}&b^{\top}\\\
b&Q^{\top}CQ\end{bmatrix}\begin{bmatrix}0&O\\\
O&Q^{\top}C^{-1}Q\end{bmatrix}\begin{bmatrix}b_{11}&b^{\top}\\\
b&Q^{\top}CQ\end{bmatrix}$ $\displaystyle=$
$\displaystyle\begin{bmatrix}0&b^{\top}Q^{\top}C^{-1}Q\\\
O&I\end{bmatrix}\begin{bmatrix}b_{11}&b^{\top}\\\ b&Q^{\top}CQ\end{bmatrix}$
$\displaystyle=$
$\displaystyle\begin{bmatrix}b^{\top}Q^{\top}C^{-1}Qb&b^{\top}\\\
b&Q^{\top}CQ\end{bmatrix}.$

Note that $\textrm{r}(L^{2})=\textrm{r}(L)=n-1=\textrm{r}(C)$. Then

$\displaystyle\textrm{r}(Q^{\top}CQ)$ $\displaystyle=$
$\displaystyle\textrm{r}(P^{\top}L^{2}P)$ $\displaystyle=$
$\displaystyle\textrm{r}\left(\begin{bmatrix}b_{11}&b^{\top}\\\
b&Q^{\top}CQ\end{bmatrix}\right)$ $\displaystyle=$
$\displaystyle\textrm{r}\left(\begin{bmatrix}b_{11}-b^{\top}(Q^{\top}CQ)^{-1}b&O\\\
b&Q^{\top}CQ\end{bmatrix}\right)$ $\displaystyle=$
$\displaystyle\textrm{r}\left(\begin{bmatrix}b_{11}-b^{\top}(Q^{\top}CQ)^{-1}b&O\\\
O&Q^{\top}CQ\end{bmatrix}\right)$ $\displaystyle=$
$\displaystyle\textrm{r}(Q^{\top}CQ)+\textrm{r}([b_{11}-b^{\top}(Q^{\top}CQ)^{-1}b]).$

Therefore, $\textrm{r}([b_{11}-b^{\top}(Q^{\top}CQ)^{-1}b])=0$ and
$b_{11}=b^{\top}(Q^{\top}CQ)^{-1}b$. Thus,
$P^{\top}L^{2}HL^{2}P=P^{\top}L^{2}P$ and so $L^{2}HL^{2}=L^{2}$.

According to the above claim, $H$ is a $g$-inverse of $L^{2}$. By Lemma 3.3,

$\displaystyle d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\delta_{uv}^{\top}H\delta_{uv}$ $\displaystyle=$ $\displaystyle
h_{uu}+h_{vv}-h_{uv}-h_{vu}$ $\displaystyle=$ $\displaystyle h_{uu}$
$\displaystyle=$ $\displaystyle(C^{-1})_{uu}$ $\displaystyle=$
$\displaystyle\dfrac{\det(L^{2}(u,v|u,v))}{\det(L^{2}(v|v))}$ $\displaystyle=$
$\displaystyle\dfrac{\det(L^{2}(u,v|u,v))}{n\tau^{2}(G)}.$

This completes the proof of Theorem 3.4. $\hfill\Box\vskip 6.0pt plus 2.0pt
minus 2.0pt$

###  3.3 The system of linear equations characterization

In this section, we will characterize biharmonic distance between any two
distinct vertices of a connected graph from the perspective of the system of
linear equations. Given a real vector $f\in\mathbf{R}^{n}$, let
$\|f\|_{2}=\sqrt{f^{\top}f}$. We call $\|f\|_{2}$ the _norm_ of $f$.

###### Theorem 3.5

Let $G$ be a connected graph with $n$ vertices. If $f$ is a solution of the
equation $Lx=\delta_{uv}$ with minimum norm, then the biharmonic distance
between two vertices $u$ and $v$ is

$d_{B}(u,v)=\|f\|_{2}.$

Proof. If $f$ is a solution of the equation $Lx=\delta_{uv}$ with minimum
norm, then $f=L^{+}\delta_{uv}$. Thus,

$\displaystyle d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\delta_{uv}^{\top}L^{2+}\delta_{uv}$ $\displaystyle=$
$\displaystyle(L^{+}\delta_{uv})^{\top}(L^{+}\delta_{uv})$ $\displaystyle=$
$\displaystyle f^{\top}f$ $\displaystyle=$ $\displaystyle\|f\|_{2}^{2}.$

Therefore, $d_{B}(u,v)=\|f\|_{2}$. $\hfill\Box\vskip 6.0pt plus 2.0pt minus
2.0pt$

##  4 The biharmonic index of a graph

In this section, we investigate the biharmonic index $B(G)$ of a connected
graph $G$. We first derive an expression for the biharmonic index of a graph
$G$ in terms of the Laplacian eigenvalues of $G$. Then we reveal a
relationship between biharmonic index and Kirchhoff index. Moreover, we study
the variation of $B(G)$ when modifications on the graph are made, in
particular when an edge is added. A lower bound for $B(G)$ is given and the
equality case is studied.

###  4.1 The spectral characterization

###### Theorem 4.1

Let $G$ be a connected graph with $n$ vertices and $\lambda_{i}$ the $i$th
smallest eigenvalue of the Laplacian matrix of $G$. Then

$B(G)=n\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}}.$

Proof. Since

$L^{2+}\mathbf{1}=L^{+}L^{+}\mathbf{1}=L^{+}L^{+}LL^{+}\mathbf{1}=L^{+}L^{+}L^{+}L\mathbf{1}=\mathbf{0},$

we have $\sum_{v\in V(G)}L_{uv}^{2+}=0$ for any $u\in V(G)$. Then

$\displaystyle B(G)$ $\displaystyle=$ $\displaystyle\dfrac{1}{2}\sum_{u\in
V(G)}\sum_{v\in V(G)}d_{B}^{2}(u,v)$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\sum_{u\in V(G)}\sum_{v\in
V(G)}(L_{uu}^{2+}+L_{vv}^{2+}-2L_{uv}^{2+})$ $\displaystyle=$ $\displaystyle
n\sum_{u\in V(G)}L_{uu}^{2+}$ $\displaystyle=$ $\displaystyle
n~{}\text{trace}(L^{2+})$ $\displaystyle=$ $\displaystyle
n\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}}.$

This completes the proof of Theorem 4.1. $\hfill\Box\vskip 6.0pt plus 2.0pt
minus 2.0pt$

###  4.2 A relationship to Kirchhoff index

For a graph $G$, its Kirchhoff index $Kf(G)$ [15] is defined as the sum of
resistance distance over all its vertex pairs. That is,
$Kf(G)=\sum_{i<j}r_{ij}$, where $r_{ij}$ is the resistance distance between
the vertices $v_{i}$ and $v_{j}$ in $G$. This index can be expressed via the
non-zero Laplacian eigenvalues as follows.

######  Lemma 4.2 ([14])

If $G$ be a connected graph with $n$ vertices, then

$Kf(G)=n\sum_{i=2}^{n}\frac{1}{\lambda_{i}}.$

Now, we establish a relationship between biharmonic index and Kirchhoff index
in the following discussion.

###### Theorem 4.3

Let $G$ be a connected graph with $n$ vertices. Then for $n\geq 2$,

$B(G)\geq\dfrac{Kf(G)^{2}}{n(n-1)},$

where the equality holds if and only if $G$ is isomorphic to $K_{n}$.

Proof. By Theorem 4.1, Cauchy-Schwarz inequality and Lemma 4.2, we have

$\displaystyle B(G)$ $\displaystyle=$ $\displaystyle
n\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}}$ $\displaystyle=$
$\displaystyle\frac{n}{n-1}\cdot(n-1)\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}}$
$\displaystyle\geq$
$\displaystyle\frac{n}{n-1}(\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}})^{2}$
$\displaystyle=$ $\displaystyle\frac{n}{n-1}(\dfrac{Kf(G)}{n})^{2}$
$\displaystyle=$ $\displaystyle\dfrac{Kf(G)^{2}}{n(n-1)}.$

Note that the equality holds if and only if $\lambda_{2}=\cdots=\lambda_{n}$.
By Lemma 2.2, we obtain the desired results. $\hfill\Box\vskip 6.0pt plus
2.0pt minus 2.0pt$

###  4.3 The biharmonic index of connected graphs with given order

######  Lemma 4.4 ([8])

If $e$ is an edge of the graph $G$ and $G^{\prime}=G-e$, then

$0=\lambda_{1}(G^{\prime})=\lambda_{1}(G)\leq\lambda_{2}(G^{\prime})\leq\lambda_{2}(G)\leq\cdots\leq\lambda_{n-1}(G)\leq\lambda_{n}(G^{\prime})\leq\lambda_{n}(G),$

where $\lambda_{1}(G),\ldots,\lambda_{n}(G)$ and
$\lambda_{1}(G^{\prime}),\ldots,\lambda_{n}(G^{\prime})$ are the Laplacian
eigenvalues of $G$ and $G^{\prime}$, respectively.

###### Theorem 4.5

Let $G$ be a connected graph with $n$ vertices and $e\in E(G^{c})$. Then

$B(G+e)<B(G).$

Proof. Note that $\sum_{i=2}^{n}\lambda_{i}(G)=\textrm{trace}~{}L(G)=2|E(G)|$
and $\sum_{i=2}^{n}\lambda_{i}(G+e)=2(|E(G)|+1)$. Then
$\sum_{i=2}^{n}\lambda_{i}(G+e)-(\sum_{i=2}^{n}\lambda_{i}(G))=2$. By Lemma
4.4, there exists some $j$ such that $\lambda_{j}(G+e)>\lambda_{j}(G)$. By
Theorem 4.1, we have

$\displaystyle B(G+e)$ $\displaystyle=$ $\displaystyle
n\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}(G+e)}$ $\displaystyle<$
$\displaystyle n\sum_{i=2}^{n}\dfrac{1}{\lambda_{i}^{2}(G)}$ $\displaystyle=$
$\displaystyle B(G).$

This completes the proof of Theorem 4.5. $\hfill\Box\vskip 6.0pt plus 2.0pt
minus 2.0pt$

###### Theorem 4.6

Let $G$ be a connected graph with $n$ vertices. Then

$B(G)\geq\dfrac{n-1}{n},$

where the equality holds if and only if $G$ is isomorphic to $K_{n}$.

Proof. By the proof of Theorem 3.2, $d_{B}(u,v)=\dfrac{\sqrt{2}}{n}$ for any
two distinct vertices $u$, $v\in V(K_{n})$. Therefore, by Theorem 4.5, we have

$B(G)\geq
B(K_{n})=(\dfrac{\sqrt{2}}{n})^{2}\times\dbinom{n}{2}=\dfrac{n-1}{n}.$

Suppose $B(G)=\dfrac{n-1}{n}$. Assume to the contrary that $G$ is not a
complete graph and $e\in E(G^{c})$. Then by Theorem 4.5, we have
$\dfrac{n-1}{n}\leq B(G+e)<B(G)=\dfrac{n-1}{n}$, a contradiction. Thus,
$B(G)=\dfrac{n-1}{n}$ if and only if $G$ is isomorphic to $K_{n}$.
$\hfill\Box\vskip 6.0pt plus 2.0pt minus 2.0pt$

###### Remark 4.1

Various bounds on the sum of powers of the Laplacian eigenvalues of graphs
have been established in terms of the number of vertices, number of edges,
maximum degree, clique number, independence number, matching number or the
number of spanning trees of graphs. We refer readers to [7, 10, 17, 26]. By
Theorem 4.1, other lower bounds on the biharmonic index of connected graphs
can be obtained via the results in [7, 10, 17, 26].

##  5 Applications to some special graphs

In this section, as applications of Theorem 3.1, we investigate the biharmonic
distance of some special graphs including the complement of a graph, the
Cartesian product of two graphs and the Cayley graph of a finite abelian
group. In particular, the combinatorial expression of biharmonic distance
between any two vertices of hypercubes are determined.

###  5.1 The complement of a graph

Suppose that $\lambda_{1}\leq\cdots\leq\lambda_{n}$ are the eigenvalues of the
Laplacian matrix of $G$ and $z_{1},\ldots,z_{n}$ are the corresponding
mutually orthogonal real unit eigenvectors.

###### Theorem 5.1

Let $G^{c}$ be the complement of a graph $G$ with $n$ vertices. If $G^{c}$ is
connected, then the biharmonic distance between two vertices $u$ and $v$ in
$G^{c}$ is

$d_{B}(u,v)=(\sum_{k=2}^{n}(n-\lambda_{k})^{-2}(z_{k}(u)-z_{k}(v))^{2})^{\frac{1}{2}}.$

Proof. Note that $L(G^{c})=nI-L(G)-J$, where $J$ denotes the all-$1$ matrix.
Then $L(G^{c})z_{1}=\mathbf{0}$ and $L(G^{c})z_{k}=(n-\lambda_{k})z_{k}$ for
$k=2,\ldots,n$. Note that $\lambda_{n}+\lambda_{2}(G^{c})=n$. Since $G^{c}$ is
connected, $\lambda_{2}(G^{c})>0$. Thus, $n-\lambda_{k}\geq
n-\lambda_{n}=\lambda_{2}(G^{c})>0$ for $2\leq k\leq n$. By Theorem 3.1, we
have
$d_{B}(u,v)=(\sum_{k=2}^{n}(n-\lambda_{k})^{-2}(z_{k}(u)-z_{k}(v))^{2})^{\frac{1}{2}}$.
$\hfill\Box\vskip 6.0pt plus 2.0pt minus 2.0pt$

###  5.2 The Cartesian product of two graphs

The _Cartesian product_ of two graphs $G_{1}$ and $G_{2}$ is the graph
$G_{1}\square G_{2}$ whose vertex set is the set $\\{u_{1}u_{2}|u_{1}\in
V(G_{1}),u_{2}\in V(G_{2})\\}$, and two vertices $u_{1}u_{2},v_{1}v_{2}$ are
adjacent if $u_{1}=v_{1}$ and $\\{u_{2},v_{2}\\}\in E(G_{2})$ or if
$u_{2}=v_{2}$ and $\\{u_{1},v_{1}\\}\in E(G_{1})$. The _Kronecker product_
$A\otimes B$ of matrices $A=(a_{ij})_{m\times n}$ and $B=(b_{ij})_{p\times q}$
is the $mp\times nq$ matrix obtained from $A$ by replacing each element
$a_{ij}$ with the block $a_{ij}B$.

Suppose that $\lambda_{1}(G_{i})\leq\cdots\leq\lambda_{n_{i}}(G_{i})$ are the
eigenvalues of the Laplacian matrix of $G_{i}$ and
$z_{1}(G_{i}),\ldots,z_{n}(G_{i})$ are the corresponding mutually orthogonal
real unit eigenvectors for $i=1,2$.

######  Lemma 5.2 ([12])

Let $G_{1}\square G_{2}$ be the Cartesian product of two graphs $G_{1}$ and
$G_{2}$. Then the eigenvalues of the Laplacian matrix of $G_{1}\square G_{2}$
are $\lambda_{i}(G_{1})+\lambda_{j}(G_{2})$ for $i=1,\ldots,n_{1}$ and
$j=1,\ldots,n_{2}$ and the corresponding eigenvectors are $z_{i}(G_{1})\otimes
z_{j}(G_{2})$ for $i=1,\ldots,n_{1}$ and $j=1,\ldots,n_{2}$.

By Theorem 3.1 and Lemma 5.2, we immediately obtain the following result.

###### Theorem 5.3

Let $G_{1}\square G_{2}$ be the Cartesian product of two graphs $G_{1}$ and
$G_{2}$. If $G_{1}\square G_{2}$ is connected, then the biharmonic distance
between two vertices $u_{1}u_{2}$ and $v_{1}v_{2}$ is

$(\sum_{1\leq i\leq n_{1},1\leq j\leq
n_{2}}(\lambda_{i}(G_{1})+\lambda_{j}(G_{2}))^{-2}(z_{i}(G_{1})\otimes
z_{j}(G_{2})(u_{1}u_{2})-z_{i}(G_{1})\otimes
z_{j}(G_{2})(v_{1}v_{2}))^{2})^{\frac{1}{2}}.$

###  5.3 The Cayley graph of a finite abelian group

Let $\Gamma$ be a group and let $S$ be a subset of $\Gamma$ that is closed
under taking inverses and does not contain the identity. The Cayley graph
$G(\Gamma,S)$, is the graph with vertex set $\Gamma$ and edge set
$\\{\\{g,h\\}\mid g^{-1}h\in S\\}$. A character $\chi$ of $\Gamma$ is a group
homomorphism $\chi:\Gamma\rightarrow C^{*}$. Let $\Gamma^{\wedge}$ denote the
set of characters of $\Gamma$. Denote the trivial character by $1_{\chi}$.

######  Lemma 5.4 ([1])

Let $\Gamma$ be a finite abelian group. Then the adjacency eigenvalues of
Cayley graph $G(\Gamma,S)$ are $\alpha_{\chi}$ with $\chi\in\Gamma^{\wedge}$,
where $\alpha_{\chi}=\sum_{s\in S}\chi(s)$. In addition, $\chi$ is an
adjacency eigenvector with eigenvalue $\alpha_{\chi}$.

Note that the Cayley graph $G(\Gamma,S)$ is a regular graph and $L=|S|I-A$,
where $L$ and $A$ are the Laplacian matrix and the adjacency matrix of
$G(\Gamma,S)$ respectively. By Theorem 3.1 and Lemma 5.4, we obtain the
following result.

###### Theorem 5.5

Let $\Gamma$ be a finite abelian group. Then the biharmonic distance between
two vertices $u$ and $v$ in Cayley graph $G(\Gamma,S)$ is

$d_{B}(u,v)=(\sum_{1_{\chi}\neq\chi\in\Gamma^{\wedge}}(|S|-\alpha_{\chi})^{-2}(\chi(u)-\chi(v))^{2})^{\frac{1}{2}}.$

The $n$-dimensional hypercube $Q_{n}$ [19] is a graph with vertex set
$\\{x_{1}x_{2}\ldots x_{n}\mid x_{i}\in\\{0,1\\},1\leq i\leq n\\}$ and two
vertices are adjacent if and only if they differ exactly in one position (see
Figure 2).

Figure 2: The structures of $Q_{2}$ and $Q_{3}$.

###### Remark 5.1

Let $(Z_{2})^{n}$ be the direct products of the cyclic group $Z_{2}=\\{0,1\\}$
with $|(Z_{2})^{n}|=2^{n}$. We use the notation $v_{i}$ to denote the $i$th
coordinate of $v\in(Z_{2})^{n}$, where $v=v_{1}v_{2}\ldots v_{n}$ and
$v_{k}\in Z_{2}$ for $1\leq k\leq n$. Note that the $n$-dimensional hypercube
$Q_{n}$ can be viewed as a Cayley graph $G(\Gamma,S)$ satisfying that
$\Gamma=(Z_{2})^{n}$ and $S=\\{e^{i}\mid 1\leq i\leq n\\}$, where
$e^{i}_{h}=1$ if and only if $h=i$. Since the characters of $(Z_{2})^{n}$ are
$\chi_{I}(v)=(-1)^{\sum_{i\in I}v_{i}}$ for $I\subseteq\\{1,\ldots,n\\}$, by
Lemma 5.4, the adjacency eigenvalue corresponding to $\chi_{I}$ is
$\alpha_{I}=\sum_{j=1}^{n}\chi_{I}(e^{j})=n-2|I|$.

###### Corollary 5.6

The biharmonic distance between two vertices $u$ and $v$ in $Q_{n}$ is

$d_{B}(u,v)=\dfrac{\sqrt{2}}{2}(\sum_{\emptyset\neq
I\subseteq\\{1,\ldots,n\\}}|I|^{-2}(1-(-1)^{\sum_{i\in
I}(u_{i}+v_{i})}))^{\frac{1}{2}}.$

Proof. By Theorem 5.5 and Remark 5.1, we have

$\displaystyle d_{B}(u,v)$ $\displaystyle=$ $\displaystyle(\sum_{\emptyset\neq
I\subseteq\\{1,\ldots,n\\}}(|S|-\alpha_{I})^{-2}(\chi_{I}(u)-\chi_{I}(v))^{2})^{\frac{1}{2}}$
$\displaystyle=$ $\displaystyle(\sum_{\emptyset\neq
I\subseteq\\{1,\ldots,n\\}}(n-(n-2|I|))^{-2}((-1)^{\sum_{i\in
I}u_{i}}-(-1)^{\sum_{i\in I}v_{i}})^{2})^{\frac{1}{2}}$ $\displaystyle=$
$\displaystyle(\sum_{\emptyset\neq
I\subseteq\\{1,\ldots,n\\}}(2|I|)^{-2}(2-2(-1)^{\sum_{i\in
I}(u_{i}+v_{i})}))^{\frac{1}{2}}$ $\displaystyle=$
$\displaystyle\dfrac{\sqrt{2}}{2}(\sum_{\emptyset\neq
I\subseteq\\{1,\ldots,n\\}}|I|^{-2}(1-(-1)^{\sum_{i\in
I}(u_{i}+v_{i})}))^{\frac{1}{2}}.$

This completes the proof of Corollary 5.6. $\hfill\Box\vskip 6.0pt plus 2.0pt
minus 2.0pt$

##  6 Conclusions

In this paper, we first give some characterizations of the biharmonic distance
of a graph and obtain some bounds for it. Second, we reveal a relationship
between biharmonic index and Kirchhoff index and determine the unique graph
having the minimum biharmonic index among the connected graphs with given
order. Finally, as applications, we investigate the biharmonic distance of
some special graphs including the complement of a graph, the Cartesian product
of two graphs and the Cayley graph of a finite abelian group. In particular,
the combinatorial expression of biharmonic distance between any two vertices
of hypercubes are determined. Note that resistance distance has some
interpretations based on random walks and electrical networks. A natural
question is that whether the biharmonic distance has some similar
interpretations. Thus, we state a few challenging open problems on the
biharmonic distance of a graph.

###### Problem 6.1

Find a probability characterization of the biharmonic distance $d_{B}(u,v)$
between two vertices $u$ and $v$ for a graph $G$.

###### Problem 6.2

Find a combinatorial (or physical) characterization of the biharmonic distance
$d_{B}(u,v)$ between two vertices $u$ and $v$ for a graph $G$.

###### Problem 6.3

Find lower and upper bounds on the biharmonic index $B(T)$ for trees $T$ of
given order and characterize the extremal graphs.

## Acknowledgement

R-H. Li’s research is supported by the National Natural Science Foundation of
China (No. 62072034 and 61772346). W. Yang’s research is supported by the
National Natural Science Foundation of China (No. 11671296).

## References

  * [1] L. Babai, Spectra of Cayley graphs, J. Combin. Theory Ser. B, 27(2) (1979), 180–189.
  * [2] B. Bamieh, M.R. Jovanovic, P. Mitra, S. Patterson, Coherence in large-scale networks: Dimension-dependent limitations of local feedback, IEEE Trans. Autom. Control, 57(9) (2012), 2235–2249.
  * [3] R.B. Bapat, Graphs and matrices, Universitext, Springer, London, 2010.
  * [4] P. Barooah, J.P. Hespanha, Graph effective resistances and distributed control: Spectral properties and applications, in Proceedings of the 45th IEEE Conference on Decision and Control, San Diego, CA, USA, 2006, 3479–3485.
  * [5] J.A. Bondy, U.S.R. Murty, Graph Theory with Applications, The Macmillan Press Ltd, New York, 1976.
  * [6] H. Chen, F. Zhang, Resistance distance and the normalized Laplacian spectrum, Discrete Appl. Math., 155(5) (2007), 654–661.
  * [7] X. Chen, K.C. Das, Characterization of extremal graphs from Laplacian eigenvalues and the sum of powers of the Laplacian eigenvalues of graphs, Discrete Math., 338(7) (2015), 1252–1263.
  * [8] D. Cvetković, P. Rowlinson, S. Simić, An introduction to the theory of graph spectra, Cambridge University Press, Cambridge, 2010.
  * [9] K.C. Das, A sharp upper bound for the number of spanning trees of a graph, Graphs Combin., 23(6) (2007), 625–632.
  * [10] K.C. Das, K. Xu, M. Liu, On sum of powers of the Laplacian eigenvalues of graphs, Linear Algebra Appl., 439(11) (2013), 3561–3575.
  * [11] P.G. Doyle, J.L. Snell, Random walks and electric networks, Mathematical Association of America, Washington, DC, 1984.
  * [12] M. Fiedler, Algebraic connectivity of graphs, Czechoslovak Math. J., 23(2) (1973), 298–305.
  * [13] K. Fitch, N.E. Leonard, Joint centrality distinguishes optimal leaders in noisy networks, IEEE Trans. Control Netw. Syst., 3(4) (2016), 366–378.
  * [14] I. Gutman, B. Mohar, The quasi-Wiener and the Kirchhoff indices coincide, J. Chem. Inf. Comput. Sci., 36(5) (1996), 982–985.
  * [15] D.J. Klein, M. Randić, Resistance distance, J. Math. Chem., 12(1) (1993), 81–95.
  * [16] Y. Lipman, R.M. Rustamov, T.A. Funkhouser, Biharmonic distance, ACM Transactions on Graphics, 29(3) (2010), 1–11.
  * [17] M. Liu, B. Liu, A note on sum of powers of the Laplacian eigenvalues of graphs, Appl. Math. Lett., 24(3) (2011), 249–252.
  * [18] Y.J. Peng, S.C. Li, On the Kirchhoff index and the number of spanning trees of linear phenylenes, MATCH Commun. Math. Comput. Chem, 77(3) (2017), 765–780.
  * [19] Y. Saad, M.H. Schultz, Topological properties of hypercubes, IEEE Trans. Comput., 37(7) (1988), 867–872.
  * [20] W. Stevenson, Elements of Power System Analysis, McGraw Hill, New York, 1975.
  * [21] W. Xu, B. Wu, Z. Zhang, Z. Zhang, H. Kan, G. Chen, Coherence scaling of noisy second-order scale-free consensus networks, IEEE Trans. Cybern., 52(7) (2022), 5923–5934.
  * [22] Y. Yang, D.J. Klein, A recursion formula for resistance distances and its applications, Discrete Appl. Math., 161(16-17) (2013), 2702–2715.
  * [23] Y. Yi, B. Yang, Z. Zhang, S. Patterson, Biharmonic distance and the performance of second-order consensus networks with stochastic disturbances, in 2018 Annual American Control Conference (ACC), Milwaukee, WI, USA, 2018, 4943–4950.
  * [24] Y. Yi, B. Yang, Z. Zhang, Z. Zhang, S. Patterson, Biharmonic distance-based performance metric for second-order noisy consensus networks, IEEE Trans. Inf. Theory, 68(2) (2022), 1220–1236.
  * [25] Z. Zhang, W. Xu, Y. Yi, Z. Zhang, Fast approximation of coherence for second-order noisy consensus networks, IEEE Trans. Cybern., 52(1) (2022), 677–686.
  * [26] B. Zhou, On sum of powers of the Laplacian eigenvalues of graphs, Linear Algebra Appl., 429(8-9) (2008), 2239–2246.