1 Introduction↩︎

The von Neumann entropy is the canonical measure of entanglement in quantum many-body physics. It admits a well-defined extensive limit in the thermodynamic limit and characterises the quantum correlation structure of the ground state as a bulk thermodynamic quantity [1]–[3]. It quantifies the degree to which one half of a many-body system can be assigned a definite quantum state independent of the other. A vanishing entanglement entropy indicates that the two halves of the system are classically correlated or product-like, while a large value of the entropy indicates that neither half possesses a well-defined state on its own. In this sense the entanglement entropy encodes how thoroughly quantum correlation is spread across in the microscopic degrees of freedom [4], [5].

For quantum systems defined on graph families with large automorphism groups, nontrivial upper bounds on entanglement entropy arise from symmetry and representation-theoretic constraints [6], [7]. Much less is known about the entanglement structure of ground states on large asymmetric graphs, where such constraints are absent. The present work investigates the limit, where asymptotically asymmetric graphs for which such symmetry constraints disappear. Erdős–Rényi random graphs possess a trivial symmetry in the thermodynamic limit [8]–[10], and consequently provide a natural setting to investigate the behaviour entanglement entropy in bulk scale. Such asymmetric interaction geometries arise naturally in complex large scale quantum networks [11]–[13], [13]–[19]. Understanding the entanglement behaviour in this setting therefore offers a basic insight into the generic entangling capacity of complex quantum connectivity beyond the constraints imposed by spatial symmetries.

While the entanglement entropy of free fermions has been extensively studied in translationally invariant systems, disordered lattices, and continuum Fermi gases [20]–[22], comparatively little is known about free fermions on large random graphs. This setting is particularly interesting for the following reason. The entanglement correlation matrix, whose eigenvalues directly determine the von Neumann entropy, inherits its statistical properties entirely from the random graph structure in the thermodynamic limit. This makes the entanglement entropy a direct probe of the spectral and eigenvector statistics of the random Hamiltonian. It raises the central question addressed in this work: to what extent does the ground-state entanglement entropy of free fermions on a large asymptotically asymmetric graphs approach the Page limit, the maximal volume-law entanglement expected for a typical random pure state [23]–[25], and how does the answer depend on the connectivity of the graph?

In this work, we investigate the ground-state entanglement entropy of free fermions on Erdős–Rényi random graphs in the dense regime, where the edge probability \(p\) remains finite as \(N\to\infty\), so that the average degree \(Np\) grows linearly with system size. For the sparse regime (\(p \sim 1/N\)), existence of a extensive entropy density is not guaranteed [26]–[28] and for the threshold regime \((p \sim (\log~N)/N\)), the graph is connected with finite probability but the thermodynamic limit is delicate. The dense regime \(Np\to\infty\) is therefore the natural regime for studying thermodynamic entanglement scaling, as the graph remains connected with high probability while the average degree diverges in the large-\(N\) limit.

2 Model↩︎

We consider a system of \(N\) spinless noninteracting fermionic modes, one attached to each vertex of a graph, which permits each mode to be either empty or occupied. The full many-body Hilbert space has dimension \(2^N\), and grows exponentially with system size. The Hamiltonian in the fermion operators is given by the form \[{h}_{N,p} = \sum_{i,j} (t_{N,p})_{ij}c_i^\dagger c_j. \label{eq:second95quantized}\tag{1}\] The ground state can be represented by the Slater determinant constructed from the single-particle eigenstates of the \(N\times N\) matrix \(t_{N,p}\). The entry \((t_{N,p})_{ij}\) quantifies the hopping amplitude between modes \(i\) and \(j\), that is, a fermion at vertex \(i\) can tunnel to vertex \(j\) if an edge connects them. This reduces the many-body problem, formally posed on a \(2^N\)-dimensional Hilbert space, to the spectral theory of a single \(N\times N\) random matrix. We take this graph to be an Erdős–Rényi random graph on \(N\) vertices, in which each of the \(C(N,2)\) edges is present independently with probability \(p\in(0,1)\). The one-body Hamiltonian is the centred and normalized adjacency matrix \[H_{N,p} = \frac{t_{N,p} - p(\mathbf{1}\mathbf{1}^\top - I_N)}{\sqrt{Np(1-p)}}, \label{eq:H}\tag{2}\] where \(t_{N,p}\) is the \(N\times N\) adjacency matrix of graph, \(I_N\) is the identity, and \(\mathbf{1} = (1,\ldots,1)^\top \in\mathbb{R}^N\), so that \(\mathbf{1}\mathbf{1}^\top - I_N\) is the all-ones matrix with zero diagonal. After this normalisation, the off-diagonal entries have mean zero and variance \(1/N\). For fixed \(p\in(0,1)\), the matrix \(H_{N,p}\) falls within the Wigner universality class. Its empirical spectral measure converges to the semicircle law in the thermodynamic limit, independently of the entry distribution, as for any Wigner-type ensemble with i.i.d. entries of mean zero and finite variance [29], [30]. The transformation from the original Hamiltonian \(t_N\) to the normalized form \(H_{N,p}\) is a rescaling that places the Fermi level at zero and brings the spectrum into the Wigner universality class. As shown in Appendix ¿sec:app:p95independence?, this transformation leaves the entanglement correlation matrix \(C_A\) unchanged in the thermodynamic limit. The bulk eigenvectors of \(t_N\) and \(H_{N,p}\) are identical, and the same set of eigenvectors is occupied at half-filling under both Hamiltonians. Consequently, the entanglement entropy density \(s_\infty\) computed from \(h_{N,p}\) is identical to that computed from \(H_{N,p}\) in the thermodynamic limit.

