Spectral density of complex networks with two species of nodes

The theory of complex networks, which has been dramatically developed since the end of the last century, is based on the observation that there are universal features in real biological and social networks [1]. One such feature is the scale-free property, meaning that the degree (the number of nodes directly connected to each node) distribution function P(Δ) obeys a power law P(Δ)∝Δ^−λ for large Δ. Barabási and Albert explained the origin of this scale-free property by focusing on the network growing process [2]. Goh, Kahng and Kim formulated a static network model which exhibits the scale-free property [3].

The connection pattern of a network is mathematically described by the adjacency matrix. When the network has the scale-free property, the spectral (eigenvalue) density ρ(μ) of the adjacency matrix is also expected to obey a power law ρ(μ)∝μ^−γ for large μ. Dorogovtsev et al presented an analytic evidence of this power law behaviour [4, 5]. Moreover, a relation γ = 2λ − 1 was found between the exponents of the power laws. Rodgers et al analysed Goh, Kahng and Kim's static model and confirmed the power law behaviour of ρ(μ) [6].

In this paper, we shall study scale-free networks with two species (A and B) of nodes. Connections are made only between the nodes of different species. We introduce a static model of such bipartite scale-free networks by applying Goh, Kahng and Kim's method, and observe that each species has its own degree distribution function obeying a power law. Suppose that the exponent of the degree distribution function is λ_A for the species A and λ_B for the species B. Using the replica method in statistical physics, we are able to analytically evaluate the spectral density ρ(μ) in the limit of large mean degree [6–8]. As a result, we find that ρ(μ) also obeys a power law and the exponent γ is associated with the exponents λ_A and λ_B. In addition, the spectral density of the Laplacian matrix is similarly analysed and the power law behaviour is confirmed.

Bipartite networks find applications in the analysis of human sexual contacts [9], and of the connections between collaborators and collaboration acts [10], such as actors and movies, scientists and papers. The adjacency matrices of bipartite networks are also interesting from the viewpoint of random matrix theory, since the Gaussian matrix model with the same structure is called the chiral Gaussian ensemble and applied to physics, such as the QCD gauge theory [11].

The outline of this paper is as follows. In section 2, a static model of bipartite scale-free networks with two species of nodes is introduced, and the adjacency and Laplacian matrices are defined. In section 3, in order to evaluate the spectral density, we apply the replica method to the network model. In the limit of large mean degree, the power law behaviour of the spectral density is analytically derived. In section 4, the effective medium approximation is briefly discussed as an attempt to treat the case with a finite mean degree.

Let us suppose that there are N nodes of type A and M nodes of type B (N ⩾ M). We are interested in the asymptotic behaviour of bipartite networks with two species of nodes A and B in the limit

$$\begin{equation} N \rightarrow \infty \ {\rm and} \ M \rightarrow \infty \ {\rm with} \ c = N/M \ {\rm fixed}. \end{equation} \tag{ 2.1 }$$

We introduce a static model of such networks with power law degree distributions by applying Goh, Kahng and Kim's method. Each node of type A is assigned a probability P_j normalized as

$$\begin{equation} \sum _{j=1}^N P_j = 1, \end{equation} \tag{ 2.2 }$$

while each node of type B has a probability Q_k with

$$\begin{equation} \sum _{k=1}^M Q_k = 1. \end{equation} \tag{ 2.3 }$$

The nodes of type A and B are connected according to the following procedure. In each step, we choose a node j of type A and a node k of type B with probabilities P_j and Q_k, respectively. Then the nodes j and k are connected, unless they are already connected. After repeating such a step pN times, a node j of type A and a node k of type B are connected with the probability

$$\begin{equation} f_{jk} = 1 -(1 - P_j Q_k)^{p N} \sim 1 - {\rm e}^{- p N P_j Q_k}. \end{equation} \tag{ 2.4 }$$

Let us consider an N × M matrix C (N ⩾ M), where C_jk = 1 if the node j of type A is directly connected to the node k of type B, and C_jk = 0 otherwise. This random matrix C describes the connection pattern of the network with two species of nodes. Each matrix element C_jk is independently distributed with the probability density function (p.d.f.)

$$\begin{equation} \mathcal{P}_{jk}(C_{jk}) = (1 - f_{jk}) \delta (C_{jk}) + f_{jk} \delta (1 - C_{jk}). \end{equation} \tag{ 2.5 }$$

We assume that P_j and Q_k are given by

$$\begin{equation} P_j = \frac{j^{-\alpha }}{ \sum _{l=1}^N l^{-\alpha }} \sim (1 - \alpha ) N^{\alpha - 1} j^{-\alpha }, \ \ \ 0 < \alpha < 1 \end{equation} \tag{ 2.6 }$$

and

$$\begin{equation} Q_k = \frac{k^{-\beta }}{ \sum _{l=1}^M l^{-\beta }} \sim (1 - \beta ) M^{\beta - 1} k^{-\beta }, \ \ \ 0 < \beta < 1. \end{equation} \tag{ 2.7 }$$

There are thus two parameters α and β controlling the p.d.f. of the matrix C.

We define the degree d_j of the type-A node j as the number of directly connected type-B nodes:

$$\begin{equation} d_j = \sum _{k=1}^M C_{jk}. \end{equation} \tag{ 2.8 }$$

Then the type-A node degree distribution function is given by

$$\begin{equation} P^{(A)}(\Delta ) = \left\langle \frac{1}{N} \sum _{j=1}^N \delta (\Delta - d_j) \right\rangle , \end{equation} \tag{ 2.9 }$$

where the brackets denote the average over the p.d.f. (2.5) and δ(x) is Dirac's delta function. We can similarly introduce the degree e_k of the type-B node k:

$$\begin{equation} e_k = \sum _{j=1}^N C_{jk} \end{equation} \tag{ 2.10 }$$

and the type-B node degree distribution function

$$\begin{equation} P^{(B)}(\Delta ) = \left\langle \frac{1}{M} \sum _{k=1}^M \delta (\Delta - e_k) \right\rangle . \end{equation} \tag{ 2.11 }$$

In appendix A, a useful asymptotic relation

$$\begin{equation} \ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle \sim p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k({\rm e}^{ - {\rm i} t_{jk}} - 1) \end{equation} \tag{ 2.12 }$$

is derived in the limit (2.1). Here t_jk, which depends neither on N nor on M, is in the neighbourhood of the origin so that $|{\rm e}^{ - {\rm i} t_{jk}} - 1| < 1$.

As special cases, we can readily derive asymptotic relations for

$$\begin{equation} F^{(A)}_j(t) \mathrel {\mathop :}= \ln \langle {\rm e}^{- {\rm i} d_j t} \rangle , \ \ \ F^{(B)}_k(t) \mathrel {\mathop :}= \ln \langle {\rm e}^{- {\rm i} e_k t} \rangle \end{equation} \tag{ 2.13 }$$

as

$$\begin{equation} F^{(A)}_j(t) \sim p N P_j ({\rm e}^{- {\rm i} t} - 1), \ \ \ F^{(B)}_k(t) \sim p N Q_k ({\rm e}^{- {\rm i} t} - 1). \end{equation} \tag{ 2.14 }$$

Then we can readily see that

$$\begin{eqnarray} \langle d_j \rangle & = & {\rm i} \left. \frac{\partial }{\partial t} F^{(A)}_j(t) \right|_{t=0} \sim p N P_j, \nonumber \\ \langle e_k \rangle & = & {\rm i} \left. \frac{\partial }{\partial t} F^{(B)}_k(t) \right|_{t=0} \sim p N Q_k, \end{eqnarray} \tag{ 2.15 }$$

so that the mean degree m^(A) of the type-A node is

$$\begin{equation} m^{(A)} = \frac{1}{N} \sum _{j=1}^N \langle d_j \rangle \sim p, \end{equation} \tag{ 2.16 }$$

while the mean degree m^(B) of the type-B node is

$$\begin{equation} m^{(B)} = \frac{1}{M} \sum _{k=1}^M \langle e_k \rangle \sim p c. \end{equation} \tag{ 2.17 }$$

It can be seen from (2.6), (2.9) and (2.14) that the type-A node degree distribution function can be written as

