Functional control of network dynamics using designed Laplacian spectra

Complex real-world phenomena across a wide range of scales, from aviation and internet traffic to signal propagation in electronic and gene regulatory circuits, can be efficiently described through dynamic network models. In many such systems, the spectrum of the underlying graph Laplacian plays a key role in controlling the matter or information flow. Spectral graph theory has traditionally prioritized unweighted networks. Here, we introduce a complementary framework, providing a mathematically rigorous weighted graph construction that exactly realizes any desired spectrum. We illustrate the broad applicability of this approach by showing how designer spectra can be used to control the dynamics of various archetypal physical systems. Specifically, we demonstrate that a strategically placed gap induces chimera states in Kuramoto-type oscillator networks, completely suppresses pattern formation in a generic Swift-Hohenberg model, and leads to persistent localization in a discrete Gross-Pitaevskii quantum network. Our approach can be generalized to design continuous band gaps through periodic extensions of finite networks.


I. INTRODUCTION
Spectral band gaps control the behavior of physical systems in areas as diverse as topological insulators [1,2], phononic crystals [3], superconductors [4], acoustic metamaterials [5], and active matter [6]. In addition to ubiquitous physical network models [7][8][9][10] ranging from aviation [11] to electronics [12], there is also considerable interest in virtual or computational networks [13] with fewer physical constraints, such as those recently used to create spiral wave chimeras in coupled chemical oscillators [14]. Often, dynamics in such systems depend on the graph Laplacian [15,16] and in particular on its spectrum of eigenvalues. Traditionally studied in periodic lattice graph models [3,5,6,17] and more recently also in hyperuniform systems [18], the targeted design of spectra of any desired shape remains a major challenge in modern materials science [5,19]. Recent breakthroughs in 3D printing [20][21][22][23] and lithography [24] make it possible now to produce and explore network structures that go beyond the traditionally considered periodic lattice geometries.
Building on such experimental and theoretical progress, we present here a mathematically rigorous solution to the longstanding question of how any desired spectrum can be realized exactly on a suitably designed positively-weighted network. Specifically, our construction of networks with specified eigenvalues allows us to place arbitrary gaps in the spectrum of the network Laplacian L = D − A, where D and A are the weighted degree and adjacency matrices respectively. These gaps, finite analogs to band gaps in continuous systems, enable precise control over the dynamics in a wide range * aforrow@mit.edu † dunkel@mit.edu of graph-based physical systems. While in a strict sense band gaps can only exist in an extended system with continuous energy bands, to follow the analogy we will name an eigenvalue-free region in our finite networks that is comparable to the range of eigenvalues a discrete band gap (DBG). That said, our construction can also be used to create continuous band gaps (Section IV). Designing a suitably weighted network topology in this way presents an alternative to control procedures based on adjusting model parameters or initial conditions on a given network [25]. The spectral approach towards functional control of network dynamics proposed here can, for example, be directly implemented with recently developed computer-coupled oscillator setups [14]. After summarizing the main mathematical result, we demonstrate its broad applicability explicitly for classical and quantum systems, by showing how suitably placed DBGs can induce chimera states [26] in oscillator networks, inhibit structural growth in generic higher-order pattern formation models, and facilitate state localization in quantum networks. In parallel, we illustrate how our construction can be combined with sparsification algorithms [27,28] to yield simplified networks preserving DBGs. This approach complements the more traditional procedure of constructing graph ensembles with predefined statistical adjacency characteristics [29][30][31][32]. Finally, we discuss periodic extensions of finite networks as a systematic procedure for designing continuous band gaps.