Let \(\{|\phi_k\rangle\}_{k=1}^N\) denote the orthonormal eigenstates with energy \(\varepsilon_k\) of \(H_{N,p}\). The eigenstates are written in the single-particle basis \(\{|j\rangle\}_{j=1}^N\) of vertex states as \(|\phi_k\rangle = \sum_{j=1}^N \phi_k(j)|j\rangle\). The fermionic creation operator for the eigenmode \(|\phi_k\rangle\) is the corresponding linear combination of site creation operators, \[d_k^\dagger = \sum_{j=1}^N \phi_k(j)\,c_j^\dagger, \label{eq:dk}\tag{3}\] which creates a fermion in the energy eigenstate \(|\phi_k\rangle\). The fermion vacuum \(|\Omega\rangle = \otimes^N_{k=1} |0\rangle_k\) is the state with all \(N\) modes unoccupied, satisfying \(c_j|\Omega\rangle = 0\) for every \(j=1,\ldots,N\). At half-filling with Fermi energy \(\mu_F = 0\), the many-body ground state is the Slater determinant \[|\Psi_0\rangle = \prod_{k:\,\varepsilon_k < 0} d_k^\dagger\,|\Omega\rangle, \label{eq:groundstate}\tag{4}\] where \(\varepsilon_k\) is the eigenvalue of \(H_{N,p}\) associated with the eigenvector \(|\phi_k\rangle\). Since the spectrum of \(H_{N,p}\) is symmetric about zero in the large-\(N\) limit, exactly \(m = N/2\) modes satisfy negative eigenvalues and are occupied with high probability.

For a balanced bipartition of the vertices \(V = A\sqcup B\) with \(|A| = |B| = N/2 \equiv m\), the entanglement between the subsystems \(A\) and \(B\) in the ground state 4 is captured by the projector onto the occupied modes \(\{|\phi_k\rangle : \varepsilon_k<0\}\), that built the ground state \(|\Psi_0\rangle\) out of the vacuum \(|\Omega\rangle\). Since \(|\Psi_0\rangle\) is a Slater determinant, Wick’s theorem guarantees that the reduced density matrix \(\rho_A = \mathrm{Tr}_B|\Psi_0\rangle\langle\Psi_0|\) is entirely determined by the one-body correlation function [20] \[(C_A)_{ij} = \langle\Psi_0|c^\dagger_jc_i|\Psi_0\rangle, \quad i,j \in A, \label{eq:CA95def}\tag{5}\] The matrix element \(\langle\Psi_0|c_j^\dagger c_i|\Psi_0\rangle\) is the projection of \(P_F = \sum_{\varepsilon_k<0}|\phi_k\rangle\langle\phi_k|\) onto sites \(i,j\in A\), written in the site basis. Since \(C_A\) is a compression of the projector \(P_F\), its eigenvalues \(\{\lambda_i\}_{i=1}^m\) lie in \([0,1]\). Physically, \(C_A\) encodes how the \(m\) occupied global modes that make up \(|\Psi_0\rangle\) are redistributed in restricted to subsystem \(A\). Each eigenvalue \(\lambda_i\) is the probability that the corresponding mode is occupied in the reduced state \(\rho_A\), with \(1-\lambda_i\) the probability that it is empty. Because \(|\Psi_0\rangle\) is Gaussian, \(\rho_A\) factorises into independent two-level systems, one per eigenmode of \(C_A\), \(\rho_A = \bigotimes_{i=1}^m \mathrm{diag}(1-\lambda_i,\lambda_i)\), so the entanglement entropy is the sum of the binary Shannon entropies of these \(m\) effective modes [20], [31], \[S = -\sum_{i=1}^{m} \bigl[\lambda_i\log\lambda_i + (1-\lambda_i)\log(1-\lambda_i)\bigr]. \label{eq:entropy}\tag{6}\]

The Fermi projector \(P_F\) can itself be written as a contour integral of the Green function of \(H_{N,p}\). For any Hermitian matrix \(H\) with eigenpairs \(\{\varepsilon_k,|\phi_k\rangle\}\), the Green function \(G(z) = (zI-H)^{-1} = \sum_k |\phi_k\rangle\langle\phi_k|/(z-\varepsilon_k)\) has a simple pole at each eigenvalue, with residue equal to the projector onto the corresponding eigenstate. By the Cauchy residue theorem, integrating \(G(z)\) around any closed contour \(\gamma\) therefore reproduces the sum of exactly those eigen-projectors whose eigenvalues lie inside \(\gamma\). This contour integral of the full \(N\times N\) Green function is the Fermi projector \(P_F\). Since \(C_A = P_F|_A\) is the \((A,A)\) block of \(P_F\), it suffices to know only the \((A,A)\) block of the Green function \((zI_N-H_{N,p})^{-1}\) at each point on \(\gamma\), rather than the full \(N\times N\) matrix.