$$\begin{eqnarray} P^{(A)}(\Delta ) & = & \frac{1}{2 \pi N} \sum _{j=1}^N \int {\rm d}t \ {\rm e}^{{\rm i} \Delta t + F^{(A)}_j(t)} \nonumber \\ && \sim \frac{1}{2 \pi } \int {\rm d}t \int _0^1 {\rm d}x \ {\rm exp}\lbrace {\rm i} \Delta t + p (1 - \alpha ) x^{-\alpha } ({\rm e}^{- {\rm i} t} - 1) \rbrace . \end{eqnarray} \tag{ 2.18 }$$

Then in the limit Δ → ∞, we find

$$\begin{equation} P^{(A)}(\Delta ) \sim \int _0^1 {\rm d}x \ \delta \lbrace \Delta - p (1 - \alpha ) x^{-\alpha } \rbrace = \frac{\lbrace p (1 - \alpha ) \rbrace ^{1/\alpha }}{\alpha } \frac{1}{\Delta ^{(1/\alpha ) + 1}}, \end{equation} \tag{ 2.19 }$$

and similarly obtain

$$\begin{equation} P^{(B)}(\Delta ) \sim \int _0^1 {\rm d}y \ \delta \lbrace \Delta - p c (1 - \beta ) y^{-\beta } \rbrace = \frac{\lbrace p c (1 - \beta ) \rbrace ^{1/\beta }}{\beta } \frac{1}{\Delta ^{(1/\beta ) + 1}}. \end{equation} \tag{ 2.20 }$$

Thus we have seen that the network has the scale-free property, as the node degree distribution functions obey power laws. The exponents of the power laws defined as

$$\begin{equation} P^{(A)}(\Delta ) \propto \Delta ^{-\lambda _A}, \ \ \ P^{(B)}(\Delta ) \propto \Delta ^{-\lambda _B}, \ \ \ \Delta \rightarrow \infty \end{equation} \tag{ 2.21 }$$

are found to be λ_A = (1/α) + 1 and λ_B = (1/β) + 1.

In this paper, we study the adjacency and Laplacian matrices of this scale-free network. The adjacency matrix $\mathcal{A}$ of this network is defined as

$$\begin{equation} \mathcal{A} = \left(\begin{array}{@{}cc@{}}O_N & C \\ C^{\rm T} & O_M\end{array} \right), \end{equation} \tag{ 2.22 }$$

where C^T is the transpose of C and O_n is an n × n matrix with zero elements. The Laplacian matrix $\mathcal{L}$ is an (N + M) × (N + M) symmetric matrix with

$$\begin{equation} \mathcal{L}_{jl} = \left\lbrace \begin{array}{@{}ll@{}}d_j, & j = l \ {\rm and} \ 1 \le j \le N, \\ e_{j-N}, & j = l \ {\rm and} \ N+1 \le j \le N+M, \\ - \mathcal{A}_{jl}, & j \ne l.\end{array} \right. \end{equation} \tag{ 2.23 }$$

Let us define that J is the adjacency matrix $\mathcal{A}$ or the Laplacian matrix $\mathcal{L}$. The spectral density of J is defined as

$$\begin{equation} \rho (\mu ) = \left\langle \frac{1}{N+M} \sum _{j=1}^{N+M} \delta (\mu - \mu _j) \right\rangle , \end{equation} \tag{ 3.1 }$$

where μ_j, j = 1, 2, ..., N + M, are the eigenvalues of J. In order to calculate ρ(μ), we introduce the partition function

$$\begin{equation} Z(\mu ) = \int ^{\infty }_{-\infty } \prod _{j=1}^{N+M} {\rm d}\Phi _j \ {\rm exp}\left( \frac{{\rm i}}{2} \mu \sum _{j=1}^{N+M} \Phi _j^2 - \frac{{\rm i}}{2} \sum _{j=1}^{N+M} \sum _{l=1}^{N+M} J_{jl} \Phi _j \Phi _l \right). \end{equation} \tag{ 3.2 }$$

Using the partition function Z, we can write the spectral density as

$$\begin{eqnarray} \rho (\mu ) & = & \frac{1}{(N+M) \pi } {\rm Im}{\rm Tr}\langle \lbrace J - (\mu + {\rm i} \epsilon ) I\rbrace ^{-1} \rangle \nonumber \\ & = & \frac{2}{(N+M) \pi } {\rm Im} \frac{\partial }{\partial \mu } \langle \ln Z(\mu + {\rm i }\epsilon ) \rangle , \end{eqnarray} \tag{ 3.3 }$$

where is an infinitesimal positive number and I is an (N + M) × (N + M) identity matrix. Then we can utilize the relation

$$\begin{equation} \lim _{n \rightarrow 0} \frac{\ln \langle Z^n \rangle }{n} = \langle \ln Z \rangle \end{equation} \tag{ 3.4 }$$

to obtain

$$\begin{equation} \rho (\mu ) = \lim _{n \rightarrow 0} \frac{2}{(N+M) n \pi } {\rm Im} \frac{\partial }{\partial \mu } \ln \langle \lbrace Z(\mu ) \rbrace ^n \rangle . \end{equation} \tag{ 3.5 }$$

Therefore, it is necessary to evaluate the average 〈Zⁿ〉.

The replica method explained in appendix B is known to be a powerful tool for that purpose. It follows in the limit (2.1) that

$$\begin{equation} \langle Z^n \rangle \sim \int \prod _{j=1}^N \mathcal{D}\xi _j({\vec{\phi }}) \prod _{k=1}^M \mathcal{D}\eta _k({\vec{\psi }}) \ {\rm e}^{S_0 + S_1 + S_2}. \end{equation} \tag{ 3.6 }$$

Here,

$$\begin{equation} S_0 = - \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) \ln \xi _j({\vec{\phi }}) - \sum _{k=1}^M \int {\rm d}{\vec{\psi }} \ \eta _k({\vec{\psi }}) \ln \eta _k({\vec{\psi }}), \end{equation} \tag{ 3.7 }$$

$$\begin{equation} S_1 = \frac{{\rm i}}{2} \mu \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) {\vec{\phi }}^2 + \frac{{\rm i}}{2} \mu \sum _{k=1}^M \int {\rm d}{\vec{\psi }} \ \eta _k({\vec{\psi }}) {\vec{\psi }}^2 \end{equation} \tag{ 3.8 }$$

and

$$\begin{equation} S_2 = p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k \int {\rm d}{\vec{\phi }} \int {\rm d}{\vec{\psi }} \ \xi _j({\vec{\phi }}) \ \eta _k({\vec{\psi }}) \lbrace f({\vec{\phi }},{\vec{\psi }}) - 1 \rbrace \end{equation} \tag{ 3.9 }$$

with

$$\begin{equation} f({\vec{\phi }},{\vec{\psi }}) = \left\lbrace \begin{array}{@{}ll@{}}{\rm e}^{- {\rm i} {\vec{\phi }} \cdot {\vec{\psi }}}, & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm adjacency} \ {\rm matrix} \ \mathcal{A}, \\ {\rm e}^{- \frac{{\rm i}}{2}({\vec{\phi }} -{\vec{\psi }})^2}, & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm Laplacian} \ {\rm matrix} \ \mathcal{L}. \end{array} \right. \end{equation} \tag{ 3.10 }$$

The functional integrations are taken over the auxiliary functions $\xi _j({\vec{\phi }})$ and $\eta _k({\vec{\psi }})$ satisfying

$$\begin{equation} \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) = \int {\rm d}{\vec{\psi }} \ \eta _k({\vec{\psi }}) = 1. \end{equation} \tag{ 3.11 }$$

In the limit (2.1), the functional integrations over $\xi _j({\vec{\phi }})$ and $\eta _k({\vec{\psi }})$ are dominated by the stationary point satisfying

$$\begin{equation} \fl\delta \left\lbrace S_0 + S_1 + S_2 + \sum _{j=1}^N \theta _j \left( \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) - 1 \right) + \sum _{k=1}^M \omega _k \left( \int {\rm d}{\vec{\psi }} \ \eta _k({\vec{\psi }}) - 1 \right) \right\rbrace = 0, \end{equation} \tag{ 3.12 }$$

where θ_j and ω_k are the Lagrange multipliers. It follows from this equation that

