††thanks: Correspondence should be addressed to Y.T and P.S.

]Department of Psychology & Tsinghua Laboratory of Brain and Intelligence,
Tsinghua University, Beijing, 100084, China Also at ]Tsien Excellence in
Engineering Program, Tsinghua University, Beijing, 100084, China

]Department of Physics, Tsinghua University, Beijing, 100084, China.

]Faculty of Health Sciences, City University of Macau, Macau, 999078, China.

]Faculty of Health Sciences, City University of Macau, Macau, 999078, China
Also at ]Faculty of Data Science, City University of Macau, Macau, 999078,
China.

# Simplex path integral and simplex renormalization group for high-order
interactions

Aohua Cheng<EMAIL_ADDRESS>[  [  Yunhui Xu<EMAIL_ADDRESS>[  Pei
Sun<EMAIL_ADDRESS>[  Yang Tian<EMAIL_ADDRESS>[  [

###### Abstract

Modern theories of phase transitions and scale-invariance are rooted in path
integral formulation and renormalization groups (RG). Despite the
applicability of these approaches on simple systems with only pairwise
interactions, they are less effective on complex systems with un-decomposable
high-order interactions (i.e., interactions among arbitrary sets of units). To
precisely characterize the universality of high-order interacting systems, we
propose simplex path integral and simplex renormalization group (SRG) as the
generalizations of classic approaches to arbitrary high-order and
heterogeneous interactions. We first formalize the trajectories of units
governed by high-order interactions to define path integrals on corresponding
simplexes based on a high-order propagator. Then we develop a method to
integrate out short-range high-order interactions in the momentum space,
accompanied by a coarse graining procedure functioning on the simplex
structure generated by high-order interactions. The proposed SRG, equipped
with a divide-and-conquer framework, can deal with the absence of ergodicity
arised from the sparse distribution of high-order interactions and renormalize
a system with intertwined high-order interactions on $p$-order according to
its properties on $q$-order ($p\leq q$). The associated scaling relation and
its corollaries support to differentiate among scale-invariant, weakly scale-
invariant, and scale-dependent systems across different orders. We have
validated our theory in multi-order scale-invariance verification, topological
invariance discovery, organizational structure identification, and information
bottleneck analysis to demonstrate its capacity for identifying intrinsic
properties of high-order interacting systems during system reduction.

††preprint: APS/123-QED

##  I Introduction

###  I.1 Unknowns about high-order interactions

Over the past decades, the studies on phase transition phenomena, especially
the non-equilibrium ones, in different interacting systems have accomplished
substantial progress [1, 2]. This progress should be credited to the
development of path integral [3, 4, 5] and renormalization group [6, 7]
theories, which significantly deepen our understanding of system dynamics and
provide a precise formulation of scaling and criticality.

However, a non-negligible theoretical vacancy can be found in the existing
path integral and renormalization group approaches if we subdivide
interactions into pairwise and high-order categories. As the name suggests, a
pairwise interaction only involves a pair of units and does not require the
participation of any other unit. This kind of interaction can be represented
by edges in networks [8, 9] and are implicitly used in the derivations of
classic path integrals and renormalization groups (e.g., see Ref. [10]). A
high-order interaction, on the contrary, is the mutual coupling among more
than two units [11, 12, 13]. While some high-order interactions are reducible
and can be decomposed into a group of pairwise interactions, most high-order
interactions are un-decomposable and in-equivalent to the direct sum of the
pairwise ones [9, 14]. For instance, the triplet collaboration among three
agents is not equivalent to the trivial sum of three pairs of individual
collaborations. These un-decomposable high-order interactions have intricate
effects on system dynamics [13, 14] and can not be trivially represented by
ordinary networks [12]. Although notable efforts have been devoted to studying
the optimal characterization of high-order interactions (e.g., using
simplicial complexes [15, 16, 17, 18] or hypergraphs [19, 20, 21]), there are
fewer works focusing on developing specialized path integral and
renormalization group theories for complex systems with high-order
interactions [22, 17]. The difficulty of proposing these specialized
frameworks arises from the intrinsic dependence on pairwise interactions when
researchers constrain the discrete replicas of interacting systems as
conventional networks. Once this constraint is relaxed, the intertwined
effects of newly included high-order interactions would make the analysis on
emerging dynamics highly non-trivial and even lead to novel results that can
not be directly derived by classic theories (e.g., see instances in diffusion
and random walks [16, 23, 20], social dynamics [24, 25, 26], and neural
dynamics [27, 28]).

###  I.2 Related works and remaining challenges

To suggest a possible direction for developing path integral and
renormalization group theories for high-order interactions, we summarize the
accomplishments and limitations of previous works.

To date, most progress of the computational implementations of path integrals
and renormalization groups is achieved on real systems with pairwise
interactions [29, 30, 10, 31, 32, 33]. Although box-covering methods [34, 35,
36, 37] may be the most straightforward ways for coarse graining while
preserving the organizational properties of a system represented by a network,
they depend on the assumption of internal fractal properties and convey no
information about system dynamics. Different from box-covering, spectral
method [38] pays more attention to inherent dynamics properties (e.g., random
walks) during coarse graining but may suffer from the obstacles of finding
super-units implied by small-world effects [39, 33, 10]. Aiming at resolving
the coexisting and correlated scales, a geometric renormalization group is
developed to identify the potential geometric scaling rather than path-
distance scaling [33]. Nevertheless, this approach essentially relies on the
_a-priori_ -knowledge of the hidden metric space of networks (e.g., the
hyperbolic space [40]). With the pursuit of realizing the metric-free coarse
graining whose iterability is not limited by small-world effects, a Laplacian
renormalization group (LRG) [10] has been derived from the diffusion on
networks, a continuum counterpart of the block-spin transformation [41], to
deal with the network topology, dynamics, and geometry encoded by the graph
Laplacian operator [42]. This approach has a natural connection with path
integrals as well [43]. Although these properties make the LRG a promising
foundation for in-depth explorations, this framework suffers from several non-
negligible limitations. First, there is no apparent correspondence between the
LRG and an appropriate approach for analyzing un-decomposable high-order
interactions because the Laplacian operator can not characterize polyadic
relations directly [11, 12, 13, 9, 14]. Second, the intrinsic dependence of
the LRG on the ergodic property of system dynamics reflected by network
connectivity limits the applicability of this approach in analyzing high-order
interactions, which are more sparsely distributed in the system than pairwise
interactions and do not necessarily ensure ergodicity.

There exist fewer pioneering works devoted to developing renormalization
groups for high-order interactions, among which, a notable framework is the
real space renormalization group proposed by Refs. [22, 17]. This approach is
initially defined on the Laplacian operators of the network skeletons of the
Apollonian [44] and the pseudo-fractal [45] simplicial complexes to calculate
spectral dimensions [22]. Then, it is generalized to the normalized up-
Laplacian operator [23, 16] of simplicial complexes to deal with high-order
spectra. Despite the effectiveness of this approach, its dependence on the
non-trivial manual derivations using the Gaussian model [46] and the
specialized application scope about spectral dimensions [22, 17] make it less
applicable to the real scenarios where a programmatic implementation is
demanded for renormalizing the complex systems generated by empirical data
automatically.

In sum, while representative works such as the LRG [10] have suggested a
promising way to renormalize complex systems, the appropriate path integral
and renormalization group theories for high-order interactions remain clouded
by numerous theoretical unknowns and technical difficulties.

Figure 1: System evolution governed by high-order interactions. (a) The
simplicial complexes of different orders are presented. A $n$-order simplex
manifests itself as a clique with $n+1$ units in the network and represents
the $n$-order interactions among these $n+1$ units (i.e., each $n$-order
interaction is formed by the simultaneous participation of $n+1$ units). As an
instance, the $1$-order simplex is a pairwise interaction. (b) An instance of
the interacting system with high-order interactions is presented, where there
are $200$ units whose pairwise interactions follow a Barabási-Albert network
(the number of edges between a new unit and existing units is $c=6$, see
definitions in Appendix A) [47, 48]. The units engaged with the first four
orders of interactions are extracted using simplicial complexes. (c) The
associated network sketches of the system shown in (b) are presented on each
order, where two units share an edge in the network sketch on $n$-order only
if they are engaged with a $n$-order interaction in the system. (d-e) The
corresponding multi-order Laplacian (d) and high-order path Laplacian (e)
operators of the system shown in (b) are illustrated in every order.

###  I.3 Our frameworks and contributions

To fill these gaps, we generalize classic path integral and renormalization
group theories from dyadic to polyadic relations. To realize such
generalizations, we are required to derive the appropriate definitions of the
canonical density operators (i.e., the functional over fields) [49, 50, 10]
and the coarse graining processes in real and moment spaces [6, 7] for high-
order interactions.

Our work suggests a way to define these concepts on simplexes, in which
simplexes can join any sub-set of units to characterize un-decomposable high-
order interactions among them (see Fig. 1(a) for instances and Appendix A for
definitions). We propose two kinds of high-order propagators as the
generalized path integral formulation of the abstract diffusion processes
driven by multi-order interactions, which further creates a system mode
description tantamount to the moment space. The coarse graining in the moment
space is implemented by progressively integrating out short-range system modes
corresponding to some certain order of interactions. Parallel to this process,
we propose a real space coarse graining procedure to reduce the simplex
structures associated with short-range high-order interactions, where we
suggest a way to overcome the dependence on ergodicity. The correspondence
between the moment and real spaces enables our framework to the system coarse
graining according to both organizational structure and internal dynamics.
Taken together, these definitions give birth to our simplex path integral and
simplex renormalization group (SRG) theories, which are applicable to
arbitrary orders of interactions.

The SRG can be used to uncover the potential scale invariance of high-order
interactions, whose non-trivial fixed point may suggest the existence of
critical phenomena. Different from conventional renormalization groups, the
SRG is designed to renormalize $p$-order interactions (each $p$-order
interaction requires the simultaneous participation of $p+1$ units) in the
system according to the structure and dynamics associated with $q$-order
interactions ($p\leq q$). Therefore, we can compare the renormalization flow
of $p$-order interactions guided by $p$-order interactions with that guided by
arbitrary higher-order interactions. This feature allows us to investigate the
impacts of higher-order interactions on lower-order ones and study the
differences between lower-order and higher-order interactions in terms of
scaling properties. In a special case with $p=q=1$ (i.e., pairwise
interactions), the two kinds of high-order propagators developed in our work
both reduce to the classic counterpart used in the LRG [10] and our proposed
SRG can reproduce the outputs of the LRG in the ergodic case.

To comprehensively demonstrate the capacity of our approaches in identifying
the latent characteristics of diverse real complex systems, we validate our
framework from the perspectives of multi-order scale-invariance detection,
topological invariance discovery, latent organizational structure
identification, and information bottleneck optimization, respectively. An
efficient code implementation of this framework is provided in [51]. One can
also see a detailed explanation about the programmatic implementation and
usage of the SRG in Appendix K.

##  II Un-decomposable high-order interactions

The Laplacian operator is a natural choice for characterizing interactions and
their spreading among units [52, 53, 54]. In a special case with pairwise
interactions, the Laplacian has a classic expression,
$\mathbf{L}_{ij}=\delta_{ij}\sum_{k}\mathbf{A}_{ik}-\mathbf{A}_{ij}$, where
$\delta$ denotes the Kronecker delta function and $\mathbf{A}$ is a weighted
adjacency matrix that describes the non-negative coupling strengths among all
units in unit set $V$ (here the non-negativity is required by the Laplacian).
Under the autonomous condition, the evolution of system dynamics given an
initial state $\mathbf{x}_{0}$ is defined by
$\mathbf{x}=\exp\left(-\tau\mathbf{L}\right)\mathbf{x}_{0}$, where
$\tau\in\left(0,\infty\right)$ denotes a time scale [55]. However, in more
general cases with un-decomposable high-order interactions, the classic
Laplacian operator is no longer applicable. To fill the gap, previous studies
have explored the generalization of the Laplacian operator on uniform
hypergraphs [56, 57] and simplicial complexes [58, 59, 23, 9], where diverse
variants of the Laplacian (e.g., the combinatorial Laplacian [60], the Hodge
Laplacian [61, 23], and the multi-order Laplacian [9]) have been proposed to
characterize high-order interactions. In our work, we first derive a multi-
order Laplacian following the spirit of Ref. [9]. Then we develop a new
operator, referred to as the high-order path Laplacian, as an alternative
description of high-order interactions from a different perspective. As shown
in Fig. 1, these two operators are defined on clique complexes, a kind of
commonly used simplicial complexes in the persistent homology analysis [62,
60] (see Appendix A).

###  II.1 High-order interactions through faces

Let us consider the $n$-order interactions that manifest as the simultaneous
participation of $n+1$ related units. The un-decomposability of these
interactions lies in that the couplings among $n+1$ units are required to
occur simultaneously and none of them is dispensable. For instance, the
triplet collaboration requires the participation of all three agents together.
As discussed in Ref. [9], these interactions can be represented by the
$n$-faces of the $n$-simplex where $n+1$ units are placed on. For the sake of
clarity, we define $\\{i\\}_{n+1}\subset V$ as an arbitrary set of $n+1$ units
in which unit $i$ is included and denote $S^{\\{i\\}}_{n+1}$ as the set of all
permutations (i.e., rearrangements) on $\\{i\\}_{n+1}$. The multi-order
Laplacian $\mathbf{L}_{M}^{\left(n\right)}$ associated with $n$-simplexes is
defined as

$\displaystyle\left[\mathbf{L}_{M}^{\left(n\right)}\right]_{ij}=n\delta_{ij}\mathbf{D}_{i}^{\left(n\right)}-\mathbf{A}_{ij}^{\left(n\right)}.$
(1)

In Eq. (1), term $\mathbf{D}_{i}^{\left(n\right)}$ counts the number of
$n$-simplexes that contains unit $i$

$\displaystyle\mathbf{D}_{i}^{\left(n\right)}=\frac{1}{\left(n+1\right)!}\sum_{\\{i\\}_{n+1}\subset
V}\prod_{\omega\in S^{\\{i\\}}_{n+1}}\mathbf{M}_{\omega},$ (2)

where $\frac{1}{\left(n+1\right)!}$ serves the normalization for high-order
simplexes (i.e., the enumeration in Eq. (2) involves repetition in high-order
cases). Term $\mathbf{A}_{ij}^{\left(n\right)}$ counts the number of
$n$-simplexes that contains units $i$ and $j$

$\displaystyle\mathbf{A}_{ij}^{\left(n\right)}=\frac{1-\delta_{ij}}{\left(n+1\right)!}\sum_{\\{i,j\\}_{n+1}\subset
V}\prod_{\omega\in S^{\\{i,j\\}}_{n+1}}\mathbf{M}_{\omega},$ (3)

where $\\{i,j\\}_{n+1}\subset V$ is an arbitrary set of $n+1$ units in which
units $i$ and $j$ are included. The general definition of
$\mathbf{M}_{\omega}$ used in Eq. (2-3) is given as a product of a series of
matrices for any permutation $\omega=\left(k_{1}\ldots k_{n}\right)$

$\displaystyle\mathbf{M}_{\omega}=\mathbf{A}_{k_{1}k_{2}}\cdots\mathbf{A}_{k_{n-1}k_{n}},$
(4)

where $\mathbf{A}=\mathbf{A}^{\left(1\right)}$ is the classic adjacency matrix
of units.

In Eq. (2), vector $\mathbf{D}^{\left(n\right)}$ is similar to the degree
vector while matrix $\mathbf{A}^{\left(n\right)}$ can be understood as the
adjacency matrix. In a special case where $n=1$, Eq. (2) is equivalent to the
classic Laplacian for pairwise interactions. In the more general cases where
$n>1$, high-order interactions can be characterized by Eq. (2). Please see
Fig. 1(d) for the illustrations of the multi-order Laplacian. Although the
definitions of Eqs. (2-4) are different from the approach proposed by Ref.
[9], they are mathematically consistent with Ref. [9].

###  II.2 High-order interactions along paths

Then, we turn to another type of $n$-order interactions, which manifest as the
sequential actions of $n+1$ related units after they all participate in. The
un-decomposability of these interactions arises from that the mutual coupling
effects emerge only after $n+1$ units have been engaged with the interplay
progressively. For example, a pipeline-like collaboration may require all
agents to communicate with each other and behave in order. Each previous
action is globally known to all agents and affects all subsequent actions.
This property can be characterized by the paths (i.e., a sequence of
$1$-simplexes) that successively pass through the $n+1$ units placed on the
$n$-simplex. To describe the interactions propagating along the paths on
$n$-simplexes, we develop an operator named as the high-order path Laplacian
$\mathbf{L}_{H}^{\left(n\right)}$

$\displaystyle\left[\mathbf{L}_{H}^{\left(n\right)}\right]_{ij}=\frac{1}{n}\left(\delta_{ij}\mathbf{P}_{i}^{\left(n\right)}-\mathbf{B}_{ij}^{\left(n\right)}\right).$
(5)

In Eq. (5), coefficient $\frac{1}{n}$ denotes the normalization of the
propagation speed of interactions along paths. Term
$\mathbf{P}_{i}^{\left(n\right)}$ counts the number of paths that traverse all
$n+1$ units on the $n$-simplex without repetition and terminate at unit $i$

$\displaystyle\mathbf{P}_{i}^{\left(n\right)}=\frac{1}{n+1}\sum_{\\{i\\}_{n+1}\subset
V}\prod_{\omega\in S^{\\{i\\}}_{n+1}}\mathbf{M}_{\omega},$ (6)

where $\frac{1}{n+1}$ denotes the normalization. Term
$\mathbf{B}_{ij}^{\left(n\right)}$ counts the number of paths that traverse
all $n+1$ units on the $n$-simplex without repetition, initiate at unit $j$,
and terminate at unit $i$

$\displaystyle\mathbf{B}_{ij}^{\left(n\right)}=\frac{1-\delta_{ij}}{n\left(n+1\right)}\sum_{\\{i,j\\}_{n+1}\subset
V}\prod_{\omega\in S^{\\{i,j\\}}_{n+1}}\mathbf{M}_{\omega}.$ (7)

Please see Fig. 1(e) for the instances of the high-order path Laplacian. In
the case with $n=1$, Eqs. (6-7) reduce to the classic degree vector and
adjacency matrix, respectively (i.e., we have
$\mathbf{P}^{\left(1\right)}=\mathbf{D}^{\left(1\right)}$ and
$\mathbf{B}^{\left(1\right)}=\mathbf{A}$). Therefore, Eq. (5) denotes the
classic Laplacian when interactions are pairwise.

In general, one can choose the multi-order Laplacian
$\mathbf{L}_{M}^{\left(n\right)}$ and the high-order path Laplacian
$\mathbf{L}_{H}^{\left(n\right)}$ according to application demands. While
$\mathbf{L}_{M}^{\left(n\right)}$ conveys information about the high-order
interactions that units are simultaneously engaged with, the definition of
$\mathbf{L}_{H}^{\left(n\right)}$ is more applicable to the cases where a
high-order interaction is formed by an irreducible sequence of the
participation of related units. We offer efficient code implementations of
these two kinds of operators, which can be seen in [51]. Meanwhile, in
Appendix B, we explain the similarities and differences between our considered
Laplacian operators and other Laplacian operators used in topology theories
(e.g., the combinatorial Laplacian [60] and the Hodge Laplacian [61, 23] in
persistent homology analysis [63, 64, 62, 65, 66, 67]), which helps us
understand the theoretical significance of designing renormalization groups on
$\mathbf{L}_{M}^{\left(n\right)}$ and $\mathbf{L}_{H}^{\left(n\right)}$. For
convenience, we uniformly denote them as $\mathbf{L}^{\left(n\right)}$ in our
subsequent derivations. One can specify the actual definition of
$\mathbf{L}^{\left(n\right)}$ in the application.

Figure 2: Illustrations of the simplex path integral. (a) An interacting
system of $100$ units is generated, where the pairwise interactions of units
follow a Barabási-Albert network ($c=6$). This system is used to show the
derivation of the time evolution operators of the first three orders, where we
use the multi-order Laplacian (MOL) for calculation and set $\tau=1$. For each
order of interactions, we illustrate the first four terms of the series
expansion of time evolution operator in Eq. (5). (b) For the system defined in
(a), we show the density operators calculated based on the MOL or the high-
order path Laplacian (HOPL) on each order. (c) We show that values in the
density operator can be related to the evolution paths of system states. The
blue path leads from one state to another while the red path starts from and
ends at the same state.

##  III Simplex path integral

Given the high-order interactions described by $\mathbf{L}^{\left(n\right)}$,
we turn to explore the path integral formulation of the trajectories of units
governed by high-order interactions. As first suggested by Ref. [42] and then
discussed in Refs. [68, 69, 10, 70, 43], the Laplacian can be understood as a
Hamiltonian if we consider the system dynamics governed by it. In the case of
high-order interactions, we can consider the eigendecomposition,
$\mathbf{L}^{\left(n\right)}=\sum_{\lambda^{\left(n\right)}}\lambda^{\left(n\right)}\big{|}\lambda^{\left(n\right)}\big{\rangle}\big{\langle}\lambda^{\left(n\right)}\big{|}$,
using the bra-ket notation and define a time evolution operator

$\displaystyle\mathbf{U}^{\left(n\right)}_{\tau}=$
$\displaystyle\exp\left(-\tau\mathbf{L}^{\left(n\right)}\right)$ (8)
$\displaystyle=$
$\displaystyle\sum_{k=0}^{\infty}\frac{1}{k!}\left(-\tau\mathbf{L}^{\left(n\right)}\right)^{k},$
(9) $\displaystyle=$
$\displaystyle\sum_{\lambda^{\left(n\right)}}\exp\left(-\tau\lambda^{\left(n\right)}\right)\big{|}\lambda^{\left(n\right)}\big{\rangle}\big{\langle}\lambda^{\left(n\right)}\big{|}.$
(10)

As suggested by Eq. (9), this operator characterizes the probability for each
unit to evolve from one state to another state within a given time scale
$\tau$, i.e.,
$\big{|}\mathbf{x}\big{\rangle}=\mathbf{U}^{\left(n\right)}_{\tau}\big{|}\mathbf{x}_{0}\big{\rangle}$.
The characterization is realized by considering the distributions of units,
$\frac{\left(-\tau\mathbf{L}^{\left(n\right)}\right)^{k}}{k!}$, generated by
all the possible ways of evolution, $k\in\left(0,\infty\right)$, during a time
period of $\tau$. See Fig. 2(a) for illustrations. As indicated by Eq. (10),
the probability amplitude is intrinsically determined by the exponential decay
rates of diffusion modes in a moment space. Short-range interactions
associated with fast decays (i.e., large eigenvalues) in the moment space have
small impacts on units.

At equilibrium, the Gibbs state of the system associated with $n$-order
interactions is

$\displaystyle\mathbf{\rho}^{\left(n\right)}=\frac{\mathbf{U}^{\left(n\right)}_{\tau}}{\operatorname{tr}\left(\mathbf{U}^{\left(n\right)}_{\tau}\right)}=\frac{\sum_{\lambda^{\left(n\right)}}\exp\left(-\tau\lambda^{\left(n\right)}\right)\big{|}\lambda^{\left(n\right)}\big{\rangle}\big{\langle}\lambda^{\left(n\right)}\big{|}}{\sum_{\lambda^{\left(n\right)}}\exp\left(-\tau\lambda^{\left(n\right)}\right)},$
(11)

where $\operatorname{tr}\left(\cdot\right)$ denotes the trace and parameter
$\tau$ serves as the inverse of a finite temperature. Please see Fig. 2(b) for
instances. The Boltzmann constant is assumed as unitary for simplicity. The
partition function,
$\sum_{\lambda^{\left(n\right)}}\exp\left(-\tau\lambda^{\left(n\right)}\right)$,
is proportional to the average return probability of random walks within a
time scale $\tau$ [42]. See Figs. 2(b-c) for illustrations.

Eq. (11) defines a density operator of system evolution characterized by high-
order interactions

$\displaystyle\big{\langle}\mathbf{x}\big{|}\mathbf{\rho}^{\left(n\right)}\big{|}\mathbf{x}_{0}\big{\rangle}=\frac{1}{\operatorname{tr}\left(\mathbf{U}^{\left(n\right)}_{\tau}\right)}\big{\langle}\mathbf{x}\big{|}\mathbf{U}^{\left(n\right)}_{\tau}\big{|}\mathbf{x}_{0}\big{\rangle},$
(12)

which naturally relates to the path integral formulation. Let us consider a
minimum time step $\varepsilon$ such that $N=\frac{\tau}{\varepsilon}$. At the
limit $N\rightarrow\infty$ (or $\varepsilon\rightarrow 0$), the numerator of
the density operator (i.e., the time evolution operator) implies a power form
of the corresponding high-order propagator
$\mathbf{K}^{\left(n\right)}\left(\mathbf{x},\mathbf{x}_{0},\tau\right)=\big{\langle}\mathbf{x}\big{|}\mathbf{U}^{\left(n\right)}_{\tau}\big{|}\mathbf{x}_{0}\big{\rangle}=\lim_{N\rightarrow\infty}\big{\langle}\mathbf{x}\big{|}\left(\mathbf{U}^{\left(n\right)}_{\varepsilon}\right)^{N}\big{|}\mathbf{x}_{0}\big{\rangle}$.
On the one hand, this power form expression directly leads to a path integral
in conventional quantum field theory

$\displaystyle\mathbf{K}^{\left(n\right)}\left(\mathbf{x},\mathbf{x}_{0},\tau\right)$
$\displaystyle=$
$\displaystyle\lim_{N\rightarrow\infty}\int\prod_{k=0}^{N-1}\big{\langle}\mathbf{x}_{\left(k+1\right)\varepsilon}\big{|}\mathbf{U}^{\left(n\right)}_{\varepsilon}\big{|}\mathbf{x}_{k\varepsilon}\big{\rangle}\mathsf{d}\mathbf{x}_{\varepsilon}\cdots\mathsf{d}\mathbf{x}_{\left(N-1\right)\varepsilon},$
(13)

where we mark $\mathbf{x}=\mathbf{x}_{N\varepsilon}$. On the other hand, this
power form expression can also be reformulated using the imaginary time and
Eqs. (8-9)

$\displaystyle\mathbf{K}^{\left(n\right)}\left(\mathbf{x},\mathbf{x}_{0},\tau\right)=$
$\displaystyle\sum_{k=0}^{\infty}\frac{1}{k!}\Bigg{\langle}\mathbf{x}\Bigg{|}\left(-\frac{i}{\hbar}\int_{0}^{\tau}-i\hbar\mathbf{L}^{\left(n\right)}\mathsf{d}\varepsilon\right)^{k}\Bigg{|}\mathbf{x}_{0}\Bigg{\rangle},$
(14)

where notion $\hbar$ denotes the reduced Planck constant (i.e., or named as
the Dirac constant) [3, 4]. Eq. (14) is closely related to the imaginary time
evolution governed by operator $\mathbf{L}^{\left(n\right)}$ in the infinite-
dimensional path space. Specifically, the probability amplitude for the system
to evolve from $\big{|}\mathbf{x}_{0}\big{\rangle}$ to
$\big{|}\mathbf{x}\big{\rangle}$ during a period of $\tau$ conveyed by the
propagator in Eq. (13) can be measured by summing over all possible evolution
paths connecting between these two states (i.e., the discrete counterpart of
path integrals) in Eq. (14). This is the reason why the $n$-order network
topology conveyed by $\mathbf{L}^{\left(n\right)}$ intrinsically shapes path
integrals. See Fig. 2(c) for illustrations.

Moreover, the denominator of the density operator (i.e., the partition
function) can be interpreted by considering the expectation
$\big{\langle}\mathbf{U}^{\left(n\right)}_{\tau}\big{\rangle}=\operatorname{tr}\left(\rho^{\left(n\right)}\mathbf{U}^{\left(n\right)}_{\tau}\right)$

$\displaystyle\operatorname{tr}\left(\rho^{\left(n\right)}\mathbf{U}^{\left(n\right)}_{\tau}\right)$
$\displaystyle=$
$\displaystyle\int\big{\langle}\mathbf{x}_{0}\big{|}\mathbf{U}^{\left(n\right)}_{\tau}\big{|}\mathbf{x}_{0}\big{\rangle}\mathsf{d}\mathbf{x}_{0},$
(15) $\displaystyle=$
$\displaystyle\int\Bigg{\langle}\mathbf{x}_{0}\Bigg{|}\left(\oint_{\Gamma}\exp\left(-\frac{i}{\hbar}\int_{0}^{\tau}-i\hbar\mathbf{L}^{\left(n\right)}\mathsf{d}\varepsilon\right)\mathsf{d}\mathbf{x}_{\varepsilon}\right)\Bigg{|}\mathbf{x}_{0}\Bigg{\rangle}\mathsf{d}\mathbf{x}_{0},$
(16)

where $\Gamma$ is an arbitrary closed curve that starts from and ends at
$\mathbf{x}_{0}$. Eq. (16) suggests why the partition function is proportional
to the average return probability of random walks within a time scale $\tau$
from the perspective of path integral formulation. See Fig. 2(c) for
illustrations.

The path integrals in Eqs. (13-16) are fully characterized by
$\mathbf{L}^{\left(n\right)}$, which describes latent simplex structures of
the system. They are integrals over all the possible evolution paths of a
system governed by high-order interactions. Consequently, they are referred to
as the simplex path integrals in our work. Although our framework is presented
using the terminologies of quantum field theory [3, 4] for convenience, it is
generally applicable to diverse classical systems (e.g., the brain in
neuroscience) in the world. In fact, we can naturally find the Maxwell-
Boltzmann statistics, a common description of classical systems, in Eqs.
(11-12).

##  IV Simplex renormalization group

The connection between Eqs. (8-16) and the moment space naturally inspires us
to consider the generalization of renormalization groups. The freedom degrees
corresponding to short-range interactions (e.g., fast diffusion within the
clusters of strongly correlated units) can be safely coarse grained without
significant information loss, which is consistent with our interpretation of
Eq. (10). Here we develop a framework to generalize renormalization groups to
high-order interactions, where we also propose possible ways to overcome the
limitations of the Laplacian renormalization group (LRG) [10].

###  IV.1 Renormalization procedure

To offer a clear vision, we first introduce our approach, referred to as the
simplex renormalization group (SRG), in an ergodic case as assumed by Ref.
[10]. The SRG renormalizes a system on the $p$-order according to the
properties of $q$-order interactions ($p\leq q$). Specifically, the
renormalization is realized as the following:

Figure 3: Illustrations of the simplex renormalization group (SRG). (a) The
conceptual illustration of the SRG is presented. For the sake of legibility,
the simplicial complexes of different orders are not elaborated and are
abstractly represented by dashed lines with different colors. Given a system
where multi-order interactions coexist, the SRG can be applied to renormalize
the system on the $p$-order based on the structure and dynamics associated
with $q$-order interactions ($p\leq q$). (b) Conceptual illustrations of our
divide-and-conquer procedure for dealing with the dependence of
renormalization flows on the ergodicity. (c) A system with $500$ units is
defined, where pairwise interactions follow a Barabási-Albert network ($c=2$).
We use this system to illustrate the selection of time scale
$\tau^{\left(q\right)}$ based on the specific heat ($q=2$).

  * (1)

In each $k$-th iteration, we have two iterative Laplacian operators,
$\mathbf{L}^{\left(p\right)}_{k}$ and $\mathbf{L}^{\left(q\right)}_{k}$, and
their associated high-order network sketches,
$\mathbf{G}^{\left(p\right)}_{k}$ and $\mathbf{G}^{\left(q\right)}_{k}$, with
$N^{\left(p\right)}_{k}=N^{\left(q\right)}_{k}$ units (there is an edge
connecting between two units in $\mathbf{G}^{\left(p\right)}_{k}$ or
$\mathbf{G}^{\left(q\right)}_{k}$ only if they share a $p$\- or $q$-order
interaction).

  * (2)

We find the eigenvalues of $\mathbf{L}^{\left(q\right)}_{k}$ that are smaller
than $\frac{1}{\widehat{\tau}^{\left(q\right)}}$, where
$\widehat{\tau}^{\left(q\right)}$ is the time scale for $q$-order simplicial
complexes (see Sec. IV.2 for the selection of
$\widehat{\tau}^{\left(q\right)}$). These eigenvalues correspond to long-range
high-order interactions because the contributions of their eigenvectors have
slow decays in Eq. (10). We denote the number of these eigenvalues as
$m_{k}^{\left(q\right)}$. Other eigenvalues correspond to short-range high-
order interactions and can be coarse grained. Moreover, according to Eq. (11),
the totality of long-range high-order interactions (i.e., corresponding to
smaller eigenvalues) is generally greater than that of short-range ones (i.e.,
associated with larger eigenvalues). Therefore, integrating out short-range
high-order interactions does not damage the effective information contained in
the system significantly.

  * (3)

In the moment space, we reduce the modes with fast decays in the system. The
current $q$-order Laplacian, $\mathbf{L}^{\left(q\right)}_{k}$, is reduced to
the re-scaled contributions of the eigenvectors associated with long-range
high-order interactions

$\displaystyle\mathbf{Q}^{\left(q\right)}_{k}=\sum_{\lambda^{\left(q\right)}_{k}}\Theta\left(\frac{1}{\widehat{\tau}^{\left(q\right)}}-\lambda^{\left(q\right)}_{k}\right)\lambda^{\left(q\right)}_{k}\big{|}\lambda^{\left(q\right)}_{k}\big{\rangle}\big{\langle}\lambda^{\left(q\right)}_{k}\big{|},$
(17)

where $\Theta\left(\cdot\right)$ is the Heaviside step function (i.e.,
$\Theta\left(x\right)=1$ if $x\geq 0$ and $\Theta\left(x\right)=0$ if $x<0$).
Meanwhile, the current $p$-order Laplacian, $\mathbf{L}^{\left(p\right)}_{k}$,
is reduced to the contributions of the $m_{k}^{\left(q\right)}$ smallest
eigenvalues (i.e., long-range $p$-order interactions)

$\displaystyle\mathbf{Q}^{\left(p\right)}_{k}=\sum_{\omega^{\left(p\right)}_{k}}I_{\Omega^{\left(p\right)}_{k}}\left(\omega^{\left(p\right)}_{k}\right)\omega^{\left(p\right)}_{k}\big{|}\omega^{\left(p\right)}_{k}\big{\rangle}\big{\langle}\omega^{\left(p\right)}_{k}\big{|},$
(18)

where $\Omega^{\left(p\right)}_{k}$ is the set of the $m_{k}^{\left(q\right)}$
smallest eigenvalues of $\mathbf{L}^{\left(p\right)}_{k}$ and
$I_{A}\left(\cdot\right)$ is the indicator function defined on an arbitrary
set $A$. Note that the eigenvalues in $\Omega^{\left(p\right)}_{k}$ are not
necessarily smaller than $\frac{1}{\tau^{\left(q\right)}}$.

  * (4)

In the real space, we search for every cluster of units sharing short-range
high-order interactions within the time scale
$\widehat{\tau}^{\left(q\right)}$, such that we can coarse grain these units
into a macro-unit. To this end, we generate a reference network,
$\mathbf{H}^{\left(q\right)}_{k}$, in the following way. We first initialize
$\mathbf{H}^{\left(q\right)}_{k}$ as a null network (i.e., with no edge). We
follow Eq. (11) to calculate the density operator associated with time scale
$\widehat{\tau}^{\left(q\right)}$ as
$\widehat{\rho}^{\left(q\right)}_{k}=\exp\left(-\widehat{\tau}^{\left(q\right)}\mathbf{L}^{\left(q\right)}_{k}\right)/\operatorname{tr}\left(\exp\left(-\widehat{\tau}^{\left(q\right)}\mathbf{L}^{\left(q\right)}_{k}\right)\right)$.
In operator $\widehat{\rho}^{\left(q\right)}_{k}$, its main diagonal elements,
related to self-interactions of units, are frequently greater than most of its
off-diagonal elements (e.g., see Fig. 2(b)). If the $\left(i,j\right)$-th
element is greater than the $\left(i,i\right)$-th and $\left(j,j\right)$-th
elements, then the $q$-order interactions between units $i$ and $j$ are
sufficiently quick such that their accumulations within time scale
$\widehat{\tau}^{\left(q\right)}$ are faster than those of self-interactions.
In this case, we add an edge between units $i$ and $j$ in network
$\mathbf{H}^{\left(q\right)}_{k}$. By progressively implement the edge adding
procedure for every pair of units, we can form $n_{k}^{\left(q\right)}$
clusters in $\mathbf{H}^{\left(q\right)}_{k}$. Each cluster contains a set of
strongly correlated units with short-range high-order interactions. In other
words, these clusters can be understood as the Kadanoff blocks (i.e., the
quasi-independent blocks whose mean size serves as the correlation length of
the system) in the real space renormalization procedure [71, 72, 73]. In
$\mathbf{G}^{\left(q\right)}_{k}$, each set of units that belong to the same
cluster in $\mathbf{H}^{\left(q\right)}_{k}$ are aggregated into a super-unit
to generate a new network $\mathbf{G}^{\left(q\right)}_{k+1}$. Consequently,
there remain $n^{\left(q\right)}_{k}$ super-units in
$\mathbf{G}^{\left(q\right)}_{k+1}$ after coarse graining. Parallel to this
process, we also aggregate the units belonging to the same cluster of
$\mathbf{H}^{\left(q\right)}_{k}$ into a super-unit in
$\mathbf{G}^{\left(p\right)}_{k}$. The generated
$\mathbf{G}^{\left(p\right)}_{k+1}$ contains $n^{\left(q\right)}_{k}$ super-
units as well.

  * (5)

The correspondence between the moment space and the real space for $q$-order
interactions is ensured by a Laplacian $\mathbf{L}^{\left(q\right)}_{k+1}$.
This operator is generated following
$\left(\mathbf{L}^{\left(q\right)}_{k+1}\right)_{\nu_{i}\nu_{j}}=\big{\langle}\nu_{i}\big{|}\mathbf{Q}^{\left(q\right)}_{k}\big{|}\nu_{j}\big{\rangle}$,
where $\nu_{i}$ and $\nu_{j}$ are super-units in
$\mathbf{G}^{\left(q\right)}_{k+1}$. Notion $\big{|}\nu_{j}\big{\rangle}$ is a
$N^{\left(q\right)}_{k}$-dimensional ket where unitary components correspond
to all the units in $\mathbf{G}^{\left(q\right)}_{k}$ that are aggregated in
$\nu_{i}$ and zero components correspond to all other units in
$\mathbf{G}^{\left(q\right)}_{k+1}$ that are not covered by $\nu_{i}$.
Mathematically, this procedure actually defines a similarity transformation
$T$ between $\mathbf{Q}^{\left(q\right)}_{k}$ and
$\mathbf{L}^{\left(q\right)}_{k+1}$

$\displaystyle
T^{-1}\mathbf{Q}^{\left(q\right)}_{k}T=\operatorname{diag}\left(\left[\mathbf{L}^{\left(q\right)}_{k+1},\mathbf{0}\right]\right),$
(19)

where $\mathbf{0}$ is a
$\left(N^{\left(q\right)}_{k}-n^{\left(q\right)}_{k}\right)$-dimensional all-
zero square matrix. It is clear that the similarity transformation is given as
$T=\left(\big{|}\nu_{1}\big{\rangle},\ldots,\big{|}\nu_{n^{\left(q\right)}_{k}}\big{\rangle},\big{|}\eta_{n^{\left(q\right)}_{k}+1}\big{\rangle},\ldots,\big{|}\eta_{N^{\left(q\right)}_{k}}\big{\rangle}\right)$,
where each ket $\eta_{j}$ is selected such that the columns of $T$ define a
group of orthonormal bases. Similarly, we can derive the Laplacian
$\mathbf{L}^{\left(p\right)}_{k+1}$ for $p$-order interactions based on
$\mathbf{Q}^{\left(p\right)}_{k}$. The derived
$\mathbf{L}^{\left(p\right)}_{k+1}$ and $\mathbf{L}^{\left(q\right)}_{k+1}$
are used in the $\left(k+1\right)$-th iteration.

The renormalization of $p$-order interactions in both the real and moment
spaces is realized under the guidance of $q$-order interactions (e.g., see
steps (3-4) for details). After renormalizing
$\mathbf{L}^{\left(p\right)}_{k}$ by repeating steps (1-5), the SRG
progressively drives the system towards an intrinsic scale of $p$-order
interactions that exceeds the microscopic scales. More specifically, the SRG
creates a $q$-order-guided iterative scale transformation that, irrespective
of starting from which local scale, always leads to the latent critical point
associated with $p$-order interactions (if there is any). It provides
opportunities to verify whether the concerned thermodynamic functions are
scale-invariant or not. See Fig. 3(a) for a summary of key steps in our
approach. Note that the using of a reference network in step (4) is not
necessary for theoretical derivations, yet it is favorable for computer
programming.

In a non-ergodic case, the presented procedure can not be directly adopted.
This is because our proposed renormalization approach is rooted in the
relations among the time evolution operators defined by Eqs. (8-10), the
exponential decay of diffusion modes in the moment space, and the evolution
paths in the path space, which implicitly depends on the assumption of the
ergodicity of system states. In the non-ergodic case where a set of states can
never evolve to or be transformed from another set of states (i.e., the
network connectivity is absent), path integrals between these state sets are
ill-posed and the diffusion over the network can be subdivided into several
irrelevant sub-processes (i.e., the diffusion processes defined on two
disconnected clusters are independent from each other). Dealing with these
properties is crucial for our work since the SRG is oriented toward analyzing
high-order interactions, which are usually sparsely distributed and may lack
ergodicity.

Figure 4: Instances of the renormalization flows of the simplex
renormalization group (SRG). The first three iterations of renormalization on
the systems whose pairwise interactions follow the random tree (a), the square
lattice (b), the triangular lattice (c), and the Watts-Strogatz networks (d-e)
are presented, where we set $\left(p,q\right)=\left(1,1\right)$,
$\left(p,q\right)=\left(1,2\right)$ and $\left(p,q\right)=\left(1,3\right)$
for (a-b), (c-d) and (e), respectively.

Here we develop a divide-and-conquer approach. Our idea arises from the fact
that we can treat the diffusion from one cluster (i.e., a connected component)
to another cluster on any order as being decelerated to a condition with zero
velocity (i.e., being infinitely slow). This infinitely long-term diffusion
process contains no short-range high-order interactions and, therefore, should
not be integrated out during renormalization. Consequently, the SRG does not
need to deal with the diffusion across any pair of clusters and can focus on
the diffusion within each cluster severally. In the $k$-iteration, we consider
the following two cases:

  * (A)

If $\mathbf{G}^{\left(q\right)}_{k}$ is connected (i.e., the ergodicity
holds), we directly apply steps (1-5) on $\mathbf{G}^{\left(q\right)}_{k}$ to
guide the renormalization of $p$-order interactions and derive the inputs of
the $\left(k+1\right)$-iteration.

  * (B)

If $\mathbf{G}^{\left(q\right)}_{k}$ is disconnected (i.e., the ergodicity
does not hold) and has $r$ clusters denoted by $C_{1},\;\ldots,C_{r}$, we
respectively apply steps (1-5) on each cluster of
$\mathbf{G}^{\left(q\right)}_{k}$ to guide the renormalization of $p$-order
interactions formed among the units that belong to this cluster. Specifically,
we treat each cluster $C_{i}$ as a network and input it into steps (1-5) with
$\mathbf{G}^{\left(p\right)}_{k}$. Please note that the above procedure does
not affect $\mathbf{G}^{\left(p\right)}_{k}$ if the selected cluster $C_{i}$
is trivial (i.e., contains only one unit). After dealing with all $r$
clusters, the obtained results are used for the $\left(k+1\right)$-iteration.

In Fig. 3(b), we conceptually illustrate our approach, which enables us to
relax the constraint on system ergodicity while implementing the SRG.

In Appendix B, we explain how the SRG is related to the persistent homology
theory [62, 74, 75, 76] and suggest its potential applicability in analyzing
the lifetime of topological properties. In Appendix C, we discuss the SRG in a
form of conventional renormalization groups, suggesting the benefits of
including high-order interactions into renormalization group analysis and
indicate the relations between the SRG and other frameworks [31, 77, 78, 79].
Specifically, we find that the ideas underlying the SRG can be clearly
interpreted by relating the SRG with a PCA-like renormalization group theory
[31].

###  IV.2 Critical time scale

After developing the renormalization procedure, we attempt to present a
detailed analysis of the time scale $\tau^{\left(q\right)}$ and its associated
scaling relation.

In the ergodic case, the setting of $\tau^{\left(q\right)}$ in the SRG (as we
have mentioned in step (2)) is finished in $1$-st iteration and required to
make the specific heat, an indicator of the transitions of diffusion
processes, maximize (see similar ideas in Ref. [10]). In Appendix C, we
explain why it is valid to define renormalization flows based on the specific
heat. Specifically, we consider the $q$-order spectral entropy in the $1$-st
iteration of the SRG (see Ref. [42] for the definition of spectral entropy)

$\displaystyle S^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)$
$\displaystyle=$
$\displaystyle-\operatorname{tr}\left[\rho^{\left(q\right)}_{1}\log\left(\rho^{\left(q\right)}_{1}\right)\right],$
(20) $\displaystyle=$
$\displaystyle\log\left(\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\right)+\tau^{\left(q\right)}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho},$
(21)

where
$\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho}=\operatorname{tr}\left(\mathbf{L}^{\left(q\right)}_{1}\rho^{\left(q\right)}_{1}\right)$.
The $q$-order specific heat is defined according to the first derivative of
the $q$-order spectral entropy

$\displaystyle X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)$
$\displaystyle=$
$\displaystyle-\frac{\mathsf{d}}{\mathsf{d}\log\left(\tau^{\left(q\right)}\right)}S^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right),$
(22) $\displaystyle=$
$\displaystyle-\left(\tau^{\left(q\right)}\right)^{2}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho}.$
(23)

