Notes on Resonant and Synchronized States in Complex Networks

Synchronization and resonance on networks are some of the most remarkable collective dynamical phenomena. The network topology, or the nature and distribution of the connections within an ensemble of coupled oscillators, plays a crucial role in shaping the local and global evolution of the two phenomena. This article further explores this relationship within a compact mathematical framework and provides new contributions on certain pivotal issues, including a closed bound for the average synchronization time in arbitrary topologies; new evidences of the effect of the coupling strength on this time; exact closed expressions for the resonance frequencies in terms of the eigenvalues of the Laplacian matrix; a measure of the effectiveness of an \textit{influencer node}'s impact on the network; and, finally, a discussion on the existence of a resonant synchronized state. Some properties of the solution of the linear swing equation are also discussed within the same setting. Numerical experiments conducted on two distinct real networks - a social network and a power grid - illustrate the significance of these results and shed light on intriguing aspects of how these processes can be interpreted within networks of this kind.

Synchronization and resonance on networks are some of the most remarkable collective dynamical phenomena. The network topology, or the nature and distribution of the connections within an ensemble of coupled oscillators, plays a crucial role in shaping the local and global evolution of the two phenomena. This article further explores this relationship within a compact mathematical framework and provides new contributions on certain pivotal issues, including a closed bound for the average synchronization time in arbitrary topologies; new evidences of the effect of the coupling strength on this time; exact closed expressions for the resonance frequencies in terms of the eigenvalues of the Laplacian matrix; a measure of the effectiveness of an influencer node's impact on the network; and, finally, a discussion on the existence of a resonant synchronized state. Some properties of the solution of the linear swing equation are also discussed within the same setting. Numerical experiments conducted on two distinct real networks -a social network and a power grid -illustrate the significance of these results and shed light on intriguing aspects of how these processes can be interpreted within networks of this kind.
Why do pink flamingos perform a pre-breeding dance in perfect unison? How does the entire Australian coral reef manage to synchronize to release millions of gametes at the same time each year? How does it happen that a number of individuals in a social network take part with increasing involvement in collective ritual events? These questions concern the problem of how multiple agents within a network can achieve an identical synchronized dynamic behavior, possibly with an amplification of the intensity of the phenomenon over time. Prompted by the ubiquity and transversality of these processes, the paper offers some novel insights into the already immense research field on the dynamics of networks of coupled identical harmonic oscillators.

