Universal hypotrochoidic law for random matrices with cyclic correlations

TL;DR

This paper generalizes the elliptic law for random matrices to include higher-order cyclic correlations, revealing that their eigenvalue spectra are bounded by hypotrochoid curves with universal applicability across different matrix types.

Contribution

It introduces a new hypotrochoid law for eigenvalue distributions in matrices with cyclic correlations, extending the classical elliptic law to higher-order motifs.

Findings

Eigenvalue spectra are bounded by hypotrochoid curves with k-fold symmetry.

The law applies to both full and sparse matrices, demonstrating universality.

Extension to matrices with competing cycle motifs leads to polytrochoid spectral boundaries.

Abstract

The celebrated elliptic law describes the distribution of eigenvalues of random matrices with correlations between off-diagonal pairs of elements, having applications to a wide range of physical and biological systems. Here, we investigate the generalization of this law to random matrices exhibiting higher-order cyclic correlations between $k$-tuples of matrix entries. We show that the eigenvalue spectrum in this ensemble is bounded by a hypotrochoid curve with $k$-fold rotational symmetry. This hypotrochoid law applies to full matrices as well as sparse ones, and thereby holds with remarkable universality. We further extend our analysis to matrices and graphs with competing cycle motifs, which are described more generally by polytrochoid spectral boundaries.

Equations

$$G(z,z^{\ast })={\left(\begin{array}{cc}z\text{openone}-M& i\lambda \text{openone}\\ i\lambda \text{openone}& z^{\ast }\text{openone}-M^{†}\end{array}\right)}^{-1}\equiv \left(\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right),$$

$$\mu (z)$$

$$g(z,z^{\ast })$$

$$G={\mathcal{Z}}^{-1}\sum _{\ell =0}^{\infty }(\mathcal{M}{\mathcal{Z}}^{-1}{)}^{\ell }$$

$$\mathcal{Z}$$

$$\overline{G}$$

$$={\mathcal{Z}}^{-1}+{\mathcal{Z}}^{-1}\dot{\mathcal{M}}\overline{G}\dot{\mathcal{M}}\overline{G},$$

$$Z=N/\overline{g}+\left(\begin{array}[]{cc}\tau_{2}^{2}\overline{g}_{11}&\sigma^{2}\overline{g}_{12}\\ \sigma^{2}\overline{g}_{21}&\tau_{2}^{2}\overline{g}_{22}\\ \end{array}\right).$$

$$z={\sigma }^2{\overline{g}}_{11}^{\ast }+{\tau }_2^2{\overline{g}}_{11}.$$

$$\overline{G}$$

$$Z=N/\overline{g}+\left(\begin{array}[]{cc}\tau_{k}^{k}\overline{g}_{11}^{k-1}&\sigma^{2}\overline{g}_{12}\\ \sigma^{2}\overline{g}_{21}&\tau_{k}^{k}\overline{g}_{22}^{k-1}\\ \end{array}\right)$$

$$\begin{array}{rl}z& ={\sigma }^2{\overline{g}}_{11}^{\ast }+{\tau }_k^k{\overline{g}}_{11}^{k-1},\\ {\overline{g}}_{12}^2& =|{\overline{g}}_{11}{|}^2-N/{\sigma }^2.\end{array}$$

$$z_b(\phi )=e^{-i\phi }+{\rho }_k{e}^{i(k-1)\phi }.$$

$$z_b(\phi )=\frac{1}{t}{e}^{-i\phi }+\hat{d}{t}^{k-1}{e}^{i(k-1)\phi },$$

$$z_b(\phi )={\hat{d}}^{-1/2}{t}^{-1}{e}^{-i\phi }-({\hat{d}}^{1/2}-{\hat{d}}^{-1/2})t^2{e}^{2i\phi }$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\overline{\text{Tr}{M}^{3l}}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\overline{\text{Tr}(M{M}^{†}{)}^l}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\overline{\text{Tr}{M}^{3l}}$$

$$\underset{N\to \infty }{\lim }\frac{1}{N}\overline{\text{Tr}(M{M}^{†}{)}^l}$$

