# A non-isomorphic decomposition analysis of percolation, correlated systems,
and renormalization groups

Yang Tian  Department of Psychology & Tsinghua Laboratory of Brain and
Intelligence, Tsinghua University, Beijing, 100084, China.  Laboratory of
Advanced Computing and Storage, Central Research Institute, 2012 Laboratories,
Huawei Technologies Co. Ltd., Beijing, 100084, China.
<EMAIL_ADDRESS>Pei Sun  Department of Psychology & Tsinghua
Laboratory of Brain and Intelligence, Tsinghua University, Beijing, 100084,
China<EMAIL_ADDRESS>
###### Abstract

Any system with internal correlations (e.g., interactions), such as spin
glasses, random media, and brains, can be studied as a network where nodes are
elements and edges denote correlated interactions. As internal correlations
increase, correlated behaviours originate and propagate among elements to
create percolation phenomena. In this work, we study percolation in correlated
systems by proposing a framework referred to as the non-isomorphic
decomposition analysis. This framework analyzes the evolution of the space of
non-isomorphic percolation configuration classes governed by internal
correlations, where a metric called as the automorphism entropy is proposed to
measure the reduction of information (i.e., freedom degrees) in the space. A
special case of our theory under the mean-field approximation is elaborated,
where all concepts can be explicitly calculated based on group theory and
percolation theory. Our framework has been applied to studying absorbing phase
transition, spreading process, and synchronization, whose results reveal the
shared properties of percolation and dynamics in different correlated systems.
First, the percolation process driven by internal correlations is
intrinsically an information loss process in the probability space of non-
isomorphic percolation configuration classes. The metric of information loss
degree, automorphism entropy, serves as a sensitive indicator of percolation
transition and related processes in diverse systems. Second, the information
loss intensity depends on the strength of internal correlations. Sufficiently
strong correlations ensure slight information losses and make specific freedom
degrees omissible, which supports the implementation of renormalization
groups. Sufficiently weak correlations, on the contrary, imply significant
information losses and make renormalization groups less effective for the
correlated system.

  * February 2023

Keywords: percolation, correlated systems, group theory

###### Contents

  1. 1 Introduction
    1. 1.1 Correlation systems and percolation
    2. 1.2 Our framework and contribution
  2. 2 Non-isomorphic decomposition
    1. 2.1 Mean-filed approximation
    2. 2.2 Non-isomorphic percolation configuration class
  3. 3 Automorphism entropy
    1. 3.1 General definition of automorphism entropy
    2. 3.2 Precise formalization of automorphism entropy
    3. 3.3 Empirical formula of automorphism entropy
  4. 4 Infinite cluster automorphism entropy
    1. 4.1 Giant cluster in Bernoulli percolation process
    2. 4.2 Infinite cluster automorphism entropy in Bernoulli percolation
  5. 5 Automorphism entropy spectrum
    1. 5.1 Cluster size distribution and automorphism entropy spectrum
    2. 5.2 Moments of automorphism entropy spectrum
  6. 6 Non-isomorphic decomposition and the correlated systems
    1. 6.1 Non-isomorphic decomposition and absorbing phase transition
    2. 6.2 Non-isomorphic decomposition and spreading processes
    3. 6.3 Non-isomorphic decomposition and synchronization
  7. 7 Physics aspects of non-isomorphic decomposition
    1. 7.1 Coarse graining guided by non-isomorphic decomposition
    2. 7.2 Manifestation in Laplacian renormalization group
    3. 7.3 Manifestation in moment space phenomenological renormalization group
    4. 7.4 Manifestation in real space phenomenological renormalization group
    5. 7.5 Non-isomorphic decomposition and renormalization groups
  8. 8 Conclusion
  9. A A practical calculation of the Pólya enumeration theorem
  10. B Computation of automorphism entropy
  11. C Empirical estimation of automorphism entropy
  12. D Computation of the Bernoulli percolation
  13. E The probability for a node to belong to a cluster
  14. F Computation of automorphism entropy spectrum
  15. G Computation of directed percolation in branching model

##  1 Introduction

###  1.1 Correlation systems and percolation

Correlated systems (or referred to as interacting systems), such as those
described by the Wilson-Cowan equation [1, 2, 3, 4], spin glasses [5], the
random cluster model [6], the Kuramoto model [7], and the cellular Potts model
[8, 9, 10], are important research objects in physics. By its nature, a
correlated system can be represented by the dynamics on a network
$G\left(V,E\right)$, where nodes in $V$ are system elements and edges in $E$
denote correlated interactions. The existence of edge between elements $i$ and
$j$ is non-trivially determined by $C_{ij}$, the $\left(i,j\right)$-th entity
in a specific correlation matrix $C$. The definition of element correlation
can be arbitrary, which is not limited to the statistical correlation.

If we seek an abstract formalization that characterizes different correlated
systems as networks in a unified manner, percolation on network [11, 12] may
be a natural choice. The key idea is to generate an all-to-all network of
elements (i.e., a complete graph) during initialization. Then each edge is
occupied if and only if the associated elements synchronously behave (i.e.,
these two elements are correlated). All unoccupied edges are removed from the
initial network, after which a filtered network is obtained to describe the
propagation of correlated behaviours. In a probabilistic manner, we can
consider an occupation probability vector

$\displaystyle\rho=f\left(C\right)=\left(\rho_{e}:\;e=\left(i,j\right)\in
E\right)\in\left[0,1\right]^{E},$ (1)

where $f\left(\cdot\right)$ is a non-trivial element-wise function that maps
correlation matrix $C$ to occupation probability vector $\rho$. A percolation
configuration

$\displaystyle\mathbf{\chi}=\left(\chi_{e}:\;e=\left(i,j\right)\in
E\right)\in\\{0,1\\}^{E}$ (2)

lives in a probability space
$\mathbb{P}_{\rho}=\left(\\{0,1\\}^{E},\sigma_{\chi},P_{\rho}\right)$, where
$1$ denotes occupation (i.e., the edge is kept) and $0$ denotes non-occupation
(i.e., the edge is deleted). Notion $\sigma_{\chi}$ denotes a $\sigma$-algebra
and $P_{\rho}$ is a probability measure characterized by vector $\rho$. The
propagation of correlated behaviours within the system is fully characterized
by $\mathbb{P}_{\rho}$, based on which we can apply percolation theory to
study the emergence of a giant cluster of correlated elements or the evolution
of the average size of clusters [11, 13, 12].

###  1.2 Our framework and contribution

In this work, we focus on the associate dynamics in the statistical ensemble
of percolation configurations during the percolation process of correlated
behaviours. We suggest a new perspective to analyze percolation and correlated
systems in terms of the variation of information (i.e., freedom degrees) in
percolation configuration ensemble.

Mathematically, we can explore this perspective by comparing
$\mathbb{P}_{\rho}$ with
$\mathbb{P}_{\infty}=\left(\\{0,1\\}^{E},\sigma_{\chi},P_{\infty}\right)$, a
space that includes the percolation configurations generated by all possible
relations among system elements (i.e., element relation can be arbitrary
rather than be limited to $\rho$) and reflects the neutral properties of the
target system when there is no constraint on element relations. We pursue to
develop a possible approach to realize this comparison and suggest how it
inspires new techniques in studying percolation and correlated systems.

In Sec. 2, we develop a framework, referred to as the non-isomorphic
decomposition analysis, to transform our concerned question into the space of
non-isomorphic percolation configuration classes. The change degree of
information in $\mathbb{P}_{\rho}$ compared with that in $\mathbb{P}_{\infty}$
can be reflected by the evolution of the probability measure of non-isomorphic
percolation configuration classes, which is more practical to analyze. In this
paper, we primarily introduce a special case of the non-isomorphic
decomposition analysis under the mean-field approximation, yet more general
versions of our framework can be freely developed without mean-field
approximation.

In Sec. 3, we develop a metric referred to as automorphism entropy to measure
the reduction of information (i.e., the loss of freedom degrees) in the space
of non-isomorphic percolation configuration classes as a function of $\rho$.
Automorphism entropy can be calculated based on the orbit-stabilizer theorem
and its lemmas [14]. A smaller automorphism entropy corresponds to a lower
degree of information loss when element relations are constrained by $\rho$.

In Sec. 4, we present a calculation instance of the automorphism entropy in an
abstract case where the propagation of correlated behaviours among the
elements of a system is characterized by Bernoulli percolation on a network
[13]. We propose the infinite cluster automorphism entropy to measure the
information loss degree associated with the infinite cluster of correlated
elements. The quantified information loss is observed to serve as an indicator
of percolation transition. This instance suggests the quantitative relation
between information losses in the space of non-isomorphic percolation
configuration classes and phase transitions in correlated systems.

In Sec. 5, we propose the automorphism entropy spectrum to measure the
information loss degree associated with the finite cluster of correlated
elements. Similar to the infinite cluster automorphism entropy, the evolution
of the automorphism entropy spectrum and its moments can be used to predict or
analyze phase transitions controlled by internal correlations.

In Sec. 6, we apply our framework to analyzing the dynamics in several
representative correlated systems, including the absorbing phase transition in
branching models [15, 16, 17], spreading processes in epidemic models [18,
19], and synchronization in oscillator models [20, 7]. The phase transition
phenomena in these correlated systems, irrespective of their distinct details,
can all be reflected and analyzed in terms of information loss in percolation
configuration ensemble.

In Sec. 7, we elaborate the physics idea underlying the non-isomorphic
decomposition analysis. We suggest how our concerned question relates to
renormalization group theories [21, 22, 23], based on which, an abstract
coarse graining approach is proposed under the guidance of the non-isomorphic
decomposition analysis. This abstract framework has concrete manifestations in
various existing renormalization groups, such as the Laplacian renormalization
group [24], the moment space phenomenological renormalization group [25, 26],
and the real space phenomenological renormalization group [27, 28]. The non-
isomorphic decomposition analysis enables us to study when and how a
correlated system can be manipulated by renormalization groups in an effective
way.

In Sec. 8, we present a summary of our work and suggests potential directions
to generalize our framework.

Although our work is primarily presented in an abstract form, we endow equal
importance to the computational practicality of our framework. Therefore, we
elaborate lots of details of computer programming in the appendixes. Moreover,
because numerous concepts across different branches of statistical physics are
used in our work, it is less feasible to assign a unique notion for every
concept. As a compromise, some notions are slightly abused if they are
contextually dependent and do not cause confusions.

##  2 Non-isomorphic decomposition

To realize our target comparison, we analyze the spaces of non-isomorphic
percolation configuration classes associated with $\mathbb{P}_{\infty}$ and
$\mathbb{P}_{\rho}$ separately. Our analysis framework is referred to as the
non-isomorphic decomposition analysis because it begins with decomposing the
percolation configuration ensemble into numerous non-isomorphic parts. We
elaborate the basic framework of non-isomorphic decomposition analysis by
formalizing a special case of it under the mean-filed approximation.

###  2.1 Mean-filed approximation

We first narrow down our target question to homogeneous systems, where a mean-
field approximation can be applied on Eq. (1). Specifically, all elements in
the correlated system are assumed to be homogeneous such that the occupation
probability vector in Eq. (1) can be reasonably replaced by a mean occupation
probability

$\displaystyle\rho=\frac{1}{|E|}\sum_{\left(i,j\right)\in
E}g\left(C_{ij}\right),$ (3)

where $g\left(\cdot\right)$ is an arbitrary mapping that normalizes all
entities in $C$ into $\left[0,1\right]$ (e.g., the min-max normalization).
Based on Eq. (3), the percolation process underlying the propagation of
correlated behaviours reduces to the Bernoulli percolation on networks [13].
Please note that the mean-filed approximation is not a necessary condition of
our theory. It only corresponds to a special case of our theory that is easy
to interpret (see a calculation instance of our theory without the mean-filed
approximation in Sec. 6).

Given a mean occupation probability $\rho$ in Eq. (3), we can define the
corresponding $\mathbb{P}_{\rho}$ and analyze $\rho_{n}$, the associated
probability for a random node to belong to the giant cluster (i.e., an
infinite connected component with endless branching process). Under the mean-
filed approximation (i.e., the probability is invariant across different
nodes), the giant cluster size is

$\displaystyle N=\lfloor\rho_{n}|V|\rfloor,$ (4)

where $\lfloor\cdot\rfloor$ is the floor function. Please note that although
the giant cluster is theoretically infinite, its embodiment within a finite
correlated system must have a finite size (i.e., $N$) as well. Meanwhile, we
can also analyze the size distribution of finite clusters (i.e., an connected
component with finite branching process) following the approach proposed by
Refs. [29, 30, 31, 32, 33] (see Sec. 4 for detailed derivations).

The difference between $\mathbb{P}_{\infty}$ and $\mathbb{P}_{\rho}$ manifests
as the sets of nodes included by specific infinite or finite clusters. In
space $\mathbb{P}_{\rho}$, the nodes included by an arbitrary cluster are
required to form a connected component following a certain wiring diagram. In
space $\mathbb{P}_{\infty}$, there is no restraint on the wiring diagram of
any set of nodes. This difference implies the distinct freedom degrees of
these two spaces and may serve as a window to explore our target question.

###  2.2 Non-isomorphic percolation configuration class

Assuming that $V_{c}\subset V$ is the node set of an arbitrary cluster (i.e.,
infinite or finite) of size $M$, we can see the arbitrariness of $V_{c}$,
i.e., we only require $|V_{c}|=M$ and have no constraint on which set of $M$
nodes should be included (nodes are homogeneous to each other). Therefore, our
analysis of the difference mentioned above always involves plentiful
duplicated cases. To prune these redundancies, we can transform our analysis
into the space of non-isomorphic percolation configuration classes.

For convenience, we consider a complete graph $G_{c}\left(V_{c},E_{c}\right)$
defined on $M$ nodes. We mark $E_{\chi}\subset E_{c}$ as the edge set
corresponding to a percolation configuration $\chi\in\\{0,1\\}^{E_{c}}$

$\displaystyle E_{\chi}=\\{e|e=\left(i,j\right)\in E_{c},\;\chi_{e}=1\\}.$ (5)

We consider a symmetric group $S_{|E_{c}|}$ that acts on the filtered edge set
$E_{\chi}$, i.e., a permutation acts on any edge $\left(i,j\right)\in
E_{\chi}$ and maps it to
$\left(\sigma\left(i\right),\sigma\left(j\right)\right)$. A permutation
$\sigma:E_{\chi}\rightarrow E^{\sigma}_{\chi}$ is treated as auto-isomorphism
if

$\displaystyle\left(i,j\right)\in
E_{\chi}\Leftrightarrow\left(\sigma\left(i\right),\sigma\left(j\right)\right)\in
E^{\sigma}_{\chi}.$ (6)

Then, we define the automorphism group $\mathsf{Aut}\left(\chi\right)$ as the
set of all automorphisms of $E_{\chi}$ under mapping composition. Based on it,
we can define an equivalence relation $\sim$ between percolation
configurations

$\displaystyle\chi\sim\chi^{\prime}\Leftrightarrow\left(E_{\chi^{\prime}}=E^{\sigma}_{\chi},\;\sigma\in\mathsf{Aut}\left(\chi\right)\right).$
(7)

The space of non-isomorphic percolation configuration classes is formally
defined as

$\displaystyle\Gamma\left(M\right)=\\{0,1\\}^{E_{c}}/\sim,$ (8)

which corresponds to a partition of set $\\{0,1\\}^{E_{c}}$. Each element in
$\gamma$ is a set of percolation configurations that are isomorphic. Each pair
of elements in $\Gamma$ are not equivalent to each other in terms of $\sim$
(i.e., they are not isomorphic).

Now, let us consider the probability space defined on $\Gamma$. We begin with
the case where the probability of each percolation configuration is defined by
$\mathbb{P}_{\infty}$. Although $\mathbb{P}_{\infty}$ can not be explicitly
defined because it is impossible to enumerate all possible element relations,
a simple idea can be applied to make the analysis possible. If one percolation
configuration class $\gamma$ is less to occur given a relation among elements,
there will be another relation that makes $\gamma$ occur with a higher
probability. Therefore, we can reasonably consider a symmetry on $\Gamma$,
i.e., a postulate of _equal a priori probabilities_

$\displaystyle
P_{\infty}\left(\gamma\right)=\frac{1}{|\Gamma\left(M\right)|}.$ (9)

Please note that Eq. (9), similar to the postulate of _equal a priori
probabilities_ for equilibrium system in statistical mechanics, has no strict
proof yet. It is used as a basic assumption in our theory.

Then we turn to the case where the probability of percolation configuration is
defined by $\mathbb{P}_{\rho}$. For convenience, we define
$L\left(\gamma\right)\leq M$ as the largest connected component size of the
network associated with each percolation configuration $\chi$ in class
$\gamma$ (note that the largest connected component size is invariant across
all configurations in $\gamma$ because these configurations are isomorphic).
The probability of percolation configuration class $\gamma$ is defined as