I. INTRODUCTION
Up to several thousand great flamingos, during the prebreeding period, gather in southern France to perform in synchrony a variety of stereotyped movements. Such group performances, involving elaborate courtship rituals and synchronized dances, are believed to play a role in enabling reproductive synchronization 1 . In the Australia's Great Barrier Reef, many coral species tune their reproductive cycles with incredible precision, synchronizing each year all on the same day, just after the full moon in November, to collectively release millions of gametes in one of the largest scale mass reproduction events on the planet 2 . In the epidemiological context, synchronization processes govern the dynamics of disease spread through complex networks in anthropic environments. Examples are the dynamics of childhood diseases such as measles, mainly determined by the annual seasonal cycle, or the pulse vaccinationinduced synchronization in simple epidemic models. It is a peculiar feature that ecological communities tend to adjust to periodic environmental or climatic driving variables. Although ecologists and epidemiologists have essentially opposite goals, the mathematical laws underlying the different population dynamics are very similar and frequently involve internal synchronization or resonance with external driving factors 3 . We find countless examples in different settings. In the biomechanical context, Brumley et al. 4 study a mechanism of self-propulsion in ancestral eukaryotes, unicellular organisms equipped with tiny mobile flagella that allow them to move in a liquid medium. Correlated changes in the directions of the flagellar beating can be explained by the raise of some kind of large scale coordination and synchronization, even in a weakly coupled system as is the case for the green alga Volvox. In the human social setting, the emergence of synchronization processes among individuals, on a physical or mental level, can be induced by collective dynamics and participation in repetitive rituals or dances, as occurs in Zikr, a Sufi mystical dance performed by Chechen Sufis and involving hundreds of men at once. However, it is well known that synchronization phenomena do not always have positive implications. Strong synchronization may indeed be related to pathological activities as is the case for the epileptic seizures and Parkinson's disease. Hammond et al. 5 have observed an abnormally synchronized oscillatory activity at multiple levels of the basal ganglia-cortical loop, in humans with Parkinson's disease. This excessive synchronization correlates with motor deficit, and its suppression by specific therapies could suggest an effective approach to improve the living conditions of patients with severe motility disorders.
These examples show that synchronization emerges, in general, from the collaboration, competition and interaction of multiple agents when they undergo periodic and cyclic dynamics. A first interpretive model, albeit oversimplified, reduces such systems to a set of coupled identical harmonic oscillators all synchronizing on the same frequency.
The two pioneering papers on synchronization in coupled systems 6 and synchronization in chaotic systems 7 have stimulated a great deal of interest in the study of complete synchronization of coupled nonlinear dynamical processes. The first paradigmatic model of diffusively coupled identical oscillators on complex network is given in Pecora et al. 8 , where the master stability function is introduced for the first time.
Such a function accomplishes two goals: decrease the computational load and provide a guidance in designing coupling configurations that conform to the required stability. This approach is currently the most widely applied one to determine the stability of the fully synchronized state in terms of the eigenstructure of the connectivity matrix.
The interest in the study of synchronization in coupled systems grew considerably after the discovery of the small-world and scale-free properties of many natural and artificial complex networks. Barahona and Pecora 9 analyze the implications of the small-world phenomenon on the synchronization of oscillator networks of arbitrary topology. Subsequent analysis provided evidence that wellconnected networks are easier to globally synchronize than regular networks. For instance, Chen et al. 10 found that synchronization processes typically start first from a small fraction of hub nodes and then spread to nodes with smaller degrees. This is replicated at the cluster level, where a typical synchronization process starts from partial synchronization within clusters to evolve to global complete synchronization.
Ren 11 analyzes coupled second-order linear harmonic oscillators, described as point masses, under rather mild network connectivity assumptions, showing that their positions and velocities can be synchronized in the presence or absence of a leader. A similar approach is used in Zhan et al. 12 where the vibration spectra of a classical mass spring model on complex networks are investigated; and, in Shang 13 , where synchronization of a leader-follower system of coupled harmonic oscillators connected by dampers is studied in presence of random noises and time delays with an interaction topology modeled by a weighted directed graph.
In general, the presence of dissipative couplings has the effect to equate the state of the nodes so that the system enters a positively invariant set, the absorbing domain, in finite time. It is well-known that if G is a constant matrix whose eigenvalues have only negative real parts, then e tG acts as a uniform contraction; that is, ||e tG || ≤ Ke −ηt for some positive real constants, K and η. This fact guarantees that the solution of the equation is uniformly asymptotically stable and the synchronization manifold remains invariant under the flow of the equation that governs the system.
For numerical estimates and model prediction purposes, Pietras and Daffertshofer 14 stress the importance of reduction of complex oscillatory systems. A number of surveys are devoted to this topic 15 . Arola-Fernández et al. 16 study the synchronized state in a population of networkcoupled heterogeneous oscillators, showing that the steadystate solution of the linearized dynamics may be written as a geometric series whose terms represent the different spatial scales of the network. They prove that this expansion converges for arbitrary frequency distributions and for both undirected and directed networks, provided that the adjacency matrix is primitive. Also Eroglu et al. 17 propose a survey on the main theory behind complete, generalized and phase synchronization phenomena in complex networks with an approach close to the one followed in the present paper.
Investigations of the transient and steady-state dynamics of synchronous machines have significantly impacted the description of power grid networks. Bhatta et al. 18 study the swing equation and its linearization on transmission networks, highlighting the role of symmetries in reducing the computational complexity of modal decomposition 19 . Bhatta et al. 20 extend this approach to multilayer networks.
A crucial issue is the time and speed of synchronization. Grabow et al. 21 focus on the speed of convergence as a collective time scale for synchronizing systems on different topologies, ranging from completely ordered and grid-like, to completely disordered and random topologies. They find, for instance, that, at fixed in-degree, topological randomness induces shorter synchronization times, whereas intermediate randomness in the small-world regime induces longer synchronization times.
Many other contributions have been given on specific issues, such as in Skardal et al. 22 and Skardal et al. 23 , where authors introduce a synchrony alignment function that encodes the interplay between the network structure and oscillator frequencies as an objective measure for the synchronization properties of a network of heterogeneous oscillators with natural frequencies or in Su et al. 24 , where authors revisit the synchronization problem for coupled oscillators in a dynamic proximity network. Dörfler et al. 25 prove a closed-form condition for synchronization that can be stated in terms of the Moore-Penrose pseudoinverse of the Laplacian matrix. Liu and Barabási 26 , in the wider context of control theory, address the problem of how networks organize themselves to balance control with functionality.
In particular, they review the process of local pinning synchronization, when controllers, applied to a subset of nodes, help synchronize the network. Also in Lu et al. 27 , pinning synchronization in the new setting of networks of networks is investigated. Slotine et al. 28 study the role, in both nature and system design, of co-existing power leaders, to which the networks synchronize, and knowledge leaders, to whose parameters the networks adapt.
One of the most important contributions to the study of synchronization in coupled systems is undoubtedly the Kuramoto model 29 . In this model the underlying topology is fully connected, even if it was later adapted to the scenario of nearest neighbor interactions. Despite the simplicity and effectiveness of the Kuramoto model and all its variants, several works have highlighted some of its shortcomings. For example, Shahal et al. 30 address the issue of the importance of synchronization in human networks for our civilization. They study, in an accurate experimental setting, the synchronization between violin players, arranged in a network framework with full control in terms of network connectivity, coupling strength and delay. Their results show that players can adjust their playing period and cancel connections by ignoring frustrating signals, to find a stable solution. Their observations reveal that the usual models for coupled networks, such as Kuramoto model, cannot always be applied to human networks, as additional degrees of freedom enable new strategies and produce better solutions than are possible within such a model. Corroborating these observations, it is argued that, during social interactions, participants are continuously active, each modifying their actions in response to the changing actions of their partners. This continuous mutual adaptation results in an interactional synchrony to which all members contribute. Dumas et al. 31 study the brain activity of the participants in this process, and found that states of interactional synchrony correlate with the emergence of an interbrain synchronization network between right centroparietal regions. These regions have been found to play a central role in social interaction within an interindividual brainweb.
Another limiting aspect of classical models is the presence of strictly local interactions. For this reason, Estrada et al. 32 study how long-range interactions can affect synchronization dynamics, proposing a model of coupled oscillators based on the d-path Laplacian matrices. Their analysis reveals that long-range interactions improve network synchronizability with an impact that depends on the original structure.
The dynamics of opinion formation within a social network is also a widely studied field. In this kind of network, each node can be assigned, for instance, a discrete or continuous score according to his/her agreement or disagreement on a given topic. By introducing a coupling between members in order to describe the dragging effect in their opinion formation process, each member in the network can change his/her opinion according to a synchronization mechanism or according to some kind of resonant phenomenon with some source/leader status.
Related to this issue, Jenkins 33 presents a very interesting analysis of the so-called self-oscillation phenomenon, which is distinct from the more familiar phenomenon of forced resonance, whereas Tönjes et al. 34 focus on networks with influencers, a group of hubs that couple strongly and connect different parts of the network. They show that, subjecting influencers to an optimal noise, can result in enhanced network synchronization. This effect, called coherence resonance, emerges from a synergy between network structure and stochasticity and is highly nonlinear, vanishing as the noise is too weak or too strong.
In the light of the previous literature review, this paper aims to analyze, in a sense from first principles, some aspects of the synchronization and resonance phenomena, and, in general, of the transient dynamics in networks of oscillators with a heterogeneous topological structure. A unified approach is maintained in the description of all phenomena. The main new contributions and results of the present work can be summarized as follows: 1. a synthetic framework in which to inscribe the description of synchronization and resonance phenomena on networks of coupled identical harmonic oscillators with arbitrary topology; 2. a closed expression of the average synchronization time that encodes information about the topological structure of the network; 3. new evidences that, for a weakly coupled system, complete and dense networks achieve the state of perfect synchronization more rapidly than, for instance, the path network and that the complete network shows a synchronization time which is shorter for weak coupling than for strong coupling; 4. closed expressions for the resonance frequencies of the network as a function of the eigenvalues of the Laplacian matrix; 5. a detailed account of how a single node in the network can serve as a leader (influencer) or a follower (receiver) depending on its topological location, when a resonance phenomenon arises inside the network; 6. a proof that the presence of dissipative terms pulls down the possible resonance frequencies to a single one, corresponding to the synchronization frequency.
7. a new perspective on the linearization of the swing equation and its solution in the proposed setting; 8. some new mathematical properties of the polar decomposition of the symplectic matrix governing the synchronization phenomenon.
The paper is organized as follows. Section II introduces the basic notations and the general model. Section III has a didactic purpose: it shortly illustrates the case of networked coupled identical linear oscillators without damping. In Section IV, existing results are framed into the broader context described, and new ones about the mathematical properties of the involved operators and the rate of synchronization are presented. Section V is devoted to resonance phenomena and, after some recalls on the case of the single driven oscillator, provides new results related to resonant frequencies on networks when an external driving force is placed at a specific node. Section VI deals with a special case of the general equation introduced, specifically the linear swing equation and its solution. Finally, section VII applies the previous findings to two real world networks, the Zachary club social network and the Syrian power grid, by illustrating a number of interesting remarks. All the proofs of the propositions in this article are given in the Appendixes A and B.