To compute the correlation matrix \(C_A\) in 5 without diagonalising the full matrix \(H_{N,p}\), we write \(H_{N,p}\) in the block form induced by the bipartition, \[H_{N,p} = \begin{pmatrix} H_{AA} & H_{AB} \\ H_{AB}^\top & H_{BB} \end{pmatrix}, \label{eq:block}\tag{7}\] where \(H_{AA}\) and \(H_{BB}\) are \(m\times m\) real symmetric matrices of intra-subsystem hoppings and \(H_{AB}\) is an \(m\times m\) real matrix of inter-subsystem hoppings. Since \(H_{AA}\), \(H_{BB}\), and \(H_{AB}\) correspond to disjoint edge sets of the graph and the edge indicators are jointly independent, the three blocks are mutually independent random matrices, each with off-diagonal entries of mean zero and variance \(N^{-1}\). The Schur complement identity gives an exact contour integral representation of entanglement correlation matrix [32] (see Appendix A) \[C_A = \frac{1}{2\pi i} \oint_{\gamma} G_{AA}(z) dz, \label{eq:CA95contour}\tag{8}\] where \[G_{AA}(z) = \bigl(zI_m - H_{AA} - \Sigma(z)\bigr)^{-1}, \label{eq:GAA}\tag{9}\] is the Green function of the joint operators \(H_{AA} + \Sigma(z)\) at spectral parameter \(z \in \mathbb{C}\), and

\[\Sigma(z) = H_{AB}\bigl(zI_m - H_{BB}\bigr)^{-1}H_{AB}^\top \label{eq:Sigma}\tag{10}\] is the self-energy from integrating out subsystem \(B\). The contour of the integral \(\gamma\) is the boundary of the rectangle \(\bigl\{z = x+iy : x_L \le x \le x_R,\;-\eta \le y \le \eta \bigr\} \subset \mathbb{C}\), with \(x_L < -2\) chosen to the left of the entire spectrum of \(H_{N,p}\), \(x_R = 0^-\) chosen immediately to the left of the Fermi energy, and \(\eta>0\) a fixed imaginary offset keeping \(\gamma\) off the real axis. The contour is traversed counter-clockwise, so that every eigenvalue \(\varepsilon_k<0\) of \(H_{N,p}\) lies in the interior of \(\gamma\) (see Fig. 1).

The operator \(H_{AA} + \Sigma(z)\) is the effective Hamiltonian of subsystem \(A\) after the \(B\) degrees of freedom have been integrated out, with the self-energy \(\Sigma(z)\) encoding the influence of \(B\) on \(A\) mediated through the inter-subsystem hopping \(H_{AB}\). Eq. 8 is exact for any finite \(N\), and the thermodynamic limit behaviour of \(\Sigma(z)\) is the key object of the analysis that follows.

3 Universality of Entanglement Entropy in the dense regime↩︎

We analyse Eq. 8 in the thermodynamic limit where the edge probability \(p\) is fixed. The blocks \(H_{AA}\) and \(H_{BB}\) are \(m\times m\) Wigner-type matrices with centred entries of variance \(1/N\). Since \(m = N/2\), the effective variance per unit of block size is \(1/2\), and by the Wigner semicircle law [30] their empirical spectral measures converge in the thermodynamic limit to the semicircle supported on \([-\sqrt{2},\sqrt{2}]\). The matrix \(W = H_{AB}H_{AB}^\top\) is an \(m\times m\) positive semidefinite square Wishart-type matrix whose empirical spectral measure converges to the Marchenko–Pastur law with unit ratio supported on \([0,2]\) [33], [34] (see Appendix ¿sec:app:B?).

For that we approximate the Green function using the asymptotic freeness between the blocks \(H_{AA}\) or \(H_{BB}\) and \(H_{AB}\) of the matrix in 7 . By the local semicircle law for \(H_{BB}\) [30], the Green function of \(H_{BB}\) in the thermodynamic limit tends to \[\begin{align} \lim_{N \rightarrow \infty}&&(zI_m - H_{BB})^{-1} = g(z)I_m, \nonumber\\ \text{where} && g(z) = z - \sqrt{z^2 - 2} \label{eq:gz} \end{align}\tag{11}\] is the Stieltjes transform of the semicircle on \([-\sqrt{2},\sqrt{2}]\) (see Appendix ¿sec:app:B?). Since \(H_{AB}\) is independent of \(H_{BB}\), the exact self-energy 10 concentrates to \[\lim_{N\to\infty} \Sigma(z) = g(z)W. \label{eq:Sigma95largeN}\tag{12}\] The asymptotic freeness assumes that the subsystem \(A\) and its environment \(B\) are independent, implying the three blocks \(H_{AA}\), \(H_{BB}\), and \(H_{AB}\) of the Hamiltonian \(H_{N,p}\) correspond to disjoint edge sets of the Erdős–Rényi graph. Therefore, the self-energy \(\Sigma(z) = H_{AB}(zI_m - H_{BB})^{-1}H_{AB}^\top\) involves the product of \(H_{AB}\) and a function of the independent matrix \(H_{BB}\), and its thermodynamic limit is computed using asymptotic freeness between independent Wigner-type random matrices [29], [35]. However, asymptotic freeness captures only the leading-order behaviour of \(\Sigma(z)\) in the thermodynamic limit, with corrections encoded in \(\Delta(z) = \Sigma(z) - g(z)W\).

