Fundamental interactions in self-organized critical dynamics on higher-order networks

In functionally complex systems, higher-order connectivity is often revealed in the underlying geometry of networked units. Furthermore, such systems often show signatures of self-organized criticality, a specific type of non-equilibrium collective behaviour associated with an attractor of internal dynamics with long-range correlations and scale invariance, which ensures the robust functioning of complex systems, such as the brain. Here, we highlight the intertwining of features of higher-order geometry and self-organized critical dynamics as a plausible mechanism for the emergence of new properties on a larger scale, representing the central paradigm of the physical notion of complexity. Considering the time scale of the structural evolution with the known separation of the time scale in self-organized criticality, i.e., internal dynamics and external driving, we distinguish three classes of geometries that can shape the self-organized dynamics on them differently. We provide an overview of current trends in the study of collective dynamics phenomena, such as the synchronization of phase oscillators and discrete spin dynamics with higher-order couplings embedded in the faces of simplicial complexes. For a representative example of self-organized critical behaviour induced by higher-order structures, we present a more detailed analysis of the dynamics of field-driven spin reversal on the hysteresis loops in simplicial complexes composed of triangles. These numerical results suggest that two fundamental interactions representing the edge-embedded and triangle-embedded couplings must be taken into account in theoretical models to describe the influence of higher-order geometry on critical dynamics.

The role of self-organised criticality (SOC) has been increasingly recognised in various complex systems, from brain functioning to social and geophysical phenomena 1 as a mechanism enabling robust functioning and emergent properties at a larger scale that reside on collective dynamic behaviours; for a recent review of general features of complex systems, see [2][3][4][5] .These are nonlinear dynamical systems repeatedly driven by environmental forces that self-organise towards an attractor with a stationary state characterised by avalanches, long-range correlations and scale invariance.These stationary out-of-equilibrium critical states enable the system's response to driving forces at all scales, thus providing robust functioning and stability.Critical dynamics out of equilibrium is currently an active research field.Two challenging issues are the impact of a dynamically changing environment and underlying complex geometry 1 .
Besides numerical simulations, the renormalisation-group (RG) theory with scaling concepts and methodologies based on quantum and statistical field theory for equilibrium phase transitions are extended to non-equilibrium dynamics and SOC; see recent reviews 6 and 7 and references there.The described models with SOC behaviours 7 demonstrate the field-theory methodologies as powerful tools to characterise collective fluctuations in stationary states of complex systems that are driven out of equilibrium.Beyond external and internal noise terms, this methodology properly considers the dynamical environment, which can critically impact the intrinsic SOC dynamics 8 .Another promising way to understand the role of the complex systems' coupling to the environment can be built by the use of quantum formalism; a representative example is the study of adaptive complex systems 9 , where the formal theory leads to a requirement that the emergent "quantum" potential needs to be effectively balanced by the environmental coupling.It should be noted that the RG theory is based on continuous field models.For example, studies 8,10,11 of SOC are based on continuous versions of the original discrete sandpile automata (SPA) model, incorporating its essential intrinsic anisotropy of diffusion on otherwise homogeneous space.However, considering more complex underlying geometries 12 within the RG theory 13,14 , in particular, those enabling geometry-embedded higher-order interactions [15][16][17][18] , remains a challenging problem.Here, we aim to highlight the interplay of higher-order geometry and emergent SOC behaviour in complex systems.
Three types of underlying geometry that can shape the collective dynamics and SOC can be recognised as networks with fixed, co-evolving, and temporal structures.The dynamical units associated with the network's nodes interact via the network-provided connections of different orders.As explained below, this separation is conditional in real-world systems, relating the time scale of the network's evolution vs the dynamics of interacting units on it.
Fixed geometry network substrates extend the idea of physics models, e.g., spins at sites of a regular 2D-or 3Dlattice to a more complex structure described by the network or Simplical complex measures; see more details below.The dynamical units associated with the network's nodes and interacting through the local or nonlocal structure may have identical dynamics, i.e., as phase oscillators 19,20 or spins on simplicial complexes 21 ; another option often used in the agent-based modelling is that each unit has individual dynamics and parameters characteristic of each agent associated with a given node 22 .The temporal evolution of each unit is subject to interactions (of different orders), driving fields, and constraints by the surrounding geometry.It can be described by a set of update rules, as in the case of cellular automata or spin kinetics or by solving differential equations with identical, for example, in the case of phase oscillators, or individual forms and parameters, as in the case of agent-based models.
Co-evolving networks represent structural patterns that appear through the dynamics of interacting units and develop over time.They can be visualised as mathematical graphs, similar to the OSN of users mentioned above, or by bipartite graphs, for example, in blog data, where indirect interactions among users (as one partition) are effectuated over another partition as posted subjects of communications; In this case, graph mapping enables advanced graph theory methods to quantitatively study the structure of these dynamical patterns and how they change over time.However, the co-evolutionary mechanisms are more subtle, involving collective dynamics with SOC signatures; see, for example, 22 and references there.The geometry emerges from the dynamics of interacting units, meanwhile, the currently existing structure facilitates the diffusion (of information, knowledge, emotion) to spread further.Unsurprisingly, in cooperative users' activities, such as collective knowledge building 23,24 , some stable network structures emerge in the other partition.For example, they represent the network of emergent knowledge (of subjects) in the case of knowledge-creation dynamics in Mathematics data or a stable "social graph" emerges among users of Ubuntu Chats that persists over the years, satisfying the social 'weak ties' hypothesis; see 25 and references there.Evolutionary processes, where the respective time scale is system characteristic, comprise one of the fundamental features of complexity 3 .In this context, the network mappings of the Brain represent different types of functional patterns rather than a fixed structure graph; for a recent survey, see 26 .A more detailed description of the human connectome is given in I A .
Time-varying geometries represent another category of structures compared to the above-discussed cases; they exhibit partial or global reconstruction, which is virtually independent of the processes of dynamical units at their nodes but occurs at the time scale of these processes.The structural changes cause nonlocal effects that can be examined directly in real and phase space or indirectly by their impact on the dynamics.This issue deserves more attention, particularly in higher-order networks.In analogy to methods of programable self-assembly of materials 27 , local reconstruction due to, for example, built-in defects in the simplicial complexes architecture 28 can cause a collapse of hierarchical structure at a larger scale.
This perspective paper focuses on the interplay of higher-order geometry features and SOC dynamics.In this context, the time scale of the structural evolution is particularly relevant, given the critical requirement of the time scale separation between the internal dynamics and driving for the occurrence of SOC.In the following, we describe the concepts of SOC and higher-order geometry with simplicial complexes and survey current research trends.As a representative example, a more detailed study is presented considering spin-reversal dynamics on a structure of geometrically assembled triangles, where the competing pairwise and triangle-embedded interactions lead to emergent SOC behaviours.The paper ends with a brief discussion with a summary of open questions and research directions.