$\displaystyle
P_{\rho}\left(\gamma\right)=\frac{1}{Z_{\gamma}}\Theta\left(L\left(\gamma\right)-M\right),$
(10)

where $\Theta\left(\cdot\right)$ is the Heaviside step function and
$Z_{\gamma}$ is a partition function

$\displaystyle
Z_{\gamma}=\sum_{\gamma\in\Gamma\left(M\right)}\Theta\left(L\left(\gamma\right)-N\right).$
(11)

Mathematically, term $\frac{1}{Z_{\gamma}}$ can be understood as a re-
normalized version of $P_{\infty}\left(\gamma\right)$, where all percolation
configurations whose largest connected component sizes are smaller than $M$
have been excluded (i.e., $V_{c}$ forms a connected component of size $N$ so
it always corresponds to the percolation configurations whose largest
connected component sizes are $M$).

As we suggest below, the above definitions enable us to implement our target
comparison in a practical way.

##  3 Automorphism entropy

Given the basic framework in Sec. 2, we can compare between
$\mathbb{P}_{\infty}$ and $\mathbb{P}_{\rho}$ by proposing a possible metric
of the variation degree of information in $\mathbb{P}_{\rho}$ controlled by
$\rho$. The developed metric, as a fundamental concept of our non-isomorphic
decomposition analysis, is referred to as the automorphism entropy.

###  3.1 General definition of automorphism entropy

Our metric arises from a natural idea. Let us measure the variability of
$\Gamma$ in terms of the Shannon entropy $H$

$\displaystyle H_{\infty}\left(M\right)$
$\displaystyle=\sum_{\gamma\in\Gamma\left(N\right)}P_{\infty}\left(\gamma\right)\log\frac{1}{P_{\infty}\left(\gamma\right)},$
(12) $\displaystyle H_{\rho}\left(M\right)$
$\displaystyle=\sum_{\gamma\in\Gamma\left(N\right)}P_{\rho}\left(\gamma\right)\log\frac{1}{P_{\rho}\left(\gamma\right)}.$
(13)

It is trivial to know

$\displaystyle H_{\rho}\left(M\right)\leq H_{\infty}\left(M\right)$ (14)

because numerous freedom degrees (i.e., the percolation configurations whose
largest connected component sizes are smaller than $M$) have been excluded
while defining $P_{\rho}\left(\cdot\right)$. The variability of
$\Gamma\left(M\right)$ is larger given the symmetric structure in
$\mathbb{P}_{\infty}$ while it is smaller given the limited freedom degrees in
$\mathbb{P}_{\rho}$.

A natural metric of the latent information loss can be defined as

$\displaystyle
A\left(M\right)=1-\frac{H_{\rho}\left(M\right)}{H_{\infty}\left(M\right)}\in\left[0,1\right],$
(15)

which is referred to as the automorphism entropy since it is defined in the
space of non-isomorphic percolation configuration classes. A smaller
automorphism entropy suggests a smaller degree of information loss.

###  3.2 Precise formalization of automorphism entropy

Below, we propose a calculation approach of Eq. (15) based on group theory.
Key steps of this approach is summarized in Fig. 1.

Figure 1:  Conceptual illustrations of the non-isomorphic decomposition
analysis. (a) summarizes the key steps of calculating $|\Gamma|$ and
$|\Gamma^{*}|$, based on which, other concepts of the non-isomorphic
decomposition analysis can be readily derived. (b) illustrates the instance of
spaces $\Gamma$ and $\Gamma^{*}$ on a correlated system with a gaint connected
component of size $4$ (i.e., $M=4$).

We first deal with the calculation of Eq. (9), which essentially requires us
to measure $|\Gamma\left(M\right)|$, the total number of non-isomorphic
percolation configuration classes that can emerge from $M$ nodes.
Mathematically, we can calculate $|\Gamma\left(M\right)|$ following the Pólya
enumeration theorem [34], a generalization of the lemma of the orbit-
stabilizer theorem [14]

$\displaystyle|\Gamma\left(M\right)|=\frac{1}{|S_{M}|}\sum_{\omega\in
S_{N}}2^{c\left[\sigma\left(\omega\right)\right]},$ (16)

where $2=|\\{0,1\\}|$ measures the number of possible occupation states and
$c\left(\cdot\right)$ measures the number of cycles in the permutation. Notion
$\sigma\left(\omega\right)\in S_{|E_{c}|}$ denotes the corresponding
permutation on edges given $\omega\in S_{M}$, the permutation on nodes. In
general, there are two ways to calculate Eq. (16). The first way depends the
precise derivation of cycle index, a multivariate generating function, of
group $S_{|E_{c}|}$ and may require more manual derivations [34, 35]. The
second way, shown in A, utilizes the properties of group orbits and is easier
to follow. In our work, the second approach is applied, which calculates
$|\Gamma\left(M\right)|$ following Eqs. (88-89) and solves Eq. (9) as

$\displaystyle
P_{\infty}\left(\gamma\right)=\frac{\sum_{O\left(\omega\right)\in\mathbf{O}_{M}}|O\left(\omega\right)|}{\sum_{O\left(\omega\right)\in\mathbf{O}_{M}}|O\left(\omega\right)|2^{c\left[\sigma\left(\omega\right)\right]}},$
(17)

where $O\left(\omega\right)$ denotes the conjugate class of $\omega$ and
$\mathbf{O}_{M}$ is the set of all conjugate classes of group $S_{M}$. Based
on Eq. (17), entropy $H_{\infty}\left(M\right)$ in Eq. (12) can be derived.

Then we turn to calculating Eq. (10), which requires us to count the size of
$\Gamma^{*}\left(M\right)\subset\Gamma\left(M\right)$, the set of all
percolation configuration class associated with connected graphs on $M$ nodes
(i.e., all nodes belong to the same connected component). Please see Fig. 1
for the illustrations of $\Gamma\left(M\right)$ and $\Gamma^{*}\left(M\right)$
in a case with $M=4$. Although the counting task of
$|\Gamma^{*}\left(M\right)|$ is hard to be implemented with the Pólya
enumeration theorem [34], previous studies have suggested that the sequence of
$\\{|\Gamma^{*}\left(1\right)|,|\Gamma^{*}\left(2\right)|,|\Gamma^{*}\left(3\right)|,\ldots\\}$
is the inverse Euler transform of the sequence of
$\\{|\Gamma\left(1\right)|,|\Gamma\left(2\right)|,|\Gamma\left(3\right)|,\ldots\\}$
[36, 34]. For convenience, we introduce a new sequence
$\\{\Lambda\left(1\right),\Lambda\left(2\right),\Lambda\left(3\right),\ldots\\}$
where $\Lambda\left(1\right)=|\Gamma\left(1\right)|$. Then we can have the
following relations

$\displaystyle\Lambda\left(M\right)$
$\displaystyle=N|\Gamma\left(M\right)|-\sum_{i=1}^{M-1}\Lambda\left(i\right)|\Gamma\left(M-i\right)|,$
(18) $\displaystyle|\Gamma^{*}\left(M\right)|$
$\displaystyle=\frac{1}{M}\sum_{i|M}\phi\left(\frac{M}{i}\right)\Lambda\left(i\right),$
(19)

where $\phi\left(\cdot\right)$ denotes the Möbius function and notion $i|M$
represents that $i$ divides $M$. As a special case of the empty graph, i.e.,
$M=0$, we define $|\Gamma\left(0\right)|=1$ and $|\Gamma^{*}\left(0\right)|=1$
because the empty graph can have any property [37]. Based on Eqs. (18-19), we
can transform Eqs. (10-11) as

$\displaystyle
P_{\rho}\left(\gamma\right)=\frac{M}{\sum_{i|M}\phi\left(\frac{M}{i}\right)\Lambda\left(i\right)},$
(20)

which enables us to calculate entropy $H_{\rho}\left(M\right)$ in Eq. (13).

Figure 2:  Calculation of the automorphism entropy as a function of $M$. (a)
shows fraction $\frac{|\Gamma^{*}|}{|\Gamma|}$ as a function of $M$ in a log-
log plot. The inserted sub-plot shows the sequences of $|\Gamma|$,
$|\Gamma^{*}|$, and $|\Lambda|$ for reference. (b) illustrates the
automorphism entropy across different $M$ in a log-log plot. The actual values
of entropies $H_{\infty}$ and $H_{\rho}$ are shown in the inserted sub-plot.

Based on the above derivations, we can precisely define the automorphism
entropy as a function of mean occupation probability $\rho$. Specifically,
given a correlated system, we can measure $\rho$ to derive $M$, based on which
we obtain $P_{\infty}\left(\cdot\right)$ and $P_{\rho}\left(\cdot\right)$ to
calculate the automorphism entropy $A\left(\cdot\right)$. In our research, we
primarily define $A\left(0\right)=A\left(1\right)=1$ because they share an
ill-posed expression of $1-\frac{0}{0}$ following Eq. (15). Other possible
definitions (e.g., $A\left(0\right)=A\left(1\right)=0$) are also acceptable
because they do not affect our analysis. In B, we offer a step-by-step
explanation of the code for computing the automorphism entropy.

In Fig. 2, we present the calculation results of the automorphism entropy with
$1\leq M\leq 50$. We first show the value of
$\frac{|\Gamma^{*}\left(M\right)|}{|\Gamma\left(M\right)|}$ as a function of
$N$. While a sharp drop of
$\frac{|\Gamma^{*}\left(M\right)|}{|\Gamma\left(M\right)|}$ occurs when $M$
changes from $1$ to $2$, a continuous and bounded increase of
$\frac{|\Gamma^{*}\left(M\right)|}{|\Gamma\left(M\right)|}$ can be observed
once $N\geq 2$. In other words, as $M$ enlarges, an increasing proportion of
non-isomorphic percolation configuration classes in space $\Gamma$ corresponds
to connected graphs. We also show the automorphism entropy across different
values of $M$. As $M$ increases, the automorphism entropy rapidly reduces to
$0$, suggesting the decrease of information loss degree according to the
definition of the automorphism entropy.

###  3.3 Empirical formula of automorphism entropy

Apart from analytic calculation, we also present a numerical estimation of the
automorphism entropy based on the least squares method. The simple variation
trend of the automorphism entropy in Fig. 2 ensures an accurate fitting of it
for $M\geq 2$ (we exclude $A\left(0\right)$ and $A\left(1\right)$ from fitting
because their definitions can be arbitrary). Specifically, an exponential
decay of $A\left(M\right)$ is obtained with high accuracy ($R^{2}=1$ and
$RMSE=0.0005618$)

$\displaystyle
A\left(M\right)=3.956^{+0.014}_{-0.014}\exp(-0.6876^{+0.0015}_{-0.0015}M),\;M\geq
2$ (21)

where $3.956^{+0.014}_{-0.014}$ (i.e., $\left[3.942,3.970\right]$) and
$-0.6876^{+0.0015}_{-0.0015}$ (i.e., $\left[-0.6891,-0.6861\right]$) denote
the intervals of regression coefficients with 95% confidence, respectively.
This empirical formula is more convenient in the case with a sufficiently
large $M$ (e.g., $M=10^{5}$), where the analytic calculation method presented
above becomes computationally costly (i.e., the analytic method requires
integer partition). Please see C for the code implementation of the empirical
formula.

##  4 Infinite cluster automorphism entropy

###  4.1 Giant cluster in Bernoulli percolation process

In this section, we present a calculation instance of the automorphism entropy
in a concrete percolation process. We implement a Bernoulli percolation
process controlled by occupation probability $\rho$ in an initial network with
$100$ nodes (note that the initial network is a complete graph). The code is
presented in D.

We primarily focus on the automorphism entropy defined on the giant cluster
(i.e., infinite cluster), which is referred to as the infinite cluster
automorphism entropy, $A_{\infty}$, in our work. To calculate this quantity,
we need to deal with the calculation of $N$, the size of the giant cluster of
correlated elements in the concerned system. Because we have reduced our
question to the Bernoulli percolation on networks [13] based on mean-field
approximation, the precise value of $N$ given a mean occupation probability
$\rho$ can be readily calculated by the generating function approach.
Specifically, the probability for a random node to belong to the giant
cluster, $\rho_{n}$, and the probability for a random edge to belong to the
giant cluster, $\rho_{e}$, obey the following relations

$\displaystyle\rho_{e}$
$\displaystyle=1-\sum_{k}P_{k}\frac{k}{\sum_{k}kP_{k}}\left(1-\rho_{e}\rho\right)^{k-1},$
(22) $\displaystyle\rho_{n}$
$\displaystyle=1-\sum_{k}P_{k}\left(1-\rho_{e}\rho\right)^{k},$ (23)

where $P_{k}$ denotes the probability of finding a node with $k$ degrees in
the corresponding graph (i.e., the initial network) of the correlated system
[11, 32, 33]. Because the initial network is set as a complete graph [11, 12],
i.e., $P_{k}=\delta\left(k,|V|-1\right)$ where
$\delta\left(\cdot,\cdot\right)$ denotes the Kronecker delta, we can
reformulate Eqs. (22-23) as

$\displaystyle\rho_{e}$ $\displaystyle=1-\left(1-\rho_{e}\rho\right)^{|V|-2},$
(24) $\displaystyle\rho_{n}$
$\displaystyle=1-\left(1-\rho_{e}\rho\right)^{|V|-1}.$ (25)

In practice, one can first solve Eq. (24) and insert $\rho_{e}$ into Eq. (25)
to derive $\rho_{n}$. Given a $\rho_{n}$, the size of the giant cluster of
correlated elements, $N$, can be directly obtained following Eq. (4).

###  4.2 Infinite cluster automorphism entropy in Bernoulli percolation

Figure 3:  Calculation of the infinite cluster automorphism entropy during the
Bernoulli percolation process. (a) shows the order parameter $\rho_{n}$ of the
Bernoulli percolation as a function of $\rho$. The inserted sub-plot presents
$\rho_{e}$ as a function of $\rho$ for reference. (b) illustrates the infinite
cluster automorphism entropy across different values of $\rho$. The inserted
sub-plot shows the size of the giant connected component as a function of
$\rho$.

In Fig. 3, we show the order parameter, $\rho_{n}$, of the Bernoulli
percolation and the associated value of $N$ as the functions of $\rho$.
Meanwhile, we calculate the associated infinite cluster automorphism entropy
$A_{\infty}\left(N\right)$ at each value of $N$ for a qualitative comparison.
As $\rho$ increases from $0$ to $1$, a percolation transition emerges at a
certain threshold, after which a system-level connected component with a size
of $\mathcal{O}\left(|V|\right)$ occurs in the network (here
$\mathcal{O}\left(\cdot\right)$ denotes the big O notion of growth rate). A
coherent variation trend can be observed in $A_{\infty}\left(N\right)$, where
the percolation transition threshold is near a change point for
$A_{\infty}\left(N\right)$ to exhibit a sharp reduction. Therefore, the
infinite cluster automorphism entropy has the potential to serve as an
indicator of percolation transition. The interval of $\rho$ in which
$A_{\infty}\left(N\right)=1$ principally corresponds to the period before
percolation transition. The change point at which $A_{\infty}\left(N\right)$
starts to depart from $1$ is close to the transition threshold. The reduction
process of $A_{\infty}\left(N\right)$ is also the process for a giant
connected component to include more nodes.

After qualitatively suggesting the sensitivity of the infinite cluster
automorphism entropy towards percolation transition, we turn to deriving the
quantitative relation between the infinite cluster automorphism entropy and
percolation transition. According to Eq. (25) and our setting, we know
$A_{\infty}\left(N\right)<1$ if and only if $N>2$ and
$A_{\infty}\left(N\right)=1$ when $N\in\\{0,1,2\\}$. Therefore, the change
point $\rho_{c}^{A}$ of $A_{\infty}\left(N\right)$ in terms of occupation
probability is $\rho_{c}^{A}=\frac{2}{N}$. Meanwhile, we know that the
percolation threshold $\rho_{c}$ of the Bernoulli percolation follows a form
[11, 32, 33]

$\displaystyle\rho_{c}=\left(\frac{\partial}{\partial
x}\sum_{k}P_{k}\frac{k}{\sum_{k}kP_{k}}x^{k-1}\Bigg{|}_{x=1}\right)^{-1},$
(26)

which can be reformulated as

$\displaystyle\rho_{c}=\frac{\mathbb{E}\left(k\right)}{\mathbb{E}\left(k^{2}\right)-\mathbb{E}\left(k\right)}=\frac{1}{N-2}$
(27)

because the initial network is a complete graph. Consequently, we can obtain
the following relation