$$\begin{eqnarray} \xi _j({\vec{\phi }}) & = & \Theta _j \ {\rm exp}\left[ \frac{{\rm i}}{2} \mu {\vec{\phi }}^2 + p N P_j \sum _{k=1}^M Q_k \int {\rm d}{\vec{\psi }} \eta _k({\vec{\psi }}) \lbrace f({\vec{\phi }},{\vec{\psi }}) - 1 \rbrace \right], \nonumber \\ \eta _k({\vec{\psi }}) & = & \Omega _k \ {\rm exp}\left[ \frac{{\rm i}}{2} \mu {\vec{\psi }}^2 + p N Q_k \sum _{j=1}^N P_j \int {\rm d}{\vec{\phi }} \xi _j({\vec{\phi }}) \lbrace f({\vec{\phi }},{\vec{\psi }}) - 1 \rbrace \right], \end{eqnarray} \tag{ 3.13 }$$

where Θ_j and Ω_k are normalization constants.

In the limit of large mean degree p → ∞, the variational equations (3.13) are satisfied by the Gaussian ansatz

$$\begin{equation} \fl\xi _j({\vec{\phi }}) = \frac{1}{(2 \pi {\rm i} \sigma _j)^{n/2}} \ {\rm exp} \left( - \frac{{\vec{\phi }}^2}{2 {\rm i} \sigma _j} \right), \ \ \ \eta _k({\vec{\psi }}) = \frac{1}{(2 \pi {\rm i} \tau _k)^{n/2}} \ {\rm exp} \left( - \frac{{\vec{\psi }}^2}{2 {\rm i} \tau _k} \right), \end{equation} \tag{ 3.14 }$$

as shown in appendix C. Here, Im σ_j ⩽ 0 and Im τ_k ⩽ 0. This property simplifies the problem and enables us to evaluate the asymptotic spectral density.

Let us consider the limit p → ∞ with a scaling variable $E = \mu /\sqrt{p}$. In appendix C, we find the asymptotic spectral density of the adjacency matrix $\mathcal{A}$ as

$$\begin{equation} \rho (\mu ) \sim \frac{2}{\sqrt{p} (1 + c)} \left\lbrace \frac{c^{1/\beta } (1 - \beta )^{1/\beta }}{\beta E^{(2/\beta )+ 1}} + \frac{c (1 - \alpha )^{1/\alpha }}{\alpha E^{(2/\alpha )+ 1}} \right\rbrace \end{equation} \tag{ 3.15 }$$

in the tail region E → ∞. The exponent γ of the spectral density defined as

$$\begin{equation} \rho (\mu ) \propto \mu ^{-\gamma }, \ \ \ \mu \rightarrow \infty , \end{equation} \tag{ 3.16 }$$

is (2/α) + 1 if α ⩾ β, and is (2/β) + 1 if β ⩾ α. Thus γ is associated with λ_A = (1/α) + 1 and λ_B = (1/β) + 1 as γ = 2min (λ_A, λ_B) − 1.

It is also explained in appendix C that the asymptotic spectral density of the Laplacian matrix $\mathcal{L}$ is given by

$$\begin{eqnarray} &&\rho (\mu ) \sim \frac{c \lbrace p (1 - \alpha )\rbrace ^{1/\alpha }}{(1 + c) \alpha } \frac{1}{\mu ^{(1/\alpha )+1}} \ H\lbrace \mu - p (1 - \alpha ) \rbrace \nonumber \\ &&\hphantom{ \rho (\mu ) } + \frac{\lbrace p c (1 - \beta )\rbrace ^{1/\beta }}{(1 + c) \beta } \frac{1}{\mu ^{(1/\beta )+1}} \ H\lbrace \mu - p c (1 - \beta ) \rbrace \end{eqnarray} \tag{ 3.17 }$$

in the region μ = O(p) with p → ∞. Here, H(x) is defined as

$$\begin{equation} H(x) = \left\lbrace \begin{array}{@{}ll@{}}0, & x < 0, \\ 1, & x > 0.\end{array} \right. \end{equation} \tag{ 3.18 }$$

The exponent γ_L of the spectral density $\rho (\mu ) \propto \mu ^{-\gamma _L}$ (μ → ∞) is (1/α) + 1 if α ⩾ β, and is (1/β) + 1 if β ⩾ α. Thus γ_L is associated with λ_A and λ_B as γ_L = min (λ_A, λ_B).

In the previous section, we have dealt with the spectral density in the limit p → ∞. The calculation of the spectral density with a finite mean degree p is a much more involved problem, for which sophisticated numerical schemes have been proposed [12–15]. In this section, we briefly discuss a simple approximation method (effective medium approximation) for that problem [8, 16–20]. In this approximation, we put the Gaussian ansatz (3.14) into formulas (3.7), (3.8) and (3.9), and solve the stationary point equations

$$\begin{equation} \frac{\partial }{\partial \sigma _j} (S_0 + S_1 + S_2) = 0 \end{equation} \tag{ 4.1 }$$

and

$$\begin{equation} \frac{\partial }{\partial \tau _k} (S_0 + S_1 + S_2) = 0. \end{equation} \tag{ 4.2 }$$

In the case of the adjacency matrix $\mathcal{A}$, the above procedure results in the effective medium approximation (EMA) equations

$$\begin{eqnarray} & & \mu - \frac{1}{\sigma _j} - p N P_j \sum _{k=1}^M \frac{Q_k \tau _k}{1 - \sigma _j \tau _k} = 0, \nonumber \\\ms & & \mu - \frac{1}{\tau _k} - p N Q_k \sum _{j=1}^N \frac{P_j \sigma _j}{1 - \sigma _j \tau _k} = 0. \end{eqnarray} \tag{ 4.3 }$$

As for scale-free networks with a single species of nodes, Nagao and Rodgers calculated the 1/p expansion of the spectral density by using the corresponding EMA equation [20]. A similar analytical treatment could also be possible in the present case. Here, however only results of numerical iterations of (4.3) are shown in figure 1 as the EMA spectral densities. They are compared with the spectral densities of positive eigenvalues calculated by numerical diagonalizations of numerically generated adjacency matrices (averaged over 100 samples). The EMA gives a better fit for a larger p, as expected from the fact that the variational equations (3.13) are satisfied by the Gaussian ansatz (3.14) in the limit p → ∞. When p = 1, the agreement significantly breaks down around the origin, although it is still fairly good in the tail region with large μ.

Zoom In Zoom Out Reset image size

In the limit α, β → 0, we obtain the adjacency matrix of a classical random graph with two species, where the connections are made only between the nodes of different species. In that case, σ_j and τ_k can be written as σ and τ, respectively, because they depend neither on j nor on k. The EMA equations become a cubic equation for σ

$$\begin{equation} c \sigma ^3 + \frac{c (p - 2) + 1}{\mu } \sigma ^2 - \frac{\mu ^2 + (c - 1)(p - 1)}{\mu ^2} \sigma + \frac{1}{\mu } = 0 \end{equation} \tag{ 4.4 }$$

and

$$\begin{equation} \tau = \frac{c (\mu \sigma - 1) + 1}{\mu }. \end{equation} \tag{ 4.5 }$$

These equations are equivalent to Nagao and Tanaka's SEMA (symmetric EMA) equations concerning the spectral density of sparse correlation matrices [19], and can be analysed in the same way.

We can similarly derive the EMA equations for the Laplacian matrix $\mathcal{L}$ as

$$\begin{eqnarray} & & \mu - \frac{1}{\sigma _j} - p N P_j \sum _{k=1}^M \frac{Q_k}{1 - \sigma _j - \tau _k} = 0, \nonumber \\ & & \mu - \frac{1}{\tau _k} - p N Q_k \sum _{j=1}^N \frac{P_j}{1 - \sigma _j - \tau _k} = 0. \end{eqnarray} \tag{ 4.6 }$$

In the limit α, β → 0, σ_j and τ_k can again be reduced to σ and τ, respectively. Then we find a cubic equation for σ

$$\begin{equation} \fl(1 - c) \sigma ^3 + \frac{2 c + (p - \mu )(1 - c)}{\mu } \sigma ^2 + \frac{\mu - 1 + c (p - 2 \mu - 1)}{\mu ^2} \sigma + \frac{c}{\mu ^2} = 0 \end{equation} \tag{ 4.7 }$$

and

$$\begin{equation} \tau = \frac{\sigma }{(1 - c)\mu \sigma + c}. \end{equation} \tag{ 4.8 }$$

The author thanks Professor G J Rodgers and Professor Toshiyuki Tanaka for valuable discussions. This work was partially supported by the Japan Society for the Promotion of Science (KAKENHI 20540372).