$$A^{(m)}(l,\hat{d})=\frac{d}{l}\sum _{j=0}^{l-1}\left(\genfrac{}{}{0ex}{}{ml}{j}\right)(l-j)(\hat{d}-1{)}^j.$$

$$z_b(\phi )=e^{-i\phi }+\sum _k{\rho }_k{e}^{i(k-1)\phi }.$$

$$\frac{z}{\bar{d}}=e^{-i\phi }+d_1{\left(\frac{w_1}{\bar{d}}\right)}^{k_1}{e}^{i(k_1-1)\phi }+d_2{\left(\frac{w_2}{\bar{d}}\right)}^{k_2}{e}^{i(k_2-1)\phi },$$

$$G_{n\tilde{n},n\tilde{n}}=(Z-{\mathcal{M}}_{n\tilde{n},\star }{G}^{(n)}{\mathcal{M}}_{\star ,n\tilde{n}}{)}^{-1},$$

$$G_{n\tilde{n},n\tilde{n}}=\bigg{(}Z-\sum_{m\sim n}G^{(n,m)}_{k\tilde{2},2\tilde{k}}\bigg{)}^{-1},$$

$$G^{(n)}_{s\tilde{s},s\tilde{s}}=\Bigg{[}Z\otimes\openone_{k-1}-\mathcal{S}-(d-1)G^{(n)}_{k\tilde{2},2\tilde{k}}\otimes\openone_{k-1}\Bigg{]}^{-1}\,,$$

$$G^{(n)}_{s\tilde{s},s\tilde{s}}=\left(\begin{array}[]{cc}\alpha-S&\beta\\ \gamma&\alpha^{*}-S^{T}\\ \end{array}\right)^{-1}$$

$$\alpha$$

$$\beta$$

$$G^{(n)}_{s\tilde{s},s\tilde{s}}=\left(\begin{array}[]{cc}A&-\beta AA^{\dagger}\\ -\gamma A^{\dagger}A&A^{\dagger}\\ \end{array}\right)+\mathcal{O}(\beta^{2},\beta\gamma,\gamma^{2}),$$

$$G_{k,\tilde{k}}^{(n)}=(d-1)G_{k,\tilde{k}}^{(n)}\sum _{\ell =1}^{k-1}|\alpha {|}^{-2\ell },$$

$$z_b(\phi )=\frac{1}{t}{e}^{-i\phi }+(d-1)t^{k-1}{e}^{i(k-1)\phi }.$$

$$G_{n\tilde{n},n\tilde{n}}=(Z-d_1{w}_1^2{G}_{k_1\tilde{2},2{\tilde{k}}_1}^{(1)}-d_2{w}_2^2{G}_{k_2\tilde{2},2{\tilde{k}}_2}^{(2)}{)}^{-1}$$

$$\displaystyle G^{(1)}_{s_{1}\tilde{s}_{1},s_{1}\tilde{s}_{1}}=\Bigg{[}Z\otimes\openone_{k_{1}-1}-\mathcal{S}^{(1)}-(d_{1}-1)w_{1}^{2}G^{(1)}_{k_{1}\tilde{2},2\tilde{k}_{1}}\otimes\openone_{k_{1}-1}-d_{2}w_{2}^{2}G^{(2)}_{k_{2}\tilde{2},2\tilde{k}_{2}}\otimes\openone_{k_{1}-1}\Bigg{]}^{-1},$$

Full text