$\displaystyle\rho_{c}=\frac{N}{2N-4}\rho_{c}^{A}.$ (28)

At a limit of $N\rightarrow\infty$, Eq. (28) reduces to
$\rho_{c}=\frac{1}{2}\rho_{c}^{A}$. In sum, the change point of
$A_{\infty}\left(N\right)$ is proportional to the percolation threshold
according to $N$. The infinite cluster automorphism entropy can be used to
derive the percolation threshold in arbitrary correlated systems.

##  5 Automorphism entropy spectrum

###  5.1 Cluster size distribution and automorphism entropy spectrum

In the previous sections, we have introduced the infinite cluster automorphism
entropy for comparing between $\mathbb{P}_{\infty}$ and $\mathbb{P}_{\rho}$.
This metric quantifies the reduction of freedom degree in $\mathbb{P}_{\rho}$
and is primarily defined based on the giant cluster (i.e., the infinite
cluster where the branching process is endless) emerged during percolation.
Although this idea is natural and comprehensible, we can not ignore the fact
that the giant cluster may not cover all nodes completely (i.e., there may
exist finite clusters in the network). To fill this gap, we generalize our
calculation to a spectrum of automorphism entropies measured for all finite
clusters.

Let us consider the distribution of finite clusters during percolation [29,
30, 31, 32, 33]. The probability for a randomly selected node to belong to a
cluster of size $r\in\mathbb{N}^{+}$ given an occupation probability $\rho$ is
calculated as (see Refs. [29, 30] for detailed derivations of Eq. (29))

$\displaystyle\pi_{r}\left(\rho\right)=\frac{\rho^{r-1}\sum_{k}kP_{k}}{\left(r-1\right)!}\left[\frac{\partial^{r-2}}{\partial
x^{r-2}}\left(\sum_{k=0}^{\infty}\frac{\left(k+1\right)P_{k+1}}{\sum_{k}kP_{k}}x^{k}\right)^{r}\right]\Bigg{|}_{x=1-\rho},\;r\geq
2,$ (29)

which reduces to

$\displaystyle\pi_{r}\left(\rho\right)$
$\displaystyle=\frac{\rho^{r-1}\left(|V|-1\right)}{\left(r-1\right)!}\left[\frac{\partial^{r-2}}{\partial
x^{r-2}}x^{r\left(|V|-2\right)}\right]\bigg{|}_{x=1-\rho},\;r\geq 2,$ (30)
$\displaystyle=\frac{\rho^{r-1}\left(|V|-1\right)}{\left(r-1\right)!}\prod_{i=0}^{r-3}\left[r\left(|V|-2\right)-i\right]\left(1-\rho\right)^{r\left(|V|-3\right)+2},\;r\geq
3$ (31)

in our case (i.e., the initial network is a complete graph with $|V|$ nodes).
For the special case with $r=1$ (i.e., the node is isolated), we have
$\pi_{1}\left(\rho\right)=\left(1-\rho\right)^{|V|-1}$ (i.e., all the
associate edges of the node in the initial network are not occupied). For the
special case with $r=2$ (i.e., the node is connected with only one node) in
Eq. (30), Eq. (31) is reformulated as
$\pi_{2}\left(\rho\right)=\rho\left(|V|-1\right)\left(1-\rho\right)^{2|V|-4}$.

In fact, we can derive a distribution of $\pi_{r}\left(\rho\right)$ in a much
simpler manner. Let us consider $Q_{k}$, the degree distribution of the
diluted network. Given an occupation probability $\rho$, it is trivial to know
that $Q_{k}$ is nothing else but the Poisson degree distribution

$\displaystyle Q_{k}$
$\displaystyle\simeq\exp\left(\rho|V|-\rho\right)\frac{\left(\rho|V|-\rho\right)^{k}}{k!}.$
(32)

The distribution of $\pi_{r}\left(\rho\right)$ in a network with a Poisson
degree distribution can be readily derived as (note that we can no longer
repeatedly analyze occupation while calculating $\pi_{r}\left(\rho\right)$
because the network has already been diluted before deriving $Q_{k}$)

$\displaystyle\pi_{r}\left(\rho\right)$
$\displaystyle=\frac{\sum_{k}kQ_{k}}{\left(r-1\right)!}\left[\frac{\partial^{r-2}}{\partial
x^{r-2}}\left(\sum_{k=0}^{\infty}\frac{\left(k+1\right)Q_{k+1}}{\sum_{k}kQ_{k}}x^{k}\right)^{r}\right]\Bigg{|}_{x=0},\;r\geq
2,$ (33) $\displaystyle=\frac{\exp\left(\rho
r-\rho|V|r\right)\left(\rho|V|r-\rho r\right)^{r-1}}{r!},\;r\geq 2.$ (34)

One can easily prove the equivalence between Eq. (31) and Eq. (34). In E, we
present our code for computing $\pi_{r}\left(\rho\right)$.

Given the probability for a randomly selected node to belong to a cluster of a
given size, we can derive $\theta_{r}\left(\rho\right)$, the distribution of
cluster size in the diluted network, according to the following relation [11,
32, 33]

$\displaystyle\pi_{r}\left(\rho\right)=\frac{r\theta_{r}\left(\rho\right)}{\sum_{r}r\theta_{r}\left(\rho\right)}.$
(35)

Eq. (35) holds because the probability for a randomly selected node to belong
to a cluster is proportional to the size of this cluster. The derivation of
$\theta_{r}\left(\rho\right)$ based on Eq. (35) is presented in F. The
obtained result, organized according to the Kronecker delta
$\delta\left(\cdot,\cdot\right)$, is