The magnitude of this correction is governed by the boundary-to-bulk ratio of the graph. In an Erdős–Rényi graph with balanced bipartition, \(|A|=|B|=N/2\), each of the \((N/2)^2\) possible boundary edges between \(A\) and \(B\) is present independently with probability \(p\). It yields an expected number of boundary edges \(|E_{AB}| \sim pN^2/4\). Similarly, each of the \(C(N/2,2)\) possible bulk edges within \(B\) is present independently with probability \(p\), giving \(|E_{BB}| \sim pN^2/8\). The boundary-to-bulk ratio is \(|E_{AB}|/|E_{BB}| = 2\), which is \(\mathcal{O}(1)\) and therefore independent of both \(N\) and \(p\). Unlike geometrically local models in \(d\) spatial dimensions, where the boundary-to-bulk ratio vanishes as \(\mathcal{O}(N^{-1/d})\) [3], the Erdős–Rényi random graph has no geometric locality, and the boundary remains extensive even in the thermodynamic limit. Consequently, the correction term \(\Delta(z)\) remains \(\mathcal{O}(1)\) in the thermodynamic limit and contributes a finite correction to the entropy computed from the approximation of asymptotic freeness.

After assuming this approximation of asymptotic freeness , the Green function 9 therefore simplifies to \[\bar{G}_{AA}(z) = \lim_{N\to\infty} G_{AA}(z) = \bigl(zI_m - H_{AA} - g(z) W\bigr)^{-1}. \label{eq:GAA95limit}\tag{13}\] Equation 13 is therefore the central formula of the dense regime. The expression is exact within the asymptotically free effective theory, in the thermodynamic limit, as it involves only the \(m\times m\) blocks \(H_{AA}\) and \(W\), and it requires no diagonalization of the full Hamiltonian and the limiting Green function \(\bar{G}_{AA}(z)\) is universal. The semicircle law governing \(H_{AA}\) and the Marchenko–Pastur law governing \(W\), as well as the scalar function \(g(z)\), all depend only on the entry variance \(N^{-1}\), and are independent of the edge probability \(p\) and the microscopic distribution of the matrix entries. Consequently, in the thermodynamic limit, \(\bar{G}_{AA}(z)\) depends only on these two universal spectral laws.

Since \(G_{AA}(z)\to\bar{G}_{AA}(z)\) at every \(z\in\gamma\), and the contour \(\gamma\) is fixed deterministic independent of \(N\), \(p\), and the graph realisation, the entanglement correlation matrix concentrates to a deterministic limit, \[\bar{C}_A = \lim_{N\to\infty} C_A = \frac{1}{2\pi i} \oint_{\gamma}\bar{G}_{AA}(z)~dz. \label{eq:CA95limit}\tag{14}\] This follows from the uniform bound \(\|G_{AA}(z)\|\leq 1/\eta\) on \(\gamma\) (where \(\eta>0\) is the imaginary offset of the contour) and the dominated convergence theorem (the thermodynamic limit may be taken inside the contour integral) applied to the contour integral. Since the integrand is universal and the contour is deterministic, the limit \(\bar{C}_A\) is universal, and its empirical spectral measure \(\bar{\rho}_{C_A}\) is a deterministic limit independent of \(p\). The universality of \(C_A\) implies an exact volume law for the entanglement entropy. Defining the empirical spectral measure of \(C_A\) as \(\rho_{C_A} = \frac{1}{m}\sum_{i=1}^m\delta(\lambda-\lambda_i)\), the entropy density 6 , in the thermodynamic limit, can be written exactly as \[\frac{S}{m} = \int_0^1 h(\lambda)\,d\bar{\rho}_{C_A}(\lambda) \equiv s_\infty, \label{eq:s95limit}\tag{15}\] where \(h(x)\) is the binary Shannon entropy function.

For the volume law to hold, \(s_\infty > 0\) is required. If \(s_\infty = 0\), then by the half-filling constraint \(\int\lambda\,d\bar{\rho}_{C_A} = 1/2\), the measure \(\bar{\rho}_{C_A}\) must equal \(\frac{1}{2}\delta(\lambda) + \frac{1}{2}\delta(\lambda-1)\), requiring \(C_A\) to be a projector and the subspace \(A\) to be invariant under \(P_F\). However, by the quantum ergodicity of Wigner matrices [30], the eigenvectors of \(H_{N,p}\) are delocalized with \(|\psi_k(i)|^2\sim 1/N\) for all \(i\), so no eigenvector is supported entirely on \(A\) or \(B\). Therefore \(C_A\) has generic nonzero off-diagonal elements and \(s_\infty > 0\). Therefore the entanglement entropy scales as the volume law in the thermodynamic limit. The entanglement entropy grows linearly with system size \(N\), with a universal coefficient \(s_\infty\) that is independent of \(p\) and the microscopic structure of the graph. This universality reflects the fact that in the thermodynamic limit the entanglement structure is governed entirely by the Wigner semicircle and Marchenko–Pastur spectral laws, which are insensitive to the microscopic details of the connectivity. The volume law arises because the ER random graph has no geometric locality; every vertex is connected to \(\mathcal{O}(Np)\) others throughout the graph, so there is no notion of a boundary whose area could bound the entanglement, in sharp contrast to geometrically local lattice models where an area law holds [3].

Since \(C_A\) is a compression of the Fermi projector onto \(m = N/2\) filled states, the half-filling condition implies \(\frac{1}{m}\sum_{i=1}^m \lambda_i = \frac{1}{m}\operatorname{Tr} C_A = \frac{1}{2}\), so the empirical spectral measure \(\bar{\rho}_{C_A}\) has mean \(1/2\) in the thermodynamic limit. Since the binary entropy function \(h(\lambda)\) is strictly concave on \([0,1]\) with unique maximum \(h(1/2) = \log 2\), by the Jensen inequality, we obtain \[\begin{align} s_\infty =\int_0^1 h(\lambda)d\bar{\rho}_{C_A}(\lambda) \leq h\left(\int_0^1 \lambda d \bar{\rho}_{C_A}(\lambda)\right) = \log 2, \label{eq:Jensen} \end{align}\tag{16}\] with equality if and only if \(\bar{\rho}_{C_A} = \delta(\lambda - \tfrac{1}{2})\).