A. Collective dynamics & Self-organised criticality in Complex Systems
As stated above, the term SOC refers to out-of-equilibrium critical behaviour occurring in a steady state near an attractor of intrinsic nonlinear dynamics without apparent phase transition.For a more recent review, see [29][30][31] and references there.Such attractors are reached by the system's response to repeated driving by external forces when the driving rate is slow compared to the time scales of the intrinsic dynamics.Self-organised dynamics possess characteristic self-similarity manifested in avalanching behaviours, temporal correlations and scaling, adequately described by theoretical concepts 8 and numerical methods 32 .
Ever since Bak and coworkers introduced it and defined the paradigmatic SPA model 33,34 , the idea of self-organised dynamics was utilised to understand mechanisms underlying complexity 1, 35,36 .Signatures of SOC are increasingly found in many complex systems across the scales in physics 30 -from nano assemblies to rainfalls 37 , geophysics and solar activity 38 , and biology 39,40 .Furthermore, collective dynamics based on social cooperation represents a specific type of SOC behaviour evidenced by empirical data analysis, for example, in human activity devoted to collective knowledge creation 24 but also in schooling fish 41 .Further examples include bio-social epidemic processes 42,43 , sociotechnological 44,45 , and socio-economical systems 46 .In addition to characteristic avalanches superimposed to cyclical trends 47 , monitoring changes in the co-evolving network geometry serves as a "blueprint" of complexity in social dynamics 25 .
Understanding Brain functions and mechanisms of Brain disorders in neurology [48][49][50] represent the most challenging issues in the science of complexity and their applications, in particular, in the era of artificial intelligence 51 .In this context, network mapping 52 and the complexity of the Brain dynamics as a SOC system of firing neurons are central questions both in theoretical and experimental research [53][54][55] .The whole Brain computational connectomics 56 and imaging techniques monitoring different Brain functions, for example, attention 57 and cognition 58 or pain processing 59 , reveal multiple scales interplay between integration and segregation processes that involve distinct Brain regions and variable communications among them.Therefore, information about the architecture of these distributed neuronal processes can be gained by modelling Brain dynamics and identifying temporal variations of local order in the corresponding Brain networks; see a recent survey in 60 .Among phenomenological models, synchronisation processes of Kuramoto phase oscillators are often studied on Brain networks.In this picture, "the human brain is a complex system comprising subregions that dynamically exchange information between its various parts through synchronisation" persisting in a natural metastable at the edge of synchrony 61 .These studies gain some importance in quantifying brain dynamics in various psychiatric illnesses and Brain disorders implicating altered metastability 61,62 .Mathematically, it was shown 63 that the structure of the Brain network does not permit a stable full synchronisation 63 .Potential mechanisms to maintain partial synchronisation with co-evolving groups of weakly synchronised nodes are demonstrated in the human connectome, as illustrated in Fig. 1, considering the core network around Brain bubs 64 .
FIG. 1: A pattern of phase evolution at negative pairwise coupling with partially synchronised groups of nodes (brain regions) in human connectome core networks, consisting of simplexes of all orders attached to the eight Brain hubs.Data from 64 .

