{'fraction_non_alphanumeric': 0.04597540525433203, 'fraction_numerical': 0.01239718917192366, 'mean_word_length': 4.259107611548556, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The factorization method, self-similar potentials and quantum algebras

This is a brief review of the Schrodinger's factorization method and its relations to supersymmetric quantum mechanics and its nonlinear (parastatistical, etc) modifications, self-similar infinite soliton potentials, quantum algebras, coherent states, Ising chains, discretized random matrices and 2D lattice Coulomb gases.


Introduction
From the very beginning quantum mechanics served as a rich source of good mathematical problems. It played a major role in the development of the theory of generalized functions, functional analysis, path integrals, to name a few. In the last two decades the "quantum" disease became so widespread in mathematics that it is difficult to guarantee that a randomly chosen mathematical term will not get its cousin with such an adjective in the foreseeable future (if it does not have already). The influence of quantum mechanics upon the theory of special functions is also indispensable. A bright example is given by the angular momentum theory which has lead to the Racah polynomials. An aim of this lecture is to outline another fruitful interplay between these two scientific fields inspired by the factorization method. This method was suggested by Schrödinger as a convenient operator language tool for working with quantum mechanical spectral problems [42]. It was reformulated as a problem of searching of the factorization chain solutions by Infeld [25]. The review [26] became a basic reference in this field. A fairly recent revival of the interest to this method occurred due to the discovery of its relation to the notion of supersymmetry (see, e.g. the review [22] and references therein). The author himself turned to this subject from the quantum chromodynamics due to an idea of a generalization of the supersymmetric quantum mechanics [40]. Although it was clear that special functions play an important role in this formalism, it took some time to recognize that, heuristically, special functions are defined in this approach as the functions appearing from similarity reductions of the factorization chain (for a more precise formulation, see the Appendix).
From the mathematical side, a related techniques for solving linear ordinary differential equations of the second order was proposed by Darboux long ago [12]. Its various generalizations are referred to in the modern theory of completely integrable systems as the Laplace, Darboux, Bäcklund, dressing, etc transformations [1]; in the analysis of isomonodromic deformations one deals with the Schlesinger transformations [29]. In the theory of special functions such transformations appear as contiguity relations. An important step in the development of the subject was performed in a series of papers by Burchnall and Chaundy [9], which contain even some parts of the operator formulation of the approach ("the transference"). For a completeness, let us mention also the terms "shift operator" and "transmutation", which are used in some other variations of the formalism. For a rigorous mathematical treatment of some aspects of the Darboux transformations technique or the factorization method, see [11,14,15,30,41].
As far as the discrete recurrence relations are concerned, actually, it is the Christoffel's theory of kernel polynomials that provides a first constructive approach of such kind to spectral problems. This theory is based upon the simplest discrete analogue of the Darboux transformations. A complementary part to this Christoffel's transformation was found by Geronimus in [23]. First applications of the Schrödinger's factorization method to finitedifference equations is given in [37]. Recent developments in this direction are reviewed in [49]. Let us remark that the same techniques was rediscovered in the works on numerical calculations of matrix eigenvalues. More precisely, the well-known numerical LR, QR, galgorithms, etc provide particular instances of the chains of discrete Darboux transformations. The literature on the taken subject is enormous, it is not possible to describe all its branches. The list of references given at the end of this manuscript is not complete, it contains mainly the papers encountered by the author during his own work (additional lists can be found in [47,49]).
Despite of a vast variety of existing constructions, part of which was just mentioned, we limit ourselves in this lecture to the simplest possible case based upon the stationary one-dimensional Schrödinger equation describing the motion of a non-relativistic particle on the line x ∈ R in the potential field u(x), which is assumed to be bounded from below. The operator L is called the Hamiltonian or the Schrödinger operator. For convenience, the particle's mass variable and the Planck's constanth are removed by rescalings of the coordinate x and the energy λ. The equation (1) describes an eigenvalue problem for L and the eigenvalues (or the permitted bound states energies of the quantum particle) are determined from the condition that the modulus of the wave function ψ(x) is square integrable, ψ(x) ∈ L 2 (R).
Depending on the physical situation, the real line R may be replaced by an interval with an appropriate boundary conditions upon ψ(x). However, if one is interested not in the spectra themselves, but in the differential operators L with some formal properties, then it is convenient to take x, λ ∈ C.

The factorization method
Let us factorize the second order differential operator L = −d 2 /dx 2 + u(x) as a product of two first order ones up to some real constant λ 0 : The function f (x) is a solution of the Riccati equation is a solution of the original Schrödinger equation (1) for λ = λ 0 .
If f (x) is a smooth function, then A + is a hermitean conjugate of A − in L 2 (R) and L is a self-adjoint operator. Under these circumstances λ 0 cannot be bigger than the smallest eigenvalue of L. Suppose that λ 0 is the smallest eigenvalue of L, and let ψ 0 (x) be the corresponding eigenfunction. It is well known that ψ 0 (x) is nodeless and may be normalized where a, b are arbitrary constants, is the general solution of the equation Lφ 0 = λ 0 φ 0 . For b = 0, the resulting f (x) is singular at some point and the operators A ± in (2) are not well defined. Let us exclude this situation, i.e. set a = 1, b = 0.
Using the lowest eigenvalue eigenfunction of L, one can always factorize L as described in (2) with the well-defined operators A ± . Vice versa, if one manages to find the factorization (2) such that the zero mode ψ 0 of the operator A − , A − ψ 0 = 0, belongs to L 2 (R) and f (x) is not singular, then λ 0 is the lowest eigenvalue and is the corresponding normalized eigenfunction.
Let us define now a new Schrödinger operatorL by the permutation of the operator factors in (2)L whose potential has the formũ(x) = f 2 (x) + f x (x) + λ 0 . Evidently, one has the intertwining relations playing the key role in the formalism. From (5) one deduces that if ψ(x) satisfies (1), then the functionsψ = A − ψ provide formal eigenfunctions ofL. Indeed, Actually, this gives general solutions of the differential equationLψ = λψ for all λ, except of the point λ = λ 0 , where the zero mode of A − is located. This problem is curable since one can find the general solution of the equationLφ 0 = λ 0φ0 separately: where g, e are arbitrary constants and ψ 0 is the lowest eigenvalue eigenfunction of L (note that A + ψ −1 0 = 0). Denote as λ n ,λ n and ψ n ,ψ n = A − ψ n discrete eigenvalues and corresponding eigenfunctions of L andL respectively. From (6) it follows that the spectra λ n andλ n almost coincidẽ λ n = λ n , n = 1, 2, . . .. The only difference that may occur in these spectra concerns the zero mode of the operator A − . In fact, the point λ = λ 0 does not belong to the spectrum ofL.
Thus, λ 1 is the lowest eigenvalue ofL and the point λ 0 was "deleted" from the spectrum of L. Repeating the same procedure once more, i.e. takingL =Ã +Ã− + λ 1 and permuting the operator factors, one can delete the point λ = λ 1 , etc. This procedure allows one to remove an arbitrary number of smallest eigenvalues of L. Often it is much easier to find the smallest eigenvalue of a given operator than the other ones. If the lowest eigenvalue of L is determined separately by some means, then, by construction, it will coincide with the second eigenvalue of L, etc. This observation is the central one in the factorization method [42] since it reduces the problem of finding complete discrete spectrum of a taken operator to the problem of finding lowest eigenvalues of a sequence of operators built from L by the "factorize and permute" algorithm.
One can invert the procedure of a deletion of the smallest eigenvalue. Namely, a givenL with known lowest eigenvalue λ 1 may be factorized asL = A − A + + λ 0 , with λ 0 < λ 1 . If the zero mode of A − is normalizable, then the operator L = A + A − + λ 0 has the same spectrum asL with an additional inserted eigenvalue at an arbitrary point λ = λ 0 . If one factorizes L (orL) in such a way that f (x) is a non-singular function, but none of the zero modes of the operators A ± are normalizable, then the discrete spectra of L andL coincide completely (an isospectral situation).
There are more complicated possibilities for changing spectral data of a given Schrödinger operator L. E.g., if the first factorization is "bad", in the sense thatL has a singular potential and considerations given above are not valid, then one may demand that after a number of additional refactorizations one gets a well defined self-adjoint Schrödinger operator. In this way one can delete not only the lowest eigenvalues of L or to insert the new ones, but delete or insert a bunch of spectral points above the smallest one [30]. In particular, two step refactorization procedure allows one to imbed an eigenvalue into the continuous spectrum, etc [14].
Let us give now a "discrete time" formulation of the construction. Denote L ≡ L j and L ≡ L j+1 and take j ∈ Z (j may be treated as a continuous parameter and we could denotẽ L ≡ L j−1 -all this is a matter of agreement). This gives an infinite sequence of Schrödinger Neighboring L j are connected to each other via the _ abstract factorization chain Intertwining relations take the form Substituting explicit forms of A ± j into (9) one gets a differential-difference equation upon f j (x): This chain was derived in [25] and the problem of searching "exactly solvable" spectral problems was formulated as a problem of the search of solutions of (10) such that the points λ j define the discrete spectrum of an operator L, say, L ≡ L 0 . E.g., one may try to find solutions of the equation (10) f j (x) in the form of power series in j (λ j are considered as unknown functions of j). As shown in [25,26] the finite term expansion occurs iff f j (x) = a(x)j + b(x) + c(x)/j, where a, b, c are some elementary functions of x. This leads to the 2 F 1 hypergeometric function and well known "old" exactly solvable potentials of quantum mechanics.
As evident from the construction, the constants λ j , j ≥ 0, determine the smallest eigenvalues of L j under the condition that the zero modes of A − j are normalizable and f j (x) are not singular. In general, the j → j + 1 transitions may describe all three possibilities -removal or insertion of an eigenvalue and isospectral transformations (sometimes it is convenient to parameterize inserted eigenvalues as λ 0 > λ 1 > . . . > λ n ).
For a positive integer n, let us introduce the operators The intertwining relations guarantee that for almost all λ solutions of the equations L j ψ (j) = λψ (j) and L j+n ψ (i+n) = λψ (j+n) are related to each other as ψ (j+n) ∝ M − j ψ (j) and ψ (j) ∝ M + j ψ (j+n) . As a result, the product M + j M − j should commute with L j and M − j M + j should commute with L j+n . A simple computation leads to the equalities Let ψ (j) (x) ∈ L 2 (R) be an eigenfunction of L j with the eigenvalue λ and some finite norm Zeroes on the r.h.s. for λ = λ k for some k indicate that the corresponding eigenvalues were deleted from the spectrum. Under the condition that zero modes of A − j , j = 0, . . . , n − 1, are normalizable and nodeless, 0 || = 1, one finds that the functions define the unit norm eigenfunctions of L 0 with the eigenvalues λ n : L 0 ψ (0) n = λ n ψ (0) n and ||ψ (0) n || = 1.
The main advantage of the factorization method consists in its pure operator language formulation. One may replace the Schrödinger operator by any other (higher order differential, finite-difference, integral, etc) operator L admitting factorizations into the well defined operator factors of a simpler nature A ± (for some applications they are not necessarily hermitean conjugates of each other). In all cases one deals with the abstract operator factorization chain (9) with an appropriate interpretation of the constants λ j . It should be noted that this method does not have straightforward generalizations to the multidimensional spectral problems, only some of its features are preservable, see e.g. [3].

Supersymmetry
Supersymmetry is a symmetry between bosonic and fermionic degrees of freedom of particular physical systems. The corresponding symmetry algebras are distinguished from the standard Lie algebras by the presence of anticommutator relations. The simplest superalgebra is realized upon two neighboring operators in the factorization chain. Define the matrix operators (supercharges) where A ± are the factorization operators defined in (2). They form the following algebra (see, e.g. [22]) where the Hamiltonian describes a particle with spin 1/2 (a fermionic variable) on the line in an external magnetic field, σ 3 is the Pauli matrix. In terms of the hermitean supercharges As a formal consequence of the algebra (12) the discrete spectrum of H is doubly degenerate (such degeneracies are characteristic to the supersymmetry) with possible exception of the lowest eigenvalue which cannot be negative (this is another constraint upon the supersymmetric systems).
This simple construction was generalized in [40] to a symmetry between particles with parastatistics. The corresponding symmetry algebras are polynomial in the generators. E.g., in the simplest case involving a parafermion of the second order (or the spin 1 particle) H is given by a 3 × 3 diagonal matrix with the entries L 0 , L 1 , L 2 , the corresponding symmetry generators being (Q + ) ij = A + i δ i,j−1 and its conjugate. For the hermitean charges Q 1,2 defined as above, one has now cubic relations The discrete spectrum of H is now triply degenerate with possible exception of two smallest eigenvalues.
One can go further and propose other modifications of the algebra (12). For instance, one can q-deform it [46]: This lifts the degeneracy of spectra, but creates a nontrivial scaling relation between spectral points. As a result, it opens a way for building various q-harmonic oscillator models. Another construction proposed in [4] refers to a polynomials generalization of (12). It appears after a reduction of the algebra (13) and similar higher order polynomial relations to a twodimensional subspace of eigenfunctions: This is a pure supersymmetry again in the sense of a symmetry between bosons and fermions.
However, the consequences are quite different from the standard case (12). In particular, there are no such severe constraints upon the smallest eigenvalue of H (the vacuum energy).
All these general constructions are useful because they contain new concepts providing nonstandard viewpoints upon physical systems and new mathematical tools for their exploration. Special functions emerge when one starts to work with particular "exactly solvable" models with such symmetries.

Darboux transformations
Let us describe the Darboux transformations technique in its modern appearance. Consider compatibility conditions of the following three linear equations where t is some additional continuous parameter ("evolution time") and u j (x, t), f j (x, t) are some free functions. The compatibility condition of (14) and (15) generates the intertwining where λ j is an integration constant. Thus one arrives again to the infinite chain of nonlinear differential-difference equations (10) and the formal factorizations (8).
The compatibility condition of the equations (14) and (16) yields the celebrated Kortewegde Vries (KdV) equation: a completely integrable Hamiltonian system with an infinite number of degrees of freedom (17), then it takes the form It appears that the compatibility conditions of this equation with the Darboux transformation (15) imposes an additional non-trivial constraint V j = 0, which is called the modified KdV equation.
A substitution of 0 and the Darboux transformation (15) takes the form: k (x, t) are formal eigenfunctions of L j with the eigenvalues λ j+k equal to constants of integration mentioned above for all k. Denote as W (W (φ 1 , . . . , φ n , ξ 1 ), W (φ 1 , . . . , φ n , ξ 2 )) = W (φ 1 , . . . , φ n )W (φ 1 , . . . , φ n , ξ 1 , ξ 2 ), one can show that [11] These formulae give an explicit representation of the transformed potentials u j+n in terms of the initial Schrödinger equation solutions φ  In the formulae given above one may permute any pair of the φ (j) k functions and the final result is not changed, i.e. such a permutation of the Darboux transformations is a symmetry of the factorization chain (10). In the simplest case, it has the following explicit form [2]: all other f j (x), λ j staying intact for j = k, k+1. This discrete symmetry may be "discretized" further by passing to the discrete Schrödinger equation (or the three term recurrence relation for orthogonal polynomials) and the corresponding discrete Darboux transformations [48].
The well known KdV tau-function is introduced as u j (x, t) = −2∂ 2 x log τ j (x, t). In its terms relations (18), (19) are rewritten as τ j+n = W (φ Introducing the variables ρ j = −∂ x log τ j , one can write u j = 2ρ j,x and f j = ρ j+1 − ρ j . As a result, the relation between u j and f j yields the equation which starts to play the role of the factorization chain. The function τ j is a very convenient object since its zeros in x correspond to poles of the potential.

Operator self-similarity and quantum algebras
Let us turn now to the problem of searching particular solutions of the operator factorization chain (9). At first glance it is not clear how to proceed, but the harmonic oscillator problem -a base model for the whole quantum mechanics and quantum field theory -provides a guiding idea. One has to try to form from the factorization operators A ± j some nontrivial symmetry algebras. In fact, the relations (11) look already as defining relations of some algebra, but they are not closed -the relations between the eigenvalue problems for L j and L j+n are too weak (they are valid for any starting potential u j (x)). In order to close the system, one has to assume that there is an additional relation between L j and L j+n , say L j+n = g(L j ), which would force the operators M ± j to map eigenfunctions of a taken operator L j onto themselves.
In the simplest case one demands that the sequence of Hamiltonians L j is periodic: . This is a remarkable fact since the existence of additional conserved quantities may simplify solution of the eigenvalue problem for L j . In the differential operator realization of A ± j this leads to the commuting differential operators [9]. There is a generalization of this pure periodicity condition to the periodicity up to a twist condition L j+N = UL j U −1 , where U is some invertible operator. This leads again to commuting operators [B ± j , L j ] = 0, where but now there is an essential additional freedom in the choice of U.
Another possible "closure" or a reduction of the sequence of operators L j consists in the requirement of their periodicity up to a constant shift and a twist, where µ is a constant. This results in the ladder relation [L j , B ± j ] = ±µB ± j and the operator identities where B ± j operators are defined as in the previous case. Denoting B ± ≡ B ± 0 , L ≡ L 0 , one can form a polynomial algebra where For a representation theory of such algebras, see e.g. [45]. Note that for N = 1 this is the Heisenberg-Weyl or the harmonic oscillator algebra, and for N = 2 it coincides with the su(1, 1) algebra.
Quantum algebras, or q-analogues of the algebras (21) appear from the following operator self-similarity constraint imposed upon the chain (9): When the numerical factor q 2 = 1, one can remove µ by the uniform shift L j → L j +µ/(1−q 2 ).
Such a q-analogue of the Heisenberg-Weyl algebra was encountered in physics long ago [10,18]. In the modern times it became quite popular due to the inspirations coming from the quantum groups, see e.g. [35]. For N = 2 relations (23) determine a particular q-analogue of the conformal algebra su(1, 1) admitting the Hopf algebra structure, etc [46,47].

Self-similar potentials
Consider the Schrödinger equation realizations of the algebraic relations described in the previous section. Let us start from the closureL = ULU −1 + µ, where L = −∂ 2 x + u(x) is the Schrödinger operator in the notations of Sect. 2 and U is a translation operator ) and one gets the Heisenberg-Weyl algebra [L, The potential entering L is defined as For a = 0 one gets from this equation the standard harmonic oscillator model f (x) ∝ x, u(x) ∝ x 2 , which is related to the Hermite polynomials. Although for a = 0 and µ = 0 the author has found a meromorphic solution of (25) in terms of the Weierstrass P-function in general it is quite difficult to build its solutions analytic in some region, especially for µ = 0. Let us remark, that the equation (25) may be derived also after imposing the constraint f j+1 (x) = f j (x + a) upon the chain (10).
Turn now directly to the relations (23). In this case one can take U as the dilation (or qdifference) operator: Uf (x) = |q| 1/2 f (qx). For real q = 0 this is a unitary operator U † = U −1 .
one realizes the identities (23), provided f j (x) satisfy the following system of nonlinear differential-delay equations: Equivalently, these equations appear from the chain (10) after imposing the following con- having a simple group-theoretical interpretation.
One easily arrives at (27) using the Lie's symmetry reduction technique (for a differential context of this theory, see e.g. [38], and for an extension to difference equations, see e.g. [31,36]). First, one notices that if f j (x), µ j are solutions of the chain (10), then their discrete scaling transformation qf j (qx), q 2 µ j gives solutions of (10) as well. Analogously, a shift in the numeration f j (x) → f j+N (x), µ j → µ j+N maps solutions of (10) into solutions. Let us consider a set of self-similar solutions of (10) which is invariant under a combination of these two symmetries. E.g., demanding that these two transformations are equivalent to one another, we arrive to (27). For N = 1, this reduction may be rewritten as f j (x) = q j f 0 (q j x), λ j = q 2j λ 0 in which form it was first met in [43]. The general q-periodic reduction (27) was found by the author [46] in an attempt to q-deform the parastatistical supersymmetry algebras (13). Let q be arbitrary and f 2 j (x) = 1 This yields L 0 = −d 2 /dx 2 or the free quantum mechanical particle, which acquires in this way a q-algebraic interpretation [47].
Suppose that f 0 (x) is not singular at x = 0, then in the crystal base limit q → 0 the potential u 0 (x) boils down to the general KdV N-soliton potential.
The N = 4 case is related to the Painlevé-V function [2], etc. This gives also a commutator representation of the Painlevé equations which was noticed first in [16]. For q = −1 one gets similar situations with the functions f j (x) obeying certain parity symmetry [47]. If N −1 m=0 µ m = 0 then for odd N and some cases of even N the potential u 0 (x) is expressed through the hyperelliptic functions [50,51].
The case when q is a primitive root of unity, q n = 1, q = ±1, is quite interesting [44].
For any odd n and, under certain conditions, for even n potentials are expressed through hyperelliptic functions with additional crystallographic or quasi-crystallographic symmetries.
For example, for N = 1 and q 3 = 1 or q 4 = 1 one recovers the Lamé equation with the equianharmonic or lemniscatic Weierstrass P-functions. Thus, this classical differential equation of the second order appears to be related to the representations of the q-harmonic oscillator algebra for q a root of unity.
The general family of self-similar potentials unifies some of the Painlevé functions with hyperelliptic functions. Due to the connections with quantum algebras, it may be taken as a new class of nonlinear q-special functions (viz. continuous q-Painlevé transcendents) defined upon the differential equations. The standard basic hypergeometric functions [21] appear in this formalism through the coherent states.
As far as the discrete spectrum of self-similar potentials is concerned, from the factorization method and (27) it follows that in the simplest case it is composed from N independent geometric progressions: λ pN +k = λ k q 2p , k = 0, . . . , N − 1, p ∈ N. The same follows from the theory of unitary representations of the algebra (23). The only condition for the validity of this formal conclusion is that the functions f j (x) are non-singular for x ∈ R and positive for x → +∞ and negative for x → −∞.

Coherent states
There are many definitions of quantum mechanical coherent states [39]. In the context of spectrum generating algebras, such as (21) and (23), they are defined as eigenfunctions of the corresponding lowering operators. Such eigenfunctions play the role of generating functions for the space elements of irreducible representations of a taken algebra. Let us describe briefly coherent states for the quantum algebra (23) (for more details, see [47]).
Suppose that the operators B ± are conjugated to each other with respect to some inner product σ|λ . We assume that the spectrum of L is not degenerate. Then the action of operators B ± upon |λ has the form: Let N be odd, 0 < q 2 < 1 and λ 0 < . . . < λ N −1 < q 2 λ 0 < 0. Then for λ < 0 the operator L may have a discrete spectrum formed only from up to N geometric progressions corresponding to the lowest weight unitary irreducible representations of the algebra (23).
This follows from the fact that B − is the lowering operator for the λ < 0 states and N −1 k=0 (λ− λ k ) becomes negative for λ < λ 0 . Since B − |λ k = 0, this problem is avoided for special values of λ, namely, for λ = λ pN +k ≡ λ k q 2p , p ∈ N. For normalizable |λ k , one gets series of normalizable eigenstates of the form |λ pN +k ∝ (B + ) p |λ k .
Coherent states of the first type are defined as eigenstates of B − : where α is a complex variable. Representing |α − as a sum of the states |λ pN +k with some coefficients, one finds where N ϕ N −1 is a standard basic hypergeometric series with the operator argument z = (−1) N αB + /λ 0 · · · λ N −1 and C (k) (α) is a normalization constant (in the set of parameters of the N ϕ N −1 function the term λ k /λ k = 1 is assumed to be absent). The superscript k simply counts the number of lowest weight irreducible representations of the algebra (23) each of which has its own coherent state (or a generating function). The states |α (k) − are normalizable for |α| 2 < |λ 0 · · · λ N −1 |. Coherent states of this type were first constructed for the q-harmonic oscillator algebra (24) in [5], for a particular explicit model of them, see e.g. [6].
It is not difficult to see, that zero modes of the operator L define a special degenerate representation of the algebra (23). Since B ± operators commute with L in the subspace of these zero modes, one can consider them as coherent states as well.
Unusual coherent states appear from the non-highest weight representations of (23) corresponding to the λ > 0 eigenstates of L. Let λ + > 0 be a discrete spectrum point of L.
Then the operators B ± generate from |λ + a part of the discrete spectrum of L in the form of a bilateral geometric progression λ + q 2n , n ∈ Z. Since N −1 k=0 (λ + − λ k ) > 0 for arbitrary N and λ + > 0, the number of such irreducible representations in the spectrum of L is not limited. Moreover, a continuous direct sum of such representations may form a continuous spectrum of L.
An important fact is that for λ > 0 the role of lowering operator is played by B + . Therefore one can define coherent states as eigenstates of this operator as well [47]: Suppose that the bilateral sequence λ + q 2n > 0 belongs to the discrete spectrum of L. Then the states |α + are expanded in the series of eigenstates |λ + q 2n : where 0 ψ N is a bilateral q-hypergeometric series with the operator argument z = αB − /(−λ + ) N and C(α) is a normalization constant. In this case the states |α + are normalizable if Suppose that the λ > 0 region is occupied by the continuous spectrum. The corresponding states |λ may be normalized as σ|λ = λδ(λ − σ). Let N = 1 and the states |λ are not degenerate. Taking ρ = 1 in (24), we get the expansion |α (s) There is a countable family of such coherent states which have a unit norm for |α| 2 > 1/(1 − q 2 ) under the following choice of the normalization constant C(α): The integral under consideration is a particular subcase of a Ramanujan q-beta integral [21].
We were considering until now the abstract coherent states. In the case of Schrödinger equations with the self-similar potentials one has the doubly degenerate continuous spectrum for λ > 0. The structure of |α ± states in this case is very complicated. Even for the zero potential case (the free particle) we get a highly nontrivial situation. Let L = −d 2 /dx 2 and where U is the dilation operator, Uf (x) = q 1/2 f (qx), 0 < q < 1. It can be checked directly coincides now with the retarded pantograph equation, which was investigated in detail in [28]. As follows from the results of this paper, in the retarded case the pantograph equation does not admit solutions belonging to L 2 (R). However, the operator B + has infinitely many normalizable eigenstates. The equation B + ψ + α (x) = αψ + α (x) coincides now with the advanced pantograph equation which has infinitely many solutions belonging to L 2 (R) [28]. For the free particle realization of the N > 1 symmetry algebras, one gets the generalized pantograph equations considered in [27]. As one of the open problems in this field, let us mention a need to characterize a minimal set of solutions of the equation (30) providing a complete basis of the Hilbert space L 2 (R).

Solitons, Ising chains, and random matrices
where τ N is a Wronskian of N different solutions of the free Schrödinger equation (18). Actually, there are many determinantal representation for this tau-function, e.g. [1]  A subclass of the self-similar potentials appearing from the q-periodic closure (27) can be viewed as a particular infinite soliton potential. Indeed, consider the KdV pM-soliton solution with the parameters k j subject to the constraint k j+M = qk j , 0 < q < 1, and take the p → ∞ limit. The limiting potential has the self-similar spectrum k pM +m = q p k m , m = 0, . . . , M − 1, and an infinite number of free parameters θ (0) j . The scaled potentialũ(x, t) = q 2 u(qx, q 3 t) has the same solitonic interpretation with j , i.e. impose the additional constraint θ j (qx, q 3 t) = θ j+M (x, t), which is seen to match with the condition (27). In this picture the discrete dilation x → qx, t → q 3 t just deletes M solitons corresponding to the smallest eigenvalues of L -this feature is the simplest physical characterization of the self-similar potentials.
Let us turn now to the statistical mechanics applications. N-soliton tau-function of the KdV equation (31) can be represented in the following Hirota form [1,24]: The coefficients A ij , describing phase shifts of the i-th and j-th solitons after their collision, have the form As noticed in [33], for θ i = θ (0) this tau-function coincides with the grand partition function of a lattice gas model on a line (for a two-dimensional Coulomb gas picture, see [34] and the next section). In this interpretation the discrete variables σ i describe filling factors of the lattice sites by molecules, θ (0) is a chemical potential, and A ij are the interaction energies of molecules [7].
It is well known that a simple change of variables in the lattice gas partition function leads to the partition function of some Ising chains: where N is the number of spins s i = ±1, J ij are the exchange constants, H i is an external magnetic field, β = 1/kT is the inverse temperature. Indeed, a replacement of filling factors in (32) by the spin variables σ i = (s i + 1)/2 yields where Similar relations with Ising chains are valid for the whole infinite Kadomtsev-Petviashvili (KP) hierarchy of equations and some other differential and finite difference nonlinear integrable evolution equations. Corresponding tau-functions have the same form (32) with a more complicated structure of the phase shifts A ij and phases θ i (a partial list of such equations can be found in [1,24]).
The general relation between N-soliton solutions of integrable hierarchies and partition functions of the lattice statistical mechanics models established in [33] brings a new point of view upon the self-similar potentials as well. Namely, as shown in [33] self-similar spectra can be deduced from the condition of translational invariance of the infinite spin chains.
E.g., let us demand that the spin chain is invariant with respect to the shift of the lattice by one site j → j + 1 which assumes that J i+1,j+1 = J ij . As a result, the exchange constants J ij (or the soliton phase shifts A ij ) depend only on the distance between the sites |i − j|.
This very simple and quite natural physical requirement forces k i to form one geometric progression (N.B.: uniquely) where k 0 and 0 < q < 1 are some free parameters. In a more complicated case one demands the translational invariance with respect to the shifts by M lattice sites, i.e. J i+M,j+M = J ij , and this results in the general self-similar spectra k j+M = qk j . The main drawback of these relations between the solitons and Ising chains is that for fixed q the temperature T (or β) is fixed too, which is clearly seen from the comparison of (36) and (37) (in the Coulomb gas interpretation one actually finds that β = 2).
For finite chains 0 ≤ j ≤ N − 1, the translational invariance is not exact. The infinite soliton limit N → ∞ corresponds to the thermodynamical limit in statistical mechanics. In this picture, all coordinates of integrable hierarchies (x, t and higher hierarchy "times") are interpreted as particular parameters of the external magnetic field H i . Since 0 < q < 1, the x and t-dependent part of H i decays exponentially fast for i → ∞. As a result, in the limit N → ∞ only the constants θ (0) i are relevant for the partition function (speaking more precisely, they determine the leading asymptotics term of the partition function at N → ∞).
Although in the thermodynamic limit the x, t-dependence is washed out and the Ising chain is periodic, presence of an infinite number of free parameters θ In order to compute the partition function one can use determinantal representations for the tau-function, e.g. the Wronskian formula (18). As shown in [33], = q θ 2 4 (ν/2, q 2 ) θ 2 1 (ν/2, q 2 ) (θ 1,4 are the standard Jacobi theta-functions). In the derivation of this formula the Ramanujan 1 ψ 1 sum was used at the boundary of convergence of the corresponding infinite bilateral basic hypergeometric series. A standard physicists' trick was used for dealing with that, namely, a small auxiliary parameter ǫ was introduced into the original expression for the density function ρ(ν) which guarantees the absolute convergence of the series for ǫ > 0. Although the limit ǫ → 0 leads to some singularity in ν on the right hand side of (39), it does not lead to divergences after taking the integral in (38) (a more precise mathematical justification of this physical "harmlessness" is desirable). The total magnetization of the lattice m(H) = −∂ H f I = lim N →∞ N −1 N −1 i=0 s i has the following appealing form In order to change the fixed value of the temperature, one may try to replace (37) by where n is some sequence of integers, and look for integrable equations admitting N-soliton solutions with such phase shifts. This is a highly nontrivial task, and in [33] only one more permitted value of the temperature was found. It corresponds to n = 2 and appears from special reduction of the N-soliton solution of the KP-equation of B-type (BKP) [13]. The general BKP τ -function generates much more complicated Ising chains than in the KdV case. The corresponding exchange constants have the form where a i , b i are some free parameters. Translational invariance of the spin lattice, J ij = J(i − j), results in the following spectral self-similarity where we set a 0 = 1 and assume that α > 0. This gives the exchange The parameter b is restricted to the unit disk |b| ≤ 1 due to the inversion symmetry b → b −1 .
For −1 < b < −q one gets now the ferromagnetic Ising chain, i.e. J ij < 0; for q < b ≤ 1 and b = e iφ = −1 the chain is antiferromagnetic, i.e. J ij > 0. In the thermodynamic limit N → ∞ the partition function Z N is given again by the determinant of a Toeplitz matrix, which is diagonalized by the discrete Fourier transformation. Similar picture holds for arbitrary M-periodic chains, see [33] for explicit expressions for the magnetization and technical details of computations.
It is well known that the n-particle lattice gas partition functions define probability distribution functions in the theory of random matrices. In a natural way, partition functions of the Ising models considered in the previous section may be related to the probability distribution functions of some random matrix models with a discretized set of eigenvalues.
Let us describe briefly this correspondence comparing the KP N-soliton solution and the Dyson's circular ensemble.
Dyson has introduced an ensemble of unitary random n×n matrices with the eigenvalues ǫ j = e iφ j , j = 1, . . . , n, such that after the integration over the auxiliary "angular variables" the probability distribution of the phases 0 ≤ φ j < 2π takes the form P dφ 1 . . . dφ n ∝ i<j |ǫ i − ǫ j | 2 dφ 1 . . . dφ n . One may relax some of the conditions used by Dyson and work with a more general set of distribution functions. In [19] Gaudin has proposed a circular ensemble with the probability distribution law containing a free continuous parameter ω. It interpolates between the Dyson's case (ω = 0) and the uniform distribution (ω = 1). Let us note in passing, that the paper [19] contains in it implicitly a specific model of the q-harmonic oscillator. The model (43) admits also an interpretation as a lattice gas on the circle with the partition function where β = 1/kT = 2 is fixed and the potential energy has the form A surprising fact is that the grand partition function of this model may be obtained in a special infinite soliton limit of the N-soliton tau-function of the KP hierarchy [33]. This means that the finite KP soliton solutions provide a discretization of the model, namely, they define a lattice gas on the circle. This leads to random matrices with a discrete set of eigenvalues. E.g., one may take unitary n × n matrices with eigenvalues equal to N-th roots of unity, i.e. ǫ j = exp(2πim j /N), m j = 0, . . . , N − 1. The probability measure is taken to be continuous in the auxiliary "angular" variables of the unitary matrices and discrete in the eigenvalue phase variables φ j . More precisely, the integrals over φ j are replaced by finite sums over m j and the continuous model is recovered for m j , N → ∞ with finite m j /N: The n-particle partition function becomes while the grand canonical ensemble partition function takes the form where η = ln(2π/N) is an excessive chemical potential, and To conclude, the grand partition function of the Gaudin's circular ensemble coincides with the particular infinite soliton KP tau-function at zero hierarchy "times". The root of unity discretization of circular ensembles has been considered by Gaudin himself [20], where a connection with Ising chains was noticed as well, but the relation with soliton solutions of integrable equations was not established. The BKP hierarchy of equations suggests a generalization of the distribution law (43) to and its discrete analogue, which were not investigated yet appropriately.
It is well known that the n-tuple Selberg integral provides an explicit evaluation of the n particle Coulomb gas partition functions for arbitrary values of the inverse temperature β (i.e. not just for β = 1, 2, 4, corresponding to orthogonal, unitary and symplectic ensembles).
It would be interesting to know whether a similar universal formula exists for the grand partition functions of lattice Coulomb gases. The energy of an electrostatic system ("plasma") consisting of N particles in a bounded domain of the plane is

Lattice Coulomb gas on the plane
where z j = x j + iy j and q j are the particles' coordinates and charges. The first term is the standard Coulomb energy, the second one describes an interaction with the boundaries (or the charge-image interaction), and the last term corresponds to the external fields contribution.
Suppose that our plasma is composed from the particles of equal charges, q j = +1 (in the two component case q j = ±1), upon a discrete lattice Γ on the plane. The grand partition function of such lattice Coulomb gas is where µ is an effective chemical potential. In general case one can rewrite G N in the form [7] where σ(z) = 0 or 1 depending on whether the site with the coordinate z is empty or occupied by a particle, and the functions q(z) and µ(z) characterize distribution of particles of different types. E.g., in the two-component case, when q(z ± ) = ±1 charges occupy the {z ± } sublattices, one has µ(z ± ) = µ ± . Note that these 2D lattice gases describe simultaneously some 2D Ising magnets.
Let us take now the N-soliton tau-function (32) and replace in it the soliton number j by a variable z taking N discrete values. As a result, it may be rewritten in the form But this is precisely (50) -one just needs to identify the phase shifts A zz ′ with the Coulomb interaction potential W (z, z ′ ), and the phases θ(z) with the function w(z). As a result of this observation made in [34], one can construct a number of new solvable models of Coulomb gases in addition to already known ones (see, e.g. [17,19,20] and references therein).
For the KP-hierarchy soliton solutions one has [13] A zz ′ = ln where t p is the p-th KP hierarchy "time" and a z , b z are some arbitrary functions of z. For the BKP-hierarchy A zz ′ are given by (41) and where t 2p−1 are the BKP evolution "times".
Let us take a z = z = x + iy, b z = −z * = −x + iy. Then in the KP case where V (z, z ′ ) is the potential created by a positive unit charge placed at z ′ over a conductor with its surface occupying the y = 0 line. In this case V (z, z ′ ) solves the Poisson equation with the boundary condition E t (y = 0) = 0. The same potential is created by a positive charge placed at the point z ′ and its image of opposite charge located at the point (z ′ ) * .
Similar to the random matrix models situation, the correspondence between 2D lattice Coulomb gases with specific boundary conditions and solitons is valid only for fixed β, which is found from the comparison of (53) with (50) to be β = 2. Since w(z) = θ(z), one finds an explicit expression for the zero time soliton phases θ (0) (z) in (52): where the middle term corresponds to the charge-image interaction, and φ(z) is the potential created by a neutralizing background of some density ρ(r), ∆φ(r) = −2πρ(r), φ x (y = 0) = 0.
The harmonic term ∞ p=0 (z p −(z * ) p )t p = −βφ ext (z) in (52) corresponds to an external electric field. One may conclude that the imaginary hierarchy times evolution of the KP soliton solutions describes the behavior of a 2D lattice Coulomb gas in a varying external electric field.
Usefulness of the conformal transformations z → f (z) is well known in the potential theory. With their help one can map plasma particles to various regions. So, the map z → z n puts them inside a corner with the π/n angle between the conducting walls. The exponential map z → e πz/L leads to a gas inside the strip 0 < ℑz < L between two parallel conductors, etc.
As to the BKP equation case, the choice a z = z, b z = z * leads to the Coulomb gas placed inside the upper right quarter of the plane with the dielectric and conductor walls along the x and y-axes respectively. The discrete temperature appears to be the same, β = 2. A curious set of solvable dipole gas models is generated by the self-similarity constraints imposed upon the spectral data of the corresponding multi-soliton systems. For further details, see [34].

Appendix. A heuristic guide to special functions of one variable
There are many handbooks and textbooks on special functions. However, none of them contains a formal list of properties which a function should have in order to be "special".
One usually talks on functions of some particular type (hypergeometric, L-functions, etc).
R. Askey has suggested one universal definition that if a function is so useful that it starts to bear some name, then it is "special". Another essential, but not so universal, characteristic refers to the asymptotic behavior. Namely, for special functions one is expected to be able to deduce asymptotics at infinity from a known local behavior, i.e. the connection problems should be solvable -such an approach is advocated by people working on the Painlevé type functions.
To the author's taste these two definitions are relying upon the secondary features of special functions. One has to have already the functions in hands in order to start to investigate their properties. If one takes as a goal the search of new special functions, then it is necessary to find a definition containing a more extended set of technical tools for work.
In this respect it should be stressed that even the term "classical special functions" appears to be not so stable. E.g., it is by now accepted that the family of classical orthogonal polynomials consists of not just the Jacobi polynomials and their descendants, but of the essentially more complicated Askey-Wilson polynomials invented just two decades ago.
The group theory provides a number of tools for building new functions, but, unfortunately, connections with the representation theory often provide only interpretations for functions already defined by other means. Still, the symmetry approach seems to be the central one in the theory of special functions. In particular, the key "old" special functions appear from separations of variables in the very simple (and, so, useful and universal) equations [38]. For the last decade the author's research was tied to the following working definition: special functions are the functions associated with self-similar reductions of spec-tral transformation chains for linear spectral problems. Speaking differently, special functions are connected with fixed points of various continuous and discrete symmetry transformations for a taken class of eigenvalue problem. This definition works well only for special functions of one independent variable (which may, however, depend on an infinite number of parameters) and even for them it does not pretend to cover all possible special functions. On the one hand, this definition comes from the theory of completely integrable systems, where a search of self-similar solutions of nonlinear evolution equations is a standard problem [1].
On the other hand, particular examples found from this approach show that it has in its heart the contiguous relations -linear or nonlinear equations connecting special functions at different values of their parameters.
Schematically, a heuristic algorithm of looking for these "spectral" special functions consists of the following steps: 1. Take a linear eigenvalue problem determined by differential or difference equations. with respect to taken symmetry transformations. As a result, some finite sets of nonlinear differential, differential-difference, difference-difference, etc equations are emerging which de-fine "nonlinear" special functions. Solutions of the original linear equations with coefficients determined by these self-similar functions define "linear" special functions.
The last two steps are heuristic since, despite of a big progress in the general theory of symmetry reductions (see, e.g. [31,36]), no completely regular method has been built yet. E.g., the reductions used in the derivation of the associated Askey-Wilson polynomials' recurrence coefficients [49] and in the discovery of new explicit systems of biorthogonal rational functions described in Zhedanov's lecture at this meeting did not find yet a purely group-theoretical setting.
Another important ingredient of the theory of special functions, which was not listed in this scheme, is the theory of transcendency. Painlevé functions are known to be transcendental over the differential fields built from a finite tower of Picard-Vessiot extensions of the field of rational functions. In solving differential (difference, or whatever) equations one has to determine eventually which differential field a taken solution belongs to (e.g., whether it belongs to the field of functions over which the differential equation is defined).
As an open problem in this field, which was discussed partially in [47], we mention a need to find a differential (or difference) Galois theory interpretation of the self-similar reductions of factorization chains.