Limiting spectral distribution of a random matrix with specific structure

Question

First, consider an $N \times N$ Hermitian random matrix $V$ from the Gaussian Unitary Ensemble (GUE). It is well known that the empirical spectral distribution of the GUE satisfies the semicircle law in the large-$N$ limit. Let $U$ be the eigenvector matrix that diagonalizes $V$. The eigenvalue distribution of $T^\dagger T$, where $T^\dagger$ is the complex conjugate of $T$, and $T$ is the $k \times m$ upper-left block of $U$, converges in distribution to the Wachter distribution:

$$ f^{\text{Wachter}}(x) = \frac{1}{2 \pi \lambda} \cdot \frac{\sqrt{(\lambda_+ - x)(x - \lambda_-)}}{x(1 - x)} \cdot \mathbf{I}_{[\lambda_-, \lambda_+]}(x) + \left( 1 - \frac{\kappa}{\lambda} \right) \Theta(\lambda - \kappa) \delta(x), $$

where $ \lambda_- = \left( \sqrt{\kappa(1 - \lambda)} - \sqrt{\lambda(1 - \kappa)} \right)^2 , \lambda_+ = \left( \sqrt{\kappa(1 - \lambda)} + \sqrt{\lambda(1 - \kappa)} \right)^2$ , with $ \lambda = \frac{m}{N}$ and $\kappa = \frac{k}{N} .$

Now, consider the empirical eigenvalue distribution of a matrix $M := WW^\dagger $, where $W$ is the $k \times m$ upper-left block of $X$, and $X$ is the eigenvector matrix that diagonalizes an $N \times N$ adjacency matrix of a $d$-regular graph. In this example, $M$ is here not normalized. $M$ is Hermitian and has real eigenvalues.

For a random d-regular graph, the adjacency matrix is not purely random but has a deterministic structure dictated by the regularity of the graph. This means that the eigenvectors contained in matrix $W$ are not independent and random as in the GUE case, but are influenced by the specific properties of the graph.

For $d \to \infty$, the empirical spectral distribution of $M$ in the large-$N$ limit and $\kappa$, $ \lambda$ $\in (0,1)$ converges in distribution again to the Wachter distribution.

How can I compute the distribution of $M$ for $d = 2$? How do I combinatorially find the $k$-th moment $\beta_k$ of the empirical spectral distribution $F^M$, written as:

$$ \beta_k(M) = \int_{-\infty}^{\infty} x^k F^M(dx) = \frac{1}{n} \text{tr}(M^k), $$

or, how do I find the probability distribution $d\mu$ such that:

$$ \lim_{N \to \infty} E[\text{tr}(M^k_N)] = \int_{-\infty}^{\infty} x^k d\mu(x)? $$

Attached is an example of a numerically generated empirical spectral distribution of the matrix $M$ for the upper-left corner of size $N/2 \times N/2$ for $W$. Is there explicit literature for combinatorics in graph theory related to structured product matrices?

# Load necessary libraries

library(ggplot2)
library(igraph)
library(gridExtra)

# Wachter density function
wachter_density <- function(x, kappa, lambda) {
  lambda_minus <- (sqrt(kappa * (1 - lambda)) - sqrt(lambda * (1 - kappa)))^2
  lambda_plus <- (sqrt(kappa * (1 - lambda)) + sqrt(lambda * (1 - kappa)))^2
  
  # Compute the Wachter density
  if (x >= lambda_minus & x <= lambda_plus) {
    density <- (1 / (2 * pi * lambda)) * sqrt((lambda_plus - x) * (x - lambda_minus)) / (x * (1 - x))
  } else {
    density <- 0
  }
  
  # Dirac-Delta term
  if (x == kappa) {
    density <- density + (1 - kappa / lambda)
  }
  
  return(density)
}

# Define kappa and lambda parameters
kappa <- 0.5
lambda <- 0.5

# Function to simulate the eigenvalue distribution for a regular graph with degree d
simulate_eigenvalue_distribution <- function(N, d, k, m) {
  # Create a d-regular graph with N nodes
  g <- sample_k_regular(n = N, k = d, directed = FALSE, multiple = FALSE)
  
  # Extract adjacency matrix
  A <- as.matrix(as_adjacency_matrix(g))
  
  # Diagonalize the adjacency matrix
  eig <- eigen(A)
  X <- eig$vectors  # Eigenvector matrix
  
  # Extract the k-by-m upper-left block of X
  W <- X[1:k, 1:m]
  
  # Form the matrix M = W^T * W
  M <- t(W) %*% W
  
  # Compute eigenvalues of M
  eigenvalues_M <- eigen(M, symmetric = TRUE, only.values = TRUE)$values
  
  return(eigenvalues_M)
}

# Set parameters for N, k, and m
N <- 1000  # Number of nodes
k <- 500  # Block size k
m <- 500  # Block size m

# Simulate for d = 2 and d = 100
eigenvalues_d2 <- simulate_eigenvalue_distribution(N, d = 2, k, m)
eigenvalues_d100 <- simulate_eigenvalue_distribution(N, d = 100, k, m)

# Create data frames for plotting
df_empirical_d2 <- data.frame(Eigenvalue = eigenvalues_d2)
df_empirical_d100 <- data.frame(Eigenvalue = eigenvalues_d100)

# Compute Wachter density for comparison
x_vals <- seq(0, 1, length.out = 1000)
density_vals <- sapply(x_vals, wachter_density, kappa = kappa, lambda = lambda)
df_wachter <- data.frame(Eigenvalue = x_vals, Density = density_vals)

# Plot for d = 2
p1 <- ggplot(df_empirical_d2, aes(x = Eigenvalue)) +
  geom_histogram(aes(y = ..density..), bins = 50, fill = "blue", alpha = 0.5, color = "black") +
  geom_line(data = df_wachter, aes(x = Eigenvalue, y = Density), color = "red", size = 1) +
  labs(title = "Empirical Density of Eigenvalues (d = 2)",
       x = "Eigenvalue", y = "Density") +
  theme_minimal()

# Plot for d = 100
p2 <- ggplot(df_empirical_d100, aes(x = Eigenvalue)) +
  geom_histogram(aes(y = ..density..), bins = 50, fill = "green", alpha = 0.5, color = "black") +
  geom_line(data = df_wachter, aes(x = Eigenvalue, y = Density), color = "red", size = 1) +
  labs(title = "Empirical Density of Eigenvalues (d = 100)",
       x = "Eigenvalue", y = "Density") +
  theme_minimal()

# Combine the two plots side by side
grid.arrange(p1, p2, ncol = 2)