II. NETWORK CONSTRUCTION
The problem of recovering a network from its eigenvalues has been studied extensively, both from an algorithmic [33][34][35] and mathematical [36,37] perspective. However, with a few limited exceptions [37], most prior work has focused only on unweighted networks [38], where  The eigenvalues for the graphs in (b). The mode on the complete DBG network with the k-th largest nonzero eigenvalue is supported on the first k + 1 vertices, counted counterclockwise from the top red vertex, and highly localized on vertex k + 1, which is colored to match in (b). Grey lines indicate the borders of the unstable region for the Swift-Hohenberg model with the parameters used in Fig. 3. (d), Sparsified networks retain a significant gap even for relatively large . Each point shows the mean number of edges and gap size at fixed between 1 (left) and 0.01 (right), starting from a graph on 200 vertices designed to have 100× eigenvalue 5 and 99× eigenvalue 20. The solid curve shows the worst-case gap estimate, reduction by a factor 1 − 5 3 . Sample size is 1000 for 0.1 and 300 for < 0.1. Error bars are ±1 standard deviation; horizontal error bars are smaller than the marker size.
there are a finite number of graphs on n vertices and thus only a finite number of possible spectra. We here construct an exact solution for the weighted case.
Our main result is that, given a set {λ i } of desired eigenvalues ordered so λ 1 . . . λ n−1 λ n = 0, there is a weighted graph G on n vertices with nonnegative edge weights whose Laplacian L has spectrum λ 1 , . . . , λ n−1 , 0. The Laplacian, which determines the graph, can be reconstructed from its eigenvalues and eigenvectors with the eigenvalue decomposition; we therefore need to find a set of eigenvectors that together with {λ i } give a graph Laplacian. In fact, the same set suffices for any spectrum. These eigenvectors are strongly localized: the inverse participation ratio (4norm) v (k) 4 4 = 1−2k −1 +O(k −2 ) indicates near-perfect localization v (k) 4 4 → 1 for almost all k, itself a desirable phenomenon [16,39]. As the v (k) are mutually orthonormal and orthogonal to the vector of all ones, the matrix L = n−1 k=1 λ k v (k) v (k) has the desired spectrum with 1 √ n 1 as the final eigenvector for k = n with eigenvalue zero. By explicitly computing the sum over k for i < j, we find that is, that the off-diagonal elements of L are all nonpositive (Appendix B); L therefore corresponds to a graph with nonnegative weight −L ij between vertices i and j. If all of the eigenvalues are nonzero, all of the off-diagonal elements of L will be nonzero and the resulting graph will be complete. Some spectra can only be realized on complete graphs. A graph G with approximately constant spectrum must be complete: if L has nonzero eigenvalues λ k = λ + k for k < n and λ n = 0, then The off-diagonal elements of λI − L, which equal the original edge weights of G, are therefore λ n + O( ). For small k , every edge has nonzero weight. More generally, our construction also shows that the spectrum of any non-complete weighted graph with no isolated vertices cannot uniquely specify that graph, in line with older results on, for example, the spectra of trees [36].
This construction allows us to create networks with precisely specified gaps. For instance, choosing λ 1 = λ 2 = . . . = λ n/2−1 and λ n/2 = λ (n/2)+1 = . . . = λ n−1 leads to a graph with edge weights −L ij = λ n−1 /n if i > n/2 or j > n/2 and −L ij = (2λ 1 −λ n−1 )/n otherwise (Appendix A); that is, there are two groups of vertices, one strongly connected within itself and one weakly connected to everything. Adding a small amount of noise to each eigenvalue then lifts the eigenvalue degeneracy while preserving the connectivity structure and retaining a gap (Fig. 1b,c).
Since complete graphs can be difficult to realize physically, we explore the effect of the sparsification-byresistances algorithm developed by Spielman and Srivastava [27]. Given an accuracy parameter , this sparsification creates a network with O(n log(n)/ 2 ) edges whose eigenvalues match the eigenvalues of the original network to within a multiplicative factor 1 ± with high probability. Sparsification by resistances aims to preserve the entire spectrum, not just a gap; future sparsification algorithms directly constructed to preserve a gap could therefore improve on its efficiency. In other applications, the networks of interest are virtual ones [14] and sparsification may not be necessary.
We can use the 1 ± multiplicative error bound to estimate the size of a discrete band gap after sparsification. Suppose we start from a network with eigenvalues λ 1 , λ 2 , and 0, with some multiplicities, where λ 1 > λ 2 . The eigenvalues {µ i } of the sparsified graph corresponding to λ 1 should be no smaller than µ i (1 − )λ 1 , while the eigenvalues {ν i } corresponding to λ 2 should be no larger than ν i (1+ )λ 2 . The sparsified graph should therefore have a gap ∆ = min i µ i − max i ν i of size That is, the gap contracts by a factor at most 1− λ1+λ2 λ1−λ2 . For the parameters used in Fig. 1d, this is (1 − 5 3 ).

