Generating directed networks with prescribed Laplacian spectra

Complex networks play a role of paramount importance for a wide range of problems, of cross-disciplinary breadth. In several cases of interest, networks define the skeleton of pairwise interaction between coupled populations, families of homogeneous constituents anchored to a given node of the collections. The nature of existing paired relationships between mutually entangled populations is encoded in the weight of the links, that bridge adjacent nodes. Cooperative and competitive interactions can, in principle, be accommodated for by allowing for the weights to take positive or negative values. Local and remote non linear couplings shape the ensuing dynamics, and possibly steer the system towards a stationary stable equilibrium which is compatible with the assigned initial condition. The stability of the fixed points can be analyzed by studying the dynamical system in its linearized version. For reaction-diffusion systems defined on networks, the stability of the inspected equilibrium is ultimately dictated by the spectrum of the discrete Laplacian matrix [1–3]. The eigenvalues of the Laplacian define, in fact, the support of the dispersion relation, the curve that sets the rate for the exponential growth of the imposed perturbation. More specifically, external disturbances can be, in general, decomposed on the basis formed by the eigenvectors of the Laplacian operator. Each eigenvector defines an independent mode, which senses the web of intricate paths made accessible across the network: the perturbation can eventually develop, or, alternatively, fade away, along the selected direction, depending on the corresponding entry of the dispersion relation, as fixed by its associated eigenvalues. Stability is an attribute of paramount importance as it relates to resilience, the ability of a given system to oppose to external perturbations that would take it away from the existing equilibrium. Similarly, synchronization, a widespread phenomenon in distributed systems, can be enforced by properly adjusting the spectrum of the matrix which encodes for intertwined pairings. Based on the above, it is therefore essential to devise suitably tailored recipes for generating networks, which display a prescribed Laplacian spectrum, compatible with the stability constraint [4, 5].

The problem of recovering a network from a set of assigned eigenvalues has been tackled in the literature both from an algorithmic [6–8] and formal [9–11] standpoints. In [10], a procedure is discussed to generate an undirected and weighted graph from its spectrum. The result extends beyond the well-known theorem of Botti and Merris [12] which states that the reconstruction of non-weighted graphs is, in general, impossible since almost all (non-weighted) trees share their spectrum with another non-isomorphic tree. In [13], a method is proposed to obtain a, directed or undirected, graph whose eigenvalues are constrained to match specific bounds, which ultimately reflect the nodes degrees, as well as the associated weights. In [14], a mathematically rigorous strategy is instead developed to yield weighted graphs which exactly realize any desired spectrum. As discussed in [14], the method translates into an efficient approach to control the dynamics of various archetypal physical systems via suitably designed Laplacian spectra. The results are however limited to undirected Laplacians, characterized by a real spectrum. The purpose of this paper is to expand beyond these lines, by proposing and testing a procedure which enables one to recover a signed Laplacian operator which displays a prescribed complex spectrum. Signed Laplacians are often used in the literature for applications which relate to social contagion, cluster synchronization or repulsive–attractive interactions [19, 20]. In engineering, they are often employed in modeling microgrids dynamics [21].

The paper is organized as follows. The first section is devoted to illustrating the devised method, focusing on the mathematical aspects. We then turn to discussing the implementation of the scheme and introducing the sparsification algorithms that are run to cut unessential links. In the subsequent section, we elaborate on the conditions that are to be met to generate a positively weighted network. This discussion is carried out with reference to a specific setting. Then, we apply the newly introduced technique to the study of an ensemble made of coupled Stuart–Landau oscillators [22–24] and to (a simplified version of) the Kuramoto model [27, 28]. Finally, in the last section, we sum up the contributions and provide concluding remarks.

Consider a network made of Ω nodes and denote by A the (weighted) adjacency matrix, where structural information is encoded. More precisely, the element A_ij is different from zero when a directed link exists from j to i. The entries of the matrix A are real numbers and their signs reflect the specificity of the interaction at play: negative signs stand for inhibitory (or antagonistic) couplings, while positive entries point to excitatory (or cooperative) interaction. From the adjacency matrix, one can define its associated Laplacian operator. This is the matrix L, whose elements are L_ij = A_ij − k_i δ_ij , where k_i = ∑_j A_ij represents the natural extension of the concept of (incoming) connectivity to the case of a weighted network and δ_ij denotes the Kronecker delta.

We shall here discuss a procedure to generate a Laplacian matrix, which displays a prescribed set of eigenvalues. As anticipated above, we will focus in particular on directed Laplacians, which yields, in general, complex spectra. Concretely, we begin by introducing a collection of Ω = 2N + 1 complex quantities defined as ⁴

$$\begin{equation}\left\{{{\Lambda}}_{i}\right\}=\left\{{{\Lambda}}_{1},{{\Lambda}}_{2},\dots ,{{\Lambda}}_{2N},{{\Lambda}}_{2N+1}=0\right\}\end{equation} \tag{ 1 }$$

The first 2N elements come in complex conjugate pairs and we set in particular ${{\Lambda}}_{i}={{\Lambda}}_{i+N}^{{\ast}}$, ∀ i = 1, ..., N, where (⋅)* stands for complex conjugate. The aim of this section is to develop a rigorous procedure to construct a directed (and weighted) graph G with 2N + 1 vertices, whose associated Laplacian has {Λ_i } for eigenvalues. Recall that {Λ_i } contains the null element, since this latter is, by definition, a Laplacian's eigenvalue.

The procedure that we are going to detail in what follows exploits the eigenvalue decomposition of the Laplacian matrix. To this end we will seek to introduce a proper eigenvector basis such that

$$\begin{equation}L=VD{V}^{-1}\end{equation} \tag{ 2 }$$

is a Laplacian. In equation (2), D is a diagonal matrix where the sought Laplacian eigenvalues are stored. More specifically, D_ii = Λ_i−1 for i = 2, ..., 2N + 1 and D₁₁ = Λ_2N+1 = 0. The problem is hence traced back to constructing V, whose columns are the right eigenvectors of L. We also recall that rows of the inverse matrix V⁻¹ are the left eigenvectors of L. As outlined in [15], Laplacian (right and left) eigenvectors should satisfy a set of conditions:

• (a)  
  The columns of V, which refer to complex conjugate eigenvalues, must be complex conjugate too.
• (b)  
  The same condition holds for the rows of V⁻¹.
• (c)  
  Moreover, the columns of V (resp. the rows of V⁻¹) corresponding to eigenvalues different from 0 should sum up to zero.
• (d)  
  Finally, the right eigenvector relative to the null eigenvalue should be uniform (i.e. display identical components).

In light of the above, we put forward for V the following structure:

$$\begin{equation}V=\left(\begin{matrix}\hfill c\hfill & \hfill {v}^{\mathrm{T}}\hfill & \hfill {{v}^{\mathrm{T}}}^{{\ast}}\hfill \\ \hfill cu\hfill & \hfill \mathrm{i}U\hfill & \hfill -\mathrm{i}U\hfill \\ \hfill cu\hfill & \hfill U\hfill & \hfill U\hfill \end{matrix}\right)\end{equation} \tag{ 3 }$$

where c is an arbitrary constant, i stands for the imaginary unit, U is an invertible N × N matrix having real entries. Further, vector u = (1...1)^T has dimension N × 1 and vector v is defined as

$$\begin{equation}{v}^{\mathrm{T}}=-\left(1+\mathrm{i}\right){u}^{\mathrm{T}}U\end{equation} \tag{ 4 }$$

The first column of V is hence a uniform vector, corresponding to the eigenvector associated to the null eigenvalue. By construction, every other column sums up to zero, that is (4) holds. In the following, we will write D_jj = α_j + iβ_j , which, in turn, implies D_j+N,j+N = α_j − iβ_j , for j = 2, ..., N + 1. Here, α_j and β_j are real quantities and respectively denote the real and imaginary parts of the jth eigenvalue. To proceed further, one needs to determine the inverse of V.