II. PRELIMINARIES
Our purpose is basically to describe the behavior of a system of n coupled harmonic oscillators within a network structure. To model the local interactions between the n harmonic oscillators, we use an undirected graph G = (V, E) without loops, where V is the set of vertices or nodes and E ⊂ V ×V is the set of edges or links. Throughout the paper we will assume that the network is connected i.e., that a path exists between each pair of distinct vertices. The actual configuration is completely encoded in the adjacency matrix A = {a i j }, i, j = 1, . . . , n, where a i j = 1 when (i, j) ∈ E, that is when nodes i and j show an actual local interaction, and 0 otherwise. We will assume a i j = a ji , which is a perfectly symmetric interaction within a pair of nodes. The degree of a node i is defined as k i = ∑ n j=1 a i j = (A 2 ) ii and K = diag {k i } is the diagonal matrix of the degrees of all the nodes. We will denote by λ i and ψ i , i = 1, . . . , n, respectively the eigenvalues and eigenvectors of the adjacency matrix A, where the eigenvalues are listed in decreasing order: λ 1 ≥ · · · ≥ λ n . The Laplacian matrix is then defined as L = K − A. We will denote by µ i and φ i , i = 1, . . . , n, respectively the eigenvalues and eigenvectors of the Laplacian matrix L, where again the eigenvalues are listed in decreasing order: µ 1 > µ 2 > · · · > µ n = 0. In particular, φ n = 1 √ n [1, 1, . . . , 1] T = 1 √ n u T , being u ∈ R n the vector of all 1's. Similarly, we will denote by 1 = uu T the n−square matrix of all 1's and by I the identical matrix.
To have a concrete picture of the system of coupled oscillators, let us look at the network in Fig. 1. This example will be used throughout the paper as a toy model. We have n = 4 objects of mass m i , i = 1, . . . , 4, attached to a rigid support by identical springs. These objects are also connected to each other according to a network scheme with adjacency matrix A by other identical springs. All . , x n ] T be the vector of the displacements of each mass from its equilibrium position, so that x i ∈ R represents the one-dimensional (vertical, in the figure) displacement of the ith object. Let F i be the force acting on mass i due to the vertical spring and F i j be the force acting on mass i due to its interaction with mass j. Then F i = −c 1 x i and F i j = −c 2 (x i − x j ) and the total force on node i is given by The forces acting on all the nodes in the network are then collected in the force vector Let us add now the damping forces, which are proportional to the velocitiesẋ = dx(t)/dt with coefficients c 1 for the vertical damping and c 2 for the coupling between nodes. As above, in a similar way, they can be collected into the vector: The most general equation of motion of the network of coupled and damped harmonic oscillators with masses M = diag {m i } is then given by which is equivalent to the system of 2n differential equations given by with initial conditions x(0) = x 0 and v(0) = v 0 . Let us define the vector y ∈ R 2n as y := x v (6) so that, in matrix form, We will focus on the case of coupled identical harmonic oscillators, so that M = I. Hence, the system (7) becomes ẏ = Gy y(0) = y 0 (8) where the matrix G is We stress that coefficients c 2 and c 2 determine the strength of the coupling between the oscillators and that in the following the coupling is meant to be weak when c Since we are interested in resonant phenomena, we add to the system an external periodic force applied to one node. Let us imagine that such a driving force f (t), f : [0, +∞) → R, acts on node h, h = 1, . . . , n; node h will be called the leader. The total force on node i is now where δ hi = 1 only if i = h, 0 otherwise. The forces acting on the nodes in the network are then described by the vector where e h = [0, 0, . . . , 1, . . . , 0] T , with 1 in position h, is in the standard basis of R n . System (8) is then generalized by Finally, in addition to the cases described above, we will study a further variant of the system (12) in order to show how this formalism can accommodate the linearization of an important nonlinear problem. In particular, by setting c 1 = 0, c 2 = 1, c 1 = 1, c 2 = 0, and by giving the factor b(t) in Eq. (12) the form b(t) = [0, p] T , where p is a constant vector in R n , we will gain the linear swing equation describing the dynamic response of a network of synchronous machines to small disturbances.
Shortly, we can summarize in the following Table I the cases we are interested in, where c 1 , c 2 , c 1 , c 2 are the coupling terms in the matrix G in Eq. (9) and the last column specifies the form of the driving term b(t) in Eq. (12).  The first case in Section III is for illustrative purposes only and it is used to properly introduce notations. The second case in Section IV will lead to synchronization phenomena on network. The third and fourth cases in Section V are the core of the paper and lead to the resonance phenomena on networks. Finally, the fifth case in Section VI will allow us to discuss some properties of the linear swing equation and its solution.

III. NETWORK OF COUPLED IDENTICAL HARMONIC OSCILLATORS
The analysis of this first case is worthwhile since it allows the exact solution of system (8) to be written in a closed form useful for later discussion. Therefore, let us go back to system (8) and set c 1 = 0 and c 2 = 0: Let us denote H = c 1 I + c 2 L and B the square root matrix of H: B 2 = H. Since, for any value of the positive coefficients c 1 and c 2 , the matrix H is positive definite, B can be easily obtained as B = ΦΛ 1/2 Φ T , where Φ is the matrix whose columns are the eigenvectors of H and L and Λ is the diagonal matrix of the positive eigenvalues of H.
It is well known that the general solution of the equationẍ + B 2 x = 0, equivalent to system (14), is given by where x (1) 0 and x (2) 0 are complex conjugate constant vectors satisfying the following conditions: Therefore and the general solution (15) takes the usual form Let us notice that solution (18) of the system (8) can be deduced equivalently, by observing that where now matrix G in Eq. (9) is G = 0 I −H 0 .
A straightforward computation shows that, in this case, where, as above, B 2 = H. Hence which is equivalent to Let us remark that, if we set c 1 = c 1 = c 2 = 0 and c 2 = 1 in system (8), we have a network of free harmonic oscillators connected by identical springs.
In this case, H = L, and the system is described bÿ x + Lx = 0. Since det B = det L = 0, B −1 does not exist, and if rank(B) = rank(B|v 0 ) = n − 1, system (16) admits ∞ 1 solutions. In fact, if we set x Therefore, there are ∞ 1 values for the imaginary part b of the initial coefficients. Let us notice that B −1 doesn't exist because of the null eigenvalue µ n = 0 of the Laplacian L and then of the matrix B. This eigenvalue is associated with the global translation of all the nodes with the same velocity, i.e. to the stationary state of the whole network, which, in this case, is free and not bounded to any fixed support. Let us also observe that, since ∑ j B i j = ∑ i B i j = 0, i.e. the sum of the elements along a row or a column in the matrix B is null, the initial velocity vector v 0 has to satisfy Then, for the system (16) to be consistent, the initial velocity vector must satisfy the condition ∑ i v 0i = 0, equivalent to the classical linear momentum conservation.