As the applications of random matrix theory have diversified, so too have the ensembles under study. Existing generalizations of the elliptic law fall broadly into three categories. There is a large body of theoretical work concerning random matrix ensembles defined by potential functions, where the quadratic case recovers the elliptic law, but more complicated potentials produce interesting spectral distributions (see e.g. Elbau and Felder (2005); Bleher and Kuijlaars (2012)). Alternatively, for some applications it is necessary to impose a system-level structure such as modularity (see e.g. Grilli et al. (2016)), often by addition or multiplication with another matrix; see Rogers (2010); Bordenave (2011) for methods and rigorous mathematical results. Lastly, growing interest in the study of complex networks has led many authors to consider sparse ensembles containing many zero matrix elements, where the location of the non-zero elements encodes the adjacency matrix of some directed graph (or digraph), and can have striking eigenvalue distributions Metz et al. (2011); Bollé et al. (2013) (see Metz et al. (2018) for a recent topical review). In none of these directions of work has the problem of high-order correlations been fully and directly addressed. This is not for lack of interest. Multi-party interactions in dense systems are important in biological applications such as ecology Levine et al. (2017), stabilization of microbial communities Guo and Boedicker (2016) or gene–gene interactions Wang et al. (2014) and can provide valuable engineering insights into machine learning Sejnowski (1986); Personnaz et al. (1987) and control theory Doyle and Stein (1981). Moreover, most of the sparse random matrix literature hinges on the assumption of a tree-like interaction structure; it remains a long-standing and important problem to allow the relaxation of this assumption, inducing higher-order correlations.

Here, we investigate the generalization of the elliptic law to ensembles, both dense and sparse, featuring higher-order correlations. In the case that there is a single dominant correlation order $k$, we find that the spectrum is bounded by a hypotrochoid curve (that is, the path of a point located on a small wheel that rolls inside a larger wheel) with $k$-fold rotational symmetry; see Fig. 1 for an illustration. Surprisingly, we thereby recover the spectral boundary of a highly structured system, a sparse regular digraph Metz et al. (2011); Bollé et al. (2013), but in a much more general and universal context. Extending this result further to the case where more than one correlation order is present, we find a more general family of polytrochoid boundary curves described by multiple linked gears 111See, for example, US patent 3037488A “Rotary hydraulic motor”, George M Barrett, (1962).

Dense matrices.—We study $N\times N$-dimensional nonhermitian matrices with real or complex entries independently drawn from a distribution with zero mean and bounded variance. For any of such matrix $M$, the density $\mu(z)$ of complex eigenvalues $z$ can be obtained from a Green’s function of the form Feinberg and Zee (1997); Janik et al. (1999); Rogers (2010)

[TABLE]

where $z^{*}$ denotes complex conjugation and $M^{\dagger}$ the adjoint; for a real matrix this is just the transpose. Thereby

[TABLE]

For a random matrix $M$, the ensemble-averaged density $\overline{\mu}(z)$ therefore follows from $\overline{g}(z,z^{*})$.

This can be used to derive the elliptic law of random matrices in the large $N$ limit, which serves as useful preparation. Let us expand

[TABLE]

into a geometric series, where

[TABLE]

When we perform the averaging, only certain products of matrix elements $M_{nm}$ survive. These can be organised into groups of terms whose number scales differently in the matrix dimension $N$. In particular, there will be many terms where we can pair $M_{nm}$ with $(M^{\dagger})_{mn}$. If no further correlations are present the leading order is

[TABLE]

where the dot denotes matrices in which we pair the elements $M_{nm}$ with $(M^{\dagger})_{mn}$. This factorization of the average is called the non-crossing or planar approximation.

The elliptic law is derived straightforwardly from the expression above. Carrying out the average of the dotted matrices, taking the partial trace on both sides, and rearranging for $Z$, one obtains

[TABLE]

Comparing terms in the off-diagonal in this equation (and noting the constraints $\bar{g}_{11}=\bar{g}^{*}_{22}$ and $g_{12}=g_{21}>0$), we find two possible solutions: (i) either $g_{12}=0$, or (ii) $|\bar{g}_{11}|^{2}-\bar{g}_{12}^{2}=N/\sigma^{2}$. The first of these yields $\bar{g}_{11}$ proportional to $z$ and therefore holds only outside of the support of the spectrum. Examining the diagonal elements of (6) in the case (ii) we obtain

[TABLE]

