Self-Organised Critical Dynamics as a Key to Fundamental Features of Complexity in Physical, Biological, and Social Networks

: Studies of many complex systems have revealed new collective behaviours that emerge through the mechanisms of self-organised critical ﬂuctuations. Subject to the external and endogenous driving forces, these collective states with long-range spatial and temporal correlations often arise from the intrinsic dynamics with the threshold nonlinearity and geometry-conditioned interactions. The self-similarity of critical ﬂuctuations enables us to describe the system using fewer parameters and universal functions that, on the other hand, can simplify the computational and information complexity. Currently, the cutting-edge research on self-organised critical systems across the scales strives to formulate a unifying mathematical framework, utilise the critical universal properties in information theory, and decipher the role of hidden geometry. As a prominent example, we study the ﬁeld-driven spin dynamics on the hysteresis loop in a network with higher-order structures described by simplicial complexes, which provides a geometric-frustration environment. While providing motivational illustrations from physical, biological, and social systems, along with their networks, we also demonstrate how the self-organised criticality occurs at the interplay of the complex topology and driving mode. This study opens up new promising routes with powerful tools to address a long-standing challenge in the theory and applications of complexity science ingrained in the efﬁcient analysis of self-organised critical states under the competing higher-order interactions embedded in complex geometries.

In this context, the sandpile automata (SPA) models serve as the paradigm of SOC, marked by the formulation of BTW sandpile automaton [3]. They have been a proper mathematical framework to define different quantities and characterise the threshold instability and the character of the emergent critical states [7]. In addition to the deterministic symmetrical BTW, other SPA models with probabilistic rules [8], directed deterministic [7] or stochastic [9] SPA all introduced new ways to SOC and different universal features (classes of criticality) [10].
Since the beginning, the idea of SOC was linked to complexity and contingency [11] in the sense of the emergent new feature through the collective dynamics. Hence, SOC provides an underlying mechanism for the physical complexity-the genesis of new collective properties from individual constituents, which is found in many complex systems across the scales (see more details below). It has been argued that the same physical principles lay the foundations for biological complexity [12]. Remarkably, the frustration caused by competing interactions at different levels could be identified as a "general driving force behind the emergence of complexity" [12].
Another link connects the information complexity [13] and SOC. In the context of data analysis, it concerns the relevance [14] of the information that can be stored in the critical state, as well as the related computational complexity [15,16], the problems of prediction of individual states [17], and mesoscopic avalanche prediction [18]. The idea that complex systems exhibit a hidden order, emerging through self-organised dynamics, links it to information that that ordering patterns may contain [19]. This idea becomes increasingly more interesting in the era of quantum information, which substantially extended the views of classical information theory introducing new concepts linked to complexity, see the main focus point in [20]. In this context, the occurrence of critical states with their self-similarity leads to the optimisation of information entropy [21], thus reducing the information complexity. Remarkably, considering the information relevance, it was shown [21] that the occurrence of a power-law frequency distribution (such as Zipf's law) "arises at the optimal trade-off between resolution (i.e., compression) and relevance". Such "maximally informative" samples, on the other hand, allow assessment of the number of parameters that can be inferred from a given dataset without prior knowledge about the underlying generative processes.
Furthermore, the information theory analysis also revealed a formal similarity of the SPA automata and several other systems with a threshold-like nonlinearity, in particular, the zero-temperature Ising spin dynamics [22]. Recent studies considered SPA dynamics on non-compact lattices and graphs [23], fractals and percolation clusters [24][25][26] and Penrose tailing [27], adapting the dynamical rules to the local connectivity of sites (nodes of the graph). It has been understood that the substrate structure may profoundly impact the dynamics yielding the stable critical states of different universality classes that can be reached. Moreover, the abelian nature of the dynamics (independence of the order of toppling) was one of the key features leading to the observed critical states [23]. Another line of research uses continuous models and connects SOC with tropical geometry in log-log scale (which can be understood as the zero-temperature limit of classical algebraic and geometrical operations; see [28] and references there). It was shown that the critical state and its evolution, where the toppling function satisfies the least action principle, can be described in terms of tropical geometry [28]. Other variants of the SPA models with disorder and finite driving rates that allow the study of "avalanche queues" lead to non-universal critical behaviours [29]. The models with the correlated disorder were studied by renormalization-group in [30][31][32].
The diversity of avalanches occurring in a SOC state and the role of open boundaries in maintaining the critical state is another relevant issue that receives renewed attention, e.g., in the context of the large avalanche prediction [18]. Again, SPA models appear as useful tools in this context too. Remarkably, in [33] the class of avalanches that spill the grains over the open boundaries are identified as responsible for the observed scaling behaviour of the SPA. The diversity of avalanches can be even higher in natural systems, and more realistic models of physical systems having multiscale structural patterns and entropy complexity are required [34]. They depend on the fractal energy landscape; for example, in the bistable spin systems under the slowly increasing external fields, the spanning avalanches in different dimensions obey different scaling properties [35]. Besides the physical geometry of the sample [36,37], the implicated portion of the phase space plays the role. Thus, the avalanches originating in the centre of the hysteresis loop, when the system is away from the two global minima, have different scaling properties compared to the ones observed at both ends of the loop [36]. In social systems, a different type of scale invariance is observed when the network that supports the dynamics grows during the avalanching dynamics [38], such that each avalanche touches the boundaries but without spill-over. Notably, in both cases, the sequence of avalanches represents a multifractal set with a non-random distribution of avalanche returns.
This feature article has two goals. Firstly, we give a brief survey of a diversity of current research trends of SOC systems across different scales and types of interactions. Secondly, we present new results on the field-driven spin dynamics in complex nanonetworks, an appearing prominent example of SOC behaviour induced by the substrate's geometry. Using several representative examples of SOC systems of different nature and interaction patterns, we highlight some fundamental aspects of the dynamic complexity. The noisy output signal, representing a fingerprint of the system's functioning, is often accessible from data, e.g., in biological, social and physical systems at different scales and often measured in the laboratory and produced by numerical simulations. We briefly define several quantities as signatures of SOC that can be determined from the output signal. Furthermore, focusing on the role of the complex-geometry substrate, we simulate the fielddriven spin dynamics on nanonetworks with simplicial-complexes architecture. With the antiferromagnetic interactions among the neighbour spins, the triangle faces of the built-in simplexes provide the conditions for the spin frustration effects. By performing a detailed analysis of the simulated magnetisation noise signals, we demonstrate SOC behaviour occurring on the hysteresis loop in these complex assemblies at the interplay of geometric frustrations and external driving. Our study allows to identify and constructively describe the promising ways to investigate self-organised critical states under the competing higherorder interactions embedded in complex geometries.

Properties of the Self-Organised Critical States
In this section, we describe those fundamental properties of SOC states that are most important in our further discussion, and in doing so, we rely on some established results known from the literature. As stated above, the fingerprints of SOC are connected with the avalanching dynamics, exhibiting spatial and temporal correlations. Consequently, the corresponding quantities that characterize the occurrence of SOC states and the related universality class can be theoretically defined and numerically determined [5]. Specifically, some of the quantities that can be inferred from empirical time series data or model simulations are introduced in the following (see more in Section 4). For the continuous models, the renormalization group for non-equilibrium dynamical systems can be used; see, for example, references [30][31][32]. Recently, a new machine learning method was introduced to assess multi-scale structural complexity of emergent patterns [34].
The temporal activity correlations, 1/ f φ -noise, and multifractal fluctuations. The system's activity {m t } time series of T max data points, cf. an example in Figure 2, exhibits self-similarity with non-trivial temporal fluctuations. They are characterised by the temporal correlation function C(δt) = m t m t+δt ∼ (δt) −γ , algebraically decaying with the time lag δt. At the same time, the correlations G(δt) of the signal's profile M t = ∑ t t m t increases, i.e., G(δt) = (M t+δt − M t ) 2 1/2 t ∼ (δt) α , where α = 1 − γ/2. These features are further manifested in the power-spectrum density exhibiting a power-law decay valid for a range of frequencies f , where theoretically the scaling law φ = 3 − γ = 2α + 1 holds. These identified attributes of the system are closely related with the power-law singularity in time series, which manifests in the multi-fractal fluctuations [39] of the time series profile Y(i) = ∑ i k=1 (m k − m t t ). It is suitably studied by the detrended multifractal analysis, see [40,41] and references therein. Specifically, dividing the signal into K = int(T max /n) non-overlapping segments of the length n, the standard deviation around the local trend y µ (i) at each segment µ = 1, 2, . . . , K is determined as The process is repeated starting from the end of the signal. Then, the r-th order fluctuation function is determined as [40,41] Here, H r defines a generalised Hurst exponent, which can be determined from the slope of a straight line in the log-log plot of the fluctuation function for each r value. The parameter r takes a range of the negative and positive values, corresponding to the amplification of the small and large fluctuations, respectively. Note that for a mono-fractal time series H r H 2 for all values of r, where H 2 is the standard Hurst exponent. Self-similarity, power-law avalanche distributions, and finite-size scaling. In the SOC systems, an avalanche of the activity may occur in response to the external perturbation (driving). The driving force can act at a single site, e.g., as in the SPA, which helps define mathematical limit known as the infinitely slow driving, or at several different sites, corresponding to a finite driving rate. Moreover, the SPA paradigm allows defining an adiabatic evolution, where the driving force is halted during the evolution of the avalanche, which is theoretically appealing, in contrast to the laboratory systems (see a more detailed discussion in [42]). The avalanche size s is then determined as the number of different sites where the activity occurred before the system is quiescent again. Similarly, the corresponding time since the avalanche is triggered until the activity ceases defines the duration T of the avalanche. Other related quantities, e.g., dissipated energy, can be defined depending on the physical system. Note that these quantities can be also determined from the activity time series, considering the avalanche as a segment of the time series between two consecutive drops to zero (or to a baseline, cf. Figure 2, discriminating random fluctuations [38].) The most striking feature of SOC states is self-similarity or the absence of the characteristic length of avalanches within a given range determined by the finite system size. Ideally, their size X ≡ s and duration X ≡ T are described by a power-law function with the corresponding scaling exponent τ X Moreover, the exponential type of cut-off is induced by the finite system size. Given that different avalanche properties stem from the same nonlinear dynamics, their characteristics, e.g., the size and duration, are statistically related. Specifically, selecting all avalanches of a given duration, their average size scales with the duration T is given as s T ∼ T γ ST , and the exponent is expressed as γ ST = (τ T − 1)/(τ s − 1). We note that other properties in the dynamics, e.g., defects, can introduce a different scale, as discussed below. The scale-invariance of the avalanches has two immediate consequences. First, in the SOC state, a system-size event (catastrophic avalanche) has a finite probability. Moreover, the system size plays a vital role in determining the features of the events inside the system. Consequently, a thorough finite-size scaling (FSS) analysis is necessary, analogously to one known at the second-order phase transitions [43].
An additional relevant scale can appear along with the finite system size, leading to a two-parameter scaling function. As typical representative examples, we mention the system's specific shape [36] and the intrinsic disorder, causing new types of critical behaviour [30,31]. The scale invariance, occurring when the linear system length L is rescaled by the factor , can be determined in the avalanche distribution P(X, p, L), e.g., for the avalanche size X ≡ s or the duration X ≡ T in the presence of disorder p and the finite system of the linear size L. It is expressed by the generalised scaling function form P(X, p, L) = λ P( λ X X, λ p p, −1 L), where λ is the generalised scaling factor of the distribution, meanwhile the corresponding variables scale with their respective factors λ X and λ p . Choosing the scaling factor, for example, such that λ X X ∼ 1, one can obtain the scaling form of the distribution of X. Specifically, we have P(X, , L) = X −λ/λ X P(1, X −λ p /λ X p, X 1/λ X L), where we can identify λ/λ X ≡ τ X is the scaling exponent measured in the power-law region of the distribution, and the two scaling variables can be identified, leading to the scaling form where D X = −λ X is the fractal dimension of the variable X. Meanwhile, ξ X (p) ∼ p −D X /λ p is the characteristic length induced by the disorder variable p. The appropriate numerical values of the exponents τ X , D X , and D X /λ p enable a successful scaling collapse. We note that D X describes how the cut-off length scales with the system size, X 0 ∼ L D X , while λ p depends on the type of defects. For example, if the random distribution of defect holes is considered, then we have [44] that λ p ∼ −1 defining the characteristic linear distance between the defects. In the absence of disorder ξ X (p) → ∞, and the standard SPA scaling function of the variable X/L D X is recovered. The scaling collapse according to Equation (4) can be attempted in two limits. First, a weak disorder ξ X (p) L does not affect the scaling collapse for different system sizes as long as this length separation applies. In the opposite limit, when the disorder is strong, the length ξ X (p) L, and the scaling collapse can be attempted considering different p values and fixed large system size L, see, for example, Ref. [44] for the avalanche distributions in the case of the directed SPA with uncorrelated site defects. In between these two limits, however, a two-variable scaling function is needed. Such scaling collapse is challenging to achieve numerically, and it is rarely performed. Moreover, a disorderinduced dynamical phase transition may occur from SOC to a strong-disorder dynamical regime. The FSS with a two-variable function can then be used to determine the critical line, for example in [36,45].
The spatial propagation of avalanches strongly depends on the dynamical rules and the underlying substrate geometry. Often, the finite-size scaling form of the type (4) with a well-defined fractal dimension does not hold. Instead, a mutifractal scaling of avalanche distributions is obtained [29,46,47] according to compatible with the spectrum Ψ X (α X ) of the fractal dimensions α X . Studies of SPA models have revealed that this type of scaling is related to multiple toppling of a site during the same avalanche, so-called toppling waves in BTW model [48]. Hence, the toppling numbers at neighbouring sites may differ considerably, exhibiting a power-law singularity, in analogy to the above mentioned time series. In the directed SPA, such conditions arise in temporally disordered flows [49], as well as in the case of dissipative sandpiles with finite driving rates [29]. It should be noted that the linear scale L in the above Expressions (4) and (5) is a good measure of the system size in the case of compact lattices and regular graphs. In the case of structured graphs, the graph's diameter scales with the number of nodes N as d ∝ logN; meanwhile, the topological inhomogeneities at different scales may impact the avalanche propagation, depending on the stochastic process in question. The power-law ranking (Zipf's and Heaps' laws) and inter-activity distributions. These properties of SOC states are closely related to the avalanche evolution. For example, in the earthquake data [50], knowledge-creation [51] and biological data [52], to mention only a few, they can be determined from the activity time series. In general, it has been shown [53] that in non-ergodic systems, the reduction of phase space due to self-adaptation of the system along an evolutionary path leads to Zipf's law.

Self-Organised Critical Systems and Their Networks at Different Scales
Studies of a vast number of complex systems at different scales, natures and types of interactions have shown signatures of SOC; see the survey in the book [2]. The examples shown in Figure 3 are representative in complex systems and are along the significant research lines indicated below.  [54]; (b) network of the self-assembled nano-materials grown by the model rules [55]; (c) dialogue-based structure of the online social network, data from [56]; (d) bipartitegraph representation of the interacting agents and viruses in the model [57], indicating the infection transmission events. Different colours of nodes in (a-c) indicate the characteristic community structural features, which have a profound impact on the nature of collective behaviours. The brain as an SOC system has been studied intensively as a prominent example of biological complexity. Both the brain imaging data and laboratory neural systems were considered; see reviews [58,59]. The main idea is that the SOC attractor maintains a stable response under ever-changing inputs, thus enabling robust brain function. It was demonstrated [60] that the self-organising dynamics of the brain waves reside on balanced excitatory-inhibitory interactions on the connectome. Another research line connects the hierarchical architecture of brain networks and SOC dynamics [61], also pointing out to the structural adaptation [62]. Recently, the occurrence of higher-order structures (simplicial complexes) in human connectomes were demonstrated [54,63]. The role of such complex topologies in brain dynamics was pointed out [64]. Even though the brain data increasingly confirm the evidence of SOC, understanding the precise role of brain criticality remains a challenging question [65,66]. The underlying dilemma, as pointed out in a recent perspective paper [65], is whether the primary cause for the brain criticality lies in the demand of its robustness or improved functionality, both of which are compatible with SOC states. Given the multiple brain functions that maintain human activity, psychology, and behaviour, there is almost no easy answer to these questions. In contrast to complex artificial systems designed to improve either robustness or functionality, for example, by using SOC as a guiding rule, the biological systems have developed through the evolutionary optimization processes. During these evolutionary developments, adapting both robustness and functional properties played their role. A deeper insight into the role of criticality in balancing the robustness and functionality can be achieved by studying empirical data regarding brain disorders [67]. For example, neurology studies confirm that several disorders, such as Parkison's disease, are related to the appearance of specific protein aggregates and their subsequent diffusion among connected brain regions. In [68], where such proteinopathy inclusions in engineered human neural networks have been studied, it has been recognized that "self-organized criticality represents the critical point between resilience against perturbation and adaptational flexibility." Physical systems exhibiting well-pronounced signatures of SOC have been investigated in the laboratory experiments with liquid foams [69], criticality at the hysteresis loop in ferromagnets and ferroelectric switching [36,70,71], crackling noise in strained and porous materials [72], and shape-memory alloys [73]. On a larger scale, SOC was demonstrated in the data of natural systems from the forest fire [74], earthquakes [75] and ocean geophysical turbulence [76] to solar flares [77] and other astrophysical phenomena [2]. In these systems, the actual degrees of freedom, the nature of interactions, and driving forces are well known, enabling the description of the system's evolution within the theoretical physics concepts. The Hamiltonian operator is readily constructed to describe the system's energy at the microscopic scale, meanwhile taking into account all interaction properties and pertinent constraints, e.g., due to the system's boundaries and geometry of the underlying grid. Furthermore, considering the global energy minimisation principle in physical systems, the frustration dynamics due to competing interactions at different scales can be formally described [78][79][80]. Continuum models of SOC are based on the stochastic PDEs describing, for example, the fluctuations of height (mass) [81], or energy correlations [82] of a SPA. Analysis of these evolution equations provides new insights into the nature and stability of the critical state and the scaling region of the implicated parameters. Moreover, continuum models enable the use of the dynamic renormalisation group to determine different classes of universal scaling behaviours and their stability depending on the identified relevant parameters [30][31][32].
Social systems with cooperative dynamics, studied using the dynamical data from online social sites [38,51,[83][84][85], appear to possess the mechanisms of SOC that enable the genesis of collective values from the individual features (knowledge, emotions, and behaviours) of the communicating participants. The same active principle seems to underscore the emergence of languages [86,87], as well as the social protests [88], wars [89], and historical processes [90]. In contrast to the above mentioned physical systems, the interactions underlying social self-organisation have different attributes that are elusive to similar formal theories. However, these interactions are elucidated by the concepts of social psychology, and cognitive sciences [91]. For example, the participants' cooperative behaviour during the collective knowledge creation manifests the meaningful interactions, more precisely, the interactions "designed to serve the needs of others" [92]. Such social endeavours, where individuals engage in sharing their abilities with the community [93] resulting in a collective value, often involve the mechanisms of SOC. A striking example is the dynamics of knowledge creation via "questions-and-answers" [38]. The growing knowledge shared by the active participants is an immediate common value quantified by different collective behaviour measures [51]. Furthermore, the network of knowledge contents used by the participants [84,94] remains as a permanent "explicit knowledge" outcome of the process.
Bio-social dynamics, such as epidemics spreading, is an example of collective phenomena that emerge from the features and behaviours of individuals. The genesis of the collective social behaviours connects the social interactions with factors related to the biol-ogy of the antigen (virus) [57,95]. It remains to be studied if the SOC signatures are present in the empirical data or if they can be incorporated to affect the evolution of the epidemic.
Designed systems utilise properties of the SOC critical state to improve their functional stability and efficiency. Examples include designed computer networks [96] and job scheduling [97], but also modelling the hospital management processes [98], empirical data of urban systems [99], and stock market dynamics [100]. New classes of functional materials, for example, frustrated spin systems [101][102][103], are being designed to improve their functional properties, for example, relevant to spintronics. We note that, in this inverse problem approach, the system's elements, interactions and parameters of the dynamics are varied such that they eventually result in the desired SOC behaviour. In the following, our main exemplification of the analysis is given for a designed nano-structured system with simplicial complexes that support SOC dynamics of spin reversal. Here, we describe the key elements of the model for the self-assembly of pre-formatted groups on nanoparticles [55], which we use to grow the nano-network that underlies the spin interactions.
In the model, [55], the groups of particles that simultaneously join the growing structure are represented by cliques (complete graphs) of different sizes. The attachment of each clique is constrained by the geometric compatibility of its faces with the current structure (see also [104] for the case of defect edges), and the chemical affinity factors, controlled by the parameter ν. More precisely, a clique of the size b ∈ [2, 10] is added by sharing one of its faces, of the size q < b, with a clique present in the growing structure. Here, the faces q = 0, 1, 2, 3 · · · q max − 1 define a single node, link, triangle, tetrahedron, and so on up to the largest sub-clique q max = b − 1. The remaining n a = q max − q new nodes are added simultaneously and mutually connected such that they form the clique with the nodes belonging to the shared face. The probability for attaching along the clique's q-face is given by [55] Here, the factor denoted by c q (θ) is the number of geometrically compatible locations for docking a simplex of the order q, found at the growth step θ on the entire network. In the most general case, the geometrical factor is weighted by the chemical affinity of the system towards the addition of n a = q max − q new nodes. Notably, a more significant number of nodes can be added when ν < 0, where the cliques share a single node, in contrast to the case with ν > 0, where the probability to share progressively more significant face increases with the increasing ν. By varying the parameter ν, the sizes of cliques that appear in the process, and the number of steps θ, the model permits the growth of a variety of simplicial complexes (see the online demo on this link [105]). Some prominent examples of the emergent simplicial complexes and their structural analysis can be found in [55,104].

Hysteresis-Loop Criticality in Nanonetworks with Simplicial Complexes
The nano-network structure with simplicial complexes that enables complex interaction patterns is grown by the self-assembly algorithm described in [55,104]. It attaches cliques of different sizes b ∈ [2,10], where their population decays algebraically with size as b −2 . Compared to the general model, we consider the case without the chemical affinity (ν = 0); hence, the assembly rules are governed strictly by the geometric compatibility of the attaching clique with the growing structure. In addition, each clique has a potential defect bond that can affect the docking, as explained in [104]. The resulting assembly has an intrinsic simplicial complexes structure, which has been described in [104]. Consequently, its topological graph, i.e., 1-skeleton of the simplicial complex, is a Gromov hyperbolic graph with a small hyperbolicity parameter δ G = 1, as shown in [55].
For the present study, we show in Figure 4 several structural features of the network, consisting of N = 10,000 nodes and E = 29,483 edges, that are relevant to the spin dynamics. Particularly, the network's heterogeneity (see Figure 3b) implies a fat-tail distribution of the node's connectivity, k, the main part of which can be approximated by Tsallis q T -exponential distribution P(k) = A(1 + k/k 0 ) 1/1−q T with q T ≈ 1.43. Its average degree is < k > = 5.897, clustering coefficient < Cc > = 0.712, and the path length is < > = 6.815; the graph's diameter is d = 22. Moreover, it possesses assortative correlations among the node's degrees, according to < k > nn ∼ k µ , where µ = 1.31 ± 0.05 for k 10. In the second inset of Figure 4, we show the distribution of the number of simplexes and faces f q of the size q, where q = 0, 1, 2, . . . , 9 indicates the nodes, edges, triangles, tetrahedrons, and so on, up to the largest clique of ten nodes. As described below, this structure can support complex interaction patterns of the dynamical variables associated with different geometrical faces. To study the field-driven spin reversal dynamics on this network substrate, an Ising spin S i = ±1 is attached to each network node and the interactions among spins are enabled by the simplex faces of different order. In the most general case, the Hamiltonian is given by where h is the external field and the interaction J i 1 ,i 2 ,...,i k differs from zero when the indexes belong to one of the simplexes faces in the set F k of the corresponding size k ; the upper limit k ≤ q max + 1 indicates the highest order of the considered interaction, which is determined by the order q max of the largest clique in the simplicial complex. Note that a precise identity of nodes that make each simplex and face of the order k is necessary for determining the contribution of the interactions of that order in (7), see, for example, studies in [106,107] for the competing pairwise and triangle-based interactions. By growing the structure, we keep track of nodes that make each attached simplex. In general, such information is obtained by Q-analysis of a given simplicial complex [104,108]. It should also be noted that the leading (pairwise) interactions chiefly promote the long-range order; meanwhile, the higher-order terms determine the type of criticality and shape the hysteresis loop, in agreement with the renormalization-group theory [30][31][32] and also recent numerical studies of the dynamics on simplicial complexes [106,107,109].
In the present study, we introduce the antiferromagnetic interactions of the leading order, i.e., J i 1 ,i 2 = −1; in conjunction with the triangle faces of simplexes, they induce the geometric frustration effects [101,102,110]. Specifically, attempting to minimize the energy, the orientations of all three spins on a triangle cannot be simultaneously satisfied. Consequently, a local interaction can have a large-scale impact. The system is driven by the external magnetic field h and applying the zero-temperature dynamics [36]. Specifically, starting from a state with all spins down, and a large negative field h = −h max , the field is increased adiabatically in small steps. We note the absence of any magnetic disorder in the present system, and the natural driving rate is determined by the changing connectivity of nodes, cf. Figure 4. Hence, it is compatible with the integer steps ∆h = 1, and h max = k max + 1 corresponds to the largest connectivity node in the network.
A parallel update of all nodes is applied. For each time step t, the local field consisting of the current value of the external field and the contribution of the nearest-neighbour spins of each node is computed. Under the applied field, each spin attempts to align with its local field to minimize global energy. Alignment occurs with a probability c 1, allowing frustrated spins to select one of the possible states. A reversed spin changes the local field at its neighbour nodes, which can trigger these spins to flip in the next time step, and so on, producing an avalanche. The avalanche stops when no more unstable spins are found. Then, the field is increased again. The process is repeated until h = +h max , completing an ascending branch of the hysteresis loop, then the field is slowly decreased along the descending branch to complete the loop. A detailed program flow can be found in [110]. For the present analysis, the temporal variation of the field h, the magnetisation  The hysteresis loop shows plateaus of fractional magnetisation, cf. Figure 6b, resembling ones observed experimentally in frustrated antiferromagnetic materials [111]. Between these plateaus, the magnetisation increases in jumps of different sizes, see Figure 5 bottom panel. These magnetisation jumps comprise the Barkhausen noise signal m t , an example is shown in Figure 6a. A detailed analysis of the noise signal reveals several quantitative measures of SOC, as shown in Figures 5 and 6. In particular, the power spectrum of the signals for both ascending and descending branches, shown Figure 5 top panel, obey the power-law decay in Equation (1) for a range of large frequencies with the exponent φ 1.23 ± 0.04. Furthermore, by computing the fluctuation function defined in Equation (2) of the magnetisation jumps m t , we find multi-fractal features with large variations of the generalised Hurst exponent for a range of time-series segments n 50; meanwhile a monofractal fluctuations with the Hurst exponent H 2 ∼ 0.69 ± 0.02 occur at smaller segments, as shown in Figure 6c,d. Furthermore, we determine the avalanche size as s = ∑ drops to zero at t b and t e points. The cumulative distributions of the avalanche size and duration, averaged over ten network samples, are shown in Figure 7. They are compatible with the expression in Equation (3), where we can extract the exponents τ s − 1 = 0.44 ± 0.06 and τ T − 1 = 0.89 ± 0.12 in the segment before the stretched exponential cut-offs; note that these values are close to the mean-field exponents. We further show that these cut-offs are due to the network size and no additional scale in the system. Consequently, the finite-size scaling form (4) with a single scaling variable applies, as shown in the inset to Figure 7 right panel. We note that we used the precise cut-off value s 0 to scale sizes in each distribution. Thus, we cannot determine the fractal dimension of these avalanches, given that the variations of the relevant linear scale are minimal for finite networks.

Discussion and Conclusions
The SOC occurs in many complex systems and networks at various scales, types of interactions, and intrinsic dynamics. They all obey some universal behaviours that can be captured by the properties of the emergent critical states. These are the long-range correlations, fractality, avalanching dynamics and scale invariance. It has been understood that these properties of the critical states can provide a deeper understanding of different aspects of complexity. In particular, recent research on various SPA models and real-world systems strives to underpin self-organised critical behaviour in the mechanisms underlying the emergence of new collective features, essential for the physical and biological complexity. They also provide a more transparent structure of information stored in the critical state and reduced computational complexity. In the context of complexity, understanding the role of various geometrical constraints in the critical dynamics and hidden geometry features that enable competing interactions at different scales are of paramount importance.
As a prominent example, we have studied the field-driven magnetisation reversal in antiferromagnetic nanonetworks with the simplicial complexes architecture. The underlying dynamics is related to the motion of domain walls separating the oppositely oriented Ising spins in this complex structure. With the state-of-the-art analysis, we have demonstrated the occurrence of SOC at the hysteresis loop in these assemblies. Remarkably, the simulated Barkhausen noise signal exhibits the multi-fractal features, the power-law decay of the power spectrum, and the scaling of avalanches. It should be stressed that the critical avalanches in our assembly appear without any magnetic disorder, in contrast to the case of random-field ferromagnets on a compact lattice where the criticality is associated with a critical disorder line [36]. The origin of SOC in our system can be attributed to the interplay between antiferromagnetic interactions and the network structure with triangle faces of simplexes providing the geometric frustrations, inducing non-local effects. Moreover, the finite driving rate (appropriate to this type of systems) may have its impact. In this context, we note that for a specific type of interaction, the driving mode can potentially change the criticality related to a dynamical phase transition to SOC [112]. The studied complex nano-assembly represents a rich laboratory to explore the geometry impact on the Ising spins dynamics, for example, by varying strength of the pairwise interactions [113] that can reduce the geometric frustration. Moreover, including the higher-order interactions of different order and sign according to Equation (7) and studying how they can modify the observed SOC behaviour remains one of the outstanding open questions. To this end, we have introduced the generalised Hamiltonian in Equation (7), indicating how these higherorder interactions can be implemented based on a concrete simplicial complex architecture. The renormalisation-group approaches complementing the numerical studies would be desired for determining the relevance of these geometrically constrained interactions of increasing orders.
Besides the popular SPA models, the motion of domain walls in complex substrates appear to be suitable physical models for the study of SOC. The domain walls are extended objects that can promote long-range impact based on local interactions. On the other hand, the substrate's topology determines multiple walls and constraints their motion, potentially leading to the phenomenon of "domain-wall jamming", for example, seen in ferroelectrics [72]. Our study has shown how the network substrates with the hierarchical architecture of simplicial complexes enable the competing interactions of different orders, potentially significant driving force of SOC [12] in these complex geometries. The presented analysis, obtained results, and motivating representative examples from physical, biological, and social systems, along with their networks, should stimulate further research in SOC dynamics, providing a deeper understanding of complexity.

Informed Consent Statement: Not applicable.
Data Availability Statement: This is a theoretical study. Figure 3a,c is based on publicly available data first collected and used in references [54,56], respectively, cf. Figure caption.

Conflicts of Interest:
The authors declare no conflict of interest.