A. Mathematical properties
Let us set in system (8): c 1 = c 2 = 1 and c 2 = c 1 = 0. This choice for the values of the parameters corresponds to a system of harmonic oscillators attached to a fixed support, free to oscillate around their rest point with constant c 1 = 1 but coupled in a damping network with coupling coefficient c 2 = 1. Basically it represents a system of n networked point-mass agents connected and interacting only by virtual dampers and can be viewed as a general distributed algorithm for the synchronization of their positions and velocities. Matrix G in Eq. (9) now becomes Let us first observe that G is a symplectic matrix, that is As a consequence, we have that if λ is an eigenvalue of G, so are its complex conjugate λ , 1 λ and, hence, 1 λ . Moreover, if λ has multiplicity m, then 1 λ has multiplicity m. Finally, both G and G −1 = −JG T J have the same set of eigenvalues. In particular, denoting by λ G± i , i = 1, . . . , n, the eigenvalues of G and ψ G± i the corresponding eigenvectors, Ren 11 proves the following proposition Proposition 1. The eigenvalues of G are the n couples of reciprocal values given by and the corresponding normalized eigenvectors are given by where µ i and φ i are the eigenvalues and the corresponding eigenvectors of the Laplacian matrix L. In Appendix B, we propose a straightforward alternative proof for the undirected case in line with the one proposed by Ren 11 for the directed case. Here, we focus on the following remarks.
Remark. Proposition 1 has the following straightforward consequence: G is a Hurwitz matrix, i.e. all its eigenvalues are negative or have a negative real part: Re(λ G± i ) < 0, ∀i = 1, . . . , n. Indeed, we may have the following four possibilities as in Table II: Let us also notice that, since µ 1 > k max ≥ 2 for any non trivial network, there is always a real negative eigenvalue. Moreover, as expected, the trace of the matrix G is equal to where k tot is the total degree of the network.
A general solution to system (8) is again of the usual form with matrix G as in Eq. (24).
Remark. Let us immediately notice that, if we apply the general time evolution in Eq. (28) to an initial state equal to an eigenvector ψ G± i as in Eq. (27), such that µ i > 0, we get that, for t → +∞: We can then identify the eigenstates of the matrix G as the initial states such that the synchronization leads to the extinction of the oscillation for all the nodes of the network.
In Appendix A, we propose some theoretical results about the polar decomposition of the matrix G in Eq. (24), which have relevance in the synchronization process. In particular, we will show that the asymptotic behavior is entirely encoded in the unitary component U of the matrix G. Since G is a real non-singular matrix, its polar decomposition can be directly obtained as where The next proposition states that, in terms of asymptotic behaviors, the evolution operators e Gt and e Ut are equivalent.
Proposition 2. For any initial condition y(0) = y 0 , both e Gt y 0 and e Ut y 0 show the same asymptotic behaviorỹ(t), for t → +∞, wherẽ and 1 is the n-square all 1's matrix.

B. Example
Let us consider again the toy model in Fig. 1. The four eigenvalues of L are: The eight eigenvalues of operators G , P , U and the four angles defined in Eq. (A7) in Appendix A are listed in Table III  In the Figures 3, 4, 5, and 6 we illustrate some different synchronization processes on the same network, with different initial conditions, as specified in the caption of the figures. For each initial condition, plots represent the position and velocity of the four nodes and the phase portrait with the corresponding limiting cycle.