Please see Appendix D for detailed derivations. The value of
$\tau^{\left(q\right)}$ is selected to ensure
$\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)=0$,
which indicates an infinite deceleration condition of diffusion processes. See
Fig. 3(c) for illustrations. Apart from the infinite deceleration condition,
determining the time scale value $\widehat{\tau}^{\left(q\right)}$ used in the
step (2) of the SRG also requires us to consider the high-order properties of
diffusion. For convenience, we denote $\widetilde{\tau}^{\left(q\right)}$ as
the value of $\tau^{\left(q\right)}$ under the infinite deceleration condition
(i.e., there is
$\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)\big{|}_{\widetilde{\tau}^{\left(q\right)}}=0$).
The optimal time scale for the SRG is set as

$\displaystyle\widehat{\tau}^{\left(q\right)}=\frac{2}{q+1}\widetilde{\tau}^{\left(q\right)},$
(24)

which is the average subdivided time scale for each $1$-simplex in a
$n$-simplex (i.e., the average number of adjacent $1$-simplexes between
arbitrary units $i$ and $j$ is $\frac{n+1}{2}$). Note that
$\widehat{\tau}^{\left(q\right)}=\widetilde{\tau}^{\left(q\right)}$ holds for
pairwise interactions.

The above procedure can be directly generalized to the non-ergodic case, where
we treat each cluster $C_{i}$ as a network and derive a
$\widetilde{\tau}^{\left(q\right)}$ for it if it has not been assigned with a
time scale yet. The time scale of a cluster does not change after it is
assigned. See Fig. 3(c) for illustrations.