In this appendix, we derive an asymptotic relation

$$\begin{equation} \ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle \sim p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k ({\rm e}^{ - {\rm i} t_{jk}} - 1), \end{equation} \tag{ A.1 }$$

where t_jk is a parameter which is independent of N and M. We moreover assume that t_jk is in the neighbourhood of the origin so that |S_jk| < 1 holds for $S_{jk} = {\rm e}^{ - {\rm i} t_{jk}} - 1$. A similar argument for Goh, Kahng and Kim's model is found in [21].

The Taylor expansion of the logarithmic function gives

$$\begin{eqnarray} &&\fl\ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle = \sum _{j=1}^N \sum _{k=1}^M \ln (1 + f_{jk} S_{jk}) \nonumber \\ &&{\hphantom{\fl\ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle}} = p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k S_{jk} + \sum _{j=1}^N \sum _{k=1}^M (f_{jk} - p N P_j Q_k) S_{jk} \nonumber \\ &&{\hphantom{\fl\ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle}}\quad+ \sum _{\ell =2}^{\infty } \frac{(-1)^{\ell +1}}{\ell } \sum _{j=1}^N \sum _{k=1}^M (f_{jk} S_{jk})^\ell . \end{eqnarray} \tag{ A.2 }$$

We show

$$\begin{equation} \ln \left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} t_{jk} \right) \right\rangle - p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k S_{jk} = o(N) \end{equation} \tag{ A.3 }$$

in two steps.

Step 1

Let us first prove that

$$\begin{equation} \sum _{j=1}^N \sum _{k=1}^M (f_{jk} - p N P_j Q_k) S_{jk} = o(N). \end{equation} \tag{ A.4 }$$

We define

$$\begin{equation} S_{\max } = \max _{jk} | S_{jk} | \end{equation} \tag{ A.5 }$$

and

$$\begin{equation} G_1(x) = x - 1 + {\rm e}^{-x}. \end{equation} \tag{ A.6 }$$

Then we see that

$$\begin{eqnarray} \left| \sum _{j=1}^N \sum _{k=1}^M (f_{jk} - p N P_j Q_k) S_{jk} \right| & \le & S_{\max } \sum _{j=1}^N \sum _{k=1}^M |f_{jk} - p N P_j Q_k| \nonumber \\ & = & S_{\max } \sum _{j=1}^N \sum _{k=1}^M G_1(p N P_j Q_k). \end{eqnarray} \tag{ A.7 }$$

A monotonously decreasing continuous function F(x) satisfies

$$\begin{equation} \sum _{k=1}^M F(k) \le \int _1^M F(y) \ {\rm d}y + F(1), \end{equation} \tag{ A.8 }$$

so that

$$\begin{equation} \sum _{k=1}^M G_1(p N P_j Q_k) \le \int _1^M \ G_1\left(p N P_j \frac{y^{-\beta }}{\zeta _M(\beta )} \right) {\rm d}y + G_1\left( p N P_j \frac{1}{\zeta _M(\beta )} \right) \end{equation} \tag{ A.9 }$$

with

$$\begin{equation} \zeta _M(\beta ) = \sum _{\ell =1}^M \ell ^{-\beta }. \end{equation} \tag{ A.10 }$$

Then one can again use (A.8) to obtain

$$\begin{equation} \sum _{j=1}^N \sum _{k=1}^M G_1(p N P_j Q_k) \le \sum _{\nu = 1}^4 I_{\nu }, \end{equation} \tag{ A.11 }$$

where

$$\begin{eqnarray} I_1 & = & \int _1^N {\rm d}x \int _1^M {\rm d}y \ G_1\left(p N \frac{x^{-\alpha }}{\zeta _N(\alpha )} \frac{y^{-\beta }}{\zeta _M(\beta )} \right), \nonumber \\ I_2 & = & \int _1^M {\rm d}y \ G_1\left( p N \frac{1}{\zeta _N(\alpha )} \frac{y^{-\beta }}{\zeta _M(\beta )} \right), \ \ \ I_3 = \int _1^N {\rm d}x \ G_1\left( p N \frac{x^{-\alpha }}{\zeta _N(\alpha )} \frac{1}{\zeta _M(\beta )} \right), \nonumber \\ I_4 & = & G_1\left( p N \frac{1}{\zeta _N(\alpha )} \frac{1}{\zeta _M(\beta )} \right). \end{eqnarray} \tag{ A.12 }$$

Using the notations

$$\begin{equation} \epsilon _1 = \frac{\sqrt{p N}}{\zeta _N(\alpha )} N^{-\alpha }, \ \ \ \epsilon _2 = \frac{\sqrt{p N}}{\zeta _M(\beta )} M^{-\beta }, \end{equation} \tag{ A.13 }$$

we see that

$$\begin{equation} I_1 = \frac{MN \epsilon _1^{1/\alpha } \epsilon _2^{1/\beta }}{\alpha \beta } \int _{\epsilon _1}^{\epsilon _1 N^{\alpha }} {\rm d}u \int _{\epsilon _2}^{\epsilon _2 M^{\beta }} {\rm d}v \ \frac{G_1(uv)}{u^{1 + (1/\alpha )} v^{1 + (1/\beta )}}. \end{equation} \tag{ A.14 }$$

Then, using the inequality

$$\begin{equation} G_1(x) \le x^2/2, \ \ \ x \ge 0, \end{equation} \tag{ A.15 }$$

we find

$$\begin{eqnarray} I_1 & \le & \frac{MN \epsilon _1^{1/\alpha } \epsilon _2^{1/\beta }}{2 \alpha \beta } \int _{\epsilon _1}^{\epsilon _1 N^{\alpha }} {\rm d}u \int _{\epsilon _2}^{\epsilon _2 M^{\beta }} {\rm d}v \ u^{1 - (1/\alpha )} v^{1 - (1/\beta )} \nonumber \\ & = & O(N^{(2 \alpha - 1)} N^{(2 \beta - 1)}), \end{eqnarray} \tag{ A.16 }$$

where

$$\begin{equation} N^{(2 \alpha - 1)} \equiv \left\lbrace \begin{array}{@{}ll@{}}1, & 0 < \alpha < 1/2, \\ \ln N, & \alpha = 1/2, \\ N^{2 \alpha - 1}, & 1/2 < \alpha < 1.\end{array} \right. \end{equation} \tag{ A.17 }$$

In the case min (α, β) > 1/2, we similarly employ

$$\begin{equation} G_1(x) \le \left\lbrace \begin{array}{@{}ll@{}}x^2/2, & 0 \le x \le 1, \\ x, & x \ge 1 \end{array} \right. \end{equation} \tag{ A.18 }$$

to obtain

$$\begin{equation} I_1 = O(N^{(\alpha ,\beta )}), \end{equation} \tag{ A.19 }$$

where

$$\begin{equation} N^{(\alpha ,\beta )} = \left\lbrace \begin{array}{@{}ll@{}}N^{(\alpha + \beta - 1)/\min (\alpha ,\beta )}, & \alpha \ne \beta , \\ N^{(2 \alpha - 1)/\alpha } \ln N, & \alpha = \beta .\end{array} \right. \end{equation} \tag{ A.20 }$$

Using inequality (A.15), we can similarly derive the estimates

$$\begin{equation} I_2 = O(N^{2 \alpha - 1} N^{(2 \beta - 1)}), \ \ \ I_3 = O(N^{(2 \alpha - 1)} N^{2 \beta - 1}). \end{equation} \tag{ A.21 }$$

If min (α, β) > 1/2, then we utilize (A.18) to find

$$\begin{equation} I_2 = O(N^{(\alpha + \beta - 1)/\beta }), \ \ \ I_3 = O(N^{(\alpha + \beta - 1)/\alpha }). \ \ \ \end{equation} \tag{ A.22 }$$

Moreover, one can readily see from the inequality G₁(x) ⩽ x (x ⩾ 0) that

$$\begin{equation} I_4 = O(N^{\alpha + \beta - 1}). \end{equation} \tag{ A.23 }$$

It follows from (A.16), (A.19), (A.21), (A.22) and (A.23) that

$$\begin{equation} \sum _{j=1}^N \sum _{k=1}^M G_1(p N P_j Q_k) \le \sum _{\nu = 1}^4 I_{\nu } = o(N), \end{equation} \tag{ A.24 }$$

which yields (A.4).

Step 2