which must be solved together with the constraint that $|\bar{g}_{11}|^{2}-N/\sigma^{2}>0$. The boundary of the spectrum is therefore found by determining the values of $z$ for which $|\bar{g}_{11}|^{2}=N/\sigma^{2}$; in the present case one finds an ellipse with foci $\pm 2\tau\sqrt{N}$. Inside the support, one can solve (7) to determine $\overline{g}_{11}$ and apply (2) to find that the spectral density $\mu(z)$ is uniform.

To generalize these results to ensembles with high-order correlations where $\overline{\textrm{Tr}M^{k}}/N$ is a fixed parameter, we reinterpret the Hermitian contributions of weight $\tau_{2}$ in the elliptic law as correlations of order $k=2$. Introducing correlations of general order $k$, we pick up additional contributions corresponding to contractions

[TABLE]

We assume for now that there is just a single extra term of these, of fixed $k$ and with weight $\tau_{k}$. This gives

[TABLE]

and results, inside the spectrum, in the equations

[TABLE]

The boundary of the support is therefore determined by the condition $|\overline{g}_{11}|^{2}=N/\sigma^{2}$. In the large $N$ limit we choose the scalings $\sigma^{2}=N^{-1}$, $\tau_{k}^{k}=\rho_{k}N^{1-k}$ to balance the contribution of terms in (10). With the parameterization $\overline{g}_{11}=Ne^{i\varphi}$, which we insert into Eq. (10), the boundary curve then becomes

[TABLE]

This equation is precisely the complex parameterization of a hypotrochoid curve in which the small and large wheels have radii in a ratio of $1:(k-1)$. In Fig. 2 we illustrate this result numerically for various values of $k$ and $\rho_{k}$; see Appendix A for the matrix generation algorithm.

Note that (although hard to determine from Fig. 2), for $k>2$ the density of eigenvalues inside the hypotrochoid support is in fact not uniform; in general solutions to (10) do not have the property that $\partial\bar{g}_{11}/\partial z^{*}$ is constant. This distribution is, however, universal in the sense that it is determined entirely by the parameters $\sigma^{2}$ and $\rho_{k}$, and other properties of the distribution of matrix elements are unimportant. We will now explore to what extent this universality extends to sparse matrices.

Sparse digraphs.— Square matrices can be seen as an alternative representation of weighted digraphs, where the entry $M_{nm}$ corresponds to the weight of the edge going from node $n$ to node $m$. This implies that for large dense graphs with random weights we can obtain the eigenvalues of their adjacency matrix by the previous methods. However, it is well-known that sparsity can change the eigenvalue distribution substantially Semerjian and Cugliandolo (2002); Biroli and Monasson (1999); Rogers and Castillo (2009). Surprisingly, the hypotrochoidic law (11) also applies to highly-structured sparse systems. Randomly generated digraphs in which each node belongs to exactly $d$ directed cycles of length $k$ were studied in Metz et al. (2011); Bollé et al. (2013), where it was shown that the adjacency matrices of these graphs hypotrochoidic spectra in the limit of large network size.

This result extends further to disordered directed random graph ensembles. In Appendix B we apply effective medium approximation (EMA) Dorogovtsev et al. (2003) to derive the following hypotrochoidic law for the spectral boundary of cyclic random digraphs:

[TABLE]

where $k$ is the length of cycles, $\hat{d}$ (degree biased) number of cycles per node, and $t$ is the unique positive real solution of $\hat{d}t^{2k}-dt^{2}+1=0$. The EMA is technically valid in the case $1\ll d\ll N$, however, numerical simulations show excellent agreement down to relatively small values of $d$, see Fig. 3. For fixed node degrees it is exact.

To understand how both sparse and dense cyclic ensembles share the same universal behaviour, we explore the asymptotic behaviour of (29) as $\hat{d}\to\infty$. For a direct comparison we rescale the adjacency matrix of the graph by $\hat{d}^{-1/2}$, which corresponds to the factor $\sigma=N^{-1/2}$ that we used before and that normalizes rows and columns. The previous result (29) for $k=3$ then reads

[TABLE]