III. APPLICATIONS
We now demonstrate the practical potential of DBGs with three generic network models. In each case, we compare the dynamics on a complete DBG network (Fig. 1b, left) both to a sparsified approximate DBG network (Fig. 1b, top) and to a random connected network (Fig. 1b, bottom) constructed to have the same weighted vertex degrees as the DBG network (Appendix C). The gap is approximately preserved in the sparsified network and vanishes entirely in the random graph (Fig. 1c,d).
Matching the degrees in the random graph to the DBG network ensures that any differences in dynamics are due to the gap and not differences in coarse features like the average connectivity. Often, the behavior of optimized networks is sensitive to small perturbations [40]; here, behaviors preserved in the sparsified graph are robust to significant changes.
Simulations were performed using a third or fourth order Adams-Bashforth linear multistep method with a time step ∆t = 10 −4 . All simulations were written in C++ using Armadillo [41].

A. Kuramoto oscillators
Our first application is the Kuramoto model of coupled oscillators [42,43].
Recent experiments coupling Belousov-Zhabotinsky reactions via a computercontrolled projector have shown the emergence of chimeras [14]; we will show this can be achieved in the Kuramoto model by designing an appropriately gapped spectrum. Here phases θ i on the vertices evolve with a natural frequency ω and a nonlinear coupling defined by the network adjacency matrix: On any connected graph with α = 0, there is a single global attractor θ i = θ 0 + ωt. The rate of convergence to this attractor is controlled by the eigenvalues of the Laplacian L [15]. Both the complete and sparsified graphs have no eigenvalues near zero, so they synchronize much faster than the random graph ( Fig. 2a-c). The gap divides the modes into two groups, one synchronizing faster than the other (Fig. 2d-e); moreover, on the complete graph, the localization of the eigenvectors causes staggered synchronization of vertices (Fig. 2a).  If α is sufficiently large, the oscillators no longer synchronize at a single frequency. On DBG networks, global coherence gives way to weak chimera states [44] where vertices synchronize into two clusters with distinct frequencies ( Fig. 2g-i, Movie 1). For the exactly-gapped network with edges of weight w 1 = λ n−1 /n + (λ 1 − λ n−1 )/(m + 1) or w 2 = λ n−1 /n there is a steady state with θ i = θ 1 for i m + 1 and θ i = θ n for i > m + 1. In this state, The two phases θ 1 and θ n can synchronize if This synchronization is possible if α is small enough that the right hand side is less than one. If the two groups do not synchronize, and nw 2 = λ 2 is not too large, the sine term in Eq. (7) will average to nearly zero giving an approximate mean frequency difference d dt (θ n − θ 1 ) ≈ [nw 2 − (m + 1)(w 1 + w 2 )] sin(α), (8) which reduces to −(λ 1 − λ n−1 ) sin(α) if m + 1 = n 2 as in Fig. 2. More general cluster synchronization [45,46] could be achieved by adjusting the number and size of the gaps. In contrast, the random graph becomes thoroughly incoherent at comparable values of α (Fig. 2i,l). The coherence can be quantified by the order parameter r = | j e iθj |, which oscillates for the complete and sparsified networks but is near zero for the random graph ( Fig. 2jl), indicating complete disorder.

B. Swift-Hohenberg pattern formation
As the second application, we study generic Swift-Hohenberg pattern formation dynamics on a network [47,48]. Consider a scalar field u i on the vertices obeying The fixed point u i = 0, which exists for any values of the parameters D 1 , D 2 and α, is linearly stable to perturbations in a Laplacian eigenmode with eigenvalue λ if the growth rate σ ≡ −α − D 1 λ − D 2 λ 2 < 0. With α and D 2 positive, σ is negative for small and large λ, but choosing D 1 < −2 √ αD 2 creates a range of unstable λ in between. This can drive pattern formation that is eventually stabilized by the nonlinearity. The patterns can only form, however, if L has eigenvalues in the unstable range. Controlling the spectrum of L therefore allows us to completely suppress pattern formation in arbitrarily large systems by placing a gap around the unstable region (Figs. 1c, 3a, Movie 2). If we sparsify the network with sufficiently small , the gap will be preserved and again no patterns will form. Eventually, though, increased sparsification will push some eigenvectors into the edges of the unstable region and bring back partial pattern formation (Fig. 3b), which becomes fully developed in the random graph (Fig. 3c). The maximum for which patterns will be fully suppressed for given parameter settings can be predicted straightforwardly from the expected changes in the eigenvalues, in a similar fashion to the post-sparsification gap size in Fig. 1d.