After determining the time scale, the SRG proposed in Sec. IV.1 can be applied
to renormalize interacting systems in a non-parametric manner. In Appendix E,
we use the SRG to process a system whose pairwise interactions exhibit small-
world properties [80] and illustrate how the renormalization flow changes
across different combinations of $\left(p,q\right)$. This result qualitatively
suggests the effects of higher-order interactions on lower-order interactions
during renormalization, whose quantification is explored in our subsequent
analysis. In Fig. 4, we present more instances of the SRG applied on synthetic
interacting systems whose pairwise interactions follow the random tree (RT),
the square lattice (SL), the triangular lattice (TL), and the Watts-Strogatz
networks (WS). Please see Appendix A for the definitions of these random
networks.

###  IV.3 Scaling relation and high-order effects

Given the SRG, a crucial task is to analyze the potential scaling relation and
quantify the effects of high-order interactions on it.

In the ergodic case, we denote $\widetilde{X}^{\left(q\right)}$ and
$\widetilde{\rho}^{\left(q\right)}$ as the specific heat and density operator
associated with $\widetilde{\tau}^{\left(q\right)}$. We can insert these
variables into Eq. (23) to derive the following relation

$\displaystyle\frac{\widetilde{X}^{\left(q\right)}}{\widetilde{\tau}^{\left(q\right)}}=\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\widetilde{\rho}}+\mu,\;\forall\mu\in\mathbb{R},$
(25)