The maximal entropy density would require all eigenvalues of \(C_A\) to converge to \(1/2\). Since every matrix element of \(C_A\) is a coherent sum over the occupied eigenvectors of the random graph Hamiltonian, the exact saturation of the maximal entropy density would therefore require all these sums to cancel asymptotically for every pair \(i\neq j\). The vanishing of all off-diagonal elements of \(C_A\) simultaneously is impossible in the Erdős–Rényi ensemble, for the following reasons. (a) The ensemble has no geometric symmetry that could enforce systematic cancellation among the coherent sums \((C_A)_{ij} = \sum_{\varepsilon_k<0}\phi_k(i)\phi_k(j)\). This implies that all \(C(N/2,2)\sim N^2/8\) off-diagonal elements vanish simultaneously imposing \(O(N^2)\) independent conditions on the \(N/2\) eigenvectors, which is a massively overdetermined system that generically has no solution. (b) By the quantum ergodicity of the ER random graph, the eigenvectors \(\{\phi_k\}\) are delocalized, with components satisfying \(|\phi_k(i)|^2 \sim 1/N\) for all \(i\) [30]. Each term \(\phi_k(i)\phi_k(j) \sim 1/N\) is therefore nonzero, and the coherent sum over \(N/2\) occupied modes gives \((C_A)_{ij} = \mathcal{O}(1)\), which does not vanish in the thermodynamic limit. Therefore the inequality 16 is strict, and the entropy deficit \[\Delta s = \log 2 - s_\infty = \int_0^1 \bigl[\log 2 - h(\lambda)\bigr] d\bar{\rho}_{C_A}(\lambda) > 0 \label{eq:deficit}\tag{17}\] is strictly positive. The quantity \(\Delta s\) measures the spread of \(\bar{\rho}_{C_A}\) around \(1/2\). It vanishes only in the maximally entangled limit \(\bar{\rho}_{C_A} = \delta(\lambda - \tfrac{1}{2})\), corresponding to the Page value \(s = \log 2\), and it increases as the spectral measure broadens away from \(1/2\). Therefore, the ground state does not achieve the Page value, and the deficit is a universal constant.

Although analytical theory establishes the existence of an extensive entanglement entropy, obtaining the universal coefficient \(s_{\infty}\) analytically requires determining the spectral measure of the restricted correlation matrix \(C_A\). While the asymptotic Green function is obtained implicitly through the Dyson equation, extracting the spectral measure of \(C_A\) entails solving a nonlinear algebraic equation for its Stieltjes transform, for which no closed-form solution is currently known. Consequently, the universal coefficient \(s_{\infty}\) should be determined numerically.

4 Numerical results↩︎

Figure 2 (a) shows the entropy density as a function of system size \(N\) for different values of the edge probability. The entropy density is independent of \(p\) for all system sizes considered, confirming the analytical result. The entropy density converges to a universal constant \(s_\infty \approx 0.386\) nats as \(N\to\infty\), establishing the exact volume law numerically. The finite-\(N\) corrections are visible for small system sizes but become negligible as \(N\) increases.

Figure 2 (b) shows the empirical spectral density \(\rho_{C_A}(\lambda)\) of the entanglement correlation matrix \(C_A\) for \(N = 500\) and \(p = 0.5\), computed from the exact many-body ground state of and from the asymptotic freeness approximation via the contour integral 8 . Both distributions exhibit a U-shaped profile. This U-shape reflects the spread of eigenvalues of \(C_A\) away from the uniform value \(1/2\), which is responsible for the positive entropy deficit. The freeness approximation places slightly less weight near the endpoints than the exact distribution. Hence we see a difference in the entropy density \(s_\infty^{\text{free}} \approx 0.426 > s_\infty \approx 0.386\) is observed.

The discrepancy between the two distributions in Fig. 2 (b) has a physical origin. The asymptotic freeness approximation replaces the exact environment response \((zI_m - H_{BB})^{-1}\) by its isotropic scalar limit \(g(z)I_m\). This approach preserves the average spectral density of subsystem \(B\) via the Stieltjes transform but discards all information encoded in the eigenvectors of \(H_{BB}\) and all higher-order correlations between the environment modes. From the perspective of subsystem \(A\), the environment is replaced by an isotropic effective medium that retains only the coarse-grain spectral information of \(H_{BB}\). This structure breaks the non-homogeneous mixing of the subsystem \(A\) to different modes of \(B\), allowing the concentration of the eigenvalues of \(C_A\) toward \(1/2\) and therefore reducing the entropy below the freeness approximation value.

5 Discussion and Conclusions↩︎