with $t$ being the same as in Eq. (29). We note that for large $\hat{d}$, $t\sim\hat{d}^{-1/2}$, so that we recover the support (11) for full matrices with an effective parameter $\rho_{3}\sim\hat{d}^{-1/2}$.

The same connection indeed appears when we compare the definition of the quantities $\sigma_{2}$ and $\rho_{\tau}$ from the traces of $M$. By applying the noncrossing approximation to the full matrix we find

[TABLE]

while carrying out the corresponding combinatorics for the graphs gives

[TABLE]

which we express in terms of the number of number of self-returning walks of length $2l$ from the root of an infinite $k$-regular tree:

[TABLE]

Asymptotically for large $\hat{d}$, $A^{(m)}(l,\hat{d})=\frac{\hat{d}^{l}}{l}\binom{ml}{l-1}=\frac{\hat{d}^{l}}{ml-l+1}\binom{ml}{l}$, so that both expressions again match up for $\rho_{3}\sim\hat{d}^{-1/2}$. It is worth mentioning that the traces of $M$, which can be computed explicitly, relate to asymptotic spectral statistics in the so-called ”trace formulas” which are relevant in random matrix theory Haake et al. (1996) as well as in semiclassical physics Haake (2013); Gutzwiller (2013).

Polytrochoid spectra.— Starting from Eq. (8), our approach generalizes to matrices with correlations of multiple orders, leading to a boundary curve

[TABLE]

The curve described by this equation is an example of the very general polytrochoid family. For the particular case of two competing correlation orders, the curve is described by the tracing the path of a point in a wheel rotating around a larger wheel, which itself is rotating in the opposite direction around the origin. Adding further correlation orders would correspond to the addition of more linked wheels with the same drawing procedure.

Polytrochoid spectral boundaries are also found in digraphs with mixed cycle lengths, where each vertex connects to $d_{k}$ $k$-cycles of weight $w_{k}$; Appendix D for a derivation in the case of two competing cycle motifs. In the limit of large degrees, we obtain the explicit formula

[TABLE]

where $\bar{d}=\sqrt{d_{1}w_{1}^{2}+d_{2}w_{2}^{2}}$. Fig. 4 shows a numerical illustration of this result, which is in excellent agreement with our derivations for relatively degrees.

Discussion.— We have shown that high-order correlations in the entries of a random matrix give a surprising and even beautiful shape to its eigenvalues, with the boundary given by a hypotrochoid, generalizing the classical result from Girko Girko (1986). Furthermore, we have uncovered a remarkable degree of universality of this result by connecting it to the known case of regular graphs with cyclic motifsMetz et al. (2011), and studied the relationship between both cases. Our derivations are in excellent agreement with numerical results.

Our results have a simple interpretation in terms of systems theory. A feedback loop is a classical control theory tool to enhance or dampen certain frequencies. A cycle of length $k$ where the product of the edge weights is $w$ is a feedback loop with delay $k$ and weight $w$, therefore a graph with abundance of cycles of length $k$ with positive feedback would resonate at a frequency $\frac{1}{k}$. On the spectral side, the presence of large positive correlations of order $k$ leads to dominant eigenvalues (or poles, in control theory terms) with phases $\frac{2\pi j}{k}$ for $j<k$.

This interpretation can also be used in the other direction. Matrices and graphs are useful models to represent the interactions between elements, and our results show that random matrix theory can account for intricate multi-element interactions. That is, the stability and resonance properties of complex systems with interconnected feedback loops can be studied by converting the loops into graph cycles and then applying our results in random matrix theory. For example, designing large networked systems such as the Internet or power grids is a challenging problem due to the amount of feedback loops present Fairley (2004); Low et al. (2002). Likewise, biological regulatory systems Thomas et al. (1995); Becskei and Serrano (2000); Csete and Doyle (2004) include many intertwined feedback loops that render their analysis difficult. In all those examples, their stability and dynamical properties can be studied through the spectrum of their adjacency matrix, which can be difficult to estimate. However, finding frequent cycles and their corresponding delays is typically easier Milo et al. (2002); Shen-Orr et al. (2002), meaning that we can leverage the simplicity of finding cycles to obtain rigorous results on the stability of those systems. We hope that the techniques and results presented here may inspire new developments in these fields.

