The asymptotic volume of diagonal subpolytopes of symmetric stochastic matrices
J. de Jong R. Wulkenhaar

Alpha Fair Coded Caching
    Apostolos Destounis^1,
		Mari Kobayashi^2, 
		Georgios Paschos ^1,
		Asma Ghorbel^2 
		

		^1 France Research Center, Huawei Technologies Co. Ltd., email: firstname.lastname@huawei.com
			
		
		
	^2Centrale-Supélec, 
		France, email: firstname.lastname@centralesupelec.fr
	

    
======================================================================================================================================================================================================================================================================================


The asymptotic volume of the polytope of symmetric stochastic matrices can be determined by asymptotic enumeration techniques as in the case of the Birkhoff polytope. These methods can be extended to polytopes of symmetric stochastic matrices with given diagonal, if this diagonal varies not too wildly. To this end, the asymptotic number of symmetric matrices with natural entries, zero diagonal and varying row sums is determined.




Keywords: Asymptotic enumeration, Polytope volumes

MSC 2010: 05A16, 52B11



§ INTRODUCTION

Convex polytopes arise naturally in various places in mathematics. A fundamental problem is the polytope's volume. Some results are known for low-dimensional setups <cit.>, polytopes with only a few vertices, or highly symmetric cases <cit.>. This work belongs to the latter category.


A convex polytope P is the convex hull of a finite set S_P={v_j∈ℝ^n} of vertices.

Stochastic matrices are square matrices with nonnegative entries, such that every row of the matrix sums to one. The symmetric stochastic N× N-matrices are an example of a convex polytope. It will be denoted by 𝒫_N. Its vertices are given by the symmetric permutation matrices. There are ∑_j=0^N/2N2j(2j-1)!! such matrices. It follows directly from the Birkhof-Von Neumann theorem that all symmetric stochastic matrices are of this form. A basis for this space is given by

    {I_N}∪{B^(jk)|1≤ j<k≤ N}  ,

where I_N is the N× N identity matrix and the matrix elements of B^(jk) are given by

    B^(jk)_lm={[         B^(jk)_lm=1 , if {l,m}={j,k}  ;;         B^(jk)_lm=1   , if j≠ l=m≠ k  ;;         B^(jk)_lm=0         , otherwise ].  .

All these vertices are linearly independent and it follows that the polytope is

. 

A convex subpolytope P' of a convex polytope P is the convex hull of a finite set {v'_j∈P} of elements in P.

Slicing a polytope yields a surface of section, which is itself a convex space and, hence, a polytope. Determining its vertices is in general very difficult.

Spaces of symmetric stochastic matrices with several diagonal entries fixed are examples of such slice subpolytopes of 𝒫_N, provided that these entries lie between zero and one. The slice subpolytope of 𝒫_N, obtained by fixing all diagonal entries h_j∈[0,1], will be called the diagonal subpolytope P_N(h_1,…,h_N) here. This is a polytope of dimension N(N-3)/2. These polytopes form the main subject of this paper.


To keep the notation light, vectors of N elements are usually written by a bold symbol. The diagonal subpolytope with entries h_1,…,h_N will thus be written by P_N(h).


The main results are the following two theorems.

Let V_N(t;λ) be the number of symmetric N× N-matrices with an empty diagonal and entries in the natural numbers such that t_j is the j-th row sum. Denote the total entry sum by x=∑_j=1^Nt_j and let λ be the average matrix entry 

    λ=x/N(N-1)>C/log N  .

If for some ω∈(0,1/4) the limit

    lim_N→∞t_j-λ(N-1)/λ N^1/2+ω=0    for all j=1,…,N  ,

then the number of such matrices is asymptotically (N→∞) given by

    V_N(t;λ)=√(2)(1+λ)^N2/(2πλ(λ+1)N)^N/2(1+1/λ)^x/2exp[14λ^2+14λ-1/12λ(λ+1)]
       ×exp[-1/2λ(λ+1)N∑_m(t_m-λ(N-1))^2]exp[-1/λ(λ+1)N^2∑_m(t_m-λ(N-1))^2]
       ×exp[2λ+1/6λ^2(λ+1)^2N^2∑_m(t_m-λ(N-1))^3]exp[-3λ^2+3λ+1/12λ^3(λ+1)^3N^3∑_m(t_m-λ(N-1))^4]
       ×exp[1/4λ^2(λ+1)^2N^4(∑_m(t_m-λ(N-1))^2)^2]×(1+𝒪(N^-1/2+6ω))  .




Let h=h_1,…,h_N with h_j∈[0,1] and χ=∑_j=1^Nh_j. If

    lim_N→∞N^1/2-ωN-1/N-χ·|h_j-χ/N|=0  for all j=1,…,N    ,

and for some ω∈(loglog N/2log N,1/4), then the asymptotic volume (N→∞) of the polytope of symmetric stochastic N× N-matrices with diagonal (h_1,…,h_N) is given by

    (P_N(h))=√(2)e^7/6(e(N-χ)/N(N-1))^N2(N(N-1)^2/2π(N-χ)^2)^N/2exp[-N(N-1)^2/2(N-χ)^2∑_j(h_j-χ/N)^2]
       ×exp[-(N-1)^2/(N-χ)^2∑_j(h_j-χ/N)^2]exp[-N(N-1)^3/3(N-χ)^3∑_j(h_j-χ/N)^3]
       ×exp[-N(N-1)^4/4(N-χ)^4∑_j(h_j-χ/N)^4]exp[(N-1)^4/4(N-χ)^4(∑_j(h_j-χ/N)^2)^2]
       ×(1+𝒪(N^-1/2+6ω))  .



The outline of this paper is as follows. In Paragraph <ref> the volume problem is formulated as a counting problem and subsequently as a contour integral. Under the assumption of a restricted region this is subsequently integrated in Paragraph <ref>. Paragraph <ref> is dedicated to a fundamental lemma to actually restrict the integration region. The volume of the diagonal subpolytopes is extracted from the counting result in Paragraph <ref>.



§ COUNTING PROBLEM

The volume of a polytope P in ℝ^n with basis {ℬ_j∈ℝ^n|1≤ j≤ d} is obtained by

    ∫_[0,1]^du 1_P(∑_j=1^du_jℬ_j)  ,

where 1_P is the indicator function for the polytope P. If the polytope is put on a lattice (aℤ)^n with lattice parameter a∈(0,1), an approximation of this volume is obtained by counting the lattice sites inside the polytope and multiplying this by the volume a^n of a single cell. This approximation becomes better as the lattice parameter shrinks. In the limit this yields

    (P)=lim_a→ 0 a^n  |{P∩ (aℤ)^n}|  .

This approach is formalized by the Ehrhart polynomial <cit.>, which counts the number of lattice sites of ℤ^n in a dilated polytope. A dilation of a polytope P by a factor a^-1>1 yields the polytope a^-1P, which is the convex hull of the dilated vertices S_a^-1P={a^-1v|v∈ S_P}. That the obtained volume is the same, follows from the observation

    |{a^-1P∩ℤ^n}|=|{P∩ (aℤ)^n}|  .


The volume integral of the diagonal subpolytope P_N(h) is

    (P_N(h))={∏_1≤ k<l≤ N∫_0^1 u_kl}  1_P_N(h)(I_N+∑_1≤ k<l≤ Nu_kl(B^(kl)-I_N))  .

To see that this integral covers the polytope, it suffices to see that the any symmetric stochastic matrix A=(a_kl) is decomposed in basis vectors as

    A=(a_kl)=I_N+∑_1≤ k<l≤ Na_kl(B^(kl)-I_N)  .


The next step is to introduce a lattice (aℤ)^N2 and count the sites inside the polytope. Each such site is a symmetric stochastic matrix with h_1,…,h_N on the diagonal.

Since the volume depends continuously on the extremal points, it can be assumed without loss of generality that all h_j are rational. This implies that a dilation factor a^-1 exists, such that all a^-1(1-h_j)=t_j∈ℕ and that the matrices that solve

    ([    0 b_12    ⋯ b_1N; b_12    0    ⋯ b_2N;    ⋮    ⋮    ⋱    ⋮; b_1N b_2N    ⋯    0 ])([ 1; 1; ⋮; 1 ])=([ t_1; t_2;   ⋮; t_N ])

with t_j,b_jk∈ℕ are to be counted. This yields a number V_N(t). The polytope volume is then given by

    (P_N(h))=lim_a→0 a^N(N-3)/2V_N(1-h_1/a,…,1-h_N/a)  ,

where 

    V_N(t)=∮_𝒞 w_1/2π i w_1^1+t_1…∮_𝒞 w_N/2π i w_N^1+t_N ∏_1≤ k<l≤ N1/1-w_kw_l  .

To see this, let the possible values m for the matrix element b_jk be given by the generating function

    1/1-w_jw_k=∑_m=0^∞(w_jw_k)^m  .

Applying this to all matrix entries shows that V_N(t) is given by the coefficient of the term w_1^t_1w_2^t_2… w_N^t_N in ∏_1≤ j<k≤ N1/1-w_jw_k. Formulating this in derivatives yields

    V_N(t)=1/t_1!/ w_1|_w_1=0^t_1…1/t_N!/ w_N|_w_N=0^t_N ∏_1≤ k<l≤ N1/1-w_kw_l  .

By Cauchy's integral formula the number of matrices (<ref>) follows from this. The contour 𝒞 encircles the origin once in the positive direction, but not the pole at w_kw_l=1.


The next step is to parametrize this contour explicitly and find a way to compute the integral for N→∞. This must be done in such a way that a combinatorial treatments is avoided. A convenient choice is

    w_j=√(λ_j/λ_j+1)e^iφ_j    , with λ_j∈ℝ_+ and φ_j∈[-π,π)  .

Later a specific value for λ_j will be chosen.

The counting problem has now been turned into an integral over the N-dimensional torus

    V_N(t)=(∏_j=1^N(1+1/λ_j)^t_j/2)(2π)^-N∫_𝕋^Nφ e^-i∑_j=1^Nφ_jt_j
       ×∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) 1/1-√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) (e^i(φ_k+φ_l)-1)  ,

where we have written φ for φ_1…φ_N.


The notations

    x=∑_jt_j=∑_j=1^Nt_j    and    ∑_k<l(φ_k+φ_l)=∑_1≤ k<l≤ N(φ_k+φ_l)

are used, when no doubt about N can exist. When no summation bounds are mentioned, these will always be 1 and N. The notation a≪ b indicates that a<b and a/b→0.


The main tool for these integrals will be the stationary phase method, also called the saddle-point method. In the form used in this paper, the exponential of a function f is integrated around its maximum x̃, so that

    lim_Λ→∞∫ x e^Λ f(x)=lim_Λ→∞exp[Λ f(x̃)]∫ x exp[Λ f^(2)(x̃)/2(x-x̃)^2+Λ f^(3)(x̃)/6(x-x̃)^3]
       ×exp[Λ f^(4)(x̃)/24(x-x̃)^4]
       =exp[Λ f(x̃)]√(-2π/Λ f^(2)(x̃))(1+15/162(f^(3)(x̃))^2/9Λ(-f^(2)(x̃))^3+3/4f^(4)(x̃)/6Λ (f^(2)(x̃))^2+𝒪(Λ^-2))  .


Many counting problems can be computed asymptotically by the saddle-point method <cit.>. Often it is assumed that all t_j are equal, but we show that it suffices to demand that they do not deviate too much from this symmetric case.



§ INTEGRATING THE CENTRAL PART


The integrals in (<ref>) are too difficult to compute in full generality. A useful approximation can be obtained from the observation that the integrand

    |1/1-μ(e^iy-1)|^2=1/1-2μ(μ+1)(cos(y)-1)    for y∈(-2π,2π)

is concentrated in a neighbourhood of the origin and the antipode y=±2π, where it takes the value 1. This is plotted in Figure <ref>. For small y and μ y the absolute value of the integrand factor can be written as

    |1/1-μ(e^iy-1)|=√(1/1+μ(μ+1)y^2)(1+𝒪(y^4))   .




It is concentrated in a small region around the origin and the antipode. The form of the region is assumed to be [δ_N,δ_N]^N with 

    δ_N=N^-αζ_N/min_j{λ_j}  ,

where α∈(0,1/2) and ζ_N tends slowly to infinity. In the remainder of this paragraph the integral inside this box will be computed.

To this end, a lower bound is introduced. Below this threshold we do not strive for accuracy. The aim is thus to find the asymptotic number V_N(t) for configurations t, such that this number is larger than the Lower bound.
Lower bound