where
$\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\widetilde{\rho}}=\operatorname{tr}\left(\mathbf{L}^{\left(q\right)}_{1}\widetilde{\rho}^{\left(q\right)}\right)$
and constant $\mu$ is obtained when we solve Eq. (23) as a differential
equation. If the Laplacian eigenvalue spectrum (i.e., the density of the
eigenvalues of $\mathbf{L}^{\left(q\right)}_{1}$) follows a general power-law
form
$\operatorname{Prob}\left(\lambda^{\left(q\right)}_{1}\right)\sim\left(\lambda^{\left(q\right)}_{1}\right)^{\beta^{\left(q\right)}}$,
we can obtain a scaling relation

$\displaystyle\beta^{\left(q\right)}=\widetilde{X}^{\left(q\right)}-1-\nu,\;\exists\nu\in\left[0,\infty\right),$
(26)

whose derivations are presented in Appendix F. Parameter $\nu$ is a specific
constant. For a system with scale-invariance property on the $q$-order, the
power-law form of the Laplacian eigenvalue spectrum holds and should be
invariant under the transformation of the SRG. Therefore, we expect to see an
approximately constant $\widetilde{X}^{\left(q\right)}$ in a certain interval
of time scale (i.e., the infinite deceleration condition), accompanied by a
fixed $\beta^{\left(q\right)}$ during renormalization. For a system without
scale-invariance, the scaling relation vanishes.

Figure 5: Scaling relation of the simplex renormalization group. (a-e) We
verify the scaling relation on the systems whose pairwise interactions follow
a Barabási-Albert network (BA, $c=10$), a Watts-Strogatz network (WS, each
unit initially has $10$ neighbors and edges are rewired according to a
probability of $0.1$), and an Erdos-Renyi network (ER, two units share an edge
with a probability of $0.2$), the triangular lattice (TL, with the
interactions of the first two orders), and the random tree (RT, with only
pairwise interactions), respectively. Scatters denote the mean values of
$\left(\widetilde{X}^{\left(q\right)},\widetilde{\tau}^{\left(q\right)}\Big{\langle}\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\Big{\rangle}_{i}+\nu\right)$
averaged across all replicas on each order and error bars denote standard
deviations. See the definitions of random networks in Appendix A. (f) The
standard deviations compared across different systems. (g-i) The high-order
effects measured on the replicas in (a-c), where each color area denotes the
distribution of high-order effects in the $k$-iteration of the SRG
($k\in\\{1,2,3\\}$).

To generalize our analysis to the non-ergodic case, we respectively derive
$\widetilde{\tau}^{\left(q\right)}\left(i\right)$ and
$\beta^{\left(q\right)}\left(i\right)$ for each cluster $C_{i}$ of
$\mathbf{G}^{\left(q\right)}_{1}$ following Eqs. (25-26). These definitions
enable us to obtain

$\displaystyle\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\widetilde{\rho}}=\Big{\langle}\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\Big{\rangle}_{i}+\nu,\;\exists\nu\in\left[0,\infty\right),$
(27)

and

$\displaystyle\widetilde{X}^{\left(q\right)}=\widetilde{\tau}^{\left(q\right)}\Big{\langle}\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\Big{\rangle}_{i}+\nu,\;\exists\nu\in\left[0,\infty\right),$
(28)

where $\langle\cdot\rangle_{i}$ denotes weighted averaging across all clusters
of $\mathbf{G}^{\left(q\right)}_{1}$ (i.e., weighted according to cluster
size). In the non-ergodic case, variables $\widetilde{X}^{\left(q\right)}$ and
$\widetilde{\tau}^{\left(q\right)}$ refer to the global specific heat and the
time scale measured on $\mathbf{L}^{\left(q\right)}_{1}$ under the infinite
deceleration condition. Eqs. (25-26) actually serve as a special case of Eqs.
(27-28) where there exists only one cluster. Please see Appendix G for
derivations of these equations. Moreover, although the scaling relations in
Eqs. (25-28) are derived in a case with $k=1$, we can also explore these
scaling relations on $\mathbf{G}^{\left(q\right)}_{k}$ for $k>1$ (i.e., we
directly replace $\mathbf{G}^{\left(q\right)}_{1}$ with
$\mathbf{G}^{\left(q\right)}_{k}$ during calculating Eqs. (25-28)). In
summary, we can verify if a system obeys scaling relations and further explore
whether the SRG preserves scaling relations.

Eqs. (25-28) are derived for $q$-order interactions during renormalization and
may divide from the original scaling relation of $p$-order interactions. These
$p$-order interactions, initially corresponding to
$\widetilde{X}^{\left(p\right)}$ and $\beta^{\left(p\right)}$, are
renormalized according to $\widetilde{X}^{\left(q\right)}$ and
$\beta^{\left(q\right)}$ in the SRG. In the case with $p<q$, we can measure
the effects of $q$-order interactions on the scaling relation of $p$-order
ones as

$\displaystyle\zeta\left(p,q\right)=|\widetilde{X}^{\left(p\right)}-\widetilde{X}^{\left(q\right)}|,$
(29)

which is applicable to both ergodic and non-ergodic cases. Meanwhile, similar
to the scaling relation analysis, the measurement of high-order effects can be
implemented in any iteration of the SRG.

Overall, the proposed SRG offers an opportunity to verify the potential
existence of scale-invariance on different orders and classify interacting
systems according to their high-order scaling properties.

##  V Application of the simplex renormalization group

To this point, we have elaborated the SRG framework and its properties. Below,
we validate its applicability in analyzing the systems with high-order
interactions.

###  V.1 Verification of the multi-order scale-invariance

We aim at verifying the multi-order scale-invariance (i.e., being invariant
under the SRG transformation on different orders). Specifically, being scale-
invariant on an order requires the system to be located at a certain fixed
point of the renormalization flow (i.e., the system behaves similarly across
all scales), leading to scale-independent characteristics of the system
described by the power-law scaling. In other words, if we describe system
states by a set of observables, their measurements are expected to maintain
invariant across all the iterations of the SRG. Due to the internal complexity
of the system on different orders, we may see a rich variety of fixed points
in empirical data. In an opposite where the system lacks a fixed point, its
behaviours significantly change across scales.

In Appendix H, we implement our verification according to the behaviours of
the Laplacian eigenvalue spectrum (the behaviours of degree distributions are
also presented as auxiliary information). Based on the existence of invariant
power-law forms of Laplacian eigenvalue spectra under the SRG transformation,
we can classify systems into scale-invariant, weakly scale-invariant, and
scale-dependent types on different orders.

Apart from the above approach, there exists a more direct way to realize the
verification. Because a scale-invariant system on $q$-order should invariably
follow the scaling relation, being born to satisfy Eq. (27) is one of the
necessary conditions of being scale-invariant. Consequently, we can first
measure the deviations from the scaling relation of concerned systems before
renormalization. Given significantly large deviations, the system can not be
scale-invariant no matter how it behaves during renormalization. In Fig. 5,
the measurement is conducted on multiple types of systems in each $k$-th
iteration ($k\in\\{1,2,3\\}$) of the SRG (see Appendix A for random network
definitions and see Appendix I for experiment details). As shown in Figs.
5(a-e), although the mean observations of
$\left(\widetilde{X}^{\left(q\right)},\widetilde{\tau}^{\left(q\right)}\Big{\langle}\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\Big{\rangle}_{i}+\nu\right)$
averaged across all replicas collapse onto the scaling relation on every order
and in each $k$-th iteration, the departures from the scaling relation
measured by the standard deviations of $\widetilde{X}^{\left(q\right)}$ and
$\widetilde{\tau}^{\left(q\right)}\Big{\langle}\frac{\beta^{\left(q\right)}\left(i\right)+1}{\widetilde{\tau}^{\left(q\right)}\left(i\right)}\Big{\rangle}_{i}+\nu$
in some systems can be non-negligible. Specifically, as suggested in Fig.
5(f), the systems whose pairwise interactions follow known scale-invariant
structures, such as the triangular lattice or the random tree, do ensure that
all their replicas follow the scaling relation with vanishing deviations in
each $k$-th iteration. The replicas of the systems whose pairwise interactions
follow the Barabási-Albert network, a structure with weak scale-invariant
property (i.e., the Barabási-Albert network can be renormalized within a
specific set of scales but is not strictly scale-invariant [10]), have
relatively small deviations on the first three orders as well. The replicas of
other kinds of systems whose pairwise interactions follow the Watts-Strogatz
or the Erdos-Renyi networks [81] exhibit larger deviations from the scaling
relation on most orders, suggesting the absence of mutli-order scale-
invariance.

In sum, the method reported in Fig. 5 can offer a convenient detection of
scale-dependent systems in practice because scale-invariant systems are expect
to satisfy Eq. (27) before and during renormalization. To seek a more precise
verification, one should apply the approach in Appendix H.

Finally, we measure the high-order effects on the scaling relation applying
Eq. (29). As shown in Figs. 5(g-i), these high-order effects are non-
negligible and generally increase with the difference between $p$ and $q$.
These results explain how the SRG differs from conventional renormalization
groups that focus on a single order.

Figure 6: The topological invariance explored by the SRG. (a-d) The $m$-order
Betti number $\mathcal{B}_{m}$ ($m\in\\{0,\ldots,3\\}$) is calculated on
$\mathbf{G}^{\left(p\right)}_{k}$ in each $k$-th iteration of the SRG
($k\in\\{1,\ldots,5\\}$), where the multi-order Laplacian (MOL) is used. (e-h)
The same analysis is repeated using a SRG defined with the high-order path
Laplacian (HOPL). Each line in (a-h) corresponds to an instance of the
evolution of the $m$-order Betti number during renormalization. Please note
that we define the $y$-axis as $\mathcal{B}_{m}+1$ to meet the demands of a
log-scale axis.

###  V.2 Exploration of the topological invariance

In Appendix B, we have explained how the SRG is theoretically related to the
persistent homology theory and other concepts, where we suggest the potential
applicability of the SRG in analyzing persistent homology. Below, we validate
this applicability by using the SRG to explore the topological invariance.

In brief, the persistent homology theory can be used to study the lifetimes of
different topological properties in a filtration process, during which we
progressively add or remove simplices from the system [62, 74, 75, 76] (i.e.,
this is similar to the network evolution process in complex network studies
[82]). Although the SRG is not exactly same as the classic filtration approach
because it modifies the system in a more complicated manner, it does offer an
opportunity to study the evolution of topological properties during
renormalization (i.e., we can measure different topological properties on
$\mathbf{G}^{\left(p\right)}_{k}$ for each $k\geq 1$). Similar to the analysis
of scale-invariance, a topological property can be treated as invariant if it
is persistent across all the iterations of the SRG (i.e., referred to as the
topological invariance). In our analysis, we mainly focus on the Betti number
on different orders, a family of important topological properties used in
biology [83, 84, 85] and physics [86, 87, 88] studies. In general, a $m$-order
Betti number, $\mathcal{B}_{m}$, counts the number of $m$-dimensional holes in
the system, which is a key topological characteristic of high-order structures
and defines the dimension of $m$-th homology [74]. In the terminology of graph
theory, a $0$-order Betti number counts the number of connected components and
a $1$-order Betti number is the number of cycles. Computationally, the
$m$-order Betti number can be calculated using the toolbox offered by Ref.
[89].

In Fig. 6, we explore the potential existence of topological invariance in
different random network models, including the Barabási-Albert network
($c=4$), the Watts-Strogatz network (each unit initially has 10 neighbors and
edges are rewired according to a probability of 0.3), the Erdos-Renyi network
(two units share an edge with a probability of 0.1), the triangular lattice,
and the random tree. One can see Appendix A for their definitions. We generate
$100$ replicas for each random network model, where every replica consists of
$200$ units. In the experiment, we design a SRG with a certain combination of
$\left(p,q\right)$ and let it function on each random network to generate a
sequence of $\mathbf{G}^{\left(p\right)}_{k}$ for $k\in\\{1,\ldots,5\\}$.
Given each $\mathbf{G}^{\left(p\right)}_{k}$, we measure the $m$-order Betti
number, $\mathcal{B}_{m}$, on it ($m\in\\{0,\ldots,3\\}$). If a $m$-order
Betti number is generally constant across all iterations, its corresponding
topological property can be treated as invariant during renormalization.