Finally, the system theory interpretation also provides relevant theoretical insights and questions. For instance, the effects of negative correlations, which in systems theory corresponds to a rotation of the poles, is also explained by Eq. (11) and (29). On the other hand, the combination of cycles with positive and negative feedbacks for the same length is not fully covered by our method: in the case of a full matrix or a dense digraph, the positive and negative feedbacks appear in a single connected graph and thus cancel each other – meaning that Eq. (11) and (29) remain valid –, but in very sparse digraphs the presence of many finite-size isolated subgraphs implies that there are components dominated by either positive or negative feedback, and thus the effective medium approximation does not hold.

Appendix A Appendix A: Generation of Matrices with high-order correlations

The algorithm to generate a random matrix with positive correlations of order $k$ proceeds by going through the nodes $n$ one by one and modifying only the weights of the edges between node $n$ and nodes $<n$. In more precise terms:

Generate a $N\times N$ matrix $M$ with random i.i.d entries. Set $n=k-1$ 2. 2.

Take the square submatrix $M_{n\times n}$ which contains the first $n$ entries of the first $n$ rows. 3. 3.

Then compute the matrix $P_{n,k-1}$ of weights of paths of length $k-1$ between all pairs of nodes that do not pass twice by the same node. 4. 4.

Compute the weights $w(e)$ of the cycles passing through node $n+1$ through all edges $e=(c,n+1)$ for $c\leq n$ by $w(e)=P_{n,k}M_{c,l+1}M_{n+1,c}$. 5. 5.

For each one of those edges, if the weight $w(e)$ is negative, flip the sign of $e$ with probability $p$. 6. 6.

Increase $n$ by 1. If $n<N$, go back to step 2.

The paths can be computed by $P_{n,k}=M_{n\times n}P_{n,k-1}-\text{diag}\left[M_{n\times n}P_{n,k-1}\right]$ where $\text{diag}\left[\cdot\right]$ is the operator that sets all non-diagonal entries to zero and the recursion is started by $P_{n,1}=M_{n\times n}$. Combining correlations of different orders can be achieved by adding two matrices generated by this algorithm and normalizing. The computational complexity of this algorithm is limited by the matrix multiplications needed to calculate the path weights, which has a complexity of $O(n^{2.373})$ and has to be called $n$ times at each of the $N$ iterations. This leaves us with $O(N^{4.373})$.

Appendix B Appendix B: Cavity method calculation for spectral density of regular cyclic graphs

We derive the hypotrochoidic law for digraphs making use of an adaption of the cavity method Rogers and Castillo (2009); Rogers et al. (2010), as used in Metz et al. (2011); Bollé et al. (2013), coupled with an effective medium approximation (EMA) Dorogovtsev et al. (2003). We start by expanding the Green’s function $G=(\mathcal{Z}-\mathcal{M})^{-1}$ around an arbitrary node $n$. Note that the matrix $\mathcal{M}$ describes the original graph $M$ and a replica with inverted connections $M^{\dagger}$. Consider then the $2\times 2$ block of $G$ corresponding a node $n$ and the corresponding node $\tilde{n}=n+N$ in the replica graph. Using the Schur complement we can write

[TABLE]

where $\mathcal{M}_{n\tilde{n},\star}$ refers to rows $n$ and $\tilde{n}$ of $\mathcal{M}$, $\mathcal{M}_{\star,n\tilde{n}}$ to the corresponding columns, and $G^{(n)}$ is the Green’s function of the graph where nodes $n$ and $\tilde{n}$ are removed.

For graphs composed of $k$-cycles, where $k\geq 3$, we develop our theory around pairing the source and drain endpoints of the segment $s$ of $k-1$ points that is obtained when removing the point $n$ from a single cycle. It is convenient to use indices $2$ to $k$ for these points. For the replica graph, $\tilde{p}$ represents the points $\tilde{2}$ (now the drain) and $\tilde{k}$ (now the source). We are interested in the case where the node degrees and the cycle length $k$ are small relative to the graph size $N$, and nodes are assigned to cycles at random, such that the local topology approaches a tree-like digraph.