We next prove

$$\begin{equation} \sum _{\ell =2}^{\infty } \frac{(-1)^{\ell +1}}{\ell} \sum _{j=1}^N \sum _{k=1}^M (f_{jk} S_{jk})^\ell = o(N). \end{equation} \tag{ A.25 }$$

Using

$$\begin{equation} G_0(x) = 1 - {\rm e}^{-x}, \end{equation} \tag{ A.26 }$$

we see that

$$\begin{eqnarray} \left| \sum _{\ell =2}^{\infty } \frac{(-1)^{\ell +1}}{\ell} \sum _{j=1}^N \sum _{k=1}^M (f_{jk} S_{jk})^\ell \right| & \le & \sum _{\ell =2}^{\infty } \frac{S_{\max }^{\ell }}{\ell } \sum _{j=1}^N \sum _{k=1}^M f_{jk}^{\ell } \nonumber \\ & = & \sum _{\ell =2}^{\infty } \frac{S_{\max }^{\ell }}{\ell } \sum _{j=1}^N \sum _{k=1}^M G_0(p N P_j Q_k)^{\ell }. \end{eqnarray} \tag{ A.27 }$$

We can again employ (A.8) to obtain

$$\begin{equation} \sum _{j=1}^N \sum _{k=1}^M G_0(p N P_j Q_k)^{\ell } \le \sum _{\nu = 1}^4 J_{\nu }, \end{equation} \tag{ A.28 }$$

where

$$\begin{eqnarray} \fl J_1 & = & \int _1^N {\rm d}x \int _1^M {\rm d}y \ G_0\left(p N \frac{x^{-\alpha }}{\zeta _N(\alpha )} \frac{y^{-\beta }}{\zeta _M(\beta )} \right)^{\ell }, \nonumber \\ \fl J_2 & = & \int _1^M {\rm d}y \ G_0\left( p N \frac{1}{\zeta _N(\alpha )} \frac{y^{-\beta }}{\zeta _M(\beta )} \right)^{\ell }, \ \ \ J_3 = \int _1^N {\rm d}x \ G_0\left( p N \frac{x^{-\alpha }}{\zeta _N(\alpha )} \frac{1}{\zeta _M(\beta )} \right)^{\ell }, \nonumber \\ \fl J_4 & = & G_0\left( p N \frac{1}{\zeta _N(\alpha )} \frac{1}{\zeta _M(\beta )} \right)^{\ell }. \end{eqnarray} \tag{ A.29 }$$

Making use of the identity

$$\begin{equation} G_0(x) \le x, \ \ \ x \ge 0, \end{equation} \tag{ A.30 }$$

we find

$$\begin{equation} J_1 = O( N^{\langle \alpha ,\ell \rangle } N^{\langle \beta ,\ell \rangle } ), \end{equation} \tag{ A.31 }$$

where

$$\begin{equation} N^{\langle \alpha ,\ell \rangle } = \left\lbrace \begin{array}{@{}ll@{}}N^{(2 - \ell )/2}, & \ell < 1/\alpha , \\ N^{(2 \alpha - 1)/(2 \alpha )} \ln N, & \ell = 1/\alpha , \\ N^{\ell (2 \alpha - 1)/2}, & \ell > 1/\alpha . \end{array} \right. \end{equation} \tag{ A.32 }$$

In the case α + β > 1, by means of

$$\begin{equation} G_0(x) \le \left\lbrace \begin{array}{@{}ll@{}}x, & 0 \le x \le 1, \\ 1, & x \ge 1,\end{array} \right. \end{equation} \tag{ A.33 }$$

we obtain

$$\begin{equation} J_1 = \left\lbrace \begin{array}{@{}ll@{}}O(N^{(\alpha ,\beta )}), & \alpha > 1/2, \ \beta > 1/2, \\ O(N^{2 \alpha - 1} \ln N), & \alpha > 1/2, \ \beta \le 1/2, \ \ell = 1/\beta , \\ O(N^{2 \alpha - 1}), & \alpha > 1/2, \ \beta \le 1/2, \ \ell \ne 1/\beta , \\ O(N^{2 \beta - 1} \ln N), & \alpha \le 1/2, \ \beta > 1/2, \ \ell = 1/\alpha , \\ O(N^{2 \beta - 1}), & \alpha \le 1/2, \ \beta > 1/2, \ \ell \ne 1/\alpha . \end{array} \right. \end{equation} \tag{ A.34 }$$

Here the symbol N^{(α, β)} is defined in (A.20). Inequality (A.33) similarly gives the estimates

$$\begin{eqnarray} J_2 & = & \left\lbrace \begin{array}{@{}ll@{}}O(N^{1 + \ell (\alpha - 1)}), & \ell < 1/\beta , \\ O(N^{(\alpha + \beta - 1)/\beta } \ln N), & \ell = 1/\beta , \\ O(N^{(\alpha + \beta - 1)/\beta }), & \ell > 1/\beta , \end{array} \right. \nonumber \\ J_3 & = & \left\lbrace \begin{array}{@{}ll@{}}O(N^{1 + \ell (\beta - 1)}), & \ell < 1/\alpha , \\ O(N^{(\alpha + \beta - 1)/\alpha } \ln N), & \ell = 1/\alpha , \\ O(N^{(\alpha + \beta - 1)/\alpha }), & \ell > 1/\alpha . \end{array} \right. \end{eqnarray} \tag{ A.35 }$$

Moreover, it is evident from the inequality G₀(x) ⩽ 1 (x ⩾ 0) that

$$\begin{equation} J_4 = O(1). \end{equation} \tag{ A.36 }$$

Now we can easily see from (A.31), (A.34), (A.35) and (A.36) that

$$\begin{equation} \sum _{j=1}^N \sum _{k=1}^M G_0(p N P_j Q_k)^{\ell } \le \sum _{\nu = 1}^4 J_{\nu } = o(N) \end{equation} \tag{ A.37 }$$

for any ℓ ⩾ 2. This relation results in the asymptotic estimate (A.25).

Let us first discuss the spectral density of the adjacency matrix $\mathcal{A}$. The eigenvalues μ_j, j = 1, 2, ..., N + M of $\mathcal{A}$ consist of M pairs ±ν_j, j = 1, 2, ..., M and N − M zeros. Note that ν²_j > 0 are identified with the eigenvalues of the M × M correlation matrix V with

$$\begin{equation} V_{kl} = \sum _{j=1}^N C_{jk} C_{jl}. \end{equation} \tag{ B.1 }$$

Using the notations

$$\begin{equation} \phi _j = \Phi _j, \ \ \ j = 1,2,\ldots ,N \end{equation} \tag{ B.2 }$$

and

$$\begin{equation} \psi _k = \Phi _{k+N}, \ \ \ k = 1,2,\ldots ,M, \end{equation} \tag{ B.3 }$$

we can rewrite the partition function Z defined in (3.2) as

$$\begin{equation} \fl Z(\mu ) = \int \prod _{j=1}^N {\rm d}\phi _j \int \prod _{k=1}^M {\rm d}\psi _k \ {\rm exp}\left( \frac{{\rm i}}{2} \mu \sum _{j=1}^N \phi _j^2 + \frac{{\rm i}}{2} \mu \sum _{k=1}^M \psi _k^2 - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} \phi _j \psi _k \right). \end{equation} \tag{ B.4 }$$

Then we introduce the replica variables

$$\begin{equation} {\vec{\phi }}_j = \big(\phi _j^{(1)},\phi _j^{(2)},\ldots ,\phi _j^{(n)}\big), \ \ \ {\vec{\psi }}_k = \big(\psi _k^{(1)},\psi _k^{(2)},\ldots ,\psi _k^{(n)}\big) \end{equation} \tag{ B.5 }$$

and

$$\begin{equation} {\rm d}{\vec{\phi }}_j = {\rm d}\phi _j^{(1)}{\rm d}\phi _j^{(2)} \cdots {\rm d}\phi _j^{(n)}, \ \ \ {\rm d}{\vec{\psi }}_k= {\rm d}\psi _k^{(1)}{\rm d}\psi _k^{(2)} \cdots {\rm d}\psi _k^{(n)} \end{equation} \tag{ B.6 }$$

to obtain

