Unification of the Nature’s Complexities via a Matrix Permanent—Critical Phenomena, Fractals, Quantum Computing, ♯P-Complexity

We reveal the analytic relations between a matrix permanent and major nature’s complexities manifested in critical phenomena, fractal structures and chaos, quantum information processes in many-body physics, number-theoretic complexity in mathematics, and ♯P-complete problems in the theory of computational complexity. They follow from a reduction of the Ising model of critical phenomena to the permanent and four integral representations of the permanent based on (i) the fractal Weierstrass-like functions, (ii) polynomials of complex variables, (iii) Laplace integral, and (iv) MacMahon master theorem.


Introduction
We find a remarkable explicit connection between the major types of complexity in nature. They represent the critical phenomena, fractal structures in the theory of chaos, quantum information processing in many-body physics, cryptography, number-theoretic complexity in mathematics, and P-complete problems in the theory of computational complexity. We show that all of them are analytically related to a well-known in mathematics matrix permanent via the fractal Weierstrass-like functions and polynomials or determinants involving complex variables.
The analysis is based on the concept of the P/NP-complexity of computations and quantum information processing and computing (Section 2) as well as on a nontrivial reduction of the critical phenomena problem to a permanent (Section 3) and new integral representations of the permanent revealing its deep explicit relation to the fractals and chaos (Section 4), complex stochastic multivariate polynomials (Section 5), number-theoretical functions (Section 6), asymptotics of a Toeplitz determinant employed in the Onsager's solution of Ising model and given by the Szegő limit theorems (Section 7).
The permanent of a n × n matrix (A pq ) is defined similar to a determinant [1][2][3][4] via a sum running over permutations σ of 1, 2, ..., n, that is, over the symmetric group S n , perA = ∑ (1) for the Ising model of N spins s = 1 2 located in a cubic lattice via a permanent of a circulant matrix A = 2 + g −1 1 . Correlation functions and other characteristics of the critical phenomena also can be expressed in a similar way. The result (2) is based on a bosonization of a many-body constrained system via a Holstein-Primakoff representation of spins. A matrix g 1rr is a correlation function g 1 (r, r ) of an auxiliary, much simpler system of the related unconstrained bosons and can be found by known methods. The submatrix [A] {r} differs from A by an absence of one row and one column which intersect at the entry A rr .
This nontrivial constructive reduction of the Ising model to computing the permanent, that is a P-problem, has been annotated in References [32,33]. It is derived in the present paper. In particular, it implies that a full analysis of the Ising model by numerical simulations alone is intractable. This fact stresses an importance of the exact general representation of the solution via a permanent, which unveils a remarkably canonical analytic structure of statistics and thermodynamics of the critical phenomena and guides to the adequate approximations and asymptotics for their computation including the known mean-field (random-phase) and renormalization-group approximations.
This novel approach can be compared with two other, more formal, plain approaches also pointing to a connection between the critical phenomena and permanent. One of them is a long known combinatorial method of obtaining all expansions and formulas of quantum field theory without diagrammatic expansions [34]. It stems from the basic representation of a many-particle wave function in a system of N Bose particles as a symmetrized product of the single-particle wave functions that is the Caianiello permanent. The latter is similar to the Slater determinant which represents a many-particle wave function in a system of N Fermi particles. The other approach is based on the graph theory. In particular, an ad hoc counting of all matchings of a bipartite graph representing a monomer-dimer model of phase transitions allows one to express its partition function via the permanent of a 0-1 matrix adjacent to the bipartite graph [35]. Importantly, the graph theory and the Markov chain Monte Carlo method provide a fully polynomial randomized approximation scheme (FPRAS) for numerical computation of the permanent of nonnegative matrices and a ferromagnetic Ising model [36][37][38][39]; for a discussion of a different scheme, see Reference [40].
For calculating the exact partition function, magnetization, Green's functions (13), and other nonpolynomial averages, we employ a nonpolynomial diagram technique and partial operator contractions [32,41]. The point is that the constrained, true Green's functions do not obey equations of a Dyson type due to a presence of the nonpolynomial functions θ(2s −n r ), and a standard diagram technique is not suited to deal with them.

The Constrained Spin Bosons in the Holstein-Primakoff Representation
Let us consider a 3d cubic lattice of N interacting quantized spins s = 1 2 with a period a in a box with a volume L 3 and periodic boundary conditions. (The method is valid for an arbitrary dimensionality of the lattice d = 1, 2, 3, ...). The lattice sites are enumerated by a position vector r. According to the Holstein-Primakoff representation [42], worked out also by Schwinger [43], each spin is a system of two spin bosons, which are constrained to have a fixed total occupation n 0r +n r = 2s;n r =â † râr ,n 0r =â † 0râ0r .
Theâ r andâ 0r are the annihilation operators obeying the Bose canonical commutation relations: [â r ,â † r ] = δ r,r , [â 0r ,â † 0r ] = δ r,r , and all (r)-operators commute with all (0r )-operators; δ r,r is a Kronecker delta. A vector spin operatorŜ r at a site r is given by its components: A proper reduction of a many-body Hilbert space ensures [41] that this system is isomorphic to a system of N spin-boson excitations, described by annihilation operatorsβ r at each site r and obeying the Bose canonical commutation relations [β r ,β † r ] = δ r,r , if we cutoff them by a step-function θ(2s −n r ); θ(x) = 1 if x ≥ 0 and θ(x) = 0 if x < 0. This isomorphism is valid on an entire physically allowed Hilbert space and is achieved by equating the annihilation operatorsβ r =β r θ(2s −n r ) of those constrained, true excitations to the cutoff Holstein-Primakoff's transition operators: β r =â † 0r (1 + 2s −n r ) −1/2â r θ(2s −n r ).
Here and thereinafter we add a prime to a symbol of an unconstrained quantity to denote its cutoff, constrained counterpart. The vector components of the spin operator arê where the spin raising and lowering operators are equal tô S + r = √ 2s −n rβ r ,Ŝ − r =β † r √ 2s −n r ;n r =β † rβ r .
The aforementioned isomorphism is not trivial since it is not valid outside the constrained, physically allowed Hilbert space and the commutation relations for the creation and annihilation operators of the true spin excitations in Equations (5) are not canonical, [β r ,β † r ] = δ r,r (1 − (2s + 1)δn r ,2s ). (8) A free Hamiltonian of a system of N spins in the lattice is determined by a Zeeman energy −gµ B B extŜ z of a spin in an external magnetic field B ext (which is assumed homogeneous and directed along the axis z) via a g-factor and a Bohr magneton µ B = eh 2Mc . We intentionally define the free Hamiltonian in Equation (9) via the unconstrained occupation operatorŝ n r =β † rβr on a full Fock space generated by a set of the creation operators {β † r }, that is on the extended many-body Hilbert space without any θ(2s −n r ) cutoff factors. This makes the free Hamiltonian purely quadratic which is necessary for a validity of the standard diagram technique. The latter is crucial for a derivation of the Dyson-type equations, like Equation (15). One is allowed to skip the θ(2s −n r ) cutoff factors in H 0 in virtue of an equalityβ † rβr =β † rβ r , valid on the physical many-body Hilbert space, and a fact that the occupation operatorn r =β † rβr leaves that space invariant. An interaction Hamiltonian of the Ising model [44] Here a coupling between spins is a symmetric function J r,r = J r−r of a vector r − r , connecting spins. For a spin at a site r 0 there are only the coordination number p of the nonzero couplings J r 0 ,r l = 0 with the neighboring spins at sites r l = r 0 + l; l = 1, ..., p. The result in Equation (11) generalizes the Holstein-Primakoff's one [42] by including the nonpolynomial operator θ(2s −n r )-cutoff functions, which add a spin-constraint nonlinear interaction and are crucially important in a critical region.
Since the Holstein-Primakoff's paper of 1940, there were many unsuccessful attempts to convert it into a rigorous and tractable microscopic theory of critical phenomena in magnetic phase transitions. Note that a well-known Dyson's theory of spin waves in a ferromagnet [45] is invalid in the critical region and restricts an analysis to just a well-formed ordered phase. Due to a lack of a proper mathematical apparatus, in particular, a technique of a partial contraction of operators and a diagram technique for the nonpolynomial averages, Dyson thought that "the Holstein-Primakoff formalism is thus essentially nonlinear and unamenable to exact calculations".
A total Hamiltonian H = H 0 + H defines, for any operatorÂ, a Matsubara operatorÃ τ = e τHÂ e −τH evolving in an imaginary time τ ∈ [0, 1 T ] in a Heisenberg representation. A symbol T denotes a temperature. A symbolÃ jτ stands for an operator itselfÃ 1τ =Ã τ at j = 1 and a Matsubara-conjugated operator A 2τ =Ã τ at j = 2.
The unconstrained and true Matsubara Green's functions for spin excitations are defined by T τ -ordering [46]: G Here an unconstrained thermal average over an equilibrium statistical operator ρ = e − H T /Tr{e − H T } of the spin-boson excitations is denoted by the angles as and a true, constrained thermal average is denoted as . . .θ /P s . A partition function P s = θ is equal to a cumulative probability of all occupations of the spin excitations in the unconstrained Fock space to be within physically allowed intervals n r ∈ [0, 2s] for all lattice sites r;θ = ∏ r θ(2s −n r ) is a product of all N cutoff factors. In the Ising model there is no coherence, β rτ = 0, and the unconstrained Green's functions obey the usual Dyson equation with a total irreducible self-energy Σ Here the integral operatorsΣ orǦ (0) , applied to any function f jx of an index j and a four-dimensional coordinate x = {τ, r}, stand for a convolution of that function f jx over the variables j, τ, r with the total irreducible self-energy Σ or the free propagator G (0) , respectively, for example, The total irreducible self-energy is defined by an equation

