Virial expansion of the nonlinear sigma model in the strong coupling limit

One of the most powerful approaches to study complex quantum systems is based on the field-theoretical formulation. The nonlinear σ-model, as an effective field theory describing long-distance physics, was first derived by Wegner [1, 2] in the context of disordered quantum systems using the replica trick. An alternative formulation of the σ-model given by Efetov [3, 4] employed the supersymmetry method making use of both commuting and anticommuting variables. The latter version of the σ-model can be successfully applied not only for problems, which can be solved on the perturbative level, but allows also for a non-perturbative treatment of the theory. Moreover, the supersymmetric σ-model can be derived as an effective field theory for various random matrix ensembles. In particular, it was used to prove that the statistical properties of the energy levels in small metallic disordered grains can be described by the Wigner–Dyson random matrix theory (RMT) [3, 4].

The partition function of the generic supersymmetric nonlinear σ-model on a lattice

$$\begin{equation} Z=\int DQ\; {\rm e}^{-\frac{1}{t}\sum _{m\ne k}U_{mk} {\mathrm{Str}}Q_mQ_k-\sum _{n}S_0[Q_n]} \end{equation} \tag{ 1 }$$

is defined in terms of supermatrices Q_m associated with the mth lattice site, a symmetric coupling matrix U_mk > 0 and a coupling constant t. S₀ comprises all the local terms of the action including source terms. An explicit parametrization of Q and of the measure DQ = ∏_mdQ_m is given in appendix A.

An important feature of the σ-model is the fact that the super-matrix field Q satisfies a nonlinear constraint Q² = 1. It is this constraint that makes generally the analysis of the σ-model highly non-trivial. One possible way to tackle the problem is to use the perturbation theory of the nonlinear σ-model in the metallic regime. The latter is characterized by a small value of the coupling constant t, which corresponds to high electrical conductivity in the case of disordered conductors. In this regime the statistical properties of energy levels and eigenfunctions are close to the corresponding ones of the Wigner–Dyson RMT. A large number of important and interesting results have been obtained following this route (see [5] and references therein).

The opposite limit of the strong coupling has received much less attention. One of the reasons for that is the fact that in many cases the mapping between original microscopic models and the σ-model can be justified only at small values of the coupling constant. At the same time it is possible to construct microscopic models (such as for example n-orbital Wegner model [6]), which can be mapped onto the σ-model at all values of the coupling constant. This situation can be realized, for example, in a cluster of coupled disordered grains or chaotic systems, provided that each element of a cluster is well described by the Wigner–Dyson RMT. In such a case the properties of the σ-model in the weak and in the strong coupling limits are equally important. Moreover, some properties of the σ-model in the strong coupling limit might be of interest independently of whether or not it can be directly derived from a microscopic model. This is especially true for critical exponents of the σ-models which remain critical at large values of the coupling constant, as critical exponents are the same for all microscopic models belonging to the same universality class [7].

Another reason that makes the strong coupling limit less explored is that technically this problem is more challenging, as the standard perturbation theory approach can not be applied there. This is in contrast to the weak coupling limit, where after separating out slow and fast modes, the latter can be integrated out order by order using the Wick theorem.

In this paper we develop a perturbation theory for the supersymmetric nonlinear σ-model defined on a lattice in the strong coupling limit. The general idea of our approach is similar to the one of the virial expansion methods introduced in [8] for studying almost-diagonal random matrices. After presenting a construction of the virial expansion for the generic σ-model, characterized by an arbitrary coupling matrix U, we calculate explicitly the first non-trivial terms of the expansion for two different quantities: the moments of the eigenfunctions I_q and the two-level correlation function R₂(ω). Then we compare the obtained results with the corresponding expressions derived for almost-diagonal random matrices. Our main conclusion is that the expressions for I_q and R₂(ω) are in one-to-one correspondence with the random matrix results. This fact can be traced back to the properties of the target space of Q and its counterpart discussed in [9]. We show that such a correspondence can be established at any order of the virial expansion under very general assumptions.

The layout of this paper is as follows. In section 2, we present the general framework for the virial expansion formulated for the generic nonlinear σ-model. Section 3 describes the calculations of the moments of the eigenfunctions I_q and the distribution function of the eigenfunction intensities P(z) in the lowest non-trivial order of the coupling constant. Similar calculations for the two-level correlation function R₂(ω) are presented in section 4. In section 5 we apply the general formulas to the power-law nonlinear σ-model and compare the obtained expression with the results known for the almost-diagonal random matrices. Finally, in section 6 we discuss the origin of the results and demonstrate the equivalence between the σ-model description and the supersymmetric approach developed in [8] in the strong coupling limit.

The main idea of the virial expansion method is to represent the 'interaction' part of the action $\frac{1}{t}\sum _{m\ne k}U_{mk} {\mathrm{Str}}Q_mQ_k$ as a sum over contributions characterized by a given number of 'interacting' matrices Q_m [8]. To this end, we denote by $S[Q_m,Q_k]\equiv \frac{2}{t}U_{mk} {\mathrm{Str}}Q_mQ_k$ the sum of two terms which couple Q_m with Q_k. Then the 'interaction' part of the partition function can be expanded as follows

$$\begin{eqnarray} &&\fl {\rm e}^{-\sum _{m>k}S[Q_m,Q_k]}=\prod _{m>k}{\rm e}^{-S[Q_m,Q_k]}=\prod _{m>k}(1+F_{mk})\nonumber \\ &&{\hphantom{{\rm e}^{-\sum _{m>k}S[Q_m,Q_k]}}} \fl =1+\sum _{m>k}F_{mk}+\sum _{m>k>p}(F_{mk}F_{kp}+F_{mp}F_{kp}+F_{mk}F_{mp}+F_{mk}F_{kp}F_{mp})+\cdots , \end{eqnarray} \tag{ 2 }$$

where F_mk stands for ${\rm e}^{-S[Q_m,Q_k]}-1$. The first term of the above expansion is trivial. In the second one the contributions of pairs of supermatrices are collected together, the third one takes into account the contributions of triplets and so on. As a result we obtain the virial expansion of the partition function:

$$\begin{equation} Z=\int DQ\; {\rm e}^{-\sum _{n}S_0[Q_n]}+\sum _{m>k}\int DQ({\rm e}^{-S[Q_m,Q_k]}-1)\, {\rm e}^{-\sum _{n}S_0[Q_n]}+\cdots . \end{equation} \tag{ 3 }$$

In the case of the disordered metallic grains the above expansion is useful, when the grains are weakly coupled.

This general scheme is elaborated in the next two sections, where it is used to calculate the moments of the eigenfunctions I_q and the two-level correlation function R₂(ω) in the unitary symmetry class (β = 2).

The expression for the moments I_q of the eigenfunctions (also known as generalized inverse participation ratios) in terms of the σ-model reads [10, 11]

$$\begin{equation} I_q=\lim _{\eta \rightarrow 0}\frac{\eta ^{q-1}}{4N}\int DQ \sum _n {\mathrm{Str}}^q(Q_n k \Lambda ){\rm e}^{-\frac{1}{t}\sum _{m\ne p}U_{mp} {\mathrm{Str}}Q_mQ_p-\eta \sum _{m} {\mathrm{Str}}\Lambda Q_m}, \end{equation} \tag{ 4 }$$