$\displaystyle\theta_{r}\left(\rho\right)$ $\displaystyle=$
$\displaystyle\left(\frac{\exp\left(\rho r-\rho|V|r\right)\left(\rho|V|r-\rho
r\right)^{r-1}}{r^{2}\left(r-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{1-\delta\left(1,r\right)}\left(1+\sum_{k=2}^{|V|}\frac{\exp\left(\rho
k-\rho|V|k\right)\left(\rho|V|k-\rho
k\right)^{k-1}}{k^{2}\left(k-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{-1}$
(36)

for any $r\geq 1$.

Figure 4:  Calculation of the automorphism entropy spectrum during the
Bernoulli percolation process. (a) shows the distribution of
$\pi_{r}\left(\rho\right)$ given different $\rho$. (b) shows the distribution
of $\theta_{r}\left(\rho\right)$ given different $\rho$. (c) illustrates the
automorphism entropy spectrum across different values of $\rho$.

Based on the non-trivial cluster size distribution presented in Eq. (36), we
can define the automorphism entropy spectrum as $\zeta_{A}\left(\rho\right)$,
the probability distribution of the automorphism entropies measured for every
possible cluster. In the case where the empirical formula of the automorphism
entropy in Eq. (25) is applied, the spectrum can be directly defined as

$\displaystyle\zeta_{A}\left(\rho\right)$
$\displaystyle=\theta_{r}\left(\rho\right),\;A<1,$ (37) $\displaystyle r$
$\displaystyle=\frac{1}{-0.6876^{+0.0015}_{-0.0015}}\ln\left(\frac{1}{3.956^{+0.014}_{-0.014}}A\right).$
(38)

Please note that
$\zeta_{1}\left(\rho\right)=\theta_{1}\left(\rho\right)+\theta_{2}\left(\rho\right)$
because $A\left(1\right)=A\left(2\right)=1$. An efficient function for
computing the automorphism entropy spectrum is proposed in F . Certainly, Eq.
(38) can also be replaced by the analytic expression of the automorphism
entropy (see Eqs. (16-20)). Here we use empirical formula only for
convenience.

Figure 5:  The evolution of the automorphism entropy spectrum near percolation
transition. (a) shows several instances of the evolution of
$\zeta_{A}\left(\rho\right)$ where $A\neq 1$. (b) illustrates the evolution of
$\zeta_{A}\left(\rho\right)$ where $A=1$ (we use $A\left(1\right)$ as an
instance). Given
$\zeta_{1}\left(\rho\right)=\theta_{1}\left(\rho\right)+\theta_{2}\left(\rho\right)$,
the inserted sub-plot also presents the evolution of
$\theta_{1}\left(\rho\right)$ and $\theta_{2}\left(\rho\right)$ for reference.

In Fig. 4, we present the instances of $\pi_{r}\left(\rho\right)$,
$\theta_{r}\left(\rho\right)$, and $\zeta_{A}\left(\rho\right)$ given
different values of $\rho$. These instances are calculated during the
Bernoulli percolation process in a network with $100$ nodes, where $\rho$
increases from $\rho_{c}$ to $10\rho_{c}$. As shown in Fig. 4, the increase of
$\rho$ on a large scale (e.g., from $\rho_{c}$ to $10\rho_{c}$) seems to
robustly imply a sharper decay of the probability for the automorphism entropy
to reach at a smaller value. However, our subsequent analysis on a small scale
(e.g., near $\rho_{c}$) suggests that this phenomenon is not universal and
only serves as a part of a latent property of the automorphism entropy
spectrum.

As demonstrated in Fig. 5, the automorphism entropy spectrum exhibits a clear
change at $\rho_{c}$, the percolation transition threshold. Given the same
experimental setting of Fig. 4, we primarily focus on the cases where $\rho$
is near $\rho_{c}$. Computation results suggest that the value of
$\zeta_{A}\left(\rho\right)$ ($A\neq 1$) increases when $\rho<\rho_{c}$ and
decreases when $\rho>\rho_{c}$. A maximum value of
$\zeta_{A}\left(\rho\right)$ ($A\neq 1$) emerges at the percolation transition
threshold. An exactly opposite variation trend can be seen on
$\zeta_{A}\left(\rho\right)$ ($A=1$), whose value is minimized at the
percolation transition threshold. Therefore, the phenomenon observed in Fig. 4
is only a special case of the evolution of the automorphism entropy spectrum
when $\rho>\rho_{c}$, which can not be treated as universal.

If we summarize the above finding in terms of the condensation (i.e., the
degree to which the probability density is condensed at a sub-set of its
support) of probability distribution, the following conclusion on the property
of the automorphism entropy spectrum can be drawn: The condensation degree of
the automorphism entropy spectrum at $A=1$ decreases before percolation
transition and is minimized at the transition threshold. After percolation
transition, the condensation degree at $A=1$ gradually increases. The
condensation degree of the automorphism entropy spectrum in the interval of
$A>1$ increases before percolation transition and is maximized at the
transition threshold. The condensation degree in this interval decreases after
percolation transition.

Figure 6:  The evolution of the moments of the automorphism entropy spectrum
near percolation transition. (a) shows several instances of the evolution of
$\mathbb{E}\left(A^{m}\right)$ where $m\in\\{1,\ldots,20\\}$. (b) illustrates
$\Big{|}\frac{\partial}{\partial\frac{\rho}{\rho_{c}}}\mathbb{E}\left(A^{m}\right)\Big{|}$,
the absolute change rate of the moments for $m\in\\{1,\ldots,20\\}$, near
percolation transition.

###  5.2 Moments of automorphism entropy spectrum

Given the evolution of the automorphism entropy spectrum near percolation
transition, we are also interested in its embodiment in the moments of the
automorphism entropy spectrum.

Let us consider the $m$-th order central moment of the spectrum

$\displaystyle\mathbb{E}\left(A^{m}\right)=\sum_{A}\zeta_{A}\left(\rho\right)\left(A-\mathbb{E}\left(A\right)\right)^{m},\;m\geq
1.$ (39)

In Fig. 6, we calculate the moments of the automorphism entropy spectrum based
on the data in Fig. 5. For convenience, we primarily derive the first twenty
moments as instances. As shown in Fig. 6, all these moments exhibit a kind of
convex-function-like evolution and are minimized at the percolation transition
threshold. Meanwhile, we also measure
$\Big{|}\frac{\partial}{\partial\frac{\rho}{\rho_{c}}}\mathbb{E}\left(A^{m}\right)\Big{|}$,
the absolute change rate with respect to $\frac{\rho}{\rho_{c}}$, for these
moments. It turns out that
$\Big{|}\frac{\partial}{\partial\frac{\rho}{\rho_{c}}}\mathbb{E}\left(A^{m}\right)\Big{|}$
is minimized at the transition threshold as well. According to these
properties, we suggest that the evolution of the automorphism entropy spectrum
and its moments may serve as indicators to reflect percolation transition.

##  6 Non-isomorphic decomposition and the correlated systems

To this point, we have introduced the infinite cluster automorphism entropy
and the automorphism entropy spectrum as two basic concepts in the non-
isomorphic decomposition analysis. Our previous analyses are presented in an
abstract manner, where all the details of correlated systems are omitted to
offer a clear vision. Below, we apply our theory on concrete examples of
correlated systems with full details to show its potential in correlated
system analysis.

###  6.1 Non-isomorphic decomposition and absorbing phase transition

Our first concerned system is the branching model. This model is of interest
because it and its variants have extensive applications across different
fields, such as neuroscience [38, 39, 40, 41, 42, 43], social science [44, 45,
46], and geophysics [47, 48].

For convenience, we consider a branching model defined on a square lattice
with periodic boundary conditions [41]. In this model, all elements have two
possible states, i.e., being active or silent. Each element only interacts
with its four neighbors. A silent element has a small probability (e.g.,
$10^{-4}$ in our research) to be spontaneously activated in every time step
and it will return to silence in the next time step (i.e., resetting). An
active element has a probability of $p$ to activate each of its neighbors,
creating a propagation process of activation. Given a small $p$, the
propagation process rapidly vanishes because of resetting. As $p$ increases,
activation may propagate among elements for a longer time. Given a
sufficiently large $p$, the propagation process of activation may even become
self-sustaining, i.e., the propagation is strong enough to suppress resetting.
A self-sustaining propagation process can last for an infinite time, which is
intrinsically different from those who are not self-sustaining and only last
for a finite time. Therefore, two phases can be distinguished from each other
in the phase space. The first phase, referred to as the active phase,
corresponds to the cases where propagation is self-sustaining. The second
phase, referred to as the absorbing phase, corresponds to the cases where
propagation vanishes in a finite time. The boundary between these two phases
corresponds to the absorbing phase transition characterized by the directed
percolation universality class, where the avalanches of activation emerge in
the system [15, 16, 17]. One can see reviews of these concepts in terms of
neuroscience in Refs. [43, 49] to acquire more detailed and pictorial
explanations. Here we no longer elaborate these explanations and directly
focus on the mathematical properties of the absorbing phase transition.

Mathematically, the branching model exhibits several scaling behaviours when
it is captured by the directed percolation universality class. Specifically,
the emerged avalanches during phase transition are exhibit the following
power-law behaviours [43, 49]

$\displaystyle P\left(t\right)$ $\displaystyle\sim t^{-\alpha},$ (40)
$\displaystyle P\left(s\right)$ $\displaystyle\sim s^{-\beta},$ (41)

where $t$ and $s$ denote the lifetime and size of an avalanche. Meanwhile, the
mean avalanche size given a certain lifetime, denoted by
$\mathbb{E}\left(s\left(t\right)\right)$, is expected to follow an exponential
relation

$\displaystyle\mathbb{E}\left(s\left(t\right)\right)$ $\displaystyle\sim
t^{\gamma}.$ (42)

In the directed percolation universality class, a scaling relation can be
derived for the above exponents [43, 49]

$\displaystyle\gamma=\frac{\alpha-1}{\beta-1}.$ (43)

Apart from the scaling relation, an universal collapse should be observed in
the scaled average temporal shape of avalanches. Specifically, there is

$\displaystyle\mathbb{E}\left(s\left(\tau|t\right)\right)$
$\displaystyle=t^{\psi}h\left(\frac{\tau}{t}\right),\;\forall\tau\in\left[0,t\right],$
(44) $\displaystyle\psi$ $\displaystyle=\gamma-1,$ (45)

where $f\left(\cdot\right)$ denotes a universal scaling function. In most
cases, a parabolic form of $f\left(\cdot\right)$ can be derived [50]. On can
see the derivation of Eqs. (40-45) in Refs. [15, 16, 17, 43].

Figure 7:  The power-law behaviours, scaling relation, and universal collapse
shape of avalanches. (a) shows the power-law behaviour of avalanche size $s$,
suggesting that $\beta=1.453$. (b) shows the power-law behaviour of avalanche
lifetime $t$, suggesting that $\alpha=1.625$. (c) illustrates the exponential
behaviour of $\mathbb{E}\left(s\left(t\right)\right)$, deriving that
$\gamma=1.4892$. The scaling relation is principally satisfied because
$\frac{\alpha-1}{\beta-1}=1.3796$. (d) illustrates the scaled temporal shapes
of different avalanches (marked by colorful solid lines). Meanwhile, a
quadratic polynomial is defined to fit $f\left(\cdot\right)$ for the scaled
average temporal shape (marked by black dashed line) to derive that
$\psi=0.49\simeq\gamma-1$.

In Fig. 7, we implement a branching model with $100$ elements and
$p=\frac{1}{4}$. Here we set $p=\frac{1}{4}$ because lattices are squares and
$p=\frac{1}{4}$ implies that each active element can activate $1$ neighbor on
average. According to Ref. [41], this setting is sufficient to create
avalanches. This property can be understood in terms of the Galton-Watson
branching process, where the critical point corresponds to the case where an
activate element gives birth to $1$ child on average [51]. In the experiment,
each element generates an activity sequence with a time length of $10^{6}$ to
support counting avalanches. Given the generated data, we verify Eqs. (40-45)
using the statistical toolbox proposed in Ref. [41]. As shown in Fig. 7, the
power-law behaviours, scaling relation, and universal collapse shape defined
in Eqs. (40-45) are principally satisfied by the avalanche data. Although some
errors are inevitably caused by data noise and finite size effects, their
magnitudes are controlled within a reasonable interval.

Now, let us consider a general case where we lack _a priori_ knowledge about
the branching model, i.e., element iteration rules and lattice properties are
unknown, and the model is involved with non-negligible noise. Although Eqs.
(40-45) are expected to be satisfied when the system is critical (i.e., the
necessary condition), we may not be able to use them to precisely determine if
an unknown branching model is critical or not (i.e., the sufficient
condition). This is because the latent power-law behaviours, scaling relation,
and universal collapse shape of the system can be easily disturbed or covered
by data noise, finite size effects, or other factors. As we demonstrate in G,
the determination of criticality based the latent power-law behaviours,
scaling relation, and universal collapse shape can be amphibolous or even
misleading, especially when we lack information about the system property.
This may serve as a reason for the controversies about criticality in diverse
correlated systems (e.g., see those in the brain [43, 49]).

Figure 8:  The evolution of the infinite cluster automorphism entropy and the
automorphism entropy spectrum as $p$ increases. (a) visualizes the
automorphism entropy measured on the giant cluster across different $p$. (b)
shows several instances of the evolution of $\mathbb{E}\left(A^{m}\right)$ as
$p$ increases, where $m\in\\{1,\ldots,20\\}$. (c) illustrates
$\Big{|}\frac{\partial}{\partial p}\mathbb{E}\left(A^{m}\right)\Big{|}$, the
absolute change rate of the moments for $m\in\\{1,\ldots,20\\}$, as $p$
increases.

To verify the latent criticality of real data in a more precise manner, we
need to ensure that the indicator is significant (i.e., change drastically at
the critical point) and robust (i.e., can not be covered by data noise)
enough. Inspired by the evolution of the automorphism entropy measured on the
giant cluster, $A_{\infty}$, and the automorphism entropy spectrum,
$\zeta_{A}$, near percolation transition, we suggest them as possible
candidates. We present a possible approach to calculate $A_{\infty}$ and
$\zeta_{A}$ in any branching model (please note that our approach only serves
as an example and other approaches are also acceptable). The branching model
is defined with $100$ elements and $p\in\left[0,0.4\right]$. Under each
condition of $p$, the model runs for $10^{6}$ time steps to generate element
activities, based on which we can estimate a Pearson correlation matrix
$C\left(p\right)$ among these elements. Then we calculate a mean absolute
correlation value, $\mathbb{E}_{i,j}|C_{ij}\left(p\right)|$, given each $p$. A
mean occupation probability $\rho$ is defined as
$\rho=g\left(\mathbb{E}_{i,j}|C_{ij}\left(p\right)|\right)$, where
$g\left(\cdot\right)$ is a linear function that normalizes
$\left[g\left(\mathbb{E}_{i,j}|C_{ij}\left(0\right)|\right),\ldots,g\left(\mathbb{E}_{i,j}|C_{ij}\left(0.4\right)|\right)\right]$
into an interval of $\left[0,1\right]$ (i.e., the min-max normalization).
Finally, our concerned $A_{\infty}$ and $\zeta_{A}$ can be derived given each
mean occupation probability $\rho$.

In Fig. 8, we show the infinite cluster automorphism entropy, $A_{\infty}$,
across different values of $\rho$, suggesting that $A_{\infty}$ is sensitive
to change of $p$ and minimized near $p=0.25$. Then we calculate the first
twenty central moments of the automorphism entropy spectrum under each
condition of $p$. As shown in our results, all these moments exhibit sharp
changes near $p=0.25$, the point where the branching model is captured by the
directed percolation universality class. Quantitatively, we measure the
absolute change rate of each moment as a function of $p$, demonstrating that
$p=0.25$ implies the most dramatic change of each moment of the spectrum.
Compared with the results in G, our calculated $A_{\infty}$ and $\zeta_{A}$
can offer a much more clear indication of the latent criticality.

An unexpected finding is that the automorphism entropy spectrum and its
moments seem to exhibit less dramatic but still observable changes near
$p=0.1$. We can not provide a clear explanation about this phenomenon yet and
suggest it as a remaining question for future explorations.

###  6.2 Non-isomorphic decomposition and spreading processes

Spreading processes, such as the epidemic spreading, are important topics in
statistical physics due to their high application values [11]. A spreading
process essentially corresponds to a correlated system because spreading
relays on element interactions.

To offer clear vision, let us begin with the simple SIR model defined with
$|V|$ elements (i.e., agents) [18, 19]. Each element has three possible
states, i.e., being susceptible (S), infectious (I), or recovered (R). For
convenience, we consider a case where there is no constraint on element
adjacency relations such that each element is free to interact all other
elements as its neighbors. This is a special case of the SIR model on
networks, where the network is a uncorrelated complete graph (see the general
case in Ref. [11]). The infection rate of the epidemic is $\xi$ such that an
infectious element has a probability of $\xi\mathsf{d}t$ to make another
susceptible element become infectious in an infinitesimal time interval
$\mathsf{d}t$. Given a recovery rate $\mu$, each infectious element has a
recovery probability of $\mu\mathsf{d}t$ to recover and become susceptible
after $\mathsf{d}t$. These definitions allow us to derive the evolution
equations

$\displaystyle\frac{\partial}{\partial t}I\left(t\right)$
$\displaystyle=\xi\left(|V|-1\right)S\left(t\right)\Theta\left(t\right)-\mu
I\left(t\right),$ (46) $\displaystyle\frac{\partial}{\partial
t}S\left(t\right)$
$\displaystyle=-\xi\left(|V|-1\right)S\left(t\right)\Theta\left(t\right),$
(47) $\displaystyle\frac{\partial}{\partial t}R\left(t\right)$
$\displaystyle=\mu I\left(t\right),$ (48)

where $I\left(t\right)$, $S\left(t\right)$, and $R\left(t\right)$ denote the
fraction of infections, susceptible, and recovered elements among elements at
moment $t$. Notion $\Theta\left(t\right)$ measures the probability for an
interaction (i.e., an edge) from an infected element links to another infected
element

$\displaystyle\Theta\left(t\right)$
$\displaystyle=\frac{\left(|V|-2\right)I\left(t\right)}{|V|-1},$ (49)

whose integral during a duration can be calculated as

$\displaystyle\Phi\left(t\right)$
$\displaystyle=\int_{0}^{t}\Theta\left(t^{\prime}\right)\mathsf{d}t^{\prime},$
(50) $\displaystyle=\frac{|V|-2}{\mu\left(|V|-1\right)}R\left(t\right).$ (51)

Figure 9:  The evolution of the infinite cluster automorphism entropy and the
automorphism entropy spectrum as $\tau$ changes. (a) presents the infinite
cluster automorphism entropy across different
$\frac{\left(|V|-2\right)\left(2|V|-4\right)}{\tau|V|}$. (b) shows the
instances of the evolution of $\mathbb{E}\left(A^{m}\right)$ across different
$\tau$, where $m\in\\{1,\ldots,20\\}$. (c) illustrates
$\Big{|}\frac{\partial}{\partial p}\mathbb{E}\left(A^{m}\right)\Big{|}$, the
absolute change rate of the moments for $m\in\\{1,\ldots,20\\}$, when $\tau$
increases.

At the time limit $t\rightarrow\infty$, infection vanishes and all elements
are either susceptible or recovery [11, 18, 19]. The steady state of this
system can be discovered according to

$\displaystyle 0$ $\displaystyle=\frac{\partial}{\partial
t}\Phi\left(t\right)\Big{|}_{t\rightarrow\infty},$ (52)
$\displaystyle=\frac{|V|-2}{|V|-1}I\left(\infty\right),$ (53)
$\displaystyle=\frac{|V|-2}{|V|-1}\left(1-S\left(\infty\right)\right)-\mu\Phi\left(\infty\right).$
(54)

Meanwhile, given the initial conditions, $S\left(0\right)\simeq 1$,
$I\left(0\right)\simeq 0$, and $R\left(0\right)=0$, of Eqs. (46-48), we can
derive a solution

$\displaystyle
S\left(t\right)=\exp[-\xi\left(|V|-1\right)\Phi\left(t\right)].$ (55)

Combining Eqs. (54-55), we have the self-consistent equation of
$\Phi\left(\infty\right)$

$\displaystyle\Phi\left(\infty\right)=\frac{|V|-2}{\mu\left(|V|-1\right)}\left(1-\exp[-\xi\left(|V|-1\right)\Phi\left(\infty\right)]\right).$
(56)

Eq. (56) only has a solution $\Phi\left(\infty\right)=0$ (i.e., there is no
epidemic process) if

$\displaystyle\frac{\partial}{\partial\Phi\left(\infty\right)}\left\\{\frac{|V|-2}{\mu\left(|V|-1\right)}\left(1-\exp[-\xi\left(|V|-1\right)\Phi\left(\infty\right)]\right)\right\\}\Bigg{|}_{\Phi\left(\infty\right)=0}<1.$
(57)

Therefore, the epidemic threshold, i.e., the threshold of $\frac{\xi}{\mu}$
that allows the epidemic to break out, can be discovered as

$\displaystyle\frac{\partial}{\partial\Phi\left(\infty\right)}\left\\{\frac{|V|-2}{\mu\left(|V|-1\right)}\left(1-\exp[-\xi\left(|V|-1\right)\Phi\left(\infty\right)]\right)\right\\}\Bigg{|}_{\Phi\left(\infty\right)=0}$
$\displaystyle=1,$ (58) $\displaystyle\frac{\xi_{c}}{\mu_{c}}$
$\displaystyle=\frac{1}{|V|-2}.$ (59)

The result in Eq. (59) is exactly the percolation threshold of a complete
graph (i.e., see Eq. (27) for explanations). Therefore, the phase transition
of the above epidemic process can readily predicted by $A_{\infty}$, the
infinite cluster automorphism entropy. Specifically, we have
$\rho_{c}^{A}=\frac{2|V|-4}{|V|}\frac{\xi_{c}}{\mu_{c}}$ according to Eq. (28)
and Eq. (59). As for $\zeta_{A}$, the automorphism entropy spectrum, we can
define occupation probability $\rho$ as the probability of effective epidemic
spreading controlled by $\frac{\xi}{\mu}\in\left[0,1\right]$. A small value of
$\rho=\frac{\xi}{\mu}$ corresponds to the case where epidemic spreading is
less effective and vanishes soon. To distinguish between the phases of this
system at the time limit, we can consider the outbreak time scale [11, 18, 19]

$\displaystyle\tau=\frac{\mu\mathbb{E}\left(k\right)}{\xi\left(\mathbb{E}\left(k^{2}\right)-\mathbb{E}\left(k\right)\right)}=\frac{\mu}{\xi\left(|V|-2\right)}=\frac{\mu}{\xi}\frac{\xi_{c}}{\mu_{c}},$
(60)

where a small $\tau$ indicates a sufficiently fast spreading of the infection
and $\tau=1$ holds at the epidemic threshold.

In Fig. 9, we show the evolution of the infinite cluster automorphism entropy
and the automorphism entropy spectrum across different values of $\tau$. The
epidemic spreading is defined on a system with $100$ elements, i.e.,
$|V|=100$. By simple calculation, we can know that the transition point of
$A_{\infty}$ in such a system is expected as
$\rho_{c}^{A}=\frac{\left(|V|-2\right)\left(2|V|-4\right)}{\tau|V|}\simeq
0.2$. Our results in Fig. 9 are consistent with this calculation. Meanwhile,
it can be seen that all the moments of the automorphism entropy spectrum are
minimized at $\tau=1$, which corresponds to the epidemic threshold. Meanwhile,
the epidemic threshold implies the local minima of the absolute change rates
of these moments. In sum, the evolution behaviours of the infinite cluster
automorphism entropy and the automorphism entropy spectrum may be used to
distinguish between and further study different phases of epidemic outbreak, a
typical instance of spreading processes.

One thing needs to be mentioned about the automorphism entropy spectrum. In
the more universal cases where there are certain constraints on element
adjacency relations, i.e., the epidemic spreading happens on an arbitrary
network, our derivation of the spectrum can be generalized as well.
Specifically, one needs to derive the degree distribution $P_{k}$ in the given
network and insert it into Eq. (29) to recalculate distribution
$\pi_{r}\left(\rho\right)$, cluster size distribution, and all subsequent
results. Note that Eqs. (30-36) can no longer be followed because they are
limited to complete graphs. Here we do not elaborate these generalizations in
details.

###  6.3 Non-isomorphic decomposition and synchronization

Synchronization is a widespread phenomenon in diverse correlated systems, such
as brains, social networks, and bird flocks [52, 53, 54, 55]. The Kuramoto
model, a model of coupled phase oscillators, is a representative framework for
analyzing synchronization [20, 7]. In our subsequent analysis, we apply the
non-isomorphic decomposition analysis to studying the change of freedom
degrees of the Kuramoto model during synchronization.

For convenience, we consider the classic Kuramoto model with $|V|$ oscillators
[20, 7]

$\displaystyle\frac{\partial}{\partial
t}\theta_{i}\left(t\right)=\omega_{i}+\frac{\kappa}{|V|}\sum_{i=1}^{|V|}\sin\left[\theta_{j}\left(t\right)-\theta_{i}\left(t\right)\right],$
(61)

where $\theta_{i}$ and $\omega_{i}$ denote the phase and natural frequency of
the $i$-th oscillator, respectively. Parameter $\kappa$ measures the coupling
strength among oscillators. The order parameter of the above correlated system
is

$\displaystyle
r\left(t\right)\exp\left[\mathsf{i}\psi\left(t\right)\right]=\frac{1}{|V|}\sum_{i=1}^{|V|}\exp\left[\mathsf{i}\theta_{j}\left(t\right)\right],$
(62)

in which $\psi\left(t\right)$ is the time-dependent average phase and
$r\left(t\right)\in\left[0,1\right]$ measures the synchronization degree of
the system. Based on the order parameter, we can reorganize Eq. (61) as

$\displaystyle\frac{\partial}{\partial
t}\theta_{i}\left(t\right)=\omega_{i}+\kappa
r\left(t\right)\sin\left[\psi\left(t\right)-\theta_{i}\left(t\right)\right].$
(63)

According to Eq. (63), a feedback loop exists between coupling strength
$\kappa$ and order parameter $r$, i.e., any increment in $r$ due to the
increasing $\kappa$ will enlarge the effective coupling among oscillators and
attract more oscillators to the synchronous populations in return [20, 7].

We are interested in the non-trivial evolution of the infinite cluster
automorphism entropy and the automorphism entropy spectrum while the system
becomes increasingly synchronous (e.g., by enlarging coupling strength). We
define a Kuramoto model with $100$ oscillators. Each oscillator has a random
natural frequency uniformly selected from $\left[1,5\right]$ and exhibits
activities for $1000$ time steps. All oscillators are initialized with a
random phase selected from $\left[-\pi,\pi\right]$. We set an increasing
coupling strength $\kappa\in\left[0,5\right]$ and repeat the experiment under
each condition of $\kappa$. In Fig. 10, we show the sequence of order
parameter $r\left(t\right)$ across different $\kappa$. Given a sufficiently
large $\kappa$, the sequence of $r\left(t\right)$ gradually converges to $1$
as the number of time step accumulates, suggesting that all the oscillators
become strongly synchronous. Given a small $\kappa$, the sequence of
$r\left(t\right)$ exhibits irregular oscillations rather than converge to a
certain value, corresponding to the case where oscillators are not
synchronous. These results can be validated according to
$\mathbb{E}\left(r\right)$, the mean order parameter derived by time-averaging
under each condition of $\kappa$, as well. A sharp increase of
$\mathbb{E}\left(r\right)$ can be seen near $\kappa\simeq 2$.

Figure 10:  The evolution of the infinite cluster automorphism entropy and the
automorphism entropy spectrum as the synchronization degree changes. (a)
illustrates the time series of $r\left(t\right)$ under each condition of
$\kappa$. As $\kappa$ increases from $1$ to $5$, the color of the line of
$r\left(t\right)$ changes from blue to green. (b) presents the mean order
parameter (upper parallel) and the infinite cluster automorphism entropy
(bottom parallel) across different $\kappa$. (c) shows the instances of the
evolution of $\mathbb{E}\left(A^{m}\right)$ across different $\kappa$, where
$m\in\\{1,\ldots,20\\}$.

To relate these results with our theory, we define a correlation matrix $C$ in
terms of synchronization

$\displaystyle C_{ab}$
$\displaystyle=\mathbb{E}\left[r_{ab}\left(t\right)\right],$ (64)
$\displaystyle
r_{ab}\left(t\right)\exp\left[\mathsf{i}\psi\left(t\right)\right]$
$\displaystyle=\frac{1}{2}\left\\{\exp\left[\mathsf{i}\theta_{a}\left(t\right)\right]+\exp\left[\mathsf{i}\theta_{b}\left(t\right)\right]\right\\},$
(65)

where the $\left(a,b\right)$-th entity measures the local synchronization
degree between oscillators $a$ and $b$. Then we generate the network of
oscillators, where there exists an edge between oscillators $a$ and $b$ only
if $C_{ab}\simeq 1$, i.e., these two oscillators are strongly synchronous.
Here we do not require $C_{ab}=1$ because numerical errors may lead to a non-
one synchronization degree even if oscillators $a$ and $b$ are strictly
synchronous (during computation, we determine $C_{ab}\simeq 1$ if $C_{ab}\geq
1-10^{-3}$). Given the network under each condition of $\kappa$, we can
calculate the infinite cluster automorphism entropy and the automorphism
entropy spectrum. Although the calculation lacks analytic expressions (i.e.,
the formulas of infinite cluster size and finite cluster size distribution
remain highly non-trivial), we can empirically count cluster sizes in each
network to support calculation. In Fig. 10, the derived infinite cluster
automorphism entropy is shown as a function of $\kappa$. It can be seen that
the infinite cluster automorphism entropy exhibits drastic changes when the
mean synchronization degree of the system, $\mathbb{E}\left(r\right)$, starts
to increase sharply. Meanwhile, the moments of the automorphism entropy
spectrum maintain generally constant before $\mathbb{E}\left(r\right)$ sharply
increases near $\kappa\simeq 2$. Once $\mathbb{E}\left(r\right)$ increases,
these moments starts to oscillate and slightly decrease. Therefore, the
evolution of the infinite cluster automorphism entropy and the automorphism
entropy spectrum can reflect the emerging synchronization within the system.

Apart from the discussion on synchronization, our analysis presented above
also serves as an instance of the non-isomorphic decomposition analysis
without mean-field approximation, i.e., our network generation approach
presented above does not require a mean occupation probability $\rho$ and
defines each edge according to its corresponding local synchronization degree.
This instance highlights the fact that mean-field approximation is not a
necessary condition of our theory. However, we note that excluding mean-field
approximation does create difficulties for analytic derivations. In the
future, more theoretical explorations are necessary for establishing a non-
mean-field version of our framework.

##  7 Physics aspects of non-isomorphic decomposition

In our previous sections, we have defined the basic concepts of the non-
isomorphic decomposition analysis and presented calculation instances in
concrete correlated systems. In this section, we return to the very beginning
and summarize the physics idea underlying the proposed framework.

###  7.1 Coarse graining guided by non-isomorphic decomposition

As we have mentioned before, our framework is motivated by analyzing
percolation and correlated systems in terms of the variation of information
(i.e., freedom degrees) in percolation configuration ensemble. This analysis
can be implemented by comparing between $\mathbb{P}_{\rho}$, a space
corresponding to the target properties of a correlated system with constrained
element relations, and $\mathbb{P}_{\infty}$, a space corresponding to the
neutral properties of the correlated system when there is no constraint on
element relations.

We first note that $\mathbb{P}_{\rho}$ and $\mathbb{P}_{\infty}$ actually
correspond to two extreme cases of a correlated system. According to the
definition, space $\mathbb{P}_{\infty}$ denotes an extreme case where all
system sates share the same probability to occur. Space $\mathbb{P}_{\rho}$
stands for an extreme case where only the system states corresponding to given
element correlations are allowed to occur (e.g., any system state in which the
clusters of correlated elements do not satisfy the desired infinite cluster
size or finite cluster size distribution is not allowed to occur). A real
correlated system always places itself at somewhere between these two extreme
cases. This is because internal correlations manifest themselves in
probabilistic forms and can not ensure the property of a real system in a
deterministic way (e.g., a real system with probability $\rho$ has a giant
cluster of size $N$ on average rather than deterministically). Therefore, if
we denote the actual probability space where the system lives in as
$\mathbb{P}_{\mathsf{ac}}$ and measure its information quantity as
$H_{\mathsf{ac}}$, we can extend Eq. (14) to

$\displaystyle H_{\rho}\leq H_{\mathsf{ac}}\leq H_{\infty}.$ (66)

Eq. (66) holds because $\mathbb{P}_{\mathsf{ac}}$ is more flexible in
maintaining freedom degrees than $\mathbb{P}_{\rho}$ while
$\mathbb{P}_{\infty}$ actually maximizes entropy (i.e., the uniform
distribution in a finite interval is the maximum entropy distribution).

In our work, we propose the infinite cluster automorphism entropy,
$A_{\infty}$, and the automorphism entropy spectrum, $\zeta_{A}$, to measure
the difference in freedom degrees between $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\infty}$ that are related to the infinite and finite clusters,
respectively. Because these clusters are pre-determined by a correlation
matrix, i.e., a matrix of element correlation strength, our proposed
$A_{\infty}$ and $\zeta_{A}$ naturally evolve as internal correlation changes.
As shown in Secs. 4-6, the increase of internal correlation, irrespective of
being related to percolation, emerged avalanches, the outbreak of spreading
process, synchronization, or other mechanisms, generally reduces $A_{\infty}$
and $\zeta_{A}$. In other words, as the system becomes increasingly
correlated, the informational difference between $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\infty}$ gradually decreases.

Given a sufficiently small difference between $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\infty}$, the difference between $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\mathsf{ac}}$ must be small as well

$\displaystyle H_{\mathsf{ac}}-H_{\rho}\leq H_{\infty}-H_{\rho}.$ (67)

Therefore, we can find omissible freedom degrees and approximate
$\mathbb{P}_{\mathsf{ac}}$ based on the properties of $\mathbb{P}_{\rho}$.
Specifically, for an arbitrary cluster of multiple nodes, we can safely omit
the possibilities where these nodes are not connected in
$\mathbb{P}_{\mathsf{ac}}$. This is because omitting these freedom degrees
does not imply a significant information loss, i.e., the derived
$\mathbb{P}_{\rho}$ is still similar to $\mathbb{P}_{\mathsf{ac}}$.

The operation described above is essentially similar to coarse graining in
renormalization group [21, 22, 23], where omissible freedom degrees (i.e.,
those with low contributions or effects) are eliminated as well. Given this
similarity, we may find the physics connotations of our proposed non-
isomorphic decomposition analysis, i.e., our framework may convey information
about when and how a real correlated system in $\mathbb{P}_{\mathsf{ac}}$ can
be transformed into $\mathbb{P}_{\rho}$ by specific coarse graining
approaches.

In this section, we propose a possible coarse graining approach based on the
non-isomorphic decomposition analysis, which enables us to discuss the
relation between our framework and diverse renormalization groups [21, 22,
23].

Let us consider an arbitrary correlated system with a correlation matrix $C$.
If element relations are weak (i.e., the absolute values of entities in $C$
are small), coarse graining on elements may be improper because every element
is unique and independent. In an opposite case, each certain group of elements
with strong correlations can be replaced by a super-element to realize coarse
graining. In more general cases, we can identify elements to coarse grain
using a function of the correlation matrix, denoted by $K=y\left(C\right)$.
Here $y\left(\cdot\right)$ is an appropriate function that satisfies
application demands and ensures

$\displaystyle\frac{\partial}{\partial|C_{ij}|}y\left(C_{ij}\right)\geq 0,$
(68)

for any $\left(i,j\right)$ pair.

Based on the above definition, we develop a possible way to coarse grain the
system:

  * (A)

The coarse graining is initialized by representing the correlated system as a
complete graph
$G^{\left(0\right)}\left(V^{\left(0\right)},E^{\left(0\right)}\right)$ whose
adjacency matrix is defined as $|C|$.

  * (B)

In the $k$-th iteration of coarse graining, we generate $H^{\left(k\right)}$,
a copy of $G^{\left(k\right)}$, as a reference. Then the edges corresponding
to small values in $K^{\left(k\right)}$, the function of the correlation
matrix in the $k$-th iteration, are removed from $H^{\left(k\right)}$ one by
one until there remain $n_{k}$ clusters (i.e., connected components) in
$H^{\left(k\right)}$. Here $n_{k}$ is a non-trivial parameter that requires a
careful definition. Each remaining cluster in $H^{\left(k\right)}$ contains a
set of $M$ elements with relatively strong correlations (i.e., such that the
associate values in $K^{\left(k\right)}$ can be large). According to the non-
isomorphic decomposition analysis, this cluster can be equivalently analyzed
as a network with $M$ nodes, which has a sufficiently small automorphism
entropy because of strong correlations (e.g., $A_{\infty}$ and $\zeta_{A}$).
Therefore, the difference between the associate spaces, $\mathbb{P}_{\rho}$
and $\mathbb{P}_{\infty}$, is small (note that here $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\infty}$ denote the probability spaces of non-isomorphic
percolation configuration classes that can form on $M$ nodes). According to
Eq. (67), the difference between $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\mathsf{ac}}$ is small as well. Therefore, we can treat the
correlated elements contained in this cluster as always coherent (i.e., they
always form a connected component as required by $\mathbb{P}_{\rho}$). The
possibilities that they may occasionally be incoherent can be eliminated
without significant information loss because $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\mathsf{ac}}$ have become sufficiently similar for these
elements. Consequently, we can analyze these elements as a whole and replace
them by a super-element.

  * (C)