The Green’s function of the graph $M$ on nodes $n$ and $\tilde{n}$ can then be written as

[TABLE]

where $G^{(n,m)}=(\mathcal{Z}-\mathcal{N})^{-1}$ is the Green’s function of the graph with $n$ removed expanded around .

Truncated to the segments $s$ and $\tilde{s}$, this Green’s function can be written as

[TABLE]

where $\mathcal{S}=\mathrm{diag}\,(S,S^{T})$ with $S_{i,j}=\delta_{i,j+1}$ represents the graph on the isolated segment $s$. This gives the condition

[TABLE]

where, utilizing that $G^{(n)}_{\tilde{2},\tilde{k}}={G^{(n)}_{k,2}}^{*}$, we find

[TABLE]

On the boundary of the support of the spectrum, the off-diagonal elements $G^{(n)}_{k,\tilde{k}}$ and $G^{(n)}_{\tilde{2},2}$ vanish, while we always need to consider the limit $\lambda\to 0$. This means that $\beta$ and $\gamma$ become small, giving the expansion

[TABLE]

where the matrix $A$ has elements $A_{nm}=\alpha^{-1-n+m}$ if $n\geq m$ and 0 otherwise. Comparing the matrix element $k,\tilde{k}$ on both sides gives the condition

[TABLE]

so that we can write $\alpha=t^{-1}e^{-i\varphi}$ where $t$ is the unique positive real solution of $(d-1)t^{2k}-dt^{2}+1=0$. Comparing the matrix element $k,2$ on both sides of (27) we furthermore have the condition $G^{(n)}_{k,2}=\alpha^{1-k}$. Given the definition of $\alpha$ in Eq. (25), we can express $z=\alpha+(d-1)G^{(n)}_{k,2}$, so that the boundary of the support is parameterized as

[TABLE]

This again describes a hypotrochoid curve, and is in agreement with numerical results as presented in Fig. 3.

Appendix C Appendix C: Derivation for graphs with multiple cycle structures

For a graph where each vertex connects to $d_{1}$ $k_{1}$-cycles of weight $w_{1}$ and $d_{2}$ $k_{2}$-cycles of weight $w_{2}$, the Green’s functions on a node $n$ and its replica $\tilde{n}$ can be written as

[TABLE]

with

[TABLE]

where the isolated segments $s_{i}$ with adjacency matrix $S^{(i)}$ now occupy points $2$ to $k_{i}$.

These matrices are still of the form (24), where now

[TABLE]

On the boundary it now is consistent to require that both $\beta_{1}$ and $\beta_{2}$ are small. Inverting the matrices perturbatively as before, we then find the following conditions,

[TABLE]

where the last two require

[TABLE]

with $\Sigma_{r}=\sum_{l=1}^{k_{r}-1}|\alpha_{r}/w_{r}|^{-2l}$. Furthermore, using the definitions of $\alpha_{1}$ and $\alpha_{2}$ in Eq. (33) and the first two conditions in Eq. (34) we obtain

[TABLE]

which have to agree, hence

[TABLE]

We now set $\alpha_{r}=(w_{r}/t_{r})\exp(-i\varphi_{r})$, so that

[TABLE]

Equation (35) then relates $t_{2}$ to $t_{1}$, while Eq. (38) further relates $\varphi_{2}$ to $\varphi_{1}$. Given these relations, the boundary curve can then be written, in symmetrised form, as

[TABLE]

Explicit expressions are obtained for $d_{1},d_{2}\gg 1$. Then, conditions (35) and (38) give in leading order $\varphi_{2}=\varphi_{1}\equiv\varphi$ and

[TABLE]

so that

[TABLE]

where $\bar{d}=\sqrt{d_{1}w_{1}^{2}+d_{2}w_{2}^{2}}$.