The Order Parameter and Correlation Functions via the True Probabilities of Spin-Boson Occupations
The magnetization at a lattice site r, being an order parameter of the Ising model, is equal to a true average of the spin z-component in Equation (6). For a spin s = 1 2 , it is that is determined by aθ-cutoff, true probability ρ n r =1 of a spin boson at site r to have one quantum of excitation. j 1 τ 1 r 1 depends on r 1 and r 2 only via r 2 − r 1 . It is a Toeplitz matrix with respect to indexes r 1 , r 2 . The matrix g −1 , which is inverse to the matrix (g I 2 I 1 ) of unconstrained correlations, can be calculated by a technique of Toeplitz matrices, known from the theory of the 2d Ising model [44,[47][48][49][50][51][52]. In Equation (19), a quantity ρ n r 1 =n 1 ,n r 2 =n 2 = δn r 1 ,n 1 δn r 2 ,n 2θ /P s (20) stands for a true probability for two spin bosons at the sites r 1 , r 2 to have n r 1 = n 1 and n r 2 = n 2 quanta of excitations. Namely, Equation (19) involves a true probability for two spin-bosons to have zero quanta of excitations n 1 = n 2 = 0 simultaneously.

The Unconstrained Probabilities of Spin-Boson Occupations via the Unconstrained Correlation Matrix
The next step is finding a joined non-cutoff probabilities of the spin-boson occupations at all N lattice sites Actually, for the reduction of the Ising model to a permanent we need to calculate just a particular joined probability ρ 1{m} ≡ ρ {n r =1:r=r 1 ,...,r m ; n r =0:r =r 1 ,...,r m } = f m , of getting unity occupations n r k = 1 for m spin bosons at a subset of sites {m} = {r k , k = 1, ..., m} and zero occupations for all other N − m spin bosons in the lattice, since the latter probability (22) determines the true joined statistics of the spin-boson occupations (calculated below in Section 3.4) and enters Equations (18) and (19) for the true magnetization and true correlation functions. Also, for the self-energy in Equations (17) and (57), we need a similar joined unconstrained distribution of the spin-boson occupations at a subset of sites {M} = {r k , k = 1, ..., M}, which admits arbitrary occupations n r = 0, 1, 2, ..., ∞, r = r 1 , ..., r M , at all other N − M lattice sites, that is, it is non-cutoff averaged over the latter occupations. Again, we need just a particular joined probability the derivatives of which yield those joined probability distributions: We find the characteristic function Θ N by means of the partial operator contraction within the nonpolynomial diagram technique [32,33] as follows Here the diagonal matrices Z Θ and z are related as Z Θ = z/(z − 1). The obtained solution in Equation (31) is normalized to unity at a point {u r = 0}, Θ N ({u r = 0}) = 1, as it should be for a characteristic function of any distribution.
The probability of the unity occupations for m spin bosons and zero occupations for all other spin bosons, Equation (22), is set by a differentiation of that characteristic function: It is a coefficient in front of the multilinear term z r 1 ...z r m in a Taylor expansion of the characteristic function Θ N over the variables {z r } at the zero point {z r = 0}.
In order to evaluate the Taylor expansion, we employ a well-known MacMahon master theorem [53,54]. It yields a Taylor expansion of a function, inversely proportional to a determinant of a matrix 1 − Ax, over the variables {x i }, where s i ≥ 0 is a non-negative integer (i = 1, ..., N), x = diag{x 1 , ..., x N } a diagonal matrix, A a N × N matrix, and per (s 1 ,...,s N ) A a generalized permanent of the matrix A. For the required by Equation (32) multilinear terms with a subset of unity integers {s i = 1; i = 1, ..., m} and the rest of integers being equal zero, the corresponding permanents are reduced to the standard permanent, of the corresponding (m × m)-submatrix A {m} . In order to get the derivatives in Equation (32), one may compute a multilinear expansion of the characteristic function in Equation (31) by taking into account (i) a square-root function, additional to the MacMahon master Equation (33), via the corresponding Bell polynomials of that Faá di Bruno's formula and (ii) an equality of the variables z 1r 1r = z 2r 2r = z r in the two adjacent columns with the same site-index r.
For simplicity's sake, we consider below a case of vanishing anomalous correlations, g 2r 1r = 0 for ∀r, r , and non-zero normal correlations g 1r 1r = g 2r 2r , which are real-valued in the homogeneous phases, when the matrix g j r jr is a circulant Toeplitz matrix in indexes r, r and the arbitrary phases of spin-bosons' annihilation operatorsβ r are calibrated properly (for the calculations in a general case, see Reference [33]). In this case, one has det[g 1r ]] 2 and Equations (31)- (34) yield Here the elements of the (N × N)-matrices g 1 and z 1 as well as (m × m)-matrix (1 + g −1 1 ) {m} are labeled solely by the site-indexes r, r . The effect of the normal cross-correlations g 1r 1r = 0 between the spin bosons at different sites on their joined unconstrained statistics, described by Equation (35), remains highly nontrivial even in that case of vanishing anomalous correlations g 2r 1r = 0 for ∀r, r . For the joined unconstrained, non-cutoff distribution of the spin-boson occupations at only a subset of lattice sites {M} = {r k , k = 1, ..., M}, M ≤ N, defined in Equation (24), a derivation is similar. We just need to restrict the (2N × 2N)-matrices g, Z Θ , z to the corresponding quasi-diagonal (2M × 2M)-block The result for the characteristic function is similar to Equation (31). Obviously, its differentiation, yields the corresponding, similar to Equation ( Finally, we present the explicit formulas for the characteristic functions of the joined non-cutoff probability distributions in the most important cases of the single-site (M = 1) and two-sites (M = 2) subsets of the spin bosons:

The Partition Function and the True Probabilities of Spin-Boson Occupations
Here we give the exact analytic formulas for the partition function and true joined probability distributions ρ n r , ρ n r 1 ,n r 2 , ρ {n r } of the physically allowable spin-boson occupations n r = 0, 1.
Those distributions are simply theθ-cutoff versions of the unconstrained distributions ρ n r , ρ n r 1 ,n r 2 , ρ {n r } , restricted to the unity-occupation ones ρ 1{m} in Equation (32). Note that the unconstrained occupation distributions, calculated in the Section 3.3, already contain all effects of the constraints and spin interaction, except theθ-cutoff only, since they were calculated for the exact, constrained andθ-cutoff, Hamiltonian (11).
We start the analysis of the true joined distribution of the occupations {n r = 0 or 1} for all N spin bosons, with an evaluation of the partition function P s from Equation (13). It is equal to the sum of the probabilities ρ 1{m} in Equation (32) over all occupation configurations {n r = 0 or 1; r = r 1 , ..., r N }, which can be written as follows The second equality in the equation for P s is due to the fact that the terms with the square, z 2 r , and higher powers of any variable z r do not contribute to the considered derivative at the zero point {z r = 0}. Note that the newly introduced function Θ differs from the characteristic function Θ N in Equation (31) only by a substitution of the matrix A = 1 + g −1 with the matrix Thus, an evaluation of the partition function P s can be done similar to the evaluation of the probability ρ 1{m} at m = N described above. In particular, in the case of the vanishing anomalous correlations g 2r 1r = 0 and non-zero normal correlations g 1r 1r = g 2r 2r , as in Equation (35), we find A result for the single-site zero occupation probability differs from P s only by an absence of one partial derivative ∂/∂z r 1 and by a normalization factor. In particular, in the case of the vanishing anomalous correlations g 2r 1r = 0 and non-zero normal correlations g 1r 1r = g 2r 2r , we have The true single-site unity occupation probability is equal The true two-sites zero occupation probability differs from the single-site one in Equation (46) only by an absence of one more partial derivative ∂/∂z r 2 . So, in the case of the vanishing anomalous correlations and non-zero normal correlations g 1r 1r = g 2r 2r , as in Equation (35), one has ρ n r 1 =0,n r 2 =0 = The true two-sites probabilities for other occupation combinations can be computed from the probabilities, presented above, as follows ρ n r 1 =0,n r 2 =1 = ρ n r 1 =0 − ρ n r 1 =0,n r 2 =0 , ρ n r 1 =1,n r 2 =1 = ρ n r 1 =1 − ρ n r 1 =1,n r 2 =0 .
These equations stem from a fact that the true single-site occupation distribution is equal to the true two-sites occupation distribution, averaged over the physically allowable occupations n r 2 = 0, 1 of a spin boson at the second site: ρ n r 1 = ρ n r 1 ,n r 2 =0 + ρ n r 1 ,n r 2 =1 .
The true three-sites and other multiple-sides occupation distributions are not required for calculating the true order parameter and correlation functions, but are necessary for the analysis of the true multiple-sides correlations and statistics. Those m-sides occupation distributions can be computed in a similar way by an induction: ρ n r 1 ,...,n r m−1 ,n rm =1 = ρ n r 1 ,...,n r m−1 − ρ n r 1 ,...,n r m−1 ,n rm =0 .
A detailed analysis of the true spin-boson occupation probability distributions obtained above as well as the true order parameter and correlation functions in Equations (18) and (19) will be given elsewhere, since they are not required for finding the self-energy in Section 3.5 and the exact self-consistency equation in Section 3.6 below.

The Exact Solution for the Total Irreducible Self-Energy via the Unconstrained Correlation Matrix
In order to find the aforementioned spin-boson occupation probabilities and the correlation matrix g J J , it is crucial to get an exact solution to Equation (17) for the total irreducible self-energy which allows one to go beyond standard second-order or ladder approximations. For a given site r 0 in a lattice and any its nearest-neighbor site r l , l = 1, ..., p, we introduce a correlation (4 × 4)-matrix Here q = q † is hermitian, R and R run over two values {r 0 , r l }, A means anti-normal ordering, (2 × 2)-matrices g j j (l) = g j r l jr 0 of basis auto-and cross-correlations are denoted as g(0) = S = S † and g(l = 0) = C(l), respectively. An exact solution for the self-energy is a matrix 2(p + 1)-banded in indexes I 0 = {j 0 , r 0 } and I = {j, r}, where the (2 × 2)-matrix blocks Σ(l) = (Σ jr l j 0 r 0 (l)) are Here ρ n r 0 = δñ r 0 ,n r 0 and ρ n r 0 ,n r l = δñ r 0 ,n r 0 δñ r l ,n r l are the non-cutoff probabilities for the spin bosons at the sites r 0 and r l to acquire n r 0 and n r l quanta of excitations. The probabilities of the zero occupations follow from Equation (39): For the single-site probability one has M = 1 and (2 × 2)-matrix g {1} = S, Equation (56), so that Equations (38) and (40) yield For the two-sites probabilities one has M = 2 and (4 × 4)-matrix g {2} = q (where hereinafter a matrix q j R jR is defined similar to Equation (56) with the R and R running over arbitrary two sites {r 1 , r 2 }, not necessarily neighboring sites), and Equations (38) and (41) yield where [q] I I stands for a I I -submatrix of q, that is, for the matrix q with the I-th row and I -th column deleted.

The Exact Closed Self-Consistency Equation for the Unconstrained Correlation Matrix
Now we can make a final, crucial step in the exact reduction of the Ising model to the matrix permanent-find an exact closed self-consistency equation for the nearest-neighbors', basis normal and anomalous auto-and cross-correlations g (56). Indeed, the total irreducible self-energy in Equation (57) and the spin-boson unconstrained occupation probabilities, Equations (60)-(64), entering formulas for the self-energy, are known exactly via the (1 + p) basis correlation (2 × 2)-matrices g(l), l = 0, 1, ..., p, Equation (56), that is, the matrix S ≡ g(0) = (g j r 0 jr 0 ) of the auto-correlations for a spin boson at the site r 0 and the coordination number p matrices C(l) ≡ g(l = 0) = (g j r l jr 0 ) of the cross-correlations of a spin boson at the site r 0 with the nearest-neighbors at the sites r l = r 0 + l. Due to the complex-conjugation relations there are only two independent, normal g 1r l 1r 0 and anomalous g 2r l 1r 0 , correlation parameters per each basis correlation (2 × 2)-matrix, that is, only 2(1 + p) numbers, which determine all details of the critical phenomena.
Thus, we can find the self-consistency equation for those 2(1 + p) basis auto-and cross-correlations in two steps. First, we solve the Dyson-type Equation (15) for the unconstrained Green's functions in terms of those basis correlations. Second, we close the loop by expressing the basis correlations themselves via those Green's functions.
For the considered stationary homogeneous phases, the Green's functions, the equal-time correlation functions, and the self-energy depend only on the differences of their arguments τ = τ 1 − τ 2 and r = r 2 − r 1 , that is, G Hence, it is straightforward to solve the Dyson-type Equation (15) by means of the Fourier transformation over the imaginary time τ ∈ [− 1 T , 1 T ] and the discrete Fourier transformation over the space. The latter has a following form where the sums run over all lattice sites r with a period a and discrete wave vectors k = {k i |i = 1, . . . , d}, k i = 2π L q with an integer q; k i ∈ [− π a , π a ]. We discern the Fourier transform and its inverse by the arguments k and r. A result for the normal and anomalous Green's functions is where the two quasiparticle eigen-energies depend on the wave vector k via the self-energies The latter Fourier transform of the self-energy consists of only 1 + p terms within a neighborhood of the nearest sites for which there are nonzero couplings J r 0 ,r l = 0. This is a consequence of the fact that the self-energy matrix is a 2(p + 1)-banded matrix. The (2 × 2)-matrix blocks Σ(l) are given explicitly in Equations (57)- (59) representing the exact solution to Equation (17).
The spatial Fourier transforms of the normal and anomalous equal-time correlation functions follow from Equations (68) and (69) in the limit τ → +0: Thus, we derive the equations for the values of the normal and anomalous correlation functions at (1 + p) difference position vectors l = r l − r 0 of the neighboring spins: Their right hand side is determined by the left hand side g 1j (l) itself via Equations (57) and (70)-(73). They constitute an exact closed system of 2(1 + p) self-consistency equations. Its finding is related to a solution of the Ising problem in the same way as finding of a self-consistency equation in the mean-field theory is related to a solution of a phase transition problem. However, now the self-consistency Equation (74) is an exact equation valid in the entire critical region, not just its mean-field approximation. One can analyze these explicit exact self-consistency equations by the well-known in the mean-field theory analytic and numerical tools. It is relatively simple for the Ising model with the zero off-diagonal self-energy Σ 12 = 0 and zero anomalous correlations, when only the (1 + p) self-consistency equations remain. Moreover, in the isotropic case, when the cross-correlations with all p nearest neighbors are the same, the system is reduced to just two equations. Note that the close exact self-consistency equations exist only for the unconstrained, auxiliary basis normal and anomalous auto-and cross-correlations. When the latter are found, the actual, observable statistical and thermodynamic quantities can be explicitly expressed in terms of those basis correlations via the matrix permanent, as is shown in the Sections 3.2-3.4 for the true, constrained partition function, order parameter, correlation functions, and joined statistics of the spin-boson occupations (in particular, see Equations (31)- (35), (45) and (55)).
This completes the exact general reduction of the Ising model to computing the matrix permanent and provides a basis for the calculation of all statistical and thermodynamic characteristics of the critical phenomena via the permanent of the sub-matrices composed from a relatively simple, unconstrained spin-boson correlation matrix g.
Next, we proceed with a discussion of other remarkable features and representations of the matrix permanent ad rem to the analysis of the nature's complexities.

The Permanent and the Fractals
Here and in Section 5 we present a remarkable finding of a direct relation between the permanent and the fractals and chaos. It is based on the two new integral representations of the n × n matrix permanent via (i) an analog of the famous Weierstrass function which is known for its fractal structure and nontrivial Hausdorff dimension and (ii) a mean value of a random multivariate polynomial.
We show that the fractals and chaos are intrinsic to evaluation of the matrix permanent.

The 1d Integral Representation of the Permanent: A Fractal Integrand and a Weierstrass Function
A move is to find an integral representation of the permanent for any n × n matrix A pq in a form of a 1d integral. The idea is as follows. Let us form a sum of quasi-harmonics c k e iπtν k at frequencies ν k with amplitudes c k determined by the matrix entries A pq in such a way that one spectral component of a known frequency ν 0 would have the amplitude equal to the permanent, c 0 = perA. Then, employ an appropriate Fourier integral to discriminate this component and find the permanent as its amplitude. So, we introduce a permanental function P A (z) as a product of the row functions B p (z), The amplitude of the spectral component at the frequency gives the permanent, c 0 = perA, for any base b > 1 if there are no coincidental resonances: ∑ n q=1 n q b q−1 = ν 0 for any partition {n q ≥ 0|q = 1, ..., n} of n = ∑ q n q except the unity one, n q = 1 ∀q. Under this condition, we find the permanent's 1d integral representation as follows A (e iπt )dt for an integer base b = 2, 3, ..., Here an integrand is a following function of a complex variable z = e iπt , For an integer base b, this function is a polynomial in z.
is exponentially broad, with big Hadamard gaps. It was designed so on purpose to make the 1d representation (76) possible. At b ≥ 2, Equation (76) remains valid even if the row function B p is extended from a finite sum of n terms to an infinite seriesB p (z) = ∑ ∞ q=1 A pq z b q−1 by adding the higher z-powers with q > n and any, unrelated to A, factors. The seriesB p is a Weierstrass-type function, like a complex extensioñ of the Weierstrass cosine function W a,b (t) = ReW a,b (e iπt ). The Weierstrass functions are famous for being continuous everywhere but differentiable nowhere. The extensionW a,b (z) is a lacunary (cf. Hadamard gaps) complex power series. A fractal (box or, as is believed, graph is determined by the Hurst, or scaling, exponent α = log b a. The dimension is greater than unity, D > 1, if a ∈ (1, b); see reviews on a fractal geometry of the Weierstrass functions [55,56]. It is known that an image of the unit circle, |z| = 1, under the complex-valued Weierstrass mapW a,b (z = e iθ ), θ ∈ [−π, π], of an integer base b covers an open subset of the complex plane, that is, forms a Peano curve, if the amplitude a is close to 1, that is if the fractal dimension, Thus, the result in Equation (76) reveals a fractal nature of the permanent discussed in detail below. For the critical phenomena, in particular, for the Ising model, the matrix A pq is determined by the correlation function of the unconstrained bosons and evolves from a fast exponential decay, for example, A pq ∼ b −αq , in the disordered phase to a slower than exponential, for example, power-law A pq ∼ q −η , η > 0, decay in the ordered phase. In terms of the permanental row functionsB p (z), such an evolution of the many-body system across the critical region means a transition from (a) the Weierstrass-type function of a large exponent α ≈ 1 and a trivial, almost non-fractal structure with the dimension D ≈ 1 through (b) a sequence of the Weierstrass-type functions of a smaller exponent α 1 and a nontrivial fractal structure with the dimension D = 2 − α larger than unity towards (c) the Weierstrass-type functions of an effectively zero exponent α ≈ 0 and a fully developed fractal structure with the maximal dimension D ≈ 2.

A Fractal Nature of the Matrix Permanent
Here we demonstrate phenomenal fractal properties of the permanental functionP A (z) forming the 1d integral representation in Equation (76). They manifest themselves already in the cases of very simple n × n matrices A pq ≡ 1, A pq = a −q , and A pq = 1 + aδ p,q who's permanents, perA, are n!, n!/a n(n+1)/2 , and e a Γ(n + 1, a) (see Equation (142) below), respectively. In the first two cases the permanental function can be replaced by the n-th power of the Weierstrass function in Equation (78) In order to illustrate a fractal nature of the permanent, let us consider the fractal properties of the integrand and the integral's accumulation in the permanent's representation (76) with the integration range increasing from zero to the ultimate value T = 1 or T = ∞ for the integer or non-integer base b, respectively: First, we elaborate on the basic case of the n × n matrix with unity entries, A pq ≡ 1, for which the asymptotics of the permanent is given by the Stirling's formula, We find that there are two qualitatively different patterns by which the integral (79) approaches the permanent's value: I A (T) → perA at T → 1 or ∞. Which of the two patterns is realized depends on whether the base b of the permanental functionP A is less or larger than the base of the exponential factor in the denominator of the permanent's asymptotics, which is e = 2.718... in the case of Equation (80). (A factor n n could be eliminated by re-scaling the matrix A, say, to a doubly stochastic one.)

Permanent's Fractal: The Case of the Integer Base
The first pattern is illustrated in Figure 1 for the binary base b = 2 < e. In this scenario, the integral in Equation (79) quickly reaches an exceedingly large maximum value ∼n n /b n perA at an exponentially small displacement T ∼ 1/b n 1 from zero and then gradually, with some oscillations, decreases by a sequence of fractal, self-similar steps to the actual value of the permanent as the integration range T tends to unity. (Hereafter, for simplicity's sake, we skip a logarithmic factor ∼log b n in all of the order-of-magnitude estimates.) The second pattern is illustrated in Figure 2 for the ternary base b = 3 > e. In this scenario, the integral in Equation (79) gradually grows from zero to the permanent's value perA all the way from T = 0 to T = 1. This pattern has a fractal, self-similar structure similar to a famous Cantor-Lebesgue function, or Devil's staircase [56]. These patterns could be unambiguously understood by taking into account a remarkable fine structure of the permanental functionP A = [W a,b (e iπt )] n which is neither a smooth function nor a structureless noise, but a fractal hierarchy of the ultrashort pulses/peaks of a width ∆t ∼ b −n and an amplitude regularly scaling from one hierarchy level to the next one. Its fractal structure stems from a more elementary fractal structure of each Weierstrass-function factor. A real part of the function involving a truncated version of the Weierstrass function relevant to a matrix of a finite size n, is shown in Figure 3. A primary series of extrema for the permanental functionP A consists of the pulses located near the points t k 1 = 2 −k 1 , k 1 = 0, 1, ..., n − 1. At the large matrix size n → ∞, all of them, except a few (logarithmic number of) pulses located close to the boundary values of the index k 1 (say, k 1 = 0, 1, 2, 3 and k 1 = n − 4, n − 3, n − 2, n − 1), have equal amplitudes and a universal profile F(∆t), localized within a narrow deviation ∆t ∼ 2 −n from the points t k 1 and described by a special function The primary series of extrema is preceded by an exceptionally large, main peak located at the origin t = 0 and holding the same universal profile of the width ∼2 −n , A secondary series of extrema for the permanental functionP A consists of two sequences of pulses located to the right (s 2 = +1) and to the left (s 2 = −1) from each primary extremum t k 1 at the points Namely these two sequences of secondary peaks are shown in Figure 3 in a vicinity of the primary extremum at the point t 1 = 1 2 in the case of a = 1, b = 2, n = 20. This fractal, self-similar hierarchy of the enclosed into each other extrema's series continues with the ternary and higher, r-order series of pulses surrounding each (r − 1)-order extremum at the points t (s 2 ,...,s r ) k 1 ,k 2 ,...,k r = At the large matrix size n → ∞, all of the r-order pulses, except a logarithmic number of pulses located close to the boundary values 1 and n − 1 − k 1 − ∑ r−1 j=2 k j of the index k r , have equal amplitudes and the universal profile:P A (e iπ(t (s 2 ,...,sr ) k 1 ,k 2 ,...,kr +∆t) ) ≈ ζ r n n F(∆t).
The extrema's amplitudes in this hierarchy of series exponentially decrease with increasing order r of the series as the r-th power of the scaling factor defined in Equation (82), At the same time, the number of pulses in the series grows with increasing order r of the series roughly as 2 r−1 C r n−2 , where a binomial coefficient C r n−2 accounts for a number of r-compositions of the integer n − 1 = ∑ r j=1 k j into a sum of r integers k j ≥ 1 and a factor 2 r−1 accounts for a presence of two, right and left, branches of extrema (s j = ±1) for each j-series in the hierarchy.
For a finite matrix size n, each pulse contributes to the integral in Equation (76), representing the matrix permanent perA, with the universal factor on the order of the pulse width ∼2 −n . Thus, the convergence of the integral in Equation (79) to the exact permanent's value, I A (T) → perA, with the increasing range of integration T → 1, shown in Figure 1, is a subtle interplay between the effects of decreasing amplitude (∼ζ r ) and increasing number (∼C r n−2 ) of pulses in the extrema's r-series as well as a phase modulation of their contribution. The I A (T) accumulates only the real part of the complex-valued permanental functionP A (e iπt ) which contains a phase shift varying with the increasing order r, that is evidenced already by the fact that the factor ζ r in Equation (86) holds the complex number (87). This interplay results in the exponentially small prefactor ∼ √ 2πn(2 n /e n ) 1 in front of the product of the pulse width, 2 −n , and normalization, n n , factors in the value of the permanental integral in Equation (79) at T = 1, I A (T = 1) = perA, as is required by the Stirling's formula in Equation (80). This observation explains why for the binary base b = 2 < e the accumulation pattern of the permanental integral in Equation (79), shown in Figure 1, involves a huge initial growth of the integral to a value ∼n n /2 n due to the main peak (84) and, then, its subsequent almost complete cancellation and fine tuning of the interference contributions from many extrema's r-series of different orders r which finally (at T = 1) lead to the actual, exponentially smaller by the factor ∼ √ 2πn(2 n /e n ) 1, value of the permanent. A similar analysis could be done for the ternary, b = 3 > e, and other integer bases larger than e. In this case, the width of the pulses in the hierarchy of the extrema's series is on the order of b −n that is much smaller than the exponential factor e −n required by the Stirling's formula (80). It yields a value that is much less, by an exponential factor ∼e n /(b n √ 2πn) 1, than the actual value of the permanent and calls for an accumulation of the contributions from many extrema's r-series. This observation explains the fractal pattern of the Cantor-Lebesgue, or Devil's staircase, type in Figure 2.
These results suggest that, starting from a wide enough range of the integer bases [2, b max ] and dividing it in halves according to the observed patterns of the permanental integral accumulation shown in Figure 1 or Figure 2, one could find, in a logarithmic number of steps, an approximation of the exponential factor n n /b n A in the permanent's asymptotics by bounding a true asymptotics' base b A between the two neighboring integer bases, b < b A < b + 1, as it was demonstrated above for the case of the Stirling's asymptotics in Equation (80)

Permanent's Fractal: The Case of the Non-Integer Base
The analysis presented above could be extended to the non-integer bases by switching to the permanent's integral representation with the non-integer base b > 1 in Equation (76). It constitutes an alternative and, probably, more efficient way of computing the true base b A and the pre-exponential factor in the permanent's asymptotics, similar to the factor √ 2πn in the Stirling's asymptotics (80). For instance, let us illustrate how the permanent's asymptotics in Equation (80) arises from the integral representation with the natural-logarithm base b = e in Equation (76).
With increasing range of the integration T, it steadily converges to the matrix permanent's value, as is shown in Figure 4. The analysis of the permanental function with the non-integer base is similar to the binary and ternary ones. Yet, it requires the infinite-limits Fourier integral in Equation (76), instead of the finite-limits integral, for the evaluation of the permanent via the integral spectral discrimination since the permanental function (81) is not a periodical function anymore.
For the case of the n × n matrix with unity entries, A pq ≡ 1, the permanental function is given by Equation (81) with the parameters a = 1 and b = e as follows The hierarchy of its extrema is illustrated in Figure 5 and can be understood in terms of a harmonics' synchronization as follows. Let us consider a differential counter-clockwise rotation of n unity-length links (q = 1, ..., n) in the chain representing the sum,W  The m = 0 extremum is located at the point t 0 = 0 where all links (harmonics) have the same orientation in the east direction: exp(iπt 0 e q−1 ) = 1; q = 1, ..., n.
The next, m = 1, extremum sits at the point t 1 where the last, q = n, link (harmonic) makes a full 2π rotation plus a little extra rotation at an angle equal to a total rotation of all other links (harmonics) q = 1, ..., n − 1. The latter condition can be approximately written as a transcendental equation for the position t 1 as follows For the case of n = 10 in Figure 5, Equation (89) gives the value t 1 ≈ 0.000257 which is very close to the exact position t 1 = 0.000271 of this extremum of the row-sum polynomial in Equation (75). At the point t = t 1 , the next to the last link (harmonic) q = n − 1 makes 90 o + 55 o degrees, that is about 2π/e, of a counter-clockwise rotation and is directed mainly to the west, while all other links (harmonics) are still directed mainly to the east at moderate to small angles above the positive real-valued axis on the complex plane: 35 o for q = n, 53 o for q = n − 2, 20 o for q = n − 3, 7 o for q = n − 4, etc.
The other extrema in the primary series located at the positions t m < 1, m = 2, ..., n − 1, can be viewed and found similarly. Their positions can be approximated as All these extrema at t = t m , m = 1, ..., n − 1, correspond to the configurations with just one link (harmonic of the row-sum polynomial), namely the one labeled by the index q = n − m, directed mainly to the west and all of the other links (harmonics) directed mainly to the east.
The next extremum, m = n, the first one located at the position, t n ≈ t 1 e n−1 > 1, further away from the origin than unity, has all links (harmonics) directed mainly to the east. The following, m = n + 1, extremum in the series (90) has only one, namely the last one with q = n, link (harmonic) directed mainly to the west, while all other links (harmonics) are east directed. The next three, m = n + 2, n + 3, n + 4, extrema in the series (90) have only two (q = 2n − m + 1 and q = 2n − m + 2) links (harmonics) oriented mainly in the west direction.
This pattern continues in the series (90) as following. The next three, m = n + 5, n + 6, n + 7, extrema have already three (q = 2n − m + 1, 2n − m + 2, and q = 2n − m + 5) links (harmonics) oriented mainly in the west direction. The next two, m = n + 8, n + 9, extrema have already four (q = 2n − m + 1, 2n − m + 2, 2n − m + 5, 2n − m + 7) links (harmonics) oriented mainly in the west direction. The next two, m = n + 10, n + 11, extrema have already five links (harmonics) oriented mainly in the west direction. Then the pattern becomes more complicated but still can be followed via the picture of differential rotation of links in the chain of harmonics. Of course, a similar picture is valid for any n, not only for n = 10.
Obviously, there are many higher-order series of extrema in the fractal hierarchy of extrema of the permanental function (88). All of them correspond to the local extrema of the row-sum polynomial in Equation (75), that is, the functionW (n−1) 1,e (t). In fact, at large n, there exist exponentially many closely located local extrema of the row polynomial that makes extremely difficult to numerically differentiate a particular m-order extremum from other series' extrema. Such a Weierstrass fractal structure is globally homogeneous, ergodic along the t axis since this property is required for the asymptotically linear growth of the permanental integral in Equation (76) with an increasing range limit T → ∞.
We find that, with increasing range of the variable t ∈ [−T, T], T → ∞, the path of the scaled first-row permanental function for the n × n matrix A pq = 1, introduced in Equation (75) and additionally scaled by its maximum value n, fully covers a finite 2d region in the complex plane enclosed by a hypocycloid, A number of its cusps is equal to the size n of the matrix. Remarkably, in the limit n → ∞, the hypocycloid perimeter, P n , remains longer than the unit circle perimeter, 2π, by a finite amount, although its area, S n , tends to the area of the unit circle, So, the related scaled permanental function entering the permanent's integral representation in Equation (76),P A /n n = [B 1 (e iπt )/n] n , fully covers a finite 2d region in the complex plane enclosed by the n-th power of the hypocycloid that acquires a teardrop shape at n 1, This fractal property of the row and entire permanental functions is reminiscent of the properties of the Peano and similar fractal curves and illustrated in Figure 7. 1 (u, v) and ρ 1 (u, v), respectively) with a support on the aforementioned finite region of the complex plane z = u + iv. For the row function of the matrix A pq = 1, it is shown in Figure 8 for the matrix size n = 5, 6, 8. The distribution function ρ 1 (u, v) for the entire permanental function is shown in Figure 9.
The analysis of the fractal properties of the permanental functions entering the permanent's integral representation in Equation (76), that was presented above for the case of the matrix A pq = 1, could be easily extended to the matrices A pq = a −q and A pq = 1 + aδ p,q , who's permanents n!/a n(n+1)/2 and e a Γ(n + 1, a), as per Equation (142), are known, as well as to other circulant matrices with still unknown permanents. We skip it here.
We only state that the related row and entire permanental functions possess the similar fractal properties and point to some of their typical modifications due to a variation of the matrix entries.
For instance, Figure 10 shows the 2d probability density function ρ (1) 1 (u, v) of the scaled fractal first-row function (91),B 1 /A 1 , for the n × n circulant matrix with a power-law decay of its first-row entries, A 1q = q −2 , and n = 7. A comparison with Figure 8 clearly proves that the aforementioned fractal behavior in the present case of the circulant matrix with varying entries qualitatively remains the same as in the case of the constant matrix A pq = 1 discussed above. Just its sharp features become smoother, that is, less pronounced, and the borders and topology of the fractal support region on the complex plane are modified. In particular, an increasing variation of the matrix entries results in an appearance of a no-support, empty region emerging in the central part of the fractal's support as is shown in Figure 10 and also in Figure 12 below (cf. Figure 7). The inner and outer borders of the permanent's fractal support could be found analytically, via a method of the Lagrange multipliers, for any complex matrix A. The related results will be presented elsewhere.    + 1, a), with increasing range of the integration T for the circulant n × n matrix A pq = 1 + aδ p,q with varying entries by plotting the relative value of the integral in Equation (79), I A (T)/perA, with the base b = e = 2.718... as a function of T in Figure 11. The convergence is stable, similar to that in Figure 4.

Multivariate Representations of the Matrix Permanent
Here we find a relation of the permanent to complex random polynomials and determinants of many variables.

The Integral Representation of the Permanent via a Multivariate Polynomial of Complex Variables:
HereB p ({z q }) = 1 z ∑ n q=1 A pq z q , z q = e iθ q ,z n = ∏ q z q . A validity of the formula (96) stems from the construction of the permanental functionP A in such a way that all the terms in Equation (96) which are present in the perA do not depend on the phases θ q in virtue of the imposed phase-locking factor ∏ n q=1 z q z = 1 and, hence, are not affected by the integration. All the other terms contain, at least, one factor z q = e iθ q and, hence, vanish after the integration.

Discrete Analogs of the Permanent's Integral Representations: BBFG Formula & Its Generalization
The method that we employed for the derivation of the permanent's integral representations in Equations (76) and (96) immediately yields also their discrete analog-a formula of Balasubramanian [57], Bax-Franklin [58] & Glynn [23] where δ q=n = 1. Namely, it suffices to replace the complex variables z q = e iθ q running over the unity circle in Equation (96) by the discrete variables δ q running over two values ±1 and then, instead of imposing the condition ∏ n q=1 z q z = 1 that selects all summands constituting the permanent to be independent on the variables z q and hence not to vanish after applying the integration, multiply the product of the matrix-row sums by the factor ∏ n k=1 δ k . The latter does the same job of selecting all summands constituting the permanent to be sign independent on the discrete variables δ q = ±1, while making any summand irrelevant to permanent the antisymmetric function of those variables δ q = ±1 the column's index q of which is missing in that summand. As a result, the summation ∑ {δ q } in Equation (97) nullifies all these irrelevant summands and yields the exact value of the permanent.
The proof of Equation (97) presented above immediately allows us to conclude that there is an entire series of different discrete-sum representations of the permanent for any integer m = 2, 3, ..., where the variables δ q run over the set R m of all m-th roots of unity, that is δ m q = 1 but ∑ δ q ∈R m δ q = 0. The BBFG formula (97) is a particular case of the result in Equation (98) at m = 2.

Permanent vs. Determinant: The MacMahon Master Theorem
Note that we employed in Section 2 another novel representation of the permanent -via a determinant and the MacMahon master theorem. Here it is given in the integral and discrete forms.
Namely, we can rewrite a well-known MacMahon master theorem [54] as an integral representation of the permanent of any matrix A via a n-dimensional Fourier integral where z = diag{z k |k = 1, ..., n} is a diagonal matrix with the entries z k = e iθ k . The discrete sum representation is where the diagonal matrix δ = diag{δ q |q = 1, ..., n} has the entries δ pq = δ q δ p,q determined by a set of the integers {δ q = ±1|q = 1, ..., n} and the Kronecker delta δ p,q . The result (100) stems from the MacMahon master theorem written in a form of the n-dimensional partial derivative taken at the zero values of the variables z k , k = 1, ..., n.

Permanent's Fractal vs. Complex Stochastic Multivariate Polynomial
The chaotic fractal behavior of the permanental function can be analyzed and understood via the multivariate polynomial of the complex variables comprising the permanent's integral representation in Equation (96). Here we illustrate this method by two examples.
First, we show that the support region of the fractal associated with the row function in Equation (75) (cf. Equation (91)) on the complex plane is given by a range of the row-sum multivariate polynomial B p ({z q |q = 1, ..., n}) in Equation (96) for any n × n matrix A pq . The row polynomial is a complex-valued stochastic variable given by the function of n random phases θ q , q = 1, ..., n, uniformly distributed over their domains θ q ∈ [−π, π], A fact is that mapping the n-dimensional domain of variation of a multivariate random vector {θ 1 , ..., θ n }, that is, the n-dimensional direct sum of the n intervals [−π, π], or [−π, π] n , to the complex plane via the row functionB p ({z q }) in Equation (102) yields the support of the fractalB p (e iπt ) in Equation (77), which is a function of just one real variable t. This is shown in Figure 12. Second, it is possible to find the ergodic measure of the permanental-function fractal in Equation (77) by calculating the 2d probability density function, ρ n (u, v), of the multivariate polynomial in Equation (96), . The latter can be found by mapping a distribution of the product of n random row functionsB p ({z q }), p = 1, ..., n, from the uniform distribution of the multivariate random phase vector {θ 1 , ..., θ n } in the n-dimensional cube [−π, π] n to the complex plane of the values ofP A ({z q }) = u + iv. For instance, for a n × n matrix A of a large size n 1 and a moderate variation of entries, we find its reasonable zeroth-order approximation analytically as follows It is obtained as a distribution function of the n-th power of one row functionB 1 ({z q }) by neglecting correlations between different row functionsB p ({z q }) as well as correlations imposed by the phase factor z n . The zeroth-order approximation for the distribution function of the row functionB p ({z q }) is Here J 0 is a Bessel function, r = √ u 2 + v 2 . Equations (104) and (103) approximate the main, independent on a polar angle, part of the distributions ρ (p) n (u, v) and ρ n (u, v) quite well. They are shown in Figures 10 and 13 along with the 2d probability density functions ρ (p) 1 (u, v) and ρ 1 (u, v) of the fractal row and entire permanental functions in Equation (76),B p (e iπt ) andP A (e iπt ), respectively. In the left three quarters of Figure 13, the approximation (103) (light gray) is a bit larger and, hence, shields the actual distribution (dark gray). In the right quarter, the actual distribution is larger and shields the surface (103).
Importantly, we verified that the fractal's distributions ρ  Thus, the matrix permanent could be calculated as a mean value of the stochastic permanental function or polynomial by averaging their values over the complex plane with the weight given by the distribution function ρ 1 (u, v) or ρ n (u, v), respectively: So, finding the next-to-zeroth and higher-order approximations for the distribution function ρ 1 (u, v) of the permanent's fractal for the 1d integral representation in Equation (76) or ρ n (u, v) for the permanent's n-dimensional integral representation in Equation (96) allows one to compute the permanent via Equation (105).

Manifestation of a Number-Theoretic Complexity in the Permanent of Schur/Fourier Matrices
Let us relate the complexity of the critical phenomena in physics of many-body systems to the number-theoretic complexity in mathematics via the permanent of a n × n Schur matrix S pq = e 2πipq/n , employed in a fast Fourier transform, or the degenerate Schur matrices S ν , which differ from S by a replication or deletion of some columns.

Permanent's Representation via Laplace Integrals
Applying the Binet-Cuachy expansion to a circulant n × n matrix A = PSΛS † P −1 /n, like the one in Equation (2), written via the diagonal matrices of phases P = diag(e 2πip/n ) and A's eigenvalues Λ = diag(λ p ), p = 1, ..., n, we find [59] For the perS ν , we get the Laplace integral representation where a coefficient η and a polynomialf Q of Q variables q k and integers m k (k = 1, ..., Q) defined via a generalized multivariate q-Pochhammer symbol [60,61] are specified right after Equation (110). The results (106) and (107) yield many links between the number theory and permanents. Say, a number of the degenerate Schur matrices with nonzero permanent and a number of terms in the permanent of the generic circulant matrix, both equal ∑ d|n ϕ( n d ) (2d−1)! d!(d−1)!n , are related to the Euler's totient function ϕ(k). The permanent of the Schur matrix of the odd order n [62], perS = ∑ d|n u n (d)¯( n d ), is related to the Möbius function µ(k) and a cardinality u n (d) of a subset of permutations {π ∈ S n | ∑ n l=1 lπ(l) = d mod n}. A contribution of the first two orders (Q = 1, 2) in Equation (107) to this permanent, perS = n(n − 2)!(H n−2 − 2 + D n ) + {terms of order Q ≥ 3}, includes the harmonic number H m = ∑ m k=1 1/k and another number-theoretic function D n involving the greatest common divisor (n, m). The latter is related to the Euler's totient function and Ramanujan's sum via a discrete Fourier transform. At last, Equation (107) is equivalent to a sum over multiset partitions, a major focus of the number theory and combinatorics.
More links between the number theory and the permanent could come from Equation (98) which, like the Möbius function and the Ramanujan's sum, involves a sum over the roots of unity.
Thus, the complexity of the major number-theoretic functions is closely related to the complexity of the permanent.

The Permanent of the Schur/Fourier and Circulant Matrices vs. the Number Theory
Here we brief on the details of the rigorous definitions and properties of the circulant and Schur/Fourier matrices, polynomials and number-theoretic functions relevant to the permanent's representations in Equations (106) and (107).
The n × n circulant matrix A with the p-th row and q-th column entries A pq is a Toeplitz matrix [51] with rows obtained via the consecutive cyclic permutations of the elements of the first row. It is given by the discrete Fourier transform of the set of its eigenvalues {λ l | l = 1, ..., n}, λ l e 2πi(q−p)(l−1)/n .
A submatrix [A] {i k |k=1,...,m} of the matrix A is a (n − m) × (n − m) matrix obtained from A by deletion of m rows and m columns which intersect at the diagonal entries A i k i k specified by m integers i k ∈ [1, n]. The Schur matrix S pj = e 2πipj/n is employed in the fast Fourier transform and sometimes is called the Fourier matrix. A n × n degenerate Schur/Fourier matrix S ν is the Schur matrix (S pj ) each, j-th column of which is replicated ν j times, (The multiplicities, ν j ≥ 0, are integers; j = 1, ..., n.) It is specified by a n-tuple ν = (ν 1 , ..., ν n ), ∑ n j=1 ν j = n; θ(x) is a step-function. All J + 1 different columns j = j(i) of S ν are enumerated in the increasing order j(0) < j(1) < ... < j(J) by an index i = 0, 1, ..., J; J > 0. S ν{i k } is the (n − m) × (n − m) matrix obtained from S ν by deleting m rows with indexes p = i k (k = 1, ..., m) and truncating the column-index range to q = 1, ..., n − m.
The Laplace integral representation of the permanent of the degenerate Schur/Fourier matrix in Equation (107) involves a symmetric homogeneous polynomial in Q variables q k of degree n − ν j(0) with Q parameters m k , -actually, its reduced counterpartf Q ({q k , m k }) which is built from the polynomial f Q ({q k , m k }) by keeping only those monomials which simultaneously depend on all Q variables q k ; η = (−1) (n+1)j(0)+n−ν j(0) (ν j(0) !). The polynomial (110) can be written in terms of a special function, that is, a multivariate ν-generalized q-Pochhammer symbol A normal q-Pochhammer symbol of q-analysis [60,61] is The number Q of variables q k in the polynomialsf Q ({q k , m k }) contributing to perS ν in Equation (107) depends on the n-tuple ν of the degenerate Schur/Fourier matrix S ν but, in any case, is bounded from above by an inequality Q ≤ Q max ≤ ∑ n j=1 jν j /n − j(0), (n − ν j(0) )/2. Equation (107) proves that only very few degenerate Schur/Fourier matrices have the nonzero permanent perS ν = 0, namely, the matrices S ν satisfying the following Diophantine equation ∑ n j=1 jν j = 0 mod n. The number-theoretic functions discussed above in conjunction with the permanent perS ν are related to the prime numbers, namely, the numbers coprime to n/d, where d is a divisor of the matrix size n. The Euler's totient function ϕ(k) counts the number of positive integers j up to a given integer k that are relatively prime: j ≤ k and the greatest common divisor (j, k) = 1 is unity. The Möbius function µ(k) = ∑ k j=1;(j,k)=1 e 2πij/k is the sum of the primitive k-th roots of unity and takes on the three values µ(k) = −1, 0, 1. Similarly, the Ramanujan's sum c k (m) = ∑ k j=1;(j,k)=1 e 2πimj/k is the sum of the m-th powers of the primitive k-th roots of unity. The Euler's totient function and the Ramanujan's sum can be calculated via the Möbius function as the following sums and are related to the greatest common divisor (j, n) via the discrete Fourier transform: These facts explain why a close relation between the number-theoretic functions and the permanent [1-4] of the degenerate Schur/Fourier matrices is so natural. It is known that a computational complexity of the Euler's totient function is that of a factoring of an integer into a product of prime numbers. It is believed that the latter problem belongs to a NP-intermediate class, that is, it is much harder than any polynomial, P-class problem, but it is not one of the hardest, NP-complete problems. The factoring of large numbers constitutes a basis for a famous RSA public-key cryptosystem and as such had been studied in detail for more than four decades. There is no efficient, polynomial-time algorithm for solving this problem, although the running time of the best known algorithm, the general number field sieve (GNFS) algorithm, for factoring a b-bit number is sub-exponential ∼O(exp{4[b(log b) 2 /9] 1/3 }).
Apparently, the computational complexity of the number-theoretic functions contributes to the permanent's computational complexity. At the same time, the latter includes also other factors. Even for the quite special Schur/Fourier matrix S, computing the permanent [62], perS = ∑ d|n u n (d)¯( n d ), requires computing a cardinality u k (d) of a subset of permutations {π ∈ S n | ∑ n l=1 lπ(l) = d mod n}, along with the Möbius function µ(n/d). Moreover, even a contribution of the first two orders (Q = 1, 2) in Equation (107) involving the greatest common divisor (n, m) via the binomial coefficients C q p = p!/(q!(p − q)!) and shown in Figure 14  For the degenerate Schur/Fourier matrices S ν specified by the n-tuples {ν 1 , ..., ν n } or the circulant matrices specified by the eigenvalues {λ 1 , ..., λ n } (see Equation (106)), the permanent's computational complexity greatly increases due to a necessity to compute perS ν for a very large number r(n) = ∑ d|n ϕ( n d ) (2d−1)! d!(d−1)!n of the different n-tuples ν generating the nonzero permanents [63] perS ν , even though this number is much less than the total number of all degenerate Schur/Fourier matrices T n = (n + 1)C n /2, where C n = (2n)!/[(n + 1)!n!] is the Catalan number. This additional complexity is encrypted into the polynomial (110) of the permanent's integral representation (107) via the other special function-the multivariate ν-generalized q-Pochhammer symbol (112).
An overall complexity of the degenerate Schur/Fourier matrix permanent can be understood as a combinatorial, number-theoretic complexity of the multiset partitions ν Q = {ν (k) |k = 1, ..., Q} of the n-tuple ν = {ν 1 , ..., ν n } which constitute the matrix's discrete representation [59] Here Q submultisets ν (k) = {ν They are naturally ordered, ν (k) ν (k+1) , in accord with the decreasing order ν of their components. The lengths d l (l = 1, ..., L) of the intervals of equal consecutive submultisets ν (k l ) = ν (k l +1) = ... = ν (k l +d l −1) are called the degeneracy factors of the multiset partition ν Q . Each submultiset The sum in Equation (117) runs over those multiset partitions ν Q for which all submultiset partitions ν (k) have a span Thus, the permanent, perS ν , of the degenerate Schur/Fourier matrix, both in the form of the Laplace integral representation in Equation (107) and in the equivalent form of the discrete representation in Equation (117), is the sum over the multiset partitions which constitute a major focus in the number theory and combinatorics.
In addition to the aforementioned functions related to perS ν , the other combinatorial and number-theoretical functions enter the scene when one calculates the permanent of the circulant matrix via the power expansion over the matrix eigenvalues in Equation (106). This happens because the latter includes the sum over the n-tuples ν = {ν 1 , ..., ν n } specifying the degenerate Schur/Fourier matrices, in addition to the aforementioned sum over the multiset partitions. In particular, the number of terms with a given value of the exponent ν 1 in λ ν 1 1 is equal to the well-known rencontres (encounter) numbers [53] D n,ν 1 = (!(n − ν 1 ))C ν 1 n = n!
x means rounding x up for even n and down for odd n.

Asymptotics of the Permanent and the Szegő Limit Theorems
The result (2) points to a fundamental open problem of finding the permanent's asymptotics for the case of the large-size circulant matrix and its analogy with the Szegő limit theorems on the Toeplitz determinant employed in the Onsager's exact solution of the 2d Ising model [47][48][49][50][51][52].
A starting point for finding this asymptotics could be the expansion (106) of perA in powers of the A's eigenvalues.
An interesting example of an asymptotic reduction of the permanent of a doubly stochastic n × n matrix A, with a moderate variation of its entries, to the determinant, was found in Reference [64]; I pq = δ p,q , J pq = 1/n, A stands for the transpose of A. It works quite well, for example, for the nonnegative circulant matrices with a power-law, a q ∼ q −ξ , or similar moderate variation of entries. Other ways to estimate the permanent have been recently reviewed in Reference [4]. For the symmetric circulant matrix A, the matrix under the determinant in Equation (119) is circulant, that is Toeplitz. Below we show that the permanent's asymptotics (119) is directly related to the Toeplitz-determinant asymptotics.

The Circulant Determinant vs. the Toeplitz Determinant
Finding the permanent's asymptotics for the large-size circulant matrix remains an open problem who's solution would be the key to a practical application of the permanents in the theory of critical phenomena and quantum many-body processes as well as related studies of the nature's complexities. An analogous problem of finding the asymptotics of the Toeplitz determinant was the key to the Onsager's solution of the 2d Ising model [47][48][49][50][51][52]. The latter problem had been solved by the first and second Szegő limit theorems on the asymptotics of the Toeplitz determinant. Here we elaborate on this analogy in view of the McCullagh permanent's asymptotics [64] in Equation (119) applied to the symmetric circulant matrix A. In this case the problem is reduced to finding the determinant's asymptotics of the circulant and, hence, Toeplitz matrix I + J − A A where A is the transpose of A. This fact sets a direct relation between the permanent's asymptotics and the Szegő limit theorems.
This relation is based on an important fact that the determinant of the n × n circulant matrix C n = (c q−p mod n ) can be explicitly computed as a product of its eigenvalues given by a discrete Fourier transform of the first row of the circulant matrix, (C n ) 1q = c q−1 , q = 1, ..., n.
Assuming that all of the eigenvalues are not zero, one has det C n = exp Here an associate polynomial λ (n) (z) of the circulant matrix C n is taken on the unit circle in the complex plane, z = e −2πix , at the discrete, homogeneously distributed values of a polar angle 2πx corresponding to the values x l = l/n, l = 0, 1, ..., n − 1 of a real variable x ∈ [0, 1].
The exact Equation (121) is a discrete analog of the first Szegő limit theorem [51]. The latter gives a leading term of the determinant's asymptotics for the Toeplitz n × n matrix T n = (t q−p ) at n → ∞ via a mean value of the logarithm of its associate polynomial, symbol t(z), on the unit circle: This leading asymptotic behavior corresponds to the homogeneous distribution of the eigenvalues on the unit circle.
The next-to-leading term in the Toeplitz-determinant asymptotics, namely, the independent-on-n pre-factor in front of the leading term (that, according to Equation (122), has an exponent growing linearly with n) is given by the second (strong) Szegő limit theorem as follows lim n→∞ det T n e n 1 where (log t) k = 1 0 e −2πikx log[t(e 2πix )]dx is the k-th Fourier coefficient of the logarithm of the symbol t(z) associated with a limit Toeplitz matrix T = lim n→∞ T n .
On the one hand, an appearance of the nontrivial factor (123) due to the strong Szegő theorem manifests a more complex structure of the Toeplitz matrices in the general case compared to the structure of a special subset of the Toeplitz matrices -the circulant matrices. Essentially, a product of two circulant matrices is always the circulant matrix, while the Toeplitz matrices in the general case lack this property of a multiplicativity, that is, do not form a multiplicative group [65]. On the other hand, among the Toeplitz matrices there is a wide subset of matrices for which the symbol t(z) is well defined, for instance, when a well-known sufficient condition for a convergence of its series in Equation (122), ∑ k |t k | < ∞, is satisfied due to a relatively fast decrease of the entries t ±k at k → ±∞. The Szegő limit theorems are directly applicable only to this subset of the Toeplitz matrices.
For the circulant matrices C n , a definition of the associate Toeplitz symbol employed in Equation (122) is not appropriate since a series t based on the associate polynomial λ (n) (z) introduced in Equation (121) could be used, but only if the sequence c k is decreasing so that, for instance, the convergence condition ∑ k |c k | < ∞ is satisfied. However, it would exclude from consideration the basic case of the periodic boundary conditions for the phase transition problem when the first row of the circulant matrix, that is, the correlation function, has the same values at the symmetrically located positions, |c k | = |c n−k |. In the latter case, let us define a circulant-matrix symbolλ(z) via a modified associate polynomial which monomials z k acquire a shift −n of the exponent k if it is greater than an integer part of n/2, The latter preserves the circulant determinant (121): The most interesting, nontrivial situation of passing the critical region of a phase transition corresponds to the case when the related circulant matrix entering the result in Equation (2) and discussed in Section 2 above ceases to have a convergent symbol, neither λ(z) norλ(z). This occurs in the critical region of parameters when the correlation function representing the first row of this circulant matrix experiences a transition from an abrupt exponential decay to spreading over an entire macroscopic volume of the system. In this case computing the asymptotics of the circulant determinant amounts to finding the asymptotics of the series in the exponent of Equation (121) or Equation (126). Note that even then this series could be represented sometimes as an integral from a smooth function over the unit interval x ∈ [0, 1]. The latter function may be different from the logarithm of the polynomial λ (n) (z) orλ (n) (z) associated with the circulant matrix as per Equation (121) or Equation (126) since this polynomial is strongly oscillating, ill-defined at n → ∞. An example is given below.
The aforementioned nontrivial, complimentary relation between the asymptotics of the circulant determinant in terms of the explicit exact formula in Equation (121) or Equation (126) and the asymptotics of the Toeplitz determinant in terms of the Szegő limit theorems could be clarified by considering a typical case of the circulant matrix being represented as a sum of two Toeplitz matrices, C n = T n +T n , who's entries t ±k ort ±k generally decrease in value with increasing k or n − k from 1 to n − 1, respectively. Then, the determinant of the circulant matrix C n differs from the determinant of the Toeplitz matrix T n by a factor: det C n = det T n det(1 +T n T −1 n ).
Let's consider a symmetric circulant matrix with the entries c k = c n−k decreasing in absolute value when k runs from 1 to [n/2]. Let the Toeplitz matrix T n contains the central (1 + 2[n/2])-band of the circulant matrix C n and zero entries in the complimentary to this band upper and lower triangular parts, that is, t k = c k for k = 0, 1, ..., [n/2], t k = c n+k for k = −1, ..., −[n/2], and t k = 0 for |k| > [n/2]. Suppose ∑ k |t k | < ∞, so that the Toeplitz symbol t(z) in Equation (122) exists and the Szegő limit theorems in Equations (122) and (123) apply. Then, since the Toeplitz-matrix symbol t(z) coincides with the circulant-matrix symbolλ(z) as per Equations (125) and (126), the circulant determinant in the left hand side of Equation (127) coincides with the leading term of the Toeplits-determinant asymptotics given by the first Szegő theorem (122). Obviously, the second (strong) Szegő theorem (123) gives the pre-factor 1/ det(1 +T n T −1 n ) that accounts for a fact that computing the Toeplitz determinant cannot be reduced to computing a mean value of the logarithm of the symbol over the homogeneous distribution of eigenvalues on the unit circle and, hence, is more complex than computing the circulant determinant. Only in the case when a contribution due to the complimentary Toeplitz matrixT n becomes negligible, the strong Szegő theorem (123) is reduced to just a trivial, unity factor and, hence, the Toeplitz determinant is reduced to a value given by the first Szegő theorem and coinciding with the exact result in Equation (121) or (126) for the circulant determinant.

McCullagh Asymptotics of the Permanent and Two Opposite Limits for the Circulant Determinant
Now we are ready to unfold the permanent's asymptotics (119) at large n for the doubly stochastic matrix A n = C n /λ 1 with a moderate variation of its entries, n n perA n n! ∼ 1 The eigenvalues of the matrix A n = C n /λ 1 differ from that of the matrix C n , Equation (120), only by the scaling factor ∑ n−1 k=0 c k = λ 1 . Each eigenvalue of the matrix I + J − A n A n entering the McCullagh asymptotics equals to the sum of the eigenvalues of matrices I, J, and −A n A n since all of these matrices are circulant in the case of the symmetric circulant matrix A considered here and, hence, can be diagonalized by the same discrete Fourier transformation. Thus, the determinant of the matrix I + J − A n A n is given by the product of its eigenvalues: Here we implemented the fact that only the first eigenvalue of the matrix J is not zero, namely, it equals unity. As a result, calculating the product in Equation (129) via the exact eigenvalues in Equation (120), we get the explicit asymptotics of the permanent by means of Equation (128). Since the entries c k and the eigenvalues λ l of a sequence of the circulant matrices {C n } constitute two counterparts of the discrete Fourier transform as per Equation (120), there are two opposite limiting cases. Namely, either (i) the matrix entries given by the function c k of the integer variable k enumerating matrix columns are well localized and do not spread over an entire range of the column index k with increasing matrix size n, for example, when the conversion condition ∑ k |c k | < ∞ is satisfied, or (ii) the eigenvalues λ l are similarly well localized and do not spread over an entire range of the eigenvalue index l with increasing matrix size n.
In the first case, the limiting distribution of the eigenvalues λ l could be described by a smooth, well-defined function, say, λ(e 2πix ) orλ(e 2πix ), of the variable x = l/n ∈ [0, 1]. In the second case, the limiting distribution of the entries c k could be described by a smooth, well-defined function of the variable x = k/n ∈ [0, 1].
This alternative is a manifestation of the uncertainty principle: The more concentrated a function is, the more spread out its Fourier transform must be. In the application to the critical phenomena in a spontaneous symmetry breaking, a transition from the first case to the second case corresponds to the transition from a disordered phase with a strongly localized correlation function to an ordered phase with the correlation function spread out over an entire, macroscopic dimension of a many-body system. The complexity of the critical phenomena revealed by the matrix permanent's complexity arises due to a simultaneous absence of such smooth, well-defined functions in the both dual domains of the matrix entries and the eigenvalues that takes place in the central part of the critical region of phase transitions.
Let us illustrate a related transition in the behavior of the permanent's asymptotics by considering a sequence of the symmetric circulant n × n matrices {C n } who's first rows c k are given by smooth functions c (n) (x) of the continuous variable x ∈ [0, 1] at the discrete points x = k/n. In accordance with the Poisson summation formula for a discrete-time Fourier transform (DTFT), the eigenvalues are given by Equation (120) as a periodic summation of Fourier coefficientsc (n) (j) of the function c (n) (x) as follows The determinant of the matrix C n equals their product: In a class of the functions c (n) (x) = h(β n x) which are given by one and the same smooth function h(x), x ∈ [0, ∞), and depend on the matrix size n only via a scaling factor β n , the first of the two opposite limiting cases discussed above can be illustrated by the following typical case of the scaling β n = n. If the function h(x) is decreasing fast enough at large values of |x|, then, in the limit of large n, the matrix C n becomes essentially a relatively narrow band matrix and the entries c k ≡ c (n) (x = k/n) = h(k) constitute the Fourier coefficients of the well-defined symbolλ(z), Equation (125). So, the associated polynomials λ (n) (z) of the matrices C n could be considered as projections of the symbolλ(z) onto a subspace of polynomials spanned by the finite-power monomials z k , k ∈ [[ n 2 ] + 1 − n, [ n 2 ]]. If the symbolλ(z) is a smooth positive function, so that the projectionsλ (n) (z) tends toλ(z) with increasing n, then the eigenvalues λ l tend toλ(e 2πil/n ) in a dense and homogeneous set of points on the unit circle. In this case, the leading contribution to the determinant in Equation (130) is provided solely by the Fourier coefficients within the central n-period of harmonics, l + nj ∈ [−n + 1 + [ n 2 ], [ n 2 ]], and the asymptotics of the determinant is given by the integral of the logarithm of the symbol over the unit circle: This formula for the determinant of circulant matrices has the same form as the first Szegő limit theorem in Equation (122). Note, however, that the analogous asymptotics of the Toeplitz determinant includes an additional nontrivial pre-exponential factor given by the strong Szegő theorem, Equation (123). The reason for this difference was discussed above, after Equation (127), and stems from the fact that the circulant and Toeplitz matrices are close to each other only in terms of the weak or Hilbert-Schmidt norm, but differ significantly in terms of the strong norm [51]. Hence, their determinants have different limiting values.
The second of the two opposite limiting cases discussed above can be illustrated by the case of the unity scaling factor β n = 1. Now, with increasing n, the entries c k ≡ c (n) (x = k/n) = h(k/n) of the circulant matrices C n tend to the values of the smooth function h(x) on a dense homogeneous set of points {k/n} within the unit interval x ∈ [0, 1]. In virtue of the uncertainty principle, it implies that the eigenvalues λ l determined by the entries' Fourier transform as per Equation (130) appear to be strongly localised in the dual, frequency domain. In the case of the symmetric doubly-stochastic matrix A n = C n /λ 1 with the entries c k decreasing from the side columns k = 0 and k = n towards the central column k = [n/2], which corresponds to the problem of phase transition in a system with the periodic boundary conditions and is relevant to the McCullagh permanent's asymptotics, the generating function h(x) is real-valued and positive. So, at n → ∞, all of the significant eigenvalues are concentrated only near the boundaries l = 1 and l = n, that is, at l n and n − l n, as the symmetric pairs λ l = λ n−l . In the leading order, their asymptotics is given by just one relevant termc (n) (l) =c (n) (n − l) in the periodic summation in Equation (130). The rest of the eigenvalues, spanning the entire central part of the range of Fourier harmonic numbers l between these boundary layers, tend to zero. As a result, the McCullagh permanent's asymptotics in Equations (128)  (133)
A straightforward calculation of the product in Equation (129) yields the explicit asymptotics of the permanent via Equation (128) for the large finite values of n. Here we present it only for n 1 + α when the determinant gets to the limit lim n→∞ det(I + J − A n A n ) = α 2 sinh 2 (α √ 2) 2 sinh 4 (α) (136) and the leading term of the permanent's asymptotics is lim n→∞ n n perA n n! = The latter result follows from Equation (133), that is, calculation of the product in Equation (129) as (∏ r k=2 ) 2 for r = [n/2] → ∞, by means of the integration over a parameter p ∈ [0, 1], Dependence of the permanent on the matrix size n (the exact values calculated numerically and shown by dots) and its McCullagh asymptotics (128) (shown by crosses) are illustrated in Figure 15 for the doubly stochastic circulant matrix A n = C n /λ 1 specified in Equation (134) in the case of the relatively large correlation decay parameter, α = 2.5. This parameter sets a scale of the entries' variation from the maximum value c k=0 = cosh(α) 1 at the main diagonal to the minimum value at a half way to the matrix edge c k=[n/2] that is 1 for even n and about 1 for odd n > α. The fact that this scale of variation is a fixed, independent on n quantity makes the matrix A n resemble the matrix J at n → ∞ and, hence, justifies McCullagh asymptotics which approximates the permanent quite well starting already from n ∼ 20. At the same time, for relatively small n, a steep, exponential decay of matrix entries with a deviation from the main diagonal makes the matrix A n resemble the narrow band matrix that validates an opposite-type approximation of the permanent as the product of the diagonal entries, associated with a random phase approximation [33]. Thus, the example shown in Figure 15 simultaneously illustrates both opposite limiting cases (i) and (ii) of the circulant permanent behavior discussed above and a nontrivial transition between them. An unambiguous agreement of the exact numeric calculations with corresponding analytic results in Equations (140) (a dashed-dotted curve) and (137) (a dashed line) at small (n < 5) and large (n > 20) matrix sizes n, respectively, is remarkable.  Figure 15. The scaled permanent, n n perA n /n!, of the doubly stochastic circulant n × n matrix A n = C n /λ 1 specified in Equation (134) as a function of the matrix size: the dots-an exact numerical calculation, the crosses-the McCullagh asymptotics in Equation (128), the dashed-dotted curve-the diagonal (random phase) approximation in Equation (140), the dashed line-the leading asymptotics in Equation (137); α = 2.5.

The Permanent's Asymptotics for the Circulant Matrices with the Power-Law Varying Entries
A similar analysis of the permanent for the circulant matrices with a power-law variation of the entries, which is characteristic for emergence of the ordered phase in the critical region within the renormalization-group approach and is significantly less steep than the exponential variation in Equation (134) discussed above, also confirms the aforementioned two opposite types of the permanent behavior. Here, for the sake of room, we skip it along with an illustrative plot of the n-dependence of the permanent which is similar to the one shown in Figure 15.
Its first eigenvalue is λ 1 = c 0 − c 1 + nc 1 . All other eigenvalues are equal to λ k = c 0 − c 1 , k = 2, ..., n. We find its permanent analytically using combinatorics and the well-known rencontres number !k of the derangements (118), that is, the permutations of the set {1, 2, ..., k} without fixed points: where Γ is an upper incomplete gamma function. Introducing a relative variance of the diagonal and off-diagonal entries, α n = c 0 c 1 − 1, which plays a part similar to that of the correlation decay parameter α in Equation (134), we find the exact analytic formula for the permanent of the related doubly stochastic circulant matrix A n = C n /λ 1 , perA n = e α n Γ(n + 1, α n ) (n + α n ) n , and its McCullagh asymptotics in Equations (128) and (129), perA n ∼ n! n n 1 − 1 (1 + n/α n ) 2 (1−n)/2 at n → ∞.
Again, there are two opposite limiting cases. If the entries' variance parameter a n is very large, then the matrix A n is close to the diagonal matrix and its permanent is given by the diagonal (random phase) approximation perA n ≈ c 0 λ 1 n = α n + 1 α n + n n .
In the opposite case of small values of the variance parameter α n , the matrix A n is of the J type and the McCullagh asymptotics in Equation (144) applies. A transition between these two limit regimes occurs at the moderate values of α n . It is described by the exact analytic result in Equation (143) and is similar to the one discussed in Figure 15. Importantly, the exact solution in Equation (143) allows one to analytically obtain all of the permanent's asymptotics and limits for the matrix (141) via the known asymptotics and limits of the incomplete gamma function.
A comparative analysis of various permanent's asymptotics will be presented elsewhere.

Conclusions
We infer that the method of a matrix permanent can be employed as a universal tool for analyzing and measuring complexity in nature. A nontrivial reduction of the critical phenomena, fractals, many-body processes in quantum systems and computing of NP-and P-problems to the permanent as well as the permanent's integral representations presented above clearly demonstrate the power and capability of the permanent-based analytic technique for the unification of the nature's complexities.
Of course, there are other functions, different from the matrix permanent, which are also involved into various P-or NP-complete problems [66]. However, in view of a number of the aforementioned results and facts, the permanent is marked by its universal and central role in describing complexity both in physics and mathematics.

Conflicts of Interest:
The authors declare no conflict of interest.