The coarse graining can be done in $G^{\left(k\right)}$ by aggregating this
set of nodes and replacing them by a super-node. After realizing this coarse
graining for every cluster in $H^{\left(k\right)}$, we can generate
$G^{\left(k+1\right)}$ to continue the iteration.

Please note that the reference graph $H^{\left(k\right)}$ in step (B) is
proposed only for convenience (i.e., for idea clarification and computational
programming). We can also realize step (B) without the help of
$H^{\left(k\right)}$.

The key steps of the above approach are to define an appropriate function
$y\left(\cdot\right)$ and a meaningful $n_{k}$ in each $k$-th iteration, which
can vary across different application demands. Below, we suggest their
manifestations in existing renormalization groups [21, 22, 23].

###  7.2 Manifestation in Laplacian renormalization group

We first consider the Laplacian renormalization group (LRG) [24]. In the
spirit of free field theories [21, 22, 23], a Laplacian is defined in the LRG

$\displaystyle L_{ij}=\delta\left(i,j\right)\sum_{k}|C_{ik}|-|C_{ij}|$ (69)

using the Kronecker delta function $\delta\left(\cdot,\cdot\right)$ and a
normalized network propagator

$\displaystyle K=y\left(C\right):=\frac{\exp\left(-\tau
L\right)}{\mathsf{Tr}\left[\exp\left(-\tau L\right)\right]}$ (70)

is considered, where $\mathsf{Tr}\left(\cdot\right)$ denotes the trace of an
operator. Given a time scale $\tau$, this propagator is known as the
counterpart of the path-integral of general diffusion processes [56].

In each $k$-th iteration, parameter $n_{k}$ is defined by the following way.
Based on the bra-ket formalism, $L^{\left(k\right)}$, the Laplacian of
$G^{\left(k\right)}$, is expressed as

$\displaystyle
L^{\left(k\right)}=\sum_{\lambda}\lambda|\lambda\rangle\langle\lambda|,$ (71)

where each $\lambda$ and $|\lambda\rangle$ denote an eigenvalue and the
associate (column) eigenvector, respectively. Parameter $n_{k}$ is defined as

$\displaystyle n_{k}=\sum_{\lambda}\Theta\left(\frac{1}{\tau}-\lambda\right),$
(72)

where $\Theta\left(\cdot\right)$ is the Heaviside step function. In the
momentum space, the eigenvalues satisfying $\lambda>\frac{1}{\tau}$ correspond
to strong correlations and fast diffusion processes, whose contributions to
$K^{\left(k\right)}$ exhibit a rapid decay because

$\displaystyle
K^{\left(k\right)}=\prod_{\lambda}\exp\left(-\left(\tau\lambda\right)|\lambda\rangle\langle\lambda|\right).$
(73)

According to the non-isomorphic decomposition analysis, this kind of
eigenvalues can be eliminated during coarse graining because of strong
correlations. After coarse graining, the remaining $n_{k}$ eigenvalues
correspond to slow diffusion processes, whose contributions to
$K^{\left(k\right)}$ are long-range and long-term. In the real space, the LRG
simultaneously makes $G^{\left(k\right)}$ coarse grained in a way similar to
the step (B) of our approach. Specifically, nodes are gradually aggregated if
they are connected by the edges that correspond to large values in
$K^{\left(k\right)}$. Each set of aggregated nodes is replaced by a super-
node. This process is repeated until there remain $n_{k}$ nodes (including
super-nodes) in the coarse grained network. In fact, the step (B) of our
approach, which includes $H^{\left(k\right)}$ as a reference, realizes the
above coarse graining in a more easy-to-program way.

In the LRG, the correspondence between the moment space and the real space is
realized based on the renormalized Laplacian $L^{\left(k+1\right)}$

$\displaystyle L^{\left(k+1\right)}=\tau L^{\left(k\right)}_{\mathsf{cg}},$
(74)

where $L^{\left(k\right)}_{\mathsf{cg}}\in\mathbb{R}^{n_{k}\times n_{k}}$ is
the Laplacian of the coarse grained network in the $k$-th iteration. This
Laplacian is expected to be consistent with the coarse graining in the moment
space

$\displaystyle
U^{-1}\left(\sum_{\lambda}\Theta\left(\frac{1}{\tau}-\lambda\right)\lambda|\lambda\rangle\langle\lambda|\right)U$
$\displaystyle=\mathsf{diag}\left(L^{\left(k\right)}_{\mathsf{cg}},\mathbf{0}\right),$
(75)

where $U$ denotes a similarity transformation, notion
$\mathsf{diag}\left(\cdot\right)$ stands for the diagonal, and notion
$\mathbf{0}$ is a all-zero square matrix whose width is $D_{k}-n_{k}$ (here we
mark
$D_{k}=|V^{\left(k\right)}|=\mathsf{dim}\left(L^{\left(k\right)}\right)$). The
left side of Eq. (75) corresponds to the moment space, where all eigenvalues
satisfying $\lambda>\frac{1}{\tau}$ are eliminated, and the right side of Eq.
(75) concerns the real space. Clearly, the similarity transformation exactly
corresponds to the coarse graining on $G^{\left(k\right)}$

$\displaystyle
U=\left(|\nu_{1}\rangle,\ldots,|\nu_{n_{k}}\rangle,|\sigma_{n_{k}+1}\rangle,\ldots,|\sigma_{D_{k}}\rangle\right),$
(76)

where each $\nu_{i}$ is a super-node in the coarse grained network. Notion
$|\nu_{i}\rangle$ is a $D_{k}$-dimensional column ket, whose unitary
components correspond to all the nodes in $G^{\left(k\right)}$ that are
aggregated in $\nu_{i}$ while zero components correspond to all other nodes in
$G^{\left(k\right)}$ that are not covered by $\nu_{i}$. Each ket $\sigma_{j}$
is designed such that
$\left(|\nu_{1}\rangle,\ldots,|\nu_{n_{k}}\rangle,|\sigma_{n_{k}+1}\rangle,\ldots,|\sigma_{D_{k}}\rangle\right)$
defines a group of orthonormal bases. Meanwhile, Eq. (76) is required to
satisfy

$\displaystyle\langle\nu_{i}|\left(\sum_{\lambda}\Theta\left(\frac{1}{\tau}-\lambda\right)\lambda|\lambda\rangle\langle\lambda|\right)|\nu_{j}\rangle$
$\displaystyle=\left(L^{\left(k\right)}_{\mathsf{cg}}\right)_{\nu_{i}\nu_{j}},$
(77)
$\displaystyle\langle\nu_{i}|\left(\sum_{\lambda}\Theta\left(\frac{1}{\tau}-\lambda\right)\lambda|\lambda\rangle\langle\lambda|\right)|\sigma_{j}\rangle$
$\displaystyle=0,$ (78)
$\displaystyle\langle\sigma_{i}|\left(\sum_{\lambda}\Theta\left(\frac{1}{\tau}-\lambda\right)\lambda|\lambda\rangle\langle\lambda|\right)|\sigma_{j}\rangle$
$\displaystyle=0.$ (79)