where $k=\mathbf {1}\otimes \mathrm{diag}(1,-1)$ and Λ = diag(1, −1) in the retarded-advanced representation and N is the total number of sites.

Applying the general scheme of the virial expansion outlined in the previous section, we obtain the following expansion for I_q:

$$\begin{equation} I_q=I_q^{(0)}+I_q^{(1)}+\cdots . \end{equation} \tag{ 5 }$$

The first term in this series corresponds to the first term of equation (3) and is given by

$$\begin{equation} I_q^{(0)}=\frac{1}{N}\sum _n\lim _{\eta \rightarrow 0}\frac{\eta ^{q-1}}{4}\int DQ {\mathrm{Str}}^q(Q_n k \Lambda )\,{\rm e}^{-\eta \sum _{m} {\mathrm{Str}}\Lambda Q_m}. \end{equation} \tag{ 6 }$$

One can notice that the integrand can be factorized and the integrals over different supermatrices Q_m can be calculated separately. Moreover all the integrals over Q_m with m ≠ n are equal to unity due to supersymmetry. The only non-trivial integral over Q_n reproduces the standard expression for the moments of eigenvectors in the Wigner–Dyson RMT [11]. Thus we find that I⁽⁰⁾_q coincides with the standard RMT result

$$\begin{equation} I_q^{(0)}=\Gamma (q+1). \end{equation} \tag{ 7 }$$

Physically it corresponds to the statistics of the eigenfunctions in a system of completely decoupled metallic grains.

The second term of the expansion (5) corresponding to the second term in equation (3) reads

$$\begin{equation} \fl I_q^{(1)}=\sum _n\sum _{m>p}\lim _{\eta \rightarrow 0}\frac{\eta ^{q-1}}{4N}\int DQ \,{\mathrm{Str}}^q(Q_n k \Lambda )\big({\rm e}^{-S[Q_m,Q_p]}-1\big)\,{\rm e}^{-\eta \sum _{m} {\mathrm{Str}}\Lambda Q_m} . \end{equation} \tag{ 8 }$$

Again one can notice that, if m ≠ n and p ≠ n, then the double integral over Q_m and Q_p is determined solely by the anomalous contributions and hence equals to zero. Thus the terms with m = n or p = n are the only ones which should be taken into account:

$$\begin{equation} \fl I_q^{(1)}=\sum _{m,n}\lim _{\eta \rightarrow 0}\frac{\eta ^{q-1}}{4N}\int {\rm d}Q_m \,{\rm d}Q_n {\mathrm{Str}}^q(Q_n k \Lambda )\big({\rm e}^{-S[Q_m,Q_n]}-1\big){\rm e}^{-\eta {\mathrm{Str}}\Lambda (Q_m+Q_n)} . \end{equation} \tag{ 9 }$$

Now we need to calculate the double integral in the above expression. This task is straightforward, but technically quite challenging appreciating a non-trivial parametrization of the supermatrices. Fortunately, exactly the same integral appears at an intermediate step of the derivation of the nonlinear σ-model associated with the ensemble of the banded random matrices [10, 12]. It turns out that the integrals over the Grassmann variables and the phases can be calculated exactly. The next step, which simplifies evaluation of the integral significantly, is to rescale the bosonic variables λ_B in such a way that the limit η → 0 can be taken explicitly in the integrand. Borrowing these results we arrive at the following expression for I⁽¹⁾_q:

$$\begin{eqnarray} &&\fl I_q^{(1)}=\frac{1}{N}\sum _{m,n}\bigg[q(q-1)\int _0^{\infty }{\rm d}z\int _0^{\infty }{\rm d}u\: F_{mn}(z,u)-\Gamma (q+1)\bigg], \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&\fl F_{mn}(z,u)=\sqrt{\frac{\gamma }{2\pi }}z^{q-2}u^{-\frac{3}{2}}{\rm e}^{-z(u+1)-\frac{\gamma }{2}\big(u+\frac{1}{u}\big)}\bigg[\cosh \gamma -\frac{\sinh \gamma }{2\gamma } +\bigg(u+\frac{1}{u}\bigg)\frac{\sinh \gamma }{2}\bigg], \end{eqnarray} \tag{ 11 }$$

where parameter γ = 4U_mn/t is introduced in order to lighten the notation. The appearance of Γ(q + 1) in equation (10) is not accidental and can be traced back to the term containing −1 in equation (9), which gives exactly the same contribution as the one calculated for I⁽⁰⁾_q.

The expansion (5) makes sense only if its different terms scale as different powers of 1/t. Therefore from the exact expression for I⁽¹⁾_q derived above, we need to extract only the first non-trivial term in the limit t → ∞ or γ → 0. At this point our calculations depart from the σ-model derivation, which assumes the opposite limit γ → ∞.

It is obvious from the form of the function F_mn that the limit γ → 0 can not be taken before at least one of the integrations in equation (10) is performed. One way to deal with this problem is to evaluate the integral over u exactly:

$$\begin{equation} \fl \int _0^{\infty }{\rm d}u\: F_{mn}(z,u)= {\rm e}^{-z-\sqrt{\gamma (2z+\gamma )}}z^{q-2}\frac{\gamma (2z+\gamma )\cosh \gamma +(z+\gamma )\sqrt{\gamma (2z+\gamma )}\sinh \gamma }{\gamma (2z+\gamma )}. \end{equation} \tag{ 12 }$$

The expansion of the last expression into the power-series in $\sqrt{\gamma }$ yields

$$\begin{equation} \int _0^{\infty }{\rm d}u\: F_{mn}(z,u)={\rm e}^{-z}z^{q-2}\bigg(1-\sqrt{\frac{\gamma z}{2}}+O(\gamma )\bigg). \end{equation} \tag{ 13 }$$

Substituting this result back into equation (10) and integrating over z we obtain

$$\begin{equation} I_q^{(1)}=-\frac{1}{N}\sum _{m,n}\sqrt{\frac{\gamma }{2}}q(q-1)\Gamma (q-1/2)+O(\gamma ),\quad q>{\textstyle \frac {1}{2}}. \end{equation} \tag{ 14 }$$

The condition q > 1/2 is necessary for convergence of the integral at z → 0. Restoring the dependence of γ on U_mn and taking into account the expression for I⁽⁰⁾_q we find the final result for the moments of the eigenfunctions

$$\begin{equation} \fl I_q=\Gamma (q+1)\bigg(1-\frac{\Gamma (q-1/2)}{\Gamma (q-1)}\frac{1}{N}\sum _{mn}\sqrt{\frac{2 U_{mn}}{t}}+O(t^{-1})\bigg),\quad q>{\textstyle \frac{1}{2}}, \end{equation} \tag{ 15 }$$

which is valid for an arbitrary coupling matrix U_mn.

The appearance of the factor Γ(q + 1) corresponding to the standard RMT result suggests that the eigenfunctions can be written as a product of two contributions. One of them describes the fluctuations of the eigenfunctions inside a grain and it is given by the Wigner–Dyson RMT. Another one, represented by the series in equation (15) describes a change of the amplitudes of the wavefunctions at different grains. The same factorization occurs for the statistics of the eigenfunctions of one-dimensional and quasi-one-dimensional disordered wires and for arbitrary dimensionality in the metallic regime [5].