C. Gross-Pitaevskii localization
Having discussed two classical applications to nonconservative systems, we now show how DBGs can control quantum dynamics with conserved energy. In a network version of the Gross-Pitaevskii model [49][50][51], (a-c), When the wavefunction in the Gross-Pitaevskii model of Eq. (10) is initialized at a weakly connected vertex with low kinetic energy, localization or delocalization (indicated by high or low potential energy, respectively) is controlled by the interplay between the graph spectrum and the rate of potential energy loss g. The random graph (purple) always delocalizes, due to its dense spectrum. However, while the sparsified graph (yellow) can delocalize for low g (a) and high g (c), again due to available eigenmodes, intermediate g (b) places the range of allowed modes inside the spectral gap, preventing delocalization. The complete graph (blue) always inhibits spreading due to the extreme localization of its eigenvectors. similar to those used in studying Bose-Einstein condensates in optical lattices [52], we find that the interplay of the total energy conservation constraint in such a model with the kinetic energy gap inhibits spreading of the wavefunction on DBG networks. Similarly to the Swift-Hohenberg example, we take the Gross-Pitaevskii equation for a complex wavefunction ψ and replace the continuous Laplacian ∇ 2 with its discrete analog −L: This can be written i dψi dt = ∂E ∂ψ * i , where the energy E is the sum of the kinetic energy T = i,j ψ * i L ij ψ j and the potential energy V = 1 2 g j (ψ * j ψ j ) 2 . The potential energy quantifies the localization of ψ: with g > 0, it is large when the probability ψ * ψ is concentrated at a single vertex and small when ψ * ψ is spread out. Delocalization is limited by the size of the network, as V g 2n , but can vary widely even on a finite network. If ψ is initialized at a single vertex j, then V = g/2, independent of j, while T = L jj equals the degree of j.
Since energy is conserved, the wavefunction can delocalize and reduce its potential energy only by converting it to kinetic energy. The rate of potential energy loss, set by g, must therefore match the rate of kinetic energy gain, set by the differences in eigenvalues among the modes involved. Suppose the wavefunction is mostly in a localized mode j with eigenvalue λ j . Spreading to a higher mode k with λ k − λ j g would increase kinetic energy by more than it would decrease potential energy, while a weak higher mode 0 < λ k − λ j g or a lower mode λ k < λ j would not increase kinetic energy by enough, if at all. Both are barred by energy conservation. The amplitude in mode j can only be reduced if there are other modes k with λ k ∼ λ j + g.
To see this in more detail, suppose we have a wavefunction comprising two modes, j , with initial complex amplitudes c 1 , c 2 . Suppose also that these eigenmodes are localized on two different vertices, with v (1) ≈ (−1, , , . . . , ) and v (2) ≈ ( , −1, , . . . , ). The system energy as a function of c 1 and c 2 is then If the squared amplitudes change slightly, to |c 1 | 2 − δ and |c 2 | 2 + δ, the change in energy to leading order in δ is Conservation of energy requires ∆E = 0, so in order to transfer a noticeable amplitude δ from the first mode to the second we must have λ 2 − λ 1 + g(|c 2 | 2 − |c 1 | 2 ) ≈ 0. In the cases considered in Fig. 4, where |c 1 | ≈ 1 and |c 2 | ≈ 0, this reduces to λ 2 − λ 1 ≈ g. Thus on a network with a spectral gap, the localization of ψ can depend nontrivially on the interplay between g and the spectrum.
Initializing ψ at a weakly-connected vertex brings out this interplay as g is varied (Movie 3). The initial state, with high potential energy and low kinetic energy, is localized on modes with eigenvalue below the spectral gap. On the sparsified network, a low value of g makes nearby modes below the gap accessible for delocalization, causing the wavefunction to spread (Fig. 4a). However, increasing g pushes the region where transfer is possible inside the spectral gap, inhibiting the spread of the wavefunction on the sparsified network (Fig. 4b). Further increase of g once again enables delocalization as the modes above the gap becomes accessible for energy transfer (Fig. 4c). In contrast, the dense spectrum of the random graph means delocalization occurs in all three instances (Fig. 4). Interestingly, the complete DBG network appears to remain localized for all values of g in our simulations (Fig. 4); this is likely due to the strong localization and near-zero overlap of the eigenmodes.