The contour integral representation of the correlation matrix and the Schur complement identity, provides a general framework for understanding the scaling behavior of the ground-state entanglement entropy. The three exact equations 8 –10 hold for any free fermion Hamiltonian with a balanced bipartition, independently of the graph structure. The entanglement entropy satisfies an exact volume law if and only if two conditions hold in the thermodynamic limit: (a) the empirical spectral measure \(\rho_{C_A}\) converges to a deterministic limit \(\bar{\rho}_{C_A}\) in the thermodynamic limit, and (b) the limiting entropy density \(s_\infty = \int_0^1 h(\lambda)\,d\bar{\rho}_{C_A}(\lambda)\) is strictly positive. The first condition holds automatically whenever the empirical spectral measures of the diagonal blocks \(H_{AA}\) and \(H_{BB}\) converge to deterministic limits and the off-diagonal block \(H_{AB}\) has bounded operator norm. The second condition follows from the boundedness of the Green function \(\|G_{AA}(z)\|\leq 1/\eta\) on the contour \(\gamma\). The second condition also requires the eigenvalues of \(C_A\) to be spread strictly across \((0,1)\) rather than concentrated at \(0\) and \(1\) in the thermodynamic limit, which holds if and only if the subspace \(A\) is not invariant under the Fermi projector \(P_F\).

For the Erdős–Rényi graph, both conditions are satisfied rigorously. The first condition follows from the Wigner semicircle law for \(H_{AA}\) and \(H_{BB}\), the independence of \(H_{AB}\) from \(H_{BB}\) by the graph construction, and the bounded operator norm \(\|H_{AB}\|= \mathcal{O}(1)\). The second condition follows from the quantum ergodicity of Wigner matrices [30], so no eigenvector is supported entirely on \(A\) or \(B\). For graph ensembles where one or both conditions fail, such as geometrically local graphs where the boundary-to-bulk ratio vanishes or scale-free graphs where \(\|H_{AB}\|\to\infty\) and the contour integral may break down, the volume law does not hold. A different scaling of entanglement entropy is expected in these cases.

We have studied the ground-state entanglement entropy of free fermions on Erdős–Rényi random graphs in the dense regime, establishing that the entropy obeys a volume law in the thermodynamic limit. We showed that the entropy density is independent of the edge probability and of the microscopic distribution of the matrix entries. The entropy density interpreted as the entanglement per fermion, an intensive thermodynamic quantity, that characterises how much quantum correlation each fermion carries when embedded in a scrambled non-local environment [23], [36]. The strictly positive entropy deficit from the maximum value \(\log 2\) achievable by a single fermionic mode is itself a universal constant, providing a sharp characterisation of the departure of the free-fermion ground state from maximal entanglement.

A natural question is whether this entropy density depends only on the divergence of the average degree, rather than on the specific scaling of the edge probability. This raises the possibility that graph ensembles with \(Np\rightarrow \infty\) but \(p \rightarrow 0\), such as \(p \sim 1/\log~N\), belong to the same entanglement universality class. Investigating this regime and identifying possible crossover behaviour remain important directions for future work.

Another central result of our analysis is that the bipartite entanglement entropy of free fermions on random graphs is not completely determined by the asymptotically free description of the block Hamiltonians. While the asymptotic-freeness approximation accurately captures the average spectral properties of the environment and reproduces the dominant contribution to the entropy, a finite but small discrepancy between the exact and free-theory entropies persists in the thermodynamic limit. This demonstrates that entanglement entropy retains sensitivity to structural information beyond the asymptotically free spectral data. We find therefore entanglement entropy as a non-free-probabilistic observable, even in the thermodynamic limit, for random graphs.

Beyond its mathematical simplicity, the Erdős–Rényi ensemble serves as an idealized model of generic connectivity in disordered quantum systems. Similar random-network descriptions arise naturally in disordered quantum-dot arrays [37] and probabilistic quantum communication networks [13], [38], [39]. Although the connectivity in real systems generally does not follow the Binomial distribution, our approach provides a technique to solve models of asymptotically asymmetric connectivity. Therefore the entropy density constitutes a theoretical prediction for the entanglement per fermion in any system within a universality class, against which experiments and numerical simulations of more complex models can be compared. In particular, the role of interactions, absent in the present free-fermion baseline, in modifying the entropy remains to be studied in strongly correlated systems in random graphs.

APPENDIX A↩︎

APPENDIX B↩︎

Spectral Distributions of \(H_{AA}\) and \(H_{BB}\)↩︎

The blocks \(H_{AA}\) and \(H_{BB}\) are each \((N/2)\times(N/2)\) Wigner-type random matrices with centred independent entries of variance \(1/N\). Since their size is \(M=N/2\) but their entry variance is \(1/N = 2/M\cdot(1/2)\), the effective variance parameter in units of the block size is \(\sigma_*^2 = 1/2\). By the Wigner semicircle law [30], the empirical spectral distributions of \(H_{AA}\) and \(H_{BB}\) each converge to the semicircle law \[\mu(x) = \frac{1}{\pi}\sqrt{2-x^2}, \qquad x\in[-\sqrt2,\sqrt2]. \label{eq:semicircle}\tag{19}\] The corresponding Stieltjes transform is \[m_{AA}(z) = m_{BB}(z) = \int_{-\sqrt2}^{\sqrt2} \frac{\mu(x)}{z-x} dx, \qquad z\in\mathbb{C}\setminus[-\sqrt2,\sqrt2]. \label{eq:m95def}\tag{20}\]

The empirical spectral distributions of \(H_{AA}\) and \(H_{BB}\) converge to the semicircle law \[\mu(x)=\frac{1}{\pi}\sqrt{2-x^2}, \qquad x\in[-\sqrt2,\sqrt2].\]