B. Higher-Order Connectivity: Dynamics on Simplicial Complexes
In networks mapping functional connections in many complex systems, hidden geometry of higher-order relations occurs 65 .They can be revealed by advanced mathematical techniques beyond standard graph theory 66 , such as the algebraic topology of graphs 67 .They are described by aggregates of simplexes (simplicial complexes) 68 and multigraphs.For example, in social graphs 17 , Human connectome [69][70][71][72] , etc., such structures naturally evolved with the self-organised dynamics.Meanwhile, in materials design 73 such complex structures often emerge from self-assembly processes, particularly those based on preformatted building blocks, e.g., groups of nanoparticles 74 .Such higher-order geometries provide a basis for multiple interactions that play their role in the dynamics and determine the system's collective behaviours.Therefore, studies of the impact of higher-order geometries on the dynamics are vital for understanding the mechanisms underlying dynamic critical behaviours and can be used to estimate the predictability limits 75 in the system's evolution.The interplay of the dynamics and higher-order structures can be also utilised to design new methodologies of network control 76 .
Recently, attention has been devoted to the dynamics of units associated with the nodes in simplicial complexes with a fixed structure 15 .Below, we describe some key features of such geometries and embedded higher-order interactions in the case of spin-reversal dynamics and synchronisation among phase oscillators.For this type of study, an underlying simplicial complex is grown, i.e., using a model; depending on the research aims, several generative models of simplicial complexes with different emergent structures are known in the literature 28,[77][78][79] .Alternatively, a real-world network, e.g., human connectome 64 , can be considered, and methods of Q-analysis 80 applied to determine its detailed architecture 69 .
Here, we briefly describe the model introduced in 77 , where the rules for the attachment of simplexes are motivated by the above-mentioned cooperative self-assembly of nanostructured materials 74 .Specifically, starting with a single simplex of the size s = q max + 1, at each growth step t g , a new simplex is added such that it shares one of its faces of the order q = 0, 1 • • • q max − 1, i.e., a node, an edge, triangle, etc., with an existing simplex, randomly selected in the structure; the attaching probability is given by P (q max , q; t g ) = c q (t g )e −ν(qmax−q) qmax q=0 c q (t g )e −ν(qmax−q) . ( where c q (t g ) stands for the number of geometrically compatible locations that are found in the existing structure at the growth step t g .The parameter ν represents the chemical affinity of the existing structure towards adding of q max − q new nodes.Thus, simplicial complexes with a sparse architecture are grown when ν < 0, representing a "tree of cliques" that predominantly share a single node, whereas the structure becomes increasingly more compact with increasing ν > 0, where sharing larger faces is more probable.When ν = 0, the process is controlled by strictly geometrical rules.We note that the spectral dimension of the Laplacian operator associated with the adjacency matrix of the underlying graph, i.e., (1-skeleton) of a simplicial complex is another measure of its architecture that determines the nature of collective dynamics on it 81,82 .The graphs of complexes grown at different chemical affinity is characterised by the spectral dimension that varies from the random-tree dimension at ν < 0 through d s ∼ 2 at ν = 0, and continuously increasing with ν > 0 and the dimension of simplexes; see detailed analysis in 83 .The size of the newly added simplex can be fixed in advance or drawn from a given distribution of sizes; see original reference 77 for examples and detailed analysis of the architecture of complexes for different chemical affinity and their Q-analysis.Note that by construction, the distances between the building simplexes are zero, and then these simplicial complexes are 1-hyperbolic 84 ; similarly, the size of the largest clique determines the dimension of the simplicial complex; see detailed discussion in 77 .For demonstration, an example of self-assembled triangles is shown in Fig. 2. It exhibits several branches of triangles attached via shared nodes and edges.These branches often appear to belong to topological communities, here indicated by different colours; they are interconnected through large hub nodes.A representative structure of the assembled triangles based on geometrical attachment rules is used in this work in section II.In this case, we have ν = 0, hence the probability of sharing faces of the size s = 1 (nodes) and s = 2 (edges) are governed by strictly geometric compatibility rules; its spectral dimension is d s ≃ 2.1; see 83 .We note that about 21% of edges and 33% of single-node faces are shared among two or more triangles in this type of structure.Moreover, as shown in 85 , assortative nodes' correlations are observed and broad distributions of the generalised degree (number of edges per node and the number of triangles per node, which behave statistically similar in this case).The degree distribution shows a power-law decay with an exponent γ ≥ 3 for a segment of intermediate degrees, excluding the nodes with extensive connectivity (hubs).During the assembly process we enumerate each edge {< i I , j I >}, I = 1, 2, • • • E and each triangle {< i T , j T , k T >}, T = 1, 2, • • • ∆ and make the lists of the nodes that make them; such lists comprise the corresponding adjacency tensors for the geometry-embedded interactions among spins, as discussed in the following section.In the remaining part of this section, we give an insight into how the structure of simplicial complexes yields the geometry-embedded interactions in the more general cases and how they are taken into account in numerical simulations of two fundamental processes leading to collective phenomena: phase synchronisation and spin dynamics.
Spin dynamics on SC: We consider the case where Ising spins S i = ±1 is attached to each node, and interactions among spins enabled by the simplexes and their faces of different sizes k ′ = 2, 3 • • • q max + 1, where q max is the order of the simplicial complex in question.Given the pairwise adjacency matrix, as is usually the case with real-world network data, all simplexes and their faces can be identified by Q-analysis, as stated above.Alternatively, in the self-assembly model described above, by monitoring all added simplexes, we keep track of the identity of the nodes that make them.In this way, we can produce a unique list of simplexes of all sizes k ′ ; then, we can write the Hamiltonian with all possible geometry-embedded interactions in that simplicial complex as follows 1 (2) Here, h is the external magnetic field and J i1,i2,•••i k ′ stands for the interaction tensor of the order k ′ , whose elements differ from zero when the indexes (i ) match one of the simplexes in the list L k ′ of simplexes of the size k ′ , and zero otherwise.For example, J i1,i2 is nonzero when the indexes correspond to the nonzero elements of the network's adjacency matrix; similarly, the nonzero values of J i1,i2,j3 correspond to all triangles in the structure, and so on.As usual in studies of spin systems with pairwise interactions, the sign and strength of interactions and potential disorder, as well as the presence of given higher-order interactions, are determined considering the physics of the problem and the study objectives.

Synchronisation on SC:
In most studied Kuramoto model that we present here, 1-dimensional phase oscillators θ i are associated with the nodes i = 1, 2 • • • N of the simplicial complex, and the geometry-embedded interactions are provided by simplexes of different orders q = 0, 1, • • • q max ; see review in 15 .For the synchronisation of topological signals associated with faces of simplexes, see recent work in 86 and references there.Then, the evolution of phases of all nodes is obtained by solving the differential equations interconnected via interaction tensors, as described below.For this purpose, for each interaction order q, we prepare a unique list L q i of the simplexes of the order q that contain the considered node i.Then, the generalised equation of motion can be written as where ω i is the node's internal frequency and the interaction tensor B i,j1,j2•••jq = 1 when the set of indexes (i, j 1 , j 2 • • • j q ) matches one of the entries in the list L q i of simplexes of the order q, and zero otherwise; the number of such simplexes is the node's i generalised degree k (i) q .Note that the coupling function in (3) satisfies general conditions for the diffusive-like, non-invasive and natural coupling 15 .The corresponding interaction constants K q are varied to explore the system's transition to the synchronised states; see for example, 87,88 , for the case of highdimensional simplicial complexes with the pairwise and triangle-embedded interactions.It has been recognised that higher-order interactions cause new collective dynamics phenomena; for example, triangle-embedded interactions induce the broadening of the hysteresis loop and an abrupt desynchronisation transition 19 .Studies 82,88 have demonstrated the relevance of the architecture of simplicial complexes to synchronisation, quantified by varied spectral dimensions.Moreover, spectral analysis and eigenvector localisation have recently been explored to predict cluster synchronisation in theoretical 20 and experimental studies 89,90 .

II. HYSTERESIS-LOOP SELF-ORGANISED CRITICALITY WITH TRIANGLE-EMBEDDED INTERACTIONS
As mentioned above, here we present a more detailed analysis of the spin-reversal dynamics driven by the slow ramping of the magnetic field along the hysteresis loop, demonstrating the emergence of SOC due to complex geometry.The spins are situated at nodes of a large complex of self-assembled triangles; cf.Fig. 2. Therefore, we have two types of geometry-embedded interactions in the Hamiltonian (4), specifically: where L 2 and L 3 stand for unique lists of the network's edges and triangles, respectively.For this study, we set J ij = −J 2 fixing the antiferromagnetic pairwise interactions among all implicated pairs of spins, and J ijk = J 3 ; a parameter κ ∈ [0, 1] is added to balance their respective contributions.As shown below, we can differentiate the ferromagnetic and J ijk = −J 3 for antiferromagnetic triangle-embedded interactions; see also 85 .As usual, the dimensionless units can apply; thus, we set J 2 = 1 and J 3 = 1.Moreover, the magnetisation (in Bohr magnetons µ B ) is determined as M = (N + − N − )/N ∈ [−1, 1] by the balance of the respective up and down oriented spins N + and N − with respect to their number N .The external field (in the units µ B ) , h ∈ [−h max , +h max ], where h max is related to the maximum number of neighbours of the vertices, is varied in a quasistatic manner, as explained below.
Starting with a uniform state of all spins {S i = −1} and h = −h max , the spin-reversal process is driven by slow field ramping h → h + δh along an ascending branch of the hysteresis and then reversing the field to close the loop.As it is a widely accepted approach in the study of Barkhausen noise in disordered magnetic systems a zero-temperature dynamics is applied; see, for example, 91 and references therein.In particular, the spin S i at the node i flips to align along the external field which can balance the local field due to the neighbouring spins, i.e., when h loc i S i < 0. Here, is the local field acting on the spin at node i at a current value of the external field h and the summation indicates the corresponding subsets of the edges L i 2 and triangles L i 3 that contain the node i.Thus, the spin-flip at node i causes changes in the local fields of its neighbours, which can satisfy the condition to flip, and so on, resulting in an avalanche.The avalanche stops when no more spins satisfy the above condition to flip.The boundary between the domain of flipped and unflipped spins defines the position of the domain wall at a given value of the external field.When the avalanche stops, the external field is changed again by δh (adiabatic driving).We note that spin frustration 92 with the antiferromagnetic pairwise interaction among spins on a triangle prevents all three spin pairs from simultaneously ordering with the field; thus, some spins may flip back even though the above condition is fulfilled.To avoid frequent back-and-forth flips of the same spin at a given field value, flips with a probability p ≲ 1 are adopted.Moreover, without the magnetic disorder, the field ramping parameter δh = 1 is the lowest value that may move the domain wall 1,85 .In the simulations, a time step t consists of a parallel update of all spins in the network.At each time step during the spin activity avalanches, we sample the value of the external field h t , the total number of flipped spins n t , and the number of spins ordered with the field (unnormalised magnetisation) Here, our focus is on the distributions of avalanches; for three representative values of the parameter κ, in particular, κ = 0, describing purely antiferromagnetic pairwise −J 2 interactions, κ = 1, corresponding to the case where only triangle-embedded interactions J 3 are present, and the case κ = 0.5, where these two interaction types are well balanced.(See ref 85 for the analysis of hysteresis loop shape and temporal correlations of the signals {n t }.)As this figure shows, the spin-activity avalanches appear in different sizes and durations that can be described by a cumulative distribution function having a power-law segment and an exponential cut-off where X stands for the size s (the number of spins undergoing dynamics) and the duration T (measured by the number of time steps from the field rump till the activity stops).In the critical dynamics where the scaling exists, the corresponding scaling exponents obey the relation γ sT = τ T −1 τs−1 , where γ sT scales the average size of the avalanches of a given duration T , i.e., < s > T ∼ T γ sT .In Fig. 3, we show the results for the distributions of size and duration of the spin-activity avalanches for the case of purely pairwise antiferromagnetic coupling, κ = 0, purely triangular ferromagnetic coupling,κ = 1, and the intermediate case with balanced pairwise and triangular interactions, κ = 0.5.The hysteresis loops showing the magnetisation vs external field are given in the corresponding lower panel.As Fig. 3 shows, the distributions of spin-activity avalanches are exhibiting a power-law decay with an exponential cut-off according to Eq. ( 5), however, two sets of scaling exponents apply.In particular, when κ = 0, corresponding to strictly pairwise antiferromagnetic interactions, large avalanches decay with the exponents τ s − 1 = 0.526 ± 0.018 and τ T − 1 = 0.993 ± 0.023, close to the mean-field SOC 93 .Whereas, the exponents found by fitting the corresponding curves in the case of κ = 1 are lower, τ s −1 = 0.326±0.013and τ T −1 = 0.49±0.021.We note that these exponents are numerically close to the ones in the universality class of directed percolation cellular automata 94 .See more discussion below.FIG.4: The fluctuation function Fq(n) vs time interval length n for q ∈ [−4.5, +4.5] (every second line is shown for better vision), for the case κ = 0.5 of balanced antiferromagnetic pairwise and positive triangle-based interactions, using the whole signal up-down loop branch.Inset: the singularity spectrum Ψ(α) vs α for three cases of the balance pairwise antiferromagnetic and triangle-based ferromagnetic interactions, the parameter κ =0.0, κ =0.5, corresponding to the fluctuation function on the left, and κ =1.0.
To demonstrate the multifractal properties of these magnetisation reversal processes, we exploit the underlying scale-invariance of the fluctuation function F q (n) as a function of the time segment n and determine the spectrum of generalised Hurst exponents H q .Here, we analyse the time series of the number of spin flips {n t } along the entire hysteresis loop.T max is the time series length.In particular, following the procedure for the de-trended multifractal analysis described in 95,96 , the profile Y (i) = i t=1 (n t − ⟨n t ⟩) of the considered time series is determined.The profile is then divided into 2N s = 2int(T max /n) non-overlapping segments of the length n, starting from the beginning of the series and then the process is repeated starting from the end.At each segment µ = 1, 2 • • • N s , the local trend y µ (i) is determined, and the standard deviation around it is computed as Then, the r-th order fluctuation function is for the segment of the length n is obtained and averaged over all such segments as For varying the parameter q ∈ [−4.5, 4.5], we explore the scale-invariant parts of the fluctuation function F q (n) against n, where the segment length is varied in the range n ∈ [2, int(T max /4)], cf.Fig. 4. In this way, the segments of the signal representing small and large fluctuations, respectively, are enhanced in different ways to maintain the self-similarity of the whole time series.Here, a generalised Hurst exponent H q depends on the amplification factor q, as shown in the inset of Fig. 4. Intuitively, this means that the small fluctuations (amplified when q < 0) along the signal have different scaling properties than the large fluctuations ( q> 0).In contrast, H q = H 2 for all q, where H 2 is the standard Hurst exponent for the monofractal signals.Having the spectrum of the generalised Hurst exponents H q , we determine other multifractal measures; see, for example, 95 .Specifically, we obtain τ q = qH q − 1, which is related to the standard (box probability) measure, and the singularity spectrum Ψ(α) = qα − τ q , where α = dτ /dq; in this context, different α values indicate variations in the power-law singularities along the considered time series.In the inset to Fig. 4, we show how the shape of the singularity spectrum changes when the balance between antiferromagnetic pairwise and triangle-embedded interactions is varied via the parameter κ.

III. DISCUSSION AND FUTURE DIRECTIONS
In this paper, we have given a survey of recent research activities aiming at understanding the role of higher-order connectivity in functional complex systems.Our leading idea concerns the intricate interdependences of the higherorder structures and the out-of-equilibrium self-organised critical dynamics as a plausible mechanism for emerging new properties at a larger scale, the central paradigm of physical complexity.Having briefly described both issues, we have mentioned an increasing number of studies that evidence the occurrence of higher-order (hidden) structures and signatures of SOC in many complex systems across the scales.We have highlighted the most striking example, mapping the brain structure and its dynamics, which has considerable influence on the science of complexity theory and practice, as well as the applications in artificial intelligence.Furthermore, we highlighted the relevance of different time scales in the light of the structure-dynamics interdependences.In this context, the time-scale separation between the intrinsic dynamics and driving is required for the nonlinear systems to evolve towards SOC attractors.In addition, we stress the relevance of the time scale of changes in the underlying geometry evolution as compared to the SOC dynamics; they can be characterised by the appropriate graph or simplicial complex measures.In this regard, we differentiate two limiting cases, particularly structures co-evolving with the intrinsic dynamics and the quenched structure, whose evolution time exceeds both time scales relevant to the SOC.We referred to several pertinent examples of both groups of systems studied in the literature.Another interesting, much less understood case concerns the underlying geometry reconstruction at an intermediate time scale between the intrinsic dynamics and driving.It can occur virtually independent from the dynamics of the actual units comprising the system, for example, reconstructions via built-in temporal defects, or can have a hidden connection with the driving forces.We have briefly described two currently most studied dynamic models-synchronisation of phase oscillators and zero-temperature spin dynamics on simplicial complexes of a fixed architecture, which provides geometry-embedded higher-order interactions.A more detailed study of the spin-reversal dynamics with the competing geometry-embedded interactions leading to SOC behaviour is given as an illustrative example.
In the studied example with Ising spins on an assembly of triangles, we have shown that the spin reversal zerotemperature dynamics with competing interactions of different orders when driven by avalanches-adapted field changes lead to self-organised critical dynamics on the hysteresis loop.We observe two different universality classes of SOC that are associated with the order of the primary spin interactions.Here, we emphasise that the system with triangleembedded interactions without background pairwise couplings evolves towards the SOC state in a new universality class, compared to the one induced by the pairwise interactions.Specifically, for κ = 1 where only triangle-based interactions are present, the exponents found by fitting the corresponding curves are τ s − 1 = 0.326 ± 0.013, and τ T − 1 = 0.49 ± 0.021, which are close to the exact exponents of SOC in the directed sandpile automata τ s = 4/3 and τ T = 3/2, derived in 94 .In contrast, numerically determined exponents τ s −1 = 0.526±0.018and τ T −1 = 0.993±0.023close to the mean-field SOC 93 , i.e., τ s = 3/2 and τ T = 2, are found in the limit κ = 0, where only the antiferromagnetic pairwise interactions are present.
To understand the nature of these dynamic critical states, we refer to the theoretical concepts regarding the equivalence between the zero-temperature Ising spin dynamics and the directed percolation cellular automata 99 .Furthermore, we recall that spin-reversal avalanches triggered by the external field change propagate as a directed branching process, where the actual architecture of triangles determines the underlying branching tree.In particular, for fractal networks, where the average branching number saturates for significant distances, the study 97 utilising spanning trees based on the betweenness centrality of edges reveals that the probability of cluster size s (mass when the box length ℓ B → ∞) scales with the exponent τ s = γ/(γ − 1); the exponent γ is related to the degree distribution, i.e., the number of edges per node).Then we can compute the corresponding duration exponent τ T from the scaling relation 98 τ s = 2 − 1/τ T , which gives τ T = (γ − 1)/(γ − 2).We note that the generalised degree distribution for the network of finite size, considered in the above numerical simulations, has an intermediate scale-free segment compatible with the exponent γ ∼ 3.8 ± 0.4, excluding the hub nodes; meanwhile, s smaller exponent applies for the low-degree nodes.It is relevant to mention that, in this particular case, the number of triangles per node scales similarly to the number of edges.Apart from the structure of potential branches, a crucial difference regarding the nature of avalanche branching occurs due to the spin frustration effects in the limit κ = 0.In this case, each triangle contains at least one frustrated spin, which prevents propagation of the local order set by the current value of the external field.The situation corresponds to a critical branching process, where the average branching number remains constant, which leads to the mean-field SOC according to the original work in 93 .On the other side, for triangle-based interactions only (κ = 1), no spin frustration occurs.Multiplicative branching processes are enabled by the network's fractality, leading to the exponents as above, but now the structure is related to the architecture of triangles 85 .Given the well-known equivalence of the (zero-temperature) Ising model dynamics with the Domany-Kinzel cellular automata for directed percolation 99 , as well as the related directed sandpile automata with probabilistic toppling 98 , we note that the numerical values of the exponents, in this case, are close to the ones of the directed percolation, which are known with more significant numerical precision, i.e., τ s + 1.32059, τ t = 1.47244, and the anisotropy exponent ζ = 0.63261); see 98 and references there.
In conclusion, the presented survey and the numerically studied example of the hysteresis-loop criticality in a simplicial complex with geometry-embedded interactions highlight the relevance of self-organised dynamics in functional complex systems and the critical role of geometry in the emergence of this type of collective behaviour.In this regard, our results suggest that triangle-embedded interactions play the primary role in critical dynamics beyond standard pairwise interactions in systems with higher-order geometry.These studies call for more theoretical investigations, i.e., by the renormalisation-group methodology, to confirm the existence of different fixed points representing the classes of emergent nonequilibrium critical dynamics and their stability.To this end, equivalent continuous models that consider features of complex geometry (at the local-to-mesoscopic scale) that are relevant to emergent critical dynamics need to be developed, which creates a challenging theoretical problem.We mention some ideas from the physics of hierarchical spin glasses 100,101 , mapping directed cellular automata onto zero-temperature Ising model 99 , SOC models with quenched disorder noise with specific correlations 8,11 , topological defects and the use of tropical geometry 102 that may serve as an inspiring point in this direction.Another interesting open problem regards the potentials of quantum formalism, such as in reference 9 , to elucidate the interplay of geometry and time-varying driving forces in the occurrence of self-organised critical states.On the other side, several open questions remain accessible to numerical modelling and empirical data analysis using the theoretical concept of self-organised critical behaviour 1 .Some examples include the geometry reconstruction associated with the avalanche driving rate, the role of extended defects, and the occurrence of different types of SOC in systems with "meaningful" interactions among social subjects in contrast to the physics laboratory systems and others.
One of the critical ingredients of the developed framework that we highlighted and detailed here is the interplay between self-organized critical dynamics and higher-order connectivity networks.Not only does it reveal new lines of theoretical research, but it also offers a fresh perspective for data analysis based on these scientific concepts.Namely, self-organized criticality as an active mechanism in the evolution of a complex system manifests itself in characteristic patterns of spatiotemporal correlations that a perceptive analysis of data can reveal and thus enable the use of appropriate theory for further predictions (in contrast to considering an ever-increasing amount of data as uncorrelated sets).Such a combined analysis would allow a deeper understanding of the functional characteristics of many complex systems, for example, the brain, which may allow a more critical approach in designing AI algorithms that utilize brain functions and various medical implications.
Acknowledgments.BT research supported by the Slovenian Research Agency the program P1-0044.RM is grateful to the NSERC and the CRC Program for their support, acknowledging also the support of the BERC 2022-2025 program and the Spanish Ministry of Science, Innovation and Universities through the Agencia Estatal de Investigation (AEI) BCAM Severo Ochoa excellence accreditation SEV-2017-0718.
Author contribution statement.Conceptualisation (BT, RM), computation and graphics (BT), data analysis (BT, RM), writing and editing (BT, RM) Data Availability Statement.The authors declare that the data supporting the findings of this study are available within the paper, and the related references.

FIG. 2 :
FIG. 2: The simplicial complex of triangles self-assembled under chemical affinity and geometric compatibility rules in (1) with ν = 5; a close-up view with hubs is shown, and colours indicate mesoscopic communities.

FIG. 3 :
FIG.3: Top panels: Cumulative distributions of the avalanche size (black circles) and duration (red squares) averaged over the entire hysteresis loop for κ = 0.0, 0.5, and 1.0, as indicated on the panels.For the case κ = 0.5, the distributions for J3 < 0 are shown; see text for a detailed description of fits.Lower panels: Hysteresis loops corresponding to the values of κ in the panel above; forward/backward branch shown in black/red; full and dotted lines for κ = 0.5 and κ = 1 are for J3 > 0 and J3 < 0, respectively.