To achieve this goal we begin by considering a generic matrix W, which satisfies the general constraints that are in place for V⁻¹. In formulae:

$$\begin{equation}W=\left(\begin{matrix}\hfill d\hfill & \hfill d{u}^{\mathrm{T}}\hfill & \hfill d{u}^{\mathrm{T}}\hfill \\ \hfill \left(1-\mathrm{i}\right)w\hfill & \hfill S\hfill & \hfill -\mathrm{i}{S}^{{\ast}}\hfill \\ \hfill \left(1+\mathrm{i}\right){w}^{{\ast}}\hfill & \hfill {S}^{{\ast}}\hfill & \hfill \mathrm{i}S\hfill \end{matrix}\right)\end{equation} \tag{ 5 }$$

where

$$\begin{equation}w=-\frac{1}{1-\mathrm{i}}\left(Su-\mathrm{i}{S}^{{\ast}}u\right)\end{equation} \tag{ 6 }$$

Note that the jth and (N + j)th rows of W, for j = 1, ..., N + 1, are complex conjugated, as required. Moreover, summing all the elements of each row (but the first) yields zero, a condition that the inverse of V should meet, as anticipated above. Building on these premises, we shall here determine the unknown S, w and d so as to match the identity WV = I, where I stands for the (2N + 1) × (2N + 1) identity matrix. This implies, in turn, that W ≡ V⁻¹ due to the uniqueness of the inverse matrix.

A straightforward manipulation yields the following conditions for, respectively, d and w

$$\begin{equation}\begin{cases}d=\frac{1}{c\left(2N+1\right)}\quad \hfill \\ w{v}^{\mathrm{T}}=\mathrm{i}\left(S-\mathrm{i}{S}^{{\ast}}\right)EU\quad \hfill \\ w{{v}^{\mathrm{T}}}^{{\ast}}=\left(S-\mathrm{i}{S}^{{\ast}}\right)EU\quad \hfill \end{cases}\end{equation} \tag{ 7 }$$

where use has been made of the identity u^T u = N and where

$$\begin{equation}E=u{u}^{\mathrm{T}}=\left(\begin{matrix}\hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \\ \hfill {\vdots}\hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill 1\hfill & \hfill \dots \hfill & \hfill 1\hfill \end{matrix}\right)\end{equation} \tag{ 8 }$$

The quantity d is completely specified by the first of equations (7) and solely depends on c and N, the size of the system. By making use of the identities (4) and (6), one can progress in the analysis of the second and third conditions (7) to eventually get:

$$\begin{equation}SB+{S}^{{\ast}}A=I\end{equation} \tag{ 9 }$$

$$\begin{equation}SA-{S}^{{\ast}}B=0\end{equation} \tag{ 10 }$$

where:

$$\begin{equation}A=EU-\mathrm{i}EU-\mathrm{i}U\end{equation} \tag{ 11 }$$

$$\begin{equation}B=\mathrm{i}EU+EU+\mathrm{i}U\end{equation} \tag{ 12 }$$

It is, therefore, immediate to conclude that

$$\begin{equation}S={\left(B+A{B}^{-1}A\right)}^{-1}\end{equation} \tag{ 13 }$$

$$\begin{equation}{S}^{{\ast}}={\left(B+A{B}^{-1}A\right)}^{-1}A{B}^{-1}\end{equation} \tag{ 14 }$$

The analysis can be pushed further to relate S to matrices U and E. The calculation, detailed in appendix A, yields

$$\begin{equation}S=-\frac{\mathrm{i}}{2}{U}^{-1}\left[I-\frac{1+\mathrm{i}}{1+2N}E\right]\end{equation} \tag{ 15 }$$

The expression for S* can be immediately obtained by taking the complex conjugate of the above equation. In conclusion, the matrix W defined in (5) is the inverse matrix of V, provided that d and S are respectively assigned as specified above.

Clearly, matrix L defined in (2) has the desired spectrum (1). We should, however, prove that L is a Laplacian. This amounts to showing that L is a real matrix, whose columns sum up to zero. The proof is given hereafter.

Proposition. Matrix L is real.

From equations (3) and (5), one can readily compute the elements of L via matrix products and, taking advantage of the block structures of V and W, prove that L is real. $\mathfrak{R}\left(\cdot \right)$ is introduced to represent the real part of (⋅). The generic element L_st can be written as:

$$\begin{equation}{L}_{st}={\left(VDW\right)}_{st}=\sum _{k=1}^{2N+1}{V}_{sk}{\left(DW\right)}_{kt}=\sum _{k=1}^{2N+1}{V}_{sk}\sum _{{k}^{\prime }=1}^{2N+1}{D}_{k{k}^{\prime }}{W}_{{k}^{\prime }t}\end{equation} \tag{ 16 }$$

$$\begin{equation}=\sum _{k=2}^{2N+1}{V}_{sk}{D}_{kk}{W}_{kt}\end{equation} \tag{ 17 }$$

due to the diagonal structure of the matrix D and recalling that D₁₁ = 0. By making use of the specific form of V and W, one gets:

$$\begin{equation}{L}_{st}=\sum _{k=2}^{N+1}{V}_{sk}{D}_{kk}{W}_{kt}+\sum _{k=N+2}^{2N+1}{V}_{sk}{D}_{kk}{W}_{kt}\end{equation} \tag{ 18 }$$

$$\begin{equation}=\sum _{k=2}^{N+1}{V}_{sk}{D}_{kk}{W}_{kt}+\sum _{k=2}^{N+1}{\left({V}_{sk}{D}_{kk}{W}_{kt}\right)}^{{\ast}}\end{equation} \tag{ 19 }$$

$$\begin{equation}=\sum _{k=2}^{N+1}2\mathfrak{R}\left({V}_{sk}{D}_{kk}{W}_{kt}\right)\end{equation} \tag{ 20 }$$

since $\alpha \beta +{\alpha }^{{\ast}}{\beta }^{{\ast}}=2\mathfrak{R}\left(\alpha \beta \right)$, for any complex numbers α and β. One can thus conclude that L, as generated by the above procedure, is real.

Proposition. Each column of L sums up to zero.

Because of the diagonal structure of D:

$$\begin{equation}\sum _{i}{L}_{ij}=\sum _{il}{V}_{il}{D}_{ll}{V}_{lj}^{-1}\end{equation} \tag{ 21 }$$

Then:

$$\begin{equation}\sum _{i}{L}_{ij}=\sum _{l}{V}_{lj}^{-1}{D}_{ll}\sum _{i}{V}_{il}=0\end{equation} \tag{ 22 }$$

since (i) D₁₁ = 0 and (ii) the components of all eigenvectors corresponding to non-null eigenvalues, sum up to zero. Notice that this result can also be proven by observing that the uniform vector d 1 is the left eigenvector corresponding to the null eigenvalue, that is

$$\begin{equation}d\mathbf{1}L=0\end{equation} \tag{ 23 }$$

From (23), it follows that ∑_i L_ij = 0, for every j.

Proposition.  L is balanced.

We can also show that ∑_j L_ij = 0 i.e. that the sum of all the elements of any given row i returns zero. According to (2), the first column of V is the right eigenvector corresponding to eigenvalue 0, namely

$$\begin{equation}Lc\mathbf{1}=cL\mathbf{1}=0\end{equation} \tag{ 24 }$$

This implies, in turn, ∑_j L_ij = 0 ∀ i, which ends the proof. The Laplacian is hence balanced, as the sums on the rows and on the (corresponding) columns return the same result.

From the generated Laplacian operator, one can readily calculate the adjacency matrix of the underlying network. In general, for any assigned spectrum, the recovered adjacency matrix is fully connected, meaning that there exist links connecting each pair of nodes. Notice that links are weighted and signed. The weights can be small or have a modest impact on the spectrum of the associated Laplacian. This motivates the implementation of a dedicated sparsification procedure, which seeks to remove unessential links, in terms of their reflection on the ensuing Laplacian spectrum. The next section is devoted to elaborating along these lines.