Introducing the rescaled variables \(x=\sqrt2\,u\) and \(w=z/\sqrt2\), one obtains \[m(z) = \sqrt2\,S(w), \qquad S(w) = \frac{1}{\pi} \int_{-1}^{1} \frac{\sqrt{1-u^2}}{w-u}\,du . \label{eq:S95def}\tag{21}\] The integral \(S(w)\) is the Stieltjes transform of the standard semicircle distribution on \([-1,1]\). Evaluating it, for example by setting \(u=\cos\theta\) and reducing the resulting expression to the classical integral \[\int_0^\pi \frac{d\theta}{w-\cos\theta} = \frac{\pi}{\sqrt{w^2-1}}, \qquad w\notin[-1,1],\] gives \[S(w)=w-\sqrt{w^2-1},\] where the branch is chosen so that \(S(w)\sim (2w)^{-1}\) as \(|w|\to\infty\). Undoing the rescaling yields \[m_{AA}(z)=m_{BB}(z) = z-\sqrt{z^2-2}, \label{eq:mAA95mBB}\tag{22}\] which is the expression used in the main text.

Spectral Distribution of \(W = H_{AB}H_{AB}^{\top}\)↩︎

We define the \((N/2)\times(N/2)\) matrix \[W = H_{AB}H_{AB}^{\top}.\] Since the entries of \(H_{AB}\) have variance \(1/N\), we may write \[H_{AB}=N^{-1/2}\widetilde{X},\] where \(\widetilde{X}\) has independent centred entries of unit variance. It follows that \[W = \frac{1}{N}\widetilde{X}\widetilde{X}^{\top} = \frac{1}{2} Y, \qquad Y = \frac{\widetilde{X}\widetilde{X}^{\top}}{N/2}.\] The matrix \(Y\) is a sample covariance matrix with aspect ratio \(c=1\). By the Marchenko–Pastur theorem [33], its empirical spectral distribution converges to \[\rho_Y(\lambda) = \frac{1}{2\pi} \sqrt{\frac{4-\lambda}{\lambda}}, \qquad \lambda\in(0,4].\]

Since \(W=\frac{1}{2}Y\), the spectral density rescales according to \[\mu_W(x) = 2\,\rho_Y(2x) = \frac{1}{\pi} \sqrt{\frac{2-x}{x}}, \qquad x\in(0,2]. \label{eq:MP95density}\tag{23}\]

The corresponding Stieltjes transform is \[m_W(z) = \int_0^2 \frac{\mu_W(x)}{z-x}\,dx = \frac{1}{\pi} \int_0^2 \frac{1}{z-x} \sqrt{\frac{2-x}{x}}\,dx . \label{eq:mW95def}\tag{24}\] Introducing the substitution \(x=1-\cos\phi\) reduces the integral to \[m_W(z) = \frac{1}{\pi} \int_0^\pi \frac{1+\cos\phi}{z-1+\cos\phi}\,d\phi,\] which may be evaluated using the standard Poisson-kernel integral. The result is \[m_W(z) = 1-\sqrt{\frac{z-2}{z}}, \label{eq:mW}\tag{25}\] where the branch of the square root is chosen such that \(m_W(z)\sim z^{-1}\) as \(|z|\to\infty\), as required for a Stieltjes transform.

APPENDIX C↩︎

Consider the normalized Erdős–Rényi Hamiltonian \[H_{N,p} = \frac{A_N-p(\mathbf{1}\mathbf{1}^\top-I_N)}{\sqrt{Np(1-p)}}. \label{eq:H95scaling}\tag{26}\]

We show that the thermodynamic entanglement entropy density is independent of \(p\in(0,1)\).

The \(i\)-th component of \(t_N\mathbf{1}\) is the degree of vertex \(i\) \[(t_N\mathbf{1})_i = \sum_{j\neq i}A_{ij} = \deg(i). \label{eq:degree}\tag{27}\] In the Erdős–Rényi graph \(G(N,p)\), the degree of each vertex satisfies \(\deg(i)\sim\mathrm{Binomial}(N-1,p)\), with \[\mathbb{E}[\deg(i)] = (N-1)p \approx Np, \qquad \mathrm{Var}[\deg(i)] = (N-1)p(1-p) = O(N).\] By the law of large numbers, \(\deg(i)/N \to p\) almost surely as \(N\to\infty\), so \[(t_N\mathbf{1})_i = Np + O(\sqrt{N}) \quad \text{almost surely, for each } i. \label{eq:degree95concentration}\tag{28}\] Therefore \[t_N\mathbf{1} = Np\cdot\mathbf{1} + O(\sqrt{N})\cdot\mathbf{1} \quad \text{almost surely,} \label{eq:tN1}\tag{29}\] and dividing by \(\sqrt{N}\) \[t_N\frac{\mathbf{1}}{\sqrt{N}} = Np\cdot\frac{\mathbf{1}}{\sqrt{N}} + O(\sqrt{N})\cdot\frac{\mathbf{1}}{\sqrt{N}}. \label{eq:tN95outlier}\tag{30}\] The relative magnitude of the error term is \[\frac{\left\|O(\sqrt{N})\cdot \frac{\mathbf{1}}{\sqrt{N}}\right\|}{Np} = \frac{O(\sqrt{N})\cdot \frac{\sqrt{N}}{\sqrt{N}}}{Np} = \frac{O(\sqrt{N})}{Np} = O\left(\frac{1}{\sqrt{N}}\right) \to 0 \label{eq:relative95error}\tag{31}\] as \(N\to\infty\). Therefore \(\mathbf{1}/\sqrt{N}\) is an asymptotically exact eigenvector of \(t_N\) with eigenvalue \[\lambda_{\mathrm{out}}^{(t)} \approx Np \to \infty \quad \text{as } N\to\infty. \label{eq:outlier95eigenvalue}\tag{32}\]