It follows from equations (7), (10), (13) that the distribution function of the eigenfunction intensities is given by

$$\begin{equation} P(z)=\frac{{\rm d}^2}{{\rm d}z^2}{\rm e}^{-z}\bigg(1-\sqrt{z}\frac{1}{N}\sum _{mn}\sqrt{\frac{2 U_{mn}}{t}}+O(t^{-1})\bigg). \end{equation} \tag{ 16 }$$

We would like to note that the first two terms of the series presented in equation (15) are in one-to-one correspondence with the expansion for the moments of the eigenfunctions found for almost diagonal random matrices characterized by an arbitrary variance matrix 〈|H_mn|²|〉 [13]. It implies, in particular, the same expressions for the fractal dimensions in critical ensembles, as we discuss in section 5.

In this section we focus on the correlation function R₂(ω) measuring the correlations between two levels separated by an energy interval ω

$$\begin{equation} R_2(\omega )=\frac{\langle \nu (E-\omega /2)\:\nu (E+\omega /2)\rangle }{\langle \nu (E)\rangle ^2}, \end{equation} \tag{ 17 }$$

where ν(E) is the density of states and the averaging over a random matrix ensemble or different disorder realizations is assumed.

It is well known that R₂(ω) can be expressed as the following correlator of Q matrices [4, 5]

$$\begin{equation} \fl R_2(\omega )=\bigg(\frac{1}{4N}\bigg)^2\mathrm{Re}\int DQ \bigg[\sum _n {\mathrm{Str}}(Q_n k \Lambda )\bigg]^2{\rm e}^{-\frac{1}{t}\sum _{m\ne p}U_{mp} {\mathrm{Str}}Q_mQ_p+\frac{{\rm i}\pi \nu \omega }{2}\sum _{m} {\mathrm{Str}}\Lambda Q_m}, \end{equation} \tag{ 18 }$$

where it is assumed that ω has an infinitesimal imaginary part to ensure the convergence of the integral. This formula allows us to find the general form of the first two terms of the virial expansion (3):

$$\begin{eqnarray} &&R_2^{(0)}(\omega )=\bigg(\frac{1}{4N}\bigg)^2\mathrm{Re}\int DQ \bigg[\sum _n {\mathrm{Str}}(Q_n k \Lambda )\bigg]^2{\rm e}^{\frac{{\rm i}\pi \nu \omega }{2}\sum _{m} {\mathrm{Str}}\Lambda Q_m}, \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} &&R_2^{(1)}(\omega )=\bigg(\frac{1}{4N}\bigg)^2\mathrm{Re}\int DQ \bigg[\sum _n {\mathrm{Str}}(Q_n k \Lambda )\bigg]^2\sum _{m>k}F_{mk}\,{\rm e}^{\frac{{\rm i}\pi \nu \omega }{2}\sum _{m} {\mathrm{Str}}\Lambda Q_m}. \end{eqnarray} \tag{ 20 }$$

In the following two subsections we evaluate R⁽⁰⁾₂(ω) and R⁽¹⁾₂(ω) respectively.

4.1. Two-level correlation function in the zeroth order approximation

Let us first calculate R⁽⁰⁾₂(ω). To this end we separate the diagonal and the non-diagonal contributions in equation (19)

$$\begin{equation} \bigg[\sum _n {\mathrm{Str}}(Q_n k \Lambda )\bigg]^2=\sum _n {\mathrm{Str}}^2(Q_n k \Lambda )+\sum _{n\ne p} {\mathrm{Str}}(Q_n k \Lambda ) {\mathrm{Str}}(Q_p k \Lambda ) . \end{equation} \tag{ 21 }$$

For the diagonal contribution the integration over any Q_m with m ≠ n gives unity due to supersymmetry. Thus we obtain N equal terms which we denote by I⁽⁰⁾_d. For the same reason, the integration over Q_m with m ≠ n or m ≠ p is trivial for the non-diagonal contribution. In this case we have N² − N equal terms I⁽⁰⁾_nd and hence R⁽⁰⁾₂(ω) can be written as

$$\begin{eqnarray} &&R_2^{(0)}(\omega )=\bigg(\frac{1}{4N}\bigg)^2\mathrm{Re}\bigg[N I^{(0)}_{d}+(N^2-N)I^{(0)}_{nd},\bigg] \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&I^{(0)}_{d}=\int {\rm d}Q \, {\mathrm{Str}}^2(Q k \Lambda )\,{\rm e}^{\frac{{\rm i}\pi \nu \omega }{2} {\mathrm{Str}}\Lambda Q}, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} &&I^{(0)}_{nd}=\bigg(\int {\rm d}Q\, {\mathrm{Str}}(Q k \Lambda )\,{\rm e}^{\frac{{\rm i}\pi \nu \omega }{2} {\mathrm{Str}}\Lambda Q}\bigg)^2. \end{eqnarray} \tag{ 24 }$$

One can note that the integral I⁽⁰⁾_d is exactly the same as the one which determines R₂(ω) in the Wigner–Dyson RMT. It can be easily calculated using the parametrization of Q given in appendix A:

$$\begin{equation} I^{(0)}_{d}=16+\frac{{\rm i}}{\Omega ^2}{\rm e}^{4{\rm i}\Omega }\sin 4\Omega ,\quad \Omega \equiv \frac{\pi \nu \omega }{4} . \end{equation} \tag{ 25 }$$

This result corresponds to the Wigner–Dyson expression for R₂(ω) [4], which is not surprising, as we assume that the statistical properties of a single grain are described by the standard RMT.

The calculation of the non-diagonal contribution is even simpler, as according to equation (A.5) the expression for Str(QkΛ) does not contain the full set of Grassmann variables and hence the integral in equation (24) is anomalous and equals to a constant:

$$\begin{equation} I^{(0)}_{nd}=16. \end{equation} \tag{ 26 }$$

Equation (22) shows that for Ω ≠ 0 it is precisely this contribution that is dominant in the thermodynamic limit N → ∞. Additionally we should take into account the self-correlation contribution contained in equation (25), as it survives in the thermodynamic limit, provided that ωN is kept constant. Thus we find that

$$\begin{equation} R_2^{(0)}(\omega )=1+\frac{1}{\nu N}\delta (\omega ),\quad \mbox{as }\; N\rightarrow \infty . \end{equation} \tag{ 27 }$$

There is a simple way to understand the above result without doing any calculations. Indeed, although the level correlations in each grain are non-trivial and described by the Wigner–Dyson formula, the levels in different grains are assumed to be completely uncorrelated in the zeroth order of the virial expansion. For this reason, when the number of grains N → ∞, all possible correlations except for the self-correlation vanish and 〈ν(E − ω/2) ν(E + ω/2)〉 = 〈ν(E − ω/2)〉〈ν(E + ω/2)〉, which leads to the trivial result (27), provided that the density of states in an energy interval [E − ω/2, E + ω/2] can be approximated by a constant.

4.2. Two-level correlation function in the first order approximation

The first step in calculation of R⁽¹⁾₂(ω) is to use equation (21) to separate the diagonal and the non-diagonal contributions in equation (20). Then one can realize that the diagonal terms give a non-zero result only, if n = m or n = k. All other terms are anomalous and vanish, as F_mk = 0 at Q_m = Q_k = Λ. In a similar way, for the non-diagonal contribution we should keep only the terms with n = m and p = k or n = k and p = m. Thus we need to deal with integrals of two types:

$$\begin{eqnarray} &&R_2^{(1)}(\omega )=\bigg(\frac{1}{4N}\bigg)^2 2\mathrm{Re}\sum _{m>k}\bigg[ I^{(1)}_{d}+I^{(1)}_{nd}\bigg], \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&I^{(1)}_{d}=\int {\rm d}Q\, {\rm d}Q^{\prime } {\mathrm{Str}}^2(Q k \Lambda )\big({\rm e}^{-b {\mathrm{Str}}QQ^{\prime }}-1\big)\, {\rm e}^{\frac{{\rm i}\pi \nu \omega }{2} {\mathrm{Str}}\Lambda (Q+Q^{\prime })}, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&I^{(1)}_{nd}=\int {\rm d}Q\, {\rm d}Q^{\prime } {\mathrm{Str}}(Q k \Lambda ) {\mathrm{Str}}(Q^{\prime }k \Lambda )\big({\rm e}^{-b {\mathrm{Str}}QQ^{\prime }}-1\big) {\rm e}^{\frac{{\rm i}\pi \nu \omega }{2} {\mathrm{Str}}\Lambda (Q+Q^{\prime })} , \end{eqnarray} \tag{ 30 }$$

where the dependence of I⁽¹⁾_d and I⁽¹⁾_nd on m and k is hidden in the parameter b = 2U_mk/t. Calculation of the above integrals is more challenging problem than the one we were faced with in the previous section, as now we are not allowed to consider the limit ω → 0. The main steps of the calculations are outlined below and the details are presented in appendix B.

We start with integration over the Grassmann variables and the phases in a similar way as was done in the previous section. Then we are left with a sum of four-fold integrals over bosonic and fermionic variables λ_B, λ'_B and λ_F, λ'_F. The only term which couples the bosonic with the fermionic variables originates from the measure of integration. It can be represented as an integral over an auxiliary variable of the exponential function. In this way the bosonic and the fermionic integrals can be completely decoupled and considered separately.

The next step is to recall that we are interested in the strong coupling limit t → ∞ or b → 0. One can show that in this limit R⁽¹⁾₂(ω) is given by the asymptotic series

$$\begin{equation} R_2^{(1)}(\omega )=\sum _{k\ge 0}(\sqrt{b})^k f_k\bigg(\frac{\omega }{\sqrt{b}}\bigg). \end{equation} \tag{ 31 }$$

The higher order terms of the virial expansion R⁽ⁱ⁾₂(ω) have similar representations, but with k ⩾ i − 1. Thus only the first term of the series (31) is relevant in the first order approximation. For this reason, we are allowed to take the limit b → 0 keeping $\omega /\sqrt{b}$ constant.

It turns out that the fermionic integrals can be easily calculated in this limit by expanding the corresponding integrands into the power series. All of them give just trivial constant contributions. The bosonic integrals are more tricky due to the fact that the bosonic variables are non-compact. Rescaling the integration variables and keeping only the leading terms of the expansion (31) we arrive at two-fold integrals which can be evaluated exactly and give the following expression for R⁽¹⁾₂(ω):

$$\begin{equation} R_2^{(1)}(\omega )=-\frac{2}{N^2}\sum _{m>k}\bigg[1-\frac{\pi ^{\frac{3}{2}}\nu \omega \sqrt{t}}{\sqrt{8U_{mk}}}\: {\rm e}^{-\frac{(\pi \nu \omega )^2t}{8U_{mk}}}\:\mathrm{erfi}\bigg(\frac{\pi \nu \omega \sqrt{t}}{\sqrt{8U_{mk}}}\bigg) \bigg]. \end{equation} \tag{ 32 }$$

Adding this contribution to the formula for R⁽⁰⁾₂(ω) we obtain the final result for R₂(ω)

$$\begin{equation} R_2(\omega )=\frac{1}{\nu N}\delta (\omega )+\frac{2}{N^2}\sum _{m>k}\frac{\pi ^{\frac{3}{2}}\nu \omega \sqrt{t}}{\sqrt{8U_{mk}}}\: {\rm e}^{-\frac{(\pi \nu \omega )^2t}{8U_{mk}}}\:\mathrm{erfi}\bigg(\frac{\pi \nu \omega \sqrt{t}}{\sqrt{8U_{mk}}}\bigg). \end{equation} \tag{ 33 }$$

This is the main result of this section. Similarly to equation (15) it can be used to compute R₂(ω) for the σ-model with an arbitrary coupling matrix U. Remarkably, the obtained expression coincides exactly with the corresponding result derived for the almost diagonal random matrices [8], if one identifies the coupling matrix U_mk/t with a variance matrix π²ν²〈|H_mk|²|〉/2〈|H_kk|²〉.

In this section we apply the general expressions for I_q and R₂(ω) derived in the previous sections to the σ-model showing critical behavior. One of the most studied examples of this kind is the model characterized by the power-law decay of the coupling matrix elements U_mk∝|m − k|^−α. Such a model was derived in [14] as an effective field theory corresponding to the ensemble of the power-law banded random matrices. It is believed that it exhibits critical properties at all values of the coupling constant provided that α = 2. The mapping between the ensemble of the random matrices and the σ-model can be justified only in the large bandwidth (small coupling) limit. For this reason the critical σ-model was studied in detail in this limit [14–16]. At the same time the model is well defined at arbitrary values of the coupling constant, in particular in the opposite regime of the strong coupling, where it can be explored by the method developed in this work. On the other hand, many analytical results are available for the original random matrix model in the case of the small bandwidth limit [15, 17–21]. This enables us to compare the predictions obtained for the microscopic model and for the σ-model in the situation, when it seems to be unfeasible to establish a direct relation between them.

Specifically, we consider the σ-model (1) with the coupling matrix

$$\begin{equation} \frac{8}{t}U_{mk}=B^2\frac{\pi ^2}{N^2 \sin ^2\big(\frac{\pi (m-k)}{N}\big)}. \end{equation} \tag{ 34 }$$

The periodic form of the coupling matrix (34) is more convenient for analytical calculations.

In order to compute I_q we need to plug the coupling matrix (34) into the general expression (15). Calculating the sum over m and n in the leading order in N we obtain the following expression

$$\begin{equation} I_q=\Gamma (q+1)\bigg(1-B\frac{\Gamma (q-1/2)}{\Gamma (q-1)}\ln N+\cdots \bigg),\quad q> \frac{1}{2}. \end{equation} \tag{ 35 }$$

Comparing it with the equation $I_q\propto N^{-d_q(q-1)}$ defining the fractal dimensions d_q we find the formula for d_q in the lowest order in B:

$$\begin{equation} d_q=B \frac{\Gamma (q-1/2)}{\Gamma (q)}, \end{equation} \tag{ 36 }$$

which is the same as the result that one would obtain for the power-law random model with 〈|H_mm|²〉 = 1 and

$$\begin{equation} \langle |H_{mk}|^2\rangle =\bigg(\frac{B}{\sqrt{2\pi }}\bigg)^2\frac{\pi ^2}{N^2 \sin ^2\big(\frac{\pi (m-k)}{N}\big)}, \quad m\ne k . \end{equation} \tag{ 37 }$$