For N, α∈(0,1/2), t_j∈ℕ and λ_j∈ℝ_+^N for j=1,…,N we define the Lower bound by

    ℰ_α = (2πλ(λ+1)N)^-N/2(∏_j=1^N(1+1/λ_j)^t_j/2)(∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       ×exp[14λ^2+14λ-1/12λ(λ+1)]exp[-N^1-2α]  ,

where λ=N^-1∑_jλ_j.


The integral in [-δ_N,δ_N]^N can now be cast into a simpler form, where the size δ_N of this box can be used as an expansion parameter. The expansion used is

    1/1-μ(exp[iy]-1)=exp[∑_j=1^kA_j(iy)^j]+𝒪(y^k+1(1+μ)^k+1)  .

The coefficients A_j(μ) (or A_j if the argument is clear) are polynomials in μ of degree j. They are obtained as the polylogarithms

    A_n(μ)=(-1)^n/n! _1-n(1+1/μ)  .

The first four coefficients are

    A_1=μ  ;   A_2=μ/2(μ+1)  ;   A_3=μ/6(μ+1)(2μ+1)  
       and   A_4=μ/24(μ+1)(6μ^2+6μ+1)  .

The value of the parameter μ in the above formules can be approximated. Assuming that ε_k is small compared to λ and writing ε=max_kε_k, this is

    √((λ+ε_k)(λ+ε_l))/√((λ+ε_k+1)(λ+ε_l+1))-√((λ+ε_k)(λ+ε_l))≈λ+ε_k+ε_l/2
       -2λ+1/8λ(λ+1)(ε_k-ε_l)^2+2λ^2+2λ+1/16λ^2(λ+1)^2(ε_k^3-ε_k^2ε_l-ε_kε_l^2+ε_l^3)+𝒪(ε^4/λ^3)  .


Applying this in combination with (<ref>) produces the combinations

    ∑_k<l(φ_k+φ_l)·(√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       =∑_j=1^Nφ_j[N-2/2λ_j+N/2λ-B_1(Nε_j^2+∑_mε_m^2)+C_1(Nε_j^3-ε_j∑_mε_m^2+∑_mε_m^3)]
       ×(1+𝒪(Nε^4/λ^4))    ;
    
       ∑_k<l(φ_k+φ_l)^2· A_2(√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       =[∑_j=1^Nφ_j^2((N-2)A_2+ε_jB_2(N-4)-(N-4)C_2ε_j^2-C_2∑_mε_m^2)
       +∑_j=1^Nφ_j(A_2∑_mφ_m+2B_2ε_j∑_mφ_m-2C_2ε_j^2∑_mφ_m+D_2ε_j∑_mε_mφ_m)]
       ×(1+𝒪(Nε^3/λ^3))    ;
       ∑_k<l(φ_k+φ_l)^3· A_3(√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       =[∑_jφ_j^3A_3(N-4)+3A_3∑_jφ_j^2∑_mφ_m]×(1+𝒪(ε/λ))    and
       ∑_k<l(φ_k+φ_l)^4· A_4(√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       =[∑_jφ_j^4A_4(N-8)+4A_4∑_jφ_j^3∑_mφ_m+3A_4(∑_jφ_j^2)^2]×(1+𝒪(ε/λ))  .


Here we used the additional combinations

    B_1 = 2λ+1/8λ(λ+1)  ;   C_1=2λ^2+2λ+1/16λ^2(λ+1)^2  ;   B_2=2λ+1/4
       C_2=2λ^2+2λ+1/16λ(λ+1)  ;   D_2=6λ^2+6λ+1/8λ(λ+1)

to simplify the notation.


The simplest way to compute this integral is to ensure that the linear part of the exponent is small.
Splitting λ_j=λ+ε_j and choosing the value 

    ε_j=2/N-2(t_j-λ(N-1))

is done therefore. Combined with the assumption that x=∑_jt_j=λ N(N-1), this implies that ∑_mε_m=0. Assuming furthermore that |t_j-λ(N-1)|≪λ N^1/2+ω, the error terms |ε/λ|≪ N^-1/2+ω follow.


The first step now is to focus on the integral inside the box [-δ_N,δ_N]^N, simplify and calculate this.

The estimates in Lemma <ref> cause the integral (<ref>) to depend non-trivially on λ. For that reason λ is explicitly mentioned as an argument.


Assume that K,N∈ℕ, ω,α∈ℝ_+ are chosen such that ω∈(0,log(Kα-2) + loglog N/4log N), α∈(0,1/4-ω) and K>2/α+1. Define

    δ_N=N^-αζ_N/min{λ_j}  ,

so that ζ_N→∞ and N^-δζ_N→0 for any δ>0, when N→∞. If x=∑_jt_j, the average matrix entry λ=x/N(N-1) and

    lim_N→∞t_j-λ(N-1)/λ N^1/2+ω=0  for j=1,…,N    ,

then the integral

    V_N(t)=(∏_j=1^N(1+1/λ_j)^t_j/2)(2π)^-N∫_[-δ_N,δ_N]^Nφ e^-i∑_j=1^Nφ_jt_j
       ×∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) 1/1-√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) (e^i(φ_k+φ_l)-1)

is given by

    V_N(t;λ)=2/(2π)^N(∏_j=1^N(1+1/λ_j)^t_j/2)·(∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       ×∫_[-δ_N,δ_N]^Nφ exp[-i∑_jφ_jt_j]exp[∑_n=1^K-1i^n∑_k<lA_n(√(λkλ_l)/√((1+λk)(1+λ_l))-√(λkλ_l))·(φ_k+φ_l)^n]+𝒟 ,

up to a difference 𝒟 that satisfies

    |𝒟|≤𝒪(N^2-Kα)√(2)(1+λ)^N2/(2πλ(λ+1)N)^N/2(1+1/λ)^x/2exp[10λ^2+10λ+1/4λ(λ+1)]
       exp[-1/2λ(λ+1)N∑_m(t_m-λ(N-1))^2]exp[3/4λ^2(λ+1)^2N^2∑_m(t_m-λ(N-1))^2]
       ×exp[2λ+1/6λ^2(λ+1)^2N^2∑_m(t_m-λ(N-1))^3]exp[6λ^2+6λ+1/24λ^3(λ+1)^3N^3∑_m(t_m-λ(N-1))^4]
       ×exp[6λ^2+6λ+1/8λ^3(λ+1)^3N^4(∑_m(t_m-λ(N-1))^2)^2]  .



To the fraction 

    (1-√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) (e^i(φ_k+φ_l)-1))^-1

in the integral (<ref>) the expansion (<ref>) in combination with (<ref>) and (<ref>) is applied. To prove that contributions in (<ref>) of K-th order or higher are irrelevant, we put these in the exponential exp[h(x)]. To estimate their contribution, the estimate

    |∫ x e^f(x)(e^h(x)-1)|≤𝒪(sup_x |e^h(x)-1|)·∫ x  |e^f(x)|

is applied to the integral. Taking the absolute value of the integrand sets the imaginary parts of the exponential to zero. In terms of (<ref>) and (<ref>) this means that A_3, B_1 and C_1 are set to zero. This integral is calculated in Lemma <ref>. Taking this result and setting these coefficients to zero completes the proof.



Assume that K,N∈ℕ, ω,α∈ℝ_+ are chosen such that ω∈(0,log(Kα-2) + loglog N/4log N), α∈(0,1/4-ω) and K>2/α+1. Define

    δ_N=N^-αζ_N/min{λ_j}  ,

so that ζ_N→∞ and N^-δζ_N→0 for any δ>0, when N→∞. If x=∑_jt_j, the average matrix entry λ=x/N(N-1)>C/log N and

    lim_N→∞t_j-λ(N-1)/λ N^1/2+ω=0  for j=1,…,N    ,

then the integral

    V_N(t;λ)=2/(2π)^N(∏_j=1^N(1+1/λ_j)^t_j/2)·(∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       ×∫_[-δ_N,δ_N]^Nφ exp[-i∑_jφ_jt_j]exp[∑_n=1^K-1i^n∑_k<lA_n(√(λkλ_l)/√((1+λk)(1+λ_l))-√(λkλ_l))·(φ_k+φ_l)^n]

is asymptotically (N→∞) given by

    V_N(t;λ)=√(2)/(2πλ(λ+1)N)^N/2[∏_n(1+1/λ_n)^t_n/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]
       ×exp[14λ^2+14λ-1/12λ(λ+1)]exp[∑_mε_m^2/16λ^2(λ+1)^2]exp[-(2λ+1)^2/128λ^3(λ+1)^3(∑_mε_m^2)^2]
       ×exp[-(2λ+1)^2N/128λ^3(λ+1)^3∑_mε_m^4]×(1+𝒪(N^-1/2+6ω+N^2+1/3C-Kαexp[N^4ω]))
       =√(2)(1+λ)^N2/(2πλ(λ+1)N)^N/2(1+1/λ)^x/2exp[14λ^2+14λ-1/12λ(λ+1)]
       exp[-1/2λ(λ+1)N∑_m(t_m-λ(N-1))^2]exp[-1/λ(λ+1)N^2∑_m(t_m-λ(N-1))^2]
       ×exp[2λ+1/6λ^2(λ+1)^2N^2∑_m(t_m-λ(N-1))^3]exp[-3λ^2+3λ+1/12λ^3(λ+1)^3N^3∑_m(t_m-λ(N-1))^4]
       ×exp[1/4λ^2(λ+1)^2N^4(∑_m(t_m-λ(N-1))^2)^2]×(1+𝒪(N^-1/2+6ω+N^2+1/3C-Kαexp[N^4ω]))  .

This is much larger than the Lower bound from Definition <ref>

    V_N(t;λ)/ℰ_α→∞  .



Define ε_j=2/N-2(t_j-λ(N-1)) and assume that |ε_j|≤λ N^-1/2+ω with 0<ω<1/14. It follows that ∑_jε_j=0.


To the integral V_N(t;λ) the expansion (<ref>) for k=4 in combination with (<ref>) and (<ref>) is applied. It will follow automatically that the higher orders (K>5) in this expansion will yield asymptotically irrelevant factors. This expansion produces the combinations (<ref>). 
Introducing δ-functions for S_1=∑_mφ_m, S_2=∑_mφ_m^2, T_3=∑_mε_mφ_m and T_4=∑_mε_m^2φ_m through their Fourier representation yields the integral

    V_N(t;λ)=2/(2π)^N[∏_n(1+1/λ_n)^t_n/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]∫τ_1∫ S_1∫τ_3
       ×∫ T_3∫ T_4∫τ_4∫ S_2∫τ_2 exp[2π i(τ_1S_1+τ_2S_2+τ_3T_3+τ_4T_4)
       -A_2S_1^2-2B_2S_1T_3+2C_2S_1T_4-2D_2T_3^2+3A_4S_2^2]
       ×{∏_j∫_-δ_N^δ_Nφ_j exp[iφ_j(-B_1Nε_j^2-B_1∑_mε_m^2-2πτ_1-3A_3S_2
       +C_1Nε_j^3-C_1ε_j∑_mε_m^2+C_1∑_mε_m^3-2πτ_3ε_j)]
       ×exp[-φ_j^2(A_2(N-2)+B_2(N-4)ε_j-(N-4)C_2ε_j^2-C_2∑_mε_m^2+2π iτ_2)]
       ×exp[-iφ_j^3(A_3(N-4)+4iA_4S_1)]
       ×exp[φ_j^4(A_4(N-8))] }  .

To ensure that that overall error consists of asymptotically irrelevant factors only, the φ_j-integral must be computed up to 𝒪(N^-1). Dividing the integration parameter φ_j by √(A_2(N-2)) shows that the φ_j-integral is of the form

    1/√(A_2(N-2))∫_-δ_N√(A_2(N-2))^δ_N√(A_2(N-2))φ exp[iφ Q_1/√(A_2(N-2))-φ^2Q_2-iφ^3Q_3/(A_2(N-2))^3/2+φ^4Q_4/(A_2(N-2))^2]
       =√(π/A_2(N-2))[Q_2+3iQ_3φ̃/(A_2(N-2))^3/2]^-1/2
       ×exp[iQ_1φ̃/√(A_2(N-2))-Q_2φ̃^2-iQ_3φ̃^3/(A_2(N-2))^3/2+Q_4φ̃^4/A_2^2(N-2)^2]
       ×{1-15Q_3^2/16(A_2(N-2))^3(Q_2+3iQ_3φ̃/(A_2(N-2))^3/2)^3+3Q_4/4A_2^2(N-2)^2(Q_2+3iQ_3φ̃/(A_2(N-2))^3/2)^2}  ,

which is calculated by the (<ref>) around the maximum φ̃ of the integrand. Observing that Q_1=𝒪(N^2ω), Q_2=𝒪(1) and Q_3,4=𝒪(N), shows that

    φ̃=iQ_1/2Q_2√(A_2(N-2))+𝒪(N^-3/2+4ω)=𝒪(N^-1/2+2ω)

is sufficient for the desired accuracy. This implies that

    exp[iQ_1φ̃/√(A_2(N-2))-Q_2φ̃^2-iQ_3φ̃^3/(A_2(N-2))^3/2]=exp[-Q_1^2/4A_2(N-2)]×(1+𝒪(N^-2+6ω))  .

The terms in square and curly brackets are then rewritten using

    1/√(1+y)≈ e^-y/2+y^2/4    and     1+z≈exp[z]

respectively. Using the same order of factors as in (<ref>), the result of the φ_j-integral is

    √(π/A_2(N-2))exp[-B_2(N-4)ε_j/2A_2(N-2)+C_2(N-4)ε_j^2/2A_2(N-2)+C_2∑_mε_m^2/2A_2(N-2)-iπτ_2/A_2(N-2)
       -3A_3B_1N(N-4)ε_j^2/4A_2^2(N-2)^2-3A_3B_1(N-4)∑_mε_m^2/4A_2^2(N-2)^2-3πτ_1A_3(N-4)/2A_2^2(N-2)^2-9A_3^2S_2(N-4)/4A_2^2(N-2)^2]
       ×exp[B_2^2(N-4)^2ε_j^2/4A_2^2(N-2)^2]exp[-B_1^2N^2ε_j^4/4A_2(N-2)-B_1^2Nε_j^2∑_mε_m^2/2A_2(N-2)
       -πτ_1B_1Nε_j^2/A_2(N-2)-3A_3B_1S_2Nε_j^2/2A_2(N-2)-B_1^2(∑_mε_m^2)^2/4A_2(N-2)-πτ_1B_1∑_mε_m^2/A_2(N-2)-3A_3B_1S_2∑_mε_m^2/2A_2(N-2)
       -π^2τ_1^2/A_2(N-2)-3πτ_1A_3S_2/A_2(N-2)-9A_3^2S_2^2/4A_2(N-2)]exp[-15A_3^2(N-4)^2/16A_2^3(N-2)^3]exp[3A_4(N-8)/4A_2^2(N-2)^2] .

Integrating now τ_2 yields a delta function that assigns the value

    S_2=N/2A_2(N-2)  .

Doing the same for τ_3 and τ_4 yields T_3=0 and T_4=0. The S_1-integral is

    ∫ S_1 exp[2π iτ_1S_1-A_2S_1^2]=√(π/A_2)exp[-π^2τ_1^2/A_2]  .

and the final integral

    ∫τ_1 exp[-2π^2(N-1)τ_1^2/A_2(N-2)-3πτ_1A_3N/A_2^2(N-2)-2πτ_1B_1N∑_mε_m^2/A_2(N-2)]
       =√(A_2(N-2)/2π (N-1))exp[9A_3^2N^2/8A_2^3(N-1)(N-2)+B_1^2N^2(∑_mε_m^2)^2/2A_2(N-1)(N-2)+3A_3B_1N^2∑_mε_m^2/2A_2^2(N-1)(N-2)]  .

Putting this all together yields

    V_N(t;λ)=√(2)/(2πλ(λ+1)N)^N/2[∏_n(1+1/λ_n)^t_n/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]
       ×exp[14λ^2+14λ-1/12λ(λ+1)]exp[∑_mε_m^2/16λ^2(λ+1)^2]exp[-(2λ+1)^2/128λ^3(λ+1)^3(∑_mε_m^2)^2]
       ×exp[-(2λ+1)^2N/128λ^3(λ+1)^3∑_mε_m^4]  .

Comparing (<ref>) to the Lower bound ℰ_α, it is immediately clear that V_n(t;λ) is much larger. Expand the products in square brackets around λ,

    [∏_j(1+1/λ_j)^λ(N-1)+(N-2)ε_j/2/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]
       =(1+1/λ)^x/2(1+λ)^N2exp[2λ(N-1)+(N-2)ε_j/4log(1+λ+ε_j/1+λλ/λ+ε_j)]
       ×exp[-∑_k<llog(1+λ-√((1+λ-1+λ/1+λ+ε_k)(1+λ-1+λ/1+λ+ε_l)))]
       =(1+1/λ)^x/2(1+λ)^N2
       ×exp[(∑_mε_m^2)·[Nλ^2/4λ^2(λ+1)^2+λ/4λ^2(λ+1)^2-(N-1)3λ^2+λ/8λ^2(λ+1)^2-1/8λ(λ+1)]
       +(∑_mε_m^3)·[-N6λ^3+3λ^2+λ/24λ^3(λ+1)^3+N14λ^3+9λ^2+3λ/48λ^3(λ+1)^3]
       +(∑_mε_m^4)·[N6λ^4+6λ^3+4λ^2+λ/24λ^4(λ+1)^4-N30λ^4+28λ^3+19λ^2+5λ/128λ^4(λ+1)^4]
       +(∑_mε_m^2)^2·[6λ^2+6λ+1/128λ^3(λ+1)^3]×(1+𝒪(N^-1/2+5ω))  ,

yields combined with (<ref>) the desired result.

To determine the error from the difference 𝒟 from Lemma <ref>, we divide it by V_N(t;λ). Assuming that |t_j-λ(N-1)|=λ N^1/2+ω takes maximal values, it follows that the relative difference is at most

    𝒪(N^2-Kα)exp[4λ^2+4λ+1/3λ(λ+1)]exp[4λ^2+4λ+3/4(λ+1)^2N^2ω]
       ×exp[(2λ+1)^2λ N^4ω/8(λ+1)^3]exp[(2λ+1)^2λ N^4ω/8(λ+1)^3]  .

Only the first exponential can become large, if λ is small. Assuming that λ > C/log(N), this factor adds an error N^1/3C.

To keep this relative error small, it is furthermore necessary that exp[N^4ω]≪ N^Kα-2. Solving this yields

    0<ω<log(Kα-2)+log(log(N))/4log(N)  .



Choosing the value of λ may seem arbitrary at first. It is not. Comparing (<ref>) to the Lower bound ℰ_1/2-r for some small r>0, the outcome is only much larger, if

    ∑_mε_m^4≪ N^2r    and    ∑_mε_m^2≪ N^r  .

It follows that λ N(N-1)=x in the limit. In <cit.> the number of matrices V_N(t;λ) has been calculated for the case that all t_j are equal. They require λ to be the average matrix entry for infinitely large matrices. Because Lemma <ref> covers this case too, the same value for λ had to be expected.




For any ω∈(0,loglog N/6log N) and α∈(0,1/4-ω), define

    δ_N=N^-αζ_N/min{λ_j}  ,

such that ζ_N→∞ and N^-δζ_N→0 for any δ>0. Assuming that x=∑_jt_j=λ N(N-1), |t_j-λ(N-1)|≪λ N^1/2+ω and λ>C/log(N), the integral

    V_N(t;λ)=2/(2π)^N(∏_j=1^N(1+1/λ_j)^t_j/2)·(∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       ×∫_[-δ_N,δ_N]^Nφ exp[-i∑_jφ_jt_j]exp[∑_n=1^K-1i^n∑_k<lA_n(√(λkλ_l)/√((1+λk)(1+λ_l))-√(λkλ_l))·(φ_k+φ_l)^n]

is then for any ⌊2/α⌋+1≤ K≤log(N) given by

    V_N(t;λ)=√(2)(1+λ)^N2/(2πλ(λ+1)N)^N/2(1+1/λ)^x/2exp[14λ^2+14λ-1/12λ(λ+1)]
       exp[-1/2λ(λ+1)N∑_m(t_m-λ(N-1))^2]
       ×exp[-64λ^6-192λ^5-160λ^4+40λ^2+8λ+1/64λ^4(λ+1)^4N^2∑_m(t_m-λ(N-1))^2]
       ×exp[2λ+1/6λ^2(λ+1)^2N^2∑_m(t_m-λ(N-1))^3]exp[-3λ^2+3λ+1/12λ^3(λ+1)^3N^3∑_m(t_m-λ(N-1))^4]
       ×exp[6λ^2+6λ+1/8λ^3(λ+1)^3N^4(∑_m(t_m-λ(N-1))^2)^2]
       ×exp[-8λ^4+16λ^3-8λ-1/48λ^5(λ+1)^5N^6(∑_m(t_m-λ(N-1))^2)^3]
       ×(1+𝒪(N^-1/2+6ω+N^2+1/3C-Kα))  .

where 𝒟 is the difference from Lemma <ref>.




Define ε_j=2/N-2(t_j-λ(N-1)) and assume that |ε_j|≤λ N^-1/2+ω with 0<ω<1/14. It follows that

    |ε_j/λ|≪ N^-1/2+ω    and    ∑_jε_j=0  .

To the integral V_N(t;λ) the expansion (<ref>) for k=4 in combination with (<ref>) and (<ref>) is applied. It will follow automatically that the higher orders in this expansion will yield smaller contributions. This expansion produces the combinations (<ref>) and (<ref>). Introducing δ-functions for S_1=∑_mφ_m, S_2=∑_mφ_m^2, T_3=∑_mε_mφ_m, T_4=∑_mε_m^2φ_m and T_5=∑_mε_mφ_m^2 through their Fourier representation yields the integral

    V_N(t;λ)=2/(2π)^N[∏_n(1+1/λ_n)^t_n/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]∫τ_1∫ S_1∫τ_3
       ×∫ T_3∫τ_4∫ T_4∫ S_2∫τ_2∫ T_5∫τ_5 exp[2π i(τ_1S_1+τ_2S_2+τ_3T_3+τ_4T_4+τ_5T_5)
       -A_2S_1^2-2B_2S_1T_3-2C_2T_3^2+3A_4S_2^2-D_2S_2∑_mε_m^3-E_2S_2∑_mε_m^4-F_2T_5∑_mε_m^3-G_2T_4^2
       +6B_4T_5S_2+3D_4T_5^2]
       ×{∏_j∫_-δ_N^δ_Nφ_j exp[iφ_j(-B_1Nε_j^2-B_1∑_mε_m^2-2πτ_1-3B_3T_5-3A_3S_2
       +C_1Nε_j^3-C_1ε_j∑_mε_m^2+C_1∑_mε_m^3-2πτ_3ε_j-3D_3ε_jT_5-2iD_2T_4ε_j-3B_3S_2ε_j
       -D_1(N-2)ε_j^4-D_1∑_mε_m^4+E_1ε_j∑_mε_m^3-2πτ_4ε_j^2+3C_3S_2ε_j^2+F_1ε_j^2∑_mε_m^2
       -2iC_2S_1ε_j^2-2E_1ε_j^4-F_1ε_j^4+2iD_2S_1ε_j^3+2iE_2S_1ε_j^4+2iF_2T_3ε_j^3)]
       ×exp[-φ_j^2(A_2(N-2)+B_2(N-4)ε_j-NC_2ε_j^2-C_2∑_mε_m^2+2π iτ_2+ND_2ε_j^3
       +2π iτ_5ε_j-D_2ε_j∑_mε_m^2+E_2(N-4)ε_j^4+G_2ε_j^2∑_mε_m^2-3iC_3S_1ε_j^2-6C_4S_2ε_j^2
       ×-4F_2ε_j^4-2G_2ε_j^4)]
       ×exp[-iφ_j^3(A_3(N-4)+B_3(N-8)ε_j-C_3(N-8)ε_j^2-C_3∑_mε_m^2
       +4iA_4S_1+4iB_4T_3+4iC_4T_4+4iD_4T_3ε_j+4iB_4S_1ε_j-4D_3ε_j^2+4iC_4S_1ε_j^2)]
       ×exp[φ_j^4(A_4(N-8)+B_4(N-16)ε_j+C_4(N-16)ε_j^2+C_4∑_mε_m^2-8D_4ε_j^2)] }  .

To ensure that that overall error is of order 𝒪(N^-1), the φ_j-integral must be computed up to 𝒪(N^-2). Dividing the integration parameter φ_j by √(A_2(N-2)) shows that the φ_j-integral is of the form

    1/√(A_2(N-2))∫_-δ_N√(A_2(N-2))^δ_N√(A_2(N-2))φ exp[iφ Q_1/√(A_2(N-2))-φ^2Q_2-iφ^3Q_3/(A_2(N-2))^3/2+φ^4Q_4/(A_2(N-2))^2]
       =√(π/A_2(N-2))[Q_2+3iQ_3φ̃/(A_2(N-2))^3/2]^-1/2exp[iQ_1φ̃/√(A_2(N-2))-Q_2φ̃^2-iQ_3φ̃^3/(A_2(N-2))^3/2]
       ×{1-15Q_3^2/16(A_2(N-2))^3(Q_2+3iQ_3φ̃/(A_2(N-2))^3/2)^3+3Q_4/4A_2^2(N-2)^2(Q_2+3iQ_3φ̃/(A_2(N-2))^3/2)^2}  ,

which is calculated by the (<ref>) around the maximum φ̃ of the integrand. Observing that Q_1=𝒪(N^2ω), Q_2=𝒪(1) and Q_3,4=𝒪(N), shows that

    φ̃=iQ_1/2Q_2√(A_2(N-2))+3iQ_1^2Q_3/8Q_2^3(A_2(N-2))^5/2+𝒪(N^-5/2+6ω)=𝒪(N^-1/2+2ω)

is sufficient for the desired accuracy. This implies that

    exp[iQ_1φ̃/√(A_2(N-2))-Q_2φ̃^2-iQ_3φ̃^3/(A_2(N-2))^3/2]=exp[-Q_1^2/4A_2(N-2)-Q_3Q_1^3/8Q_2^3(A_2(N-2))^3]×(1+𝒪(N^-2))  .

The terms in square and curly brackets are then rewritten using

    1/√(1+y)=e^-y/2+y^2/4-y^3/6+y^4/8×(1+𝒪(N^-2))  and   1+z=exp[z]×(1+𝒪(N^-2))

respectively. 

It is a tedious but straightforward exercise to show that[It serves the calculations and checks to write to this formula down.]

    φ̃=-iB_1Nε_j^2/2√(A_2(N-2))-iB_1∑_mε_m^2/2√(A_2(N-2))-π iτ_1/√(A_2(N-2))-3iB_3T_5/2√(A_2(N-2))-3iA_3S_2/2√(A_2(N-2))
       +iC_1Nε_j^3/2√(A_2(N-2))-iC_1ε_j∑_mε_m^2/2√(A_2(N-2))+iC_1∑_mε_m^3/2√(A_2(N-2))-iπτ_3ε_j/√(A_2(N-2))
       -3iD_3T_5ε_j/2√(A_2(N-2))+D_2T_4ε_j/√(A_2(N-2))-3iB_3S_2ε_j/2√(A_2(N-2))-iD_1(N-2)ε_j^4/2√(A_2(N-2))-iD_1∑_mε_m^4/2√(A_2(N-2))
       +iE_1ε_j∑_mε_m^3/2√(A_2(N-2))-π iτ_4ε_j^2/√(A_2(N-2))+3iC_3S_2ε_j^2/2√(A_2(N-2))+iF_1ε_j^2∑_mε_m^2/2√(A_2(N-2))+C_2S_1ε_j^2/√(A_2(N-2))
       +iB_1B_2N(N-4)ε_j^3/2(A_2(N-2))^3/2+iB_1B_2(N-4)ε_j∑_mε_m^2/2(A_2(N-2))^3/2+iπτ_1B_2(N-4)ε_j/2(A_2(N-2))^3/2
       +3iB_2B_3T_5(N-4)ε_j/2(A_2(N-2))^3/2+3iA_3B_2S_2(N-4)ε_j/2(A_2(N-2))^3/2-iB_2C_1N(N-4)ε_j^4/2(A_2(N-2))^3/2
       +iB_2C_1(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^3/2-iB_2C_1(N-4)ε_j∑_mε_m^3/2(A_2(N-2))^3/2+iπτ_3B_2(N-4)ε_j^2/(A_2(N-2))^3/2
       +3iB_2D_3T_5(N-4)ε_j^2/2(A_2(N-2))^3/2-B_2D_2T_4(N-4)ε_j^2/(A_2(N-2))^3/2+3iB_2B_3S_2(N-4)ε_j^2/2(A_2(N-2))^3/2
       -iB_1C_2N^2ε_j^4/2(A_2(N-2))^3/2-iB_1C_2Nε_j^2∑_mε_m^2/2(A_2(N-2))^3/2-iπτ_1NC_2ε_j^2/(A_2(N-2))^3/2-3iB_3C_2T_5Nε_j^2/2(A_2(N-2))^3/2
       -3iA_3C_2S_2Nε_j^2/2(A_2(N-2))^3/2-iB_1C_2Nε_j^2∑_mε_m^2/2(A_2(N-2))^3/2-iB_1C_2(∑_mε_m^2)^2/2(A_2(N-2))^3/2-iπτ_1C_2∑_mε_m^2/(A_2(N-2))^3/2
       -3iB_3C_2T_5∑_mε_m^2/2(A_2(N-2))^3/2-3iA_3C_2S_2∑_mε_m^2/2(A_2(N-2))^3/2-πτ_2B_1Nε_j^2/(A_2(N-2))^3/2-πτ_2B_1∑_mε_m^2/(A_2(N-2))^3/2
       -2π^2τ_1τ_2/(A_2(N-2))^3/2-3πτ_2B_3T_5/(A_2(N-2))^3/2-3πτ_2A_3S_2/(A_2(N-2))^3/2
       -iB_1B_2^2N(N-4)^2ε_j^4/2(A_2(N-2))^5/2-iB_1B_2^2(N-4)^2ε_j^2∑_mε_m^2/2(A_2(N-2))^5/2-iπτ_1B_2^2(N-4)^2ε_j^2/(A_2(N-2))^5/2


    -3iB_2^2B_3T_5(N-4)^2ε_j^2/2(A_2(N-2))^5/2-3iA_3B_2^2S_2(N-4)^2ε_j^2/2(A_2(N-2))^5/2+3iA_3B_1^2N^2(N-4)ε_j^4/8(A_2(N-2))^5/2
       +3iA_3B_1^2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^5/2+3iπτ_1A_3B_1N(N-4)ε_j^2/2(A_2(N-2))^5/2+9iA_3B_1B_3T_5N(N-4)ε_j^2/4(A_2(N-2))^5/2
       +9iA_3^2B_1S_2N(N-4)ε_j^2/4(A_2(N-2))^5/2+3iA_3B_1^2(N-4)(∑_mε_m^2)^2/8(A_2(N-2))^5/2+3iπτ_1A_3B_1(N-4)∑_mε_m^2/2(A_2(N-2))^5/2
       +9iA_3B_1B_3T_5(N-4)∑_mε_m^2/4(A_2(N-2))^5/2+9iA_3^2B_1S_2(N-4)∑_mε_m^2/4(A_2(N-2))^5/2+3iπ^2τ_1^2A_3(N-4)/2(A_2(N-2))^5/2
       +9iπτ_1A_3B_3T_5(N-4)/2(A_2(N-2))^5/2+9iπτ_1A_3^2S_2(N-4)/2(A_2(N-2))^5/2+27iA_3B_3^2T_5^2(N-4)/8(A_2(N-2))^5/2
       +27iA_3^2B_3S_2T_5(N-4)/4(A_2(N-2))^5/2+27iA_3^3S_2^2(N-4)/8(A_2(N-2))^5/2+𝒪(N^-3/2-2ω)  .

It is furthermore helpful to expand[It serves the calculations and checks to write to this formula down.]

    Q_2+3Q_3iφ̃/(A_2(N-2))^3/2=1+B_2(N-4)ε_j/A_2(N-2)-C_2Nε_j^2/A_2(N-2)-C_2∑_mε_m^2/A_2(N-2)
       +2iπτ_2/A_2(N-2)+D_2Nε_j^3/A_2(N-2)+2iπτ_5ε_j/A_2(N-2)-D_2ε_j∑_mε_m^2/A_2(N-2)+E_2(N-4)ε_j^4/A_2(N-2)
       +G_2ε_j^2∑_mε_m^2/A_2(N-2)-3iC_3S_1ε_j^2/A_2(N-2)-6C_4S_2ε_j^2/A_2(N-2)
    +3A_3B_1N(N-4)ε_j^2/2(A_2(N-2))^2
       +3A_3B_1(N-4)∑_mε_m^2/2(A_2(N-2))^2+3πτ_1A_3(N-4)/(A_2(N-2))^2+9A_3B_3T_5(N-4)/2(A_2(N-2))^2
       +9A_3^2S_2(N-4)/2(A_2(N-2))^2-3A_3C_1N(N-4)ε_j^3/2(A_2(N-2))^2+3A_3C_1(N-4)ε_j∑_mε_m^2/2(A_2(N-2))^2
       -3A_3C_1(N-4)∑_mε_m^3/2(A_2(N-2))^2+3πτ_3A_3(N-4)ε_j/(A_2(N-2))^2+9A_3D_3T_5(N-4)ε_j/2(A_2(N-2))^2
       +3iA_3D_2T_4(N-4)ε_j/(A_2(N-2))^2+9A_3B_3S_2(N-4)ε_j/2(A_2(N-2))^2+3A_3D_1(N-2)(N-4)ε_j^4/2(A_2(N-2))^2
       +3A_3D_1(N-4)∑_mε_m^4/2(A_2(N-2))^2-3A_3E_1(N-4)ε_j∑_mε_m^3/2(A_2(N-2))^2+3πτ_4A_3(N-4)ε_j^2/(A_2(N-2))^2
       -9A_3C_3S_2(N-4)ε_j^2/2(A_2(N-2))^2-3A_3F_1(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^2+3iA_3C_2S_1(N-4)ε_j^2/(A_2(N-2))^2


    -3A_3B_1B_2N(N-4)^2ε_j^3/2(A_2(N-2))^3-3A_3B_1B_2(N-4)^2ε_j∑_mε_m^2/2(A_2(N-2))^3-3πτ_1A_3B_2(N-4)^2ε_j/2(A_2(N-2))^3
       -9A_3B_2B_3T_5(N-4)^2ε_j/2(A_2(N-2))^3-9A_3^2B_2S_2(N-4)^2ε_j/2(A_2(N-2))^3+3A_3B_2C_1N(N-4)^2ε_j^4/2(A_2(N-2))^3
       -3A_3B_2C_1(N-4)^2ε_j^2∑_mε_m^2/2(A_2(N-2))^3+3A_3B_2C_1(N-4)^2ε_j∑_mε_m^3/2(A_2(N-2))^3-3πτ_3A_3B_2(N-4)^2ε_j^2/2(A_2(N-2))^3
       -9A_3B_2D_3T_5(N-4)^2ε_j^2/2(A_2(N-2))^3-3iA_3B_2D_2T_4(N-4)^2ε_j^2/(A_2(N-2))^3-9A_3B_2B_3S_2(N-4)^2ε_j^2/2(A_2(N-2))^3
       +3A_3B_1C_2N^2(N-4)ε_j^4/2(A_2(N-2))^3+3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^3+3πτ_1A_3C_2N(N-4)ε_j^2/(A_2(N-2))^3
       +9A_3B_3C_2T_5N(N-4)ε_j^2/2(A_2(N-2))^3+9A_3^2C_2S_2N(N-4)ε_j^2/2(A_2(N-2))^3+3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^3
       +3A_3B_1C_2(N-4)(∑_mε_m^2)^2/2(A_2(N-2))^3+3πτ_1A_3C_2(N-4)∑_mε_m^2/(A_2(N-2))^3+9A_3B_3C_2T_5(N-4)∑_mε_m^2/2(A_2(N-2))^3
       +9A_3^2C_2S_2(N-4)∑_mε_m^2/2(A_2(N-2))^3-3iπτ_2A_3B_1N(N-4)ε_j^2/(A_2(N-2))^3-3iπτ_2A_3B_1(N-4)∑_mε_m^2/(A_2(N-2))^3
       +3A_3B_1B_2^2N(N-4)^3ε_j^4/2(A_2(N-2))^4+3A_3B_1B_2^2(N-4)^3ε_j^2∑_mε_m^2/2(A_2(N-2))^4+3πτ_1A_3B_2^2(N-4)^3ε_j^2/(A_2(N-2))^4
       +9A_3B_2^2B_3T_5(N-4)^3ε_j^2/2(A_2(N-2))^4+9A_3^2B_2^2S_2(N-4)^3ε_j^2/2(A_2(N-2))^4-9A_3^2B_1^2N^2(N-4)^2ε_j^4/8(A_2(N-2))^4
       -9A_3^2B_1^2N(N-4)^2ε_j^2∑_mε_m^2/4(A_2(N-2))^4-9πτ_1A_3^2B_1N(N-4)^2ε_j^2/2(A_2(N-2))^4-27A_3^2B_1B_3T_5N(N-4)^2ε_j^2/4(A_2(N-2))^4
       -27A_3^3B_1S_2N(N-4)^2ε_j^2/4(A_2(N-2))^4-9A_3^2B_1^2(N-4)^2(∑_mε_m^2)^2/8(A_2(N-2))^4-9πτ_1A_3^2B_1(N-4)^2∑_mε_m^2/2(A_2(N-2))^4
       -27A_3^2B_1B_3T_5(N-4)^2∑_mε_m^2/4(A_2(N-2))^4-27A_3^3B_1S_2(N-4)^2∑_mε_m^2/4(A_2(N-2))^4+3B_1B_3N(N-8)ε_j^3/2(A_2(N-2))^2
       +3B_1B_3(N-8)ε_j∑_mε_m^2/2(A_2(N-2))^2+3πτ_1B_3(N-8)ε_j/(A_2(N-2))^2+9B_3^2T_5(N-8)ε_j/2(A_2(N-2))^2
       +9A_3B_3S_2(N-8)ε_j/2(A_2(N-2))^2-3B_3C_1N(N-8)ε_j^4/2(A_2(N-2))^2+3B_3C_1(N-8)ε_j^2∑_mε_m^2/2(A_2(N-2))^2
       -3B_3C_1(N-8)ε_j∑_mε_m^3/2(A_2(N-2))^2+3πτ_3B_3(N-8)ε_j^2/(A_2(N-2))^2+9B_3D_3T_5(N-8)ε_j^2/2(A_2(N-2))^2
       +3iB_3D_2T_4(N-8)ε_j^2/(A_2(N-2))^2+9B_3^2S_2(N-8)ε_j^2/2(A_2(N-2))^2-3B_1B_2B_3N(N-4)(N-8)ε_j^4/2(A_2(N-2))^3
       -3B_1B_2B_3(N-4)(N-8)ε_j^2∑_mε_m^2/2(A_2(N-2))^3-3πτ_1B_2B_3(N-4)(N-8)ε_j^2/2(A_2(N-2))^3
       -9B_2B_3^2T_5(N-4)(N-8)ε_j^2/2(A_2(N-2))^3-9A_3B_2B_3S_2(N-4)(N-8)ε_j^2/2(A_2(N-2))^3-3B_1C_3N(N-8)ε_j^4/2(A_2(N-2))^2
       -3B_1C_3(N-8)ε_j^2∑_mε_m^2/2(A_2(N-2))^2-3πτ_1C_3(N-8)ε_j^2/(A_2(N-2))^2-9B_3C_3T_5(N-8)ε_j^2/2(A_2(N-2))^2


    -9A_3C_3S_2(N-8)ε_j^2/2(A_2(N-2))^2-3B_1C_3Nε_j^2∑_mε_m^2/2(A_2(N-2))^2-3B_1C_3(∑_mε_m^2)^2/2(A_2(N-2))^2
       -3πτ_1C_3∑_mε_m^2/(A_2(N-2))^2-9B_3C_3T_5∑_mε_m^2/2(A_2(N-2))^2-9A_3C_3S_2∑_mε_m^2/2(A_2(N-2))^2
       +6iA_4B_1S_1Nε_j^2/(A_2(N-2))^2+6iA_4B_1S_1∑_mε_m^2/(A_2(N-2))^2+6iB_1B_4T_3Nε_j^2/(A_2(N-2))^2+6iB_1B_4T_3∑_mε_m^2/(A_2(N-2))^2
       +6iB_1C_4T_4Nε_j^2/(A_2(N-2))^2+6iB_1C_4T_4∑_mε_m^2/(A_2(N-2))^2+𝒪(N^-2)  .


Writing out (<ref>), where the order of factors is preserved, yields

    √(π/A_2(N-2))exp[-B_2(N-4)ε_j/2A_2(N-2)+C_2Nε_j^2/2A_2(N-2)+C_2∑_mε_m^2/2A_2(N-2)
       -iπτ_2/A_2(N-2)-D_2Nε_j^3/2A_2(N-2)-iπτ_5ε_j/A_2(N-2)+D_2ε_j∑_mε_m^2/2A_2(N-2)-E_2(N-4)ε_j^4/2A_2(N-2)
       -G_2ε_j^2∑_mε_m^2/2A_2(N-2)+3iC_3S_1ε_j^2/2A_2(N-2)+3C_4S_2ε_j^2/A_2(N-2)
    -3A_3B_1N(N-4)ε_j^2/4(A_2(N-2))^2
       -3A_3B_1(N-4)∑_mε_m^2/4(A_2(N-2))^2-3πτ_1A_3(N-4)/2(A_2(N-2))^2-9A_3B_3T_5(N-4)/4(A_2(N-2))^2
       -9A_3^2S_2(N-4)/4(A_2(N-2))^2+3A_3C_1N(N-4)ε_j^3/4(A_2(N-2))^2-3A_3C_1(N-4)ε_j∑_mε_m^2/4(A_2(N-2))^2
       +3A_3C_1(N-4)∑_mε_m^3/4(A_2(N-2))^2-3πτ_3A_3(N-4)ε_j/2(A_2(N-2))^2-9A_3D_3T_5(N-4)ε_j/4(A_2(N-2))^2
       -3iA_3D_2T_4(N-4)ε_j/2(A_2(N-2))^2-9A_3B_3S_2(N-4)ε_j/4(A_2(N-2))^2-3A_3D_1(N-2)(N-4)ε_j^4/4(A_2(N-2))^2
       -3A_3D_1(N-4)∑_mε_m^4/4(A_2(N-2))^2+3A_3E_1(N-4)ε_j∑_mε_m^2/4(A_2(N-2))^2-3πτ_4A_3(N-4)ε_j^2/2(A_2(N-2))^2
       +9A_3C_3S_2(N-4)ε_j^2/4(A_2(N-2))^2+3A_3F_1(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^2-3iA_3C_2S_1(N-4)ε_j^2/2(A_2(N-2))^2
       +3A_3B_1B_2N(N-4)^2ε_j^3/4(A_2(N-2))^3+3A_3B_1B_2(N-4)^2ε_j∑_mε_m^2/4(A_2(N-2))^3+3πτ_1A_3B_2(N-4)^2ε_j/4(A_2(N-2))^3
       +9A_3B_2B_3T_5(N-4)^2ε_j/4(A_2(N-2))^3+9A_3^2B_2S_2(N-4)^2ε_j/4(A_2(N-2))^3-3iA_3B_2^2T_3(N-4)^2ε_j/2(A_2(N-2))^3
       -3A_3B_2C_1N(N-4)^2ε_j^4/4(A_2(N-2))^3+3A_3B_2C_1(N-4)^2ε_j^2∑_mε_m^2/4(A_2(N-2))^3
       -3A_3B_2C_1(N-4)^2ε_j∑_mε_m^3/4(A_2(N-2))^3+3πτ_3A_3B_2(N-4)^2ε_j^2/4(A_2(N-2))^3+9A_3B_2D_3T_5(N-4)^2ε_j^2/4(A_2(N-2))^3
       +3iA_3B_2D_2T_4(N-4)^2ε_j^2/2(A_2(N-2))^3+9A_3B_2B_3S_2(N-4)^2ε_j^2/4(A_2(N-2))^3-3A_3B_1C_2N^2(N-4)ε_j^4/4(A_2(N-2))^3
       -3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3-3πτ_1A_3C_2N(N-4)ε_j^2/2(A_2(N-2))^3-9A_3B_3C_2T_5N(N-4)ε_j^2/4(A_2(N-2))^3
       -9A_3^2C_2S_2N(N-4)ε_j^2/4(A_2(N-2))^3-3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3-3A_3B_1C_2(N-4)(∑_mε_m^2)^2/4(A_2(N-2))^3
       -3πτ_1A_3C_2(N-4)∑_mε_m^2/2(A_2(N-2))^3-9A_3B_3C_2T_5(N-4)∑_mε_m^2/4(A_2(N-2))^3-9A_3^2C_2S_2(N-4)∑_mε_m^2/4(A_2(N-2))^3
       +3iπτ_2A_3B_1N(N-4)ε_j^2/2(A_2(N-2))^3+3iπτ_2A_3B_1(N-4)∑_mε_m^2/2(A_2(N-2))^3-3A_3B_1B_2^2N(N-4)^3ε_j^4/4(A_2(N-2))^4
       -3A_3B_1B_2^2(N-4)^3ε_j^2∑_mε_m^2/4(A_2(N-2))^4-3πτ_1A_3B_2^2(N-4)^3ε_j^2/2(A_2(N-2))^4-9A_3B_2^2B_3T_5(N-4)^3ε_j^2/4(A_2(N-2))^4


    -9A_3^2B_2^2S_2(N-4)^3ε_j^2/4(A_2(N-2))^4+9A_3^2B_1^2N^2(N-4)^2ε_j^4/16(A_2(N-2))^4+9A_3^2B_1^2N(N-4)^2ε_j^2∑_mε_m^2/8(A_2(N-2))^4
       +9πτ_1A_3^2B_1N(N-4)^2ε_j^2/4(A_2(N-2))^4+27A_3^2B_1B_3T_5N(N-4)^2ε_j^2/8(A_2(N-2))^4+27A_3^3B_1S_2N(N-4)^2ε_j^2/8(A_2(N-2))^4
       +9A_3^2B_1^2(N-4)^2(∑_mε_m^2)^2/16(A_2(N-2))^4+9πτ_1A_3^2B_1(N-4)^2∑_mε_m^2/4(A_2(N-2))^4+27A_3^2B_1B_3T_5(N-4)^2∑_mε_m^2/8(A_2(N-2))^4
       +27A_3^3B_1S_2(N-4)^2∑_mε_m^2/8(A_2(N-2))^4-3B_1B_3N(N-8)ε_j^3/4(A_2(N-2))^2-3B_1B_3(N-8)ε_j∑_mε_m^2/4(A_2(N-2))^2
       -3πτ_1B_3(N-8)ε_j/2(A_2(N-2))^2-9B_3^2T_5(N-8)ε_j/4(A_2(N-2))^2-9A_3B_3S_2(N-8)ε_j/4(A_2(N-2))^2+3B_3C_1N(N-8)ε_j^4/4(A_2(N-2))^2
       -3B_3C_1(N-8)ε_j^2∑_mε_m^2/4(A_2(N-2))^2+3B_3C_1(N-8)ε_j∑_mε_m^3/4(A_2(N-2))^2-3πτ_3B_3(N-8)ε_j^2/2(A_2(N-2))^2
       -9B_3D_3T_5(N-8)ε_j^2/4(A_2(N-2))^2-3iB_3D_2T_4(N-8)ε_j^2/2(A_2(N-2))^2-9B_3^2S_2(N-8)ε_j^2/4(A_2(N-2))^2
       +3B_1B_2B_3N(N-4)(N-8)ε_j^4/4(A_2(N-2))^3+3B_1B_2B_3(N-4)(N-8)ε_j^2∑_mε_m^2/4(A_2(N-2))^3
       +3πτ_1B_2B_3(N-4)(N-8)ε_j^2/4(A_2(N-2))^3+9B_2B_3^2T_5(N-4)(N-8)ε_j^2/4(A_2(N-2))^3
       +9A_3B_2B_3S_2(N-4)(N-8)ε_j^2/4(A_2(N-2))^3+3B_1C_3N(N-8)ε_j^4/4(A_2(N-2))^2+3B_1C_3(N-8)ε_j^2∑_mε_m^2/4(A_2(N-2))^2
       +3πτ_1C_3(N-8)ε_j^2/2(A_2(N-2))^2+9B_3C_3T_5(N-8)ε_j^2/4(A_2(N-2))^2+9A_3C_3S_2(N-8)ε_j^2/4(A_2(N-2))^2
       +3B_1C_3Nε_j^2∑_mε_m^2/4(A_2(N-2))^2+3B_1C_3(∑_mε_m^2)^2/4(A_2(N-2))^2+3πτ_1C_3∑_mε_m^2/2(A_2(N-2))^2+9B_3C_3T_5∑_mε_m^2/4(A_2(N-2))^2
       +9A_3C_3S_2∑_mε_m^2/4(A_2(N-2))^2-3iA_4B_1S_1Nε_j^2/(A_2(N-2))^2-3iA_4B_1S_1∑_mε_m^2/(A_2(N-2))^2-3iB_1B_4T_3Nε_j^2/(A_2(N-2))^2
       -3iB_1B_4T_3∑_mε_m^2/(A_2(N-2))^2-3iB_1C_4T_4Nε_j^2/(A_2(N-2))^2-3iB_1C_4T_4∑_mε_m^2/(A_2(N-2))^2]
       ×exp[B_2^2(N-4)^2ε_j^2/4(A_2(N-2))^2-B_2C_2N(N-4)ε_j^3/2(A_2(N-2))^2-B_2C_2(N-4)ε_j∑_mε_m^2/2(A_2(N-2))^2+iπτ_2B_2(N-4)ε_j/(A_2(N-2))^2
       +B_2D_2N(N-4)ε_j^4/2(A_2(N-2))^2+iπτ_5B_2(N-4)ε_j^2/(A_2(N-2))^2-B_2D_2(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^2
       +3A_3B_1B_2N(N-4)^2ε_j^3/4(A_2(N-2))^3+3A_3B_1B_2(N-4)^2ε_j∑_mε_m^2/4(A_2(N-2))^3+3πτ_1A_3B_2(N-4)^2ε_j/2(A_2(N-2))^3
       +9A_3B_2B_3T_5(N-4)^2ε_j/4(A_2(N-2))^3+9A_3^2B_2S_2(N-4)^2ε_j/4(A_2(N-2))^3-3A_3B_2C_1N(N-4)^2ε_j^4/4(A_2(N-2))^3
       +3A_3B_2C_1(N-4)^2ε_j^2∑_mε_m^2/4(A_2(N-2))^3-3A_3B_2C_1(N-4)^2ε_j∑_mε_m^3/4(A_2(N-2))^3+3πτ_3A_3B_2(N-4)^2ε_j^2/2(A_2(N-2))^3
       +9A_3B_2D_3T_5(N-4)^2ε_j^2/4(A_2(N-2))^3+3iA_3B_2D_2T_4(N-4)^2ε_j^2/2(A_2(N-2))^3+9A_3B_2B_3S_2(N-4)^2ε_j^2/4(A_2(N-2))^3


    -3A_3B_1B_2^2N(N-4)^3ε_j^4/4(A_2(N-2))^4-3A_3B_1B_2^2(N-4)^3ε_j^2∑_mε_m^2/4(A_2(N-2))^4-3πτ_1A_3B_2^2(N-4)^3ε_j^2/4(A_2(N-2))^4
       -9A_3B_2^2B_3T_5(N-4)^3ε_j^2/4(A_2(N-2))^4-9A_3^2B_2^2S_2(N-4)^3ε_j^2/4(A_2(N-2))^4+3B_1B_2B_3N(N-4)(N-8)ε_j^4/4(A_2(N-2))^3
       +3B_1B_2B_3(N-4)(N-8)ε_j^2∑_mε_m^2/4(A_2(N-2))^3+3πτ_1B_2B_3(N-4)(N-8)ε_j^2/2(A_2(N-2))^3
       +9B_2B_3^2T_5(N-4)(N-8)ε_j^2/4(A_2(N-2))^3+9A_3B_2B_3S_2(N-4)(N-8)ε_j^2/4(A_2(N-2))^3+C_2^2N^2ε_j^4/4(A_2(N-2))^2
       +C_2^2Nε_j^2∑_mε_m^2/2(A_2(N-2))^2-iπτ_2C_2Nε_j^2/(A_2(N-2))^2-3A_3B_1C_2N^2(N-4)ε_j^4/4(A_2(N-2))^3
       -3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3-3πτ_1A_3C_2N(N-4)ε_j^2/2(A_2(N-2))^3-9A_3B_3C_2T_5N(N-4)ε_j^2/4(A_2(N-2))^3
       -9A_3^2C_2S_2N(N-4)ε_j^2/4(A_2(N-2))^3+C_2^2(∑_mε_m^2)^2/4(A_2(N-2))^2-iπτ_2C_2∑_mε_m^2/(A_2(N-2))^2
       -3A_3B_1C_2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3-3A_3B_1C_2(N-4)(∑_mε_m^2)^2/4(A_2(N-2))^3
       -3πτ_1A_3C_2(N-4)∑_mε_m^2/2(A_2(N-2))^3-9A_3B_3C_2T_5(N-4)∑_mε_m^2/4(A_2(N-2))^3-9A_3^2C_2S_2(N-4)∑_mε_m^2/4(A_2(N-2))^3
       +3iπτ_2A_3B_1N(N-4)ε_j^2/2(A_2(N-2))^3+3iπτ_2A_3B_1(N-4)∑_mε_m^2/2(A_2(N-2))^3+9A_3^2B_1^2N^2(N-4)^2ε_j^4/16(A_2(N-2))^4
       +9A_3^2B_1^2N(N-4)^2ε_j^2∑_mε_m^2/8(A_2(N-2))^4+9πτ_1A_3^2B_1N(N-4)^2ε_j^2/4(A_2(N-2))^4+27A_3^2B_1B_3T_5N(N-4)^2ε_j^2/8(A_2(N-2))^4
       +27A_3^3B_1S_2N(N-4)^2ε_j^2/8(A_2(N-2))^4+9A_3^2B_1^2(N-4)^2(∑_mε_m^2)^2/16(A_2(N-2))^4+9πτ_1A_3^2B_1(N-4)^2∑_mε_m^2/4(A_2(N-2))^4
       +27A_3^2B_1B_3T_5(N-4)^2∑_mε_m^2/8(A_2(N-2))^4+27A_3^3B_1S_2(N-4)^2∑_mε_m^2/8(A_2(N-2))^4]
       ×exp[-B_2^3(N-4)^3ε_j^3/6(A_2(N-2))^3+B_2^2C_2N(N-4)^2ε_j^4/2(A_2(N-2))^3+B_2^2C_2(N-4)^2ε_j^2∑_mε_m^2/2(A_2(N-2))^3
       -iπτ_2B_2^2(N-4)^2ε_j^2/(A_2(N-2))^3-3A_3B_1B_2^2N(N-4)^3ε_j^4/4(A_2(N-2))^4-3A_3B_1B_2^2(N-4)^3ε_j^2∑_mε_m^2/4(A_2(N-2))^4
       -3πτ_1A_3B_2^2(N-4)^3ε_j^2/2(A_2(N-2))^4-9A_3B_2^2B_3T_5(N-4)^3ε_j^2/4(A_2(N-2))^4-9A_3^2B_2^2S_2(N-4)^3ε_j^2/4(A_2(N-2))^4]
       exp[B_2^4(N-4)^4ε_j^4/8A_2^4(N-2)^4]×exp[-B_1^2N^2ε_j^4/4A_2(N-2)-B_1^2Nε_j^2∑_mε_m^2/2A_2(N-2)
       -πτ_1B_1Nε_j^2/A_2(N-2)-3B_1B_3T_5Nε_j^2/2A_2(N-2)-3A_3B_1S_2Nε_j^2/2A_2(N-2)
       +B_1C_1N^2ε_j^5/2A_2(N-2)-B_1C_1Nε_j^3∑_mε_m^2/2A_2(N-2)+B_1C_1Nε_j^2∑_mε_m^3/2A_2(N-2)-πτ_3B_1Nε_j^3/A_2(N-2)
       -3B_1D_3T_5Nε_j^3/2A_2(N-2)-iB_1D_2T_4Nε_j^3/A_2(N-2)-3B_1B_3S_2Nε_j^3/2A_2(N-2)-B_1D_1N(N-2)ε_j^6/2A_2(N-2)


    -B_1D_1Nε_j^2∑_mε_m^4/2A_2(N-2)+B_1E_1Nε_j^3∑_mε_m^3/2A_2(N-2)-πτ_4B_1Nε_j^4/A_2(N-2)+3B_1C_3S_2Nε_j^4/2A_2(N-2)
       +B_1F_1Nε_j^4∑_mε_m^2/2A_2(N-2)-iB_1C_2S_1Nε_j^4/A_2(N-2)-B_1^2(∑_mε_m^2)^2/4A_2(N-2)-πτ_1B_1∑_mε_m^2/A_2(N-2)
       -3B_1B_3T_5∑_mε_m^2/2A_2(N-2)-3A_3B_1S_2∑_mε_m^2/2A_2(N-2)+B_1C_1Nε_j^3∑_mε_m^2/2A_2(N-2)
       -B_1C_1ε_j(∑_mε_m^2)^2/2A_2(N-2)+B_1C_1∑_mε_m^2∑_nε_n^3/2A_2(N-2)-πτ_3B_1ε_j∑_mε_m^2/A_2(N-2)
       -3B_1D_3T_5ε_j∑_mε_m^2/2A_2(N-2)-iB_1D_2T_4ε_j∑_mε_m^2/A_2(N-2)-3B_1B_3S_2ε_j∑_mε_m^2/2A_2(N-2)
       -B_1D_1(N-2)ε_j^4∑_mε_m^2/2A_2(N-2)-B_1D_1∑_mε_m^2∑_nε_n^4/2A_2(N-2)+B_1E_1ε_j∑_mε_m^2∑_nε_n^3/2A_2(N-2)
       -πτ_4B_1ε_j^2∑_mε_m^2/A_2(N-2)+3B_1C_3S_2ε_j^2∑_mε_m^2/2A_2(N-2)+B_1F_1ε_j^2(∑_mε_m^2)^2/2A_2(N-2)
       -iB_1C_2S_1ε_j^2∑_mε_m^2/A_2(N-2)-π^2τ_1^2/A_2(N-2)-3πτ_1B_3T_5/A_2(N-2)-3πτ_1A_3S_2/A_2(N-2)
       +πτ_1C_1Nε_j^3/A_2(N-2)-πτ_1C_1ε_j∑_mε_m^2/A_2(N-2)+πτ_1C_1∑_mε_m^3/A_2(N-2)-2π^2τ_1τ_3ε_j/A_2(N-2)-3πτ_1D_3T_5ε_j/A_2(N-2)
       -2iπτ_1D_2T_4ε_j/A_2(N-2)-3πτ_1B_3S_2ε_j/A_2(N-2)-πτ_1D_1(N-2)ε_j^4/A_2(N-2)-πτ_1D_1∑_mε_m^4/A_2(N-2)
       +πτ_1E_1ε_j∑_mε_m^3/A_2(N-2)-2π^2τ_1τ_4ε_j^2/A_2(N-2)+3πτ_1C_3S_2ε_j^2/A_2(N-2)
       +πτ_1F_1ε_j^2∑_mε_m^2/A_2(N-2)-2iπτ_1C_2S_1ε_j^2/A_2(N-2)-9B_3^2T_5^2/4A_2(N-2)-9A_3B_3S_2T_5/2A_2(N-2)
       +3B_3C_1T_5Nε_j^3/2A_2(N-2)-3B_3C_1T_5ε_j∑_mε_m^2/2A_2(N-2)+3B_3C_1T_5∑_mε_m^3/2A_2(N-2)-3πτ_3B_3T_5ε_j/A_2(N-2)
       -9B_3D_3T_5^2ε_j/2A_2(N-2)-3iB_3D_2T_4T_5ε_j/A_2(N-2)-9B_3^2S_2T_5ε_j/2A_2(N-2)-3B_3D_1T_5(N-2)ε_j^4/2A_2(N-2)
       -3B_3D_1T_5∑_mε_m^4/2A_2(N-2)+3B_3E_1T_5ε_j∑_mε_m^3/2A_2(N-2)-3πτ_4B_3T_5ε_j^2/A_2(N-2)+9B_3C_3S_2T_5ε_j^2/2A_2(N-2)
       +3B_3F_1T_5ε_j^2∑_mε_m^2/2A_2(N-2)-3iB_3C_2S_1T_5ε_j^2/A_2(N-2)-9A_3^2S_2^2/4A_2(N-2)
       +3A_3C_1S_2Nε_j^3/2A_2(N-2)-3A_3C_1S_2ε_j∑_mε_m^2/2A_2(N-2)+3A_3C_1S_2∑_mε_m^3/2A_2(N-2)-3πτ_3A_3S_2ε_j/A_2(N-2)
       -9A_3D_3S_2T_5ε_j/2A_2(N-2)-3iA_3D_2S_2T_4ε_j/A_2(N-2)-9A_3B_3S_2^2ε_j/2A_2(N-2)-3A_3D_1S_2(N-2)ε_j^4/2A_2(N-2)
       -3A_3D_1S_2∑_mε_m^4/2A_2(N-2)+3A_3E_1S_2ε_j∑_mε_m^3/2A_2(N-2)-3πτ_4A_3S_2ε_j^2/A_2(N-2)+9A_3C_3S_2^2ε_j^2/2A_2(N-2)


    +3A_3F_1S_2ε_j^2∑_mε_m^2/2A_2(N-2)-3iA_3C_2S_1S_2ε_j^2/A_2(N-2)-C_1^2N^2ε_j^6/4A_2(N-2)+C_1^2Nε_j^4∑_mε_m^2/2A_2(N-2)
       -C_1^2Nε_j^3∑_mε_m^3/2A_2(N-2)+πτ_3C_1Nε_j^4/A_2(N-2)+3C_1D_3T_5Nε_j^4/2A_2(N-2)+iC_1D_2T_4Nε_j^4/A_2(N-2)
       +3B_3C_1S_2Nε_j^4/2A_2(N-2)-C_1^2ε_j^2(∑_mε_m^2)^2/4A_2(N-2)+C_1^2ε_j∑_mε_m^2∑_nε_n^3/2A_2(N-2)-πτ_3C_1ε_j^2∑_mε_m^2/A_2(N-2)
       -3C_1D_3T_5ε_j^2∑_mε_m^2/2A_2(N-2)-iC_1D_2T_4ε_j^2∑_mε_m^2/A_2(N-2)-3B_3C_1S_2ε_j^2∑_mε_m^2/2A_2(N-2)-C_1^2(∑_mε_m^3)^2/4A_2(N-2)
       +πτ_3C_1ε_j∑_mε_m^3/A_2(N-2)+3C_1D_3T_5ε_j∑_mε_m^3/2A_2(N-2)+iC_1D_2T_4ε_j∑_mε_m^3/A_2(N-2)
       +3B_3C_1S_2ε_j∑_mε_m^3/2A_2(N-2)-π^2τ_3^2ε_j^2/A_2(N-2)-3πτ_3D_3T_5ε_j^2/A_2(N-2)
       -2iπτ_3D_2T_4ε_j^2/A_2(N-2)-3πτ_3B_3S_2ε_j^2/A_2(N-2)-9D_3^2T_5^2ε_j^2/4A_2(N-2)-3iD_2D_3T_4T_5ε_j^2/A_2(N-2)
       -9B_3D_3S_2T_5ε_j^2/2A_2(N-2)+D_2^2T_4^2ε_j^2/A_2(N-2)-3iB_3D_2S_2T_4ε_j^2/A_2(N-2)-9B_3^2S_2^2ε_j^2/4A_2(N-2)
       +A_3B_1^3N^3(N-4)ε_j^6/8(A_2(N-2))^3+3A_3B_1^3N^2(N-4)ε_j^4∑_mε_m^2/8(A_2(N-2))^3+3πτ_1A_3B_1^2N^2(N-4)ε_j^4/4(A_2(N-2))^3
    
       +9A_3B_1^2B_3T_5N^2(N-4)ε_j^4/8(A_2(N-2))^3+9A_3^2B_1^2S_2N^2(N-4)ε_j^2/8(A_2(N-2))^3
    
       +3A_3B_1^3N(N-4)ε_j^2(∑_mε_m^2)^2/8(A_2(N-2))^3+3πτ_1A_3B_1^2N(N-4)ε_j^2∑_mε_m^2/2(A_2(N-2))^3
       +9A_3B_1^2B_3T_5N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3+9A_3^2B_1^2S_2N(N-4)ε_j^2∑_mε_m^2/4(A_2(N-2))^3
       +3π^2τ_1^2A_3B_1N(N-4)ε_j^2/2(A_2(N-2))^3+9πτ_1A_3B_1B_3T_5N(N-4)ε_j^2/2(A_2(N-2))^3
       +9πτ_1A_3^2B_1S_2N(N-4)ε_j^2/2(A_2(N-2))^3+27A_3B_1B_3^2T_5^2N(N-4)ε_j^2/8(A_2(N-2))^3
       +27A_3^2B_1B_3S_2T_5N(N-4)ε_j^2/4(A_2(N-2))^3+27A_3^3B_1S_2^2N(N-4)ε_j^2/8(A_2(N-2))^3
       +A_3B_1^3(N-4)(∑_mε_m^2)^3/8(A_2(N-2))^3+3πτ_1A_3B_1^2(N-4)(∑_mε_m^2)^2/4(A_2(N-2))^3
       +9A_3B_1^2B_3T_5(N-4)(∑_mε_m^2)^2/8(A_2(N-2))^3+9A_3^2B_1^2S_2(N-4)(∑_mε_m^2)^2/8(A_2(N-2))^3
       +3π^2τ_1^2A_3B_1(N-4)∑_mε_m^2/2(A_2(N-2))^3+9πτ_1A_3B_1B_3T_5(N-4)∑_mε_m^2/2(A_2(N-2))^3
       +9πτ_1A_3^2B_1S_2(N-4)∑_mε_m^2/2(A_2(N-2))^3+27A_3B_1B_3^2T_5^2(N-4)∑_mε_m^2/8(A_2(N-2))^3
       +27A_3^2B_1B_3S_2T_5(N-4)∑_mε_m^2/4(A_2(N-2))^3+27A_3^3B_1S_2^2(N-4)∑_mε_m^2/8(A_2(N-2))^3]


    ×exp[-15A_3^2(N-4)^2/16(A_2(N-2))^3+45A_3^2B_2(N-4)^3ε_j/16(A_2(N-2))^4-45A_3^2C_2N(N-4)^2ε_j^2/16(A_2(N-2))^4
       -45A_3^2C_2(N-4)^2∑_mε_m^2/16(A_2(N-2))^4+135A_3^3B_1N(N-4)^3ε_j^2/32(A_2(N-2))^5+135A_3^3B_1(N-4)^3∑_mε_m^2/32(A_2(N-2))^5
       -45A_3^2B_2^2(N-4)^4ε_j^2/8(A_2(N-2))^5-15A_3B_3(N-4)(N-8)ε_j/8(A_2(N-2))^3
       +45A_3B_2B_3(N-4)^2(N-8)ε_j^2/8(A_2(N-2))^4+15A_3C_3(N-4)(N-8)ε_j^2/8(A_2(N-2))^3
       +15A_3C_3(N-4)∑_mε_m^2/8(A_2(N-2))^3-15B_3^2(N-8)^2ε_j^2/16(A_2(N-2))^3]
       ×exp[3A_4(N-8)/4(A_2(N-2))^2-3A_4B_2(N-4)(N-8)ε_j/2(A_2(N-2))^3+3A_4C_2N(N-8)ε_j^2/2(A_2(N-2))^3
       +3A_4C_2(N-8)∑_mε_m^2/2(A_2(N-2))^3-9A_3A_4B_1N(N-4)(N-8)ε_j^2/4(A_2(N-2))^4
       -9A_3A_4B_1(N-4)(N-8)∑_mε_m^2/4(A_2(N-2))^4+9A_4B_2^2(N-4)^2(N-8)ε_j^2/4(A_2(N-2))^4
       +3B_4(N-16)ε_j/4(A_2(N-2))^2-3B_2B_4(N-4)(N-16)ε_j^2/2(A_2(N-2))^3+3C_4(N-16)ε_j^2/4(A_2(N-2))^2
       +3C_4∑_mε_m^2/4(A_2(N-2))^2]×(1+𝒪(N^-2))  .

In the product of N of such exponentials yields various sums. These expression simplify somewhat because ∑_mε_m=0. Integrating τ_5 yields a delta function, from which it follows that

    T_5=-B_2(N-4)∑_mε_m^2/2A_2^2(N-2)^2+𝒪(N^-1)  .

This means that T_5=𝒪(N^-1+2ω), so that many of the T_5-terms may be neglected. Using the same strategy it follows that the τ_2-integral yields a delta function, so that

    S_2=N/2A_2(N-2)
       +2A_2C_2N(N-2)+B_2^2(N-4)^2-6A_3B_1N(N-4)/2A_2^3(N-2)^3∑_mε_m^2+𝒪(N^-1) .

Also the τ_4-integral yields a delta function, so that

    T_4=-3iA_3(N-4)∑_mε_m^2/4A_2^2(N-2)^2-iB_1N∑_mε_m^4/2A_2(N-2)-iB_1(∑_mε_m^2)^2/2A_2(N-2)
       -iπτ_1∑_mε_m^2/A_2(N-2)-3iA_3S_2∑_mε_m^2/2A_2(N-2)+𝒪(N^-1)  ,

which implies that all T_4-dependence may be neglected. Alternatively, the T_4-integral and afterwards the τ_4-integral may be computed using the saddle point method (<ref>), yielding the same result. Using this for the T_3-integral yields

    ∫ T_3 exp[2π iτ_3T_3-2B_2S_1T_3-2C_2T_3^2-3iB_1B_4T_3N∑_mε_m^2/A_2^2(N-2)^2]
       =√(π/2C_2)exp[-π^2τ_3^2/2C_2+B_2^2S_1^2/2C_2-iπτ_3B_2S_1/C_2+3πτ_3B_1B_4N∑_mε_m^2/2C_2A_2^2(N-2)^2
       +3iB_1B_2B_4S_1N∑_mε_m^2/2C_2A_2^2(N-2)^2]×(1+𝒪(N^-1))  .

The τ_3-integral is very similar, although a bit longer. We obtain

    ∫τ_3 exp[-π^2τ_3^2/2C_2-iπτ_3B_2S_1/C_2+3πτ_3B_1B_4N∑_mε_m^2/2C_2A_2^2(N-2)^2+9πτ_3A_3B_2(N-4)^2∑_mε_m^2/4A_2^3(N-2)^3
       -3πτ_3B_3(N-8)∑_mε_m^2/2A_2^2(N-2)^2-πτ_3B_1N∑_mε_m^3/A_2(N-2)+πτ_3C_1N∑_mε_m^4/A_2(N-2)-πτ_3C_1(∑_mε_m^2)^2/A_2(N-2)
       -π^2τ_3^2∑_mε_m^2/A_2(N-2)-3πτ_3B_3S_2∑_mε_m^2/A_2(N-2)]
       =√(2C_2/π)√(1/1+2C_2∑_mε_m^2/A_2(N-2))exp[-B_2^2S_1^2/2C_2+B_2^2S_1^2∑_mε_m^2/A_2(N-2)-3iB_1B_2B_4S_1N∑_mε_m^2/2C_2A_2^2(N-2)^2
       -9iA_3B_2^2S_1(N-4)^2∑_mε_m^2/4A_2^3(N-2)^3+3iB_2B_3S_1(N-8)∑_mε_m^2/2A_2^2(N-2)^2+iB_1B_2S_1N∑_mε_m^3/A_2(N-2)
       +B_1^2C_2N^2(∑_mε_m^3)^2/2A_2^2(N-2)^2-iB_2C_1S_1N∑_mε_m^4/A_2(N-2)+iB_2C_1S_1(∑_mε_m^2)^2/A_2(N-2)
       +3iB_2B_3S_1S_2∑_mε_m^2/A_2(N-2)]×(1+𝒪(N^-1))  .

The next intgral:

    ∫ S_1 exp[2π iτ_1S_1-A_2S_1^2+B_2^2S_1^2∑_mε_m^2/A_2(N-2)-9iA_3B_2^2S_1(N-4)^2∑_mε_m^2/4A_2^3(N-2)^3
       +3iB_2B_3S_1(N-8)∑_mε_m^2/2A_2^2(N-2)^2+iB_1B_2S_1N∑_mε_m^3/A_2(N-2)-iB_2C_1S_1N∑_mε_m^4/A_2(N-2)
       +iB_2C_1S_1(∑_mε_m^2)^2/A_2(N-2)+3iB_2B_3S_1S_2∑_mε_m^2/A_2(N-2)+3iC_3S_1∑_mε_m^2/2A_2(N-2)
       -3iA_3C_2S_1(N-4)∑_mε_m^2/2A_2^2(N-2)^2-6iA_4B_1S_1N∑_mε_m^2/A_2^2(N-2)^2-iB_1C_2S_1N∑_mε_m^4/A_2(N-2)
       -iB_1C_2S_1(∑_mε_m^2)^2/A_2(N-2)-2iπτ_1C_2S_1∑_mε_m^2/A_2(N-2)-3iA_3C_2S_1S_2∑_mε_m^2/A_2(N-2)]
       =√(π/A_2)√(1/1-B_2^2∑_mε_m^2/A_2^2(N-2))exp[-π^2τ_1^2/A_2-π^2τ_1^2B_2^2∑_mε_m^2/A_2^3(N-2)-πτ_1B_1B_2N∑_mε_m^3/A_2^2(N-2)
       -B_1^2B_2^2N^2(∑_mε_m^3)^2/4A_2^3(N-2)^2+9πτ_1A_3B_2^2(N-4)^2∑_mε_m^2/4A_2^4(N-2)^3-3πτ_1B_2B_3(N-8)∑_mε_m^2/2A_2^3(N-2)^2
       +πτ_1B_2C_1N∑_mε_m^4/A_2^2(N-2)-πτ_1B_2C_1(∑_mε_m^2)^2/A_2^2(N-2)-3πτ_1B_2B_3S_2∑_mε_m^2/A_2^2(N-2)-3πτ_1C_3∑_mε_m^2/2A_2^2(N-2)
       +3πτ_1A_3C_2(N-4)∑_mε_m^2/2A_2^3(N-2)^2+6πτ_1A_4B_1N∑_mε_m^2/A_2^3(N-2)^2+πτ_1B_1C_2N∑_mε_m^4/A_2^2(N-2)
       +πτ_1B_1C_2(∑_mε_m^2)^2/A_2^2(N-2)+2π^2τ_1^2C_2∑_mε_m^2/A_2^2(N-2)+3πτ_1A_3C_2S_2∑_mε_m^2/A_2^2(N-2)]×(1+𝒪(N^-1))  .

For the final integral we repeat the proces once more. This yields

    ∫τ_1 exp[-π^2τ_1^2/A_2-π^2τ_1^2B_2^2∑_mε_m^2/A_2^3(N-2)-πτ_1B_1B_2N∑_mε_m^3/A_2^2(N-2)
       +9πτ_1A_3B_2^2(N-4)^2∑_mε_m^2/4A_2^4(N-2)^3-3πτ_1B_2B_3(N-8)∑_mε_m^2/2A_2^3(N-2)^2
       +πτ_1B_2C_1N∑_mε_m^4/A_2^2(N-2)-πτ_1B_2C_1(∑_mε_m^2)^2/A_2^2(N-2)-3πτ_1B_2B_3S_2∑_mε_m^2/A_2^2(N-2)-3πτ_1C_3∑_mε_m^2/2A_2^2(N-2)
       +3πτ_1A_3C_2(N-4)∑_mε_m^2/2A_2^3(N-2)^2+6πτ_1A_4B_1N∑_mε_m^2/A_2^3(N-2)^2+πτ_1B_1C_2N∑_mε_m^4/A_2^2(N-2)
       +πτ_1B_1C_2(∑_mε_m^2)^2/A_2^2(N-2)+2π^2τ_1^2C_2∑_mε_m^2/A_2^2(N-2)+3πτ_1A_3C_2S_2∑_mε_m^2/A_2^2(N-2)]
       ×exp[-3πτ_1A_3N(N-4)/2A_2^2(N-2)^2-6πτ_1A_3C_2N(N-4)∑_mε_m^2/A_2^3(N-2)^3
       -15πτ_1A_3B_2^2(N-4)^3∑_mε_m^2/4A_2^4(N-2)^4+9πτ_1A_3^2B_1N(N-4)^2∑_mε_m^2/A_2^4(N-2)^4


    +9πτ_1B_2B_3(N-4)(N-8)∑_mε_m^2/4A_2^3(N-2)^3+3πτ_1C_3(N-4)∑_mε_m^2/A_2^2(N-2)^2
       -2πτ_1B_1N∑_mε_m^2/A_2(N-2)-π^2τ_1^2N/A_2(N-2)-3πτ_1B_3T_5N/A_2(N-2)-3πτ_1A_3S_2N/A_2(N-2)+2πτ_1C_1N∑_mε_m^3/A_2(N-2)
       -2πτ_1D_1(N-1)∑_mε_m^4/A_2(N-2)+3πτ_1C_3S_2∑_mε_m^2/A_2(N-2)+πτ_1F_1(∑_mε_m^2)^2/A_2(N-2)
       +3πτ_1A_3B_1^2N^2(N-4)∑_mε_m^4/4A_2^3(N-2)^3+3π^2τ_1^2A_3B_1N(N-4)∑_mε_m^2/A_2^3(N-2)^3
       +9πτ_1A_3B_1^2N(N-4)(∑_mε_m^2)^2/4A_2^3(N-2)^3+9πτ_1A_3^2B_1S_2N(N-4)∑_mε_m^2/A_2^3(N-2)^3]
       =√(A_2(N-2)/2π(N-1))√(1/1-(-B_2^2/2A_2^2(N-1)+C_2/A_2(N-1)+3A_3B_1N(N-4)/2A_2^2(N-1)(N-2)^2)∑_mε_m^2)×exp[3A_3B_1N^2∑_mε_m^2/2A_2^2(N-1)(N-2)
       +3A_3B_1N^2(∑_mε_m^2)^2/2A_2^2(N-1)(N-2)(-B_2^2/2A_2^2(N-1)+C_2/A_2(N-1)+3A_3B_1N(N-4)/2A_2^2(N-2)^2(N-1))+B_1^2N^2(∑_mε_m^2)^2/2A_2(N-1)(N-2)
       +B_1^2N^2(∑_mε_m^2)^3/2A_2(N-1)(N-2)(-B_2^2/2A_2^2(N-1)+C_2/A_2(N-1)+3A_3B_1N(N-4)/2A_2^2(N-2)^2(N-1))+9A_3^2N^2/8A_2^3(N-1)(N-2)
    
       +B_1^2B_2N^2(∑_mε_m^2)(∑_nε_n^3)/2A_2^2(N-1)(N-2)+3A_3B_1B_2N^2∑_mε_m^3/4A_2^3(N-1)(N-2)+B_1^2B_2^2N^2(∑_mε_m^3)^2/8A_2^3(N-1)(N-2)
       +C_1^2N^2(∑_mε_m^3)^2/2A_2(N-1)(N-2)-B_1C_1N^2(∑_mε_m^2)(∑_nε_n^3)/A_2(N-1)(N-2)-3A_3C_1N^2∑_mε_m^3/2A_2^2(N-1)(N-2)
       +3A_3B_1N^2(∑_mε_m^2)^2/4A_2^4(N-1)(N-2)^4(2A_2C_2N(N-2)+B_2^2(N-4)^2-6A_3B_1N(N-4))
       +9A_3^2N^2∑_mε_m^2/8A_2^5(N-1)(N-2)^4(2A_2C_2N(N-2)+B_2^2(N-4)^2-6A_3B_1N(N-4))
       -9A_3B_1B_2^2N(N-4)^2(∑_mε_m^2)^2/8A_2^4(N-1)(N-2)^3-27A_3^2B_2^2N(N-4)^2∑_mε_m^2/16A_2^5(N-1)(N-2)^3+3B_1B_2B_3N(N-8)(∑_mε_m^2)^2/4A_2^3(N-1)(N-2)^2
       +9A_3B_2B_3N(N-8)∑_mε_m^2/8A_2^4(N-1)(N-2)^2-B_1B_2C_1N^2(∑_mε_m^2)(∑_nε_n^4)/2A_2^2(N-1)(N-2)-3A_3B_2C_1N^2∑_mε_m^4/4A_2^3(N-1)(N-2)
       +B_1B_2C_1N(∑_mε_m^2)^3/2A_2^2(N-1)(N-2)+3A_3B_2C_1N(∑_mε_m^2)^2/4A_2^3(N-1)(N-2)+3B_1B_2B_3S_2N(∑_mε_m^2)^2/2A_2^2(N-1)(N-2)
       +9A_3B_2B_3S_2N∑_mε_m^2/4A_2^3(N-1)(N-2)+3B_1C_3N(∑_mε_m^2)^2/4A_2^2(N-1)(N-2)+9A_3C_3N∑_mε_m^2/8A_2^3(N-1)(N-2)
       -3A_3B_1C_2N(N-4)(∑_mε_m^2)^2/4A_2^3(N-1)(N-2)^2-9A_3^2C_2N(N-4)∑_mε_m^2/8A_2^4(N-1)(N-2)^2-3A_4B_1^2N^2(∑_mε_m^2)^2/A_2^3(N-1)(N-2)^2
       -9A_3A_4B_1N^2∑_mε_m^2/2A_2^4(N-1)(N-2)^2-B_1^2C_2N^2(∑_mε_m^2)(∑_nε_n^4)/2A_2^2(N-1)(N-2)-3A_3B_1C_2N^2∑_mε_m^4/4A_2^3(N-1)(N-2)


    -B_1^2C_2N(∑_mε_m^2)^3/2A_2^2(N-1)(N-2)-3A_3B_1C_2N(∑_mε_m^2)^2/4A_2^3(N-1)(N-2)-3A_3B_1C_2S_2N(∑_mε_m^2)^2/2A_2^2(N-1)(N-2)
       -9A_3^2C_2S_2N∑_mε_m^2/4A_2^3(N-1)(N-2)+3A_3B_1C_2N^2(N-4)(∑_mε_m^2)^2/A_2^3(N-1)(N-2)^3+9A_3^2C_2N^2(N-4)∑_mε_m^2/2A_2^4(N-1)(N-2)^3
       +15A_3B_1B_2^2N(N-4)^3(∑_mε_m^2)^2/8A_2^4(N-1)(N-2)^4+45A_3^2B_2^2N(N-4)^3∑_mε_m^2/16A_2^5(N-1)(N-2)^4-9A_3^2B_1^2N^2(N-4)^2(∑_mε_m^2)^2/2A_2^4(N-1)(N-2)^4
       -27A_3^3B_1N^2(N-4)^2∑_mε_m^2/4A_2^5(N-1)(N-2)^4-9B_1B_2B_3N(N-4)(N-8)(∑_mε_m^2)^2/8A_2^3(N-1)(N-2)^3
       -27A_3B_2B_3N(N-4)(N-8)∑_mε_m^2/16A_2^4(N-1)(N-2)^3-3B_1C_3N(N-4)(∑_mε_m^2)^2/2A_2^2(N-1)(N-2)^2-9A_3C_3N(N-4)∑_mε_m^2/4A_2^3(N-1)(N-2)^2
       +3B_1B_3T_5N^2∑_mε_m^2/2A_2(N-1)(N-2)+9A_3B_3T_5N^2/4A_2^2(N-1)(N-2)+B_1D_1N(N-1)(∑_mε_m^2)(∑_nε_n^4)/A_2(N-1)(N-2)
       +3A_3D_1N(N-1)∑_mε_m^4/2A_2^2(N-1)(N-2)-3B_1C_3S_2N(∑_mε_m^2)^2/2A_2(N-1)(N-2)-9A_3C_3S_2N∑_mε_m^2/8A_2^2(N-1)(N-2)
       -B_1F_1N(∑_mε_m^2)^3/2A_2(N-1)(N-2)-3A_3F_1N(∑_mε_m^2)^2/4A_2^2(N-1)(N-2)-3A_3B_1^3N^3(N-4)(∑_mε_m^2)(∑_nε_n^4)/8A_2^3(N-1)(N-2)^3
       -9A_3^2B_1^2N^3(N-4)∑_mε_m^4/16A_2^4(N-1)(N-2)^3-9A_3B_1^3N^2(N-4)(∑_mε_m^2)^3/8A_2^3(N-1)(N-2)^3-27A_3^2B_1^2N^2(N-4)(∑_mε_m^2)^2/16A_2^4(N-1)(N-2)^3
       -9A_3^2B_1^2S_2N^2(N-4)(∑_mε_m^2)^2/2A_2^3(N-1)(N-2)^3-27A_3^3B_1S_2N^2(N-4)∑_mε_m^2/4A_2^4(N-1)(N-2)^3]×(1+𝒪(N^-1))  ,

where (<ref>) is applied to to exp[-3πτ_1A_3S_2N/A_2(N-2)] to simplify the notation a bit.


Besides this, the expansion

    [∏_j(1+1/λ_j)^λ(N-1)+(N-2)ε_j/2/2][∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)]
       =(1+1/λ)^x/2(1+λ)^N2exp[N-1/2∑_j(λ+ε_j)log(1+λ+ε_j/1+λλ/λ+ε_j)]
       ×exp[-∑_k<llog(1+λ-√((1+λ-1+λ/1+λ+ε_k)(1+λ-1+λ/1+λ+ε_l)))]
       =(1+1/λ)^x/2(1+λ)^N2
       ×exp[(∑_mε_m^2)·[Nλ^2/4λ^2(λ+1)^2+λ/4λ^2(λ+1)^2-(N-1)3λ^2+λ/8λ^2(λ+1)^2-1/8λ(λ+1)]
       +(∑_mε_m^3)·[-N6λ^3+3λ^2+λ/24λ^3(λ+1)^3-3λ^2+λ/12λ^3(λ+1)^3
       +(N-1)14λ^3+9λ^2+3λ/48λ^3(λ+1)^3+2λ+1/16λ^2(λ+1)^2]
       +(∑_mε_m^4)·[N6λ^4+6λ^3+4λ^2+λ/24λ^4(λ+1)^4+6λ^3+4λ^2+λ/24λ^4(λ+1)^4
       -(N-1)30λ^4+28λ^3+19λ^2+5λ/128λ^4(λ+1)^4-2λ^2+2λ+1/32λ^3(λ+1)^3-6λ^2+6λ+1/128λ^3(λ+1)^3]
       +N(∑_mε_m^5)[-20λ^5+30λ^4+30λ^3+15λ^2+3λ/80λ^5(λ+1)^5+248λ^5+300λ^4+310λ^3+165λ^2+35λ/1280λ^5(λ+1)^5]
       +N(∑_mε_m^6)[15λ^6+30λ^5+40λ^4+30λ^3+12λ^2+2λ/60λ^6(λ+1)^6-504λ^6+744λ^5+1038λ^4+836λ^3+357λ^2+63λ/3072λ^6(λ+1)^6]
       +(∑_mε_m^2)^2·[6λ^2+6λ+1/128λ^3(λ+1)^3]
       +(∑_mε_m^2)(∑_nε_n^3)·[-(2λ+1)^3/128λ^4(λ+1)^4]
       +(∑_mε_m^2)(∑_nε_n^4)·[40λ^4+80λ^3+74λ^2+34λ+5/1024λ^5(λ+1)^5]
       +(∑_mε_m^3)^2·[40λ^4+80λ^3+58λ^2+18λ+3/768λ^5(λ+1)^5]]×(1+𝒪(N^-1))

is needed. Together with all the exponents the final answer

    V_N(t;λ)=√(2)(1+λ)^N2/(2πλ(λ+1)N)^N/2(1+1/λ)^x/2exp[14λ^2+14λ-1/12λ(λ+1)]
       ×exp[(∑_mε_m^2)·[Nλ^2/4λ^2(λ+1)^2+λ/4λ^2(λ+1)^2-(N-1)3λ^2+λ/8λ^2(λ+1)^2-1/8λ(λ+1)
       +32λ^4+64λ^3+40λ^2+8λ+1/256λ^4(λ+1)^4-264λ^6+792λ^5+784λ^4+248λ^3-33λ^2-25λ-3/192λ^4(λ+1)^4N]]
       ×exp[(∑_mε_m^3)·[-N6λ^3+3λ^2+λ/24λ^3(λ+1)^3-3λ^2+λ/12λ^3(λ+1)^3
       +(N-1)14λ^3+9λ^2+3λ/48λ^3(λ+1)^3+2λ+1/16λ^2(λ+1)^2
       -2λ^3+3λ^2+3λ+1/48λ^3(λ+1)^3]]
       ×exp[(∑_mε_m^4)·[N6λ^4+6λ^3+4λ^2+λ/24λ^4(λ+1)^4+6λ^3+4λ^2+λ/24λ^4(λ+1)^4
       -(N-1)30λ^4+28λ^3+19λ^2+5λ/128λ^4(λ+1)^4-2λ^2+2λ+1/32λ^3(λ+1)^3-6λ^2+6λ+1/128λ^3(λ+1)^3
       -(2λ+1)^2N/128λ^3(λ+1)^3-6λ^4+12λ^3+20λ^2+14λ+3/128λ^4(λ+1)^4]]
       ×exp[N(∑_mε_m^5)·[-20λ^5+30λ^4+30λ^3+15λ^2+3λ/80λ^5(λ+1)^5+248λ^5+300λ^4+310λ^3+165λ^2+35λ/1280λ^5(λ+1)^5
       +4λ^3+6λ^2+4λ+1/128λ^4(λ+1)^4]]
       ×exp[N(∑_mε_m^6)[15λ^6+30λ^5+40λ^4+30λ^3+12λ^2+2λ/60λ^6(λ+1)^6
       -504λ^6+744λ^5+1038λ^4+836λ^3+357λ^2+63λ/3072λ^6(λ+1)^6
       -14λ^4+28λ^3+36λ^2+22λ+5/768λ^5(λ+1)^5]]
       ×exp[(∑_mε_m^2)^2·[6λ^2+6λ+1/128λ^3(λ+1)^3+72λ^6+216λ^5+220λ^4+80λ^3-2λ^2-6λ-1/512λ^5(λ+1)^5N]]
       ×exp[(∑_mε_m^2)(∑_nε_n^3)·[-(2λ+1)^3/128λ^4(λ+1)^4+(2λ+1)^3/128λ^4(λ+1)^4]]
       ×exp[(∑_mε_m^2)(∑_nε_n^4)·[40λ^4+80λ^3+74λ^2+34λ+5/1024λ^5(λ+1)^5
       -16λ^4+32λ^3+44λ^2+28λ+5/1024λ^5(λ+1)^5]]
       ×exp[(∑_mε_m^3)^2·[40λ^4+80λ^3+58λ^2+18λ+3/768λ^5(λ+1)^5-8λ^4+16λ^3-4λ^2-12λ-3/1024λ^5(λ+1)^5]]
       ×exp[(∑_mε_m^2)^3-8λ^4-16λ^3+8λ+1/3072λ^5(λ+1)^5]×(1+𝒪(N^-1))

is obtained. Applying to this that

    e_j=2/N-2(t_j-λ(N-1))

yields the claimed result.

To determine the error from the difference 𝒟 from Lemma <ref>, we divide it by V_N(t;λ). Assuming that t_j-λ(N-1)=λ N^1/2+ω takes maximal values, it follows that the relative difference is given by

    𝒪(N^2-Kα)exp[4λ^2+4λ+1/3λ(λ+1)]
       ×exp[N^2ω/64λ^2(λ+1)^4(16(λ+1)^2λ^2(4λ^2+4λ+3)
       +64λ^6+192λ^5+160λ^4-40λ^2-8λ-1)]
       ×exp[(4λ^2+4λ+1)λ N^4ω/8(λ+1)^3]exp[λ N^6ω/48(λ+1)^5(8λ^4+16λ^3-8λ-1)]  .

Only the first exponential can become large, if λ is small. Assuming that λ > C/log(N), this factor adds an error N^1/3C.

To keep this relative error small, it is furthermore necessary that exp[N^6ω]<N^Kα-2. Solving this yields

    0<ω<log(Kα-2)+log(log(N))/6log(N)  .

Fixing the value ω=ρlog(Kα-2)+log(log(N))/6log(N) for some ρ∈(0,1) shows that the relative error is 𝒪(N^2+1/3C-KαN^(-2+Kα)^ρ)=𝒪(N^2+1/3C-Kα).



Methods to treat such multi-dimensional combinatorical Gaussian integrals in more generality have been discussed in <cit.>.



§ REDUCTION OF THE INTEGRATION REGION


In the previous paragraph the result of the integral (<ref>) in a small box around the origin was obtained. Knowing this makes it much easier to compare the contribution inside and outside of this box. This is the main aim of Lemma <ref>.


For a∈[0,1] and n∈ℕ the estimates

    exp[nalog (2)]≤ (1+a)^n≤exp[na]

hold.


The right-hand side follows from

    (1+a)^n=∑_j=0^na^jnj=∑_j=0^n(na)^j/j!n!/n^j(n-j)!≤∑_j=0^n(na)^j/j!≤exp[na]  .

For the left-hand side it suffices to show that log(1+a)≥ alog(2). Because equality holds at one and zero, this follows from the concavity of the logarithm.



For any ω∈(0,1/4) and α∈(0,1/4-ω), define

    δ_N=N^-αζ_N/min{λ_j}  ,

such that ζ_N→∞ and N^-δζ_N→0 for any δ>0. Assuming that x=∑_jt_j=λ N(N-1), |t_j-λ(N-1)|≪λ N^1/2+ω and λ>C/log(N), the integral

    V_N(t)=(∏_j=1^N(1+1/λ_j)^t_j/2)(2π)^-N∫_𝕋^Nφ e^-i∑_j=1^Nφ_jt_j
       ×∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) 1/1-√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) (e^i(φ_k+φ_l)-1)

can be restricted to

    V_N(t;λ)=2/(2π)^N(∏_j=1^N(1+1/λ_j)^t_j/2)∫_[-δ_N,δ_N]^Nφ exp[-i∑_jφ_j(t_j-λ(N-1))]
       ×∏_1≤ k<l≤ N√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) 1/1-√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l) (e^i(φ_k+φ_l)-1)
       ×(1+𝒪(N^3/2exp[-N^1-2αζ_N^2]))  .




The idea of the proof is to consider the integrand in a small box [-δ_N,δ_N]^N and see what happens to it if some of the angles φ lie outside of it.

Because x is even, it follows that the integrand takes the same value at φ and φ+π=(φ_1+π,…,φ_N+π). This means that only half of the space has to be considered and the result must be multiplied by 2.


This estimate follows directly from application of (<ref>-<ref>) to the integrand and a computation like the one in the proof of Lemma <ref>. Writing

    μ_kl=√(λ_kλ_l)/√((1+λ_k)(1+λ_l))-√(λ_kλ_l)    and    ε_j=λ_j-λ

with |ε_j|≪λ N^-1/2+ω this yields

    |∫_[-δ_N,δ_N]^Nφ ∏_1≤ k<l≤ N1/1-μ_kl(exp[i(φ_k+φ_l)]-1)|
       ≤∫_[-δ_N/2,δ_N/2]^Nφ |exp[∑_m=1i^m∑_k<lA_m(μ_kl)(φ_k+φ_l)^m]|
       ≤√(2)(2π/λ(λ+1)N)^N/2exp[10λ^2+10λ+1/4λ(λ+1)]exp[N^1/2+2ω]  .

The final exponent exp[N^1/2+2ω] here comes from the estimate μ_kl≥λ(1-N^-1/2+ω).


Now we argue case by case why other configurations of the angles φ_j are asymptotically suppressed.

Case 1. All but finitely many angles lie in the box [-δ_N,δ_N]^N. A finite number of m angles lies outside of it. We label these angles {φ_1,…,φ_m}. The maximum of the integrand

    f:(φ_m+1,…,φ_N)↦∏_1≤ k<l≤ N1/1-μ_kl(exp[i(φ_k+φ_l)]-1)

in absolute value is given by the equations

    0=∂_φ_j|f|=∑_k≠ jsin(φ_j+φ_k)/1-2μ_kl(μ_kl+1)(cos(φ_j+φ_k)-1)     for j=m+1,…,N  .

It is clear that the maximum is found for φ̃=φ_m+1=…=φ_N. The first order solution to this is then

    φ̃=-1/2(N-m-1)∑_k=1^msin(φ_k)/1+2μ_kj(μ_kj+1)(1-cosφ_k)  .

This shows that the maximum will lie in the box [-δ_N/2,δ_N/2]^N. This implies that |φ_j-φ_k|>δ_N/2, when 1≤ j≤ m and m+1≤ k≤ N. Applying the estimate (<ref>) to pairs of such angles and afterwards (<ref>) to the remaining N-m angles in the box [-δ_N,δ_N]^N-m gives us an upper bound of

    2√(2)/(2πλ(λ+1)(N-m))^N-m/2(∏_j(1+1/λ_j)^t_j/2)·(∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))
       ×Nmexp[30λ^2+30λ+3/12λ(λ+1)]exp[N^1/2+2ω](1+λ(λ+1)δ_N^2/4)^-Nm/2

on the part of the integral in the small box [-δ_N,δ_N]^N. There are Nm ways to select the m angles. Applying Lemma <ref> to the final factor and comparing the result with the Lower bound, shows that this may be neglected if

    2√(2)e^16λ^2+16λ+4/12λ(λ+1)e^m/2N^3m/2(2πλ(λ+1))^m/2exp[N^1-2α+N^1/2+2ω]e^-Nmlog(2)/8λ(λ+1)δ_N^2→0  .

The condition 0<α<1/4-ω and the sequence ζ_N→∞ guarantee this. In fact, the same argument works for all m such that m/N→0.


Case 2. If the number m=ρ N of angles outside the integration box [-δ_N,δ_N]^N increases faster, another estimate is needed, because the maximum φ̃ may lie outside of [-δ_N/2,δ_N/2]. It is clear that 0<ρ<1 in the limit.

Estimate the location φ_j=φ̃ of the maximum is much trickier now. Regardless of its precise location, we will take the maximum value as the estimate for the integrand in the entire integration box. The smaller box [-δ_N/2,δ_N/2]^N is considered once more. We distinguish two options.

-Case 2a. The maximum lies in [-δ_N/2,δ_N/2]^N, thus φ̃∈[-δ_N/2,δ_N/2].

Applying the estimate (<ref>) to this yields an upper bound

    Nρ N(2δ_N)^N(1-ρ)(2π)^ρ N(∏_j(1+1/λ_j)^t_j/2)
       ×(∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))(1+1/4λ(λ+1)δ_N^2)^-N^2ρ(1-ρ)/4  .

Applying Lemma <ref> to the last factor and dividing this by ℰ_α shows that

    (2^1/ρπδ^1-ρ/ρ(2πλ(λ+1)N)^2/ρNexp[N^-2α/ρ+N^-1/2+2ω(1-ρ)/ρ-λ(λ+1)δ_N^2N(1-ρ)log(2)/16])^ρ N→0

is a sufficient and satisfied condition.

-Case 2b. The maximum lies not in [-δ_N/2,δ_N/2]^N. This is the same as δ_N/2<|φ̃|≤δ_N.

Applying (<ref>) only to the angles φ_ρ N+1,…,φ_ρ N in the integration box gives an upper bound

    Nρ N(2δ_N)^N(1-ρ)(2π)^ρ N(∏_j(1+1/λ_j)^t_j/2)
       ×(∏_k<l√((1+λ_k)(1+λ_l))/√((1+λ_k)(1+λ_l))-√(λ_kλ_l))(1+1/4λ(λ+1)δ_N^2)^-N^2(1-ρ)^2/4  .

The same steps as in Case 2a. will do.

This shows that the integration can be restricted to the box [-δ_N,δ_N]^N. The error terms follow from Case 1., since convergence there is much slower.


Lemma <ref> shows that for every α∈(0,1/4-ω) and N∈ℕ there is a box that contains most of the integral's mass. As N increases, this box shrinks and the approximation becomes better. The parameter α determines how fast this box shrinks. Smaller values of α lower the Lower bound and, hence, increase the number of configurations within reach at the price of more intricate integrals and less accuracy.


The observation that ζ_N=log(N) and K≥α^-1(3/2+1/3C-6ω) satisfies all the demands proves Theorem <ref>.

An idea of the accuracy of these formulas can be obtained from Table <ref> and <ref>, where the reference values

    y_k=∑_j=1^N(t_j-λ(N-1))^k    for k≥ 2

are defined to compare configurations with the reference values 2^-kλ^kN^1+k/2 for k≥ 2.


 









§ POLYTOPES

In the previous paragraphs the asymptotic counting of symmetric matrices with zero diagonal and entries in the natural numbers was discussed. This allows us to return to the polytopes. The first step is to count the total number of symmetric matrices with zero diagonal and integer entries summing up to x to see which fraction of such matrices are covered by Theorem <ref>. This is easily done by a line of N2+x/2 elements, for example unit elements 1, and N2-1 semicolons. Putting the semicolons between the elements, such that the line begins and ends with a unit element and no semicolons stand next to each other, creates such a matrix. The number of elements before the first semicolon minus one is the first matrix element b_12. The number of elements minus one between the first and second semicolon yields the second matrix element b_13. In this way, we obtain the N2 elements of the upper triangular matrix. There are N2-1+x/2 positions to put N2-1 semicolons and thus

    N2-1+x/2N2-1≈1/N√(1/πλ(λ+1))(1+λ)^N2(1+1/λ)^x/2(1+𝒪(N^-1))

such matrices, where we have used Stirling's approximation and the average matrix entry condition λ=x/(N(N-1)) for the approximation.


The next step is to estimate the number of matrices within reach of Theorem <ref>. Using only the leading order, the number of covered matrices is given by

    ∫t V_N(t;λ) δ(λ N(N-1)-∑_jt_j)
       =(1+λ)^N2(1+1/λ)^x/2/π^N/2√(λ (λ+1)N)exp[14λ^2+14λ-1/12λ(λ+1)]∫ S∫σ∫τ exp[2π iσ S+S^2]
       ×{∏_j=1^N∫_-N^ω^N^ω y_j exp[2π iτ y_j]exp[-y_j^2(1+2/N-2π iσ/N)]exp[√(2)(2λ+1)y_j^3/3√(λ(λ+1)N)]
       ×exp[-3λ^2+3λ+1/3λ(λ+1)Ny_j^4]}
       =(1+λ)^N2(1+1/λ)^x/2/N√(πλ (λ+1))exp[14λ^2+14λ-1/12λ(λ+1)]∫ S∫σ∫τ exp[2π iσ S+S^2]
       ×exp[-π^2τ^2(1-2/N+2iπσ/N)]exp[-1+π i σ+π iτ(2λ+1)/√(2λ(λ+1))]exp[5(2λ+1)^2/24λ (λ+1)]
       ×exp[-3λ^2+3λ+1/4λ(λ+1)]×(1-𝒪(exp[-N^2ω]/N^2ω))^N
    
       =(1+λ)^N2(1+1/λ)^x/2/N√(πλ (λ+1))exp[-1/4λ(λ+1)]×(1-𝒪(exp[-N^2ω]/N^2ω))^N×(1+𝒪(N^-1))  .


Here we have used that
beginalign*
   |hspace-4mm∫y exp[ay-by^2+cy^3+dy^4]=exp[a^2/4b+a^3c/8b^3+9a^4c^2/64b^5+a^4d/16b^4]√(π/b-3ac/2b-9a^2c^2/8b^3-3da^2/2b^2)×(1+15c^2/16b^3+3d/4b^2)  ,

assuming that a,b=𝒪(1), c=𝒪(N^-1/2) and d=𝒪(N^-1).


A fraction exp[-1/4λ(λ)] of the matrices is covered, provided that ω is large enough. A sufficient condition is that 

    ω≥loglog N/2log N  .

Combining this with the condition 

    ω≤log (Kα-2) + log(log N)/4log (N)

shows that K≥log(N)/α+2  is necessary to satisfy both demands. However, such large values of K remain without consequences, because higher values of K only influence the the error term in Lemma <ref>.

As λ→∞, the fraction of covered matrices tends to one and the volume of the diagonal subpolytopes of symmetric stochastic matrices can be determined by (<ref>). In terms of the variables

    t_j=1-h_j/a  and  χ=∑_jh_j

the volume of the diagonal subpolytope is calculated by

    (P_N(h))=lim_a→0a^N(N-3)/2V_N(1-h/a;N-χ/aN(N-1))
       =√(2)e^7/6(e(N-χ)/N(N-1))^N2(N(N-1)^2/2π(N-χ)^2)^N/2exp[-N(N-1)^2/2(N-χ)^2∑_j(h_j-χ/N)^2]
       ×exp[-(N-1)^2/(N-χ)^2∑_j(h_j-χ/N)^2]exp[-N(N-1)^3/3(N-χ)^3∑_j(h_j-χ/N)^3]
       ×exp[-N(N-1)^4/4(N-χ)^4∑_j(h_j-χ/N)^4]exp[(N-1)^4/4(N-χ)^4(∑_j(h_j-χ/N)^2)^2]  .


The convergence criterion becomes

    |t_j-λ (N-1)|/λ N^1/2+ω=N^1/2-ω(N-1)/N-χ|h_j-χ/N|→0  .

This is the same as

    ∑_j|h_j-χ/N|^k≪ (N-χ/N-1)^kN^1-k/2+kω  for all k≥2  .


This means that we only have accuracy in a small neighbourhood around χ/N. However, the calculation (<ref>) shows that this corresponds to almost all matrices asymptotically, so that outside of this region the polytopes will have very small volumes. There, not all relevant factors are known, but missing factors will be small compared to the dominant factor. This means that for diagonals that satisfy

    lim_N→∞(N-1)^2∑_j(h_j-χ/N)^2/(N-χ)^2log(N)=0

qualitatively reasonable results are expected.


Since we are calculating a N2-dimensional volume with only one length scale, it follows that no correction can become large in this limit. It inherits the relative error from Theorem <ref>. This proves Theorem <ref>. Examples of this formula at work are given in Figure <ref> and <ref>. 



 §.§ Acknowledgments

This work was supported by the Deutsche Forschungsgemeinschaft (SFB 878 - groups, geometry & actions). We thank N. Broomhead and L. Hille for valuable discussions. In particular, we thank B.D. McKay and M. Isaev for clarifications of various aspects in asymptotic enumeration, pounting out a gap in the previous version and their suggestions to fill it.


unsrt