In this section we discuss a sparsification procedure, which aims at a posteriori simplifying the structure of the recovered network [16-18]. To this end, we begin by generating a network following the strategy outlined in the preceding section, and which yields an assigned spectrum for the associated Laplacian. The Laplacian spectrum that we seek to recover consists of Ω = 2N + 1 complex entries, the eigenvalues, which are here confined in the left portion of the complex plane, by setting $\mathfrak{R}\left({{\Lambda}}_{j}\right)={\alpha }_{j}{< }0$ for j ⩾ 2, see blue crosses in figure 1(a). This choice is somehow arbitrary, and ultimately amounts to enforce stability into a linear system of the form:

$$\begin{equation}\frac{\mathrm{d}{x}_{i}}{\mathrm{d}t}=\sum _{j}{L}_{ij}{x}_{j}\end{equation} \tag{ 25 }$$

where x_i is the ith entry of the Ω-dimensional state vector x. In the final part of the paper, we will turn to considering more complex scenarios where the stability of the examined dynamics is also influenced by local reaction terms.

Zoom In Zoom Out Reset image size

The network that we obtain, following the scheme outlined in the preceding sections yielding a Laplacian with the prescribed spectrum, is in general fully connected. In other words, a weighted link exists between any pair of nodes. The weights of the link can be, in principle, very small and, as such, bear a modest imprint on the ensuing Laplacian spectrum. Motivated by this observation, we perform an a posteriori sparsification of the obtained network: this aims to identifying and then removing the finite subset of links that appear to have a modest impact on the eigenvalues of the associated Laplacian.

The first sparsification procedure that we have considered, aims at removing unessential links while confining the spectrum of the Laplacian operator within a bounded region of the complex plane. More precisely, we focus on the links which display a weight in the range (−σ, σ), where σ is a small, arbitrarily chosen, cut off. All links whose weight is smaller that σ in absolute value are selected, in a random order. The selected link is removed and the modified Laplacian spectrum computed. Denote by ${\tilde {{\Lambda}}}_{j}$, for j = 2, ..., 2N + 1, the Laplacian eigenvalues obtained upon removal of the link. The change to the network arrangement becomes permanent, if $\vert {\mathrm{m}\mathrm{i}\mathrm{n}}_{j}\left[\mathfrak{R}\left({\tilde {{\Lambda}}}_{j}\right)\right]-{\mathrm{m}\mathrm{i}\mathrm{n}}_{j}\left[\mathfrak{R}\left({{\Lambda}}_{j}\right)\right]\vert {< }\delta$ and $\vert {\mathrm{m}\mathrm{a}\mathrm{x}}_{j}\left[\mathfrak{I}\left({\tilde {{\Lambda}}}_{j}\right)\right]-{\mathrm{m}\mathrm{a}\mathrm{x}}_{j}\left[\mathfrak{I}\left({{\Lambda}}_{j}\right)\right]\vert {< }\delta$, for j = 2, ..., N. Here $\mathfrak{I}\left(\cdot \right)$ stands for the imaginary part of (⋅) and δ is an arbitrary threshold which quantifies the amount of perturbation that is deemed acceptable for the problem at hand. As a further condition, we check that $\mathfrak{R}\left({\tilde {{\Lambda}}}_{j}\right){< }0$ for j = 2, ..., 2N + 1, which, in turn, corresponds to preserving the stability of the linear system (25). Clearly, the order of extraction of the links, which are candidate to be trimmed, matters. Different realizations of the procedure of progressive sparsification illustrated above might hence result in distinct final outcomes. In figure 1(a), the eigenvalues obtained after the sparsification algorithm are plotted (red circles) for one choice of the cutoff δ. The sparsity pattern of the adjacency matrix obtained at the end of the above procedure is displayed in panel (b) of figure 1. In panels (c) and (d) of figure 1 we plot, with an appropriate color code, the entries of the adjacency matrices, before and after the sparsification. Only weights which are significantly different from zero (see annexed colorbars) are displayed. As appreciated by visual inspection, the skeleton of the network is not altered by the applied sparsification. To monitor how the eigenvalues get redistributed within the bounded domain to which they belong, we introduce the following indicators:

$$\begin{equation}{I}_{x}=\sum _{i=2}^{2N+1}{\left({\beta }_{i}\right)}^{2}\end{equation} \tag{ 26 }$$

$$\begin{equation}{I}_{y}=\sum _{i=2}^{2N+1}{\left({\alpha }_{i}-\frac{1}{2N+1}\sum _{j=2}^{2N+1}{\alpha }_{j}\right)}^{2}\end{equation} \tag{ 27 }$$

The quantity I_x measures the dispersion along the imaginary axis, by weighting the squared distance of each eigenvalue from the horizontal axis. Conversely, I_y reflects the scattering of the eigenvalues about their mean, in the direction of the real axis. In figure 2, I_x and I_y , normalized to their respective values obtained before application of the sparsification algorithm, are shown versus N, an indicator of the size of the generated networks. The sparsification procedure shrinks the eigenvalues in the x-direction, while the opposite tendency is observed for the distribution along the y-direction.

Zoom In Zoom Out Reset image size

The second sparsification method implements a more stringent constraint. Just like before, we select the links with weights in the range (−σ, σ), where σ acts as a small threshold amount. Unlike the former case, we now eliminate the selected link only if the change produced on the modulus of each of the N eigenvalues is smaller than δ, namely if $\vert {{\Lambda}}_{j}-{\tilde {{\Lambda}}}_{j}\vert {< }\delta$, for j = 2, ..., 2N + 1. In figure 3, the eigenvalues obtained after the sparsification algorithm are plotted (red circles) for two choices of the cutoff δ. The number of links that can be effectively removed grows with δ, the size of the allowed perturbation, as clearly demonstrated in figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned above, this latter method implements a local, hence stringent, constraint. By allowing modest re-arrangements of individual eigenvalues yields a spectral distribution in the complex plane which is by construction close to the original one. This method is better suited when the room for maneuver is limited, due to the specific constraints that the reactive model implements. On the other hand the number of pruned links can be scanty, when small local perturbations are solely allowed for. A global re-organization of the operator's spectrum, which entails a significant reduction in the associated number of active links, might instead be accomplished when the eigenvalues are constrained to fall within a compact domain in the complex plane, as follows the former sparsification procedure.

Summing up, we have developed and tested a procedure to generate a network which returns an associated Laplacian matrix with a prescribed complex spectrum. The weighted network obtained following the above procedure is, in general, fully connected. Dedicated sparsification strategies can be applied to remove the links which carry a small weight, and bear a modest imprint on the ensuing Laplacian spectrum. In the following, we will consider a specific setting of the aforementioned generation scheme, which makes it possible for the Laplacian elements to be computed analytically.

In the previous sections we described a general method to generate a Laplacian matrix which displays a designated spectrum. The method assumes a generic matrix U, which can be randomly assigned. In the following, we will focus on the specific case where U is proportional to the identity matrix and progress with the analytic characterization of the obtained Laplacian. As we shall argue in the following, working in this framework allows us to derive a set of closed conditions for constraining the weights of the underlying network to strictly positive values. To proceed in this direction we set:

$$\begin{equation}U=qI\end{equation} \tag{ 28 }$$

where I stands for the identity matrix and q is scalar.

A straightforward calculation returns the following expression for matrix S:

$$\begin{equation}S=\left(\begin{matrix}\hfill a+\mathrm{i}b\hfill & \hfill a-\mathrm{i}a\hfill & \hfill \dots \hfill & \hfill a-\mathrm{i}a\hfill \\ \hfill a-\mathrm{i}a\hfill & \hfill a+\mathrm{i}b\hfill & \hfill \ddots \hfill & \hfill a-\mathrm{i}a\hfill \\ \hfill \dots \hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill {\vdots}\hfill \\ \hfill a-\mathrm{i}a\hfill & \hfill \ddots \hfill & \hfill \ddots \hfill & \hfill a+\mathrm{i}b\hfill \end{matrix}\right)\end{equation} \tag{ 29 }$$