IV. BAND STRUCTURE IN PERIODIC NETWORKS
We can construct infinite periodic networks in a standard way from any base network G by tiling periodically and rewiring edges (Fig. 5a,d). Starting from the original vertex set {j} for 1 j n and edge weights −L jk , we make an infinite string G ∞ of copies of G with vertices indexed by j, the label in G, and c ∈ Z, the unit cell. This will give a new, infinite Laplacian L ∞ . For the edges that will not be rewired, we set L ∞ jc,kc = L jk for all c. Doing this for all edges would leave the copies of G disconnected. To connect them, we choose a subset of edges {(j, k)} and rewire them to cross between unit cells; for example, if (j, k) is an edge to be rewired to have k in a unit cell to the left of j we can set L ∞ jc,k(c−1) = L jk for all c and symmetrically set L ∞ k(c−1),jc = L kj . The remainder of the entries of L ∞ are set to zero.
Since L ∞ is periodic, Bloch's Theorem allows us to write the eigenvectors as where q is a wavenumber in the first Brilloun zone −π < q < π. TheŨ then satisfy which reduces to a new eigenvalue equation for a matrix of size n: where the matrix elements ofL(q) are the same as those of L for edges within a single unit cell and differ by a factor e iq(c−d) for edges that cross between unit cells c and d.
Using these transformations, which are standard in the study of lattice systems [17], we can find the continuous spectra of periodic tilings of our designed networks. Even without any optimization of which edges to rewire, the spectral characteristics persist in the infinite system. If we rewire all edges with |j − k| > n/2, for example, a spectrum of equally-spaced eigenvalues leads to a density of states with corresponding equally-spaced large spikes (Fig. 5a-c), while a discrete band gap is almost entirely preserved ( Fig. 5d-f). In both cases, only the bottom few bands vary significantly with q. Note that, because we moved edges incident to the first vertex, all of the eigenvectors do change and are not localized for nonzero q.

V. CONCLUSIONS
Controlling dynamics on a network typically requires detailed understanding of its spectral properties. Here we have reversed the conventional approach by starting from a desired spectrum and providing a mathematically rigorous construction of a matching network. This enabled us to induce chimera states, suppress pattern formation, and control wavefunction localization [53] using suitably designed gapped spectra. Our method, which starts from global properties, complements traditional approaches using small-scale local rules to build and analyze networks [29,[54][55][56]. In the future, the above results may also prove useful as a standard of comparison for other networks. Contrasting the dynamics on an important class of networks with the dynamics on networks designed to have identical spectra can help identify the important features of that class. Moreover, as dynamics are often related to matrices other than the Laplacian [57], it will be interesting to investigate control of their spectra for weighted networks as well. Although our construction works optimally with fully-connected graphs, one can expect that improved sparsification algorithms together with recent progress in 3D printing and lithography [18,24] may soon lead to physicallyrealizable networks with arbitrary gaps; since any graph can be embedded in 3D [58], the framework introduced here lays a conceptual foundation for the targeted design of complex non-periodic metamaterials with desired spectral properties. Currently, the approach can be applied directly to interfacing biochemical oscillators with computational networks [14]. Such hybrid networks can also be naturally extended to periodic systems, where the spectral properties are preserved well without any further optimization.
Acknowledgements. The authors would like to thank Jon Kelner and Philippe Rigollet for helpful discussions. This work was supported by Trinity College There are two types of edges: edges with both endpoints in the first m + 1 vertices have weight λ n−1 /n + (λ 1 − λ n−1 )/(m + 1), while other edges have weight λ n−1 /n.