We conclude that identifying the coupling matrix U_mk/t in the σ-model with the expression π〈|H_mk|²〉/4 in the random matrix model, we obtain a one-to-one correspondence between two results.

Next we calculate R₂(ω) for the same model. To this end we need to evaluate the double sum in equation (32) for U_mk given by equation (34) in the limit N → ∞. In this limit one can first introduce rescaled variables s = Nνω, x = sin (π|m − k|/N), ϕ = πm/N and then replace the sums by the corresponding integrals:

$$\begin{equation} R_2^{(1)}(s)=-1+\frac{2s}{B\pi ^{\frac{3}{2}}}\int _0^1 {\rm d}\phi \int _0^{\sin \phi }{\rm d}x \frac{x}{\sqrt{1-x^2}}{\rm e}^{-\frac{s^2x^2}{b^2}}\mathrm{erfi}\bigg(\frac{sx}{b}\bigg). \end{equation} \tag{ 38 }$$

Integrating by parts over ϕ we are left with an integral over x, which can be evaluated exactly:

$$\begin{equation} R_2^{(1)}(s)=-{\rm e}^{-\frac{s^2}{B^2}}. \end{equation} \tag{ 39 }$$

Finally, taking into account the contribution of R⁽⁰⁾₂(s) we obtain

$$\begin{equation} R_2(s)=\delta (s)+1-{\rm e}^{-\frac{s^2}{B^2}}. \end{equation} \tag{ 40 }$$

Again we note that this expression is in one-to-one correspondence with the random matrix result, provided that we use the same identification rule U_mk/t → π〈|H_mk|²〉/4 as before.

The explicit calculations performed in the previous sections show that the first order results for both I_q and R₂(ω) are equivalent to the corresponding expressions derived for random matrix models using the virial expansion method. The equivalence is very general and holds true for an arbitrary coupling matrix U_mk. In particular, we demonstrated that there is one-to-one correspondence between these two quantities calculated in the ensemble of the power-law banded random matrices at b ≪ 1 and their counterparts calculated in the σ-model. Such a correspondence is rather unexpected taking into account the fact that the standard mapping of the power-law banded random matrix ensemble onto the σ-model can be justified only in the opposite limit b ≫ 1.

A natural question, which appears in this context, is whether the observed equivalence can be extended to higher order terms and if it is the case, then what is the reason behind it? Below we show that under rather general assumptions both the σ-models and the random matrix models can be mapped onto the same effective field theory, so that such an equivalence is preserved at any order of the virial expansion.

In [9], it was noted that one can formulate a supersymmetric field-theoretical description of disordered quantum systems or random matrix models, which is similar to the nonlinear σ-model, but is exact and free from the saddle-point approximation. The form of the action of the new theory is similar to the action of the σ-model (1), but the target space of the new supermatrix $\hat{Q}$ is different from the target space of Q. This approach was subsequently applied to the almost diagonal random matrices in [8]. In particular, it was shown in [8] that in the strong coupling limit the target space of $\hat{Q}$ can be reduced by the constraint ${\mathrm{Str}}\,\hat{Q}=0$, provided that the large scale approximation can be used. The latter assumes that an exact form of the coupling matrix U_km at small distances |k − m| plays no role. More precisely it requires that the dominant contribution to the sum $\sum _{m}\sqrt{U_{km}}$ is determined by large distances |k − m| ≳ X_c, at which U_km is small U_km ≪ 1. In particular this condition is automatically fulfilled in the thermodynamic limit N → ∞ , if $\sum _{m}\sqrt{U_{km}}$ diverges and U_km → 0 in the limit |k − m| → ∞.

The constraint ${\mathrm{Str}}\:\hat{Q}=0$ implies that only two real and four Grassmann variables are required for parametrization of $\hat{Q}$. Its block structure in the retarded-advanced representation is similar to the block structure of Q given by equation (A.1):

$$\begin{equation} \hat{Q}=\left(\begin{array}{@{}cc@{}}U&0\\ 0&V\end{array}\right)\left(\begin{array}{@{}cc@{}}\tilde{D}_{RR}&\tilde{D}_{RA}\\ \tilde{D}_{AR}&\tilde{D}_{AA}\end{array}\right)\left(\begin{array}{@{}cc@{}}U^{-1}&0\\ 0&V^{-1}\end{array}\right), \end{equation} \tag{ 41 }$$

where U and V are exactly the same matrices as those that involved in the Efetov parametrization of Q and defined by equation (A.2) [9]. Each of the matrix $\tilde{D}_{\alpha \beta }$ has only one non-zero matrix element:

$$\begin{eqnarray} &&\tilde{D}_{RR}=\mathrm{diag}(\tilde{\lambda },0),\quad \tilde{D}_{RA}=\mathrm{diag}(-\tilde{\lambda }{\rm e}^{-{\rm i}\phi },0),\nonumber \\ &&\tilde{D}_{AR}=\mathrm{diag}(\tilde{\lambda }{\rm e}^{{\rm i}\phi },0),\quad \tilde{D}_{RR}=\mathrm{diag}(-\tilde{\lambda },0),\nonumber \\ &&\tilde{\lambda }\in [0,\infty ),\quad \phi \in [0,2\pi ). \end{eqnarray} \tag{ 42 }$$

Now let us return to the virial expansion of the σ-model. To be specific we focus on R₂(ω), but our consideration can be applied for any other correlation function. The nth term of the virial expansion involves integration over n supermatrices Q_k and it can be written as the following asymptotic series:

$$\begin{equation} R_2^{(n)}(\omega )=\sum _{k\ge 0}(\sqrt{t})^{-(k+n-1)} f_k^{(n)}(\sqrt{t}\omega ). \end{equation} \tag{ 43 }$$

Using the same arguments as those, which were employed in the discussion of the validity of the approximation ${\mathrm{Str}}\:\hat{Q}=0$ (see [8] for details), one can show that actually all the terms with k ⩾ 1 can be neglected, provided that the large scale approximation is valid. In appendix C we discuss the validity of this approximation for the power-law σ-model.

In order to find the first term of the series, we scale the variables $\tilde{\lambda }=\lambda _B/\sqrt{t}$ and $\tilde{\omega }=\sqrt{t}\omega$ and take the limit t → ∞. In this limit the integrand preserves its form, but the matrices Q_k should be replaced by

$$\begin{equation} \tilde{Q}_k=\lim _{t\rightarrow \infty }\frac{1}{\sqrt{t}}Q_k\big.\big|_{\tilde{\lambda }\; \mathrm{is\; constant}}. \end{equation} \tag{ 44 }$$

Looking at the parametrization of Q ((A.1)–(A.3)) we can immediately notice that $\tilde{Q}$ can be obtained from Q by setting the variables λ_F and μ_F equal to zero and by replacing μ_B by λ_B. Identifying ϕ_B with ϕ from equation (42) we conclude that $\tilde{Q}=\hat{Q}$. The integrands for both models are the same and the measure can be written as

$$\begin{equation} {\rm d}\hat{Q}={\rm d}\tilde{Q}=2\frac{{\rm d}\tilde{\lambda }}{\tilde{\lambda }^2}\frac{{\rm d}\phi }{2\pi } {\rm d}\rho \: {\rm d}\rho ^{\ast }{\rm d} {\rm i}\sigma \: {\rm d} {\rm i}\sigma ^{\ast }. \end{equation} \tag{ 45 }$$