while w is a uniform vector with identical entries equal to Na + ((N − 1)a − b)i and the quantities d, a, b are specified by

$$\begin{equation}d=\frac{1}{\left(2N+1\right)c}\end{equation} \tag{ 30 }$$

$$\begin{equation}a=\frac{-1}{2\left(2N+1\right)q}\end{equation} \tag{ 31 }$$

$$\begin{equation}b=2Na\end{equation} \tag{ 32 }$$

From equations (16) and (18) one can obtain a closed expression for each element of the Laplacian, as function of the eigenvalues. The interested reader can find the detailed computations in appendix B. In the following the final formulae are reported.

The diagonal elements satisfy:

$$\begin{equation}{L}_{11}=\frac{2}{2N+1}\sum _{k=2}^{N+1}{\alpha }_{k}\end{equation} \tag{ 33 }$$

$$\begin{equation}{L}_{ss}=\frac{2N{\alpha }_{s}+{\beta }_{s}}{2N+1}\quad s=2,\dots ,N+1\end{equation} \tag{ 34 }$$

$$\begin{equation}{L}_{ss}=\frac{2N{\alpha }_{s}-{\beta }_{s}}{2N+1}\quad s=N+2,\dots ,2N+1,\end{equation} \tag{ 35 }$$

where use has been made of the identity $\mathfrak{R}\left({D}_{kk}\right)={\alpha }_{k}$. The first row and column are given by:

$$\begin{equation}{L}_{1t}=\frac{1-2N}{1+2N}{\alpha }_{t}-{\beta }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad t=2,\dots ,N+1\end{equation} \tag{ 36 }$$

$$\begin{equation}{L}_{1t}=\frac{1-2N}{1+2N}{\alpha }_{t}+{\beta }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t-N}^{N+1}{\alpha }_{k}\quad t=N+2,\dots ,2N+1\end{equation} \tag{ 37 }$$

$$\begin{equation}{L}_{t1}=\frac{-{\alpha }_{t}+{\beta }_{t}}{2N+1}\quad t=2,\dots ,N+1\end{equation} \tag{ 38 }$$

$$\begin{equation}{L}_{t1}=\frac{-{\alpha }_{t}-{\beta }_{t}}{2N+1}\quad t=N+2,\dots ,2N+1.\end{equation} \tag{ 39 }$$

while the remaining elements are:

$$\begin{equation}{L}_{st}=\frac{-{\alpha }_{s}+{\beta }_{s}}{2N+1}\quad s=2,\dots ,N+1\enspace \;\text{and}\enspace t=N+2,\dots ,2N+1\enspace \text{with}\enspace t\ne s-N\end{equation} \tag{ 40 }$$

$$\begin{equation}{L}_{st}=\frac{-{\alpha }_{s}-{\beta }_{s}}{2N+1}\quad s=N+2,\dots ,2N+1\;\text{and}\;t=2,\dots ,N+2\enspace \text{with}\enspace t\ne s-N\end{equation} \tag{ 41 }$$

$$\begin{equation}{L}_{s,s+N}=\frac{-{\alpha }_{s}-2N{\beta }_{s}}{2N+1}\quad s=2,\dots ,N+1\end{equation} \tag{ 42 }$$

$$\begin{equation}{L}_{s,s-N}=\frac{-{\alpha }_{s}+2N{\beta }_{s}}{2N+1}\quad s=N+2,\dots ,2N+1\end{equation} \tag{ 43 }$$

$$\begin{equation}{L}_{st}=\frac{-{\alpha }_{s}+{\beta }_{s}}{2N+1}\quad t,s=2,\dots ,N+1\;\text{and}\;s\ne t\end{equation} \tag{ 44 }$$

$$\begin{equation}{L}_{st}=\frac{-{\alpha }_{s}-{\beta }_{s}}{2N+1}\quad s,t=N+2,\dots ,2N+1\;\text{and}\;s\ne t\end{equation} \tag{ 45 }$$

4.1. Controlling the sign of non-diagonal Laplacian entries

The Laplacian matrix obtained with the procedure illustrated above has both positive and negative entries. Signed Laplacians are often used in consensus problems, where negative weights model antagonistic interactions. In other contexts, when e.g. the Laplacian is stemming from diffusive interactions, non diagonal entries are constrained to positive values. In the following, we will provide a set of necessary conditions for the assigned spectrum to eventually yield a Laplacian with positive extra diagonal elements, L_ij > 0 for i ≠ j. Clearly, L_ii < 0, as summing on the rows should return zero. The underlying network, therefore, displays positive weights, as its adjacency matrix is basically obtained from the Laplacian matrix by replacing the diagonal elements with zeros.

Further, we will set $\mathfrak{R}\left({{\Lambda}}_{j}\right)={\alpha }_{j}{< }0$ for j ⩾ 2, an assumption which corresponds to dealing with a stable linear system of the type given in equation (25). This requirement immediately yields L₁₁ < 0, as it follows from relation (33). We will also operate in the setting analyzed above, i.e. assuming U = qI. The obtained expressions for the Laplacian elements allow to recast the sought conditions on their signs as:

$$\begin{equation}\begin{cases}2N{\alpha }_{t}{< }{\beta }_{t}{< }-2N{\alpha }_{t}\quad \hfill \\ {\beta }_{t}{ >}-\frac{1-2N}{2N+1}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\beta }_{t}{< }\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\alpha }_{t}{< }{\beta }_{t}{< }-{\alpha }_{t}\quad \hfill \\ \frac{{\alpha }_{t}}{2N}{< }{\beta }_{t}{< }-\frac{{\alpha }_{t}}{2N}\quad \hfill \end{cases}\end{equation} \tag{ 46 }$$

where the inequalities hold for t = 2, ..., N + 1.

The above conditions can be simplified, after some algebraic manipulations, to return

$$\begin{equation}\frac{4N}{4{N}^{2}-1}\sum _{k}{\alpha }_{k}{< }{\alpha }_{t}{< }\frac{2}{2N+1}\sum _{k}{\alpha }_{k}\end{equation} \tag{ 47 }$$

and

$$\begin{equation}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k}{\alpha }_{k}{< }{\beta }_{t}{< }-{\alpha }_{t}+\frac{2}{2N+1}\sum _{k}{\alpha }_{k}.\end{equation} \tag{ 48 }$$

The interested reader can access the detailed steps in appendix C.

Assume that the assigned Laplacian spectrum matches the above conditions, while having α_k < 0 for k > 2. Then the Laplacian matrix obtained with the above procedure, with U = qI, displays positive non diagonal entries.

Let us explore the consequences of conditions (47) and (48). To this end, introduce $\bar{\alpha }$, the average of the non-negative real parts of the Laplacian eigenvalues, i.e. $\bar{\alpha }=\left({\sum }_{k}{\alpha }_{k}\right)/\left(2N\right)$. A straightforward analysis allows us to conclude that the generated Laplacian returns positive non-diagonal elements, if the assigned non trivial eigenvalues (i.e. Λ_k , with k > 1) fall in a bounded rectangular domain of the complex plane. More specifically, the rectangular region is symmetric, with respect to the horizontal (real) axis, and extends along the vertical direction (imaginary axis) of ${\pm}2N\bar{\alpha }/\left(4{N}^{2}-1\right)$. The rectangle is completed by two vertical sides, positioned at $\frac{4{N}^{2}}{4{N}^{2}-1}\bar{\alpha }$  and  $\frac{2N}{2N+1}\bar{\alpha }$. Working at fixed N, the larger is $\vert \bar{\alpha }\vert$ the more extended is the rectangle along the vertical direction. Conversely, when making N larger, the rectangle shrinks in the horizontal direction and becomes eventually degenerate for N → ∞. In other words, for large values of N, eigenvalues should align on a vertical segment positioned at $\bar{\alpha }$, and whose extension increases linearly with $\vert \bar{\alpha }\vert$ (while decreasing with N). This is shown in figure 5, where three different spectra are depicted (for different choices of $\bar{\alpha }$) which yield a Laplacian with positive off-diagonal elements.