$$\begin{eqnarray} &&\fl \langle Z^n \rangle = \int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \int \prod _{k=1}^M {\rm d}{\vec{\psi }}_k \ {\rm exp}\left( \frac{{\rm i}}{2} \mu \sum _{j=1}^N {\vec{\phi }}_j^2 + \frac{{\rm i}}{2} \mu \sum _{k=1}^M {\vec{\psi }}_k^2 \right) \nonumber \\ &&\times \left\langle {\rm exp}\left( -{\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} {\vec{\phi }}_j \cdot {\vec{\psi }}_k \right) \right\rangle . \end{eqnarray} \tag{ B.7 }$$

Now we can see from (2.12) that

$$\begin{equation} \fl\left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} {\vec{\phi }}_j \cdot {\vec{\psi }}_k \right) \right\rangle \sim {\rm exp}\left\lbrace p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k ( {\rm e}^{-{\rm i} {\vec{\phi }}_j \cdot {\vec{\psi }}_k} - 1 ) \right\rbrace . \end{equation} \tag{ B.8 }$$

It should be noted that this asymptotic relation holds if ${\vec{\phi }}_j \cdot {\vec{\psi }}_k$ is in the neighbourhood of the origin. This condition is justified in the limit of large mean degree p → ∞, since ${\vec{\phi }}_j^2$ and ${\vec{\psi }}_k^2$ are scaled as O(p^−1/2) or O(p⁻¹) (see equations (C.4) and (C.23)).

Using the notation

$$\begin{equation} {\tilde{\xi }}_j({\vec{\phi }}) = \delta ({\vec{\phi }}- {\vec{\phi }}_j), \ \ \ {\tilde{\eta }}_k({\vec{\psi }}) = \delta ({\vec{\psi }}- {\vec{\psi }}_k), \end{equation} \tag{ B.9 }$$

we obtain

$$\begin{eqnarray} &&\fl\left\langle {\rm exp}\left( - {\rm i} \sum _{j=1}^N \sum _{k=1}^M C_{jk} {\vec{\phi }}_j \cdot {\vec{\psi }}_k \right) \right\rangle \nonumber\\ &&\sim {\rm exp}\left\lbrace p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k \int {\rm d}{\vec{\phi }} \int {\rm d}{\vec{\psi }} \ {\tilde{\xi }}_j({\vec{\phi }}) {\tilde{\eta }}_k({\vec{\psi }}) ( {\rm e}^{- {\rm i} {\vec{\phi }} \cdot {\vec{\psi }}} - 1 ) \right\rbrace , \end{eqnarray} \tag{ B.10 }$$

so that we find

$$\begin{eqnarray} &&\fl \langle Z^n \rangle \sim \int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \int \prod _{k=1}^M {\rm d}{\vec{\psi }}_k \ {\rm exp}\left\lbrace \frac{{\rm i}}{2} \mu \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \ {\tilde{\xi }}_j({\vec{\phi }}) {\vec{\phi }}^2 + \frac{{\rm i}}{2} \mu \sum _{k=1}^M \int {\rm d}{\vec{\psi }} \ {\tilde{\eta }}_k({\vec{\psi }}) {\vec{\psi }}^2 \right\rbrace \nonumber \\ &&{\hphantom{\fl \langle Z^n \rangle =}} \times \ {\rm exp}\left[ p N \sum _{j=1}^N \sum _{k=1}^M P_j Q_k \int {\rm d}{\vec{\phi }} \int {\rm d}{\vec{\psi }} \ {\tilde{\xi }}_j({\vec{\phi }}) {\tilde{\eta }}_k({\vec{\psi }}) ( { \rm e}^{- {\rm i} {\vec{\phi }} \cdot {\vec{\psi }}} - 1 ) \right] \nonumber \\ &&{\hphantom{\fl \langle Z^n \rangle}} = \int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \int \prod _{k=1}^M {\rm d}{\vec{\psi }}_k \int \prod _{j=1}^N \mathcal{D}\xi _j({\vec{\phi }}) \int \prod _{k=1}^M \mathcal{D}\eta _k({\vec{\psi }}) \nonumber \\ &&{\hphantom{\fl \langle Z^n \rangle =}} \times \prod _{j=1}^N \prod _{\vec{\phi }} \delta (\xi _j({\vec{\phi }}) - {\tilde{\xi }}_j({\vec{\phi }})) \prod _{k=1}^M \prod _{\vec{\psi }} \delta (\eta _k({\vec{\psi }}) - {\tilde{\eta }}_k({\vec{\psi }})) \ {\rm e}^{S_1 + S_2}. \end{eqnarray} \tag{ B.11 }$$

Here, S₁ and S₂ are defined in (3.8) and (3.9). The auxiliary functions $\xi _j({\vec{\phi }})$ and $\eta _k({\vec{\psi }})$ satisfy (3.11). If J is the Laplacian matrix $\mathcal{L}$, we can similarly derive the same formula (B.11) for 〈Zⁿ〉, except the change of S₂ according to (3.10).

It can be readily seen that

$$\begin{eqnarray} &&\fl \int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \prod _{j=1}^N \prod _{\vec{\phi }} \delta (\xi _j({\vec{\phi }}) - {\tilde{\xi }}_j({\vec{\phi }})) \nonumber \\ &&= \int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \int \prod _{j=1}^N \mathcal{D}a_j({\vec{\phi }}) \ {\rm exp} \left[ 2 \pi {\rm i} \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \ a_j({\vec{\phi }}) \lbrace \xi _j({\vec{\phi }})-{\tilde{\xi }}_j({\vec{\phi }}) \rbrace \right] \nonumber \\ &&= \int \prod _{j=1}^N \mathcal{D}a_j({\vec{\phi }}) \ {\rm exp}\left[ \sum _{j=1}^N \left\lbrace 2 \pi {\rm i} \int {\rm d}{\vec{\phi }} \ a_j({\vec{\phi }}) \xi _j({\vec{\phi }}) - W_j \right\rbrace \right], \end{eqnarray} \tag{ B.12 }$$

where

$$\begin{eqnarray} W_j & = & - \ln \int {\rm d}{\vec{\phi }}_j \ {\rm exp}\left\lbrace - 2 \pi {\rm i} \int {\rm d}{\vec{\phi }} \ a_j({\vec{\phi }}) {\tilde{\xi }}_j({\vec{\phi }}) \right\rbrace \nonumber \\ & = & - \ln \int {\rm d}{\vec{\phi }}_j \ {\rm exp}\left\lbrace - 2 \pi {\rm i} a_j({\vec{\phi }}_j) \right\rbrace . \end{eqnarray} \tag{ B.13 }$$

In the limit N → ∞, the dominant contribution comes from the stationary point satisfying

$$\begin{equation} \fl\frac{\delta }{\delta a_j({\vec{\phi }})} \left\lbrace 2 \pi {\rm i} \int {\rm d}{\vec{\phi }} \ a_j({\vec{\phi }}) \xi _j({\vec{\phi }}) - W_j \right\rbrace = 2 \pi {\rm i} \xi _j({\vec{\phi }}) - 2 \pi {\rm i} {\rm e}^{ - 2 \pi {\rm i} a_j({\vec{\phi }}) + W_j } = 0, \end{equation} \tag{ B.14 }$$

which means

$$\begin{equation} - \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) \ln \xi _j({\vec{\phi }}) = 2 \pi {\rm i} \int {\rm d}{\vec{\phi }} \ a_j({\vec{\phi }}) \xi _j({\vec{\phi }}) - W_j. \end{equation} \tag{ B.15 }$$

Therefore, we find an asymptotic estimate

$$\begin{eqnarray} &&\fl\int \prod _{j=1}^N {\rm d}{\vec{\phi }}_j \prod _{j=1}^N \prod _{\vec{\phi }} \delta (\xi _j({\vec{\phi }}) - {\tilde{\xi }}_j({\vec{\phi }})) \sim {\rm exp}\left\lbrace - \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \xi _j({\vec{\phi }}) \ln \xi _j({\vec{\phi }}) \right\rbrace . \end{eqnarray} \tag{ B.16 }$$

One can similarly derive another estimate

$$\begin{eqnarray} &&\fl\int \prod _{k=1}^M {\rm d}{\vec{\psi }}_k \prod _{k=1}^M \prod _{\vec{\psi }} \delta (\eta _k({\vec{\psi }}) - {\tilde{\eta }}_k({\vec{\psi }})) \sim {\rm exp}\left\lbrace - \sum _{k=1}^M \int {\rm d}{\vec{\psi }} \eta _k({\vec{\psi }}) \ln \eta _k({\vec{\psi }}) \right\rbrace \end{eqnarray} \tag{ B.17 }$$