Since the above argument holds true at any order of the virial expansion we arrive at the conclusion that in the strong coupling limit and under the assumption of the large scale approximation both models can be mapped onto the same supersymmetric field theory. The action of this theory is the same as the one of the σ-model and the field variables $\hat{Q}_k$ are parametrized according to equations (41), (42). In particular, it implies that all the results derived for the almost diagonal random matrix models are applicable for the corresponding σ-models in the strong coupling limit and vice versa.

I have benefitted from discussions with Yan Fyodorov, Vladimir Kravtsov and Oleg Yevtushenko. I acknowledge support from the Engineering and Physical Sciences Research Council, grant number EP/G055769/1 and hospitality of the Abdus Salam ICTP.

We employ the standard parametrization of Q introduced by Efetov [3] and adopt the notation used in [12]. In the unitary symmetry class (β = 2) Q has the following block structure in the retarded-advanced representation:

$$\begin{equation} Q=\left(\begin{array}{@{}cc@{}}U&0\\ 0&V\end{array}\right)\left(\begin{array}{@{}cc@{}}D_{RR}&D_{RA}\\ D_{AR}&D_{AA}\end{array}\right)\left(\begin{array}{@{}cc@{}}U^{-1}&0\\ 0&V^{-1}\end{array}\right). \end{equation} \tag{ A.1 }$$

The matrices U and V contain only Grassmann variables ρ, ρ* and σ, σ*:

$$\begin{equation} U=\left(\begin{array}{@{}cc@{}}1+\rho \rho ^{\ast }/2 &\rho \\ \rho ^{\ast }& 1+\rho ^{\ast }\rho /2\end{array}\right)\quad V=\left(\begin{array}{@{}cc@{}}1-\sigma \sigma ^{\ast }/2 &{\rm i}\sigma \\ {\rm i}\sigma ^{\ast }& 1-\sigma ^{\ast }\sigma /2\end{array}\right). \end{equation} \tag{ A.2 }$$

The matrices D_αβ contain other four real variables λ_B, λ_F and ϕ_B, ϕ_F:

$$\begin{eqnarray} D_{RR}&=& -D_{AA}=\left(\begin{array}{@{}cc@{}}\lambda _B& 0\\ 0 & \lambda _F\end{array}\right),\;\; \lambda _B\in [1,\infty ),\: \lambda _B\in [-1,1]\nonumber \\ D_{RA} &=& \left(\begin{array}{@{}cc@{}}-\mu _B{\rm e}^{-{\rm i}\phi _B} & 0\\ 0 & \mu _F{\rm e}^{-{\rm i}\phi _F}\end{array}\right),\: D_{AR}=\left(\begin{array}{@{}cc@{}}\mu _B{\rm e}^{{\rm i}\phi _B} & 0\\ 0 & \mu _F{\rm e}^{{\rm i}\phi _F}\end{array}\right),\nonumber \\ \phi _B,\:\phi _F &\in& [0,2\pi ),\quad \mu _B=\sqrt{\lambda _B^2-1},\:\mu _F=\sqrt{1-\lambda _F^2}. \end{eqnarray} \tag{ A.3 }$$

The measure in this parametrization is given by

$$\begin{equation} {\rm d}Q=\frac{{\rm d}\lambda _B\: {\rm d}\lambda _F}{(\lambda _B-\lambda _F)^2}\frac{{\rm d}\phi _B\: {\rm d}\phi _F}{(2\pi )^2} {\rm d}\rho \: {\rm d}\rho ^{\ast }{\rm d} {\rm i}\sigma \: {\rm d} {\rm i}\sigma ^{\ast }. \end{equation} \tag{ A.4 }$$

Using the expressions for the matrices k = diag(1, −1, 1, −1) and Λ = diag(1, 1, −1, −1) and the above representation of Q, one can easily calculate Str(ΛQ) and Str(kΛQ):

$$\begin{eqnarray} {\mathrm{Str}}(\Lambda Q)&=&2(\lambda _B-\lambda _F),\nonumber \\ {\mathrm{Str}}(k\Lambda Q)&=&2\big[\lambda _B+\lambda _F+(\lambda _B-\lambda _F)(\rho \rho ^{\ast }-\sigma \sigma ^{\ast })\big]. \end{eqnarray} \tag{ A.5 }$$

According to equation (29) the diagonal contribution in the first order of the virial expansion can be written as

$$\begin{eqnarray} &&I^{(1)}_{d}=T(b)-T(0), \end{eqnarray} \tag{ B.1 }$$

$$\begin{eqnarray} &&T(b)=\int {\rm d}Q \, {\mathrm{Str}}^2(Q k \Lambda )\, {\rm e}^{2 {\rm i}\Omega {\mathrm{Str}}\Lambda Q}Y(Q),\quad \Omega \equiv \frac{\pi \nu \omega }{4}, \end{eqnarray} \tag{ B.2 }$$

$$\begin{eqnarray} &&Y(Q)=\int {\rm d}Q^{\prime }{\rm e}^{-b {\mathrm{Str}}QQ^{\prime }+2 {\rm i}\Omega {\mathrm{Str}}\Lambda Q^{\prime }}. \end{eqnarray} \tag{ B.3 }$$

Changing the integration variable Q' by U⁻¹Q'U in the integral defining the function Y(Q), one can show that Y(Q) = Y(λ_B, λ_F) [10, 12]. Using this fact and the expressions for Str ΛQ and Str QkΛ from appendix A, we obtain

$$\begin{eqnarray} &&\fl T=4\int _1^{\infty }{\rm d}\lambda _B\int _{-1}^1{\rm d}\lambda _F\frac{1}{(\lambda _B-\lambda _F)^2}\int _0^{2\pi }\frac{{\rm d}\phi _B {\rm d}\phi _F}{(2\pi )^2} \int {\rm d}\rho ^{\ast }{\rm d}\rho {\rm d}{\rm i}\sigma ^{\ast }{\rm d}{\rm i}\sigma \:{\rm e}^{4{\rm i}\Omega (\lambda _B-\lambda _F)}\nonumber \\ &&{\hphantom{T=}} \fl \times \big[(\lambda _B+\lambda _F)^2+2\big(\lambda _B^2-\lambda _F^2\big)(\rho \rho ^{\ast }-\sigma \sigma ^{\ast })-2(\lambda _B-\lambda _F)^2\rho \rho ^{\ast }\sigma \sigma ^{\ast }\big]Y(\lambda _B,\lambda _F). \end{eqnarray} \tag{ B.4 }$$

The three terms in the brackets in the above expression generate three contributions to T. The first one T₁ corresponding to (λ_B + λ_F)² is anomalous and therefore is given by

$$\begin{equation} T_1=16Y(1,1). \end{equation} \tag{ B.5 }$$

The second one corresponding to 2(λ²_B − λ_F²)(ρρ* − σσ*) is zero, as the bosonic part of the integral converges and the fermionic one is zero. Finally, the third one T₂ generated by the last term in the brackets has the form

$$\begin{equation} T_2=8\int _1^{\infty }{\rm d}\lambda _B\int _{-1}^1{\rm d}\lambda _F\:{\rm e}^{4{\rm i}\Omega (\lambda _B-\lambda _F)}Y(\lambda _B,\lambda _F). \end{equation} \tag{ B.6 }$$