Zoom In Zoom Out Reset image size

In the following section, we will apply the method of Laplacian generation here discussed to the study of two prototypical examples of dynamical systems on networks.

In this section we consider two different models of interacting oscillators, defined on a network. In both cases, the coupling between individual oscillators is implemented via a discrete Laplacian operator, which reflects the specific network arrangement. It will be shown that a suitable network arrangement can be a priori established, building on the procedure illustrated above, so as to make the inspected systems stable against external perturbations.

5.1. Coupled Stuart–Landau oscillators

Consider an ensemble made of 2N + 1 nonlinear oscillators and denote by W_i their associated complex amplitude. We assume the oscillators to be mutually coupled via a diffusive-like interaction which is mathematically modeled by a discrete Laplacian operator. Each oscillator obeys a complex Stuart–Landau equation [25]. The dynamics of the system can be cast in the form:

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{W}_{j}={W}_{j}-\left(1+\mathrm{i}{c}_{2}\right)\vert {W}_{j}{\vert }^{2}{W}_{j}+\left(1+\mathrm{i}{c}_{1}\right)K\sum _{k}{L}_{jk}{W}_{k}\end{equation} \tag{ 49 }$$

where c₁ and c₂ are real parameters. The index j runs from 1 to 2N + 1, the total number of oscillators. Here, K is a suitable parameter setting the coupling strength. Without loss of generality, in what follows it is assumed that K = 1. L_ij = A_ij − k_i δ_ij is the Laplacian, A_ij is the generic entry of the directed and weighted adjacency matrix A and k_i = ∑_j A_ij .

The system admits a homogeneous limit cycle solution in the form ${W}_{LC}\left(t\right)={\mathrm{e}}^{-\mathrm{i}{c}_{2}t}$. To characterize the stability of the cycle [26], one can introduce a non homogeneous perturbation in polar coordinates as:

$$\begin{equation}{W}_{i}\left(t\right)={W}_{LC}\left[1+{\rho }_{i}\left(t\right)\right]{\text{e}}^{\text{i}{\theta }_{i}\left(t\right)}\end{equation} \tag{ 50 }$$

By linearizing around the limit cycle solution (ρ_i (t) = 0, θ_i (t) = 0), one gets:

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix}\hfill {\rho }_{j}\hfill \\ \hfill {\theta }_{j}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill -2\hfill & \hfill 0\hfill \\ \hfill -2{c}_{2}\hfill & \hfill 0\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\rho }_{j}\hfill \\ \hfill {\theta }_{j}\hfill \end{matrix}\right)+\left(\begin{matrix}\hfill 1\hfill & \hfill -{c}_{1}\hfill \\ \hfill {c}_{1}\hfill & \hfill 1\hfill \end{matrix}\right)\sum _{k}{L}_{jk}\left(\begin{matrix}\hfill {\rho }_{k}\hfill \\ \hfill {\theta }_{k}\hfill \end{matrix}\right)\end{equation} \tag{ 51 }$$

To proceed further, expand the perturbations ρ_j and θ_j on the Laplacian eigenvectors v_(α), that is

$$\begin{equation}\left(\begin{matrix}\hfill {\rho }_{j}\hfill \\ \hfill {\theta }_{j}\hfill \end{matrix}\right)=\sum _{\alpha =1}^{2N+1}\left(\begin{matrix}\hfill {\rho }^{\left(\alpha \right)}\hfill \\ \hfill {\theta }^{\left(\alpha \right)}\hfill \end{matrix}\right){\mathrm{e}}^{{\lambda }_{\alpha }t}{\left({v}_{\left(\alpha \right)}\right)}_{i}\end{equation} \tag{ 52 }$$

By inserting this expansion in (51) and using the relation

$$\begin{equation}\sum _{j}{L}_{ij}{\left({v}_{\left(\alpha \right)}\right)}_{j}={{\Lambda}}^{\left(\alpha \right)}{\left({v}_{\left(\alpha \right)}\right)}_{i}\end{equation} \tag{ 53 }$$

for α = 1, ..., 2N + 1, we obtain a condition formally equivalent to the expression of the continuous dispersion relation

$$\begin{equation}{\lambda }_{\mathrm{max}}\left({{\Lambda}}^{\left(\alpha \right)}\right)=-{{\Lambda}}^{\left(\alpha \right)}-1+\sqrt{-{c}_{1}^{2}{\left({{\Lambda}}^{\left(\alpha \right)}\right)}^{2}-2{c}_{1}{c}_{2}{{\Lambda}}^{\left(\alpha \right)}+1}\end{equation} \tag{ 54 }$$

If ${\lambda }_{\mathfrak{R}}=\mathfrak{R}\left({\lambda }_{\mathrm{max}}\right)$ is positive for some Λ^(α), the perturbation grows exponentially in time, and the initial homogeneous state turns out to be unstable. Conversely, if ${\lambda }_{\mathfrak{R}}=\mathfrak{R}\left({\lambda }_{\mathrm{max}}\right){< }0$, for every Λ^(α), the perturbation gets re-absorbed and the system converges back to the fully synchronized state. The condition ${\lambda }_{\mathfrak{R}}{< }0$ can be further processed analytically, as discussed in [15]. In particular, it can be shown that the latter condition is fulfilled, if the Laplacian eigenvalues fall in a specific portion of the parameter plane, which reflects the choice made for the reaction parameters c₁, c₂ and K. The region of interest is the one enclosed between the two solid lines, displayed in figure 6(a) for the specific selection of the parameters. The blue symbols depicted in figure 6(a) are randomly generated so as to fall in the region of the complex plane where stability holds. They represent the spectrum of the Laplacian that we seek to recover following the method illustrated above. In figure 6(b), ${\lambda }_{\mathfrak{R}}$ is plotted against $-{{\Lambda}}_{\mathfrak{R}}=-\mathfrak{R}\left({\Lambda}\right)$, confirming the stability of the homogeneous solution.

Zoom In Zoom Out Reset image size

We now proceed by generating a Laplacian matrix, which is constructed so as to yield the spectrum depicted in figure 6(a), see crosses (in blue). Circles (in red) refer to the spectrum obtained after the sparsification. In this case we operated with the more stringent sparsification protocol, the second of those introduced above. This latter appears more suited to the case at hand because of the shape of the domain where the eigenvalues should eventually fall to enforce stability. As an additional condition to the algorithm, we accepted only the sparsification moves which prevent the eigenvalues to invade the region of the complex plane associated to instability.

The sparsified adjacency matrix A defines the network of interactions between coupled oscillators, as follows equation (49). We hence integrate numerically the governing equations, assuming the initial state to be a perturbation of the homogenous synchronized equilibrium. As expected, the perturbation fades away and the system regains its unperturbed, fully synchronized, equilibrium.

5.2. Coupled Kuramoto oscillators

As a second example we set to study the Kuramoto model. Consider a system made of 2N + 1 oscillators, denote by θ_i the phase of the ith oscillator, and ω_i its natural frequency. The oscillators evolve as dictated by the following system of 2N + 1 coupled differential equations:

$$\begin{equation}\dot {{\theta }_{i}}={\omega }_{i}+\sum _{j=1}^{2N+1}{A}_{ij}\enspace \mathrm{sin}\left({\theta }_{j}-{\theta }_{i}\right)\quad i=1,{\ldotp\ldotp}\;,2N+1\end{equation} \tag{ 55 }$$

Here, A_ij stands for the entries of the adjacency matrix A which sets the interactions between pairs of oscillators. The matrix is, in principle, weighted, and may display positive and negative entries as reflecting the specific interaction (excitatory or inhibitory) being at play.

As an additional assumption, we will here focus in the simplified setting where ω_i = ω ∀ i. We can then introduce the new variable ψ_i = θ_i − ωt, and write the governing equation in the equivalent form:

$$\begin{equation}\dot {{\psi }_{i}}=\sum _{j=1}^{2N+1}{A}_{ij}\enspace \mathrm{sin}\left({\psi }_{j}-{\psi }_{i}\right)\quad i=1,{\ldotp\ldotp}\;,2N+1\end{equation} \tag{ 56 }$$

A homogeneous solution always exists with ψ_i = Ψ ∀ i, and for any constant Ψ ∈ [0, 2π), as it can be immediately checked by substitution. To assess the stability of the solution, one sets ψ_i = Ψ + δ_i , and expands (56) at the leading order in the δ_i . In this way, one gets:

$$\begin{equation}\dot {{\delta }_{i}}=\sum _{j=1}^{2N+1}{A}_{ij}\left({\delta }_{j}-{\delta }_{i}\right)=\sum _{j=1}^{2N+1}{L}_{ij}{\delta }_{j}\end{equation} \tag{ 57 }$$

where L_ij = A_ij − k_i δ_ij is the Laplacian operator which stems from A_ij . The stability of the simplified Kuramoto model here considered is controlled by a linear system of the type introduced in (25), with the obvious replacement of x_i with δ_i . The system proves hence stable if the (non trivial) eigenvalues of the Laplacian operator display negative real parts. Our aim, here, is to generate a Laplacian (and therefore a matrix of binary weighted connections among oscillators) which warrants the stability of the system. To this end we assign the eigenvalues (which appear in conjugate pairs) to belong to the negative portion of the complex plane, see (blue) crosses in figure 7(a). The null eigenvalue is clearly included into the spectrum. Running the procedure discussed in the first part of the paper, we obtain the corresponding Laplacian and compute the associated adjacency matrix. The first of the two sparsification procedures is then applied to remove unessential links, yielding the spectrum of the associated Laplacian depicted with (red) circles in figure 7(a). The Kuramoto model (56) is then integrated numerically by assuming the recovered expression for A. As predicted, the system is stable to external perturbations as one can clearly appreciate by inspection of figure 7(b).

Zoom In Zoom Out Reset image size

Studying the dynamics of an ensemble made of interacting units on a network is central for a large plethora of applications. In many cases of interests, individual units evolve under the influence of homologous constituents, the interaction stemming in general from binary exchanges. Distinct fundamental units are assigned to different nodes of the collection, paired via physical or virtual links. For a relevant subclass of problems, the stability of the ensuing equilibrium can be traced back to the spectrum of the Laplacian operator, computed from the adjacency matrix which defines the network arrangement. Symmetric networks yield Laplacian operators with real spectrum, while directionality in the couplings reflects in an imaginary spectrum. Methods exist which allow one to generate a symmetric network, hence Laplacian, with a prescribed real spectrum. Starting from these premises, we have proposed and tested a novel procedure to generate a (signed and directed) Laplacian which returns an a priori assigned (complex) spectrum. A special case has also been considered, which enables one to recover closed analytical expressions for the entries of the sought Laplacian. Working in this setting, we can elaborate on the conditions that have to be matched for the ensuing Laplacian to solely display positive non diagonal elements. Dedicated sparsification procedures are also discussed to help removing unessential links in terms of their impact on the associated spectrum. The algorithm for Laplacian generation has been successfully tested with reference to two prototypical examples of coupled oscillators. Taken together, our work explores possible strategies for network generation with the emphasis placed on dynamical, rather than structural features. The dynamics is indirectly modulated by the spectrum of the Laplacian operator, which is here constraining the generation algorithm.

All data that support the findings of this study are included within the article (and any supplementary information files).

Given two arbitrary constants α and β, the following identity holds

$$\begin{equation}\left(I+\alpha E\right)\left(I+\beta E\right)=I+\left(\alpha +\beta \right)E+\alpha \beta {E}^{2}\end{equation} \tag{ A1 }$$

for any matrix E. If E is the matrix defined in (8) one gets:

$$\begin{equation}{E}^{2}=NE\end{equation} \tag{ A2 }$$

and then

$$\begin{equation}\left(I+\alpha E\right)\left(I+\beta E\right)=I+\left(\alpha +\beta +\alpha \beta N\right)E\end{equation} \tag{ A3 }$$

If $\beta =-\frac{\alpha }{1+\alpha N}$, from (A3) we get

$$\begin{equation}{\left(I+\alpha E\right)}^{-1}=\left(I+\beta E\right)\end{equation} \tag{ A4 }$$

The above result can be used to derive the structure of matrix S, as reported in the main body of the paper. To this end, we begin by rewriting the matrices A and B as

$$\begin{equation}A=\left[\left(1-\mathrm{i}\right)E-\mathrm{i}I\right]U=-\mathrm{i}\left[I+\left(\mathrm{i}+1\right)E\right]U\end{equation} \tag{ A5 }$$

$$\begin{equation}B=\left[\left(\mathrm{i}+1\right)E+\mathrm{i}I\right]U=\mathrm{i}\left[I+\left(1-\mathrm{i}\right)E\right]U\end{equation} \tag{ A6 }$$

From (A4) we obtain

$$\begin{equation}{B}^{-1}=-\mathrm{i}{U}^{-1}{\left[I+\left(\mathrm{i}+1\right)E\right]}^{-1}=-\mathrm{i}{U}^{-1}\left[I+\frac{\mathrm{i}-1}{1+\left(1-\mathrm{i}\right)N}E\right]\end{equation} \tag{ A7 }$$

and then, after some calculations, we get

$$\begin{equation*}B+A{B}^{-1}A=\mathrm{i}\left[I+\left(1-\mathrm{i}\right)E+I+\frac{3\mathrm{i}+2\mathrm{i}N+1}{1+\left(1-\mathrm{i}\right)N}E\right]U=2\mathrm{i}\left[I+\frac{1+\mathrm{i}}{1+\left(1-\mathrm{i}\right)N}E\right]U\end{equation*}$$

In conclusion,

$$\begin{equation*}S={\left(B+A{B}^{-1}A\right)}^{-1}=-\frac{\mathrm{i}}{2}{U}^{-1}\left[I-\frac{1+\mathrm{i}}{1+2N}E\right]\end{equation*}$$

because of

$$\begin{equation}{\left[I+\frac{1+\mathrm{i}}{1+\left(1-\mathrm{i}\right)N}E\right]}^{-1}=\left[I-\frac{1+\mathrm{i}}{1+2N}E\right]\end{equation} \tag{ A8 }$$

which proves the results.

The aim of this section is to detail the computations needed to obtain explicitly the entries of the Laplace matrix L under the assumption U = qI starting thus from equations (16) and (18).

Let us begin by computing the diagonal elements of L, namely L_ii for i = 1, ..., N. For i = 1, one gets:

$$\begin{equation}{L}_{11}=\sum _{k=2}^{N+1}{V}_{1k}{D}_{kk}{W}_{k1}\end{equation} \tag{ B1 }$$

From (4) and (6) we obtain for k = 2, ..., N + 1:

$$\begin{equation}{V}_{1k}=-1-\mathrm{i}\end{equation} \tag{ B2 }$$

$$\begin{equation}{W}_{k1}=\frac{-1+\mathrm{i}}{2\left(2N+1\right)}\end{equation} \tag{ B3 }$$

Then

$$\begin{equation}{L}_{11}=2\sum _{k=2}^{N+1}\mathfrak{R}\left[\left(-1-\mathrm{i}\right){D}_{kk}\frac{-1+\mathrm{i}}{2\left(2N+1\right)}\right]\end{equation} \tag{ B4 }$$

$$\begin{equation}=2\sum _{k=2}^{N+1}\mathfrak{R}\left[{D}_{kk}\left(-1-\mathrm{i}\right)\frac{-1+\mathrm{i}}{2\left(2N+1\right)}\right]\end{equation} \tag{ B5 }$$

$$\begin{equation}=\frac{2}{2N+1}\sum _{k=2}^{N+1}{\alpha }_{k}\end{equation} \tag{ B6 }$$

where use has been made of the identity $\mathfrak{R}\left({D}_{kk}\right)={\alpha }_{k}$.