C. Synchronization time
The preceding analysis gives us the opportunity to add some insights into the widely discussed topic of synchronization time and speed. What follows contributes to this topic by highlighting the role of network topology in the synchronization process in order to compare the behavior of different systems and to discriminate which systems converge more or less rapidly to the synchronized state.
In the literature, a measure of synchronization is defined as is the synchronization solution in Eq. (31), ε is a threshold and Θ is the step function (see, for instance, Grabow et al. 21 ). This index provides the percentage of nodes whose distance from the synchronization solution is less than a given threshold. As the network becomes completely synchronized on this solution, C (t) → 1.
There is actually no univocal definition of an index that measures the speed of synchronization. For example, it could also be taken as an indicator of the time at which the maximum of the absolute difference between couples of the oscillators positions, that is max i j x i (t) − x j (t) , or alternatively, the mean of all the differences between couples of positions, drops below a given threshold.
In light of the closed expression of the asymptotic solution in Eq. (31), we will refer directly to the mean absolute difference between positions of nodes and such a limiting solution.
We aim here at defining a proper synchronization time as a function of the topological properties of the underlying network. Since, in the singular value decomposition, the more negative are the eigenvalues, the faster the corresponding exponentials go to zero, we have to select exponentials with less negative exponents if we are to identify the dominant behavior that slows synchronization and to set a bound on the synchronizability of the global process. Let us consider an initial state given by y 0 = (x 0 , v 0 ) T and let us write the expansion of the solution to problem (8) on the basis of the eigenvectors of matrix G: The first term in the last step is equal toỹ(t) in Eq. (31) and it is the synchronization solution, so that All the 2n − 2 terms on the right hand side contain an exponential whose coefficient is λ G± i , i = n − 1, . . . , 2, 1. These coefficients can be complex or real but in both cases their real part is always negative, according to Table ??. Therefore the process is governed by the term that contains the least negative eigenvalue real part. The eigenvalues of G are expressed by Eq.
as an increasing function of µ i . Since µ 1 > k max ≥ 2 for any nontrivial network, there is always a real negative eigenvalue and the maximum real value is attained for i = 1, so that the less negative real eigenvalue is λ G+ . On the other hand, when λ G± i ∈ C, the less negative real part is gained for the minimum non zero µ i , i.e. for i = n − 1, corresponding to the algebraic connectivity, and it is given by 2 . Therefore, the less negative exponential coefficient is given by If we look at Eq. (32) by components then we can say that, as t → +∞, where λ S is the eigenvalue in Eq. (33) and ψ S (i) is the component i of the corresponding eigenvector (as far as this component is not zero). Now, we want to find out the time after which, for any node i, the difference between the actual solution y i (t) and the synchronization solutionỹ i (t) becomes less than a given threshold ε > 0: e tλ S (ψ S · y 0 ) ψ S (i) < ε. By solving for t, we get: Finally, the average synchronization time can be bounded by Remark. When |λ S | = µ n−1 2 , the greater the algebraic connectivity, the smaller is the synchronization time. Moreover, the algebraic connectivity is bounded by the inequality µ n−1 ≤ n n−1 k min with k min minimum degree in the network, and since the density of the network is given by δ = k mean n−1 , we have that so that, in this case, we have For n fixed then, as δ grows then the synchronization time decreases. Then we expect that network D would synchronize faster than E, which in turn would synchronize faster than F. In figure 8, panel (b), we plot again the mean deviation of the node positions at time t and the synchronization solution, as in Eq. (32), as a function of time t: we observe that red line (network D) decays on average more rapidly than blue line (network E), which decays more rapidly than green line (network F). Again, if we compute the mean synchronization times in formula (37) with ε = 10 −3 , we get: t D mean = 11.849, t E mean = 19.868 and t F mean = 26.510, whereas the empirical values are t D emp = 11.140, t E emp = 21.760 and t F emp = 26.940. This observation could be taken as a hint that a path will synchronize faster than a more integrated or even complete network. Things are actually more elaborate as argued in the following important remark.
Remark. So far we have assumed a constant value equal to 1 for the coupling coefficient c 2 between nodes. In this remark we want to show that different networks may exhibit different rankings in the speed and time of synchronization according to the values of the coupling coefficient. Let us consider again the triplet of networks D, E and F in figure 7 and let us assume as synchronization time the lower bound in formula (37) with a threshold equal to ε = 10 −6 . When coefficients in system (7) are c 1 = c 2 = 1 and c 2 = c 1 = 0, formula (37)  These values show again that the path network gains synchronization more rapidly than the complete network. However, if we lower the coupling coefficient to c 2 = 0.1, the same times become: with empirical values The different behavior is illustrated in figure 9, which is the analogous of figure 8, panel (b), for the weak coupled networks D, E and F. This example shows that, for a weakly coupled system, the complete network achieves the state of perfect synchronization far more rapidly than the path network. This is an outcome opposite to the previous one. More surprisingly, the complete network shows a synchronization time which is shorter for weak coupling than for strong coupling. This fact makes it interesting to analyze how this time can depend on the coupling coefficients through the eigenvalues of the involved matrices.
If we set c 2 = α, with α > 0, the eigenvalues of G become As α decreases more eigenvalues become complex and if α < 2 µ 1 they are all complex. Therefore, what matters is the real part which is just − α µ i 2 . The maximum nonnull value is then for all − α µ n−1 2 . For the path graph µ n−1 = 2 [1 − cos(π/n)], for the complete graph µ n−1 = n. Then we have to compare −α (1 − cos(π/n)) and − αn 2 . Since n > 2 (1 − cos(π/n)) for every n > 2, we have that − αn 2 < −α (1 − cos(π/n)), which shows that, for α < 2 µ 1 , the high connectivity of the complete network enables it to synchronize faster than other topologies. Similar results for network A, B and C in figure 7 confirm these conclusions. We move now to the third case in Table I by introducing a driving periodic force acting on or even produced by a single node. We will analyze the effect of a standard periodic sinusoidal force with frequency ω. Our purpose is to illustrate some new results about resonant phenomena on networks. In order to introduce method and notations, we quickly revise some ideas about the single driven resonant oscillator.
where e At = cos ω 0 t 1 ω 0 sin ω 0 t −ω 0 sin ω 0 t cos ω 0 t Let us consider a driving force F(t) = F 0 · sin(ωt) so that b(t) = f (t) e with f (t) = F 0 m sin(ωt). From Eq. (46), we have and, by a straightforward computation of the integral in Eq. (45), we get Let us note that, for ω → ω 0 , the asymptotic behaviors show that the amplitude of both position and velocity grows linearly in time.
Incidentally, it might be useful to mention that in the case of a periodical delta-type force F(t) = F 0 ∑ n k=1 δ (t − kT ), where δ is now a Dirac delta function and T = t n = 2π ω , the previous solution turns out to be very different. Indeed, from Eq. (46) we have and, by computing the integral in Eq. (45), we get By the Dirichlet kernels and replacing 2n + 1 = π+ωt π , we finally obtain Let us observe that the limit for ω → ω 0 in Eq. (52) shows a different behavior than in Eq. (48): when ω approaches ω 0 , only the term sin ω 0 ω π in the denominator approaches 0, while other terms represent an oscillation of frequency ω 0 . Therefore, the resonance produced by this kind of force preserves the proper frequency ω 0 while the amplitude goes to infinity uniformly for all t.

B. Network case
Let us now move to the network case. The general solution of system (12) is given by In the absence of dissipative terms, by setting c 1 = c 2 = 1, matrix G in Eq. (9) becomes where now H = I + L. In Eq. (20), we proved that Let us assume a driving force f (t) = F 0 sin(ωt) applied on node h, h = 1, . . . , n. The general solution and the corresponding resonant states are described by the next proposition.
The reader can find the proof of proposition (3) in the Appendix B. Here we want to focus on its consequences. Eq. (56) can be expressed by components as In particular, since µ 1 > µ 2 > · · · > µ n = 0, ω 1 > ω 2 > · · · > ω n = 1 and φ n (k) = 1 √ n , ∀k. Let us analyze first what happens when ω → ωĩ for some specificĩ = 1, . . . , n, that is when the frequency of the external force approaches one of the proper frequencies of the network. The diverging term in the position components in Eq. (57) is Since, in general, x sintx 0 −x 0 sintx In a similar vein, the diverging term in the velocity components in Eq.
The first conclusion is that the amplitudes of both x k (t) and v k (t), in general, grow linearly with t as ω approaches one of the proper frequencies of the network. However, things are slightly more involved than this, as discussed in the next remark.
Remark. Let us imagine that, for i =ĩ, the product φĩ(h)φĩ(k) is equal to 0.
Let us begin by examining the case in which φĩ(k) = 0. Then, even if ω → ωĩ, the node k cannot resonate with the source node h at such a frequency ωĩ. In other terms, if we make the source node h oscillate with frequency ω = ωĩ and φĩ(k) = 0, then node k is not affected by the oscillation in h. It is as if node k is transparent to the propagation of the oscillation of frequency ωĩ from the source node h. The network as a whole cannot transmit such a frequency to node k and this node does not participate to the global resonant process.
On the other hand, if it is φĩ(h) = 0, then an oscillation of frequency ω = ωĩ in the source node h cannot be spread throughout the network to any node and no resonance phenomenon can be induced from such a node with that frequency.
The only frequency that is able to always put in resonance the entire network is ωĩ = ω n = 1. In fact, if we look at Eq. (59), being φ n (h) = 0 and φ n (k) = 0, ∀h, k, when ω → 1, any receiver node k resonates with any source node h oscillating at such a ground frequency.
More, specifically, let us observe that, for ω → 1, we have represents the hk-component of the Moore-Penrose pseudoinverse of the Laplacian operator. Sometimes it is also called vibrational communicability and denoted by G (v) hk := L + hk . It can be interpreted as the correlation function between the displacements x h (t) and x k (t). Formula (63) contains a term that is exactly the product of the vibrational communicability between source and receiver nodes and the driving force F 0 sint. This term tunes the extent of the interaction between source and receiver node when ω tends to the ground resonance frequency. More generally, the ratio in Eqs. (57) and (58) contains the information about the local communicability between nodes and it conveys the role of the network topology in the resonant process from the source node h to the receiver node k at the ground frequency.
We also said that, if the driving force is placed in a node h such that φ i (h) = 0, then if ω → ω i , we have x k (t) = 0, ∀k = 1, . . . , n. For instance, if the force is placed in the node 3, being φ 2 (3) = 0, this implies that for ω → ω 2 = 2, all the nodes stay at rest. The network is not able to transmit such a frequency from such a node. Other frequencies may have effect but there will be no resonance with frequency ω 2 = 2 placed in node number 3.
In figure 11, we illustrate the position plots of the four nodes in the toy network when an external force f (t) = sin(ωt) is applied to node number 1 starting from initial conditions x 0 = [0, 0, 0, 0] T and v 0 = [0, 0, 0, 0] T . As ω approaches ω 4 = 1 from below, all the nodes are activated and enter a global resonant state, as shown in figure 9 panel (d). In figure 12, we illustrate the position plots of the four nodes under the same conditions as before but as ω approaches ω 3 = √ 2 from below; let us notice that all nodes are activated, except node number 3 (green line in the plot), which does not respond to node 1 stimuli at this frequency, as shown in figure 9 panel (c).