Given Eqs. (77-79), the validity of Eq. (75) is ensured.

To summarize, the LRG implements a special case of our proposed coarse
graining process described in steps (A-C), where $y\left(\cdot\right)$ serves
as a mapping from the Laplacian to the normalized network propagator and
$n_{k}$ in each $k$-th iteration is calculated according to the number of
modes corresponding to slow diffusion processes.

###  7.3 Manifestation in moment space phenomenological renormalization group

Then we turn to considering the moment space phenomenological renormalization
group (MSPRG) proposed by William Bialek _et al._ [25] (also see discussions
in Ref. [26]).

Following the idea of perturbative theory, the MSPRG considers the Gaussian-
like (i.e., may not be strictly Gaussian) truncated fluctuations of the
probability of the system state $\\{x_{i}\\}$

$\displaystyle
P\left(\\{x_{i}\\}\right)=\frac{1}{Z}\exp\left[-\frac{1}{2}\sum_{i,j}x_{i}K_{ij}x_{j}-\frac{1}{4!}q\sum_{i}x_{i}^{4}\right],$
(80)

where $Z$ serves as a normalization term and $q$ describes the kurtosis (i.e.,
the fluctuations become Gaussian if $q=0$). Matrix $K$ may lack a explicit
expression in most cases. However, once $q\rightarrow 0$ is satisfied, we can
observe an convergence relation between $K$ and the inverse of covariance
matrix $\Sigma$ [25, 26]. Therefore, eliminating the short-term contributions
of large eigenvalues and the associate eigenvectors of $K$ in the moment space
is similar to, but not necessarily equivalent to, eliminating the
contributions of small eigenvalues and the associate eigenvectors of $\Sigma$
(the equivalence holds only when $q=0$ is strictly ensured). This property can
be understood in terms of the correlation function. Specifically, the
eigenvalues of the covariance matrix in the Fourier space is fully
characterized by $L\left(\cdot\right)$, the Fourier transform of the
correlation function

$\displaystyle\sum_{\mathbf{x}_{j}}\Sigma\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)\exp\left(\mathsf{i}\mathbf{k}\mathbf{x}_{j}\right)=\exp\left(\mathsf{i}\mathbf{k}\mathbf{x}_{i}\right)L\left(\mathbf{k}\right),$
(81)

where we have slightly abused the notion to mark
$\Sigma_{ij}=\Sigma\left(\mathbf{x}_{i}-\mathbf{x}_{j}\right)$ in the Fourier
space as $q\rightarrow 0$. According to Eq. (81), we can generally treat the
Fourier transform of the correlation function as the eigenvalue spectrum of
$\Sigma$. Therefore, eliminating the contributions of small eigenvalues and
the associate eigenvectors of $\Sigma$ is similar to eliminating the short-
term or short-range correlation in the Fourier space.

Inspired by these facts, the MSPRG is proposed as the following way. In the
$k$-th iteration, the coarse graining projection $W^{\left(k\right)}$ is
defined as

$\displaystyle
W^{\left(k\right)}=\sum_{i=1}^{n_{k}}|\lambda_{i}\rangle\langle\lambda_{i}|,$
(82)

where $\\{\lambda_{1}\geq\ldots\geq\lambda_{D_{k}}\\}$ denotes the sorted
eigenvalues of $\Sigma^{\left(k\right)}$ (here $\Sigma^{\left(k\right)}$ is
the covariance matrix estimated in the system in the $k$-th iteration) and
$n_{k}$ is simply set as a constant for all iterations in Ref. [25]. Given Eq.
(82), the coarse graining is realized as

$\displaystyle x_{i}^{\left(k\right)}\rightarrow
x_{i}^{\left(k+1\right)}=\sum_{j}W^{\left(k\right)}_{ij}x_{j}^{\left(k\right)}.$
(83)

To summarize, the MSPRG implements a procedure similar to our proposed coarse
graining process in steps (A-C) in the moment space, where
$y\left(\cdot\right)$ is chosen as a mapping from an arbitrary correlation
matrix to the inverse of covariance matrix of the system under the assumption
of translational invariance (i.e., $q=0$). Parameter $n_{k}$ is simply defined
as a constant according to application demands.

###  7.4 Manifestation in real space phenomenological renormalization group

Finally, we consider the real space phenomenological renormalization group
(RSPRG), a real space counterpart of the MSPRG [27, 28]. Compared with the LRG
and the MSPRG, the implementation of RSPRG is rather simple.

Given the system described by Eq. (80), one can define $K$ as the Pearson
correlation matrix

$\displaystyle
K_{ij}=y\left(C_{ij}\right):=\frac{\Sigma_{ij}}{\sqrt{\Sigma_{ii}\Sigma_{jj}}}.$
(84)

If correlation matrix $C$ has been derived following the Pearson correlation,
then $y\left(\cdot\right)$ is nothing more but an identity mapping. The coarse
graining is simply implemented by

$\displaystyle x_{i}^{\left(k\right)}\rightarrow
x_{i^{\prime}}^{\left(k+1\right)}:=x_{i}^{\left(k\right)}+\sum_{j}\delta\left(K_{ij},\max_{r}K_{ir}\right)x_{j}^{\left(k\right)},$
(85)

where $\delta\left(\cdot,\cdot\right)$ is the Kronecker delta function. In the
case where there exists more than one possibilities of $j$ that satisfy
$\delta\left(K_{ij},\max_{r}K_{ir}\right)=1$, the RSPRG only chooses one of
these possibilities while calculating Eq. (85). In other words, parameter
$n_{k}$ is set as $n_{k}=2^{-k}|V^{\left(0\right)}|$ for the $k$-th iteration
(e.g., about half of elements are coarse grained in each iteration).

Therefore, the RSPRG realizes our coarse graining approach in the real space
by designing $y\left(\cdot\right)$ as a mapping from an arbitrary correlation
matrix to the Pearson correlation matrix and defining
$n_{k}=2^{-k}|V^{\left(0\right)}|$.

###  7.5 Non-isomorphic decomposition and renormalization groups

Let us summarize our analyses presented above. We have proposed an abstract
coarse graining approach under the guidance of the non-isomorphic
decomposition analysis. As we have suggested, this general framework has
concrete manifestations in various existing renormalization groups, such as
the LRG [24], the MSPRG [25, 26], and the RSPRG [27, 28].

Our abstract coarse graining process is valid when the internal correlations
inside the investigated system are not weak. Under this condition, each set of
$M$ nodes that need to be aggregated together in step (B) can have a small
value of the automorphism entropy (e.g., a small infinite cluster automorphism
entropy, $A_{\infty}$, or a small value of the first moment of the
automorphism entropy spectrum, $\zeta_{A}$). In other words, for this set of
nodes (i.e., elements), two extreme probability spaces of non-isomorphic
percolation configuration classes, $\mathbb{P}_{\rho}$ and
$\mathbb{P}_{\infty}$, are enough similar to each other. Consequently, the
actual probability space of non-isomorphic percolation configuration classes,
$\mathbb{P}_{\mathsf{ac}}$, must be similar to $\mathbb{P}_{\rho}$ as well
(i.e., according to Eq. (67)). Although in real cases we may never know the
precise expression of $\mathbb{P}_{\mathsf{ac}}$, we can treat
$\mathbb{P}_{\mathsf{ac}}$ as $\mathbb{P}_{\rho}$ because of the small
difference between these two spaces. Therefore, for the concerned set of $M$
elements, we can safely treat them as perfectly and unconditionally coherent
to study them as a whole, even though these elements are only correlated in a
probabilistic manner. The implied information loss, reflected by the
automorphism entropy, is small and acceptable. The way that we analyze these
elements as a whole is exactly the abstract coarse graining process. Various
renormalization groups that realize this abstract coarse graining approach can
help us find the underlying patterns of correlated systems (i.e., potential
criticality at the non-trivial fixed point of renormalization flows).

One may be curious about the consequences if we insist to renormalize the
correlated system when internal correlations are not sufficiently strong,
i.e., we insist to eliminate some freedom degrees in
$\mathbb{P}_{\mathsf{ac}}$ to make $\mathbb{P}_{\mathsf{ac}}$ approach to
$\mathbb{P}_{\rho}$ even though the informational difference,
$H_{\mathsf{ac}}-H_{\rho}$, is not small. In the LRG [24], sufficiently weak
internal correlations (i.e., small values of $|C|$) imply a small trace of the
Laplacian $L^{\left(k\right)}$ according to Eq. (69). Based on the relation
between trace and eigenvalues and the non-negativity of Laplacian eigenvalues,
we know that more Laplacian eigenvalues approach to zero. Consequently, we
derive a larger value of $n_{k}$ in Eq. (72) given a constant $\tau$. In
certain extreme cases, we may derive $n_{k}=|V^{\left(k\right)}|$ such that
the coarse graining process has no effect on the system at all (i.e., there is
no coarse grained element in each iteration). In the MSPRG [25, 26], and the
RSPRG [27, 28], the renormalization flows drive the correlated system with
weak internal correlations towards a Gaussian distribution [25, 26, 27, 28].
The emerged Gaussian distribution at the fixed points of the renormalization
flows, however, are not equivalent to the Gaussian fixed points on the
critical surface when the correlated system is embedded in a space whose
dimension is higher than the upper critical dimension [57]. On the contrary,
this Gaussian distribution is the consequence of the central limit theorem
because system elements are weakly correlated or independent [25, 26, 27, 28].
In a critical system, strongly correlated elements make the central limit
theorem fail, which implies non-Gaussian fixed points.

##  8 Conclusion

In this work, we study the dynamics of the information contained in the
statistical ensemble of percolation configurations during the percolation
process of correlated behaviours. In Secs. 2-5, we develop a framework
referred to as the non-isomorphic decomposition analysis to support a
quantitative analysis of our concerned question. We discover that the
informational difference (i.e., the proposed automorphism entropy and its
variants) between $\mathbb{P}_{\rho}$, a space of the extreme cases where the
correlated system has strictly constrained element relations, and
$\mathbb{P}_{\infty}$, a space of the extreme cases where there exists no
constraint on element relations, decreases as the internal correlations inside
the system become increasingly strong. On the one hand, this property enables
the automorphism entropy and its variants to serve as general indicators of
phase transition phenomena in diverse correlated systems, such as the
absorbing phase transition in branching models [15, 16, 17], spreading
processes in epidemic models [18, 19], and synchronization in oscillator
models [20, 7] shown in Sec. 6. On the other hand, this property also reduces
the informational difference between $\mathbb{P}_{\mathsf{ac}}$, the actual
space where the correlated system lives in, and $\mathbb{P}_{\rho}$. In other
words, this property suggests the possibility to approximate
$\mathbb{P}_{\mathsf{ac}}$ based on $\mathbb{P}_{\rho}$ with small information
loss. In Sec. 7, this possibility is demonstrated as related to
renormalization group theories [21, 22, 23]. We propose an abstract coarse
graining approach based on the non-isomorphic decomposition analysis, which
features concrete manifestations in diverse renormalization groups, including
the Laplacian renormalization group [24], the moment space phenomenological
renormalization group [25, 26], and the real space phenomenological
renormalization group [27, 28]. Our non-isomorphic decomposition analysis is
suggested to convey key information about when and how renormalization groups
can process correlated systems effectively. In sum, this work suggests a new
perspective to uncover the underlying connections between percolation,
correlated systems, and renormalization groups.

As a preliminary exploration, our work certainly has limitations. The present
version of the non-isomorphic decomposition analysis is rooted in the
developed automorphism entropy and lacks other mathematical concepts to enrich
its frameworks. The analytic calculation of the automorphism entropy and its
variants, as we have suggested, may become highly non-trivial until we apply
the mean-field approximation. Therefore, the potential direction for future
explorations may be perfecting the details of the non-isomorphic decomposition
analysis and developing more concepts built on non-isomorphic decomposition,
especially the analytic method to compute the automorphism entropy without the
mean-field approximation. Based on these improvements, we may be enabled to
uniformly analyze percolation, correlated systems, and renormalization groups
in an information-theoretical manner.

## Acknowledgements

This project is supported by the Huawei Innovation Research Program
(TC20221109044), 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 UFR de Mathématiques,
Université de Paris, for the inspiring discussions. Aohua Cheng, who studies
at the Tsien Excellence in Engineering Program, Tsinghua University, is
appreciated for proof reading.

##  Appendix A A practical calculation of the Pólya enumeration theorem

In this section, we present a calculation approach of Eq. (16) parallel to the
generating function method [34, 35]. This approach is more friendly to
computer programming.

Let us represent each permutation $\omega\in S_{M}$ as a product of disjoint
cycles (note that $S_{M}$ denotes the symmetric group). We cluster all
permutations according to their cycle types to obtain all conjugate classes.
The total number of conjugate classes, $\lambda\left(\cdot\right)$, of $S_{N}$
satisfies

$\displaystyle\lambda\left(S_{M}\right)=\frac{1}{|S_{M}|}\sum_{\omega\in
S_{M}}|O\left(\omega\right)|,$ (86)

where $O\left(\cdot\right)$ denotes the conjugate class. The size of each
conjugate class is fully determined by the corresponding cycle type [58].
Assuming $\omega$ is the product of $k\in\mathbb{N}^{+}$ types of cycles,
where the $i$-th ($i\in\\{1,\ldots,k\\}$) type of cycle has an order (i.e.,
length) of $l_{i}$ and occurs $m_{i}$ times, we have

$\displaystyle|O\left(\omega\right)|=\frac{M!}{\prod_{i=1}^{k}m_{i}!l_{i}^{m_{i}}}.$
(87)

For each conjugate class $O\left(\omega\right)$, we can analyze the property
of $\sigma\left(\omega\right)\in S_{|E_{c}|}$, the corresponding permutation
on edges.

First, given every cycle of order $l_{i}$ in $\omega$ (we assume the node set
involved with this cycle is $V_{i}$), there exist
$\lfloor\frac{l_{i}}{2}\rfloor$ associate cycles in
$\sigma\left(\omega\right)$ that consist of the edges formed among nodes in
$V_{i}$. For instance, if there exists a cycle $\left(12345\right)$ in
$\omega$, then there are $\lfloor\frac{5}{2}\rfloor=2$ associate cycles in
$\sigma\left(\omega\right)$, i.e.,
$\left[\left(1,2\right)\left(2,3\right)\left(3,4\right)\left(4,5\right)\right]$
and
$\left[\left(1,3\right)\left(2,4\right)\left(3,5\right)\left(4,1\right)\right]$
(note that notion $\left(a,b\right)$ denotes an edge between nodes $a$ and
$b$).

Second, given two cycles of orders $l_{i}$ and $l_{j}$ in $\omega$ (we assume
the node sets involved with these two cycles are $V_{i}$ and $V_{j}$), there
would be $\frac{l_{i}l_{j}}{\mathsf{lcm}\left(l_{i},l_{j}\right)}$ associate
cycles in $\sigma\left(\omega\right)$ that consist of edges between node sets
$V_{i}$ and $V_{j}$. Here notion $\mathsf{lcm}\left(\cdot,\cdot\right)$
denotes the least common multiple and $\mathsf{lcm}\left(l_{i},l_{j}\right)$
exactly derives the least common period length of these two cycles.

Based on these two properties analyzed above, we have enumerated all the
possible cycles in $\sigma\left(\omega\right)$ that are created by cycles in
$\omega$. Therefore, we know

$\displaystyle
c\left[\sigma\left(\omega\right)\right]=\sum_{i=1}^{k}\lfloor\frac{l_{i}}{2}\rfloor+\sum_{j>i\in\\{1,\ldots,k\\}}\frac{l_{i}l_{j}}{\mathsf{lcm}\left(l_{i},l_{j}\right)}.$
(88)

We can equivalently calculate Eq. (16) as

$\displaystyle|\Gamma\left(M\right)|=\frac{1}{\sum_{O\left(\omega\right)\in\mathbf{O}_{M}}|O\left(\omega\right)|}\sum_{O\left(\omega\right)\in\mathbf{O}_{M}}|O\left(\omega\right)|2^{c\left[\sigma\left(\omega\right)\right]},$
(89)

where $\mathbf{O}_{M}$ is the set of all conjugate classes of group $S_{M}$.
The elements of $\mathbf{O}_{M}$ can be enumerated in a manner of the integer
partition of $M$.

##  Appendix B Computation of automorphism entropy