Now we need to use the explicit form of the function Y, which can be found by performing integration in equation (B.3). The result is available in [10, 12]:

$$\begin{eqnarray} &&\fl Y(\lambda _B,\lambda _F)={\rm e}^{-2b(\lambda _B-\lambda _F)}\nonumber \\ &&{\hphantom{Y(\lambda _B,\lambda _F)}} \fl +\int _1^{\infty }{\rm d}\lambda _B^{\prime }\int _{-1}^1{\rm d}\lambda _F^{\prime }\:\frac{\lambda _B-\lambda _F}{\lambda _B^{\prime }-\lambda _F^{\prime }}{\rm e}^{4{\rm i}\Omega (\lambda _B^{\prime }-\lambda _F^{\prime })-2b(\lambda _B\lambda _B^{\prime }-\lambda _F\lambda _F^{\prime })} L(\lambda _B^{\prime },\lambda _F^{\prime },\lambda _B,\lambda _F), \end{eqnarray} \tag{ B.7 }$$

$$\begin{eqnarray} &&\fl L(\lambda _B^{\prime },\lambda _F^{\prime },\lambda _B,\lambda _F)=2b^2\big((\lambda _B^{\prime }\lambda _B+\lambda _F^{\prime }\lambda _F)I_0(2b\mu _B^{\prime }\mu _B)I_0(2b\mu _F^{\prime }\mu _F))\big.\nonumber \\ &&{\hphantom{L(\lambda _B^{\prime },\lambda _F^{\prime },\lambda _B,\lambda _F)=}} \fl \big.-\big[\mu _B\mu _B^{\prime }I_1(2b\mu _B^{\prime }\mu _B)I_0(2b\mu _F^{\prime }\mu _F))-\mu _F\mu _F^{\prime }I_1(2b\mu _F^{\prime }\mu _F)I_0(2b\mu _B^{\prime }\mu _B))\big]\big), \nonumber \\ \end{eqnarray} \tag{ B.8 }$$

where I₀ and I₁ are the Bessel functions. One can see that the above integral is zero at λ_B = λ_F = 1 and hence Y(1, 1) = 1. Thus we find the answer for T₁:

$$\begin{equation} T_1=16. \end{equation} \tag{ B.9 }$$

Another observation which can be made here is that the exponential function and the integral in equation (B.7) generate two different contributions in equation (B.6), which we denote by T_2a and T_2b respectively. The first one can be easily calculated:

$$\begin{eqnarray} &&T_{2a}=8\int _1^{\infty }{\rm d}\lambda _B\int _{-1}^1{\rm d}\lambda _F\:{\rm e}^{4{\rm i}\Omega (\lambda _B-\lambda _F)}{\rm e}^{-2b(\lambda _B-\lambda _F)}=f(\Omega ,b),\nonumber \\ &&f(\Omega ,b)=\frac{4{\rm i}\sin (4\Omega +2b{\rm i})}{(2\Omega +b{\rm i})^2}{\rm e}^{4{\rm i}\Omega -2b}. \end{eqnarray} \tag{ B.10 }$$

This exact result can be represented as a sum of two different contributions: the b independent part and the rest, which can be expanded into the power series:

$$\begin{equation} f(\Omega ,b)=f(\Omega ,0)+\sum _{k\ge 0}(\sqrt{b})^k f_k\bigg(\frac{\Omega }{\sqrt{b}}\bigg). \end{equation} \tag{ B.11 }$$

The explicit form of the b independent part is not important as it will be canceled out later according to equation (B.1). Moreover we are interested only in first term of the power series, as the higher order terms exceed the accuracy of the first order approximation of the virial expansion:

$$\begin{equation} T_{2a}=\frac{{\rm i}\sin (4\Omega )}{\Omega ^2}{\rm e}^{4{\rm i}\Omega }+2\frac{b}{\Omega ^2}+\cdots \; . \end{equation} \tag{ B.12 }$$

Calculation of T_2b requires much more efforts, as T_2b is given by a four-fold integral. However, one can notice that all the terms in equation (B.7) except from the first one are factorized into the fermionic and bosonic contributions. To exploit this observation we introduce the following function

$$\begin{eqnarray} &&\fl F(\Omega ,\Omega ^{\prime })=8\int _1^{\infty }{\rm d}\lambda _B\int _{-1}^1{\rm d}\lambda _F\int _1^{\infty }{\rm d}\lambda _B^{\prime }\int _{-1}^1{\rm d}\lambda _F^{\prime }\:{\rm e}^{4{\rm i}\Omega (\lambda _B-\lambda _F)+4{\rm i}\Omega ^{\prime }(\lambda _B^{\prime }-\lambda _F^{\prime })}\nonumber \\ &&{\hphantom{F(\Omega ,\Omega ^{\prime })=}} \fl \times {\rm e}^{-2b(\lambda _B\lambda _B^{\prime }-\lambda _F\lambda _F^{\prime })}L(\lambda _B^{\prime },\lambda _F^{\prime },\lambda _B,\lambda _F), \end{eqnarray} \tag{ B.13 }$$

which can be used to find T_2b

$$\begin{equation} T_{2b}=\frac{1}{4{\rm i}}\int _0^{\infty }{\rm d}t \frac{\partial }{\partial \Omega }F\big.\bigg(\Omega ,\Omega ^{\prime }+\frac{{\rm i}}{4}t\bigg)\bigg|_{\Omega ^{\prime }=\Omega }. \end{equation} \tag{ B.14 }$$

It has an advantage that now the fermionic and the bosonic integrals do not mix up:

$$\begin{eqnarray} &&\fl F(\Omega ,\Omega ^{\prime })=16 b^2\sum _{\alpha =1}^4\int _1^{\infty }{\rm d}\lambda _B\int _{-1}^1{\rm d}\lambda _F\int _1^{\infty }{\rm d}\lambda _B^{\prime }\int _{-1}^1{\rm d}\lambda _F^{\prime }\:{\rm e}^{4{\rm i}\Omega (\lambda _B-\lambda _F)+4{\rm i}\Omega ^{\prime }(\lambda _B^{\prime }-\lambda _F^{\prime })}\nonumber \\ &&{\hphantom{F(\Omega ,\Omega ^{\prime })=}} \fl \times {\rm e}^{-2b(\lambda _B\lambda _B^{\prime }-\lambda _F\lambda _F^{\prime })}f_{\alpha }(\lambda _B,\lambda _B^{\prime })g_{\alpha }(\lambda _F,\lambda _F^{\prime }), \end{eqnarray} \tag{ B.15 }$$

$$\begin{eqnarray} &&f_1=\lambda _B^{\prime }\lambda _BI_0(2b\mu _B^{\prime }\mu _B),\quad g_1=I_0(2b\mu _F^{\prime }\mu _F),\nonumber \\ &&f_2=I_0(2b\mu _B^{\prime }\mu _B),\quad g_2=\lambda _F^{\prime }\lambda _FI_0(2b\mu _F^{\prime }\mu _F),\nonumber \\ &&f_3=-\mu _B^{\prime }\mu _BI_0(2b\mu _B^{\prime }\mu _B),\quad g_3=I_0(2b\mu _F^{\prime }\mu _F),\nonumber \\ &&f_4=I_0(2b\mu _B^{\prime }\mu _B),\quad g_4=\mu _F^{\prime }\mu _FI_1(2b\mu _F^{\prime }\mu _F). \end{eqnarray} \tag{ B.16 }$$