C. Synchronization with resonance
We conclude with some remarks on the case in which the coupling between nodes is due to dissipative terms in the presence of a periodic external force of the same type discussed in the previous subsection. In this case, the fourth in Table I, matrix G in Eq. (9) takes the form A straightforward computation shows that now 36 Let us now suppose again a driving force acting on node h of the form f (t) = F 0 sin ωt. The general solution and the corresponding resonant states are described by the next proposition.
Notice that, with respect to Eq. (26), we dropped the index G in the eigenvalues to make the expressions lighter.
Leaving the proof of this proposition to Appendix B, let us now focus on the x component and its behavior for t large. After some further steps, it can be written as and, by components, formula in which we made the final term for i = n explicit.
The question now is what frequency or frequencies are able to induce resonance phenomena in the network as t grows to infinity, bearing in mind the fact that the presence of dissipative terms still drives synchronization between nodes. This amounts to asking what are the differences between formula (59) and formula (68).
First, let us note that the ratio before the square bracket in formula (68) is always a real value. This is due to the fact that λ − i and λ + i are reciprocal and to the equality Moreover, as previously discussed, the exponential terms e λ ± i t inside the brackets go to zero, as t grows. Therefore, in general, for t → +∞ we have coefficients in the sum depending on k which prevents nodes from synchronizing.
However, when ω → 1, the last term in the sum becomes dominant and this term is common to all the nodes. Indeed, the only frequency ω that can prompt resonance phenomena is the one such that ω 2 + (λ − i ) 2 = 0 and ω 2 + (λ + i ) 2 = 0 and the only value that can make these quantities vanishing is the pure imaginary value λ + n = λ − n = i.
We conclude that, under these particular conditions, the only resonant state is exactly the one involving complete synchronization between the nodes. Therefore, in general, a driving force does not favor the synchronization process and, vice versa, a resonant state is, in general, prevented by the presence of this kind of coupling for frequencies other than ω = 1. There is, thus, only one fully synchronized and resonant state for ω → 1 and t → +∞. This state is described by the following asymptotic behavior, still linearly growing in time: In figure 13, we illustrate the position plots of the four nodes in the toy network when an external force f (t) = sin ωt is applied to node number 1 starting from initial conditions x 0 = [0, 0, 0, 0] T and v 0 = [0, 0, 0, 0] T . As ω approaches 1 from below, all the nodes are activated and enter a global resonant state, as proven above. In figure 14, we illustrate the position plots of the four nodes under the same conditions as before but as ω approaches, for instance, √ 2 from below; let us notice that, as expected, no resonant state arises within the network, in contrast with what shown in figure 12 under the same conditions.

VI. THE LINEAR SWING EQUATION IN NETWORKS
Let us now return to system (12) and consider the case in which we set c 1 = 0, c 2 = 1, c 1 = γ ∈ R + , c 2 = 0 in the matrix G. Instead of considering a driving force f (t) acting on a single node as before, we now give the factor b(t) the form of a constant vector. Therefore, let us consider the system Typical examples are small perturbations in the scheduled generation of the machine, which result in a change in its rotor angle, or small loads added to the network. The response of a power system to impacts is oscillatory, and if the oscillations are damped, the system can reach in time a steady state. Transient stability studies determine whether or not synchronization is maintained after the machine has been subjected to a transient disturbance. Each variable x i in Eq. (71) represents the corresponding node displacement and v i =ẋ i its nodal velocity relative to a rotational reference value; γ, which we assume common to all the nodes, is the damping coefficient and the vector p has components equal (or, more in general, proportional) to the small power deviation generated (if positive) or consumed (if negative) by nodes. The general solution of Eq. (69) can be written as where y(0) = y 0 (and typically, in this context, y 0 = 0) and y is such that b = −Gỹ or, equivalently, such thatỹ = [x,ṽ] T with p = Lx andṽ = 0. We now limit the discussion to balanced power grids, that is we assume that ∑ n i=1 p i = 0. Furthermore, we suppose that the power grid is connected so that the multiplicity of the null eigenvalue of L is equal to 1. Under these conditions, the system p = Lx is consistent since rank(L|p) = rank(L). The family of solutions is characterized by the fact that, ifx 1 andx 2 are two of them, thenx 1 −x 2 = wu, where w ∈ R and u ∈ R n is the vector of all 1's. Since G[u, 0] T = [0, 0] T and e Gt [u, 0] T = [u, 0] T , then the general solution y(t) = e Gt y 0 + y − e Gtỹ is the same for any w ∈ R. This means that we can use any solution of the linear system p = Lx to build the general solution in Eq. (72).
Let us illustrate these ideas on the toy network in figure  1, postponing the study of a real power grid to the next section. Let us suppose that node 3 is a generator and nodes