In this section, we elaborate the details of calculating the non-isomorphic
decomposition analysis in Python. We primarily focus on computing the
automorphism entropy as a function of $M$.

We first enumerate all possible conjugate classes of symmetric group $S_{N}$
to derive $\mathbf{O}_{N}$. The code implementation is developed as

⬇

1import numpy as np

2import scipy.special

3

4def ConjugateClasses(M):

5 # This function enumerate all possible conjugate classes of group S_M

6 # CC is the array of conjugate classes.

7 # For instance, CC= [[4], [3, 1], [2, 2], [2, 1, 1], [1, 1, 1, 1]] given
group S_4

8 # [4] denotes a conjugate class with a form (****), where each element is 1
cycle of length 4

9 # [3,1] denotes a conjugate class with a form (***)(*), where each element
is the product of 1 cycle of length 3 and 1 cycle of length 1

10 # [2, 2] denotes a conjugate class with a form (**)(**), where each element
is the product of 2 cycles of length 2

11 # [2, 1, 1] denotes a conjugate class with a form (**)(*)(*), where each
element is the product of 1 ccle of length 2 and 2 cycles of length 1

12 Sol = Solution()

13 CC = Sol.UniquePartitions(M)

14 return CC

The enumeration is equivalent to integer partition in modern number theory,
which can be derived by

⬇

1class Solution:

2 # This function enumerate all possible partitions of a given integer

3 def __init__(self):

4 self.temp = []

5

6 def Solve(self, A, V, ID, S, N):

7 if S == N:

8 V.append(self.temp.copy())

9 return

10 if ID < 0:

11 return

12 self.Solve(A, V, ID-1, S, N)

13 if S < N:

14 self.temp.append(A[ID])

15 self.Solve(A, V, ID, S+A[ID], N)

16 self.temp.pop()

17

18 def UniquePartitions(self, N):

19 A = [ID for ID in range(1, N+1)]

20 V = []

21 self.Solve(A, V, N-1, 0, N)

22 V.reverse()

23 return V

Given each conjugate class of symmetric group $S_{M}$, we can calculate its
size based on the following code.

⬇

1def Size_ConjugateClasses(CC,M):

2 # This function calculates the size of each conjugate class of group S_M

3 # Each element in CCSize is the size of the corresponding conjugate class in
CC

4 CCSize=np.zeros(len(CC))

5 for ID in range(len(CC)):

6 CClass=np.array(CC[ID])

7 Uni,Count=np.unique(CClass, return_counts=True)

8 Denominator=np.prod(scipy.special.factorial(Count)*np.power(Uni,Count))

9 CCSize[ID]=np.math.factorial(M)/Denominator

10 return CCSize

Given each conjugate class $O\left(\omega\right)$, we can count the number of
all the possible cycles in $\sigma\left(\omega\right)$ that are created by
cycles in $\omega$

⬇

1def Num_Cycles(CC):

2 # This function calculates the number of cycles in the permutation
\sigma(\omega)

3 # Each element in CNum is the number of cycles in the corresponding
conjugate class in CC

4 CNum=np.zeros(len(CC))

5 for ID in range(len(CC)):

6 CClass=np.array(CC[ID])

7 Term1=np.sum(np.floor(CClass/2))

8 Term2=0

9 for ID1 in range(len(CClass)):

10 for ID2 in range(ID1+1,len(CClass)):

11 Term2=Term2+np.gcd(CClass[ID1], CClass[ID2])

12 CNum[ID]=Term1+Term2

13 return CNum

Based on these functions, we can calculate $|\Gamma\left(M\right)|$ with the
following code.

⬇

1def Gamma_M(CCSize,CNum):

2 ## This function calculates |Gamma(M)|

3 PartitionF=np.sum(CCSize)

4 Numerator=np.sum(CCSize*np.power(2,CNum))

5 GM=Numerator/PartitionF

6 return GM

Meanwhile, probability measure $P_{\infty}\left(\cdot\right)$ can be readily
calculated by the following function.

⬇

1def Prob_Inf(M):

2 ## This function calculates the probability measure P_inf

3 CC=ConjugateClasses(M)

4 CCSize=Size_ConjugateClasses(CC,M)

5 CNum=Num_Cycles(CC)

6 GM=Gamma_M(CCSize,CNum)

7 ProbI=1/GM

8 return GM,ProbI

After obtaining $|\Gamma\left(M\right)|$, we turn to deriving
$|\Gamma\left(M\right)|$ based on the inverse Euler transform. In our code
implementation, several necessary functions are proposed to support the
transform. For a given integer, we can enumerate all its divisors

⬇

1def AllDivisor(N):

2 # This function enumerate all possible divisors of a given integer

3 Divisors = set()

4 for ID in range(1, int(N**0.5)+1):

5 if N % ID == 0:

6 Divisors.add(ID)