As shown in Fig. 6, the Betti numbers of scale-dependent systems (e.g., the
Watts-Strogatz and the Erdos-Renyi networks) exhibit non-negligible changes on
most orders when the SRG transforms systems across scales. Compared with these
systems, weakly scale-invariant systems (e.g., the Barabási-Albert network)
have less variable Betti numbers during renormalization. When systems become
scale-invariant (e.g., the triangular lattice and the random tree), the Betti
numbers on all orders are generally invariant except that some fluctuations
may occur in certain replicas due to random effects. This finding is
interesting because, to our best knowledge, there is no theoretical guarantee
that the multi-order scale-invariance is correlated with the topological
invariance yet. According to our observation, we have the following
speculation: there may exist specific topological invariance properties near
the critical point that make the SRG, a framework dealing with high-order
topological structures, have less effects on the information contained in the
system. As a result, the behaviours of the system seem to be approximately
invariant under the scale transformation of the SRG, leading to the inference
that the system is born at a fixed point of the renormalization flow.
Exploring this speculation may deepen our understanding of the relation
between the persistent homology [62, 74, 75, 76] and the criticality studied
by the renormalization group theory [6, 7], which remains as a valuable
direction of future works.

###  V.3 Identification of organizational structure

Next, we use the SRG to explore the roles of different order of interactions
in preserving or defining the organizational structures (e.g., latent
communities) of complex systems, where we also validate the applicability of
the SRG in organizational structure identification.

Figure 7: The organizational structure discovered by renormalization flows.
(a-c) The Louvain communities underlying the Genetics (Ge), WI2007 (WI), and
Inflamome (In) data sets [90] are extracted using the Louvain community
detection algorithm [91]. The units of different communities are distinguished
by colors. For comparison, the SRG is applied to renormalize these systems
with $\left(p,q\right)=\left(1,2\right)$ and
$\left(p,q\right)=\left(1,1\right)$. Units are determined as belonging to the
same SRG community if they are aggregated into the same super-unit. Meanwhile,
the corrected SRG communities are presented. (d-f) The correctness of
communities and disconnectivity between communities are shown. (g-i) The
adjusted mutual information, completeness, and homogeneity are measured on SRG
communities. Notion “Ge (1,1)” denotes the SRG communities with
$\left(p,q\right)=\left(1,1\right)$ on the Genetics data set. Notion “Ge Co”
denotes the corrected SRG communities on the Genetics data set. Other notions
can be understood in a similar way. Please see Appendix A for the definitions
of all measures used in (d-i).

When organizational structures are given, we explore the conditions under
which these structures can be maintained by the SRG. In Appendix J, we
implement the verification using protein-protein interactions (e.g.,
orthologous and genetic interactions) in Caenorhabditis elegans [90]. We first
extract the initial community structures of protein-protein interactions by
applying the Louvain community detection algorithm [91] (see Appendix A for
introduction). Then, we use the SRG to renormalize protein-protein
interactions on a high order (i.e., with $\left(p,q\right)=\left(1,2\right)$)
or a low order (i.e., with $\left(p,q\right)=\left(1,1\right)$). We define
latent community structures according to unit aggregation during
renormalization (see Appendix J for details). Consistent with Ref. [38], our
results in Appendix J highlight that renormalization can never be treated as a
trivial counterpart of clustering. Meanwhile, the SRG guided by high-order
interactions (e.g., $q=2$) can generally preserve the properties of initial
community structures while the SRG guided by pairwise interactions does not.
Based on these results, we speculate that high-order interactions are
essential in forming and characterizing organizational structures. Communities
can not be treated as the trivial collections of pairwise interactions and it
may be inappropriate to partition a system into sub-systems only based on
pairwise relations. Given this observation, we naturally wonder if high-order
interactions alone are sufficient to define ideal community structures (i.e.,
do not consider pairwise relations while determining system partition).

To tackle this question, we explore whether a SRG defined with different
combination of $\left(p,q\right)$ can be applied to identify latent
communities when system partition is unknown. In Fig. 7(a-c), we define
communities according to the renormalization flows of the SRG and compare them
with the Louvain communities (i.e., the communities detected by the Louvain
algorithm [91], see Appendix A). Specifically, we classify the units into the
same SRG community if they are aggregated into the same super-unit during
renormalization (e.g., we generate SRG communities after the first iteration
of renormalization flows in Fig. 7(a-c)).

To quantitatively compare the SRG communities with the Louvain communities, we
analyze the correctness of communities (i.e., or referred to as the
performance) [92] and the disconnectivity between communities (i.e., or
referred to as the coverage) [92]. Note that both these concepts reflect the
properties of communities on $1$-order because they are defined on pairwise
interactions [92] (see Appendix A for precise definitions). As shown in Fig.
7(d-f), the behaviours of SRG are different from the Louvain community
detection algorithm, which is designed to maximize the modularity of
communities [91, 93, 94]. When the SRG is guided by pairwise interactions, the
generated SRG communities tend to be disconnected with each other (i.e., lack
pairwise interactions across communities) and have low accuracy (i.e., units
in the same community may lack pairwise interactions because lots of short-
range pairwise interactions are reduced). When the SRG is guided by high-order
interactions, the generated SRG communities tend to be precise (i.e., units in
the same community always involve in high-order interactions) and non-isolated
(i.e., there exist across-community pairwise interactions since the SRG mainly
reduces high-order interactions). We are inspired to generate the SRG
communities by combining the properties of these two cases. We suggest
considering pairwise interactions to correct the SRG communities generated by
the SRG guided by high-order interactions. Specifically, for two SRG
communities generated with $\left(p,q\right)=\left(1,2\right)$, if at least
half of the units in one community have pairwise interactions with the units
in another community, then all units in these two SRG communities are merged
into the same community. The corrected SRG communities are presented in Fig.
7(d-f), exhibiting better capability in balancing between the correctness and
the disconnectivity.

From another perspective, we quantify the consistency between the SRG
communities and the Louvain communities in defining system partition in Fig.
7(g-i). As measured by the adjusted mutual information, completeness, and
homogeneity (see Appendix A for the definitions of these metrics, whose larger
values suggest higher consistency extents), the partitions defined by the
corrected SRG communities are generally consistent with, but not exactly
equivalent to, those formed by the Louvain communities in all data sets. On
the other hand, the SRG communities defined only based on high-order (i.e.,
with $\left(p,q\right)=\left(1,2\right)$) or low-order (i.e., with
$\left(p,q\right)=\left(1,1\right)$) interactions do not necessarily keep
consistency with the Louvain communities in these data sets. Because the
Louvain communities are defined by maximizing the modularity [91, 93, 94], we
suggest that the optimization of modular structures can not be realized only
based on pairwise or high-order interactions. Instead, modular structures
depend on multiple orders of interactions.

Taken together, high-order interactions are more effective than pairwise
interactions in preserving the principal characteristics of organizational
structures. However, merely high-order interactions alone are not sufficient
to define ideal system partitions with optimal modular structures and clear
separability among communities. For real systems, their organizational
structures are defined by the intricate coexistence of various orders of
interactions, exhibiting diverse behaviours across different orders. Although
the SRG does not behave in a manner equivalent to the optimized algorithms for
community detection (e.g., the Louvain algorithm [91]), it offers a flexible
way to analyze the roles of the interactions of every order in characterizing
organizational structures since we can freely select different combinations of
$\left(p,q\right)$ in analysis. Moreover, after considering the interactions
of multiple orders, the corrected SRG communities have the potential to offer
a system partition scheme competitive to optimized algorithms (see Fig.
7(d-i)), which suggesting the applicability of the SRG in identifying latent
community structures of real systems.

###  V.4 Optimization of information bottleneck

Finally, we analyze the SRG in terms of informational properties. In previous
studies, the analysis using information theory has suggested that
renormalization groups essentially maximize the mutual information between
relevant features and the environment [95, 96] or reduce the mutual
information among irrelevant features [97]. In this work, we suggest exploring
the SRG from the perspective of information bottleneck [98] because
renormalization groups, similar to dimensionality reduction or representation
learning approaches in machine learning [99, 100, 101], essentially deal with
information encoding during data compression [102].

Figure 8: The information curve of the SRG. (a) The SRG is applied on four
synthetic interacting systems, whose pairwise interactions follow the Erdos-
Renyi network (ER, each pair of units share an edge with a probability of
$0.01$), the Watts-Strogatz network (WS, each unit initially has $5$ neighbors
and edges are rewired according to a probability of $0.3$), the Barabási-
Albert network (BA, $c=1$), and the random tree (RT), respectively. Each kind
of system consists of $500$ units and has $10$ replicas. The SRG is defined
with the multi-order Laplacian and $\left(p,q\right)=\left(1,1\right)$ for RT
while
$\left(p,q\right)\in\\{\left(1,1\right),\left(1,2\right),\left(1,3\right)\\}$
is set for all other systems. The generated information flows are used to
derive information curves. (b) The complexity per preserved information bit is
shown after averaging across all replicas.

To realize this analysis, we follow one of our earlier works [103] to
represent the associated high-order network sketch in the $k$-th iteration of
the renormalization flow, $\mathbf{G}^{\left(p\right)}_{k}$, by a Gaussian
variable

$\displaystyle\mathbf{X}^{\left(p\right)}_{k}\sim\mathcal{N}\left(\mathbf{0},\mathbf{L}^{\left(p\right)}_{k}+\frac{1}{N^{\left(p\right)}_{k}}\mathbf{J}\right),$
(30)

where $\mathbf{J}$ is an all-one matrix and $N^{\left(p\right)}_{k}$ measures
the number of units in $\mathbf{G}^{\left(p\right)}_{k}$. The derived Gaussian
variable offers a optimal representation of $\mathbf{G}^{\left(p\right)}_{k}$
with network-topology-dependent smoothness and maximum entropy properties and
enable us to compare $\mathbf{G}^{\left(p\right)}_{k}$ across different
iterations (see Ref. [103] for explanations). If necessary, one can further
add a scaled unitary matrix, $a\mathbf{I}$ ($a>0$), to the covariance matrix
of $\mathbf{X}^{\left(p\right)}_{k}$ to ensure its semi-positive property in
the non-ergodic case (i.e., the case where $\mathbf{G}^{\left(p\right)}_{k}$
is disconnected).

We suggest to consider the following information bottleneck

$\displaystyle\mathcal{L}_{IB}=\underbrace{\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)}_{\text{Preserved
information}}-\psi\underbrace{\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i}\right)}_{\text{Complexity}},\;\forall\psi>0$
(31)

where $i=k-1$ if $k>1$ and $i=k$ if $k=1$. Parameter $\psi$ denotes the
regularization strength. Notion $\mathcal{I}\left(\cdot;\cdot\right)$ denotes
the mutual information (see Appendix A for its definition), which can be
estimated by the mutual information estimator with local non-uniformity
corrections (we set the correction strength as $10^{-3}$ such that corrections
only happen when necessary) [104]. This approach has been included in the non-
parametric entropy estimation toolbox [105]. Note that Eq. (31) is different
from the information bottleneck considered in Ref. [102] due to the
discrepancy between our concerned question and Ref. [102].

Intuitively, the objective function in Eq. (31) evaluates whether the
renormalized system in the $k$-th iteration,
$\mathbf{G}^{\left(p\right)}_{k}$, preserves the information of the original
system, $\mathbf{G}^{\left(p\right)}_{1}$, by maximizing
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)$
and reduces the trivial dependence on its previous state,
$\mathbf{G}^{\left(p\right)}_{i}$, by minimizing
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i}\right)$.
For a valid renormalization group, being able to preserve some information of
the original system is an essential capability. Because a trivial solution of
maximizing
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)$
is to make $\mathbf{G}^{\left(p\right)}_{k}=\mathbf{G}^{\left(p\right)}_{k-1}$
for any $k>1$ (i.e., there is no renormalization at all), we need to include
the complexity term in Eq. (31) to avoid this trivial result.

Meanwhile, the possibility of maximizing
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)$
depends on the properties of the system as well. If the system is scale-
invariant, the preserved information does not rapidly vanish when we reduce
the complexity since the information of $\mathbf{G}^{\left(p\right)}_{1}$
holds across scales (note that
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)$
still reduces slightly because the decreasing dimensionality of
$\mathbf{X}^{\left(p\right)}_{k}$ affects the numerical behaviours of the
estimator). If the system is scale-dependent, the preserved information
crucially relies on the complexity and vanishes unless the complexity is high.
Therefore, we can divide
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{1}\right)$
by
$\mathcal{I}\left(\mathbf{X}^{\left(p\right)}_{k};\mathbf{X}^{\left(p\right)}_{i}\right)$
to measure the complexity required by a single bit of preserved information in
the $k$-th iteration. A high value of complexity per information bit suggests
that the system is scale-dependent because we can not find a low-complexity
coarse-grained system to represent it efficiently.

In Fig. 8(a), we show the information curves of the renormalization flows
generated on four kinds of interacting systems. As suggested by our results,
the renormalization flows of weakly scale-invariant (e.g., BA) and scale-
invariant (e.g., RT) systems lead to sufficiently small complexity values and
numerous preserved information bits hile the renormalization flows of scale-
dependent systems (e.g., ER and WS) do not always achieve so (although WS
realizes complexity reductions during renormalization, it is still not
competitive with BA and RT). The information bits of scale-dependent systems
are preserved at the cost of maintaining high complexity while scale-invariant
and weakly scale-invariant systems progressively become more complexity-saving
in preserving information during renormalization. These results can be
quantitatively validated based on the complexity per information bit measured
in Fig. 8(b). These findings suggest the possibility of verifying scale-
invariance from an informational perspective. Moreover, we observe that the
SRG guided by pairwise interactions (i.e., $q=1$) frequently achieve more
significant complexity reductions while preserving information. This
phenomenon may arise from the fact that high-order interactions are
distributed in systems in a sparser manner.

##  VI Discussions

Various real interacting systems share universal characteristics that govern
their dynamics and phase transition behaviours [106]. These characteristics
remain elusive because fundamental physics tools, such as path integrals and
renormalization group theories, are not well-established when _a priori_
knowledge about system mechanisms is absent [10]. This gap may be accountable
for diverse controversies concerning the critical phenomena in biological or
social systems (e.g., the controversies about brain criticality [107]) and has
induced a booming field that is devoted to developing statistical physics
approaches for the empirical data generated by unknown mechanisms [10].
Notable progress is accomplished in the computational implementations of path
integrals and renormalization groups on real systems with pairwise
interactions [29, 30, 10, 31, 32, 33], supporting remarkable applications in
analyzing neural dynamics [108] and swarm behaviours [109]. Nevertheless,
there exists no apparent equivalence between these pioneering works and an
appropriate approach for analyzing un-decomposable high-order interactions
because classic network or geometric representations are invalid in
characterizing polyadic relations [11, 12, 13, 9, 14].

In this work, we contribute to the field by developing the natural
generalizations of path integrals and renormalization groups on un-
decomposable high-order interactions. Our main contributions are summarized
below:

  * (1)

Our theoretical derivations lead us to an intriguing perspective that the
system evolution governed by high-order interactions can be fully formalized
by path integrals on simplicial complexes (e.g., the return probability of an
arbitrary system state is proportional to the path integral of all closed
curves that start from and end with this state in high-order networks), which
suggests the possibility for studying high-order interactions applying the
tools of quantum field theory [3, 4].

  * (2)