in the limit M → ∞. Then we arrive at

$$\begin{equation} \langle Z^n \rangle \sim \int \prod _{j=1}^N \mathcal{D}\xi _j({\vec{\phi }}) \prod _{k=1}^M \mathcal{D}\eta _k({\vec{\psi }}) {\rm e}^{S_0 + S_1 + S_2}, \end{equation} \tag{ B.18 }$$

where S₀ is defined in (3.7).

Putting the Gaussian ansatz (3.14) into (3.13), we see in the limit n → 0 that

$$\begin{eqnarray} \xi _j({\vec{\phi }}) & = & \Theta _j \,{\rm exp}\left[ \frac{{\rm i}}{2} \mu {\vec{\phi }}^2 + p N P_j \sum _{k=1}^M Q_k \lbrace g_k({\vec{\phi }}) - 1 \rbrace \right], \nonumber \\ \eta _k({\vec{\psi }}) & = & \Omega _k \,{\rm exp}\left[ \frac{{\rm i}}{2} \mu {\vec{\psi }}^2 + p N Q_k \sum _{j=1}^N P_j \big\lbrace h_j({\vec{\psi }})- 1 \big\rbrace \right], \end{eqnarray} \tag{ C.1 }$$

where

$$\begin{equation} g_k({\vec{\phi }}) = \left\lbrace \begin{array}{@{}ll@{}}\displaystyle {\rm exp}\left( - \frac{{\rm i} \tau _k}{2} {\vec{\phi }}^2 \right), & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm adjacency} \ {\rm matrix} \ \mathcal{A}, \\ \displaystyle {\rm exp}\left( - \frac{{\rm i}}{2 (1 - \tau _k)} {\vec{\phi }}^2 \right), & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm Laplacian} \ {\rm matrix} \ \mathcal{L} \end{array} \right. \end{equation} \tag{ C.2 }$$

and

$$\begin{equation} h_j({\vec{\psi }}) = \left\lbrace \begin{array}{@{}ll@{}}\displaystyle {\rm exp}\left( - \frac{{\rm i} \sigma _j}{2} {\vec{\psi }}^2 \right), & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm adjacency} \ {\rm matrix} \ \mathcal{A}, \\ \displaystyle {\rm exp}\left( - \frac{{\rm i}}{2 (1 - \sigma _j)} {\vec{\psi }}^2 \right), & {\rm if} \ J \ {\rm is} \ {\rm the} \ {\rm Laplacian} \ {\rm matrix} \ \mathcal{L}. \end{array} \right. \end{equation} \tag{ C.3 }$$

Let us first consider the adjacency matrix $\mathcal{A}$. We are in a position to take the limit p → ∞ with the scalings

$$\begin{eqnarray} &&\fl\mu = O(p^{1/2}), \ \ \ {\vec{\phi }}^2 = O(p^{-1/2}), \ \ \ \sigma _j = O(p^{-1/2}), \ \ \ {\vec{\psi }}^2 = O(p^{-1/2}), \ \ \ \tau _k = O(p^{-1/2}).\nonumber\\ \end{eqnarray} \tag{ C.4 }$$

Then, we obtain

$$\begin{equation} \mu - \frac{1}{\sigma _j} - p N P_j \sum _{k=1}^M Q_k \tau _k = 0, \ \ \ \mu - \frac{1}{\tau _k} - p N Q_k \sum _{j=1}^N P_j \sigma _j = 0. \end{equation} \tag{ C.5 }$$

The variational equations (3.13) are satisfied by the Gaussian ansatz (3.14), if σ_j and τ_k are determined by these equations.

In order to analytically treat (C.5), we define the scaling variables

$$\begin{equation} E = \mu /\sqrt{p}, \ \ \ s(x) = \sqrt{p} \ \sigma _j, \ \ \ x = j/N, \ \ \ t(y) = \sqrt{p} \ \tau _k, \ \ \ y = k/M. \end{equation} \tag{ C.6 }$$

Then it is straightforward to find

$$\begin{equation} \frac{x^{\alpha }}{s(x)} = E x^{\alpha } - (1 - \alpha ) (1 - \beta ) \int _0^1 y^{-\beta } t(y) \ {\rm d}y \end{equation} \tag{ C.7 }$$

and

$$\begin{equation} \frac{y^{\beta }}{t(y)} = E y^{\beta } - c (1 - \alpha ) (1 - \beta ) \int _0^1 x^{-\alpha } s(x) \ {\rm d}x \end{equation} \tag{ C.8 }$$

in the limit (2.1). Using the notations

$$\begin{equation} \fl S = c (1 - \alpha ) (1 - \beta ) \int _0^1 x^{-\alpha } s(x) \ {\rm d}x, \ \ \ T = (1 - \alpha ) (1 - \beta ) \int _0^1 y^{-\beta } t(y) \ {\rm d}y, \end{equation} \tag{ C.9 }$$

we obtain

$$\begin{equation} \fl\frac{S}{c (1 - \alpha ) (1 - \beta )} = \int _0^1 \frac{1}{E x^{\alpha } - T} \,{\rm d}x, \ \ \ \frac{T}{(1 - \alpha ) (1 - \beta )} = \int _0^1 \frac{1}{E y^{\beta } - S} \,{\rm d}y. \end{equation} \tag{ C.10 }$$

In order to evaluate the behaviour of S and T in the tail region E → ∞, we write

$$\begin{equation} S = E (s^{(R)} - {\rm i} s^{(I)}), \ \ \ T = E (t^{(R)} - {\rm i} t^{(I)}) \end{equation} \tag{ C.11 }$$

with real s^(R), s^(I), t^(R) and t^(I). Then it can be seen that

$$\begin{eqnarray} &&\fl \frac{E (s^{(R)} - {\rm i} s^{(I)})}{c (1 - \alpha ) (1 - \beta )} = \frac{1}{E} \int _0^1 \frac{1}{x^{\alpha } - t^{(R)} + {\rm i} t^{(I)}} \,{\rm d}x = \frac{1}{\alpha E} \int _0^1 \frac{s^{(1 - \alpha )/\alpha }}{s - t^{(R)} + {\rm i} t^{(I)}} \,{\rm d}s \nonumber \\ &&{\hphantom{\fl \frac{E (s^{(R)} - {\rm i} s^{(I)})}{c (1 - \alpha ) (1 - \beta )} }}= \frac{1}{\alpha E} \int _0^1 \frac{s^{(1 - \alpha )/\alpha } \ (s - t^{(R)})}{(s - t^{(R)})^2 +(t^{(I)})^2} \,{\rm d}s - \frac{{\rm i}}{\alpha E} \int _0^1 \frac{s^{(1 - \alpha )/\alpha } \ t^{(I)}}{(s - t^{(R)})^2 +(t^{(I)})^2} \,{\rm d}s.\nonumber\\ \end{eqnarray} \tag{ C.12 }$$

Let us employ an asymptotic formula [6]

$$\begin{equation} \frac{\epsilon }{(u - a)^2 + \epsilon ^2} \sim \pi \delta (u - a), \ \ \ \epsilon \downarrow 0 \end{equation} \tag{ C.13 }$$

and obtain an estimate

$$\begin{eqnarray} \frac{E (s^{(R)} - {\rm i} s^{(I)})}{c (1 - \alpha ) (1 - \beta )} & \sim & \frac{1}{\alpha E} \int _0^1 s^{(1/\alpha )-2} \,{\rm d}s - \frac{{\rm i} \pi }{\alpha E} \int _0^1 s^{(1 - \alpha )/\alpha } \delta (s - t^{(R)}) \,{\rm d}s \nonumber \\ & = & \frac{1}{E (1 - \alpha )} - \frac{{\rm i} \pi }{\alpha E} (t^{(R)})^{(1 - \alpha )/\alpha }, \ \ \ E \rightarrow \infty , \end{eqnarray} \tag{ C.14 }$$

so that

$$\begin{equation} s^{(R)} \sim \frac{c (1 - \beta )}{E^2}, \ \ \ s^{(I)} \sim \frac{c (1 - \alpha )(1 - \beta )\pi }{\alpha E^2} (t^{(R)})^{(1 - \alpha )/\alpha }. \end{equation} \tag{ C.15 }$$