VII. NUMERICAL EXPERIMENTS
A. Application to a social network Let us imagine that in a social network there is a debate around a certain topic and different members have different opinions on that topic. Take, for example, the case in which there are different opinions about the qualities of a particular sportsman or of a given member within the same network.
Each member can shape a positive or negative assessment with all possible gradations from full agreement to full disagreement and he or she can share this opinion with his or her neighbors. We can model this static situation at a fixed time, for instance, by assigning +1 to perfect agreement and −1 to total disagreement and providing an intermediate graduation of the level of disagreement or agreement continuously from −1 to +1. The value 0 would represent the totally neutral position.
However, each member is allowed to change his or her opinion over time, adjusting the assessment made at a given instant. Our goal here is to figure out what entrainment action some nodes may play with respect to others in this opinion formation process.
Following, for instance, Cartwright and Harary 37 and Altafini 38 , it is reasonable to assume that, if a member has a close friend, he is much more likely to be influenced by the latter's opinion than by others' and therefore he will move in accordance with his close neighbors' opinion dynamics. In other words, he will be induced to synchronize the direction and the intensity of the velocity at which he changes his opinion according to that of his friends.
The mapping between the previous theoretical model and the social problem under examination is based on the interpretation of the continuous variables x i as a measure of the level of agreement or disagreement of each member with the subject. The opinion entrainment effect induced by a node on its neighbors -and, thus, on the network -is then modeled by the couplings described above. What we intend to focus on here is not the numerical value of each individual variable at a given instant. Rather we are interested in the collective behavior of the set of variables x i over time. Such behavior makes clear, as will be seen, the effect of the topological structure of the network in shaping processes of synchronization of members' ideas and opinions.
A celebrated example of social network in network theory is provided by the Zachary Club network, see figure 16, in which 34 members found themselves having to decide whether to agree with the position of the instructor (+1) or with the position of the president (−1). We can start by assigning an initial 'position' vector in such a way that node number 1 in red, the instructor, has position +1 and node number 34 in blue, the president, has position −1. At the beginning, at time t = 0, all the other nodes are supposed to be totally neutral, since they still have to build an opinion, and so they are assigned a 'position' equal to 0.

FIG. 16: Zachary Club network partition
We start by assuming that the social system evolves according to Eq. (28) with a matrix G as in Eq. (24). The plot in figure 17, panel (a), illustrates the position vs time of two different groups of nodes with initial conditions equal to +1 for node number 1 (the instructor), −1 for node number 34 (the president), and 0 for all the other nodes. Specifically, we choose four nodes for each one of the two communities illustrated in figure 16 and colored in red and blue, respectively for the group of the instructor and the group of the president. The four nodes in red are nodes number 4, 11, 13, 18 and the four nodes in blue are nodes number 23, 24, 30, 33. As can be seen nodes in the two different groups tend to synchronize each other from the very beginning. As time goes on, both converge to zero due to the balanced initial conditions that see neither leader prevail over the other. What is interesting and unexpected is the fact that each group seems to assume a position which is opposite to that of their respective leader. However, it could be easily explained by saying that if node 1 is initially in position +1 then it can move only in the negative verse toward 0, dragging in the negative verse also nodes that are under its main influence. Similarly for node 34. The plot in figure 17, panel (b), illustrates the position vs time of two groups, including the two respective leaders. In this case, initial conditions see all the nodes having the same initial position equal to 0 (including nodes 1 and 34); but we imagine now that at time t = 0 node 1 starts to move with velocity +1 (in the positive verse) and node 34 starts to move with velocity −1 (in the negative verse). In the plot, nodes in red are numbers 1, 4, 11, 18 and nodes in blue numbers 34, 30, 32, 33; dashed lines represent the two leader. As can be seen, nodes in the two different groups tend to synchronize each other, to synchronize with their leader and to stay out-of-phase with members of the other group for all the time, until both converge to 0 due to the balanced initial conditions. Two synchronization processes with unbalanced initial conditions are illustrated in figure 18 for the same set of nodes, as in figure  We compare now the synchronization times given by Eq. (37) under different conditions. In detail, we are interested in understanding which of the two leaders is more effective in terms of speed of synchronization. Keeping in mind the discussion in the last remark in subsection IV C, we start with a strong coupling between nodes c 2 = 1. If only the instructor (node number 1) is endowed with an initial velocity (in our simulation equal to 4) and all other nodes (including number 34) are initially held at zero, then the time at which all nodes synchronize with the asymptotic solution in Eq. (31) below a threshold ε = 0.001 is equal to t = 91.70 with a mean time equal to t = 5.29. If we now assign the same initial velocity to the president (node 34) while keeping all others (including node number 1) at zero, then the same time becomes t = 96.20 with a mean time equal to t = 5.47. The two times are basically comparable with slightly higher effectiveness for the instructor than for the president, confirming the fact that the network is divided into two balanced communities. The same conclusion is validated using other initial conditions on both the velocity and the position of the two leaders.
Let us now move to the resonant phenomena. We assume now that the social system evolves according to Eq. (56) with a matrix G as in Eq. (54). The resonant frequencies given in proposition 3 range from 1 to 4.37 and they are listed in the Table IV. Let us imagine now that an external force with frequency ω is applied to node number 1, or, in other words, that the first leader acts as an influencer by applying the same driving force in that point of the network. Let us consider two nodes belonging to the two different clusters previously considered: node number 11 in the first cluster (red colored) and node number 30 in the second cluster (blue colored). Let us consider their behavior as ω approaches two resonant frequencies, specifically ω 30 = 1 and ω 22 = 1.719026. Results are illustrated in the next figures 19 and 20.
As discussed in the remark after proposition 3, all nodes resonate at the ground frequency ω = 1. Node number 11, which belongs to the sphere of influence of the first leader, certainly resonates at the ground frequency and so does at the second frequency under consideration; on the contrary, node number 30, which falls within the sphere of influence of the second leader, while behaving similarly to node 11 at the ground frequency, does not respond to the second frequency from node 1. This is easily explained if we look at the component of the corresponding eigenvector φ 22 . Indeed, for node number 11, it is φ 22 (11) = −0.540899038, whereas, for node number 30, it is φ 22 (30) = 0.002692861. The first component is about 200 times larger than the second one, which justifies the very low influence of node 1 on node 30 at that frequency. Let us consider now the frequency ω 9 = 2.572554. Since φ 9 (1) = 0, we expect that the oscillation of node 1 with that frequency does not affect the network. This is confirmed by figure 21, where we selected two arbitrary nodes (numbers 13 and 29) in the two different clusters, which show no amplification of their own oscillation when the frequency of node 1 approaches the indicated value.
Further confirmation comes from the frequency ω 6 = 2.827755. Since φ 6 (1) = −0.0027808652 and φ 6 (4) = −0.8231724260, we expect that the network globally responds much more if that oscillation starts from node 4 than from node 1. Figure 22 shows the different behavior of two randomly selected nodes (again numbers 13 and 29) when such a frequency is applied on the network starting from different source nodes, specifically nodes 1, 4 and 34. The amplitude of the oscillation of nodes 13 and 29 grows much more quickly when the source is placed in node 4 than when it is placed in node 34 or, worse, in node 1. The two leaders prove to be not so effective at this specific frequency.
Finally, let us consider the case of the resonant synchronization described by proposition 4, by assuming now that the social system evolves according to Eq. (68) with a matrix G as in Eq. (66). Figure 23 shows the behavior of the two leaders when the first one approaches frequency ω = 1 from below. All the other nodes display a similar behavior and enter a synchronized resonant state with the two leaders. No other frequency is able to produce a state of synchronized resonance on all nodes in the network.