7 Divisors.add(N//ID)

8 return list(Divisors)

and determine if it is a prime number.

⬇

1def isPrime(N) :

2 # This function determines if a given integer is a prime

3 if (N < 2) :

4 return False

5 for ID in range(2, N + 1) :

6 if (ID * ID <= N and N % ID == 0) :

7 return False

8 return True

Based on these two functions, we can computationally implement the Möbius
function.

⬇

1def Mobius(N) :

2 # This function determines if a given integer is a prime

3 if (N == 1) :

4 return 1

5 P = 0

6 for ID in range(1, N + 1) :

7 if (N % ID == 0 and

8 isPrime(ID)) :

9 if (N % (ID * ID) == 0) :

10 return 0

11 else :

12 P = P + 1

13 if(P % 2 != 0) :

14 return -1

15 else :

16 return 1

These functions enable us to define the inverse Euler transform of the
sequence of
$\\{|\Gamma\left(1\right)|,|\Gamma\left(2\right)|,|\Gamma\left(3\right)|,\ldots\\}$
to derive
$\\{|\Gamma^{*}\left(1\right)|,|\Gamma^{*}\left(2\right)|,|\Gamma^{*}\left(3\right)|,\ldots\\}$.

⬇

1def Prob_Rho(M):

2 ## This function calculates the probability measure P_rho

3 ## GMVec is the sequence of |Gamma(0)|, |Gamma(1)|, |Gamma(2)|, ......,
|Gamma(M)|

4 ## ProbIVec is the sequence of P_inf(0), P_inf(1), P_inf(2), ......,
P_inf(M)

5 ## GStarMVec is the sequence of |Gamma*(0)|, |Gamma*(1)|, |Gamma*(2)|,
......, |Gamma*(M)|

6 ## ProbRhoVec is the sequence of P_rho (0), P_rho (1), P_rho (2), ......,
P_rho(M)

7 ## LambdaVec is the sequence of Lambda(1), Lambda(2), ......, Lambda(M)

8 GMVec=np.zeros(M+1)

9 ProbIVec=np.zeros(M+1)

10 for ID in range(M+1):

11 print(’Processing P_inf: %f’ % (ID/(M+1)))

12 GM,ProbI=Prob_Inf(ID)

13 GMVec[ID]=GM

14 ProbIVec[ID]=ProbI

15

16 LambdaVec=np.zeros(M+1)

17 LambdaVec[1]=GMVec[1]

18 for ID in range(2,len(LambdaVec)):

19 print(’Processing Lambda: %f’ % (ID/(M+1)))

20 SumTerm=0

21 for ID2 in range(1,ID):

22 SumTerm=SumTerm+LambdaVec[ID2]*GMVec[ID-ID2]

23 LambdaVec[ID]=ID*GMVec[ID]-SumTerm

24

25 GStarMVec=np.zeros(M+1)

26 ProbRhoVec=np.zeros(M+1)

27 GStarMVec[0]=GMVec[0]

28 ProbRhoVec[0]=ProbIVec[0]

29 for ID in range(1,M+1):

30 print(’Processing P_rho: %f’ % (ID/(M+1)))

31 Divisors=AllDivisor(ID)

32 MobiusResult=0

33 for Divisor in Divisors:

34 MobiusResult=MobiusResult+Mobius(int(ID/Divisor))*LambdaVec[Divisor]

35 GStarMVec[ID]=MobiusResult/ID

36 ProbRhoVec[ID]=ID/MobiusResult

37 return GMVec,ProbIVec,GStarMVec,ProbRhoVec,LambdaVec

Finally, we calculate the automorphism entropy as a function of
$i\in\\{0,\dots,M\\}$.

⬇

1def AutomorphismEntropyFunction(M):

2 ## GMVec is the sequence of |Gamma(0)|, |Gamma(1)|, |Gamma(2)|, ......,
|Gamma(M)|

3 ## ProbIVec is the sequence of P_inf(0), P_inf(1), P_inf(2), ......,
P_inf(M)

4 ## GStarMVec is the sequence of |Gamma*(0)|, |Gamma*(1)|, |Gamma*(2)|,
......, |Gamma*(M)|

5 ## ProbRhoVec is the sequence of P_rho (0), P_rho (1), P_rho (2), ......,
P_rho(M)

6 ## LambdaVec is the sequence of Lambda(1), Lambda(2), ......, Lambda(M)

7 ## SEIVec is the sequence of Shannon entropy H_inf(0), H_inf(1), H_inf(2),
......, H_inf(M)

8 ## SERhoVec is the sequence of Shannon entropy H_rho(0), H_rho(1), H_rho(2),
......, H_rho(M)

9 ## AEVector is the sequence of Automorphism entropy A(0), A(1), A(2),
......, A(M)

10 ## Please note that we define A(0)=A(1)=1 in this work, yet other
definitions are acceptable

11 GMVec,ProbIVec,GStarMVec,ProbRhoVec,LambdaVec=Prob_Rho(M)

12 SEIVec=ShannonEntropyFunction(ProbIVec)

13 SERhoVec=ShannonEntropyFunction(ProbRhoVec)

14 AEVector=1-SERhoVec/SEIVec

15 AEVector[0]=1

16 if len(AEVector)>1:

17 AEVector[1]=1

18 return
GMVec,ProbIVec,GStarMVec,ProbRhoVec,LambdaVec,SEIVec,SERhoVec,AEVector

The Shannon entropy terms used in the above function are calculated based on
the following codes.

⬇

1def ShannonEntropy(P):

2 ## This function calculates the Shannon entropy of a given probability
measure

3 S=1/P*(-1*P*np.log(P))

4 return S

5

6def ShannonEntropyFunction(ProbVec):

7 ## This function calculates the Shannon entropy as a function of M

8 SEVec=np.zeros(len(ProbVec))

9 for ID in range(len(ProbVec)):

10 SEVec[ID]=ShannonEntropy(ProbVec[ID])

11 return SEVec

To run the above program in a unified manner, one can consider the following
instance where $M=100$.

⬇

1M=100

2GNVec,ProbIVec,GStarNVec,ProbRhoVec,LambdaVec,SEIVec,SERhoVec,AEVector=AutomorphismEntropyFunction(M)

##  Appendix C Empirical estimation of automorphism entropy

Based on the estimated coefficients of exponential decay in Eq. (25), we can
estimate the automorphism entropy based on the following code.

⬇

1def AutomorphismEntropyEstimation(M):

2 ## This function estimates the exponential decay of A(M)

3 ## AE is the estimated value of A(M)

4 AE=3.956*np.exp(-0.6876*M)

5 if M<2:

6 AE=1

7 return AE

##  Appendix D Computation of the Bernoulli percolation

In this section, we propose the Python functions for implementing the
Bernoulli percolation.

Let us begin with network generation and dilution. Given a target network size
and an occupation probability $\rho$, we first generate a complete graph of
the corresponding size and randomly delete each edge according to a
probability of $1-\rho$.

⬇

1import numpy as np

2import networkx as nx

3

4def NetworkGeneration(Size,rho):

5 ## Size is the size of the initial network

6 ## rho is the occupation probability

7 ## AMatrix is the adjacency matrix

8 ## UniDeg is the vector of unique values of node degree

9 ## ProbDegree is the vector of the probability of each unique value of node
degree

10 ## MeanDeg is the average node degree

11 G = nx.complete_graph(Size)

12 AMatrix=nx.to_numpy_matrix(G)

13 for ID in range(Size):

14 for ID1 in range(ID,Size):

15 AMatrix[ID,ID1]=np.random.binomial(1,rho,1)

16 AMatrix[ID1,ID]=AMatrix[ID,ID1]

17 DegreeV=np.sum(np.asarray(AMatrix), axis=0)

18 UniDeg, Freq=np.unique(DegreeV, return_counts=True)

19 ProbDegree=Freq/Size

20 MeanDeg=np.mean(DegreeV)

21 return UniDeg,ProbDegree,MeanDeg

Then we calculate the order parameter $\rho_{n}$ of the Bernoulli percolation
as a function of $\rho$.

⬇

1def BernoulliPercolationProcess(Size):

2 ## Size is the size of the initial network

3 ## Rho_eVector is the value of Rho_e defined as a function of
rho=0,0.01,0.02,......1

4 ## Rho_nVector is the value of Rho_n defined as a function of
rho=0,0.01,0.02,......1

5 ## Note that Rho_n is the order parameter of the Bernoulli percolation

6 rhoVector=np.arange(0.005,1.005,0.005)

7 Rho_eVector=np.zeros(len(rhoVector))

8 Rho_nVector=np.zeros(len(rhoVector))

9 for ID in range(len(rhoVector)):

10 rho=rhoVector[ID]

11 UniDeg,ProbDegree,MeanDeg=NetworkGeneration(Size,rho)

12 print(’Processing Percolation: %f’ % (ID/len(rhoVector)))

13 Rho_e_Vec=np.arange(0.005,1.005,0.005)

14 Rho_e_Eq=np.zeros(len(rhoVector))

15 for ID1 in range(len(Rho_e_Vec)):

16 PRho_e=Rho_e_Vec[ID1]

17
Rho_e_Eq[ID1]=1-np.sum(ProbDegree*UniDeg/MeanDeg*np.power(1-rho*PRho_e,UniDeg-1))

18 Rho_e=Rho_e_Vec[np.argmin(np.abs(Rho_e_Eq-Rho_e_Vec))]

19 Rho_n=1-np.sum(ProbDegree*np.power(1-rho*Rho_e,UniDeg))

20 Rho_eVector[ID]=Rho_e

21 Rho_nVector[ID]=Rho_n

22 return Rho_eVector,Rho_nVector

A sample instance of the Bernoulli percolation process in a network with
$1000$ nodes can be defined by the following code.

⬇

1Size=1000

2Rho_eVector,Rho_nVector=BernoulliPercolationProcess(Size)

##  Appendix E The probability for a node to belong to a cluster

In this section, we present the Python function to calculate
$\pi_{r}\left(\rho\right)$, the probability for a randomly selected node to
belong to a cluster whose size is $r$. The function is defined as the
following.

⬇

1def PiDistribution(Size,rho,Type):

2 ## Size is the size of the initial network

3 ## rho is the occupation probability

4 ## Type=0 or Type=1 indicates the method used to calculate the pi
distribution

5 ## Type=0 use the method in Eq.(31)

6 ## Type=1 use the method in Eq.(34)

7 ## PiProbSize is a vector of probability for a node to belong to a cluster
of a given size

8 ## The i-th element in ProbSize is the probability for a node to belong to a
cluster whose size is i

9 if Type==0:

10 ClusterSize=np.arange(1,Size,1)

11 PiProbSize=np.zeros(len(ClusterSize))

12 for ID in range(len(ClusterSize)):

13 CSize=ClusterSize[ID]

14 if CSize==1:

15 PiProbSize[ID]=np.power(1-rho,Size-1)

16 elif CSize==2:

17 PiProbSize[ID]=rho*(Size-1)*np.power(1-rho,2*Size-4)

18 elif CSize>2:

19 Term1=np.power(rho,CSize-1)*(Size-1)/scipy.special.factorial(CSize-1)

20 Term2=np.prod(CSize*(Size-2)-np.arange(0,CSize-2,1))

21 Term3=np.power(1-rho,CSize*(Size-3)+2)

22 PiProbSize[ID]=Term1*Term2*Term3

23 elif Type==1:

24 ClusterSize=np.arange(1,Size,1)

25 PiProbSize=np.zeros(len(ClusterSize))

26 for ID in range(len(ClusterSize)):

27 CSize=ClusterSize[ID]

28 Mean=(Size-1)*rho

29
PiProbSize[ID]=1/scipy.special.factorial(CSize)*np.exp(-1*Mean*CSize)*np.power(Mean*CSize,CSize-1)

30 PiProbSize=np.nan_to_num(PiProbSize, nan=0, posinf=0)

31 return PiProbSize

##  Appendix F Computation of automorphism entropy spectrum

In this section, we elaborate the details of calculating the automorphism
entropy spectrum.

We begin with deriving the probability distribution of cluster size. According
to Eq. (35), the derivation is equivalent to solving the following linear
equations

$\displaystyle\left[\begin{matrix}\pi_{1}\left(\rho\right)-1&2\pi_{1}\left(\rho\right)&\ldots&|V|\pi_{1}\left(\rho\right)\\\
\pi_{2}\left(\rho\right)&2\pi_{2}\left(\rho\right)-2&\ldots&|V|\pi_{2}\left(\rho\right)\\\
\vdots&\vdots&\vdots&\vdots\\\
\pi_{|V|}\left(\rho\right)&2\pi_{|V|}\left(\rho\right)&\ldots&|V|\pi_{|V|}\left(\rho\right)-|V|\end{matrix}\right]\left[\begin{matrix}\theta_{1}\left(\rho\right)\\\
\theta_{2}\left(\rho\right)\\\ \vdots\\\
\theta_{|V|}\left(\rho\right)\end{matrix}\right]=\left[\begin{matrix}0\\\ 0\\\
\vdots\\\ 0\end{matrix}\right],$ (90)

which are subject to $\sum_{r}\theta_{r}\left(\rho\right)=1$. Unfortunately,
the first matrix in Eq. (90) is not row full rank (i.e., one can readily find
linear correlations among rows after adding the first row with all the rest
rows). To find a non-trivial solution, we need to replace the first rows of
the first and the third matrices in Eq. (90) by the following rows

$\displaystyle\left[\begin{matrix}1&1&\ldots&1\\\
\pi_{2}\left(\rho\right)&2\pi_{2}\left(\rho\right)-2&\ldots&|V|\pi_{2}\left(\rho\right)\\\
\vdots&\vdots&\vdots&\vdots\\\
\pi_{|V|}\left(\rho\right)&2\pi_{|V|}\left(\rho\right)&\ldots&|V|\pi_{|V|}\left(\rho\right)-|V|\end{matrix}\right]\left[\begin{matrix}\theta_{1}\left(\rho\right)\\\
\theta_{2}\left(\rho\right)\\\ \vdots\\\
\theta_{|V|}\left(\rho\right)\end{matrix}\right]=\left[\begin{matrix}1\\\ 0\\\
\vdots\\\ 0\end{matrix}\right].$ (91)

Solving Eq. (91), we can derive

$\displaystyle\theta_{1}\left(\rho\right)$
$\displaystyle=\frac{1}{\sum_{r}\frac{\pi_{r}\left(\rho\right)}{r\pi_{1}\left(\rho\right)}},$
(92) $\displaystyle=\left(1+\sum_{k=2}^{|V|}\frac{\exp\left(\rho
k-\rho|V|k\right)\left(\rho|V|k-\rho
k\right)^{k-1}}{k^{2}\left(k-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{-1},$
(93)

and

$\displaystyle\theta_{r}\left(\rho\right)$
$\displaystyle=\frac{\frac{\pi_{r}\left(\rho\right)}{r\pi_{1}\left(\rho\right)}}{\sum_{r}\frac{\pi_{r}\left(\rho\right)}{r\pi_{1}\left(\rho\right)}},$
(94) $\displaystyle=\frac{\exp\left(\rho r-\rho|V|r\right)\left(\rho|V|r-\rho
r\right)^{r-1}}{r^{2}\left(r-1\right)!\left(1-\rho\right)^{|V|-1}}\left(1+\sum_{k=2}^{|V|}\frac{\exp\left(\rho
k-\rho|V|k\right)\left(\rho|V|k-\rho
k\right)^{k-1}}{k^{2}\left(k-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{-1},$
(95)

where $r\geq 2$ in Eq. (95). For convenience, we use the Kronecker delta
$\delta\left(\cdot,\cdot\right)$ to reorganize Eqs. (92-95) as

$\displaystyle\theta_{r}\left(\rho\right)$ $\displaystyle=$
$\displaystyle\left(\frac{\exp\left(\rho r-\rho|V|r\right)\left(\rho|V|r-\rho
r\right)^{r-1}}{r^{2}\left(r-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{1-\delta\left(1,r\right)}\left(1+\sum_{k=2}^{|V|}\frac{\exp\left(\rho
k-\rho|V|k\right)\left(\rho|V|k-\rho
k\right)^{k-1}}{k^{2}\left(k-1\right)!\left(1-\rho\right)^{|V|-1}}\right)^{-1},$
(96)

where $r\geq 1$.

Based on Eq. (96), a Python function can be proposed to calculate the cluster
size distribution.

⬇

1def ThetaDistribution(Size,rho):

2 ## Size is the size of the initial network

3 ## rho is the occupation probability

4 ## ThetaProbSize is a vector of probability for a cluster with a given size
to occur in the network (i.e., the cluster size distribution)

5 ## The i-th element in ThetaSize is the probability for a cluster with a
size of i to occur

6

7 ClusterSize=np.arange(1,Size,1)

8 ThetaProbSize=np.zeros(len(ClusterSize))

9 for ID in range(len(ClusterSize)):

10 CSize=ClusterSize[ID]

11 if CSize==1:

12
Term0=np.exp(rho*np.arange(2,Size,1)-rho*Size*np.arange(2,Size,1))*np.power((rho*Size*np.arange(2,Size,1)-rho*np.arange(2,Size,1)),np.arange(1,Size-1,1))

13
Term1=np.power(np.arange(2,Size,1),2)*scipy.special.factorial(np.arange(1,Size-1,1))*np.power(1-rho,Size-1)

14 ThetaProbSize[ID]=1/(1+np.sum(Term0/Term1))

15 elif CSize>1:

16
Term0=np.exp(rho*np.arange(2,Size,1)-rho*Size*np.arange(2,Size,1))*np.power((rho*Size*np.arange(2,Size,1)-rho*np.arange(2,Size,1)),np.arange(1,Size-1,1))

17
Term1=np.power(np.arange(2,Size,1),2)*scipy.special.factorial(np.arange(1,Size-1,1))*np.power(1-rho,Size-1)

18 Term2=np.exp(rho*CSize-rho*Size*CSize)*np.power((rho*Size*CSize-
rho*CSize),CSize-1)

19
Term3=np.power(CSize,2)*scipy.special.factorial(CSize-1)*np.power(1-rho,Size-1)

20 ThetaProbSize[ID]=Term2/Term3*1/(1+np.sum(Term0/Term1))

21 ThetaProbSize=np.nan_to_num(ThetaProbSize, nan=0, posinf=0)

22 return ThetaProbSize

Finally, we can calculate the automorphism entropy spectrum based on Eqs.
(37-38). A Python function is proposed to realize this calculation.

⬇

1def AutomorphismEntropySpectrum(Size,rho,Type):

2 ## Size is the size of the initial network

3 ## rho is the occupation probability

4 ## Type=0 or Type=1 indicates the method used to calculate the pi
distribution

5 ## Type=0 use the method in Eq.(31)

6 ## Type=1 use the method in Eq.(34)

7 ## AEValue is a vector of automorphism entropies calculated in all possible
clusters (i.e., clusters whose sizes are 1,.....,Size)

8 ## AEProbability is a vector of the probability of each automorphism entropy

9 PiProbSize=PiDistribution(Size,rho,Type)

10 ThetaProbSize=ThetaDistribution(Size,rho)

11 ClusterSize=np.arange(1,Size,1)

12 AEValue=np.zeros(len(ClusterSize))

13 for ID in range(len(ClusterSize)):

14 CSize=ClusterSize[ID]

15 AEValue[ID]=AutomorphismEntropyEstimation(CSize)

16 ## Note that AEValue[1]=AEValue[2]=1 (i.e., the automorphism entropy
satisfies A(1)=A(2)=1)

17 ## Therefore, we need to deal with AEValue[1:2] and AEProbability[1:2] to
ensure that the automorphism entropy spectrum is a probability distribution

18 AEValue=np.delete(AEValue, 0)

19 AEProbability=np.delete(ThetaProbSize, 0)

20 AEProbability[1]=ThetaProbSize[1]+ThetaProbSize[0]

21 return PiProbSize,ThetaProbSize,AEValue,AEProbability

##  Appendix G Computation of directed percolation in branching model

In our experiment, a branching model with $100$ elements and
$p\in\left[0,0.4\right]$ is proposed. Given each value of
$p\in\left[0,0.4\right]$, the experiment runs for $10^{6}$ time steps to
generate the data. Under each condition of $p$, we can implement an
independent analysis of the power-law behaviours, scaling relation, and
universal collapse shape to verify if Eqs. (40-45) hold.

Figure 11:  The power-law behaviours, scaling relation, and universal collapse
shape analyzed under different conditions of $p$. (a) shows the plane of
$\left(\frac{\alpha-1}{\beta-1},\gamma\right)$, where each data point
corresponds to a specific value of $p$. The inserted sub-plot visualize the
absolute difference between $\frac{\alpha-1}{\beta-1}$ and $\gamma$ as a
function of $p$, where a low-difference interval (i.e.,
$|\gamma-\frac{\alpha-1}{\beta-1}|\leq 0.1$) is marked by a blue region. (b)
illustrates $|\gamma-1-\psi|$ and $|\frac{\alpha-1}{\beta-1}-1-\psi|$ as the
functions of $p$.

In Fig. 11, we first verify if the scaling relation in Eq. (43) holds. We
illustrate the plane of $\left(\frac{\alpha-1}{\beta-1},\gamma\right)$ to see
whether $\left(\frac{\alpha-1}{\beta-1},\gamma\right)$ is close to the line of
$y=x$. As shown in our result, although an increasing value of
$p\in\left[0,0.4\right]$ makes $\left(\frac{\alpha-1}{\beta-1},\gamma\right)$
approach to the line of $y=x$, it is hard to find a piece of strict evidence
to demonstrate $p=0.25$ as the critical point. Meanwhile, we also measure
$|\gamma-1-\psi|$ and $|\frac{\alpha-1}{\beta-1}-1-\psi|$ as the functions of
$p$. When the power-law behaviours, scaling relation, and universal collapse
shape of the system are strictly satisfied, we expect to see small values of
these functions. As shown in Fig. 11, although the increase of $p$ reduces
$|\gamma-1-\psi|$ and $|\frac{\alpha-1}{\beta-1}-1-\psi|$ gradually, data
noise significantly fuzzifies the variation trends of these functions.
Consequently, one may not be able to precisely confirm the critical point
based on Fig. 11, especially when the _a priori_ knowledge about the branching
model is limited.

## References

  * [1] Wilson H R and Cowan J D 1973 Kybernetik 13 55–80
  * [2] Amari S i 1977 Biological cybernetics 27 77–87
  * [3] Bressloff P C 2010 Physical Review E 82 051903
  * [4] Bressloff P C 2010 SIAM Journal on Applied Mathematics 70 1488–1521
  * [5] Binder K and Young A P 1986 Reviews of Modern physics 58 801
  * [6] Grimmett G R 2006 The random-cluster model vol 333 (Springer Science & Business Media)
  * [7] Acebrón J A, Bonilla L L, Vicente C J P, Ritort F and Spigler R 2005 Reviews of modern physics 77 137
  * [8] Scianna M and Preziosi L 2012 Multiscale Modeling & Simulation 10 342–382
  * [9] Chen N, Glazier J A, Izaguirre J A and Alber M S 2007 Computer physics communications 176 670–681
  * [10] Chiang M and Marenduzzo D 2016 EPL (Europhysics Letters) 116 28009
  * [11] Li M, Liu R R, Lü L, Hu M B, Xu S and Zhang Y C 2021 Physics Reports 907 1–68
  * [12] Duminil-Copin H 2018 Sixty years of percolation Proceedings of the International Congress of Mathematicians: Rio de Janeiro 2018 (World Scientific) pp 2829–2856
  * [13] Georgakopoulos A and Panagiotis C 2018 arXiv preprint arXiv:1811.07404
  * [14] Hamermesh M 2005 Mathematical Tools for Physicists 189–212
  * [15] Hinrichsen H 2000 Advances in physics 49 815–958
  * [16] Lübeck S 2004 International Journal of Modern Physics B 18 3977–4118
  * [17] Henkel M, Hinrichsen H, Lübeck S and Pleimling M 2008 Non-equilibrium phase transitions vol 1 (Springer)
  * [18] Moreno Y, Pastor-Satorras R and Vespignani A 2002 The European Physical Journal B-Condensed Matter and Complex Systems 26 521–529
  * [19] Newman M E 2002 Physical review E 66 016128
  * [20] Rodrigues F A, Peron T K D, Ji P and Kurths J 2016 Physics Reports 610 1–98
  * [21] Zinn-Justin J 2007 Phase transitions and renormalization group (Oxford University Press on Demand)
  * [22] Wilson K G 1983 Reviews of Modern Physics 55 583
  * [23] Goldenfeld N 2018 Lectures on phase transitions and the renormalization group (CRC Press)
  * [24] Villegas P, Gili T, Caldarelli G and Gabrielli A 2023 Nature Physics 1–6
  * [25] Bradde S and Bialek W 2017 Journal of statistical physics 167 462–475
  * [26] Lahoche V, Samary D O and Tamaazousti M 2022 Journal of Statistical Mechanics: Theory and Experiment 2022 033101
  * [27] Meshulam L, Gauthier J L, Brody C D, Tank D W and Bialek W 2018 arXiv preprint arXiv:1812.11904
  * [28] Meshulam L, Gauthier J L, Brody C D, Tank D W and Bialek W 2019 Physical review letters 123 178103
  * [29] Newman M E 2007 Physical review e 76 045101
  * [30] Kryven I 2017 Physical Review E 95 052303
  * [31] Kryven I 2017 Physical Review E 96 052304
  * [32] Moore C and Newman M E 2000 Physical Review E 62 7059
  * [33] Newman M E, Strogatz S H and Watts D J 2001 Physical review E 64 026118
  * [34] Harary F and Palmer E M 2014 Graphical enumeration (Elsevier)
  * [35] Bóna M 2015 Handbook of enumerative combinatorics vol 87 (CRC Press)
  * [36] Sloane N J et al. 2003 The on-line encyclopedia of integer sequences
  * [37] Harary F and Read R C 2006 Is the null-graph a pointless concept? Graphs and Combinatorics: Proceedings of the Capital Conference on Graph Theory and Combinatorics at the George Washington University June 18–22, 1973 (Springer) pp 37–44
  * [38] Haldeman C and Beggs J M 2005 Physical review letters 94 058101
  * [39] Williams-García R V, Moore M, Beggs J M and Ortiz G 2014 Physical Review E 90 062714
  * [40] Fosque L J, Williams-García R V, Beggs J M and Ortiz G 2021 Physical Review Letters 126 098101
  * [41] Marshall N, Timme N M, Bennett N, Ripp M, Lautzenhiser E and Beggs J M 2016 Frontiers in physiology 7 250
  * [42] Dalla Porta L and Copelli M 2019 PLoS computational biology 15 e1006924
  * [43] Tian Y, Tan Z, Hou H, Li G, Cheng A, Qiu Y, Weng K, Chen C and Sun P 2022 Network Neuroscience 6 1148–1185
  * [44] Sood V, Mathieu M, Shreim A, Grassberger P and Paczuski M 2010 Physical review letters 105 178701
  * [45] Ermakova E, Makhmutova P and Yarovaya E 2019 Stochastic Models 35 300–317
  * [46] Panaretos V M 2007 Journal of Mathematical Biology 54 645–668
  * [47] Kovchegov Y, Zaliapin I and Ben-Zion Y 2022 Geophysical Journal International 231 567–583
  * [48] Turcotte D, Abaimov S, Dobson I and Rundle J 2009 Physical Review E 79 016101
  * [49] Girardi-Schappo M 2021 Journal of Physics: Complexity 2 031003
  * [50] Laurson L, Illa X, Santucci S, Tore Tallakstad K, Måløy K J and Alava M J 2013 Nature communications 4 2927
  * [51] Athreya K B, Ney P E and Ney P 2004 Branching processes (Courier Corporation)
  * [52] Boccaletti S, Kurths J, Osipov G, Valladares D and Zhou C 2002 Physics reports 366 1–101
  * [53] Arenas A, Díaz-Guilera A, Kurths J, Moreno Y and Zhou C 2008 Physics reports 469 93–153
  * [54] Boccaletti S, Almendral J, Guan S, Leyva I, Liu Z, Sendiña-Nadal I, Wang Z and Zou Y 2016 Physics Reports 660 1–94
  * [55] Ghosh D, Frasca M, Rizzo A, Majhi S, Rakshit S, Alfaro-Bittner K and Boccaletti S 2022 Physics Reports 949 1–63
  * [56] Feynman R P, Hibbs A R and Styer D F 2010 Quantum mechanics and path integrals (Courier Corporation)
  * [57] Jona-Lasinio G 2001 Physics Reports 352 439–458
  * [58] Scott W R 2012 Group theory (Courier Corporation)