We can proceed in analogy for the others diagonal elements of the Laplacian matrix. In particular, for s = 2, ..., N + 1, we have that

$$\begin{equation}{W}_{ss}={S}_{ss}=\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)},\end{equation} \tag{ B7 }$$

and then

$$\begin{equation}{L}_{ss}=2\sum _{k=2}^{N+1}\mathfrak{R}\left[{D}_{kk}{V}_{sk}{W}_{ks}\right]\end{equation} \tag{ B8 }$$

$$\begin{equation}=2\mathfrak{R}\left({D}_{ss}{V}_{ss}{W}_{ss}\right)\end{equation} \tag{ B9 }$$

$$\begin{equation}=2\mathfrak{R}\left[\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B10 }$$

$$\begin{equation}=\frac{2N{\alpha }_{s}+{\beta }_{s}}{2N+1}\end{equation} \tag{ B11 }$$

Conversely, for s = N + 2, ..., 2N + 1, we have

$$\begin{equation}{W}_{ss}=\mathrm{i}{S}_{ss}=\mathrm{i}\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)},\end{equation} \tag{ B12 }$$

and consequently:

$$\begin{equation}{L}_{ss}=2\mathfrak{R}\left[\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B13 }$$

$$\begin{equation}=\frac{2N{\alpha }_{s}-{\beta }_{s}}{2N+1}\end{equation} \tag{ B14 }$$

Let us proceed now with the other elements of the first row of L. For t = 2, ..., N + 1 we can write

$$\begin{equation}{L}_{1t}=2\sum _{k=2}^{N+1}\mathfrak{R}\left({V}_{1k}{D}_{kk}{W}_{kt}\right)\end{equation} \tag{ B15 }$$

$$\begin{equation}=2\left[\mathfrak{R}\left({V}_{1t}{D}_{tt}{W}_{tt}\right)+\sum _{k=2,k\enspace =t}^{N+1}\mathfrak{R}\left({V}_{1k}{D}_{kk}{W}_{kt}\right)\right]\end{equation} \tag{ B16 }$$

$$\begin{equation}=2\left\{\mathfrak{R}\left[\left(-1-\mathrm{i}\right)\left({\alpha }_{t}+\mathrm{i}{\beta }_{t}\right)\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\right]\right.\end{equation} \tag{ B17 }$$

$$\begin{equation}\quad +\left.\sum _{k=2,k\ne t}^{N+1}\mathfrak{R}\left[\left(-1-\mathrm{i}\right)\left({\alpha }_{k}+\mathrm{i}{\beta }_{k}\right)\frac{\mathrm{i}-1}{2\left(2N+1\right)}\right]\right\}\end{equation} \tag{ B18 }$$

$$\begin{equation}=\frac{1-2N}{1+2N}{\alpha }_{t}-{\beta }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\end{equation} \tag{ B19 }$$

and for t = N + 2, ..., 2N + 1 we get

$$\begin{equation}{L}_{1t}=2\left\{\mathfrak{R}\left[\left(-1+\mathrm{i}\right)\left({\alpha }_{t}-\mathrm{i}{\beta }_{t}\right)i\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\right]\right.\end{equation} \tag{ B20 }$$

$$\begin{equation}\quad +\left.\sum _{k=N+2,k\ne t}^{2N+1}\mathfrak{R}\left[\left(-1+\mathrm{i}\right)\left({\alpha }_{k}-\mathrm{i}{\beta }_{k}\right)\mathrm{i}\frac{\mathrm{i}-1}{2\left(2N+1\right)}\right]\right\}\end{equation} \tag{ B21 }$$

$$\begin{equation}=\frac{1-2N}{1+2N}{\alpha }_{t}+{\beta }_{t}+\frac{2}{2N+1}\sum _{k=N+2,k\ne t}^{2N+1}{\alpha }_{k}\end{equation} \tag{ B22 }$$

$$\begin{equation}=\frac{1-2N}{1+2N}{\alpha }_{t}+{\beta }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t-N}^{N+1}{\alpha }_{k}\end{equation} \tag{ B23 }$$

For s = 2, ..., N + 1 and t = N + 2, ..., 2N + 1 with t ≠ s − N, we get:

$$\begin{equation}{L}_{st}=\sum _{k=2}^{2N+1}{V}_{sk}{D}_{kk}{W}_{kt}\end{equation} \tag{ B24 }$$

$$\begin{equation}={V}_{ss}{D}_{ss}{W}_{st}+{V}_{s,s+N}{D}_{s+N,s+N}{W}_{s+N,t}\end{equation} \tag{ B25 }$$

$$\begin{equation}=\mathrm{i}\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\left(-\mathrm{i}\right)\frac{-\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B26 }$$

$$\begin{equation}\quad +\left(-\mathrm{i}\right)\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B27 }$$

$$\begin{equation}=\frac{-{\alpha }_{s}+{\beta }_{s}}{2N+1}\end{equation} \tag{ B28 }$$

while, for s = N + 2, ..., 2N + 1 and t = 2, ..., N + 2 with t ≠ s − N, the following expression holds:

$$\begin{equation}{L}_{st}={V}_{ss}{D}_{ss}{W}_{st}+{V}_{s,s-N}{D}_{s-N,s-N}{W}_{s-N,t}\end{equation} \tag{ B29 }$$

$$\begin{equation}=\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\frac{-\mathrm{i}-1}{2\left(2N+1\right)}+\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\frac{\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B30 }$$

$$\begin{equation}=-\frac{{\alpha }_{s}+{\beta }_{s}}{2N+1}\end{equation} \tag{ B31 }$$

For s = 2, ..., N + 1

$$\begin{equation}{L}_{s,s+N}=\end{equation} \tag{ B32 }$$

$$\begin{equation}{V}_{ss}{D}_{ss}{W}_{s,s+N}+{V}_{s,s+N}{D}_{s+N,s+N}{W}_{s+N,s+N}\end{equation} \tag{ B33 }$$

$$\begin{equation}=\mathrm{i}\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\left(-\mathrm{i}\right)\frac{2N\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B34 }$$

$$\begin{equation}\quad +\left(-\mathrm{i}\right)\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B35 }$$

$$\begin{equation}=-\frac{{\alpha }_{s}+2N{\beta }_{s}}{2N+1}\end{equation} \tag{ B36 }$$

while, for s = N + 2, ..., 2N + 1, we obtain:

$$\begin{equation}{L}_{s,s-N}=\end{equation} \tag{ B37 }$$

$$\begin{equation}{V}_{ss}{D}_{ss}{W}_{s,s+N}+{V}_{s,s-N}{D}_{s-N,s-N}{W}_{s-N,s-N}\end{equation} \tag{ B38 }$$

$$\begin{equation}=\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\frac{2N\mathrm{i}-1}{2\left(2N+1\right)}+\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{-2N\mathrm{i}-1}{2\left(2N+1\right)}\end{equation} \tag{ B39 }$$

$$\begin{equation}=-\frac{{\alpha }_{s}-2N{\beta }_{s}}{2N+1}\end{equation} \tag{ B40 }$$

For t = 2, ..., N + 1

$$\begin{equation}{L}_{t1}=2\mathfrak{R}\left({V}_{tt}{D}_{tt}{W}_{t1}\right)\end{equation} \tag{ B41 }$$