The bulk spectrum of \(t_N\) lies in the interval \([-2\sqrt{Np(1-p)}, 2\sqrt{Np(1-p)}]\) by the Wigner semicircle law applied to the fluctuation matrix \(t_N - \mathbb{E}[t_N]\) [30]. The spectral gap between the outlier eigenvalue \(\lambda_{\mathrm{out}}^{(t)} \approx Np\) and the top of the bulk spectrum \(2\sqrt{Np(1-p)}\) is \[\lambda_{\mathrm{out}}^{(t)} - 2\sqrt{Np(1-p)} \approx Np - 2\sqrt{Np(1-p)} = O(N) \to \infty. \label{eq:spectral95gap}\tag{33}\] The perturbation of each bulk eigenvector \(\psi_k\) due to the outlier direction \(\mathbf{1}/\sqrt{N}\) is bounded by the ratio of the perturbation norm to the spectral gap \[\left|\left\langle\frac{\mathbf{1}}{\sqrt{N}} \bigg|\psi_k\right\rangle\right| \leq \frac{O(\sqrt{N})}{O(N)} = O\!\left(\frac{1}{\sqrt{N}}\right) \to 0 \label{eq:DavisKahan}\tag{34}\] as \(N\to\infty\). Therefore in the thermodynamic limit \[\left\langle\frac{\mathbf{1}}{\sqrt{N}} \bigg|\psi_k\right\rangle = 0, \qquad k = 1,\ldots,N-1. \label{eq:orthogonality}\tag{35}\] For any bulk eigenvector \(\psi_k\) satisfying 35 \[p\mathbf{1}\mathbf{1}^\top\psi_k = p\mathbf{1}(\mathbf{1}^\top\psi_k) = p\mathbf{1}\cdot\sqrt{N} \left\langle\frac{\mathbf{1}}{\sqrt{N}} \bigg|\psi_k\right\rangle = p\mathbf{1}\cdot\sqrt{N}\cdot 0 = 0. \label{eq:rank195vanishes}\tag{36}\]

The normalized Hamiltonian is \[H_{N,p} = \frac{t_N - p(\mathbf{1}\mathbf{1}^\top - I_N)}{\sqrt{Np(1-p)}} = \frac{t_N - p\mathbf{1}\mathbf{1}^\top + pI_N}{\sqrt{Np(1-p)}}. \label{eq:HNp95decomposed}\tag{37}\] For any bulk eigenvector \(\psi_k\) of \(t_N\) with eigenvalue \(\lambda_k^{(t)}\), using 36 \[(t_N - p\mathbf{1}\mathbf{1}^\top + pI_N)\psi_k = t_N\psi_k - p\mathbf{1}\mathbf{1}^\top\psi_k + p\psi_k \nonumber\\ = \lambda_k^{(t)}\psi_k - 0 + p\psi_k \nonumber\\ = (\lambda_k^{(t)} + p)\psi_k. \label{eq:shifted95eig}\tag{38}\] Applying the scalar scaling \(1/\sqrt{Np(1-p)} > 0\) \[H_{N,p}\psi_k = \frac{\lambda_k^{(t)} + p}{\sqrt{Np(1-p)}} \psi_k \equiv \varepsilon_k^{(H)}\psi_k, \label{eq:H95eig}\tag{39}\] where \[\varepsilon_k^{(H)} = \frac{\lambda_k^{(t)} + p}{\sqrt{Np(1-p)}} = \frac{\lambda_k^{(t)}}{\sqrt{Np(1-p)}} + \frac{p}{\sqrt{Np(1-p)}}. \label{eq:eigenvalue95relation}\tag{40}\] The second term vanishes in the thermodynamic limit \[\frac{p}{\sqrt{Np(1-p)}} = \sqrt{\frac{p}{N(1-p)}} = O\!\left(\frac{1}{\sqrt{N}}\right) \to 0 \quad \text{as } N\to\infty. \label{eq:shift95vanishes}\tag{41}\] Therefore in the thermodynamic limit \[\varepsilon_k^{(H)} \xrightarrow{N\to\infty} \frac{\lambda_k^{(t)}}{\sqrt{Np(1-p)}}, \label{eq:eigenvalue95limit}\tag{42}\] and the bulk eigenvectors of \(t_N\) and \(H_{N,p}\) are identical in the thermodynamic limit. Since the bulk eigenvectors of \(t_N\) and \(H_{N,p}\) are identical and the occupied sets are identical, the Fermi projectors are equal in the thermodynamic limit \[P_F^{(t)} = \sum_{\lambda_k^{(t)} < \mu_{\mathrm{orig}}} |\psi_k\rangle\langle\psi_k| = \sum_{\varepsilon_k^{(H)} < 0} |\psi_k\rangle\langle\psi_k| = P_F^{(H)}. \label{eq:projectors95equal}\tag{43}\] The entanglement correlation matrices are therefore identical \[C_A^{(t)} = P_F^{(t)}\big|_A = P_F^{(H)}\big|_A = C_A^{(H)}. \label{eq:CA95equal}\tag{44}\] This confirms that the \(p\)-independence of \(s_\infty\) is a fundamental property of the Erdős–Rényi random graph ensemble, independent of the normalization convention used for the Hamiltonian.

References↩︎