title: Global eigenvalue fluctuations of random biregular bipartite graphs
authors: Ioana Dumitriu Yizhe Zhu
text: 
  
  
  
  
  
  
  
  
  Department of Mathematics, University of California, San Diego, La Jolla, CA 92093
  
  idumitriu@ucsd.edu
  
  
  
  
  
  Department of Mathematics, University of California Irvine, Irvine, CA 92697
  yizhe.zhu@uci.edu
  This work was partially supported by NSF DMS-1949617
  
    
  [2000]Primary 60C05, 60B20; Secondary 05C50
  
  
  
  Global eigenvalue fluctuations  of random biregular bipartite graphs
      Yizhe Zhu
      December 30, 2023
  ====================================================================
  
  
  
  
  
  We compute the eigenvalue fluctuations of uniformly distributed random biregular bipartite graphs with fixed and growing degrees for a large class of analytic functions. As a key step in the proof, we  obtain a total variation distance bound for the Poisson approximation of the number of cycles and cyclically non-backtracking walks in random biregular bipartite graphs, which might be of independent interest. We also prove a semicircle law for random (d_1,d_2)-biregular bipartite graphs when d_1/d_2→∞.
  As an application, we translate the results to adjacency matrices of uniformly distributed  random regular hypergraphs.
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
   §.§ Eigenvalue fluctuations of random matrices
  
  The study of fluctuations from the limiting empirical spectral distributions (ESDs) for random matrices is a well-established topic of interest in random matrix theory, originated in <cit.>, see also <cit.> and all references therein. More recently, it has been extended to sparse random matrices and random graph-related matrices in various regimes of sparsity and independence (<cit.>), and the natural next target is hypergraphs (<cit.>). 
  
  The ultimate goal in these studies is to see the equivalent of the one-dimensional Central Limit Theorem (CLT) emerge, when examining linear statistics of the spectra of random matrices and random graphs. More precisely, denote by λ_1, …, λ_n the eigenvalues of the random matrix, suitably scaled to put them with high probability on a compact set, and let f be a suitably smooth function.  When the matrices in question are not extremely sparse, one can almost invariably prove that the linear statistic
  
      ℒ(f) = ∑_i=1^n f(λ_i)
  
  has the property that, when centered, it converges to a normal distribution whose variance depends on f:
  
      ℒ(f) - 𝔼 (ℒ(f))  →  N(0, σ_f^2) .
  
  
  
  
    §.§.§ Dense and not-too-sparse Wigner cases.
  
  There is an interesting phenomenon taking place with respect to sparsity; the variance σ_f^2 is the same in the case of Gaussian Orthogonal Ensembles (GOEs) as in the case of the random regular graph under the  permutation model with growing degrees <cit.>:
  
      σ_f^2     =     2 ∑_k=1^∞ ka_k^2 ,
  
  where a_k is the k-th coefficient of f in the Chebyshev polynomial basis expansion. Small variations of this expression also occur in dense Wigner variants and the uniform regular graph model, as follows. For real  Wigner and generalized Wigner matrices in the dense case <cit.>, σ_f also depends on the fourth moments of the off-diagonal entries and the variance of the diagonal entries, which yields corrections to the constants in front of a_1^2, a_2^2 (see Theorem 1.1 in <cit.> for an explicit expression). Similarly, in the uniform regular graph model <cit.>, a correction must be introduced as there are no one- or two-cycles (as the graph is simple), and so the terms corresponding to k=1 and k=2 in the sum (<ref>) are not present.  
  
  However, in the case of sparse Wigner matrices (corresponding to Erdős-Rényi graphs G(n,p) with p→ 0, np→∞, <cit.>), the fluctuations are impacted by the fact that the number of nonzero entries in each row (i.e., the degree of each vertex) fluctuates, and the 4th moment of the scaled adjacency matrix entries grows. The variance σ_f^2 blows up, necessitating another multiplicative scaling of the linear statistic ℒ(f) - 𝔼(ℒ(f)) by √(p), and extracting only part of the expression (<ref>) (see Theorem 1 in <cit.>). 
  
  
  
    §.§.§ Dense Wishart cases.
  
  A similar phenomenon occurs in the Wishart case, i.e., for sample covariance matrices (corresponding to bipartite graphs); in the case of dense matrices with converging aspect ratio, the variance is given in different forms in <cit.>. Although these expressions are not explicit in terms of a Chebyshev polynomial expansion, in <cit.>, it is shown that the covariance between two linear statistics is diagonalized by shifted Chebyshev polynomials.  When the aspect ratio goes to ∞, <cit.> computes the variance which is consistent with the Wigner case in <cit.>. So far, we are not aware of any CLT results for sparse bipartite Erdős-Rényi graphs, but a similar argument as in <cit.> should apply. 
  
  
  For dependent entries (biregular bipartite graphs), we obtain here the variance of the eigenvalue fluctuation in Theorem <ref>, and it matches the one in <cit.>, except for the first coefficient. 
  
  
  
  
    §.§.§ Constant (expected or deterministic) degree.
   
  When p=c/n, the explicit limiting spectral distribution for Erdős-Rényi graphs G(n,p) is not known, although it is known that the measure μ_c exists for every c (given, e.g., by a Stieltjes transform equation as in <cit.>), and if c>1 it consists of a continuous part and an atomic part <cit.>. 
  Convergence of μ_c to the semicircular distribution is studied in  <cit.>, where asymptotic expressions for the moments of μ_c with an o(1/c) term are computed (as c →∞, μ_c converges to the semicircle law).
  
  However, a CLT for Erdős-Rényi graphs G(n,c/n) still holds <cit.> with a more complicated variance that does not follow the same expression as in (<ref>), see Theorem 2.2 in <cit.>. 
  By contrast, in the random d-regular graph case with d finite, the fluctuations are no longer Gaussian. Instead, they are modeled by an infinitely divisible distribution, expressed as a sum of Poisson variables (see <cit.> for the permutation model, and <cit.> for the uniform model). Notably, in the case when the matrix is not symmetric and corresponds to the (directed) cycle structure of a random permutation, <cit.> showed that the global fluctuations could be computed, and whether or not the limiting distribution is Gaussian depends on how smooth the test function is.  For random regular graphs with fixed degree d, we will not see the effect of degree fluctuation. But for Erdős-Rényi graphs G(n,d/n) with expected degree d, the degree fluctuation contributes to the fluctuation of linear statistics, and there is an extra n^-1/2 normalization in the CLT (see <cit.>). Such a difference shows that eigenvalues of random d-regular graphs for fixed d are more rigid than the corresponding Erdős-Rényi graph G(n,d/n). 
  
  For the bipartite Erdős-Rényi case, once again, the limiting distribution is not known, but results that are similar to <cit.> can be found in <cit.>.  We are not aware of any CLT-like results for the fluctuations in this case.
  We compute the fluctuations for the uniformly random biregular bipartite with fixed degrees. Just like in the regular case <cit.>, we see that the fluctuations are modeled by a sum of Poisson variables (Theorem <ref>). 
  
  
  
  
  
  
  Another important  class of random graphs  is the configuration model. The ESD of the configuration model with a large mean degree is not generally given by the  semicircle law <cit.>, but no linear statistics result is known. When the degrees grow with the size of the graph, similar to Erdős-Rényi graphs and random regular graphs, we expect a central limit theorem for linear statistics holds. There are linear statistics results for inhomogeneous matrix models beyond Wigner matrices <cit.>, and it might be possible to apply their techniques, together with a coupling argument introduced in <cit.>, to study the linear statistics for the configuration model. 
  
  
  
  
  
  
   §.§ Random biregular bipartite graphs
  
  Biregular bipartite graphs have found applications in  error correcting codes, matrix completion, and community detection, see for example <cit.>.
  An (n,m,d_1,d_2)-biregular bipartite  graph is a bipartite graph G=(V_1,V_2, E) where |V_1|=n, |V_2|=m and every vertex in V_1 has degree d_1 and every vertex in V_2 has degree d_2. Here we must have nd_1=md_2= |E|. When the number of vertices is clear, we call it a (d_1,d_2)-biregular bipartite graph for simplicity. 
  
  
  
  Let X∈{0,1}^n× m be a matrix indexed by V_1× V_2 such that 
  X_ij=1 if and only if (i,j)∈ E. The adjacency matrix of a (d_1,d_2)-biregular bipartite graph with V_1=[n], V_2=[m] can be written as 
  
      A=[   0   X; X^⊤   0 ].
  
  All eigenvalues of A come in pairs as {-λ, λ}, where |λ| is a singular value of X, along with extra |n-m| zero eigenvalues. It's easy to see  λ_1(A)=-λ_n+m(A)=√(d_1d_2). 
  
  The empirical spectral distribution for uniformly distributed random biregular bipartite graphs (RBBGs), which is the equivalent of the Kesten-McKay law, was first computed in <cit.> for the case of fixed d_1, d_2. For growing degrees, when d_1/d_2 converges to a positive constant, the analogue to the Marčenko-Pastur law was proved in <cit.>.  
  
  
  In this paper, instead of examining the spectrum of A, we will be looking at the spectrum of the matrix XX^⊤ - d_1I. This serves two purposes: one, it allows for an immediate parallel to the sample covariance matrix (Wishart) case, and two, it allows us to deal with all regimes in a unitary fashion. The eigenvalues of XX^⊤-d_1I are the shifted squares of the eigenvalues of A. Any result on global fluctuations for linear statistics of the spectrum of XX^⊤-d_1I is automatically converted into an equivalent result for the spectrum of A. However, because any result of fluctuations must necessarily put most of the eigenvalues (with the exception of the deterministic outliers) on a compact interval, scaling must be involved. This works perfectly fine when the ratio d_1/d_2 is bounded, but it becomes tricky when it is not, and the matrix XX^⊤-d_1I allows us to do the scaling in a more natural way, similarly to the sample covariance (Wishart) matrix with unbounded aspect ratio in <cit.>. 
  
  
  To prove a result on eigenvalue fluctuations, we need two special ingredients: eigenvalue confinement on a compact interval and asymptotic behavior of cycle counts. For the former, we make use of the spectral gap shown in <cit.> for the fixed degree case and <cit.> for the growing degree case. Previous results  of this kind were obtained for random regular graphs <cit.> for a fixed degree, and <cit.> for growing degrees.
  
  For the latter, we use Stein's method to approximate cycle counts as Poisson random variables by bounding the total variation distance (Theorem <ref>) and obtain a Poisson approximation of the number of cyclically non-backtracking walks (Corollary <ref>). 
  Note that computing cycle counts is a fundamental problem in the study of random graphs, ever since the seminal papers of <cit.> and more general <cit.>. 
  
  To prove our results, we follow the recipe of <cit.>  by using switching to construct exchangeable pairs of graphs that allow us to estimate cycle counts. The switching we use here differs from <cit.> and is suitable for  biregular bipartite graphs. In the analysis of switchings, a new challenge is an imbalance between the parameters d_1,d_2 when the aspect ratio is unbounded. Our  results on cycle counts hold for a large range of d_1,d_2, and are notably independent of the aspect ratio as long as the cycle length is small. It is also worth noting that the method of switching has been applied to other problems on random  biregular  bipartite graphs, for example,  <cit.>.
  
  Finally, we also obtain an algebraic relation between linear eigenvalue statistics on modified Chebyshev polynomials and cyclically non-backtracking walks (Theorem <ref>). Then based on the spectral gap results in <cit.> and approximation theory for Chebyshev polynomials <cit.>, we extend the eigenvalue fluctuation results to a general class of analytic functions.
  
  
  
  
  
   §.§ Main results
   Our main contributions are represented by Theorems <ref>, <ref>, establishing the behavior of the global fluctuations for the linear statistics of eigenvalues of RBBGs in the fixed d_1, d_2, respectively, in the d_1· d_2 →∞ cases.  Note that Theorem <ref> describes the behavior of the fluctuations even in the case when the limiting ESD does not exist since it merely requires d_1/d_2 to be bounded, rather than to converge to a number in [1, ∞) (which would be the necessary condition for the ESD to converge). In 
  addition, we  show that the covariance between two linear statistics with different test functions is given by the coefficients in their Chebyshev expansions.
  
  As part of the proofs for our main results, we also describe the asymptotic behavior of the cycle counts (Theorem <ref>). Based on the cycle counts estimates, we then use the locally tree-like structure of RBBGs to prove a global semicircle law in the case when the degree goes slowly (d_1=n^o(1)) and  d_1/d_2 is unbounded (Theorem <ref>).
  
  Finally, as an important application, we obtain  equivalent results for uniformly distributed random regular  hypergraphs, including cycle counts, global laws, spectral gaps, and eigenvalue fluctuations.
  
  
  
  
   §.§ Organization of the paper
  
  In Section <ref> we prove our results on cycle counts in random biregular bipartite graphs. Section <ref> collects relevant results for the spectral gap and eigenvalue confinement on a compact interval from the literature. Section <ref> proves our main results, Theorems <ref> and <ref>. Section <ref> proves a global semicircle law for RBBGs when d_1/d_2 is unbounded. In  Section <ref>, we use the connections established in <cit.> to prove several results on uniformly distributed regular hypergraphs. 
  
  
  
  
  
  
  § CYCLE COUNTS
  
  
  
   §.§ Counting switchings
  
  In this section, we estimate the number of switchings that create or delete a cycle in a biregular bipartite graph. The precise definitions of switchings for our purposes are given in Definition <ref> and Definition <ref>.  These estimates will be used in Section <ref> to show that cycle counts converge in distribution to Poisson random variables.
  
  
  Throughout the paper, when we say a  cycle, we mean a simple cycle, i.e., all vertices in a cycle are distinct.
  
  
  Let K_n,m be the complete bipartite graph on n+m vertices with V_1=[n], V_2=[m].
  Let H⊆ K_n,m be a subgraph with v vertices. For any i∈ K_n,m, let g_i,h_i denote the degree of i considered as a vertex in a biregular bipartite graph G=(V_1,V_2, E) and the subgraph H, respectively. Let h_max be the largest value of h_i and  |H| be the number of edges of H. Denote by 
  
      [x]_a=x(x-1)⋯(x-a+1)
   the falling factorial.  The following estimate is given in   <cit.>.
   
  	Assume d_1≥ d_2 and nd_1≥ 2d_1(d_1+h_max-2)+|H|+1. Then
  	
      ℙ(H⊆ G)≤∏_i=1^v [g_i]_h_i/[nd_1-4d_1^2-1]_|H|.
  
  
  We first prove several estimates on random biregular bipartite graphs based on Proposition <ref>. 
  
  
  
  Let G be a random (d_1,d_2)-biregular bipartite graph with d_2≤ d_1≤ n^1/3.
  
  	
    * Suppose H is a subgraph of the complete graph K_n,m in which every vertex has degree at least 2. Let e be the number of edges  in H. Suppose e=o(n^1/3). Then
  	
      ℙ(H⊆ G)≤ c_1((d_1-1)(d_2-1)/nm)^e/2.
  
  
    * Let α be a cycle of length 2k in the complete bipartite graph K_n,m. Suppose   k≤ n^1/10, then
   
      ℙ(α⊆ G)≤ c_1((d_1-1)(d_2-1)/nm)^k.
   
   
    * Let β be another cycle of length 2j≤ 2n^1/10 in the complete bipartite graph K_n,m. Suppose α,β share f edges. Then
   
      ℙ(α∪β⊆ G)≤ c_1((d_1-1)(d_2-1)/nm)^j+k-f/2.
   
  
  
  
  It suffices to prove (<ref>). Then (<ref>) and (<ref>) follow as special cases. Since H has e edges, and H is bipartite, it satisfies 
      ∑_i∈ V_1h_i=∑_i∈ V_2h_i=e.
  
  Since h_i≥ 2 for all i∈ V(H), we know [g_i]_h_i≤ (g_i(g_i-1))^h_i/2. Therefore from Proposition <ref>,
  
      ℙ(H⊆ G)   ≤(d_1(d_1-1))^e/2(d_2(d_2-1))^e/2/[nd_1-4d_1^2-1]_e
         =((d_1-1)(d_2-1)/nm)^e/2(nd_1)^e/[nd_1-4d_1^2-1]_e.
  
  
  Recall d_1≤ n^1/3,e=o(n^1/3), and (1+x)^r=1+O(rx) if rx→ 0. We have for some absolute constant c_1>0,
  
      (nd_1)^e/[nd_1-4d_1^2-1]_e≤(nd_1/nd_1-4d_1^2-e )^e=( 1+4d_1^2+e/nd_1-4d_1^2-e)^e≤ c_1.
  
   This proves (<ref>).
  
  
   Let G be a (d_1,d_2)-biregular bipartite graph.  Let C_j be the number of cycles of length 2j in G. We will always represent a cycle by a vertex sequence starting from a vertex in V_1.  Suppose α=(x_1,y_1,⋯, x_k,y_k) is a cycle of length 2k in G with x_i∈ V_1, y_i∈ V_2, 1≤ i≤ k, where   y_k is connected to x_1 in the cycle α. 
   
   
    Let e_i=u_iv_i, e_i'=u_i'v_i'  be the  edges with  with u_i, u_i'∈ V_1,v_i, v_i'∈ V_2, 1≤ i≤ k  such that neither u_i,u_i' is adjacent to y_i  for 1≤ i≤ k and neither v_i, v_i' is adjacent to x_i. See the left part of Figure <ref> for an example.
   
   
   We now introduce our definitions of switching for biregular bipartite graphs.
   
    Consider the action of deleting all  4k  edges  ẽ_i, 1≤ i≤ 2k and  e_i,e_i', 1≤ i≤ k, and replacing them by the edges x_iv_i, x_iv_i', y_iu_i, y_iu_i' for 1≤ i≤ k.   We obtain a new biregular bipartite graph G' with the cycle α deleted. We call this action induced by the 6 sequences (x_i),(y_i), (u_i),(u_i'), (v_i),(v_i'),  1≤ i≤ k a forward α-switching.
    See Figure <ref> for an example. We will consider forward α-switchings
  only up to cyclic rotation and inversion of indices in [k]; that is, we identify the 2k different forward α-switchings obtained by applying the same cyclic rotation or inversion on [k] to the 6 sequences (x_i),(y_i), (u_i),(u_i'), (v_i),(v_i'), 1≤ i≤ k. 
   
   
    
  
  
  Suppose G contains  paths v_ix_iv_i' and u_iy_iu_i' for 1≤ i≤ k, where x_i,u_i, u_i'∈ V_1,  y_i,v_i,v_i'∈ V_2. Consider  deleting all 4k  edges v_ix_i, v_i'x_i, u_iy_i, u_i'y_i  for 1≤ i≤ k,  and replacing them with u_iv_i, u_i'v_i', x_iy_i, y_ix_i+1 for 1≤ i≤ k. We obtain a new graph G' with a cycle α=(x_1,y_1,⋯, x_k,y_k). Such action is called a backward α-switching induced by the sequences (x_i),(y_i), (u_i),(u_i'), (v_i),(v_i'), 1≤ i≤ k. We also identify the 2k different backward α-switchings obtained by applying the same cyclic rotation or inversion on the index set [k].
  
   
   
   Let r be an integer; we say that a cycle is short if its length is less than or equal to 2r. 
  
   We call  a α-switching valid if α is the only short cycle created or destroyed by the switching. For each forward α-switching from G to G', there is a corresponding backward α-switching from G' to G by simply reversing the operation (i.e. from right to left  in Figure <ref>). 
   
   
   Let F_α  be the number of all valid forward α-switchings from G to some G' and let B_α be the number of all valid backward α-switchings from some G' to G. In the following two lemmas, we estimate F_α and  B_α for biregular bipartite graphs.
  
  
   
   
  Let G be a deterministic (d_1,d_2)-biregular bipartite graph with d_1≥ d_2   and  cycle counts C_k, 2≤ k≤ r. For any short cycle α⊆ G of length 2k, we have
  
      F_α≤ [n]_k[m]_kd_1^kd_2^k.
  
  If α does not share an edge with another short cycle, then for an absolute constant c_1>0, we have
  
      F_α≥ [n]_k[m]_k	d_1^kd_2^k ( 1-4k∑_j=2^rj C_j+c_1k(d_1-1)^r(d_2-1)^r/nd_1).
  	
   
  
  
  
  
  
  Consider a cycle denoted by α=(x_1,y_1,⋯, x_k,y_k). Denote  edges 
   
      ẽ_i=x_iy_i,    ẽ_i+k=y_ix_i+1, 1≤ i≤ k,
   
   where x_k+1:=x_1. There are at most [n]_kd_1^k[m]_kd_2^k
  many ways to choose  edges e_i=u_iv_i and e_i'=u_i'v_i' for 1≤ i≤ k, which gives the upper bound (<ref>).  For the  k edges e_i, 1≤ i≤ k, we require distinct u_i∈ V_1, 1≤ i≤ k, and we have d_1 choices for each v_i, given the degree constraint on u_i. This gives [n]_kd_1^k many choices altogether. For the remaining  edges e_i', 1≤ i≤ k we require distinct v_i'∈ V_2, 1≤ i≤ k and each for each u_i' we have d_2 choices, giving us a factor of [m]_k d_2^k.  Therefore (<ref>) holds.
  
   For the rest of the proof,  we always use the same way to count α-switchings by counting the choices of  e_i,e_i'. We use the parameter d_1 to control the choices from e_i, 1≤ i≤ k and the parameter d_2 for the choices from e_i',1≤ i≤ k.
  
   
  
  
  To prove the lower bound in (<ref>), we choose a subset of configurations that are guaranteed to have a valid forward α-switching. Consider e_i, e_i', 1≤ i≤ k such that the following holds: 
  
  	
    * e_i and e_i' are not contained in any short cycle in G for 1≤ i≤ k.
  
    * The distance from any vertex in  {e_i,e_i'}  to any vertex in ẽ_i is at least 2r for any 1≤ i≤ k.
  	
    * The distance between any two different edges among the 2k edges {e_i,e_i', 1≤ i≤ k} is at least r. 
  	
    * For all 1≤ i≤ k,
  	 the distance between  v_i and v_i' is at least 2r, and the distance between  u_i and u_i' is at least 2r.
  	 
  
  Recall the definition of ẽ_i in (<ref>). By Condition (<ref>), for all 1≤  i≤ k, u_i,u_i' are not adjacent to y_i, also v_i,v_i' are not adjacent to x_i, which satisfies the definition of a forward α-switching. Let G' be the  graph obtained  by applying the forward α-switching from G. We need to check that α is the only cycle deleted in G by this switching, and no other short cycles are created in G'. 
  
  Since α shares no edges with other short cycles by our assumptions, deleting α will not destroy other short cycles. From Condition (<ref>), deleting e_i,e_i' will not destroy any short cycles either.  
  
  Next, we show no other short cycles are created in G'. Suppose there exists a new short cycle β in G' created by the switching. Then β contains paths in G∩ G' separated by edges created in the forward switching   in G' (β must contain at least such edge because it is created). Any such path in G∩ G' must have a length at least r, because
  
  	
    * if it starts and ends at vertices in  α and has length less than r, then combining this path with a path in α gives a short cycle in G that intersects α, which is a contradiction to our assumption on α;
  	
    * if it starts in α and ends in {u_i,v_i, u_i', v_i'} for some i and has length less than r, then combining this path with a path in α gives a path between ẽ_̃ĩ, e_i or between ẽ_̃ĩ, e_i' of length less than 2r, which violates Condition (<ref>);
  	
    * if it starts and ends in different edges among {e_i, e_i',1≤ i≤ k}, then it must have length at least r by Condition (<ref>); 
  	
    * if it starts at some vertex in e_i and ends at some vertex in e_i, then the path must start and end at different vertices in e_i. Otherwise,  β is not a cycle in the sense of Definition <ref>. Then the path combined with e_i is a cycle.  By Condition (<ref>),  it has a length at least r, a contradiction. In the same way, it cannot start at some vertex in e_i' and end at some vertex in e_i'.
  
  This implies β contains exactly one path in G∩ G'. If not, the two separated paths together with new edges in G' have lengths greater than 2r, a contradiction to the condition that β is a short cycle.
  Given the path in G∩ G', the remainder of β has two cases: 
  
  	
    * a single edge that can be x_iv_i, x_iv_i', y_iu_i, or y_iu_i'  for some 1≤ i≤ k, then by Condition (<ref>), the path in G∩ G' connecting the two vertices in the edge has length at least 2r, which is a contradiction to the fact that β is a short cycle; 
  	
    * a single path v_ix_iv_i' or u_iy_iu_i', which is impossible by Condition (<ref>).
  
  
  From the analysis above, no such β can exist, hence any α-switching satisfying Conditions (<ref>)-(<ref>) is valid. 
  
  Next, we find the number of all switchings  satisfying Conditions (<ref>) to (<ref>) to have a lower bound on F_α. We will do this by bounding from above the number of switchings out of the [n]_k[m]_kd_1^kd_2^k many choices counted in (<ref>) that fail one of the Conditions (<ref>)-(<ref>). We treat the 4 conditions in the following (a)-(d) parts. 
  
   (a) There are a total of at most ∑_i=2^r 2jC_j edges in all short cycles of G.  For some 1≤ i≤ k, if we choose  one edge e_i  from a short cycle  and the other (2k-1)    edges arbitrarily,  we obtain a forward α-switching that fails Condition (<ref>). The number of all possible choices is at most
   
      k∑_j=2^r2jC_j · 	[n-1]_k-1[m]_kd_1^k-1d_2^k.
  
   And if we choose e_i' from a short cycle and the other (2k-1)  edges arbitrarily, the number of all possible choices is at most
   
      k∑_j=2^r2jC_j · 	[n]_k[m-1]_k-1d_1^kd_2^k-1.
  
   Altogether the number of  choices is at most 
  	
      4/nd_1[n]_k[m]_k(d_1d_2)^kk∑_j=2^rjC_j.
    
   
  
  
  	(b) To fail Condition (<ref>), we can obtain α-forward switchings by choosing   (2k-1) edges arbitrarily, and then  choose one  edge e_i or e_i' that is at most of distance 2r-1 from ẽ_i for some 1≤ i≤ k. From the degree constraints, the number of edges of distance less than 2r from some edge is at most O((d_1-1)^r(d_2-1)^r).  Similar to Part (a), by considering whether e_i or e_i' is chosen for 1≤ i≤ k,
  	the number of such switchings  is at most 
  	
      ([n-1]_k-1[m]_kd_1^k-1d_2^k +[n]_k[m-1]_kd_1^kd_2^k-1)· k ·  O((d_1-1)^r(d_2-1)^r) 
      
      =   1/nd_1	[n]_k[m]_k(d_1d_2)^k k O((d_1-1)^r(d_2-1)^r).
  
  	
  	(c)  For Condition (<ref>), there are three cases to consider depending on  whether the pair is (e_i,e_j), (e_i',e_j') or (e_i,e_j').
  	
  	Suppose the pair e_i,e_j violates Condition (<ref>). We pick the pair of edges  that are within distance r-1 and pick the remaining (2k-2) edges arbitrarily. There are (nd_1) many ways to choose e_i. When e_i is fixed,  there are at most 
  	O((d_1-1)^(r+1)/2 (d_2-1)^(r+1)/2) choices for e_j. Hence the number of switchings that fail Condition (<ref>) is at most 
    
      (nd_1)· [    O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2)] · k(k-1)	· ([n-2]_k-2 [m]_k d_1^k-2 d_2^k)
      
        =   1/nd_1[n]_k[m]_kd_1^kd_2^k ·  k^2·  O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2).
  
     By the same argument, if the pair is (e_i',e_j'), the number of switchings that fail Condition (<ref>) is at most
   
      1/nd_1[n]_k[m]_kd_1^kd_2^k k^2· O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2).
  
    When the two edges of the pair violating Condition (<ref>) are e_i,e_j' for some i,j, the number is at most  
  	
      (nd_1)· [  O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2)]· (2k^2)· ([n-1]_k-1[m-1]_k-1d_1^k-1 d_2^k-1) 
      
      	=   1/nd_1[n]_k[m]_kd_1^kd_2^k k^2O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2).
  
  
  
  	Combining the three cases in Part (c),  the number of switchings that violate Condition (<ref>) is at most
  	
      1/nd_1[n]_k[m]_k(d_1d_2)^kk^2· O((d_1-1)^(r+1)/2(d_2-1)^(r+1)/2).
  
   
    
  	(d)  Since the distance between a pair of vertices in V_1 or V_2 must be even, to violate Condition (<ref>), we can choose a pair u_i,u_i'∈ V_1 or v_i,v_i'∈ V_2 that are within distance 2r-2 first, then choose other edges arbitrarily. Similar to the cases above, the number of  switchings that fail  Condition (<ref>) is  at most   
  	
      1/nd_1[n]_k[m]_kd_1^kd_2^kkO((d_1-1)^r(d_2-1)^r).
  
  
  
  Combining the 4 Cases (a)-(d) above, from (<ref>), (<ref>),(<ref>), and (<ref>), we have at most
  
      4/nd_1[n]_k[m]_k (d_1d_2)^k(k∑_j=2^r j C_j+ O(k(d_1-1)^r(d_2-1)^r)	)
  
  many switchings that fail one of the Conditions (<ref>)-(<ref>) among the  [n]_k[m]_k(d_1d_2)^k possible switchings. Then for an absolute constant c_1>0, 
  
      F_α≥ [n]_k[m]_kd_1^kd_2^k(1-4k∑_j=2^r jC_j+c_1k(d_1-1)^r(d_2-1)^r/nd_1).
  
  Therefore (<ref>) holds.
  
  
  For the number of backward switchings, we obtain a similar upper bound, but the lower bound is only in expectation.
  
  Let G be a random (d_1,d_2)-biregular bipartite graph	and let α be a cycle of length 2k≤ 2r in the complete bipartite graph K_n,m. Let B_α be the number of valid backward switchings from G that create α. Then 
  
      B_α≤ (d_1(d_1-1))^k(d_2(d_2-1))^k,
  
  and there is an absolute constant c_2>0 such that
  
      𝔼 B_α≥ 	(d_1(d_1-1))^k(d_2(d_2-1))^k(1- c_2k(d_1-1)^r(d_2-1)^r/nd_1).
  
  
   
  
  
  	Given α, from the degree constraints, the number of choices for u_i, u_i',v_i, v_i', 1≤ i≤ k that yield a valid backward α switching is at most (d_1(d_1-1))^k(d_2(d_2-1))^k, which gives (<ref>).
  	
  
  	 For  the lower bound, we consider the quantity B:=∑_β B_β, where β is summing over all possible  cycles of length 2k in the complete bipartite graph K_n,m. 
  	As in the proof of Lemma <ref>, we give conditions that guarantee a valid backward switching.
  	
  	Assume β=(x_1,y_1,⋯, x_k,y_k). We first consider backward switchings that create β.	 
  		Suppose the paths v_ix_iv_i', u_iy_iu_i', 1≤ i≤ k in G satisfy the following conditions:
  	
  		
    * The edges x_iv_i, x_iv_i', y_iu_i, and y_iu_i' are not contained in any short cycles.
  		
    * For 1≤ i≤ k, the distance between any vertex in the path v_ix_iv_i' and any vertex in the path u_iy_iu_i' is at least 2r.
  		
    * For all 1≤ i≤ k and 1≤ j≤ k/2, the distance between the paths v_ix_iv_i' and v_i+jx_i+jv_i+j' (the index i+j is calculated modulo k) and  the distance between  u_iy_iu_i' and u_i+jy_i+ju_i+j' are at least 2r-2j+1. 
  		
    * For 1≤ i≤ k, 1≤ j≤ k/2, the distance between v_ix_iv_i' and u_i+jy_i+ju_i+j', and the distance between u_iy_iu_i' and v_i+jx_i+jv_i+j' are at least 2r-2j+2.  
  	 	
  We will show the four conditions above guarantee a valid backward β-switching.
  
  By Condition (<ref>), no short cycles are deleted. We denote x≁y if two vertices x,y are not connected in G.
  An immediate consequence of Condition (<ref>) ensures that x_i≁y_i and u_i≁v_i, u_i'≁v_i', and Condition (<ref>) ensures that y_i≁x_i+1. Therefore such  switching can be applied. 
  
  Let G' be the graph obtained by applying the backward β-switching. We need to check that no short cycles other than β are created in G'.  
  
  
  Suppose a short cycle β'≠β is created. Then β' possibly consists of paths in G∩ G', portions of β, and edges u_iv_i,u_i'v_i' for some 1≤ i≤ k. 
  Any such path in G∩ G' must have length at least r because
  	
  		
    * if it starts in one of the  sets {x_i,v_i,v_i'} or {y_i, u_i,u_i'} for 1≤ i≤ k, and ends at a different set {x_j,v_j,v_j'} or {y_j, u_j,u_j'} for 1≤ j≤ k, then Conditions (<ref>), (<ref>) and (<ref>) imply this;
  		
    * if it starts and ends in the same set {x_i,v_i,v_i'} or {y_i,u_i,u_i'}, then it follows from Condition (<ref>) that the path must have length at least r.
  	
  	It follows that β' must contain exactly one such path, otherwise, if two such paths are included in β', the length of β'  is greater than 2r, a contradiction to the fact that β' is a short cycle. 
  	
  	Besides this path in G∩ G', the remainder of β' must either be an edge u_iv_i or u_i'v_i', or a portion of β. 
  		If the remainder is some u_iv_i, then the distance between u_i and v_i in G is at most 2r-1, a contradiction to Condition (<ref>). The same holds if the remainder is some u_i'v_i'. 
  	
  	
  	If the remainder is a portion of β, then there exist two vertices in β connected by the path in G∩ G' contained in β'. If the two vertices are x_i, x_i+j for some 1≤ i≤ k, 1≤ j≤ k/2, then from Condition (<ref>), the path in G∩ G' contained in  β' that connects the two vertices has length at least 2r-2j+1. Since the path in β connecting x_i,x_i+j has length 2j, this implies β' has length at least 
  	(2r-2j+1)+2j=2r+1, a contradiction. In the same way, if the two vertices are y_i, y_i+j for some 1≤ i≤ k, 1≤ j≤ k/2, we can find a contradiction for β' from Condition (<ref>). 
  	
  	If the two vertices connected by the path  are x_i,y_i+j with 1≤ i≤ k, 1≤ j≤ k/2, then the path in β connecting the two vertices has length at least 2j-1. Combining the path in G∩ G' contained in β', from Condition (<ref>), we conclude that β' has length at least (2r-2j+2)+(2j-1)=2r+1, a contradiction. By the same argument, if the two vertices connected by the path are y_i,x_i+j for some 1≤ i≤ k, 1≤ j≤ k/2, we can find a contradiction that β' is not a short cycle.
  
  	Therefore such β' does not exist, and all backward switchings satisfying Conditions (<ref>)-(<ref>) are valid. 
  	
  
  There are [n]_k[m]_k/(2k) choices for the 2k-cycle β in the complete bipartite graph K_n,m, and at most (d_1(d_1-1))^k(d_2(d_2-1))^k choices for u_i,u_i',v_i, v_i', 1≤ i≤ k given β. We now count how many possible backward switchings  violate one of the four Conditions (<ref>)-(<ref>) to get a lower bound on B.  We treat the  Conditions (<ref>)-(<ref>) in 4 parts.
  
  (a) Suppose Condition (<ref>) is violated. We estimate the number of switchings by choosing one edge from the set of  edges in short cycles and the other edges arbitrarily. Note that by our definition of switchings, we identify 2k different switchings by applying the cyclic rotation or inversion on [k]. Suppose we choose an edge x_iv_i or x_iv_i' from short cycles, similar to the analysis in Lemma (<ref>), the number of switchings is at most 
  	
      (2k∑_j=2^r 2jC_j
      · (d_1-1))·(1/2k[n-1]_k-1 (d_1(d_1-1))^k-1· [m]_k(d_2(d_2-1))^k)
      
      =   2/nd_1[n]_k[m]_k[d_1(d_1-1)d_2(d_2-1)]^2k∑_j=2^r jC_j.
  
  Similarly, if we choose an edge y_iu_i or y_iu_i' from short cycles,  the number of switchings is at most 
  
      (2k∑_j=2^r 2jC_j
      · (d_2-1))·(1/2k[m-1]_k-1(d_2(d_2-1))^k-1· [n]_k (d_1(d_1-1))^k)
      
      =   2/nd_1[n]_k[m]_k[d_1(d_1-1)d_2(d_2-1)]^2k∑_j=2^r jC_j.
  
    Combining two parts, the number of switchings that violate Condition (<ref>) is at most
    
      8k/nd_1[n]_k [m]_k[d_1(d_1-1) d_2(d_2-1)]^k∑_j=2^r jC_j.
  
  
  
  (b) Suppose for some 1≤ i≤ k, two paths v_ix_iv_i' and u_iy_iu_i' are within distance 2r-1. The number of switching is at most
  
      [n-1]_k-1[m-1]_k-1/2k[d_1(d_1-1)d_2(d_2-1)]^k-1· (knd_1(d_1-1))· O((d_1-1)^r(d_2-1)^r+1) 
      
          =   1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r).
  
  
  
  (c) Suppose for some 1≤ i≤ k,1≤ j≤ k/2, two paths  {v_ix_iv_i', v_i+jx_i+jv_i+j'}  are within distance 2r-2j. The number of switchings is at most
   
      [n-2]_k-2[m]_k/2k(d_1(d_1-1))^k-2(d_2(d_2-1))^k· nd_1(d_1-1) ∑_i=1^k∑_j=1^⌊ k/2⌋O((d_1-1)^r-j+2(d_2-1)^r-j)
      
       =   1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r-1).
  
  Suppose for some 1≤ i≤ k, 1≤ j≤ k/2, two paths {u_iy_iu_i', u_i+jy_i+ju_i+j'} are within distance 2r-2j. Similarly, the number of switchings is bounded by
   
      1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r-1).
  
  Therefore the number of switchings that violate Condition (<ref>) is at most
    
      1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r-1).
  
   
   
  (d) Suppose two paths v_ix_iv_i', u_i+jy_i+ju_i+j' for some 1≤ i≤ k, 1≤ j≤ k/2  are within distance  2r-2j+1. The number of choices  is at most 
   
      [n]_k[m-1]_k-1/2k(d_1(d_1-1))^k(d_2(d_2-1))^k-1∑_i=1^k∑_j=1^⌊ k/2 ⌋O((d_1-1)^r-j+1(d_2-1)^r-j+2
      
       =   1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r).
  
  Suppose two paths u_iy_iu_i', v_i+jx_i+jv_i+j' for some 1≤ i≤ k, 1≤ j≤ k/2  are within distance  2r-2j+1. By the same argument, the number of choices  is at most
   
      1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r).
  
   Then the number of switchings that violate Condition (<ref>) is at most
   
      1/nd_1[n]_k[m]_k(d_1(d_1-1))^k(d_2(d_2-1))^kO((d_1-1)^r(d_2-1)^r).
  
   
  From (<ref>), (<ref>), (<ref>) and (<ref>), the lower bound of B is given by
   
      B≥	[n]_k[m]_k/2k(d_1(d_1-1))^k(d_2(d_2-1))^k( 1-8k  ∑_j=2^r jC_j+O(k(d_1-1)^r(d_2-1)^r)/nd_1).
  
   
   By Lemma <ref> (b),
   
      𝔼C_k≤[n]_k[m]_k/2kc_1(d_1-1)^k(d_2-1)^k/n^k m^k≤c_1(d_1-1)^k(d_2-1)^k/2k.
  
   Applying the inequality above to (<ref>), we obtain
   
      𝔼B≥[n]_k[m]_k/2k(d_1(d_1-1))^k(d_2(d_2-1))^k( 1-O(k(d_1-1)^r(d_2-1)^r)/nd_1).
  
   
  By the exchangeability of the vertex labels in the uniformly distributed RBBG model, the law of B_β is the same for any 2k-cycle β. Then
  
      𝔼B_α=2k/[n]_k[m]_k𝔼 B≥ (d_1(d_1-1)d_2(d_2-1))^k(1-c_2k(d_1-1)^r(d_2-1)^r/nd_1),
   for an absolute constant c_2>0. This completes the proof.
  
  
  
  
  
   §.§ Poisson approximation of cycle counts
  
  
  In this section, we prove the cycle counts in RBBGs are asymptotically distributed as Poisson random variables. The main tool we will use is the following  total variation distance bound from <cit.>.
  
  
  Let W=(W_1,…,W_r) be a random vector taking values in ℕ^r, and let the coordinates of Z=(Z_1,…,Z_r) be independent Poisson random variables with 𝔼Z_k=μ_k. Let W'=(W_1',…, W_r') be defined on the same space as W, with (W,W') an exchangeable pair. For any choice of σ-algebra ℱ with respect to which W is measurable and any choice of constants c_k, we have
  
      d_TV(W,Z)≤∑_k=1^rξ_k( 𝔼|μ_k-c_kℙ(Δ_k^+|ℱ)|+𝔼|W_k-c_kℙ(Δ_k^-|ℱ)|),
  
  where ξ_k:=min{ 1,1.4μ_k^-1/2} and 
  
      Δ_k^+:=   {W_k'=W_k+1, W_j=W_j', k<j≤ r},
      Δ_k^-:=   {W_k'=W_k-1, W_j=W_j', k<j≤ r}.
  
  
  
  We apply Stein's method to obtain the following Poisson  approximation in total variation distance.
  
  
  Let G be a random (d_1,d_2)-biregular bipartite graph with cycle counts (C_k, k≥ 2). Let (Z_k, k≥ 2) be independent Poisson random variables with 
      μ_k:=𝔼Z_k=(d_1-1)^k(d_2-1)^k/2k.
  	 For any n,m≥ 1 and r≥ 2, d_1≥ 3, there exists an absolute constant c_6>0 such that 
  
      d_TV((C_2,…,C_r),(Z_2,…,Z_r))≤c_6√(r)(d_1-1)^3r/2(d_2-1)^3r/2/nd_1.
  
  
  
   
  If d_1>n^1/3 or r>n^1/10, then 
      c_6√(r)(d_1-1)^3r/2(d_2-1)^3r/2/nd_1
      >1
   for a sufficiently large choice of c_6 and the theorem holds trivially. Thus we assume d_1≤ n^1/3 and r≤ n^1/10.
  We now construct an exchangeable pair of random biregular bipartite graphs by taking a step in a reversible Markov chain. 
  
  
  Define a graph 𝒢 whose vertex set consists of all (d_1,d_2)-biregular bipartite graphs.	
   If there is a  valid forward  or backward α-switching  from a (d_1,d_2)-biregular bipartite graph G_0 to another graph G_1 with the length of α being 2k, we make an undirected edge in 𝒢 between G_0,G_1 and place a weight of 
      1/[n]_k[m]_k(d_1d_2)^k
   on each such edge.  Define the degree of a vertex in 𝒢 to be the sum of weights from all adjacent edges. Let d_0 be the largest degree in 𝒢. To make 𝒢 regular, we add a weighted loop to each vertex if necessary to increase the degree of all vertices to d_0. 
  
  Now consider the simple random walk on 𝒢. This is a reversible Markov chain with respect to the uniform distribution on (d_1,d_2)-biregular bipartite graphs. Thus suppose G is a uniformly chosen random biregular bipartite graph, we can obtain another random biregular bipartite graph G' by taking an extra step in the random walk from G, and the pair (G,G') is exchangeable.
  
  
  Let 𝒥_k be the collection of cycles of length 2k in K_n,m with k≤ r. We have  |𝒥_k|= [n]_k[m]_k/2k. Define I_α=1{α⊆ G}.
  Then  C_k=∑_α∈𝒥_k1_α. Let I_α', C_k' be defined  on G' in the same way. Since G and G' are exchangeable, the vectors (C_2,…,C_r) and (C_2',…,C_r') are also exchangeable. We can then apply Lemma <ref> to this exchangeable pair of vectors. Now define two events
  
      Δ_k^+:=   {C_k'=C_k+1, C_j=C_j', k<  j≤ r},
      Δ_k^-:=   {C_k=C_k'+1, C_j=C_j', k<  j≤ r}.
  
  
  Through our construction of the exchangeable pair,
  
      ℙ(Δ_k^+| G)   =∑_α∈𝒥_kB_α/d_0[n]_k[m]_k(d_1d_2)^k	,
      ℙ(Δ_k^-| G)   =∑_α∈𝒥_kF_α/d_0[n]_k[m]_k(d_1d_2)^k.
  
  
  
  Applying Lemma <ref> with all c_k=d_0, 1≤ k≤ r, we have 
  
      d_TV((C_3,…,C_r),(Z_3,…,Z_r)) 
      ≤   ∑_k=2^rξ_k𝔼| μ_k-∑_α∈𝒥_kB_α/[n]_k[m]_k(d_1d_2)^k|+	∑_k=2^rξ_k 𝔼| C_k-∑_α∈𝒥_kF_α/[n]_k[m]_k(d_1d_2)^k| 
      
        =    ∑_k=2^rξ_k𝔼| ∑_α∈𝒥_k(d_1-1)^k(d_2-1)^k/[n]_k[m]_k-B_α/[n]_k[m]_k(d_1d_2)^k|+	∑_k=2^rξ_k 𝔼| ∑_α∈𝒥_kI_α-F_α/[n]_k[m]_k(d_1d_2)^k| 
      ≤   ∑_k=2^rξ_k (∑_α∈𝒥_k𝔼| (d_1-1)^k(d_2-1)^k/[n]_k[m]_k-B_α/[n]_k[m]_k(d_1d_2)^k|+ ∑_α∈𝒥_k𝔼| I_α-F_α/[n]_k[m]_k(d_1d_2)^k| ).
  
    
    For the rest of the proof, we estimate the following two sums
    
      ∑_α∈𝒥_k𝔼| (d_1-1)^k(d_2-1)^k/[n]_k[m]_k-B_α/[n]_k[m]_k(d_1d_2)^k|,
         ∑_α∈𝒥_k𝔼| I_α-F_α/[n]_k[m]_k(d_1d_2)^k|
  
  from (<ref>) in different ways. 
  
  (1) The upper bound on (<ref>).
  From Lemma <ref>, for all α∈𝒥_k, 
  
  
  
      𝔼| (d_1-1)^k(d_2-1)^k/[n]_k[m]_k-B_α/[n]_k[m]_k(d_1d_2)^k|   =	(d_1-1)^k(d_2-1)^k/[n]_k[m]_k-𝔼B_α/[n]_k[m]_k(d_1d_2)^k
         ≤c_2k(d_1-1)^r+k(d_2-1)^r+k/nd_1[n]_k[m]_k,
  
  where the first line is from (<ref>) and the second line is from (<ref>).
  Therefore   (<ref>) satisfies
  
      ∑_α∈𝒥_k𝔼| (d_1-1)^k(d_2-1)^k/[n]_k[m]_k-B_α/[n]_k[m]_k(d_1d_2)^k|	≤c_2[(d_1-1)(d_2-1)]^r+k/2nd_1.
  
  
  (2) The upper bound on (<ref>).
  To bound the  summation in (<ref>), for a given short cycle α, we consider a partition of 𝒢 in the following way:
  
      A_1^α   = {G does not contain α},
      
      A_2^α   ={G  contains α, which does not share an edge with another short cycle in G},
      
      A_3^α   ={G contains α, which shares an edge with another short cycle in G}.
  
  
  Conditioned on A_1^α, we have  I_α=F_α=0. Conditioned on A_2^α,  both the upper and lower bounds in Lemma <ref> can apply, which yield the following inequality:
  
      |I_α-F_α/[n]_k[m]_k(d_1d_2)^k|   ≤4k∑_j=2^rj C_j+c_1k(d_1-1)^r(d_2-1)^r/nd_1.
  
  Conditioned on A_3^α, we have I_α=1,F_α=0. 
  
  With the partition of 𝒢, the following inequality holds:
  
      𝔼|I_α-F_α/[n]_k[m]_k(d_1d_2)^k|    =  𝔼[1_A_2^α|I_α-F_α/[n]_k[m]_k(d_1d_2)^k|]+ℙ(A_3^α)
         ≤2k/nd_1𝔼[1_A_2^α∑_j=2^r2jC_j]+c_1k(d_1-1)^r(d_2-1)^r/nd_1ℙ(A_2^α)+ℙ(A_3^α).
  
  
  
  Let 𝒥_α be the set of all short cycles in K_n,m that share no edges with α. On the event A_2^α, the graph G contains no short cycles outside 𝒥_α except for α. Define |β| be the length of the cycle β. Then 
  
      ∑_j=2^r 2jC_j=2k+∑_β∈𝒥_α |β|I_β.
  
  Therefore the right-hand side of (<ref>) can be bounded by
  
      4k^2/nd_1ℙ (A_2^α)+2k/nd_1∑_β∈𝒥_α |β|𝔼I_αI_β+c_1k(d_1-1)^r(d_2-1)^r/nd_1ℙ(A_2^α)+ℙ(A_3^α) 
      ≤   4k^2/nd_1ℙ (α⊆ G)+c_1k(d_1-1)^r(d_2-1)^r/nd_1ℙ(α⊆ G)+2k/nd_1∑_β∈𝒥_α |β|𝔼I_αI_β+ℙ(A_3^α).
  
  By Lemma <ref>(1),
  
      4k^2/nd_1ℙ (α⊆ G)    =O( k^2[(d_1-1)(d_2-1)]^k/nd_1(nm)^k),
      c_1k(d_1-1)^r(d_2-1)^r/nd_1ℙ(α⊆ G)    =O( k[(d_1-1)(d_2-1)]^k+r/nd_1(nm)^k).
  
  Hence the first and the second term in (<ref>) combine to yield a corresponding upper bound in (<ref>) of
  
      ∑_α∈𝒥_k(4k^2/nd_1ℙ (α⊆ G)+c_1k(d_1-1)^r(d_2-1)^r/nd_1ℙ(α⊆ G) )= O( [(d_1-1)(d_2-1)]^k+r/nd_1).
  
  
  From Lemma <ref> (3), we have for any β∈𝒥_α with |β|=2j, 
      𝔼I_αI_β=ℙ(α∪β∈ G)≤c_1[(d_1-1)(d_2-1)]^j+k/(nm)^j+k.
  
  For 2≤ j≤ r, there are at most [n]_j[m]_j/(2j) cycles in 𝒥_α of length 2j. The third term in (<ref>)  then satisfies
  
      2k/nd_1∑_β∈𝒥_α |β|𝔼I_αI_β   ≤2k/nd_1∑_j=2^r[n]_j[m]_j/2j· 2j·c_1[(d_1-1)(d_2-1)]^j+k/(nm)^j+k
         =O( k[(d_1-1)(d_2-1)]^r+k/nd_1(nm)^k).
  
  Summing over all possible α∈𝒥_k, we obtain a corresponding term in (<ref>) of
  
      ∑_α∈𝒥_k2k/nd_1∑_β∈𝒥_α |β|𝔼I_αI_β   =O( [(d_1-1)(d_2-1)]^r+k/nd_1).
  
  
  Now given (<ref>) and (<ref>), to control (<ref>),   it remains to  estimate ∑_α∈𝒥_kℙ(A_3^α). Let 𝒦_α be the set of all short cycles in K_n,m that share an edge with α, not including α itself. By a union bound,
  
      ∑_α∈𝒥_kℙ(A_3^α)≤∑_α∈𝒥_k∑_β∈𝒦_αℙ(α∪β⊂ G).
  
  
  From  (<ref>) in Lemma <ref>, the upper bound for ℙ(α∪β⊂ G) depends on the lengths of α, β, and the number of edges that α,β share. To get an upper bound on (<ref>), we will classify and count the number of pairs (α,β) based on the structure of α∪β.
  
  Recall α has length 2k. Suppose β has length 2j. Let H=(V(α)∩ V(β), E(α)∩ E(β))  be the intersection of α and β.
   Suppose H has p components and f edges.  Since H is the intersection of two different cycles, H must be a forest with p+f vertices. So α∪β has 2j+2k-p-f vertices and 2j+2k-f edges. 
   Let a,b be the number of vertices in α∪β that are  from V_1 and V_2, respectively.  Then 
      a+b=2j+2k-p-f.
  
   
   
    Let v_1,v_2 be the number of vertices in V_1 and V_2 for H, respectively. Then we have
  a =j+k-v_1, 
    b =j+k-v_2, and  |a-b|=|v_1-v_2|. 
    Note that each component in H is a path. For each path, the difference between the number of vertices from V_1  and V_2 is at most 1. This implies 
   
      |a-b|=|v_1-v_2|≤ p.
  
   
   
   From the proof of Corollary 21 in  <cit.>, the number of all possible isomorphism types of   α∪β given |α|, |β|≤ 2r  and p,f≤ 2r is at most 
   
      (16r^3)^p-1/((p-1)!)^2.
   
   For each  isomorphism type, as a subgraph in K_n,m, the number of ways to label it is at most 
  [n]_a[m]_b+[n]_b[m]_a,
   where the two terms  come from assigning vertices in V_1,V_2 in two ways (pick an arbitrary starting vertex, decide whether it is from V_1 or V_2, then choose labels accordingly).
   
   
   From  (<ref>), (<ref>), and the assumption that n≤ m, we have that when f is even,
   
      [n]_a[m]_b+[n]_b[m]_a ≤  2n^j+k-p-f/2m^j+k-f/2=2n^-p(nm)^j+k-f/2.
  
   And when f is odd,
  
      [n]_a[m]_b+[n]_b[m]_a ≤ 2n^-p+1(nm)^j+k-f/2-1/2.
  
  
  By  (<ref>) in Lemma <ref>, the probability of any realization of the isomorphism type as a subgraph in G is bounded by 
   
      c_1[(d_1-1)(d_2-1)]^j+k-f/2/(nm)^j+k-f/2.
  
   With all the estimates above, the right-hand side of (<ref>) is now bounded by  
  
      ∑_j=2^r∑_1≤ p,f≤ 2r(16r^3)^p-1/((p-1)!)^2·( [n]_a[m]_b+[n]_b[m]_a )·c_1[(d_1-1)(d_2-1)]^j+k-f/2/(nm)^j+k-f/2
      ≤   ∑_j=2^r∑_1≤ p,f≤ 2r(16r^3)^p-1/((p-1)!)^2·  (2n^-p (nm)^j+k-f/2)
      ·c_1[(d_1-1)(d_2-1)]^j+k-f/2/(nm)^j+k-f/21{f is even}
         +∑_j=2^r∑_1≤ p,f≤ 2r(16r^3)^p-1/((p-1)!)^2·  (2n^-p+1 (nm)^j+k-f/2-1/2)
      ·c_1[(d_1-1)(d_2-1)]^j+k-f/2/(nm)^j+k-f/21{f is odd}
      
      =   O( [(d_1-1)(d_2-1)]^k+r-1/n)+O( [(d_1-1)(d_2-1)]^r+k-1/2/(nm)^1/2)
      
      =   O( [(d_1-1)(d_2-1)]^k+r/nd_1).
  
  
  Combining all estimates from (<ref>), (<ref>) and (<ref>), we finally obtain  
  
      ∑_α∈𝒥_k𝔼| I_α-F_α/[n]_k[m]_k(d_1d_2)^k|   =O( [(d_1-1)(d_2-1)]^r+k/nd_1).
  
  This provides an upper bound for (<ref>).
  
  
  (3) The upper bound on (<ref>).
  Now the upper bounds on (<ref>) and (<ref>) have been provided in  (<ref>) and  (<ref>), respectively. We are ready to estimate (<ref>). Recall
  
      ξ_k=min{ 1,1.4μ_k^-1/2}=2.8√(k)/[(d_1-1)(d_2-1)]^k/2.
  
  Then from (<ref>) and (<ref>), there is an absolute constant c_7>0 such that (<ref>) is bounded by 
  
      ∑_k=2^r c_7√(k)[(d_1-1)(d_2-1)]^r+k/2/nd_1=O( √(r)[(d_1-1)(d_2-1)]^3r/2/nd_1).
   This completes the proof.
  
  
  
  
  
  
  
  
  
  
   §.§ Cyclically non-backtracking walks and the Chebyshev polynomials
  
  
  
  In this section, we study non-backtracking walks in  biregular bipartite graphs and relate them to the Chebyshev polynomials. The relation will be used in Section <ref> to study eigenvalue fluctuations for random biregular bipartite graphs.
  
  
  
  We define a non-backtracking walk of length 2k in a biregular bipartite graph to be a walk (u_1,v_1,…, u_k,v_k,u_k+1) such that u_i∈ V_1, v_i∈ V_2, u_i+1≠u_i, for all 1≤ i≤ k and v_i+1≠v_i for all 1≤ i≤ k-1. Note that in our definition, all such walks start  and end at some vertices from V_1.
  
  
  
  
  
  A walk of length 2k denoted by 
      (u_1,v_1,…, u_k,v_k,u_k+1)
   is closed if u_k+1=u_1.
  A cyclically non-backtracking walk is a closed non-backtracking walk such that its last two steps are not the reverse of its first two steps. Namely, (u_1,v_1,u_2)≠(u_k+1,v_k,u_k).  Figure <ref> gives an example of a closed non-backtracking walk that is not cyclic non-backtracking.
  
  
  
  
  
  Let G_n be a random (d_1,d_2)-biregular bipartite graph and C_k^(n) be the number of cycles of length 2k in G_n. Denote NBW_k^(n) to be the number of non-backtracking walk of length 2k, and CNBW_k^(n) to be the number of cyclically non-backtracking walks of length 2k in G_n. 
  Let (C_k^(∞),k≥ 2)  be independent Poisson random variables with mean 
      μ_k=[(d_1-1)(d_2-1)]^k/2k.
   We also define C_1^(∞)=C_1^(n)=0. For k≥ 1,
  denote 
  
      _k^(∞)=∑_j| k 2jC_j^(∞).
    	
  
  For any cycle of length 2j in G_n with j| k, we can obtain 2j cyclically non-backtracking walks by choosing a starting point from V_1, fixing a direction and then walking around the cycle of length 2k repeatedly. The next lemma shows that CNBW_k^(n) can be approximated by the count of those repeated walks around cycles.
  
  
  Let G_n be a random (d_1,d_2)-biregular bipartite graph.
  	Suppose d_1≤ n^1/3,k≤ n^1/10, define
  	
      B_k^(n)=_k^(n)-∑_j| k 2jC_j^(n)
   
  to be the number of cyclically non-backtracking walks of length 2k in G_n that are not repeated walks around cycles. Then 
  
      𝔼B_k^(n)≤c_7k^7[(d_1-1)(d_2-1)]^k/n.
  
  
  
  
  	We call a cyclically non-backtracking walk bad if it's not a repeated walk on a cycle. Then from (<ref>), B_k^(n) counts the number of bad cyclically non-backtracking walks of length 2k.  
  	
  	Let (w_0,w_1,…, w_2k) with w_2k=w_0∈ V_1 be a bad cyclically non-backtracking walk in K_n,m of length 2k. For any 1≤ i≤ 2k, we say that the i-th step of the walk is  
  	
  	    
    * free if w_i did not previously occur in the walk;
  	    
    * a coincidence if w_i previously occurred in the walk, but the edge w_i-1w_i didn't;
  	    
    * forced if the edge w_i-1w_i previously occurred in the walk.
  	
  	
  Let χ+1 be the number of coincidences and f be the number of forced steps in the walk. Let χ_1+1 and χ_2 be the number of coincidence steps ending at a vertex from V_1 and V_2, respectively. Let f_1,f_2 be the number of forced steps ending at a vertex from V_1 and V_2, respectively. Denote v,e  the number of distinct vertices and edges in the cyclically non-backtracking walk, respectively.  We now have the following relations:
  
      χ+1   =χ_1+χ_2+1,
      
          f   =f_1+f_2,
      
          v   =(2k+1)-(χ+1)-f=2k-χ-f,
      
          e   =2k-f.
  
  For any repeated walk on a cycle, the number of coincidences is 1 and χ=0. Therefore if the walk is bad, we must have χ≥ 1.
  
  The following lemma bounds the number of cyclically non-backtracking walks with given parameters χ_1,χ_2,f_1, and f_2.
  
      Consider cyclically non-backtracking walks of length 2k on K_n,m such that in the subgraph spanned by this walk, all vertices from V_1 have degrees at most d_1 and vertices from V_2 have degrees at most d_2. Then the number of such walks with given χ_1,χ_2, f_1,f_2 satisfying χ≥ 1  is at most 
      (2k)^3(χ_1+χ_2)+2(d_1-1)^f_2(d_2-1)^f_1 n^k-χ_1-f_1m^k-χ_2-f_2.
   Moreover, we must have |f_1-f_2|≤χ+1.
  
   
  We count the number of such cyclically non-backtracking walks by choosing the coincidences, forced steps, and free steps separately. 
  Given that there are χ+1 coincidences, there are 2kχ+1 many possible subsets of indices in {1,…,2k} where coincidences can happen. The vertices at a coincidence have already occurred in the walk, so there are at most 2k choices for each of them, giving us a total of 
  2kχ+1(2k)^χ+1≤ (2k)^2χ+2 many choices. 
  
  For forced steps, they can only occur after a coincidence or another forced step. After each coincidence, imagine assigning some number of steps to be forced. The number of ways to do this is at most the number of weak compositions of f elements into χ+1 parts, which is f+χχ≤ (2k)^χ. For each forced step ending at a vertex from V_1, the walk can only move along an edge that has already been traversed, so there are at most (d_2-1) possible choices 
  at every step due to the non-backtracking property. Similarly, for each forced step ending at a vertex from V_2 there are at most d_1-1 possible choices. Altogether this gives us at most (2k)^χ(d_1-1)^f_2(d_2-1)^f_1 choices for all forced steps.
  
  There are  k-χ-1-f_1 many free steps ending at a vertex from V_1, we have at most n choices for the next vertex, and we have an additional n choices for w_0∈ V_1, which gives a total of at most n^k-χ_1-f_1 many choices. Similarly, the number of free steps ending at a vertex from V_2 is at most m^k-χ_2-f_2. Multiplying together every part from coincidences, forced steps, and free steps give us at most 
  
      (2k)^3χ+2(d_1-1)^f_2(d_2-1)^f_1n^k-χ_1-f_1m^k-χ_2-f_2
  
  many such cyclically non-backtracking walks.
  
  Next, we bound |f_1-f_2|. Recall that forced steps can only occur after a coincidence or another forced step. Then there are at most χ+1 many consecutive forced steps starting from a certain coincidence step. In each consecutive forced step, the number of vertices from V_1 and V_2 differ by at most 1, since the subgraph spanned by any consecutive forced steps is a path. Hence we have |f_1-f_2|≤χ+1.
  
  
  Equipped with Lemma <ref>, we continue to prove Lemma <ref>.
  
  By Part (a) in Lemma <ref>, the probability that a given bad walk appears in G_n is at most 
      c_1((d_1-1)(d_2-1)/nm)^k-f/2.
   From the upper bound on the number of such walks in Lemma <ref>, summing over all possibilities of χ_1,χ_2,f_1,f_2, we  have 
      𝔼B_k^(n)
      ≤   ∑_χ_1,χ_2:
      χ_1+χ_2≥ 1∑_0≤ f_1,f_2≤ k-1
       |f_1-f_2|≤χ+1 (2k)^3χ+2(d_1-1)^f_2(d_2-1)^f_1n^k-χ_1-f_1m^k-χ_2-f_2c_1 ((d_1-1)(d_2-1)/nm)^k-f/2
      
          =   c_1[(d_1-1)(d_2-1)]^k∑_χ_1+χ_2≥ 1n^-χ_1m^-χ_2 (2k)^3(χ_1+χ_2)+2∑_0≤ f_1,f_2≤ k-1
       |f_1-f_2|≤χ+1((d_2-1)m/(d_1-1)n)^(f_1-f_2)/2.
  
  
  Since  (d_2-1)d_1≤ (d_1-1)d_2, the following inequality holds:
  
      ∑_0≤ f_1,f_2≤ k-1
       |f_1-f_2|≤χ+1((d_2-1)m/(d_1-1)n)^(f_1-f_2)/2=   ∑_0≤ f_1,f_2≤ k-1
       |f_1-f_2|≤χ+1((d_2-1)d_1/(d_1-1)d_2)^(f_1-f_2)/2
      ≤    k^2((d_1-1)d_2/(d_2-1)d_1)^(χ+1)/2.
  
  Since d_1≤ n^1/3, k≤ n^1/10,  (<ref>) implies 
  
      𝔼B_k^(n)   ≤ c_1 k^2[(d_1-1)(d_2-1)]^k∑_χ_1+χ_2≥ 1n^-χ_1m^-χ_2 (2k)^3(χ_1+χ_2)+2((d_1-1)d_2/(d_2-1)d_1)^(χ+1)/2
         = k^2[(d_1-1)(d_2-1)]^k O( (2k)^5(d_1-1)d_2/n(d_2-1)d_1)
         =O( k^7(d_1-1)^k(d_2-1)^k/n).
  
  This completes the proof of Lemma <ref>.
  
  
  Recall the definition of _k^(∞) from (<ref>). The following corollary holds.
  
   
  Suppose d_1≤ n^1/3 and r≤ n^1/10. There exists a constant c_8>0 such that
  
      d_TV( (_k^(n),2≤ k≤ r),(_k^(∞),2≤ k≤ r))≤c_8√(r)[(d_1-1)(d_2-1)]^3r/2/nd_1.
  
  
  
  
  	By the definition of total variation distance,  for any measurable map f and random variable X,Y, we have 
  	
      d_TV(f(X),f(Y))≤ d_TV(X,Y).
  
  It follows from Theorem  <ref> that 
  	
      d_TV( (  ∑_j| k 2jC_j^(n),2≤ k≤ r),(_k^(∞ ),2≤ k≤ r ))≤c_6√(r)[(d_1-1)(d_2-1)]^3r/2/nd_1.
  
  By Markov's inequality and Lemma <ref>,
  
      ℙ(B_k^(n)≥ 1)≤c_7k^7[(d_1-1)(d_2-1)]^k/n.
   
  Summing these probabilities for k=2,…,r implies
  
      (  ∑_j| k 2jC_j^(n),2≤ k≤ r)=(_k^(n), 2≤ k≤ r)
  
  with probability 1-O(r^7[(d_1-1)(d_2-1)]^r/n).  Therefore by the coupling inequality,
  
      d_TV((  ∑_j| k 2jC_j^(n),2≤ k≤ r),(_k^(n), 2≤ k≤ r))=O(r^7[(d_1-1)(d_2-1)]^r/n).
  
  From (<ref>) and (<ref>), 
  
      d_TV( (_k^(n),2≤ k≤ r),(_k^(∞),2≤ k≤ r))
      ≤   c_6√(r)[(d_1-1)(d_2-1)]^3r/2/nd_1+ O(r^7[(d_1-1)(d_2-1)]^r/n)=O(√(r)[(d_1-1)(d_2-1)]^3r/2/nd_1).
  
  
  
  Let λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)). For the rest of this section, we  connect the spectrum of XX^⊤-d_1I/√((d_1-1)(d_2-1)) with Chebyshev polynomials and cyclically non-backtracking walks.
  Define 
  
      Γ_0(x)   =1, 
      Γ_2k(x)   =2T_2k(x/2)+d_1-2/(d_1-1)^k,
      Γ_2k+1(x)   =2T_2k+1(x/2).
  
  Here {T_k(x)} are  the  Chebyshev polynomials of the first kind on [-1,1] which satisfy
  
      T_0(x)   =1,    
          T_1(x) =x, 
      
          T_k+1(x)    =2xT_k(x)-T_k-1(x).
  
  Let {U_k(x)} be the Chebyshev polynomials of the second kind on [-1,1] such that
  
      U_-1(x)    =0,   U_0(x) = 1, 
      
      U_k+1(x)    =2xU_k(x)-U_k-1(x).
  
  Define
  
      p_k(x)=U_k(x/2)-1/d_1-1U_k-2(x/2).
  
  We begin with representing closed non-backtracking walks with p_k(x). The following lemma gives a deterministic identity. Recall in our Definition <ref>, all closed non-backtracking walks start and end at vertices in V_1.
  
  Let NBW_k^(n) be the number of closed non-backtracking walks of length 2k in a (d_1,d_2)-biregular bipartite graph G.  Let λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)).
  We have
  
      ∑_i=1^n p_k(λ_i)=(d_1-1)^-k/2(d_2-1)^-k/2NBW_k^(n).
  
  
  
  Let A^(k) be the n× n matrix such that A_ij^(k) is the number of non-backtracking walks of length 2k from i to j, where i,j∈ V_1. We have the following relations:
  
      A^(1)   =XX^⊤-d_1I,    A^(2)=(A^(1))^2-d_1(d_2-1)I, 
      
          A^(k+1)   =A^(1) A^(k)-(d_1-1)(d_2-1)A^(k-1),   ∀ k≥ 2.
  
  
  The expressions of A^(1) and A^(2) follow from the definition of non-backtracking walks. Since a non-backtracking walk of length 2k+2 can be decomposed as a non-backtracking walk of length 2k and a non-backtracking walk of length 2 which avoid backtracking at the 2k-th step, the expression (<ref>) holds. 
  We now claim that for k≥ 1,
  
      p_k(XX^⊤-d_1I/√((d_1-1)(d_2-1)))=[(d_1-1)(d_2-1)]^-k/2A^(k),
  
  and prove it by induction. Note that from (<ref>),
  
      p_1(x)=x,    p_2(x)=x^2-1-1/d_1-1.
  
  It is easy to check (<ref>) holds for k=1,2. 
  Since p_k(x) is a linear combination of U_k(x/2) and U_k-2(x/2), it satisfies the recursive relation for U_k(x/2), which is 
  
      p_k(x)=xp_k(x)-p_k-1(x).
  
  
  
  Assume (<ref>) holds for k≤ s. Let M=XX^⊤-d_1I. Then
  
      p_s+1(M/√((d_1-1)(d_2-1)))
          =   M[(d_1-1)(d_2-1)]^-(s+1)/2A^(s)-[(d_1-1)(d_2-1)]^-(s-1)/2A^(s-1)
      
          =   [(d_1-1)(d_2-1)]^-(s+1)/2(MA^(s)-(d_1-1)(d_2-1)A^(s-1))
      
          =   [(d_1-1)(d_2-1)]^-(s+1)/2 A^(s+1),
  
  where the last equality is from (<ref>). Therefore (<ref>) holds. Taking  trace on both sides in (<ref>), we obtain (<ref>).
  
  
  The next theorem  is an algebraic relation  between Γ_k and the number of cyclic non-backtracking walks. Together with Lemma <ref>, it implies the  polynomials Γ_k(x) of the eigenvalues for RBBGs converges in distribution to a sum of Poisson random variables.
  
  Let G be a (d_1,d_2)-biregular bipartite graph and λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)). Then  for any k≥ 1, we have
  	
      ∑_i=1^n Γ_k(λ_i)=(d_1-1)^-k/2(d_2-1)^-k/2CNBW_k^(n).
  
  
  
  
  
  
  
  
  
  
  
  We first relate the number of cyclically non-backtracking closed walks _k^(n) to the number of  closed non-backtracking walks  NBW_k^(n). 
  
  
  A closed non-backtracking walk of length 2k is either cyclically non-backtracking, or it can be obtained from a closed non-backtracking walk of length 2(k-2) by adding a new walk of length 2 (which we call a tail) to the beginning of the walk and its reverse to the end (see Figure <ref> for an example). For any cyclically non-backtracking walk of length 2(k-2), we can add a tail in (d_1-2)(d_2-1) many ways. For any closed non-backtracking walk of length 2(k-2) that is not cyclically non-backtracking, we can add a tail in (d_1-1)(d_2-1) many ways. 
  Therefore for k≥ 3, we have the following equation
  
      NBW_k^(n)   =CNBW_k^(n)	+(d_1-2)(d_2-1)CNBW_k-2^(n)+(d_1-1)(d_2-1)(NBW_k-2^(n)-CNBW_k-2^(n))
         =CNBW_k^(n)+(d_1-1)(d_2-1)NBW_k-2^(n)-(d_2-1)CNBW_k-2^(n),
  
  which can be written as 
  
      CNBW_k^(n)-(d_2-1)CNBW_k-2^(n)=NBW_k^(n)-(d_1-1)(d_2-1)NBW_k-2^(n).
  
  
  
  Note that CNBW_k^(n)=NBW_k^(n) for k=1,2. Applying (<ref>) recursively,  we have when k is even,
  
      CNBW_k^(n)
      
      =   NBW_k^(n)-(d_1-2) [(d_2-1)NBW_k-2^(n)	+(d_2-1)^2NBW_k-4^(n)	+⋯ +(d_2-1)^k-2/2NBW_2^(n))].
  
  And when k is odd,
  
      CNBW_k^(n)
      
      =   NBW_k^(n)-(d_1-2)[(d_2-1)NBW_k-2^(n)	+(d_2-1)^2NBW_k-4^(n)	+⋯ +(d_2-1)^k-3/2NBW_3^(n)	].
  
  
  Denote 
  
      NBW_k^(n):=   (d_2-1)^-k/2NBW_k^(n)  ,   CNBW_k^(n):=(d_2-1)^-k/2CNBW_k^(n).
  
   We can simplify the above equations (<ref>) and (<ref>) as
  
      CNBW_k^(n)=NBW_k^(n)-(d_1-2)( NBW_k-2^(n)+NBW_k-4^(n)+⋯+NBW_a^(n)),
  
  where a=2 if k is even and a=1 if k is odd. Also (<ref>) can be written as 
  
      ∑_i=1^n p_k(λ_i)=(d_1-1)^-k/2NBW_k^(n).
  
  
  
  From the proof of  Proposition 32 in <cit.>, we have the following relation between Γ_k(x) and p_k(x) for k≥ 1:
  
      Γ_2k(x)   =p_2k(x)-(d_1-2)( p_2k-2(x)/d_1-1+p_2k-4(x)/(d_1-1)^2+⋯+p_2(x)/(d_1-1)^k-1), 
      Γ_2k-1(x)   =p_2k-1(x)-(d_1-2)( p_2k-3(x)/d_1-1+p_2k-5(x)/(d_1-1)^2+⋯+p_1(x)/(d_1-1)^k-1).
  
  Then from (<ref>) and (<ref>),
  
      (d_1-1)^-kCNBW_2k^(n)   =∑_i=1^n(  p_2k(λ_i)-(d_1-2) ( p_2k-2(λ_i)/d_1-1+⋯+p_2(λ_i)/(d_1-1)^k-1))=∑_i=1^n Γ_2k(λ_i),
  
  where the last equality is from (<ref>). Similarly, from (<ref>),
  
      (d_1-1)^(2k-1)/2CNBW_2k-1^(n)=∑_i=1^n Γ_2k-1(λ_i).
  
  
  Therefore for all k≥ 1,
  
      ∑_i=1^n Γ_k(λ_i)=(d_1-1)^-k/2CNBW_k^(n)=[(d_1-1)(d_2-1)]^-k/2CNBW_k^(n).
  
  This completes the proof of Theorem <ref>.
  
  
  
  
  
  § SPECTRAL GAP
  
  
  In this section, we  provide some estimates on the second largest eigenvalue of the random biregular bipartite graphs that will be used to study eigenvalue fluctuations in Section <ref>.
  Note that the largest eigenvalue of XX^⊤-d_1I is λ_1=d_1(d_2-1). In the next theorem, we provide upper bounds on |λ| for all eigenvalues λ≠λ_1.
  
  Let	G be a (d_1,d_2)-random biregular bipartite graph with d_1≥ d_2. Let λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I. 
  
  	
    * For fixed d_1,d_2, there exists a sequence ϵ_n→ 0 such that for any eigenvalue λ≠λ_1,
  	 
      ℙ( |λ-(d_2-2)|≥ 2√((d_1-1)(d_2-1))+ϵ_n)→ 0
   
  	as n→∞.
  	
    * Suppose d_2≤1/2n^2/3, d_1≥ d_2≥ cd_1 for some constant c∈ (0,1). Then for some constant α_1>0 depending on c and any eigenvalue λ≠λ_1,
  	
      ℙ( |λ |≥α_1√((d_1-1)(d_2-1)))≤1/n^2.
   
  
    * Suppose d_2≤ C_1, d_1≤ n^2, there exists a constant α_2 depending on C_1 such that for any eigenvalue λ≠λ_1,
  	
      ℙ( |λ|≥α_2√((d_1-1)(d_2-1)))≤1/n^2.
   
  
  
  
  The probability estimates in (<ref>) and (<ref>) can be improved, see <cit.>. In order to prove the main theorems in Section <ref>, we only include a weaker version for simplicity.
  
  
  
  
  Theorem 4 in <cit.> states that for a random biregular bipartite graph with d_1≥ d_2, the eigenvalues of the adjacency matrix A satisfy the following estimates with high probability:
  
      
    * the second eigenvalue of A satisfies λ_2(A)≤√(d_1-1)+√(d_2-1)+o(1),
      
    * the smallest positive eigenvalue of A satisfies λ_min^+(A)≥√(d_1-1)-√(d_2-1)-o(1).
  
  Since eigenvalues of XX^⊤ are the squares of the eigenvalues for A, we have with high probability,
  
      λ_2(XX^⊤)-d_1-(d_2-2)    ≤ 2√((d_1-1)(d_2-1))+o(1),
      λ_n(XX^⊤)-d_1-(d_2-2)    ≥ -2√((d_1-1)(d_2-1))-o(1),
  
  therefore
   (<ref>) holds.
  
  Theorem 1.1 in <cit.>
  states that if d_2≤1/2n^2/3 and d_1≥ d_2, there exists a constant α>0 such that λ_2(A)≤α√(d_1) with probability at least 1-m^-2.  This implies for any eigenvalue λ of XX^⊤-d_1I with λ≠d_1(d_2-1), we have 
  
      ℙ( -d_1≤λ≤α^2 d_1-d_1)≥ 1-m^-2≥ 1-n^-2.
  
  Since d_1≥ d_2≥ cd_1, we can find a constant α_1>0 depending on α and c such that 
  
      ℙ( |λ|≤α_1√((d_1-1)(d_2-1))) ≥ 1-n^-2.
  
  Therefore (<ref>) holds. Theorem 1.5 in <cit.> states that if d_2≤ C_1, d_1≤ n^2, there exists a constant α_2 depending on C_1 such that
  
      ℙ(max_2≤ i≤ m+n-1|λ_i^2(A)-d_1|≥α_2 √((d_1-1)(d_2-1)))≤ n^-2.
  
  Then (<ref>) follows from the algebraic relation between the spectra of A and XX^⊤-d_1I.
  
  
  
  
  § EIGENVALUE FLUCTUATIONS
  
  
  Lemma <ref> and Corollary <ref> in Section <ref> imply the limiting laws for ∑_i=1^n Γ_k(λ_i) are given by a sum of Poisson random variables. In this section, we extend the results to a more general class of function f and study the behavior of ∑_i=1^n f(λ_i) for RBBGs with fixed and growing degrees.
  
  
  
  The following set-up for weak convergence will be used in Section <ref> to prove Theorem <ref>. We will closely follow the definitions and notations used in <cit.>. See Section 2 in <cit.> for more details.
  
  
  Denote ℕ:={1,2,…}.
  Let w⃗=(w_m)_m∈ℕ be a sequence of positive weights. Let L^2(w⃗) be the space of sequences (x_m)_m∈ℕ that are square-integrable with respect to w⃗, i.e., ∑_m=1^∞ x_m^2w_m<∞. We define a complete separable metric space 𝒳=(L^2(w⃗), ·), where for any  sequence (x_m)_m∈ℕ, 
      x=(∑_m=1^∞ x_m^2w_m)^1/2.
  
  
  Denote the space of probability measures on the Borel σ-algebra of 𝒳 by 𝒫(𝒳). We use the Prokhorov metric for weak convergence as the metric on 𝒫(𝒳). The following results are proved in Section 2 of <cit.>.
    The following holds for the complete separable metric space 𝒳.
      
      
    * Let (a_m)_m∈ℕ∈ L^2(w⃗) be such that a_m≥ 0 for every m. Then the set
      
      { (b_m)_m∈ℕ∈ L^2(w⃗): 0≤ |b_m|≤ a_m, ∀ m∈ℕ}
  
      is compact in (L^2(w⃗), ·).
          
    * Suppose {X_n} and X are random sequences taking values in L^2(w⃗) such that X_n converges in distribution to X. Then for any b∈ L^2(w⃗), the random variables ⟨ b, X_n⟩ converges in distribution to ⟨ b, X⟩.
          
    * Let x∈𝒳 and P,Q be two probability measures in 𝒫(𝒳). Suppose for any finite collection of indices (i_1,…,i_k), the law of random vector (x_i_1,…, x_i_k) is the same under both P and Q. Then P=Q on the entire Borel σ-algebra of 𝒳.
      
  
  
  
  
  
  
  
  
  
  We  also need the following results from the approximation theory.
  
  For ρ>1, let E_ρ be the image of the map z↦ (z+z^-1)/2 of the open disc of radius ρ in the complex plain centered at the origin. We can E_ρ the Bernstein ellipse of radius ρ.
  The ellipse has foci at ± 1, and the sum of the major semi-axis and minor semi-axis is exactly ρ.
  
  
  
   
      Suppose f: [-1,1]→ℝ can be analytically extended to E_ρ and is bounded by M on E_ρ. Then f has a unique expansion on [-1,1] as
      
      f(x)=∑_k=0^∞ a_kT_k(x),
  
      where T_k(x) is the Chebyshev polynomial of the first kind defined in  (<ref>),
      and the coefficients of this expansion satisfy
      
      |a_0|≤ M,    |a_k|≤2M/ρ^k.
  
  
  Define f_k(x)=∑_i=0^k a_kT_k(x).
  Applying the bound |T_k(x)|≤ 1 when x∈ [-1,1]  and Proposition <ref>, we obtain for all x∈ [-1,1],
  
      |f(x)-f_k(x)|≤2M/ρ^k(ρ-1).
  
  
  
   §.§ Poisson fluctuations with fixed degrees
  
  Now fix d_1 and d_2 as constants. We are ready to extend our results in Section <ref> to a more general class of functions as follows. Note that the following theorem is given for a sequence of RBBGs with growing n. For ease of notation, we drop the dependence on n when writing the matrix X and eigenvalues λ_1,…,λ_n.
  
  
  For fixed d_1≥ d_2≥ 2 and (d_1,d_2)≠(2,2), let G_n be a sequence of  random (d_1,d_2)-biregular bipartite graph. Let λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)). Suppose f is a function such that f(2z) is analytic on E_ρ, where ρ=[ (d_1-1)(d_2-1)]^α for some α>7/2. Then f(x) can be expanded on [-2,2] as
  
      f(x)=∑_k=0^∞ a_kΓ_k(x),
  
  and the random variable
  
      Y_f^(n):=∑_i=1^n f(λ_i)-n a_0
   
  converges in distribution as n→∞ to the infinitely divisible random variable
  
      Y_f:=∑_k=2^∞a_k/[(d_1-1)(d_2-1)]^k/2_k^(∞),
   
  where _k^(∞) is defined in (<ref>).
  
  
  
  Define
  
      f_k(x):=∑_i=0^k a_iΓ_i(x).
  
  We first show that f_k(x) is a good approximation of f(x). Applying Proposition <ref> to f(2x) gives an expansion (<ref>) with 
  
      |a_k|≤ C[(d_1-1)(d_2-1)]^-α k
   
  for some constant C that depends only on d_1,d_2 and the constant M given in Proposition <ref>. 
  
  
  By the proprieties of Chebyshev polynomials, on any interval [-K,K], we have 
  
      max_|x|≤ K|T_k(x)|=(K-√(K^2-1))^k+(K+√(K^2-1))^k/2.
  
  From (<ref>) and (<ref>), we have Γ_1(x)=x, and for any k≥ 2,
  
      | Γ_k(x)|≤ 2 | T_k(x/2)|+d_1-2/(d_1-1)^k/2≤ 2| T_k(x/2)|+1.
  
  From (<ref>),
  
      max_|x|≤ 3 2|T_k(x/2)|=(3/2-√(5)/2)^k+(3/2+√(5)/2)^k,
  
  Then for all  k≥ 2, with (<ref>) we obtain
  
      sup_|x|≤ 3| Γ_k(x)|≤(3/2+√(5)/2)^k+2≤  3^k+1,
  
  and the same bound holds when k=1. From  (<ref>) and  (<ref>), for all  x∈ [-3,3],
  
      ∑_k=0^∞|a_kΓ_k(x)|≤ 3C∑_k=0^∞ [3((d_1-1)(d_2-1))^-α]^k<∞,
  
  where the last inequality comes from the fact that (d_1-1)(d_2-1)≥ 2 and α>7/2. Hence the series ∑_k=0^∞ a_kΓ_k(x) is absolutely convergent on [-3,3], which implies the expansion of f in (<ref>) is valid on [-3,3].
  Then we have for a constant C_1>0 depending on C,
  
      sup_|x|≤ 3|f(x)-f_k(x)|    ≤sup_|x|≤ 3∑_i=k+1^∞|a_iΓ_i(x)|≤ C_1[3(d_1-1)(d_2-1)^-α]^k+1.
   
   Denote 
      K_1:=λ_1=d_1√((d_2-1))/√(d_1-1).
   For sufficiently large k, from (<ref>) and (<ref>),
  
      sup_|x|≤ K_1| Γ_k(x)|≤ (2K_1)^k.
  
  And from (<ref>) and the assumption α>7/2,
  
      ∑_k=0^∞|a_kΓ_k(x)|≤ C∑_k=0^∞ [2K_1((d_1-1)(d_2-1))^-α]^k<∞.
  
  It implies the series ∑_k=0^∞ a_kΓ_k(x) is also absolutely convergent on [-K_1,K_1], and the expansion of f in (<ref>) is  valid on [-K_1,K_1]. 
  
  
  Since 2K_1≤ 4[(d_1-1)(d_2-1)]^1/2 and (d_1-1)(d_2-1)≥ 2, for a constant C_2>0,
  
      sup_|x|≤ K_1|f(x)-f_k(x)|   ≤ C_2[2K_1((d_1-1)(d_2-1))^-α]^k+1
         ≤ C_2[4((d_1-1)(d_2-1))^-α+1/2]^k+1≤ C_2[((d_1-1)(d_2-1))^-α+5/2]^k+1.
   
  Therefore f_k converges to f uniformly on [-K_1,K_1], and the interval [-K_1,K_1] deterministically contains all eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)). 
  
  
  By the definition of _k^(∞) in (<ref>), Equation (<ref>) can be written as 
  
      Y_f:=∑_j=1^∞∑_i=1^∞a_ij/[(d_1-1)(d_2-1)]^ij/22jC_j^(∞),
  
  where Y_f is a sum of independent random variables, and
  𝔼|Y_f|^2<∞ by (<ref>). 
  
  Denote α':=α-2>3/2. 
  Choose β satisfying 1/α'<β<2/3 and define
  
      r_n    =⌊βlog n/log[(d_1-1)(d_2-1)]⌋,
      
          X_f^(n)   =∑_k=1^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(n),
      Ỹ_f^(n)   =∑_k=1^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(∞).
  
  Note that _1^(n)=0, from (<ref>),
  
      X_f^(n)   =∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(n)=∑_i=1^n f_r_n(λ_i)-na_0.
  
  By Corollary <ref>, 
  
  
      d_TV( (_k^(n),2≤ k≤ r_n),(_k^(∞),2≤ k≤ r_n))≤c_8√(r_n)[(d_1-1)(d_2-1)]^3r_n/2/nd_1=o(1).
  
  
  Since X_f^(n) and Ỹ_f^(n) are measurable functions of 
      (_k^(n),2≤ k≤ r_n)   and    (_k^(∞),2≤ k≤ r_n),
   respectively, we have
  
      d_TV(X_f^(n),Ỹ_f^(n))≤ d_TV( (_k^(n),2≤ k≤ r_n),(_k^(∞),2≤ k≤ r_n))=o(1).
  
  Note that Ỹ_f^(n) converges  almost surely to Y_f by (<ref>), so X_f^(n) converges in distribution to Y_f.  
  
  By Slutsky's theorem, to show Y_f^(n) defined in (<ref>) converges in distribution to Y_f, it remains to  show that Y_f^(n)-X_f^(n) converges to zero in probability.  The largest eigenvalue of XX^⊤-d_1I/√((d_1-1)(d_2-1)) is K_1, so from (<ref>) we have 
      lim_k→∞f_k(λ_1)= f(λ_1).
   
  Then for any δ>0 and sufficiently large n,
  
      |f(λ_1)-f_r_n(λ_1)|≤δ/2.
  
   From (<ref>), (<ref>) and (<ref>), we  have for sufficiently large n,
  
      |Y_f^(n)-X_f^(n)|   ≤∑_i=1^n |f(λ_i)-f_r_n(λ_i)|≤δ/2+∑_i=2^n |f(λ_i)-f_r_n(λ_i)|.
  
  Suppose that all the non-trivial eigenvalues λ≠λ_1 are contained in [-3,3], from (<ref>), 
  
      ∑_i=2^n |f(λ_i)-f_r_n(λ_i)|   ≤ C_1(n-1)[3(d_1-1)(d_2-1)^-α]^r_n+1
         ≤ C_1n[(d_1-1)(d_2-1)^-α+2]^r_n≤ C_1n^1-α'β=o(1),
  
  which combining (<ref>) implies for sufficiently large n,
  
      |Y_f^(n)-X_f^(n)|≤δ.
  
  
  
  
  Recall (<ref>) and the assumption d_1≥ d_2. With high probability, for a sequence ϵ_n→ 0, we have the nontrivial eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)) is contained in 
  
      [ -2-ϵ_n+d_2-2/√((d_1-1)(d_2-1)),2+ϵ_n+d_2-2/√((d_1-1)(d_2-1))]⊆ [-3,3]
   
  for sufficiently large n. Therefore
  
      ℙ( |Y_f^(n)-X_f^(n)|≥δ)≤ℙ( max_2≤ i≤ n|λ_i|≥ 3 )=o(1).
  
  This finishes the proof.
  
  
  
  As a corollary of Theorem <ref>, we obtain   eigenvalue fluctuations  for the adjacency matrices of RBBGs as follows.
  
  For fixed d_1≥ d_2≥ 2 and (d_1,d_2)≠(2,2), let G_n be a sequence of  random (d_1,d_2)-biregular bipartite graph. Let λ_1≥⋯≥λ_n+m be the  eigenvalues of its adjacency matrix A. Suppose f satisfies the same conditions as in Theorem <ref>. Then the random variable
  
      Y_f^(n):=1/2[∑_i=1^n+m f(λ_i^2-d_1/√((d_1-1)(d_2-1)))-(m-n)f(-d_1/√((d_1-1)(d_2-1))) ]-n a_0
  
  converges in distribution as n→∞ to the infinitely divisible random variable
  
      Y_f:=∑_k=2^∞a_k/[(d_1-1)(d_2-1)]^k/2_k^(∞),
  
  where _k^(∞) is defined in (<ref>).
  
  
  Recall from Section <ref> that all eigenvalues of A consist of two parts. There are 2n eigenvalues  in pair as {-λ,λ} where λ is a singular value of X. In addition, there are (m-n) extra zero eigenvalues. the result then follows from the algebraic relation between eigenvalues of A and eigenvalues of XX^⊤-d_1I.
  
  
  
  
  
  
   §.§ Gaussian fluctuations with growing degrees
  
  In this section, we consider the eigenvalue fluctuations of RBBGs when d_1· d_2→∞.
  
  We first prove the following weak convergence result for a normalized and centered version of _k^(∞).
  
  
  Suppose that d_1· d_2→∞, r_n→∞ as n→∞. For k≥ 2,
  define 
  
      N_k^(n):=1/[d_1-1)(d_2-1)]^k/2(_k^(∞)-𝔼_k^(∞)) 1_{ k≤ r_n}.
  
  Let {Z_k}_k≥ 2 be independent Gaussian random variables with 𝔼 Z_k=0 and 𝔼 Z_k^2=2k.   Define the weight w_k=b_k/(k^2log (k+1)), where (b_k)_k∈ℕ is any fixed positive summable sequence. 
  
  Let P_n be the law of the sequence (N_k^(n))_k≥ 2. Then 
  as an element in 𝒫(𝒳), P_n converges weakly to the law of the random vector (Z_k)_k≥ 2. 
  
   We first prove the following 
  Claim (1): for any fixed r, (N_k^(n))_2≤ k≤ r converges in distribution to (Z_k)_2≤ k≤ r.
  
  
  For any fixed k, when n is sufficiently large,  we can write  (<ref>) as 
  
      N_k^(n)=   1/[d_1-1)(d_2-1)]^k/2(2kC_k^(∞)-[(d_1-1)(d_2-1)]^k)
         +1/[d_1-1)(d_2-1)]^k/2∑_j| k, j<k(2jC_j^(∞)-[(d_1-1)(d_2-1)]^j).
  
  
  
  Recall C_k^(∞) is a Poisson random variable with mean (d_1-1)^k(d_2-1)^k/2k.
  The first term in (<ref>) converges in distribution to a centered Gaussian random variable Z_k with variance 2k as n→∞ from the Gaussian approximation of Poisson distribution. 
  
  To show the convergence of N_k^(n) for a fixed k, it remains to show the second term in (<ref>) converges to zero in probability. Note that  the second term in (<ref>) has a zero mean and its variance is given by
  
      Var[ 1/[d_1-1)(d_2-1)]^k/2∑_j| k, j<k(2jC_j^(∞)-[(d_1-1)(d_2-1)]^j)]=∑_j| k, j<k 2j[(d_1-1)(d_2-1)]^j-k,
  
  which goes to 0 as n→∞. Then by Chebyshev's inequality, this term converges to 0 in probability. Therefore Claim (1) holds.
  
  
  
  We further define N_1^(n)=0, Z_1=0, and consider the weak convergence of (N_k^(n))_k∈ℕ as an element in L^2(w⃗). Since
  
      𝔼∑_k=1^∞(Z_k)^2w_k=∑_k=2^∞b_k/klog (k+1)<∞,
  
  (Z_k)_k∈ℕ∈ L^2(w⃗) almost surely. From Claim (1), every sub-sequential limit of P_n has the same finite-dimensional distributions as (Z_k)_k∈ℕ. From Proposition <ref> (3),  every sub-sequential weak limit of P_n in 𝒫(𝒳)  is equal to the law of (Z_k)_k∈ℕ. 
  
  By Prokhorov's Theorem (see for example <cit.>), if {P_n}_n∈ℕ is tight, and every weakly convergent sub-sequence has the same limit μ in 𝒫(𝒳), then the sequence {P_n}_n∈ℕ converges weakly to μ. 
  Since we have already shown  every sub-sequential weak limit of  P_n is the law of (Z_k)_k∈ℕ in 𝒫(𝒳), to finish the proof, it remains to show {P_n}_n∈ℕ is tight. 
  
  From the description of compact sets in L^2(w⃗) given in Proposition <ref> (1), it suffices to show for any ε>0, there exists an element (a_k)_k∈ℕ∈ L^2(w⃗) with a_k>0, ∀ k∈ℕ, such that
  
      sup_nℙ[⋃_k∈ℕ{|N_k^(n)|>a_k}]=sup_nℙ[⋃_k=1^r_n{|N_k^(n)|>a_k}] <ε,
  
  where ⋃_k∈ℕ{|N_k^(n)|>a_k} is the complement of a compact set in L^2(w⃗).
  
  
  
  
  For any fixed ϵ>0, choose a_k=α k √(log (k+1)) for a constant α^2>32 depending on ε, then (a_k)_k∈ℕ∈ L^2(w⃗) and a_k>0, ∀ k∈ℕ.
  According to the definition of N_k^(n) in (<ref>), the above Condition (<ref>) is equivalent to 
  
      sup_nℙ[⋃_k=1^r_n{|_k^(∞)-𝔼_k^(∞)|>a_k[(d_1-1)(d_2-1)]^k/2}] <ε.
  
  
  From the proof of Theorem 22 in <cit.>, _k^(∞), as a sum of independent Poisson random variables, satisfies the following concentration inequality: for any t>0,
  
      ℙ( |_k^(∞)-𝔼_k^(∞)|>t)≤ 2exp(-t/8klog(1+t/2k[(d_1-1)(d_2-1)]^k) ).
  
   
  
  Since log (1+x)≥ x/2 for x∈ [0,1], we have from (<ref>), for sufficiently large n and  all k≤ r_n,
  
      ℙ( |_k^(∞)-𝔼_k^(∞)|>a_k[(d_1-1)(d_2-1)]^k/2)
      ≤   2exp(-a_k[(d_1-1)(d_2-1)]^k/2/8klog(1+a_k/2k[(d_1-1)(d_2-1)]^k/2) ) 
      ≤   2 exp(-a_k^2/32k^2)=2(k+1)^-α^2/32.
  
  With the assumption α^2>32, we can  make
  ∑_k=1^∞2(k+1)^-α^2/32<ϵ by choosing a sufficiently large constant α depending on ε, which guarantees (<ref>). Hence {P_n}_n∈ℕ is tight. This completes the proof. 
  
  
  
  We now continue to study the eigenvalue fluctuation for XX^⊤-d_1I/√((d_1-1)(d_2-1)) when d_1· d_2→∞. Before stating the main result, we make several assumptions on the test function f.
  Define 
  
      Φ_0(x)    =1,  Φ_k(x)=2T_k(x/2),   ∀ k≥ 1.
  
  Assume f is an entire function on ℂ.  Let K_1=max{α_1,α_2 }, where α_1 and α_2 are the constants in (<ref>), (<ref>), respectively.
  Then from Proposition <ref>, f has the expansion
  
      f(x)=∑_i=0^∞ a_i Φ_i(x)
  
  on [-K_1,K_1]. Denote 
  
      f_k(x):=∑_i=0^k a_iΦ_i(x).
  
  
  
  Suppose the following conditions hold for f:
  
      
    * For some α>3/2 and M>0, 
      
      sup_|x|≤ K_1|f(x)-f_k(x)|≤ Mexp(-α k h(k)),
  
      where h is a function such that h(r_n)≥log [(d_1-1)(d_2-1)] for a sequence 
       
      r_n    =⌊βlog n/log[(d_1-1)(d_2-1)]⌋
  
  with a constant β<1/α.
  
    * 
      lim_n→∞| f_r_n(d_1(d_2-1)/√((d_1-1)(d_2-1)))-f(d_1(d_2-1)/√((d_1-1)(d_2-1)))|=0.
  
  
  
  
  Let 
  μ_k(d_1,d_2):=𝔼_k^(∞).
  We define the following sequence: 
  
      m_f^(n):=na_0+∑_k=1^r_na_k/[(d_1-1)(d_2-1)]^k/2(μ_k(d_1,d_2)-n(d_1-2)· (d_2-1)^k/21{k is even}).
  
  
  
  
  Now we are ready to state our results for  eigenvalue fluctuations when d_1· d_2→∞. Here d_1,d_2, (λ_i)_1≤ i≤ n and the matrix X are quantities  depending on n, but for simplicity of notations, we drop the dependence on n. 
  
  
  Let G_n be a sequence of random (d_1,d_2)-biregular bipartite graphs with 
  
      d_1d_2→∞,    d_1d_2=n^o(1).
   Let λ_1≥⋯≥λ_n be the eigenvalues of XX^⊤-d_1I/√((d_1-1)(d_2-1)). 
  Suppose one of the following two assumptions holds:
  
      
    * There exists a constant c≥ 1 such that 1≤d_1/d_2≤ c.
      
    * There exists a constant c_1 such that d_2≤ c_1 for all n.
   
  Let f be an entire function on ℂ satisfying (<ref>) and (<ref>).
  Then as n→∞, the random variable
  
      Y_f^(n)=∑_i=1^n f(λ_i)- m_f^(n)
  
  converges in distribution to a centered Gaussian random variable with  variance σ_f=2∑_k=2^∞ka_k^2.
  
  Moreover, for any fixed t, consider the entire functions g_1,…, g_t   satisfying (<ref>) and (<ref>). The corresponding random vector 
  (Y_g_1^(n),…, Y_g_t^(n)) converges in distribution to a centered Gaussian random vector (Z_g_1,…, Z_g_t) with covariance
  
      Cov(Z_g_i,Z_g_j)=2∑_k=2^∞ka_k(g_i) a_k(g_j)
  
  for 1≤ i,j≤ t,
  where a_k(g_i), a_k(g_j) are the k-th coefficients in the expansion (<ref>) for g_i,g_j, respectively.
  
  
  
  
  
  We first prove the CLT for a single test function f.
  Define 
  
      X_f^(n)   :=∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(n) -𝔼∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(∞) ,
      X̃_f^(n)   :=∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(∞) -𝔼∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2_k^(∞).
  
  
  
  Recall the definition of m_f(n) in (<ref>). 
  From  (<ref>), (<ref>), and  (<ref>), X_f^(n) can be written as
  
      X_f^(n)   =∑_k=2^r_n∑_i=1^n a_kΓ_k(λ_i)-∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2μ_k(d_1,d_2)
         =∑_k=2^r_n∑_i=1^n (2a_k T_k(λ_i/2)+a_k(d_1-2)/(d_1-1)^k/21_{k is even})-∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2μ_k(d_1,d_2)
         =∑_i=1^n f_r_n(λ_i)-na_0+ ∑_k=2^r_nna_k(d_1-2)/(d_1-1)^k/21_{k is even}-∑_k=2^r_na_k/[(d_1-1)(d_2-1)]^k/2μ_k(d_1,d_2)
         =∑_i=1^n f_r_n(λ_i)-m_f^(n),
  
  where in the third line we use the fact given in (<ref>) that 
  
      ∑_i=1^n 2a_1T_1(λ_i/2)=∑_i=1^n a_1Γ_1(λ_i)=a_1[(d_1-1)(d_2-1)]^-1/2_1^(n)=0.
  
  From the definition of N_k^(n) in (<ref>),
  
      X̃_f^(n)   =∑_k=2^r_na_kN_k^(n).
  
  By Lemma <ref> and Proposition <ref> (2), X̃_f^(n) converges in distribution to a centered Gaussian random variable with variance  σ_f=∑_k=2^∞2ka_k^2. 
  
  From Corollary <ref>, the total variation distance between X_f^(n) and X̃_f^(n) satisfies 
  
      d_TV(X_f^(n), X̃_f^(n))   ≤ 	d_TV( (_k^(n),2≤ k≤ r_n),(_k^(∞),2≤ k≤ r_n))
         ≤c_8√(r_n)[(d_1-1)(d_2-1)]^3r_n/2/nd_1,
  
  which converges to 0 as n→∞ from the assumption (<ref>). Therefore X_f^(n) and  X̃_f^(n) converge to the same limit.
  
  It remains to show  Y_f^(n) and X_f^(n) converge in distribution to the same limit. We have f_r_n(λ_1)→ f(λ_1) as n→∞ from (<ref>). Then for any δ>0, 
  |f(λ_1)-f_r_n(λ_1)|≤δ/2 for sufficiently large n.
  
  Suppose that all the non-trivial eigenvalues are contained in [-K_1,K_1]. From Condition (<ref>),  we have for sufficiently large n,
  
      |Y_f^(n)-X_f^(n)|   ≤∑_i=1^n |f(λ_i)-f_r_n(λ_i)|≤δ/2+ (n-1)Mexp(-α r_nh(r_n))≤δ/2+Mn^1-αβ<δ.
  
  Therefore 
  
      ℙ( |Y_f^(n)-X_f^(n)|≥δ)≤ℙ( max_2≤ i≤ n|λ_i|≥ K_1 )=o(1),
  
  where the last inequality is from parts (2) and (3) in Theorem <ref>. Hence Y_f^(n) and X_f^(n) converge in distribution to the same limit.  This proves the CLT for (<ref>).
  
  
  
  We now extend the results to a  random vector (Y_g_1,…, Y_g_t). By Lemma <ref> and part (2) in Proposition <ref>, the random vector (X̃_g_1^(n),…, X̃_g_t^(n)) converges in distribution to the Gaussian random  vector (Z_g_1,…, Z_g_t) with covariance given in  (<ref>).
   
  Note that each entry in the vector (X_g_1^(n),…, X_g_t^(n) ) is a measurable function of (_k^(n) )_2≤ k≤ r_n, and we can find a measurable map ψ   such that 
  
      ψ( (_k^(n) )_2≤ k≤ r_n)   =(X_g_1^(n),…, X_g_t^(n)),   ψ( (_k^(∞) )_2≤ k≤ r_n)=(X̃_g_1^(n),…, X̃_g_t^(n)).
  
  Since any measurable map reduces the total variation distance between two random variables, we obtain
  from (<ref>),
  
      d_TV( (X_g_1^(n),…, X_g_t^(n)),(X̃_g_1^(n),…, X̃_g_t^(n)))
      ≤    d_TV( (_k^(n),2≤ k≤ r_n),(_k^(∞),2≤ k≤ r_n))
      ≤   c_8√(r_n)[(d_1-1)(d_2-1)]^3r_n/2/nd_1=o(1).
  
  Therefore (X_g_1^(n),…, X_g_t^(n)) converges in distribution to (Z_g_1,…, Z_g_t). 
  Finally, according to (<ref>), (X_g_1^(n),…, X_g_t^(n)) and (Y_g_1^(n),…, Y_g_t^(n)) converge in distribution to the same limit. This finishes the proof.
   
  
  
   In <cit.>, the authors proved a CLT for linear spectral statistics for normalized sample covariance matrices A=1/√(np)(XX^⊤-pI), where p/n→∞ and X=(X_ij)_n× p has i.i.d. entries with mean 0 variance 1. It is shown in Theorem 1 of <cit.> that the fluctuations of linear statistics for two analytic functions g_1,g_2 converge in distribution to a centered Gaussian vector with covariance given by (ν_4-3)a_1(g_1)a_1(g_2)+2∑_k=1^∞ka_k(g_1)a_k(g_2),
   where  ν_4=𝔼X_11^4.  
  The covariance  given in (<ref>) is  the same, except for the fact that the  coefficient in front of a_1(f_1)a_1(f_2) is 0. This can be explained by the fact that the number of 2-cycles  is 0 in RBBGs, whereas in the model used in <cit.> it is not.
   The same phenomenon was also observed in uniform  random regular graphs <cit.>, where the limiting variance  is the same as the eigenvalue fluctuations for the GOE  except for the first two terms, see Remark 22 in <cit.>.
  
  
  
  
  
  § GLOBAL SEMICIRCLE LAW
  
  
  Consider a random (n,m,d_1,d_2)-biregular bipartite graph with d_1≥ d_2. We assume d_1,d_2 satisfy the following: 
  
      lim_n→∞d_1    =∞,
      
      	d_1    =o(n^ϵ),   ∀ϵ>0, 
      d_1/d_2   →∞.
  
  Here d_2 can be fixed or a parameter depending on n. 
   In this section, we prove a semicircle law for the matrix XX^⊤-d_1I/√((d_1-1)(d_2-1)) under the assumptions (<ref>)-(<ref>).
   
   For RBBGs in this regime, we have the locally tree-like structure in the following sense. 
  Let R be fixed and τ_1 be the set of vertices in V_1 without any cycles in the R-neighborhood. The following lemma holds.
  
  	Then under Condition
  	(<ref>),
  	
      ℙ( n-|τ_1|/n>n^-1/4)=o(n^-5/4).
  
  
  To prove Lemma <ref>, the following estimates on the expectation and variance of the cycle counts  of RBBGs given in <cit.> are needed.
   Let C_k be the number of cycles of length 2k in a random  (d_1,d_2)-biregular bipartite graph. Denote μ_k=[(d_1-1)(d_2-1)]^k/2k. If d_1=o(n), k=O(log n) and kd_1=o(n), then 
      
      𝔼C_k   =μ_k(1+O( k(k+d_1)/n)),
      Var[C_k]   =μ_k(1+O( d_1^2k(k(d_1/d_2)^2k-1+(d_1/d_2)^-kd_2)/n)).
  
  
  
  
   
  
  From Lemma <ref>, for each fixed k, under Condition (<ref>),
  
      𝔼C_k=(1+o(1))μ_k,  Var[C_k]=(1+o(1))μ_k.
   
  
  If a vertex v_1∈ V_1 is not in τ_1, then for some s with 2≤ s≤ R, there exists a 2s-cycle within (R-s)-neighborhood of v_1. Hence the size of all (R-s)-neighborhoods of 2s-cycles from V_1 gives an upper bound on (n-|τ_1|). 
  
  For any 2s-cycle, the size of its (R-s)-neighborhood from V_1 is bounded by 
      c_1s[(d_1-1)(d_2-1)]^(R-s)/2+1
   with an absolute constant c_1. Define
  
      N_R:=c_1∑_s=2^Rs[(d_1-1)(d_2-1)]^(R-s)/2+1C_s.
  
  We then have n-τ_1≤ N_R. From (<ref>),
  𝔼N_R=O([(d_1-1)(d_2-1)]^R+1).
  
  Recall R is fixed. By Cauchy inequality,
  
      Var[N_R]≤ c_1^2R∑_s=2^Rs^2[(d_1-1)(d_2-1)]^R-s+2Var[C_s]=O([(d_1-1)(d_2-1)]^R+2).
  
  Then from Markov's inequality, together with our assumptions (<ref>)-(<ref>),
  
      ℙ(n-|τ_1|/n>n^-1/4)   =ℙ( n-|τ_1| >n^3/4) ≤ℙ(N_R≥ n^3/4)
         ≤𝔼[N_R^2]/n^3/2=O([(d_1-1)(d_2-1)]^2R+2n^-3/2)=o(n^-5/4).
  
  
  
  We now state our main result in this section. The proof is based on the moment method and the tree approximation of local neighborhoods, which were previously applied to random regular graphs  in <cit.>. 
  
  Let G_n be a sequence of  random (d_1,d_2)-biregular bipartite graph. Under assumptions (<ref>)-(<ref>),  the empirical spectral distribution of XX^⊤-d_1I/√((d_1-1)(d_2-1)) converges weakly to the semicircle law almost surely.
  
  
  Recall in <cit.>, when the ratio d_1/d_2≥ 1 converges to a positive constant, the ESD of XX^⊤/d_1 converges to Marčenko-Pastur law. With different scaling parameters, we obtain a different semicircle law when d_1/d_2→∞. This can be seen as an analog of the semicircle law for sample covariance matrices proved in <cit.> when the aspect ratio is unbounded.
  
  
  
  Note that for all i∈ V_1, by the degree constraint,
  
      (XX^⊤)_ii=∑_jX_ijX_ji=∑_jX_ij= (i)=d_1.
  
  	Denote M=XX^⊤-d_1I/√((d_1-1) (d_2-1)). We start with the trace expansion of M. 
  	
      1/ntr M^k   =1/n((d_1-1)(d_2-1))^k/2tr(XX^⊤-d_1I)^k 
         =	1/n((d_1-1)(d_2-1))^k/2∑_ i_1,…, i_k∈ [n]
       i_1≠i_2,…, i_k≠i_1
       j_1,…, j_k∈ [m]X_i_1j_1X_i_2j_1⋯ X_i_kj_kX_i_1j_k.
  
  	From (<ref>), the diagonal entries of XX^T-d_1I are 0, therefore we have the constraint that i_1≠i_2,…, i_k≠i_1 in (<ref>).
  	
  	Let A_k^r,c(v,v) be the number  of all
  	closed  walks of  length 2k in G starting from v∈ V_1 that use r distinct vertices from V_1, c distinct vertices from V_2, with the restriction that i_1≠i_2,…, i_k≠i_1. We have  r≤ k+1 and c≤ k, since there are at most  k+1 vertices in V_1 and k vertices in V_2 that are visited in one closed walk of length 2k.  From (<ref>), the k-th moment of the empirical spectral distribution μ_n  satisfies
  
      ∫ x^k dμ_n(x)=1/ntr M^k=1/n((d_1-1)(d_2-1))^k/2∑_v∈ V_1∑_r=1^k+1∑_c=1^k A_k^r,c(v,v).
  
  Since d_1≥ d_2, for any fixed v∈ V_1, we have  
  
      ∑_r≤ k+1,c≤ k A_k^r,c(v,v)≤ d_1^2k.
  
  For ease of notation,  in the following equations, we often omit the range of r,c in the summation.
  
  
  We may decompose the sum in (<ref>) into two parts depending on whether v∈τ_1 or not.
  For any  v∈τ_1, we write
  A_k^r,c=:A_k^r,c(v,v) since all neighborhoods of v∈τ_1 of radius k looks the same and the number of such closed walks is independent of v. 
  Now we have the following upper bound on (<ref>):  
  
      ∫ x^k dμ_n(x)    ≤1/n((d_1-1)(d_2-1))^k/2∑_v∈τ_1∑_r,cA_k^r,c(v,v)+ (n-|τ_1|) d_1^2k/n((d_1-1)(d_2-1))^k/2
         =|τ_1|/n((d_1-1)(d_2-1))^k/2∑_r,cA_k^r,c+(n-|τ_1|) d_1^2k/n((d_1-1)(d_2-1))^k/2
         ≤1/((d_1-1)(d_2-1))^k/2∑_r,cA_k^r,c+(n-|τ_1|) d_1^2k/n((d_1-1)(d_2-1))^k/2.
  
  Similarly, a lower bound holds by only counting closed walks starting with vertices in τ_1:
  
      ∫ x^k dμ_n(x)≥   1/n((d_1-1)(d_2-1))^k/2∑_v∈τ_1∑_r∑_c A_k^r,c(v,v) 
      	= |τ_1|/n((d_1-1)(d_2-1))^k/2∑_r∑_c A_k^r,c.
  
  
  From Lemma <ref> and assumption (<ref>), with probability at least 1-o(n^-5/4), for any fixed k≥ 0,
  
      (n-|τ_1|)/n((d_1-1)(d_2-1))^k/2 d_1^2k=o(1),   and  |τ_1|/n=1-o(n^-1/4).
   To show the almost sure convergence of the empirical measure to semicircle law, by the upper and lower bounds above, it suffices to show 
  
      lim_n→∞1/((d_1-1)(d_2-1))^k/2∑_r,c A_k^r,c =
      	0     if k is odd,
      
      	C_k/2    if k is even,
  
  where C_k:=1/k+12kk is the k-th Catalan number.
  
  Recall A_k^r,c counts the closed walks of length 2k on a rooted (d_1,d_2)-biregular tree starting from a root with degree d_1, ending at the same  root. Now we consider the quantity 
  
  
      1/((d_1-1)(d_2-1))^k/2∑_r,c A_k^r,c
  
    more carefully. We first consider possible ranges of r and c in the expression above. 
    
    The walk (i_1,j_1,i_2,j_2,…,i_k,j_k,i_1)  in the summation satisfies i_1≠i_2,⋯,i_k-1≠i_k, i_k≠i_1.
  This implies when a walk goes from i_t to j_t for some t, it cannot backtrack immediately to i_t. Namely, any such walk is not allowed to backtrack at even depths (here, we define the depth of the root in a tree as 1). To have a closed walk of length 2k on a tree, each edge is repeated at least twice, so the number of distinct edges is at most k. Therefore the number of distinct vertices satisfies 
  
      r+c≤ k+1.
    
   
   For fixed r and c, the number of such unlabeled rooted trees with r+c-1 distinct edges is  C_r+c-1.  Let I be the set of vertices in the odd depths of the biregular tree and J be the set of vertices in the even depths. Since the first vertex of the walk is fixed (we always start from the fixed root), for any closed walk, there are at most d_1^c many ways to choose distinct vertices from J and d_2^r-1 many ways to choose distinct vertices from I.
   Therefore we have
  
      A_k^r,c≤ d_1^c d_2^r-1C_r+c-1≤ d_1^cd_2^r-1C_k,
  
  where the last inequality is from (<ref>).
  We also know that r-1≥ c, because whenever  a new vertex in J is reached by the walk, the walk cannot backtrack, so it must reach a new vertex in I. Therefore we have 
  
      c≤ r-1   and     r+c≤ k+1,
   which implies the following conditions on c and r:
  
      c≤ k/2   and    r-1≤ k-c.
  
  
  From (<ref>), for any (r,c) satisfying (<ref>), the following holds:
  
      A_k^r,c/((d_1-1)(d_2-1))^k/2≤d_1^c/(d_1-1)^k/2d_2^r-1/(d_2-1)^k/2C_k≤d_1^c/(d_1-1)^k/2d_2^k-c/(d_2-1)^k/2C_k.
  
  	
  	
   Now we discuss two cases depending on the parity of k.
   When k is odd, from (<ref>), c≤k-1/2.  
   Since d_1/d_2→∞, we obtain 
   
      1/((d_1-1)(d_2-1))^k/2∑_r,c A_k^r,c≤   (d_1/d_2)^cd_2^kC_k/[(d_1-1)(d_2-1)]^k/2≤(d_1/d_2)^c-k/2 2^k C_k=o(1).
  
  
   When k is even, to have a non-vanishing term in the limit for A_k^r,c, we must have c=k/2 and r=k/2+1. Then we have 
   
      1/((d_1-1)(d_2-1))^k/2∑_r,c A_k^r,c=1/((d_1-1)(d_2-1))^k/2	A_k^k/2+1, k/2+o(1).
   
  
  We  continue our proof with a  more refined estimate on A_k^k/2+1,k/2. Since every  edge is repeated exactly twice in the closed walk, it's a depth-first search on the biregular tree. 
  
  If the root is at level 1, and subsequent vertices are at a level i+1 where i is the distance from the root, then all leaves must be at odd levels, since we can never backtrack at an even level. This implies that every vertex at an even level has at least one child, which means r ≥ c+1, with equality if and only if every vertex at an even level has exactly one child. Thus, one can see the tree as a subdivision of a smaller tree, where a vertex has been introduced on each edge (the “new" vertices being the vertices on an even level in the bigger tree). This is a bijection between the kind of planar rooted tree on k+1 vertices we are trying to count and the set of all planar rooted trees on k/2+1 vertices. There are C_k/2 of the latter. See Figure  <ref> for an example of a valid closed walk and an illustration of the aforementioned bijection.
  
  
  
  Moreover, given a fixed root with a  vertex label, the number of all possible ways to label the tree with vertices in a biregular bipartite graph is between d_1^k/2(d_2-1)^k/2 and (d_1-k/2)^k/2 (d_2-1)^k/2, so  the following inequality for A_k^k/2+1, k/2 holds:
  
      (d_1-k/2-1)^k/2 (d_2-1)^k/2C_k/2≤ A_k^k/2+1, k/2≤ d_1^k/2(d_2-1)^k/2C_k/2.
  
  From (<ref>) and (<ref>), we obtain for even k,
  
      lim_n→∞1/((d_1-1)(d_2-1))^k/2∑_r,c A_k^r,c=C_k/2.
  
   With (<ref>) and (<ref>),  the asymptotic behavior of moments given in  (<ref>) holds. This completes the proof of Theorem <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § RANDOM REGULAR HYPERGRAPHS
  
  
  We first include some definitions for  hypergraphs and describe a bijection between a subset of biregular bipartite graphs and the set of regular hypergraphs studied in <cit.>.
  We will use the map given in Definition <ref> to apply some of our results for RBBGs to random regular hypergraphs, see <cit.> for more details. 
  
  
  	A hypergraph H consists of a set V of vertices and a set E of hyperedges such that each hyperedge is a nonempty set of V.  A hypergraph H is k-uniform for an integer k≥ 2 if every hyperedge e∈ E contains exactly k vertices. The degree of i, denoted   (i), is the number of all hyperedges incident to i.
   A hypergraph is d-regular if all of its vertices have degree d.
  A hypergraph is (d,k)-regular if it is both d-regular and k-uniform. 
  
  
  
   A vertex i is incident to a hyperedge e if and only v is an element of e. We can define the incidence matrix X of a hypergraph H=(V,E) to be a |V|× |E| matrix indexed by elements in V and E such that X_i,e=1 if i∈ e and 0 otherwise. Moreover, if we regard X as the adjacency matrix of a graph, it defines a bipartite graph G with two vertex sets  V and E. We call G  the  bipartite graph associated to H, given by a map Φ (so Φ(H)=G).
  See Figure <ref> for an example.
   
   
   
  
  
  For  a hypergraph H with n vertices, we associate a n× n symmetric  matrix A called the adjacency matrix of H. For i≠j, we define A_ij as the number of hyperedges containing both i and j; we define A_ii=0 for all 1≤ i≤ n. When the hypergraph is 2-uniform (i.e., it is a graph), this is the usual definition for the adjacency matrix of a graph.
  
   	
  The following lemma connects the adjacency matrix of a regular hypergraph with its associated biregular bipartite graph. It formally appears in <cit.>, and it is also informally mentioned in <cit.>.
  
   
  Let H be a (d_1,d_2)-regular hypergraph, and let G be the corresponding (d_1,d_2)-biregular bipartite graph. Let A_H be the adjacency matrix of H and A_G be the adjacency matrix of G given by
  
      A_G=[   0   X; X^⊤   0 ].
   Then
  	A_H=XX^⊤-d_1I.
  
  
  
  
  A walk of length l on a hypergraph H is a vertex-hyperedge sequence (i_0,e_1,i_1,⋯ ,e_l, i_l) such that i_j-1≠i_j and {i_j-1,i_j}⊂ e_j for all 1≤ j≤ l. A walk is closed if i_0=i_l. A cycle of length l in a hypergraph H is a closed walk (v_0,e_1,…, v_l-1, e_l,v_l+1) such that 
  all edges are distinct and all vertices are distinct subject to v_l+1=v_0. In the associated bipartite graph G, a cycle of length 2l corresponds to a cycle of length l in H.
  
  
  
  Let 𝒢(n,m,d_1,d_2) be the set of all  simple biregular bipartite random graphs  with vertex set V=V_1∪ V_2 such that  |V_1|=n, |V_2|=m, and every vertex in V_i has degree d_i for i=1,2.  Without loss of generality, we assume d_1≥ d_2. Let ℋ (n,d_1,d_2) be the set of all simple (without multiple hyperedges) (d_1,d_2)-regular hypergraphs with labeled  vertex set [n] and nd_1/d_2 many labeled hyperedges denoted by {e_1,…,e_nd_1/d_2}. 
  
  
  We  can also consider all (d_1,d_2)-regular hypergraphs with labeled vertices and unlabeled hyperedges. Since all hyperedges are distinct, any such regular hypergraph with unlabeled hyperedges corresponds to (nd_1/d_2)! regular hypergraphs with labeled hyperedges. 
  
  
  
  
  It is well known (see for example <cit.>) that the map Φ defined can be extended to a bijection Φ̃ between labeled regular multi-hypergraphs and biregular bipartite graphs.  See  Figure <ref> as an example of the bijection. For a given biregular bipartite graph, if there are two vertices in V_2 that have the same set of neighbors in V_1, the corresponding regular hypergraph will have  multiple hyperedges, see Figure <ref>.  Let 𝒢'(n,m,d_1,d_2) be a subset of 𝒢(n,m,d_1,d_2) such that for any G∈𝒢'(n,m,d_1,d_2), any two vertices in V_2 have different  neighborhoods in V_1. The following lemma holds.
  
  
  
  	Φ is the restriction of the bijection Φ̃ to  ℋ(n,d_1,d_2)  and its image is 𝒢'(n, m,d_1,d_2). Hence |ℋ(n,d_1,d_2)|=|𝒢'(n, m,d_1,d_2)|.
  
   From Lemma <ref>, the uniform distribution on 𝒢' (n,m,d_1,d_2) for biregular bipartite graphs  induces the uniform distribution  on ℋ (n, d_1,d_2) for regular hypergraphs. With this observation, we can translate some of the results for spectra of random biregular bipartite graphs into results for spectra of random regular hypergraphs. A similar approach was applied in <cit.> to enumerate uniform hypergraphs with given degrees.
   
   
  
  
  
  
  
  Let G be a random biregular bipartite graph sampled uniformly from 𝒢(n,m,d_1,d_2) such that 3≤ d_2≤ d_1≤n/32.  Let 𝒢'(n,m,d_1,d_2) be the set of biregular bipartite graphs corresponding to simple regular hypergraphs.
  Then   
  
      ℙ(G∈𝒢'(n,m,d_1,d_2) )≥ 1-(nd_1/d_2)^2(4ed_2/n)^d_2.
   
  In particular, 
  
      ℙ(G∈𝒢'(n,m,d_1,d_2) )=1-O( d_1^2/nd_2^2).
  
  
  
  
  
  Lemma <ref> implies the following total variation bound.
  
  Let μ_n be the probability measure of the random (d_1,d_2)-regular hypergraph with n vertices induced on the set of all (n,m,d_1,d_2)-biregular bipartite graphs,  and let μ_n' be the uniform measure on the set of all (n,m,d_1,d_2)-biregular bipartite graphs. We have 
  
      d_TV(μ_n,μ_n')≤(nd_1/d_2)^2(4ed_2/n)^d_2.
  
  
   Since 𝒢'(n,m,d_1,d_2) is the set of all biregular bipartite graphs that are bijective to regular hypergraphs. We have 
  μ_n(𝒢'(n,m,d_1,d_2) )=1 and μ_n'(𝒢'(n,m,d_1,d_2))=|G'(n,m,d_1,d_2)|/|𝒢(n,m, d_1,d_2)|.  Let ℱ be the power set of 𝒢(n,m,d_1,d_2).
  Taking into account the fact that both μ_n and μ_n' are uniform measures, we obtain that
  
      d_TV(μ_n,μ_n')   =sup_A∈ℱ|μ_n(A)-μ_n'(A)|=|1-μ_n'(G'(n,m,d_1,d_2)) |
         =ℙ(G∉𝒢'(n,m,d_1,d_2) )≤(nd_1/d_2)^2(4ed_2/n)^d_2,
  
  where the last inequality is from  Lemma <ref>.
  
  Equipped with Lemma <ref>, we obtain several corollaries for random regular hypergraphs in the following subsections.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Cycle counts
  
  Recall the definition of cycles in a hypergraph given in Definition <ref>.
  Let C_k be the number of cycles of length k in a (d_1,d_2)-regular hypergraph. The following result holds.
  
  Let H be a (d_1,d_2)-random regular  hypergraph with cycle counts (C_k, k≥ 2). Let (Z_k, k≥ 2) be independent Poisson random variables with 𝔼Z_k=(d_1-1)^k(d_2-1)^k/2k.	 For any n,m≥ 1, r≥ 3, and  3≤ d_2≤ d_1≤n/32,
  
      d_TV((C_2,…,C_r),(Z_2,…,Z_r))≤c_6√(r)(d_1-1)^3r/2(d_2-1)^3r/2/nd_1+(nd_1/d_2)^2(4ed_2/n)^d_2.
  
  
  
  Let C̃_k be the number of cycles with length 2k in a uniform random (d_1,d_2)-biregular bipartite graph. From Lemma <ref>,
  
      d_TV((C_2,…,C_r),(C̃_2,…,C̃_r))≤ d_TV(μ_n,μ_n')≤(nd_1/d_2)^2(4ed_2/n)^d_2.
  
  Then the conclusion follows from Theorem <ref> and the triangle inequality.
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Global laws
   The limiting spectral distributions for the adjacency matrix of a random regular hypergraph can be summarized in the following corollary.
  
  Let H be a random (d_1,d_2)-regular hypergraph. 
  
  
    * If d_1,d_2 are fixed,   the empirical spectral distribution of A-(d_2-2)/√((d_1-1)(d_2-1)) converges in probability to a measure μ with density function given by 
      f(x):=1+d_2-1/q/(1+1/q-x/√(q))(1+(d_2-1)^2/q+(d_2-1)x/√(q))1/π√(1-x^2/4) dx,
  
  	where q=(d_1-1)(d_2-1).
  	 
    * For d_1,d_2→∞ with d_1/d_2→α≥ 1 and d_1≤n/32,  the empirical spectral distribution of A-(d_2-2)/√((d_1-1)(d_2-1)) converges in probability to  a measure supported on [-2,2] with a density function given by	
  	 
      g(x)=α/1+α+√(α)x1/π√(1-x^2/4) .
  
      
    * If d_1→∞, d_1=o(n^ϵ) for any ϵ>0 and d_1/d_2→∞,  the ESD of A/√((d_1-1)(d_2-1)) converges to the semicircle law in probability.
  
  
  
  
  Claim (1) is proved in Theorem 6.4 of <cit.> based on a result for deterministic regular hypergraphs in Theorem 5 of <cit.>. 
  
  
  
  Claim (2) is a combination of several results. When d_1=o(n^1/2), it is proved in Theorem 6.6 of <cit.> based on the global law for random biregular bipartite graphs in <cit.> and <cit.>.  When d_1=ω(log^4 n), the optimal local law for RBBGs  was recently  proved in <cit.>, which also implies the global law for RBBGs. 
  When d_1≤n/32 and d_1/d_2→α, from Lemma <ref>, 
  
      ℙ(G∈𝒢'(n,m,d_1,d_2) )→ 1.
  
  Therefore by the same proof of Theorem 6.6 in <cit.>, the ESD of A-(d_2-2)/√((d_1-1)(d_2-1)) for random regular hypergraphs converges in probability.
  
  
  
  
  Under the assumptions d_1→∞, d_1=o(n^ϵ) for any ϵ>0 and d_1/d_2→∞, from  Lemma <ref>, we have again
  ℙ(G∈𝒢'(n,m,d_1,d_2) )→ 1.
  Then Claim (3) follows from  Theorem <ref> and Lemma <ref>. 
  
  
  
  
  The ESD in Corollary <ref> (2) is a shifted and scaled Marčenko-Pastur law. Taking α→∞, g(x) converges to the density function of the semicircle law. The  transition from Marčenko-Pastur law to the semicircle law was also proved for sample covariance matrices in <cit.> when the aspect ratio goes to infinity.
  
  
  A semicircle law for the adjacency matrix of d_2-uniform Erdős-Rényi random hypergraphs with growing expected degrees was proved in Theorem 5 of <cit.> when d_2 is a constant. Part (3) of Corollary <ref> proves a corresponding semicircle law for random d_2-uniform d_1-regular hypergraphs where d_2 can be a parameter depending on n.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  	
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Spectral gaps
  
  The spectral gap for random regular hypergraphs  with fixed d_1,d_2 was studied in <cit.>. Here we include the results for the case when  d_1,d_2 are growing with n.
  
  
   
  
  Let	H be a random (d_1,d_2)-regular hypergraph with d_1≥ d_2. Let λ_1≥⋯≥λ_n be the eigenvalues of A. Let λ=max_2≤ i≤ n|λ_i|.
  
  	
    * Suppose d_1≥ d_2≥ 3 is fixed. There exists a sequence ϵ_n→ 0 such that
  	 
      ℙ( |λ-(d_2-2)|≥ 2√((d_1-1)(d_2-1))+ϵ_n)→ 0
   
  	as n→∞.
  	
    * Suppose 3≤ d_2≤1/2n^2/3, d_1≥ d_2≥ cd_1 for some c∈ (0,1). Then for some constant K>0 depending on c, for all n≥ 1,
  	
      ℙ( λ≥ K√((d_1-1)(d_2-1)))=O(1/n).
   
  
    * Suppose 3≤ d_2≤ C_1 for a constant C_1, and d_1=o(n^1/2). There exists a constant C depending on C_1 such that
  	
      ℙ( λ≥ C√((d_1-1)(d_2-1)))=O(d_1^2/n^2).
   
  
  
  
  Claim (1) is proved in Theorem 4.3 in <cit.>. Claim (2) and (3) follow from part (2) and (3) in Theorem <ref> with Lemma <ref>.
  
  
  
  Results in <cit.> that Claim (2) and (3) are based on have stronger probability estimates. However, Lemma <ref> we used here yields a weaker failure probability. 
  
  
  
  
   §.§ Eigenvalue fluctuations
  
  
  The following eigenvalue fluctuation results for random regular hypergraphs can be derived from  Lemma <ref>, Lemma <ref>, and the eigenvalue fluctuations results for random biregular bipartite graphs in Section <ref>.
  
  
  
  For fixed d_1≥ d_2≥ 3, let H be a random (d_1,d_2)-regular hypergraph with adjacency matrix A. Let λ_1≥⋯≥λ_n be the eigenvalues of A/√((d_1-1)(d_2-1)). Suppose f is a function satisfying the same conditions in Theorem <ref>. Then  Y_f^(n):=∑_i=1^n f(λ_i)-n a_0 converges in distribution as n→∞ to the infinitely divisible random variable
  
      Y_f:=∑_k=2^∞a_k/[(d_1-1)(d_2-1)]^k/2_k^(∞),
  
  where _k^(∞) is defined in (<ref>).
  
  
  
  
  Let Ỹ_f^(n) be the corresponding random variable of Y_f^(n) for the uniform random biregular bipartite graphs  considered in Theorem <ref>. From the total variation distance bound in Lemma <ref>, we have 
  
      d_TV(Y_f^(n), Ỹ_f^(n))≤  d_TV(μ_n,μ_n')≤(nd_1/d_2)^2(4ed_2/n)^d_2=o(1).
  
  Therefore Ỹ_f^(n) and Y_f^(n) converge in distribution to the same law.
  
  
  
  
  
  
  Let H be a random (d_1,d_2)-regular hypergraph with d_1d_2→∞ as n→∞ and d_1d_2=n^o(1). Let λ_1≥⋯≥λ_n be the eigenvalues of A/√((d_1-1)(d_2-1)). Let f be a function satisfying (<ref>) and (<ref>).
  Suppose one of the two assumptions holds:
  
      
    * there exists a constant c≥ 1 such that 1≤d_1/d_2≤ c,
      
    * 3≤ d_2≤ c_1 for a constant c_1≥ 3.
   
  Then the random variable
  
      Y_f^(n)=∑_i=1^n f(λ_i)- m_f^(n)
  
  converges in law to a Gaussian random variable with mean zero and variance σ_f=∑_k=2^∞2ka_k^2. Moreover, for any fixed t, consider the entire functions g_1,…, g_t   satisfying (<ref>) and (<ref>). The corresponding random vector 
  (Y_g_1^(n),…, Y_g_t^(n)) converges in distribution to a centered Gaussian random vector (Z_g_1,…, Z_g_t) with covariance
  
      Cov(Z_g_i,Z_g_j)=2∑_k=2^∞ka_k(g_i) a_k(g_j)
  
  for 1≤ i,j≤ t,
  where a_k(g_i), a_k(g_j) are the k-th coefficients in the expansion (<ref>) for g_i,g_j, respectively.
  
  
  
  Recall Lemma <ref> and our assumption d_1d_2=n^o(1).
  Under Case (1), we have d_2→∞ and  
  
      d_TV(μ_n,μ_n')≤(nd_1/d_2)^2(4ed_2/n)^d_2=O(n^2)(n^(-1+o(1))d_2)=o(1).
  
  Under Case (2), we have 
  
      d_TV(μ_n,μ_n')≤(nd_1/d_2)^2(4ed_2/n)^d_2=O(n^2d_1^2)( 4ec_1/n)^3=o(1).
  
  
  Then with Lemma <ref>, in both cases Y_f^(n) converges in distribution to the same limiting random variable defined in Theorem <ref>. The proof of the covariance part follows in the same way.
  
  
  plain
  
  