We suggest that the contributions of microscopic fluctuations to macroscopic
states in the moment space are precisely characterized by a function of the
multi-order Laplacian [9] or our propsed high-order path Laplacian.
Consequently, a renormalization group, the SRG, can be directly developed
based on the diffusion on simplicial complexes, which does not require any
assumption on latent metric for mapping units into target spaces (e.g., the
hyperbolic space [30, 33]). This property ensures the general applicability of
our theory on real complex systems, where _a priori_ knowledge is rare and any
assumption can be unreliable. Different from classic renormalization groups,
the SRG can renormalize $p$-order interactions in a system based on the
structure and dynamics associated with $q$-order interactions ($p\leq q$),
enabling us to analyze the effects of $q$-order interactions on $p$-order ones
and study the scaling properties across different orders.

  * (3)

We propose a divide-and-conquer approach to deal with non-ergodic cases, where
some systems states are never reachable by the evolution started from other
systems states (i.e., without global connectivity). The developed framework
improves the validity of our theory for high-order interactions, which usually
feature sparse distributions that break system ergodicity on a high order.

  * (4)

We seek for a comprehensive analysis of the scaling relation in both ergodic
and non-ergodic cases. Our theoretical derivations and computational
experiments suggest that the scaling relation can help differentiate among
scale-invariant, weakly scale-invariant, and scale-dependent systems on an
arbitrary order (consistent with verification using the behaviours of the
Laplacian eigenvalue spectrum). Meanwhile, the effects of high-order
interactions on the scaling relation can be measured as well.

  * (5)

We propose a way to use the SRG to verify the existence of topological
invariance, suggesting the possibility to relate the SRG with the persistent
homology analysis [62, 74, 75, 76]. As empirically demonstrated in our
results, the SRG can serve as a tool similar to the filtration process [62,
74, 75, 76] during studying the lifetimes of different topological properties.
In our experiments, we observe a correlated relation between the multi-order
scale-invariance and the topological invariance (i.e., multi-order scale-
invariant systems are more likely to satisfy the topological invariance),
which suggests the potential role of topology in shaping the statistical
physics properties of complex systems.

  * (6)

Apart from validating the essential difference between the SRG and clustering
[38], our experiments also suggest the applicability of the SRG in identifying
organizational structures (e.g., latent communities) of real complex systems
and analyzing the crucial roles of the interactions of different orders in
characterizing these structures.

  * (7)

We suggest that the SRG can be analyzed from the perspective of information
bottleneck, which quantifies how the renormalized system preserves the
information of its original properties while reducing complexity. We discover
that the complexity required for maintaining one information bit can
distinguish between scale-invariant and scale-dependent systems in an
informational aspect.

In sum, by extending classic path integrals and renormalization groups to
simplicial complexes, our research reveals a novel route to studying the
universality classes of complex systems with intertwined high-order
interactions. The proposed simplex path integral and simplex renormalization
group can serve as precise tools for characterizing system dynamics,
discovering intrinsic scales, and verifying potential scale invariance on
different orders. The revealed information by our theory may shed light on the
intricate effects of the interplay among multi-order interactions on system
dynamics properties (e.g., phase transitions or scale invariance). To ensure
the capacity of our theory in analyzing the real data sets that are governed
by unknown mechanisms or lack clear network structures, we suggest to consider
_a-priori_ -knowledge-free framework in future studies. This framework begins
with applying specific non-negative non-parametric metrics (e.g., distance
correlation [110] or co-mutual information [111]) to evaluate the coherence
between units and define the adjacency matrix of pairwise interactions. Then,
the high-order representation and the simplex renormalization group can be
progressively calculated following our theory.

## Data availability statement

Any data that support the findings of this study are included within the
article.

## Acknowledgements

This project is supported by the Artificial and General Intelligence Research
Program of Guo Qiang Research Institute at Tsinghua University (2020GQG1017)
as well as the Tsinghua University Initiative Scientific Research Program.
Authors appreciate Hedong Hou, who studies at the Institut de Mathématiques
d’Orsay, for his inspiring discussions.

##  Appendix A Auxiliary concept definitions

There are numerous auxiliary concepts used to support our analysis. In this
section, we elaborate their basic definitions and offer necessary references.

#### Simplex and simplicial complex.

For clarity, we first recall the elementary definitions of simplex and
simplicial complex [60, 61]. A $q$-simplex, as visualized in Fig. 1, is a
convex hull of $q+1$ affinely independent points in $\mathbb{R}^{h}$ ($h\geq
q$). A simplicial complex, $K$, is a collection of simplices in
$\mathbb{R}^{h}$ (see Fig. 1) that satisfies the following properties:

  * (1)

Each face $\sigma^{\left(p\right)}$ of simplex $\sigma^{\left(q\right)}\in K$
also belongs to $K$.

  * (2)

The non-empty intersection of each pair simplices, $\sigma^{\left(p\right)}$
and $\sigma^{\left(q\right)}$, in $K$ must be a face of both simplices
$\sigma^{\left(p\right)}$ and $\sigma^{\left(q\right)}$ (here $p$ and $q$ are
not necessarily the same).

The simplicial complex considered in our work, as shown in Fig. 1, is usually
referred to as the clique complex in the persistent homology analysis [62,
60]. In general, a $q$-simplex is a clique of $q+1$ units (i.e., a $q$-order
interaction requires the simultaneous participation of $q+1$ units).

#### Random network models.

In our computational experiments, we use the following random networks:

  * (1)

The Erdos-Renyi network (ER) [81] is an elementary type of random network in
which two units share an edge with a probability of $p$ or are disconnected
with a probability of $1-p$. We have used the ER networks in Fig. 5, Fig. 8,
and Figs. 10-11.

  * (2)

The Watts-Strogatz network (WS) [81] is a typical model of small-world
phenomenon. It is generated by an edge rewiring process, during which we first
initialize a regular network (i.e., a network where all units have the same
number of neighbors) and rewire each edge according to a probability of $p$
(i.e., the edge remains unchanged with a probability of $1-p$). Several WS
networks are implemented in Figs. 4-5, Fig. 8, and Figs. 9-11.

  * (3)

The Barabási-Albert network (BA) [47, 48] is a representative example of
weakly scale-invariant system [10]. It is generated by the network growth
process (i.e., a process for progressively adding new units and edges into a
null network), where the number of random edges to attach from a new unit to
existing units is denoted as $c$. One can see the BA networks in Figs. 1-3,
Fig. 5, Fig. 8, and Figs. 10-11.

  * (4)

The random tree (RT) [112] is an instance of scale-invariance, which is
generated from the uniformly random Prüfer sequence in computational
implementation. The instances of RT are implemented in Figs. 4-5, Fig. 8, and
Figs. 10-11.

  * (5)

The triangular lattice (TL) [113] is another instance of scale-invariance when
system size approaches to infinity. Under the periodic boundary condition, a
finite triangular lattice system is also a regular network where each unit has
six neighbors. When the periodic boundary condition is absent, only the units
within boundary form a regular network while the units on boundary may exhibit
different behaviours under finite size effects. One can see the instances of
TL in Figs. 4-5 and Figs. 10-11.

  * (6)

The square lattice (SL) [114] is also an instance of scale-invariance when
system size approaches to infinity. It shares similar properties with the
triangular lattice except that each unit in a finite square lattice system has
four rather than six neighbors under the the periodic boundary condition. The
SL has been analyzed in Fig. 4.

Among these random networks, the ER, WS, and BA networks may contain arbitrary
orders of interactions, depending on the settings of network parameters. The
RT only contains the interactions of the first order due to its tree
structure. The TL only includes the interactions of the first two orders
(i.e., we can not find any tetrahedron or any other simplicial complex of a
higher order in the TL system). The SL is formed by the interactions of the
first order (i.e., there is no triangle or any other simplicial complex of a
higher order in the SL system).

#### Organizational structure identification, preserving, and evaluation.

In our experiments, we tackle the identification and preserving tasks of
organizational structures of interacting systems, where we compare our SRG
framework with the Louvain community detection algorithm [91, 93, 94]. The
Louvain community detection algorithm searches for an ideal partition by
maximizing the modularity of communities, which is essentially different from
renormalization groups because it is purely optimization-based [91, 93, 94].
For convenience, we refer to the system partition identified by the Louvain
community detection algorithm [91, 93, 94] as the Louvain communities.

When an organizational structure has been given (e.g., if the Louvain
communities have been derived), we can explore whether the SRG can preserve
this organizational structure during renormalization. The corresponding
experiment is shown in Fig. 12 and Appendix J. We do not repeatedly elaborate
the experiment in this section.

If the organizational structure is unknown, we can explore if the SRG can
identify an appropriate structure. To evaluate the effectiveness of
organizational structure identification, we consider the following metrics:

  * (1)

The correctness of communities (i.e., or referred to as the performance)
measures the proportion of correctly interpreted pairs of units among all
pairs of units, i.e., two units within the same community are expected to
share an edge or two units belonging to different communities are expected to
be disconnected [92]. An ideal organizational structure should make all pairs
of units correctly interpreted, implying a correctness score of $1$.

  * (2)

The disconnectivity between communities (i.e., or referred to as the coverage)
measures the proportion of intra-community edges among all edges [92]. An
ideal organizational structure is expected to make all edges fall within
communities (i.e., all communities are perfectly separable from each other),
leading to a disconnectivity score of $1$.

One can see the corresponding experiment in Fig. 7, where the Louvain
community detection algorithm [91, 93, 94] is implemented for comparison.
Apart from the identification effectiveness mentioned above, the quantitative
comparison between the SRG communities and the Louvain communities can be also
realized by the following metrics:

  * (1)

The adjusted mutual information [115] is a mutual information score (i.e.,
similarity degree) between the SRG and the Louvain communities, which is
regulated by chance and is independent of permutation [115]. The adjusted
mutual information approaches to $1$ when the compared communities are
identical.

  * (2)

The completeness score [116] measures the consistency between the SRG and the
Louvain communities, which approaches to $1$ if all units within each SRG
community also belong to the same Louvain community.

  * (3)

The homogeneity score [116], similar to the completeness score, measures the
consistency between the SRG and the Louvain communities, which approaches to
$1$ if all units within each Louvain community also belong to the same SRG
community.

#### Information quantity measurement.

In our analysis, we use the mutual information to quantify the information of
one random variable, $A$, successfully encoded by another random variable,
$B$. Mathematically, the mutual information is defined as a Kullback-Leibler
divergence
$\mathcal{I}\left(A;B\right)=\mathcal{D}_{\text{KL}}\left[\text{Prob}\left(A,B\right)\|\text{Prob}\left(A\right)\text{Prob}\left(B\right)\right]$
[104]. Computationally, the mutual information is estimated by the mutual
information estimator with local non-uniformity corrections [104], which is an
approach included in the non-parametric entropy estimation toolbox [105] and
proven as efficient in dealing with high-dimensional variables. The
corresponding results are used in the information bottleneck analysis shown in
Fig. 8.

##  Appendix B Theoretical relations with the persistent homology theory and
related concepts

There are two closely related mathematical questions concerning the SRG:

  * (1)

What are the relations among the multi-order Laplacian, the high-order path
Laplacian, the combinatorial Laplacian [60], and the Hodge Laplacian [61, 23]?

  * (2)

How is the SRG related to the persistent homology theory [62]?

Although these questions are not the essential parts of our theory, discussing
them may help understand the SRG from other perspectives. Below, we sketch our
main ideas and offer necessary references.

#### Similarities and differences among four kinds of Laplacian operators.

Before we summarize key similarities and differences among four kinds of
Laplacian operators, we first present the basic definitions of the
combinatorial Laplacian [60] and the Hodge Laplacian [61, 23].

Generalized from the graph theory, a $q$-simplex also has its own adjacency
relations within a simplicial complex, whose number is referred to as the
degree of this $q$-simplex. The non-trivial part lies in that a $q$-simplex
can be adjacent to both $\left(q-1\right)$-simplices (i.e., lower adjacency)
and $\left(q+1\right)$-simplices (i.e., upper adjacency). Therefore, the
degree of a $q$-simplex is the sum of its lower and upper degrees (i.e., the
number of adjacent simplices of lower and higher orders) [60, 61]

$\displaystyle\operatorname{deg}\left(\sigma^{\left(q\right)}\right)$
$\displaystyle=\operatorname{deg}_{L}\left(\sigma^{\left(q\right)}\right)+\operatorname{deg}_{U}\left(\sigma^{\left(q\right)}\right),$
(32)
$\displaystyle=q+1+\operatorname{deg}_{U}\left(\sigma^{\left(q\right)}\right).$
(33)

Here we can derive Eq. (33) because
$\operatorname{deg}_{L}\left(\sigma^{\left(q\right)}\right)$ measures the
number of non-empty $\left(q-1\right)$-simplices in $K$ that are the faces of
$\sigma^{\left(q\right)}$, which always equals $q+1$.

To construct the combinatorial Laplacian [60] and the Hodge Laplacian [61,
23], we need to define orderings on the units of simplex and simplicial
complex [60, 61, 23]. Let us consider a $q$-simplex
$\sigma^{\left(q\right)}=\left[v_{1},\ldots,v_{q+1}\right]$, an ordering of
units $v_{1},\ldots,v_{q+1}$ defines an orientation of
$\sigma^{\left(q\right)}$. We treat two orderings on
$\left[v_{1},\ldots,v_{q+1}\right]$ are equivalent if and only if they differ
from each other by an even permutation. Simplex $\sigma^{\left(q\right)}$ is
oriented if an orientation is given. Apart from the orientation, another
concept to define is the boundary operator, which is closely related to the
chain complex in topology [60, 61, 23]. Let us consider the $q$-chain, a
formal weighted sum of $q$-simplices in $K$ whose weight coefficients are
selected from field $\mathbb{Z}_{2}$. A $q$-chain with a zero boundary is
referred to as the $q$-cycle in topology. Given the addition operation in
field $\mathbb{Z}_{2}$, the set of all $q$-chains form a chain group,
$C^{\left(q\right)}\left(K\right)$. The boundary operator is a mapping between
different chain groups,
$\partial^{\left(q\right)}:C^{\left(q\right)}\left(K\right)\rightarrow
C^{\left(q-1\right)}\left(K\right)$. Specifically, it maps a $q$-chain, a
weighted sum of $q$-simplices, to a new sum of $\left(q-1\right)$-simplices
with the same weights

$\displaystyle\partial^{\left(q\right)}\sigma^{\left(q\right)}=\sum_{i=1}^{q+1}\left(-1\right)^{i}f\left(\sigma^{\left(q\right)},i\right),$
(34)

where $f\left(\cdot,\cdot\right)$ denotes a removing operation that creates a
$\left(q-1\right)$-simplex from $\sigma^{\left(q\right)}$. Specifically, we
let
$f\left(\sigma^{\left(q\right)},i\right)=\left[v_{1},\ldots,v_{i-1},v_{i+1},\ldots,v_{q+1}\right]$
be a $\left(q-1\right)$-simplex generated from $\sigma^{\left(q\right)}$ by
omitting unit $v_{i}$. The chain complex mentioned above is a series of chain
groups transmitted by boundary operators [60, 61, 23]

$\displaystyle\cdots\xrightarrow{\partial^{\left(q+2\right)}}C^{\left(q+1\right)}\left(K\right)\xrightarrow{\partial^{\left(q+1\right)}}C^{\left(q\right)}\left(K\right)\xrightarrow{\partial^{\left(q\right)}}C^{\left(q-1\right)}\left(K\right)\xrightarrow{\partial^{\left(q-1\right)}}\cdots$
(35)