B. Application to a power grid
We apply now results discussed in Section VI to a small power grid, namely to the main electricity transmission network of Syria, represented in figure 24 panel (a). 39 A power grid consists of a network whose nodes are rotating machines and edges are transmission lines. Nodes can produce power (generating plants) or consume power (loads). The network topology of the Syrian grid is depicted in figure 24 panel (b). It consists of 19 nodes with 6 generators (nodes 1-6, in red), 13 loads (nodes 7-19, in blue), and 24 edges.     figure 25.
For instance, as can be observed, the plant corresponding to node 4 typically has the highest peak value, while node 3 has the shortest peak time among the examined nodes. For the sake of brevity, we do not consider a power distribution that respects the symmetries of the network, which is such that all nodes within the same symmetry cluster have the same power. In this case, it has been shown that if the vector p respects the symmetries, the response of the complete network is completely encoded in that of the quotient network 19 . Similarly, we do not consider the effects of heterogeneity in the values of the damping coefficients.

VIII. CONCLUSIONS AND FUTURE PERSPECTIVES
This paper proposes a comprehensive analysis of the potential responses of a complex network when it is host to coupled oscillatory phenomena. It adds several contributions to the existing literature, by providing, for instance, a closed expression for the average synchronization times and the resonance frequencies of the whole network. Interesting extensions for future works may be the following: (a) a study of the effect of symmetries in the topological structure of the network and a comparison between the response of the full network and of the so-called quotient network under network automorphisms; (b) the analysis of the effect of topologies described by non-normal matrices on the transient dynamics and the relationship between the degree of non-normality and the resilience properties of the networked dynamical system, or, finally, (c) the introduction of non-local interactions through generalizations of the Laplacian and of an external stress factor acting on the network, in order to assess their impact on resonance and synchronization outcomes.

ACKNOWLEDGMENTS
The author would like to thank the two anonymous reviewers for their helpful advice and insightful suggestions. We propose here some theoretical results about the polar decomposition of the matrix G in Eq. (24), with implications in the synchronization process. In particular, we prove that the asymptotic behavior is entirely encoded in the unitary component U of the matrix G, as stated in Proposition 2 in the main text. Since G is a real non-singular matrix, its polar decomposition can be directly obtained as G = G(G T G) −1/2 (G T G) 1/2 . Therefore, we set Matrices U and P have the following immediate properties: U: (a) U is a symplectic matrix, U T JU = J; U ∈ Sp(2n); (b) U is an orthogonal matrix, UU T = U T U = I; (c) U ∈ SO(2n).
P: (a) P is a symplectic matrix, P T JP = J; (b) P is a symmetric matrix, P T = P; (c) P is positive definite.
In the following propositions, we provide the expressions for both the eigenvalues and eigenvectors of the two matrices U and P. They play an important role in proposition 2.
Proposition 5. The eigenvalues of P are the n couples of reciprocal values given by and the corresponding normalized eigenvectors are given by Remark. All the eigenvalues λ P± Let us analyze now the operator U defined in Eq. (A2). The first step is to prove that it can be given the form stated in the following proposition. Proposition 6. The unitary operator U can be expressed as where A and B are the n−square matrices Coefficients in formulae (A6) suggest introducing suitable angles which will play a significant role in the interpretation of matrices A and B.
Let us define angles θ i as By introducing angles θ i , matrices A and B in (A6) can be expressed as Expression (A9) makes it clear the SO(2n)-nature of the operator U. Now we can focus on the eigenvalues and eigenvectors of U, in the next proposition.
Proposition 7. The eigenvalues of U are the n couples where θ j = 2arctan λ P+ j and the corresponding normalized eigenvectors are given by Remark. The angles θ j = 2 arctan λ P+ j , introduced in formula (A7), are the arguments of the unitary complex eigenvalues of U. More precisely, λ P+ j > 1 determine the angles π 2 < θ j < π of n eigenvalues of U and 0 < λ P− j < 1, which are reciprocals of λ P+ j , determine the angles π < θ j <

Proof. The equation Gψ
Then x i are eigenvectors of the Laplacian L and the eigenvalues λ G i in Eq. (26) are the roots of the equation λ G i 2 + µ i λ G i + 1 = 0. It is straightforward to show that the eigenvalues λ G± i are reciprocal: 1/λ G− i = λ G+ i . Moreover Normalization follows trivially.

Proof of Proposition 2
Proof. Since U = Ψ U Λ U Ψ † U and e tU = Ψ U e tΛ U Ψ † U and, according to the remark before the proposition statement, Re(λ U+ j ) = cos θ j < 0 for ∀ j = 1, . . . , n − 1; and Re(λ U+ n ) = cos θ n = 0, we have, for t → +∞: Similarly, for the operator G. In fact, since ψ G± i is given by Eq. (27) then, for µ n = 0, Now, let Ψ G = ψ G+ n , ψ G− n , . . . . Since Ψ −1 G Ψ G = I and ψ G± n † ψ G± n = 1 and ψ G± n † ψ G∓ n = 0, we have and, through the same steps as before, for t → +∞: We compute separately the two integrals in Eq. (B2): and, after some algebraic manipulations, it can be expressed as Finally, by setting ξ i = µ 2 i − 4: Finally, by substituting the explicit expressions for a i = a i (t) and b i = b i (t) in Eqs. (B4) and (B5), we get the final expressions for x(t) and v(t):

Proof of Proposition 5
Proof. Since P = (G T G) 1/2 , we look for the eigenvalues of G T G which can be expressed as 40 : The equation Then v i is an eigenvector of the squared Laplacian L 2 and the values λ i are roots of the equation λ i 2 − (2 + µ i 2 )λ i + 1 = 0: 3 + (1 + µ i 2 ) µ i 2 + 4 and this ends the proof.

Proof of Proposition 6
Proof. As a matter of simplicity and without loss of generality, we conduct the proof in the case n = 2. According to its definition, U is given by U = G(G T G) −1/2 = GP −1 = GΨ P Λ −1 P Ψ T P where Ψ P is the matrix whose columns are eigenvectors of P and Λ P is the diagonal matrix of the eigenvalues of P: Then P −1 = Ψ P Λ −1 P Ψ T P is given by