$$\begin{equation}=2\mathfrak{R}\left[\mathrm{i}\left({\alpha }_{t}+\mathrm{i}{\beta }_{t}\right)\frac{\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B42 }$$

$$\begin{equation}=-\frac{{\alpha }_{t}-{\beta }_{t}}{2N+1}\end{equation} \tag{ B43 }$$

and for t = N + 2, ..., 2N + 1

$$\begin{equation}{L}_{t1}=2\mathfrak{R}\left({V}_{tt}{D}_{tt}{W}_{t1}\right)\end{equation} \tag{ B44 }$$

$$\begin{equation}=2\mathfrak{R}\left[\left({\alpha }_{t}-\mathrm{i}{\beta }_{t}\right)\frac{-\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B45 }$$

$$\begin{equation}=\frac{-{\alpha }_{t}-{\beta }_{t}}{2N+1}\end{equation} \tag{ B46 }$$

For t, s = 2, ..., N + 1 and s ≠ t

$$\begin{equation}{L}_{st}=2\mathfrak{R}\left({V}_{ss}{D}_{ss}{W}_{st}\right)\end{equation} \tag{ B47 }$$

$$\begin{equation}=2\mathfrak{R}\left[\mathrm{i}\left({\alpha }_{s}+\mathrm{i}{\beta }_{s}\right)\frac{\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B48 }$$

$$\begin{equation}=\frac{-{\alpha }_{s}+{\beta }_{s}}{2N+1}\end{equation} \tag{ B49 }$$

and, finally, for s, t = N + 2, ..., 2N + 1 and s ≠ t, one gets:

$$\begin{equation}{L}_{st}=2\mathfrak{R}\left({V}_{ss}{D}_{ss}{W}_{st}\right)\end{equation} \tag{ B50 }$$

$$\begin{equation}=2\mathfrak{R}\left[\left({\alpha }_{s}-\mathrm{i}{\beta }_{s}\right)\mathrm{i}\frac{\mathrm{i}-1}{2\left(2N+1\right)}\right]\end{equation} \tag{ B51 }$$

$$\begin{equation}=-\frac{{\alpha }_{s}+{\beta }_{s}}{2N+1}\end{equation} \tag{ B52 }$$

Summing up, we have here provided closed-form analytical expressions for all the entries of the Laplacian matrix, as a function of the assigned spectrum.

The aim of this section is to work out the algebraic steps needed to rewrite (46) in the simpler form given by (47) and (48).

Let us thus rewrite (46)

$$\begin{equation*}\begin{cases}2N{\alpha }_{t}{< }{\beta }_{t}{< }-2N{\alpha }_{t}\quad \hfill \\ {\beta }_{t}{ >}-\frac{1-2N}{2N+1}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\beta }_{t}{< }\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\alpha }_{t}{< }{\beta }_{t}{< }-{\alpha }_{t}\quad \hfill \\ \frac{{\alpha }_{t}}{2N}{< }{\beta }_{t}{< }-\frac{{\alpha }_{t}}{2N}\quad \hfill \end{cases}\end{equation*}$$

where the inequalities hold for t = 2, ..., N + 1. Then the second and the third conditions of system (46) can be matched simultaneously provided that:

$$\begin{equation}{\alpha }_{t}{< }\frac{2}{2N-1}\sum _{k\ne t}{\alpha }_{k}\quad \forall \enspace t=2,\dots ,N+1\end{equation} \tag{ C1 }$$

Remark that:

$$\begin{equation}2N{\alpha }_{t}{< }{\alpha }_{t}{< }\frac{{\alpha }_{t}}{2N}\end{equation} \tag{ C2 }$$

holds for arbitrary values of α_t < 0 and N. Hence, system (46) simplifies as follows:

$$\begin{equation}\begin{cases}\frac{{\alpha }_{t}}{2N}{< }{\beta }_{t}{< }-\frac{\alpha }{2N}\quad \hfill \\ {\alpha }_{t}{< }\frac{2}{2N-1}\sum _{k\ne t}{\alpha }_{k}\quad \hfill \\ {\beta }_{t}{ >}-\frac{1-2N}{2N+1}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\beta }_{t}{< }\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \end{cases}\end{equation} \tag{ C3 }$$

for t = 2, ..., N + 1.

Let us now focus on the conditions for β_t . We assume that the following condition holds:

$$\begin{equation}\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}{ >}-\frac{{\alpha }_{t}}{2N}\end{equation} \tag{ C4 }$$

and we set to explore its consequences. Equation (C4) yields:

$$\begin{equation}{\alpha }_{t}{< }\frac{4N}{4{N}^{2}-4N-1}\sum _{k\ne t}{\alpha }_{k}=\frac{1}{N-1-\frac{1}{4N}}\sum _{k\ne t}{\alpha }_{k}\end{equation} \tag{ C5 }$$

For N > 1, conditions (C3) maps therefore in the following equivalent system:

$$\begin{equation}\begin{cases}\frac{{\alpha }_{t}}{2N}{< }{\beta }_{t}{< }-\frac{{\alpha }_{t}}{2N}\quad \hfill \\ {\alpha }_{t}{< }\frac{2}{2N-1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\alpha }_{t}{< }\frac{1}{N-1-\frac{1}{4N}}\sum _{k\ne t}{\alpha }_{k}\quad \hfill \end{cases}\end{equation} \tag{ C6 }$$

Remark that:

$$\begin{equation}\frac{1}{N-1-\frac{1}{4N}}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}{< }\frac{2}{2N-1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\end{equation} \tag{ C7 }$$

due to the inequality

$$\begin{equation}\frac{1}{N-1-\frac{1}{4N}}{ >}\frac{2}{2N-1}\end{equation} \tag{ C8 }$$

which holds for each N > 1. Hence, system (C6) takes the form:

$$\begin{equation}\begin{cases}\frac{{\alpha }_{t}}{2N}{< }{\beta }_{t}{< }-\frac{{\alpha }_{t}}{2N}\quad \hfill \\ {\alpha }_{t}{< }\frac{1}{N-1-\frac{1}{4N}}\sum _{k\ne t}{\alpha }_{k}\quad \hfill \end{cases}\end{equation} \tag{ C9 }$$

Then, (C9) has no solutions, under the working hypothesis that we have put forward to deriving it. In fact:

$$\begin{equation}{\alpha }_{t}+\frac{1}{N-1-\frac{1}{4N}}{\alpha }_{t}{< }\frac{1}{N-1-\frac{1}{4N}}\sum _{k}{\alpha }_{k}\end{equation} \tag{ C10 }$$

that is

$$\begin{equation}{\alpha }_{t}{< }\frac{1}{N-\frac{1}{4N}}\sum _{k}{\alpha }_{k}=\frac{4N}{4{N}^{2}-1}\sum _{k}{\alpha }_{k}\end{equation} \tag{ C11 }$$

and summing on every t we get

$$\begin{equation}\sum _{t}{\alpha }_{t}{< }\frac{4{N}^{2}}{4{N}^{2}-1}\sum _{k}{\alpha }_{k}\end{equation} \tag{ C12 }$$

which, in turn, implies $\frac{4{N}^{2}}{4{N}^{2}-1}{< }1$, a condition that is obviously never met. We now go back to revise ansatz (C4), and consider the alternative scenario:

$$\begin{equation}\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}{< }-\frac{{\alpha }_{t}}{2N}\end{equation} \tag{ C13 }$$

Then, (C3) becomes

$$\begin{equation}\begin{cases}_{t}{ >}-\frac{1-2N}{2N+1}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\beta }_{t}{< }\frac{1-2N}{2N+1}{\alpha }_{t}+\frac{2}{2N+1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\alpha }_{t}{< }\frac{2}{2N-1}\sum _{k=2,k\ne t}^{N+1}{\alpha }_{k}\quad \hfill \\ {\alpha }_{t}{ >}\frac{1}{N-1-\frac{1}{4N}}\sum _{k\ne t}{\alpha }_{k}\quad \hfill \end{cases}\end{equation} \tag{ C14 }$$

Following a path analogous to the one discussed above, we get:

$$\begin{equation}\frac{4N}{4{N}^{2}-1}\sum _{k}{\alpha }_{k}{< }{\alpha }_{t}{< }\frac{2}{2N+1}\sum _{k}{\alpha }_{k}\end{equation} \tag{ C15 }$$

and

$$\begin{equation}{\alpha }_{t}-\frac{2}{2N+1}\sum _{k}{\alpha }_{k}{< }{\beta }_{t}{< }-{\alpha }_{t}+\frac{2}{2N+1}\sum _{k}{\alpha }_{k}\end{equation} \tag{ C16 }$$