Thus the problem is reduced now to the calculation of the two-fold bosonic and fermionic integrals, which we denote by A_α and B_α respectively:

$$\begin{equation} F(\Omega ,\Omega ^{\prime })=16 b^2\sum _{\alpha =1}^4 A_{\alpha }B_{\alpha }. \end{equation} \tag{ B.17 }$$

Again we are interested only in the leading term of the expansion (B.11) and therefore we are allowed to consider the limit b → 0 keeping the ratios $\Omega /\sqrt{b}$ and $\Omega ^{\prime }/\sqrt{b}$ constant. In this limit the fermionic integrals can be easily evaluated by expanding the integrands into the power series and the results read

$$\begin{equation} B_1=B_3=4+O(b),\quad B_2=B_4=O(b). \end{equation} \tag{ B.18 }$$

The bosonic integrals are less trivial due to the non-compact character of the bosonic variables. In order to expand the integrals into the power series we introduce new variables u = 2bλ_Bλ'_B, v = λ_B/λ'_B, $\eta =\Omega /\sqrt{b}$ and $\eta ^{\prime }=\Omega ^{\prime }/\sqrt{b}$. Then keeping only the leading order in b term in the integrands we obtain

$$\begin{eqnarray} A_1&=&\frac{1}{8b^2}\int _0^{\infty }{\rm d}u\int _0^{\infty }{\rm d}v \:\frac{u}{v}\: {\rm e}^{4{\rm i}\big(\eta \sqrt{\frac{uv}{2}}+\eta ^{\prime }\sqrt{\frac{u}{2v}}\big)}\:{\rm e}^{-u}I_0(u).\nonumber \\ A_3&=&-\frac{1}{8b^2}\int _0^{\infty }{\rm d}u\int _0^{\infty }{\rm d}v \:\frac{u}{v}\: {\rm e}^{4{\rm i}\big(\eta \sqrt{\frac{uv}{2}}+\eta ^{\prime }\sqrt{\frac{u}{2v}}\big)}\:{\rm e}^{-u}I_1(u). \end{eqnarray} \tag{ B.19 }$$

In a similar way, one can show that A₄ converges while A₂ logarithmically diverges as b → 0. According to equation (B.18) this implies that their contributions to F(Ω, Ω') can be neglected. Substituting the results of equations (B.19), (B.18), (B.17) into equation (B.14) and integrating over t and v we find

$$\begin{equation} T_{2b}=-32\int _0^{\infty }{\rm d}u \:u K_2\bigg(-8{\rm i} \Omega \sqrt{\frac{u}{2b}}\bigg){\rm e}^{-u}\big[I_1(u)-I_0(u)\big], \end{equation} \tag{ B.20 }$$

where K₂ is the modified Bessel function. Collecting the results of equations (B.9), (B.12), (B.14) we obtain

$$\begin{eqnarray} &&\fl T(b)=16+\frac{{\rm i}\sin (4\Omega )}{\Omega ^2}{\rm e}^{4{\rm i}\Omega }+2\frac{b}{\Omega ^2}-32\int _0^{\infty }{\rm d}u \:u K_2\bigg(-8{\rm i} \Omega \sqrt{\frac{u}{2b}}\bigg){\rm e}^{-u}\big[I_1(u)-I_0(u)\big]. \nonumber\\ \end{eqnarray} \tag{ B.21 }$$

In order to compute the non-diagonal contribution equation (30) exactly the same approach can be used. Skipping the details of the calculations, we present here the final result for S(b), which plays the same role as T(b) for the diagonal contribution:

$$\begin{equation} S(b)=16-32\int _0^{\infty }{\rm d}u \:u K_0\bigg(-8{\rm i} \Omega \sqrt{\frac{u}{2b}}\bigg){\rm e}^{-u}\big[I_1(u)-I_0(u)\big]. \end{equation} \tag{ B.22 }$$

It turns out that both integrals in equations (B.21) and (B.22) can be computed analytically and the real part of their sum is given by

$$\begin{eqnarray} &&\fl \mathrm{Re}\int _0^{\infty }{\rm d}u \:u \bigg[K_0\bigg(-8{\rm i} \Omega \sqrt{\frac{u}{2b}}\bigg)+K_2\bigg(-8{\rm i} \Omega \sqrt{\frac{u}{2b}}\bigg)\bigg]{\rm e}^{-u}\big[I_1(u)-I_0(u)\big]\nonumber \\ &&=\frac{b}{16 \Omega ^2}+\frac{1}{2}-\frac{\sqrt{\pi }\Omega }{\sqrt{b}}\:{\rm e}^{-\frac{4\Omega ^2}{b}}\:\mathrm{erfi}\bigg(\frac{2\Omega }{\sqrt{b}}\bigg), \end{eqnarray} \tag{ B.23 }$$

where $\mathrm{erfi}(z)=(2/\sqrt{\pi })\int _0^z{\rm e}^{x^2}{\rm d}x$. Therefore the final expression for the sum of the diagonal and non-diagonal contributions reads

$$\begin{equation} \mathrm{Re}\big(I_d^{(1)}+I_{nd}^{(1)}\big)=-16\bigg[1-\frac{\sqrt{\pi }\Omega }{\sqrt{b}}\:{\rm e}^{-\frac{4\Omega ^2}{b}}\:\mathrm{erfi}\bigg(\frac{2\Omega }{\sqrt{b}}\bigg)\bigg]. \end{equation} \tag{ B.24 }$$

Substituting this expression into equation (28) and recalling the definition of b and Ω we arrive at equation (32).

In this appendix we show, on the example of the power-law σ-model, why the term of the order of 1/t in the expansion of I⁽¹⁾_q can be neglected in comparison with the term of the same order originating from the expansion of I⁽²⁾_q.

It follows from equations (10), (13) that the next term in equation (15) has the following structure:

$$\begin{equation} G(q)\frac{1}{N}\sum _{mn}\frac{U_{mn}}{t}, \end{equation} \tag{ C.1 }$$

where G(q) is some function of q, which does not depend on N. Substituting the expression for U_mn (34) into this formula we obtain that in the thermodynamic limit N → ∞ the result is CG(q)B² with some constant C.

On the other hand the leading order in B term for I⁽²⁾_q will be the same as the one calculated for the almost diagonal random matrices and it is given according to [18] by the following expression

$$\begin{equation} B^2(\tilde{G}_1(q)\ln ^2 N +\tilde{G}_2(q)\ln N+\cdots ), \end{equation} \tag{ C.2 }$$

with some N-independent functions $\tilde{G}_1(q)$ and $\tilde{G}_2(q)$. One can see that although this term has the same order in B it will be dominant in the thermodynamic limit N → ∞. The reason for a different N-dependence appearing here is related to the fact, that in I⁽²⁾_q there is an additional summation over the indices of the coupling matrix $\sum _{mnp}\sqrt{\frac{U_{mn}}{t}}\sqrt{\frac{U_{np}}{t}}$. The importance of such an additional summation for the justification of the large scale approximation was emphasized in [8].