which plays an important role in the persistent homology theory [62]. Given a
$q$-boundary operator $\partial^{\left(q\right)}$, we can also define the
corresponding $q$-adjoint boundary operator,
$\partial^{\left(q,*\right)}:C^{\left(q-1\right)}\left(K\right)\rightarrow
C^{\left(q\right)}\left(K\right)$. The above definitions can be directly
applied on simplex, simplicial complex, and their oriented versions. Based on
these definitions, the combinatorial Laplacian [60] is given as

$\displaystyle\Delta^{\left(q\right)}_{C}=\partial^{\left(q+1\right)}\partial^{\left(q+1,*\right)}+\partial^{\left(q,*\right)}\partial^{\left(q\right)}.$
(36)

The matrix representation of the combinatorial Laplacian,
$\mathbf{L}^{\left(q\right)}_{C}$, when $q>0$ is [60]

$\displaystyle\left[\mathbf{L}^{\left(q\right)}_{C}\right]_{ij}=\begin{cases}\operatorname{deg}\left(\sigma^{\left(q\right)}_{i}\right),&i=j\\\
1,&i\neq
j\;\;\text{and}\;\;\sigma^{\left(q\right)}_{i}\overset{U}{\nsim}\sigma^{\left(q\right)}_{j}\;\;\text{and}\;\;\sigma^{\left(q\right)}_{i}\overset{L}{\sim}\sigma^{\left(q\right)}_{j}\;\;\text{with
similar orientation}\\\ -1,&i\neq
j\;\;\text{and}\;\;\sigma^{\left(q\right)}_{i}\overset{U}{\nsim}\sigma^{\left(q\right)}_{j}\;\;\text{and}\;\;\sigma^{\left(q\right)}_{i}\overset{L}{\sim}\sigma^{\left(q\right)}_{j}\;\;\text{with
dissimilar orientations}\\\ 0,&i\neq
j\;\;\text{either}\;\;\sigma^{\left(q\right)}_{i}\overset{U}{\sim}\sigma^{\left(q\right)}_{j}\;\;\text{or}\;\;\sigma^{\left(q\right)}_{i}\overset{L}{\nsim}\sigma^{\left(q\right)}_{j}\end{cases}.$
(37)

When $q=0$, the combinatorial Laplacian matrix reduces to the classic
Laplacian matrix [117]

$\displaystyle\left[\mathbf{L}^{\left(q\right)}_{C}\right]_{ij}=\begin{cases}\operatorname{deg}\left(\sigma^{\left(q\right)}_{i}\right),&i=j\\\
1,&\sigma^{\left(q\right)}_{i}\overset{U}{\sim}\sigma^{\left(q\right)}_{j}\\\
-1,&\text{otherwise}\end{cases}.$ (38)

In Eqs. (37-38), notion $\overset{U}{\sim}$ denotes the upper adjacent
relation, i.e., two $q$-simplices are incident to the same
$\left(q+1\right)$-simplex, and $\overset{L}{\sim}$ stands for the lower
adjacent relation, i.e., two $q$-simplices are incident to the same
$\left(q-1\right)$-simplex. If two $q$-simplices are upper adjacent (e.g.,
incident to simplex $\sigma^{\left(q+1\right)}$), we say they have similar
orientations if they share the same sign in
$\partial^{\left(q+1\right)}\sigma^{\left(q+1\right)}$ and they have
dissimilar orientations if their signs are opposite.

The derivation of the Hodge Laplacian requires more knowledge about the Rham-
Hodge theory [61, 23], which will not be elaborated here due to its
theoretical complexity. We refer to Ref. [61] for more details. In physics and
applied mathematics, the discretized Hodge Laplacian is more frequent to see,
which is defined as [60]

$\displaystyle\mathbf{L}^{\left(q\right)}_{D}=\left[\mathbf{Z}^{\left(q\right)}\right]^{\top}\mathbf{S}^{\left(q+1\right)}\mathbf{Z}^{\left(q\right)}+\mathbf{S}^{\left(q\right)}\mathbf{Z}^{\left(q-1\right)}\left[\mathbf{S}^{\left(q-1\right)}\right]^{-1}\left[\mathbf{Z}^{\left(q-1\right)}\right]^{\top}\mathbf{S}^{\left(q\right)}.$
(39)

In Eq. (39), matrix $\mathbf{Z}^{\left(q\right)}$ is defined as
$\mathbf{Z}^{\left(q\right)}=\left[\mathbf{F}^{\left(q+1\right)}\right]^{\top}$,
where boundary matrix, $\mathbf{F}^{\left(q+1\right)}$, and its transpose,
$\left[\mathbf{F}^{\left(q+1\right)}\right]^{\top}$, are the discrete
counterparts of $\partial^{\left(q+1\right)}$ and
$\partial^{\left(q+1,*\right)}$ such that the combinatorial Laplacian matrix
in Eqs. (37-38) can be expressed as
$\mathbf{L}^{\left(q\right)}_{C}=\mathbf{F}^{\left(q+1\right)}\left[\mathbf{F}^{\left(q+1\right)}\right]^{\top}+\left[\mathbf{F}^{\left(q\right)}\right]^{\top}\mathbf{F}^{\left(q\right)}$.
Matrix $\mathbf{S}^{\left(q\right)}$ is a diagonal matrix that discretizes the
Hodge star [60] on simplicial complexes. The Hodge star is an operator that
transforms $q$-forms to $\left(h-q\right)$-forms in differential geometry
theories [60]. Meanwhile, the original Hodge Laplacian,
$\Delta^{\left(q\right)}_{D}$, can be approximated by a discrete counterpart,
$\left[\mathbf{S}^{\left(q\right)}\right]^{-1}\mathbf{L}^{\left(q\right)}_{C}$
[60].

Although the (discretized) Hodge Laplacian and the combinatorial Laplacian
seem to share multiple similarities, they are not necessarily equivalent in
mathematics. On most types of simplicial complexes, the spectra of these two
Laplacian operators significantly differ from each other, which may arise from
the fact that the combinatorial Laplacian is completely determined by network
connectivity and has no geometrical constraint [60]. Here we emphasize a
common property shared by the (discretized) Hodge Laplacian and the
combinatorial Laplacian, which is also one of the intrinsic differences
between these two Laplacian operators and our proposed Laplacian operators
(i.e., the multi-order Laplacian the high-order path Laplacian in the SRG).
This difference lies in the following aspects:

  * (1)

The $\left(i,j\right)$-th element of the combinatorial Laplacian,
$\mathbf{L}^{\left(q\right)}_{C}$, and the (discretized) Hodge Laplacian,
$\mathbf{L}^{\left(q\right)}_{D}$, reflects the relations between simplex
$\sigma^{\left(q\right)}_{i}$ and simplex $\sigma^{\left(q\right)}_{j}$ [60,
61, 23], which can involve with multiple system units. The dimensions of these
two operators are determined by the number of simplices, which can differ from
the number of system units. Because a $q$-simplex denotes the un-decomposable
high-order interaction among $q+1$ units, each element of the (discretized)
Hodge Laplacian and the combinatorial Laplacian actually describes the
relations between two high-order interactions. Matrix
$\mathbf{L}^{\left(q\right)}_{C}$ reduces to the classic Laplacian in graph
theory when $q=0$ because each system unit is also a $0$-simplex. On most
kinds of simplicial complexes, matrix $\mathbf{L}^{\left(q\right)}_{D}$ is not
same as the classic Laplacian even when $q=0$.

  * (2)

The $\left(i,j\right)$-th element of the multi-order Laplacian,
$\mathbf{L}^{\left(q\right)}_{M}$, and the high-order path Laplacian,
$\mathbf{L}^{\left(q\right)}_{H}$, describes the $q$-order interactions
between unit $i$ and unit $j$. In other words, these operators characterize
the un-decomposable interactions among system units on different orders.
Therefore, the dimensions of these two operators are determined by the number
of system units. When $q=1$, matrices $\mathbf{L}^{\left(q\right)}_{M}$ and
$\mathbf{L}^{\left(q\right)}_{H}$ reduce to the classic graph Laplacian
because $1$-simplex is same as edge in networks.

Consequently, if we establish a renormalization group (RG) framework on the
(discretized) Hodge Laplacian or the combinatorial Laplacian, the targets to
renormalize are high-order interactions and the criteria used to guide coarse
graining processes are the relations among high-order interactions. When we
define a RG framework on the multi-order Laplacian or the high-order path
Laplacian, as described in the main text, the objects to renormalize are
system units and the criteria of renormalization are determined by high-order
correlations among units. Please note that the difference between these two
perspectives is neither superior nor inferior to each other. We follow the
second perspective in the SRG theory only because it is normal to choose
system units as renormalization objects in practice.

The similarity shared by the combinatorial Laplacian and the (discretized)
Hodge Laplacian is that the multiplicity of their zero eigenvalues equals the
Betti number, i.e., the number of $q$-dimensional holes [60, 61, 23]. For the
multi-order Laplacian and the high-order path Laplacian, the multiplicity of
their zero eigenvalues equals the number of connected components in the high-
order network sketches of the system. Here a $q$-order network sketch of the
system is a network where two units are connected if they share at least one
$q$-order interaction, which is useful in visualizing high-order systems of
large sizes.

#### Relations with the persistent homology theory.

A common objective of the persistent homology theory is to identify specific
topological properties of simplicial complexes that have long lifetime when we
add or remove simplices from the system [62, 74, 75, 76]. Although the multi-
order Laplacian and the high-order path Laplacian have not been used in the
persistent homology theory [62, 74, 75, 76] yet, the combinatorial Laplacian
and the (discretized) Hodge Laplacian have been extensively included in
persistent homology analysis [63, 64, 62, 65, 66, 67]. Below, we suggest a
possible idea to relate the SRG with the persistent homology theory.

Our idea arises from a simple fact. No matter we implement a SRG framework
with the multi-order Laplacian or the high-order path Laplacian, when this SRG
renormalizes the system on a $p$-order according to the properties of
$q$-order interactions ($p\leq q$), we can always record the set of units in
$\mathbf{G}^{\left(p\right)}_{k}$ that are aggregated into each super-unit in
$\mathbf{G}^{\left(p\right)}_{k+1}$. By repeating this recording in every
iteration, we can recursively find the set of initial units in
$\mathbf{G}^{\left(p\right)}_{1}$ aggregated into each super-unit in
$\mathbf{G}^{\left(p\right)}_{k+1}$ for each $k\geq 1$. In other words, we can
know how the initial structure formed by system units is changed by the SRG
after $k$ iterations. If we set $p=1$ and search simplices in the initial
system, $\mathbf{G}^{\left(1\right)}_{1}$, we can use the recursively tracked
information to identify how many simplices in
$\mathbf{G}^{\left(1\right)}_{1}$ are added, removed, or modified by the SRG
after $k$ iterations. Therefore, similar to the filtration process in the
persistent homology theory [62, 76], the SRG also offers an opportunity to
study the lifetime of different topological properties (e.g., the Betti
number). Given a simplicial complex $K$, the SRG also creates a filtration
sequence of sub-complexes $\left(K_{k}\right)_{k=0}^{m}$ (here $m$ denotes the
number of iterations in the SRG). Certainly, the SRG is different from the
common filtration processes that persistently add or remove simplices
following deterministic rules [62, 76]. It has more intricate effects on the
topological properties of the system.

##  Appendix C The SRG in a form of conventional renormalization group
theories

In this section, we contextualize the SRG in the background of conventional
renormalization group (RG) theories.

Before presenting a more in-depth analysis, we first review the Gaussian model
used in complex system analysis (e.g., see Refs. [118, 31] for instances). Let
us consider a system state, $\mathbf{x}$, whose component associated with each
system unit $i$ is denoted as $\mathbf{x}\left(i\right)$. For convenience, we
assume that the system size is infinitely large and satisfies continuity. For
an interacting system, the conventional Gaussian model can be established by
omitting the fourth and higher coupling terms of the Landau model, whose
functional Hamiltonian can be expressed as [118, 31]

$\displaystyle\mathcal{H}=\int\left[\frac{1}{2}\left(|\nabla\mathbf{x}\left(i\right)|^{2}+r\mathbf{x}\left(i\right)^{2}\right)-J\mathbf{x}\left(i\right)\right]\mathsf{d}i,$
(40)

where $r$ and $J$ denote some parameters concerning the temperature and the
external field. We assume that the system described in Eq. (40) is distributed
on a network whose adjacency matrix is $\mathbf{A}$. Given this assumption, we
can rewrite Eq. (40) as

$\displaystyle\mathcal{H}$
$\displaystyle\simeq\frac{1}{2}\sum_{i,j}\left[\frac{1}{2}\mathbf{A}_{ij}\left(\mathbf{x}\left(i\right)-\mathbf{x}\left(j\right)\right)^{2}+r\mathbf{x}\left(i\right)^{2}\delta_{ij}\right]-J\sum_{i}\mathbf{x}\left(i\right).$
(41)

We notice that Eq. (41) shares similar terms with the Laplacian dynamics
analyzed in one of our earlier works [119], which inspires us to reformulate
Eq. (41) into

$\displaystyle\mathcal{H}$
$\displaystyle\simeq\frac{1}{2}\sum_{i,j}\Bigg{[}\mathbf{x}\left(i\right)\underbrace{\Bigg{(}\delta_{i,j}\sum_{j}\mathbf{A}_{ij}-\mathbf{A}_{ij}\Bigg{)}}_{\text{A
Laplacian
operator}}\mathbf{x}\left(j\right)+r\mathbf{x}\left(i\right)^{2}\delta_{ij}\Bigg{]}-J\sum_{i}\mathbf{x}\left(i\right).$
(42)

According to Eq. (41), we can naturally obtain a Gaussian model governed by
the Laplacian operator if we place the system on a network. In fact, this same
idea has been proposed and validated in existing studies extensively (e.g.,
see Refs. [118, 120] for instances).

Now, let us relate the above analysis with the SRG framework. We can consider
high-order interactions and insert the multi-order Laplacian,
$\mathbf{L}^{\left(q\right)}_{M}$, or the high-order path Laplacian,
$\mathbf{L}^{\left(q\right)}_{H}$, into Eq. (42) to derive

$\displaystyle\mathcal{H}^{\left(q\right)}$
$\displaystyle\simeq\frac{1}{2}\sum_{i,j}\Bigg{[}\mathbf{x}\left(i\right)\left[\mathbf{L}^{\left(q\right)}\right]_{ij}\mathbf{x}\left(j\right)+r\mathbf{x}\left(i\right)^{2}\delta_{ij}\Bigg{]}-J\sum_{i}\mathbf{x}\left(i\right),$
(43)
$\displaystyle\simeq\frac{1}{2}\sum_{i,j}\Bigg{[}\mathbf{x}\left(i\right)\left(\left[\mathbf{L}^{\left(q\right)}\right]_{ij}+r\delta_{ij}\right)\mathbf{x}\left(j\right)\Bigg{]}-J\sum_{i}\mathbf{x}\left(i\right),$
(44)