One can similarly derive

$$\begin{eqnarray} &&\fl\frac{E (t^{(R)} - {\rm i} t^{(I)})}{(1 - \alpha ) (1 - \beta )} = \frac{1}{E} \int _0^1 \frac{1}{y^{\beta } - s^{(R)} + {\rm i} s^{(I)}} \,{\rm d}x = \frac{1}{\beta E} \int _0^1 \frac{t^{(1 - \beta )/\beta }}{t - s^{(R)} + {\rm i} s^{(I)}} \,{\rm d}t \nonumber \\ &&{\hphantom{\fl\frac{E (t^{(R)} - {\rm i} t^{(I)})}{(1 - \alpha ) (1 - \beta )}}} \sim \frac{1}{E (1 - \beta )} - \frac{{\rm i} \pi }{\beta E} (s^{(R)})^{(1 - \beta )/\beta }, \ \ \ E \rightarrow \infty , \end{eqnarray} \tag{ C.16 }$$

so that

$$\begin{equation} t^{(R)} \sim \frac{1 - \alpha }{E^2}, \ \ \ t^{(I)} \sim \frac{(1 - \alpha )(1 - \beta )\pi }{\beta E^2} (s^{(R)})^{(1 - \beta )/\beta }. \end{equation} \tag{ C.17 }$$

It follows from (C.15) and (C.17) that

$$\begin{equation} s^{(I)} \sim \frac{1 - \beta }{\alpha } \pi c \left( \frac{1 - \alpha }{E^2} \right)^{1/\alpha }, \ \ \ t^{(I)} \sim \frac{1 - \alpha }{\beta } \pi c^{(1 - \beta )/\beta } \left( \frac{1 - \beta }{E^2} \right)^{1/\beta }. \end{equation} \tag{ C.18 }$$

Now we can evaluate the asymptotic behaviour of the spectral density ρ(μ) in the tail region E → ∞. Equations (3.5) and (3.6) can be utilized as

$$\begin{eqnarray} \rho (\mu ) & = & \lim _{n \rightarrow 0} \frac{2}{(N+M) n \pi } {\rm Im} \frac{\partial }{\partial \mu }(S_0 + S_1 + S_2) \nonumber \\ & = & \lim _{n \rightarrow 0} \frac{1}{(N+M) n \pi } {\rm Re} \left( \sum _{j=1}^N \int {\rm d}{\vec{\phi }} \ \xi _j({\vec{\phi }}) {\vec{\phi }}^2 + \sum _{k=1}^M \int {\rm d}{\vec{\psi }} \ \eta _k({\vec{\psi }}) {\vec{\psi }}^2 \right) \nonumber \\ & = & - \frac{1}{(N+M) \pi } {\rm Im} \left( \sum _{j=1}^N \sigma _j + \sum _{k=1}^M \tau _k \right) \nonumber \\ & \sim & - \frac{1}{\sqrt{p} (1 + c) \pi } {\rm Im} \left(c \int _0^1 s(x) \,{\rm d}x + \int _0^1 t(y) \,{\rm d}y \right) \end{eqnarray} \tag{ C.19 }$$

in the limit (2.1). Here,

$$\begin{eqnarray} \int _0^1 s(x) \,{\rm d}x & = & \int _0^1 \frac{1}{E - T x^{-\alpha }} \,{\rm d}x = \frac{1}{E} + \frac{T}{E} \int _0^1 \frac{1}{E x^{\alpha } - T} \,{\rm d}x \nonumber \\ & = & \frac{1}{E} + \frac{ST}{c E (1 - \alpha ) (1 - \beta )} \end{eqnarray} \tag{ C.20 }$$

and we can similarly obtain

$$\begin{equation} \int _0^1 t(y) \,{\rm d}y = \frac{1}{E} + \frac{ST}{E (1 - \alpha ) (1 - \beta )}. \end{equation} \tag{ C.21 }$$

Then it can be seen from (C.15), (C.17) and (C.18) that

$$\begin{eqnarray} \rho (\mu ) & \sim & - \frac{2}{\sqrt{p} E (1 + c) \pi (1 - \alpha ) (1 - \beta )} {\rm Im}( ST ) \nonumber \\ & = & \frac{2 E}{\sqrt{p}(1 + c) \pi (1 - \alpha ) (1 - \beta )} \lbrace s^{(R)} t^{(I)} + s^{(I)} t^{(R)} \rbrace \nonumber \\ & \sim & \frac{2}{\sqrt{p} (1 + c)} \left\lbrace \frac{c^{1/\beta } (1 - \beta )^{1/\beta }}{\beta E^{(2/\beta )+ 1}} + \frac{c (1 - \alpha )^{1/\alpha }}{\alpha E^{(2/\alpha )+ 1}} \right\rbrace . \end{eqnarray} \tag{ C.22 }$$

This gives the asymptotic spectral density of the adjacency matrix $\mathcal{A}$ in the tail region E → ∞.

We next compute the spectral density of the Laplacian matrix $\mathcal{L}$. Using the scalings

$$\begin{equation} \fl\mu = O(p), \ \ \ {\vec{\phi }}^2 = O(p^{-1}), \ \ \ \sigma _j = O(p^{-1}), \ \ \ {\vec{\psi }}^2 = O(p^{-1}), \ \ \ \tau _k = O(p^{-1}) \end{equation} \tag{ C.23 }$$

and taking the limit p → ∞, we find

$$\begin{equation} \mu - \frac{1}{\sigma _j} - p N P_j = 0, \ \ \ \mu - \frac{1}{\tau _k} - p N Q_k = 0, \end{equation} \tag{ C.24 }$$

so that

$$\begin{eqnarray} {\rm Im}\, \sigma _j & = & {\rm Im}\frac{1}{\mu + {\rm i} \epsilon - p N P_j} = - \pi \delta ( \mu - p N P_j), \nonumber \\ {\rm Im}\, \tau _k & = & {\rm Im}\frac{1}{\mu + {\rm i} \epsilon - p N Q_k} = - \pi \delta ( \mu - p N Q_k), \end{eqnarray} \tag{ C.25 }$$

where is an infinitesimal positive number. Then it follows in the limit (2.1) that

$$\begin{eqnarray} {\rm Im}\sum _{j=1}^N \sigma _j & \sim & - N \pi \int _0^1 {\rm d}x \ \delta (\mu - p (1 - \alpha ) x^{-\alpha }) \nonumber \\ & = & - N \pi \frac{\lbrace p (1 - \alpha )\rbrace ^{1/\alpha }}{\alpha } \frac{1}{\mu ^{(1/\alpha )+1}} \ H\left\lbrace \mu - p (1 - \alpha ) \right\rbrace \end{eqnarray} \tag{ C.26 }$$

and

$$\begin{eqnarray} {\rm Im}\sum _{k=1}^M \tau _k & \sim & - M \pi \int _0^1 {\rm d}y \ \delta (\mu - p c (1 - \beta ) y^{-\beta }) \nonumber \\ & = & - M \pi \frac{\lbrace p c(1 - \beta )\rbrace ^{1/\beta }}{\beta } \frac{1}{\mu ^{(1/\beta )+1}} \ H\left\lbrace \mu - p c (1 - \beta ) \right\rbrace , \end{eqnarray} \tag{ C.27 }$$

where H(x) is defined in (3.18). Therefore, we arrive at

$$\begin{eqnarray} &&\rho (\mu ) = - \frac{1}{(N+M) \pi } {\rm Im} \left( \sum _{j=1}^N \sigma _j + \sum _{k=1}^M \tau _k \right)\nonumber \\ &&\hphantom{\rho (\mu )=} \sim \frac{c \lbrace p (1 - \alpha )\rbrace ^{1/\alpha }}{(1 + c) \alpha } \frac{1}{\mu ^{(1/\alpha )+1}} \ H\left\lbrace \mu - p (1 - \alpha ) \right\rbrace \nonumber \\ &&\hphantom{\rho (\mu )=} + \frac{\lbrace p c (1 - \beta )\rbrace ^{1/\beta }}{(1 + c) \beta } \frac{1}{\mu ^{(1/\beta )+1}} \ H\left\lbrace \mu - p c (1 - \beta ) \right\rbrace . \end{eqnarray} \tag{ C.28 }$$

This gives the asymptotic spectral density of the Laplacian matrix $\mathcal{L}$ in the region μ = O(p).