where $\left[\mathbf{L}^{\left(q\right)}\right]_{ij}$ can be selected as
$\left[\mathbf{L}^{\left(q\right)}_{M}\right]_{ij}$ or
$\left[\mathbf{L}^{\left(q\right)}_{H}\right]_{ij}$. We refer to
$\mathcal{H}^{\left(q\right)}$ as the $q$-order functional Hamiltonian, which
characterizes the system behaviours on the $q$-order. If we write
$\mathbf{K}_{ij}^{\left(q\right)}=\left[\mathbf{L}^{\left(q\right)}\right]_{ij}+r\delta_{ij}$,
we notice that $\mathbf{K}_{ij}^{\left(q\right)}$ is essentially equivalent to
the term of $\phi^{4}$ theory analyzed in Ref. [31]. The only difference lies
in that the PCA-like RG proposed by Ref. [31] does not study high-order
interactions. Given this connection, we can have following ideas after
learning Ref. [31]:

  * (1)

The operator $\mathbf{L}^{\left(q\right)}$ used in the SRG mainly serves as a
foundation for determining relevant scales based on its spectrum, whose role
is similar to the precision matrix (i.e., the inverse of covariance matrix)
studied in the PCA-like RG under the condition of translation invariance [31].
To certain extents, integrating out interactions associated with large
eigenvalues of $\mathbf{L}^{\left(q\right)}$ in the SRG shares the same idea
with eliminating the modes with small variances in the PCA-like RG [31]. This
validates our idea of designing the SRG based on the spectrum of
$\mathbf{L}^{\left(q\right)}$. Apart from the relation with the PCA-like RG,
the usage of the spectrum of $\mathbf{L}^{\left(q\right)}$ is also similar
with the ideas of functional RG theories, where a regulator is used to govern
the renormalization [77, 78, 79]. Then main difference between the SRG and
functional RGs lies in that there is no explicit functional or related
formalism given in the SRG because it is hard to explicitly derive these
concepts for complex systems.

  * (2)

The non-trivial point is how to define an appropriate cutoff of the spectrum
[31]. Different from a flexible and empirically selected cutoff in the PCA-
like RG [31], we generalize the idea of Ref. [10] to define the specific heat
in both ergodic and non-ergodic cases and search for a cutoff based on the
maxima of specific heat. We choose the specific heat because it is
demonstrated as an indicator of second-order phase transitions in diffusion
processes and is informative about characteristic diffusion scales [43, 10].
Meanwhile, as suggested in Ref. [43], the specific heat is proportional to the
amplitudes of entropy fluctuations if we consider the thermal fluctuation-
dissipation theorem [121]. Moreover, the specific heat can be directly
calculated using $\mathbf{L}^{\left(q\right)}$ since diffusion on networks is
fully characterized by the Laplacian operator. In sum, the specific heat is a
valid and convenient criterion to determine a cutoff. Certainly, one can
consider other criteria in practice (e.g., the critical slowing-down of order
parameter fluctuations and the ensuing time scale separation [121]).

  * (3)

By considering a Gaussian model, we omit the fourth and higher coupling terms
among units [118, 31]. Why can we still obtain non-trivial results under this
condition? Below, we present possible explanations:

    * (a)

As shown by the perturbation analysis of the PCA-like RG [31], analyzing a
model with Gaussian fluctuations does not un-conditionally imply a trivial
consequence of the central limit theorem. Although a RG is initially proposed
for a translation invariant system with local interactions, it may also be
applicable to more complex systems and identify non-Gaussian behaviours [31].
For instance, the initial state to renormalize in diverse complex systems can
be non-Gaussian. Meanwhile, the departures from a Gaussian state can become
either more or less important under the scale transformation by the RG. These
conditions are all possible to create non-trivial phenomena.

    * (b)

Most importantly, perhaps a frequently neglected point is that
$\mathbf{K}^{\left(q\right)}$ can be defined with $q>1$ (in fact, this is the
most common case in the SRG). Under this condition, the remaining terms in Eq.
(44) can not be treated as the descriptions of conventional pairwise
couplings. Instead, they have implicitly include $q$-order interactions among
units, which require the simultaneous participation of more than two units
when $q>1$. For instance, although $\mathbf{K}_{ij}^{\left(4\right)}$ focuses
on the $4$-order relation between units $i$ and $j$, the relation itself
actually corresponds to the $4$-order interactions involved with unit $i$,
unit $j$, and three other units, which can happen over a large region and are
beyond pairwise local interactions. Consequently, by considering a high-order
functional Hamiltonian, we can study more complex correlations among units,
providing opportunities for non-trivial fixed points to occur in the SRG.

##  Appendix D The behaviour of the specific heat

As mentioned in our main text and explained in Appendix C, the specific heat
is used to define the cutoff of the spectrum of operator
$\mathbf{L}^{\left(q\right)}$, which helps determine relevant scales in the
system. Specifically, the $q$-order specific heat in Eq. (23), derived using
the first derivative of the $q$-order spectral entropy, is used to define the
optimal time scale, $\widehat{\tau}^{\left(q\right)}$, of the SRG in Eq. (24).
Then, the inverse of time scale $\widehat{\tau}^{\left(q\right)}$ is utilized
to define the number of characteristic modes for describing the system near
criticality (see step (2) in the SRG).

Figure 9: The simplex renormalization group (SRG) with different
$\left(p,q\right)$. (a) We generate a system with $200$ units whose pairwise
interactions follow a Watts-Strogatz network (each unit initially has $5$
neighbors and edges are rewired according to a probability of $0.05$) [80].
Then, the SRG is applied to renormalize the system on the $p$-order according
to the properties of $q$-order interactions ($p=1$ and $q=2$). The first five
iterations of the renormalization flows on $q$\- and $p$-orders are presented,
accompanied by their corresponding high-order path Laplacian operators. (b)
For comparison, we show the renormalization flow of the system on the
$p$-order guided by the $q$-order when $p=q=1$.

In this section, we mainly show how to derive the specific heat,
$X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)$, in Eq. (23) from
Eq. (22) when $\tau^{\left(q\right)}$ is sufficiently large. Formally, we can
solve the differential equation,
$X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)=-\frac{\mathsf{d}}{\mathsf{d}\log\left(\tau^{\left(q\right)}\right)}S^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)$,
in Eq. (22) as

$\displaystyle
X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)=-\tau^{\left(q\right)}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho}-\left(\tau^{\left(q\right)}\right)^{2}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho}-\tau^{\left(q\right)}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\log\left(\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\right),$
(45)

where the third term can be reformulated as

$\displaystyle-\tau^{\left(q\right)}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\log\left(\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\right)=$
$\displaystyle-\tau^{\left(q\right)}\frac{1}{\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right),$
(46) $\displaystyle=$
$\displaystyle\tau^{\left(q\right)}\frac{\sum_{\lambda^{\left(q\right)}_{1}}\lambda^{\left(q\right)}_{1}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)}{\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\tau^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)},$
(47) $\displaystyle=$
$\displaystyle\tau^{\left(q\right)}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho}.$
(48)

Inserting Eq. (48) into Eq. (45), we can readily derive

$\displaystyle X^{\left(q\right)}_{1}\left(\tau^{\left(q\right)}\right)$
$\displaystyle=-\left(\tau^{\left(q\right)}\right)^{2}\frac{\mathsf{d}}{\mathsf{d}\tau^{\left(q\right)}}\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\rho},$
(49)

which is the specific heat shown in our main text.

##  Appendix E Comparison of the SRG across different settings of
$\left(p,q\right)$

In Fig. 9, we use the SRG to process a system whose pairwise interactions
exhibit small-world properties [80].

We first renormalize the system on the $p$-order under the guidance of
$q$-order interactions ($p=1$ and $q=2$). As shown in Fig. 9(a), integrating
out short-range $q$-order interactions lead to an organizational-structure-
preserving renormalization flow of $p$-order interactions. Because high-order
interactions are usually distributed in a sparser manner than pairwise ones
(e.g., see the flow of $\mathbf{G}^{\left(q\right)}$), the renormalization
flow guided by $q$-order interactions reduce the system less drastically.

For comparison, we renormalize the system on the $p$-order according to the
properties of $p$-order interactions themselves (i.e., this is equivalent to
setting $p=q=1$ in the SRG) in Fig. 9(b). The generated renormalization flow
reduces the system more significantly because $p$-order interactions,
including those short-range ones, are distributed in a dense manner. The
differences between Fig. 9(a) and Fig. 9(b) qualitatively suggest the effects
of higher-order interactions on lower-order interactions during
renormalization.

Figure 10: The existence of the multi-order scale-invariance (Part I). (a-e)
The degree distributions and the Laplacian eigenvalue spectra of different
synthetic interacting systems under the transformation of the SRG are
illustrated after being averaged across all replicas. The multi-order
Laplacian operator in Eq. (1) is applied to implement the SRG.

##  Appendix F Scaling relation of the Laplacian eigenvalue spectrum in the
ergodic case

Given the general form of the Laplacian eigenvalue spectrum, i.e.,
$\operatorname{Prob}\left(\lambda^{\left(q\right)}_{1}\right)\sim\left(\lambda^{\left(q\right)}_{1}\right)^{\beta^{\left(q\right)}}=\alpha\left(\lambda^{\left(q\right)}_{1}\right)^{\beta^{\left(q\right)}}$,
we can first reformulate Eq. (25) in the main text as

$\displaystyle\frac{\widetilde{X}^{\left(q\right)}}{\widetilde{\tau}^{\left(q\right)}}$
$\displaystyle=\big{\langle}\lambda^{\left(q\right)}_{1}\big{\rangle}_{\widetilde{\rho}}+\mu,$
(50)
$\displaystyle=\frac{\frac{1}{N^{\left(q\right)}_{k}}\sum_{\lambda^{\left(q\right)}_{1}}\lambda^{\left(q\right)}_{1}\exp\left(-\widetilde{\tau}^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)}{\frac{1}{N^{\left(q\right)}_{k}}\sum_{\lambda^{\left(q\right)}_{1}}\exp\left(-\widetilde{\tau}^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)}+\mu.$
(51)

When the system is near criticality, the predominant domain of each eigenvalue
tends to be similar, making it reasonable to assume that the eigenvalue
spectrum is quasi-continuous

$\displaystyle\frac{\widetilde{X}^{\left(q\right)}}{\widetilde{\tau}^{\left(q\right)}}$
$\displaystyle=\frac{\int_{0}^{\infty}\alpha\left(\lambda^{\left(q\right)}_{1}\right)^{\beta^{\left(q\right)}+1}\exp\left(-\widetilde{\tau}^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\mathsf{d}\lambda^{\left(q\right)}_{1}}{\int_{0}^{\infty}\alpha\left(\lambda^{\left(q\right)}_{1}\right)^{\beta^{\left(q\right)}}\exp\left(-\widetilde{\tau}^{\left(q\right)}\lambda^{\left(q\right)}_{1}\right)\mathsf{d}\lambda^{\left(q\right)}_{1}}+\mu,\;\exists\mu\in\mathbb{R}.$
(52)

In Eq. (52), we can also understand $\mathsf{d}\lambda^{\left(q\right)}_{1}$
as the continuous limit of the weight associated with each eigenvalue. Using
the Euler’s gamma function, i.e.,
$\Gamma\left(z\right)=\int_{0}^{\infty}x^{z-1}\exp\left(-x\right)\mathsf{d}x$,
we can transform Eq. (52) as

$\displaystyle\frac{\widetilde{X}^{\left(q\right)}}{\widetilde{\tau}^{\left(q\right)}}$
$\displaystyle=\frac{\frac{\alpha}{\left(\widetilde{\tau}^{\left(q\right)}\right)^{\beta^{\left(q\right)}+2}}\Gamma\left(\beta^{\left(q\right)}+2\right)}{\frac{\alpha}{\left(\widetilde{\tau}^{\left(q\right)}\right)^{\beta^{\left(q\right)}+1}}\Gamma\left(\beta^{\left(q\right)}+1\right)}+\mu,$
(53)

which directly leads to

$\displaystyle\widetilde{X}^{\left(q\right)}$
$\displaystyle=\frac{\Gamma\left(\beta^{\left(q\right)}+2\right)}{\Gamma\left(\beta^{\left(q\right)}+1\right)}+\mu\widetilde{\tau}^{\left(q\right)},$
(54)
$\displaystyle=\beta^{\left(q\right)}+1+\nu,\;\exists\nu\in\left[0,\infty\right).$
(55)

Note that we have reduce $\mu\widetilde{\tau}^{\left(q\right)}$ in Eq. (54) to
$\nu$ in Eq. (55) for convenience, where $\nu$ denotes an arbitrary non-
negative number. The scope of $\nu\in\left[0,\infty\right)$ is derived from
the fact that $\beta^{\left(q\right)}\leq-1$ (i.e., as required by the power-
law distribution) and $\widetilde{X}^{\left(q\right)}\geq 0$ (i.e., the
specific heat is non-negative).

##  Appendix G Scaling relation of the Laplacian eigenvalue spectrum in the
non-ergodic case

Here we present the derivations of the scaling relation in the non-ergodic
case. For convenience, we denote
$\mathbf{L}^{\left(q\right)}_{k}\left(i\right)$ as the Laplacian of cluster
$C_{i}$ of $\mathbf{G}^{\left(q\right)}_{k}$. It is clear that the Laplacian
of $\mathbf{G}^{\left(q\right)}_{k}$ can be formulated in a block diagonal
form

$\displaystyle\mathbf{L}^{\left(q\right)}_{k}=\operatorname{diag}\left(\left[\mathbf{L}^{\left(q\right)}_{k}\left(1\right),\ldots,\mathbf{L}^{\left(q\right)}_{k}\left(r\right)\right]\right)$
(56)

because there exists no interaction between each pair of clusters. Eq. (56)
leads to an important property of the determinant

$\displaystyle\operatorname{det}\left(\mathbf{L}^{\left(q\right)}_{k}\right)$
$\displaystyle=\prod_{i=1}^{r}\operatorname{det}\left(\mathbf{L}^{\left(q\right)}_{k}\left(i\right)\right),$
(57)
$\displaystyle\prod_{\lambda^{\left(q\right)}_{k}\in\Lambda^{\left(q\right)}_{k}}\lambda^{\left(q\right)}_{k}$
$\displaystyle=\prod_{i=1}^{r}\prod_{\lambda^{\left(q\right)}_{k}\in\Lambda^{\left(q\right)}_{k}\left(i\right)}\lambda^{\left(q\right)}_{k},$
(58)

where $\Lambda^{\left(q\right)}_{k}$ and
$\Lambda^{\left(q\right)}_{k}\left(i\right)$ denote the eigenvalue sets of
$\mathbf{L}^{\left(q\right)}_{k}$ and
$\mathbf{L}^{\left(q\right)}_{k}\left(i\right)$, respectively. According to
Eqs. (57-58), the characteristic polynomial roots of
$\mathbf{L}^{\left(q\right)}_{k}$ can be constructed using the combination of
the characteristic polynomial roots of
$\\{\mathbf{L}^{\left(q\right)}_{k}\left(1\right),\ldots,\mathbf{L}^{\left(q\right)}_{k}\left(r\right)\\}$

$\displaystyle\Lambda^{\left(q\right)}_{k}=\bigcup_{i=1}^{r}\Lambda^{\left(q\right)}_{k}\left(i\right).$
(59)

Figure 11: The existence of the multi-order scale-invariance (Part II). (a-e)
The degree distributions and the Laplacian eigenvalue spectra of different
synthetic interacting systems under the transformation of the SRG are
illustrated after being averaged across all replicas. The high-order path
Laplacian operator in Eq. (5) is applied to implement the SRG.  Figure 12: The
organizational structure transformed by renormalization flows. The Louvain