Generalized Devil’s staircase and RG flows

,

The Devil's staircase is the fingerprint of the "incommensurability phenomena" in a variety of physical systems [1,2].The geometric signature of the incommensurability is the Riemann-Thomae (RT) function which often emerges in spectra of sparse systems of various physical origin.Meanwhile, the Riemann-Thomae function also appears in a plethora of fundamental problems, such as stability diagram in fractional quantum Hall effect [3,4], interactions of non-relativistic ideal anyons with rational statistics in the "magnetic gauge" approach [5], quantum 1/f noise and Frenel-Landau shift [6], distribution of quotients of reads in DNA sequencing experiment [7], frequency of specific subgraphs counting in the protein-protein network of a Drosophilla [8].Though the degree of similarity with the original RT function could vary, and experimental profiles may drastically depend on the peculiarities of each particular physical system, a general probabilistic scheme resulting in emergence of the fractal hierarchical distribution can be considered as the manifestation of number-theoretic laws in nature.
One possible pattern behind the Riemann-Thomae function and the Devil's staircase is as follows.Consider a physical problem, for example the fractional quantum Hall effect (FQHE), and push the system into the particular limit in the parameter space.For FQHE this is the so-called "thin torus limit" -see for example [4].The system hosts some defects, and in the limit under consideration defects form a lattice which is a Wigner crystal in the thin torus limit of FQHE.Consider now the propagation of a probe particle through the sample which can be studied, for instance, by analysing the spectral density.The modular SL(2, Z) group acts in the parameter space of this system.The imaginary part Im τ of the modular parameter τ gets identified with some function of disorder, while the real part Re τ corresponds to the chemical potential for the topological charge relevant for the studied problem.The motion of the probe particle in the crystal of defects can be mapped onto the motion in the fundamental domain of SL(2, Z), and the rearrangements of the lattice can be treated by analyzing the RG flow in the vicinity of transition points which are identified with points of lattice bifurcations.Generally speaking, from the probe particle perspective, the rearrangement of the lattice can be studied by varying the chemical potential of defects (or of their number).
Another view on the Devil's staircase deals with a general classification of quantum systems spectra which usually involve discrete, continuous, as well as more tricky singularcontinuous supports.Recent discussion of the latter case can be found in [9] and it incorporates the Devil's staircase as an intrinsic ingredient.The mechanism behind the appearance of the Devil's staircase in the spectrum can be illustrated in the context of the Peierls model [10] of electrons interacting with the lattice of ions.Starting with the integrable version of the Peierls model described by the Toda lattice, authors of [11] have shown that the spectral curve of the Toda system can be identified with the dispersion law of fermions.At particular values of Toda integrals of motion there very tiny bands in the fermion spectrum emerge and if one adds an arbitrary small non-integrable perturbation, the Devil's staircase structure gets formed nearby these bands at rational fermionic densities indicating the emergence of incommensurability transitions.In our work the similar viewpoint is used for studying the spectral statistics of a probe particle propagating along a weakly disordered crystal, which we interpret as a two-parametric renormalization group (RG) flows in a particular limit (as it is explained below).
The general classification of RG flows rhymes with the development of bifurcations ("catastrophes") over time in the theory of dynamical systems -see, for instance [12].In the catastrophe theory there are focuses, saddles, limits cycles and other attributes of the singularity theory, with corresponding fixed points, RG cycles and more exotic RG behavior.For instance, recently the RG counterparts of homoclinic orbits in the theory of dynamical systems have been found in the field theory [13], they also provide examples of chaotic RG flows [14].The incommensurability phenomenon is also known in the theory of dynamical systems.Hence, following the same logic, one could expect the existence of RG counterpart of the incommensurability.Indeed, the RG approach was successful in describing the Devil's staircase pattern in a Harper equation for the electron in a crystal in presence of a magnetic field [15,16] where it was argued that the tunneling in the phase space is the crucial ingredient.
In many situations it is convenient to combine two real parameters of a 2D RG flow into the single complex parameter, τ , which can be interpreted as the modulus of the complex structure for an auxiliary elliptic curve.The familiar examples are: the Anderson localization problem with the time symmetry breaking (TSB) term [17], the integer quantum Hall effect (IQHE) [18,19], and the Yang-Mills theory with the TSB θ-term [20,21].In all these examples the real part of the complex parameter is the TSB parameter.We suggest a bit more general perspective and propose to consider the following generic complex (modular) parameter: hence the RG flow unites the topology and the disorder.Let us provide some known examples supporting this perspective: (i) In the integer quantum Hall effect (IQHE) the complex parameter τ is built of two conductivities, σ xx and σ xy [18]: (ii) In gauge (Yang-Mills) theories the modular parameter τ involves the coupling constant, g Y M , and the θ-term in the following combination [20,21]: (iii) In the case of three-diagonal random matrices with the off-diagonal disorder the parameter τ enters in the combination (see [22,23] and Section III for detail): where ϵ is the function of the spectral parameter, λ, and f depends on the strength of the disorder; (iv) In the Anderson model with the TSB term the parameter τ reads [17]: where θ counts windings and D is the diffusion coefficient; (v) For the polymer propagating in the lattice of obstacles, the entanglement complexity is fully characterized by the parameter τ known as a "primitive path configuration" (see for detail [24][25][26]: τ = (chemical potential of winding) + i (complexity of entanglement).(6) where the "primitive path" was defined in [24,27] and has the meaning of a geodesic path on some hyperbolic manifold [26,28].
At any fixed value τ = χ + iξ the partition functions of considered systems fully enjoy symmetries of the SL(2, Z) modular group and hence are the modular functions.However when χ and/or ξ run over time and depend on a scale, µ, the situation is more subtle.It general, the RG flow involves two β-functions and is described by the set of equations Typically, the disorder term enjoys both the perturbative and non-perturbative remormalizations, while the topological parameter is renormalized only non-perturbatively.
There are some known patterns of β-functions with familiar properties: • For the Integer Quantum Hall Effect (IQHE) case the β-function can be expressed via the Grassmanian U (n+m) U (n)×U (m) σ-model: where D m,n and C m,n are not completely universal.
• For the Russian Doll model which is the toy example of the system with the cyclic RG flows (see, [29] for review), the RG flow is discrete • For the Berezinskii-Kosterlitz-Thouless (BKT) transition one has: We focus our attention on a specific limit of RG flows when the non-perturbative renormalization coming from instanton-like contributions dominates -see, for example, [16].This happens in all examples when ξ = Im τ → 0, which means that we are looking at the limit of a weak disorder in some frame, and the modular parameter is mainly governed by the "θ" (i.e.winding-like) terms.The details are model-dependent, however in all cases the θ-term has one and the same physical sense: it serves for counting topological defects.The weak disorder limit (ξ → 0) in some cases could mean the strong coupling.For example, in the Yang-Mills theory, since ξ = 4π , in the limit ξ → 0 the theory is strongly coupled in the "magnetic frame", while in the "electric frame" the theory is weakly coupled.We are searching for some universality in the Im τ → 0 regime.Throughout the work we argue that gRT function and generalized Devil's staircase emerge naturally in the ξ → 0 limit and are universal.
In all cases (i)-(v) mentioned above, the β-functions can be expressed in terms of elliptic functions on some Riemann surface.However, when couplings χ and ξ run in time, the construction of corresponding Riemann surfaces is a nontrivial issue.The Riemann surface is bundled over some manifold and the RG flow is identified with the dependence of the modular parameter on the point of the fibration base.The benchmark example is provided by the Seiberg-Witten solution for the low-energy effective action of the N = 2 SYM theory [30].The renormalization of both (Re τ, Im τ ) can be found directly from the partition function of ensemble of instantons which can be evaluated via the localization approach [31].
Remarkably, two more ways to handle RG flows are available.The first one can be formulated quite generally.We add the probe object into the ensemble of defects, study its dynamics and derive from the induced dynamics the behavior of β-functions.This logic works well in the Seiberg-Witten solution when the surface defect inserted into the instanton ensemble plays the role of a probe.The prepotential of the low-energy N = 2 super Yang-Mills (SYM) theory, F, yields the exact β-function [30].This β-function can be derived from the semiclassical wave function of the probe whose dynamics is governed by the integrable system of Calogero-Toda type [32][33][34].The derivation of the potential governing the dynamics of the probe from the first principle is not a simple issue and only recently this problem has been fully elucidated [35].In our study we use a similar probe analysis.
The second way to treat RG flows is more sophisticated and involves the nontrivial dynamics rooted in the so-called "vertex realization" of the infinite-dimensional algebras [36].In brief, the partition function of the instanton ensemble gets mapped onto the particular conformal block in the Liouville theory, or its extension to higher spins (the so-called "W n theory").The modular parameter becomes the position of the vertex operator insertion on the sphere.Hence the flow in the modular domain in the asymptotic regime gets reduced to the investigation of a particular asymptotic regime of a conformal block.
The paper is organized as follows.In Section II we recall the main facts concerning the Riemann-Thomae (RT) function and provide its analytic regularization in terms of the Dedekind η-function.In Section III we consider the non-perturbative problem of a probe particle propagating in a weakly disordered 1D lattice and pay attention to the spectral density of the probe whose Hamiltonian is given by the tridiagonal matrix with the offdiagonal bimodal disorder.In the same section we formulate the version of the RG flow in the fundamental domain of the modular group yielding the RT function.In Section IV we consider the propagation of the probe in 2D lattice of defects and investigate the dependence of probe dynamics on the lattice pattern.We show that a lattice with a specific asymmetry yields the Silver ratio as the asymptotic value of the flow in the fundamental domain of the group SL(2, Z) in the weak disorder limit.In the same Section we discuss the RG flows in vicinity of bifurcation (transition) points and conjecture the BKT-like divergence of the correlation length.In Section V we discuss the spin system with long-range interactions exhibiting the Devil's staircase behavior and focus on its dependence on the form of the longrange potential.We find that gRT emerges in considered systems and exponents involved in gRT function provide a kind of universality.Using duality properties between the integrable systems we comment on the similarity between the incommensurability observed in the longrange Hubbard model and in FQHE in the thin torus limit.In Discussion we overview related issues which require more profound analysis.In Conclusion we summarize our findings and formulate open questions.

II. RIEMANN-THOMAE FUNCTION, EUCLID ORCHARD AND DEVIL'S STAIRCASE
A. Reminder on the standard Riemann-Thomae function The Riemann function, g(x), also known as the Thomae function [37], has many other names: the popcorn function, the raindrop function, the countable cloud function, the ruler function, the modified Dirichlet function.It is one of the simplest number-theoretic functions possessing a nontrivial fractal structure (another famous example is the everywhere continuous but nowhere differentiable Weierstrass function).The Riemann-Thomae (RT) function is defined in the open interval x ∈ (0, 1) according to the following rule: The function g is discontinuous at every rational point: irrationals, where g vanishes, come infinitely close to any rational number.At the same time, g is continuous at irrationalssee Fig. 1a.To provide basic properties of the RT function, consider some irrational number, t, at which g(t) = 0, and take some 0 < ε < 1.Without the loss of generality, ε is assumed to be rational, otherwise one may replace ε with any smaller rational ε ′ = k s < ε such that gcd(k, s) = 1.Thus, there is a finite set of rational numbers Ω ε = {n = i j , 1 < j ≤ s, 1 ≤ i < j}, that are not smaller than g(ε).Now assign δ(ε) = inf{|n − t|, n ∈ Ω ε } that defines the vicinity of t, in which the values of g are smaller than ε: |t − y| < δ(ε) → |g(t) − g(y)| = |g(y)| < ε.These inequalities prove the continuity of the Riemann function g(x) at any irrational value of x.
An elegant representation of the Riemann-Thomae function emerges in a so-called "Euclid orchard" construction -see Fig. 1b [38].Consider an orchard of trees of unit heights located at every point (am, ak) of the two-dimensional square lattice, where m and k are nonnegative integers defining the lattice, and a is the lattice spacing, which is convenient to choose as a = 1 √ 2 .Suppose that the observer stays on the line m = 1 − k between the points A(0, a) and B(a, 0), and watches the trees in the first quadrant along the rays emitted from the origin O(0, 0).Along these rays, we see only the first (non-shadowed) tree with coprime coordinates, M (am, ak), while all other trees are shadowed.Introduce the rotated coordinate system (x, y) with the axis 0x along the segment AB and the axis 0y normal to the orchard's plane, as shown in Fig. 1a.We set the origin of the 0x axis at the point A, then the point B has the coordinate x = a.Having the focus located at the origin, the tree M N at the point M (am, ak) is projected to the tree M ′ N ′ located at M ′ (x, y = 0), where x = m k+m and the height |M ′ N ′ | of this tree is 1 k+m -see Fig. 1b.Denoting k + m by n, we immediately conclude that the "visibility diagram" in the Euclid orchard is exactly the Riemann-Thomae function for the variable x = m n .
The emergence of the distribution (12) can be understood on the basis of the generalized Euclid orchard construction if one considers a (1 + 1)-dimensional directed walk on the lattice starting at the origin and making ϕ steps along one axis, followed by ψ steps along the other axis.At every lattice site, the walk survives with the probability p and dies with the probability q = 1 − p.Having an ensemble of such walks, one arrives at the model of "hooked walks" in the Euclid orchard.Thus, for some point, ν, and at a small death probability, q, (q → 0), a fraction of survived walkers, P (ν), computed in ( 12) is described by the Riemann function (11).

B. Relation of the Riemann-Thomae function to Eisenstein series and Dedekind η-function
Define the generalized Riemann-Thomae (gRT) function, g(x), as follows: where h(n) = n −α (α > 0).The sample plots of g(x) for two arbitrary chosen values, α = 0.41 and α = 2.76 (for n = 100) are shown in Fig. 2a,b.A natural, physically justified analytic regularization of the Riemann-Thomae function is highly demanded.Below, we briefly describe such a regularization for gRT function g(x) where h(n) = n −2 considered firstly in [23] and extend our construction to any f (n) = n −α (α > 0) and even to non-algebraically decaying potentials.Namely, we demonstrate that the analytic approximation of the function g(x) for h(n) = n −2 involves the Dedekind η-function, η(x + iy), defined in the halfplane y > 0.
Let us rewrite g(x) for h(n) = n −2 on the interval x ∈ (0, 1) as follows: The function g(x) assigns zero to all irrational points.Now, using the identity: the function g(x) in ( 14) can be regularized by 0 < ε ≪ 1 as follows: where The latter term in ( 17) is proportional to δ(x) and can be neglected.The coefficient C(ε) in ( 16) can be computed from the normalization condition and it is instructive to fix it at the end of derivation.Since n and m in (17) run over all integer points except 0, it is convenient to change a sign in front of m: m → −m and rewrite (17) as follows Recall now the definition of the non-holomorphic Eisenstein series, E(z, s), [39]: where E(z, s) is a function of z = x + iy and is defined in the upper half-plane y > 0 for all Re (s) > 1. Comparing ( 18) and ( 19), we can straightforwardly conclude that the function ḡ(x, ε) matches the Eisenstein series E(z, s) at s = 1 upon the identification ε = y, i.e.
x + iε = z.Thus, The non-holomorphic Eisenstein series of weight 0 and level 1 can be analytically continued to the whole complex s-plane with one simple pole at s = 1.Notably E(z, s), as function of z, is the SL(2, Z)-automorphic solution of the hyperbolic Laplace equation: The residue of E(z, s) at s = 1 is known as the first Kronecker limit formula [40][41][42].Explicitly, it reads at s → 1: where γ is the Euler constant and η(z) is the Dedekind η-function.Equation (22) establishes the important connection between the Eisenstein series and the Dedekind η-function, which we exploit below.The Dedekind η-function is defined as follows: The argument z = x + iy is called the modular parameter, and η(z) is defined for all y > 0.
The function η(z) is invariant with respect to the action of the modular group SL(2, Z): In general, where ad − bc = 1 and ω(a, b, c, d) is a 24th degree root of unity [43].
Collecting together (20) and (22), taking into account that ε = y, and omitting the divergent constant at s → 1, we get the following asymptotic analytic expression for the Riemann-Thomae function g(x) for h(n) = n −2 : where The coefficient C(y) in (26) we compute from the condition (see the Appendix A for more detail).Specifically, the value of C(y) up to logarithmic corrections reads In Fig. 3a we have depicted the function g(x) for h(n) = n −2 with the maximal denominator n max = 10 2 together with the function − 12y π ln f (x, y) at the fixed value y = 5 × 10 −4 , where f (x, y) is defined in ( 27)- (28), while in Fig. 3b we have plotted the cumulative (integrated) function G(x), where The function G(x) has a typical "Devil's staircase" shape [44] reflected horizontally and rotated by π/2.Later on we shall consider the "generalized Devil's staircase", G α (x) defined as in (29) for h(n) = n −α (see (13)) at any values of α > 0 (G(x) ≡ G 2 (x)).

III. RIEMANN-THOMAE FUNCTION AND THE SPECTRAL DENSITY OF RANDOM TRIDIAGONAL OPERATORS
Consider an ensemble of random operators represented by tridiagonal N × N (N ≫ 1) random symmetric matrices A N with the bimodal (Bernoulli) distribution of sub-diagonal matrix elements: We are interested in spectral properties of an ensemble of such matrices.Namely, we compute the density of eigenvalues, ρ(λ), in the limit N → ∞ and demonstrate its connection to the Riemann-Thomae function g(x) defined in (11).For the first time this question was addressed in [22] and below we present slightly more extended version of our construction.
To proceed, note that at any x k = 0, one can split the matrix A N into independent blocks along the diagonal.So, it is instructive to consider subsequences of gapless sets of with x k = 1.The set of eigenvalues of a symmetric gapless n × n three-diagonal block A n with x k = 1 for all k = 1, ..., n, is The probability of having a gapless subsequence with n consecutive "1" is Q n = p n , since all x k are independently distributed Bernoulli variables.The sample plots ρ(λ) for two different values of p, namely for p = 0.9 and p = 0.5 computed numerically for N = 500 over 500 different matrix realizations, are shown in Fig. 4.
The spectral density ρ(λ) can be written in the limit N → ∞ in terms of the Riemann-Thomae function g(τ ) = g(x + iy).The resulting expression reads:  where The parameter y in (32) has a sense of a "resolution cutoff" of the Dedekind relief -see Fig. 9 in Section IIA.The relation between the strength of the disorder, p, and the cutoff, y, can be established using the following qualitative arguments.On one hand, the maximal denominator, n max , of the Thomae function (see ( 13)) defines the total number of peaks that can be resolved upon n max ln p ∼ 1, as it follows from (32).On the other hand, the cutoff y can be estimated as y ∼ 1/n max .Thus, in the limit p → 1 one has y ≈ 1 − p.In next Section we provide the outline of the derivation of (32).Let us emphasize that throughout the paper we pay attention to the non-perturbative limit p → 1 of a weak disorder.

A. Analogy with the Peierls model
Let us discuss how this toy 1D example matches the general perspective on the Devil's staircase.To this aim, we compare our example with the Peierls model of 1D superconductivity where the Devil's staircase has been recognized as well [10].One starts with the integrable version of the Peierls model [11] in which the Hamiltonian can be written as follows The Hamiltonian (34) describes the interactions between phonon-like degrees of freedom, θ n , of fluctuations around the regular lattice x n = na + θ n , n = 1 . . .N and fermions propagating at the top of the lattice.The Hamiltonian of phonons in the integrable case involves a few lowest Toda chain Hamiltonians, H k = Tr L k f , where the N × N Lax operator for the Toda chain plays the role of the Hamiltonian for the fermions.The quasimomentum, η, corresponds to the periodic lattice, p k are the momenta of phonons and c k = exp(θ k+1 − θ k ).The Lax representation for the Toda chain provides the consistency of the phonon and fermion dynamics.The solution at the generic values of the Toda Hamiltonians and some value of the fermionic density, ρ = q N , is expressed in terms of the hyperelliptic Riemann surface whose moduli are defined by the Hamiltonians and fermionic density.The Riemann surface simultaneously plays the role of the dispersion law for fermions.
The transition from commensurability to incommensurability in the fermionic spectrum occurs if two conditions are fulfilled simultaneously: • Select the values of Toda Hamiltonians in such way that some bands in the spectrum became very tiny; • Select the rational value of the fermionic density and add a very weak non-integrable deformation of the Hamiltonian.
Upon these two conditions the Devil's staircase emerges [10].It can also be reformulated in terms of the Whitham dynamics for the Toda chain describing the flow of the integrable systems solutions in the moduli space.Whitham dynamics was interpreted as a version of a RG flow for many systems -see, for instance, discussion in [45].
Let us compare the Peierls model setup and the toy model described in the previous Section.We have the 1D lattice of N sites and the particle hopping on this lattice, which can be identified with the particle Hamiltonian (30).The two-step procedure similar to the one described above for the Peirles system is as follows: • We see that in our random operator (30) all diagonal matrix elements are equal to zero.In terms of the Peierls model it means that p i = 0 and the lattice is frozen.
Freezing the lattice indeed means a deep degeneration of the Riemann surface; • Instead of adding a weak non-integrable deformation to the Peierls Hamiltonian, we introduce the weak randomness in the lattice which does the same job.
Hence, the emergence of the Devil's staircase in the spectrum of a tridiagonal matrix and in the Peierls model is of similar origin.The probe particle in our model plays the same role as Lax fermions in the Peierls model for the effective Whitham RG dynamics.The key point ensuring the full matching of the Peierls model with our scheme is the formulation of the RG procedure that can be translated into the flow in the fundamental domain of SL(2, Z).This will be explained in detail below.

B. Riemann-Thomae function and hyperbolic geometry
Let us remind a textbook definition: sequences of coprime fractions constructed via the ⊕ addition constitute the "Farey sequence" [46] A simple geometric model behind the Farey sequence, known as the Ford circles [47], is shown in The connection of the Ford construction with the spectrum of the ensemble of random operators (30) goes as follows.Write the position of the central peak at λ 1 = 0 as λ 1 = 0 = −2 cos 1 2 π .Write the spectral edge at λ ∞ = 2 as λ ∞ = −2 cos 1 1 π .Now consider the "enveloping" sequence of monotonically decreasing peaks which is denoted in Fig. 4 as S 1 .One can show that positions of eigenvalues λ 2 , λ 3 , ... shown in red in Fig. 5b are determined similarly to (37), namely Positions of resonances in other monotonic sequences, for example in the sequence S 2 in Fig. 4, one can again find recursively using the Farey construction (36).Corresponding eigenvalues for the sequence S 2 are denoted as λ 2 , λ ′ 1 , λ ′ 2 , ... and they are shown in blue in etc.The same Farey sequences can be sequentially generated by fractional-linear transformations (reflections with respect to the arcs) of the fundamental domain of the modular group SL(2, Z) -the triangle lying in the upper halfplane Im z > 0 of the complex plane z as shown in Fig. 5c.
C. Spectral density ρ(λ) and the "visibility diagram" The spectral density ρ(λ) of the ensemble of N ×N random matrices A N with the bimodal Bernullian distribution of matrix elements can be written in a form of a resolvent: where ⟨...⟩ means averaging over the distribution Q n = p n , and the identity (15) has been used to regularize the δ-function.Substituting ( 31) into (40), we find the following expression for ρ(λ): The sum in (41) looks rather complicated, however by exchanging the orders of summations (first running the summation in n and then -in k) one can advance in understanding the structure of ( 41) and its relation to the Riemann-Thomae function g(x) as stated in (32).
To proceed, note that the function ρ(λ) is nonzero only at λ k,n .In Fig. 6 we have plotted eigenvalues λ k,n = −2 cos πk n+1 as a function of k for a set of 20 fixed values of n: n = 1, ..., 20.Every point in Fig. 5 designates some eigenvalue λ k,n ; points along each solid gray curve correspond to different values of k for one and the same value of n.Each horizontal dashed line represents the set of the same eigenvalues λ coming from different n.Let us exchange orders of summation in k and in n in (41) and first sum weights of points along each horizontal dashed line over all n.In such a way we account of the degeneracy of a corresponding eigenvalue λ k and, respectively, the height of a peak in the spectral density ρ(λ) at λ k .Thus, by summing corresponding powers of p we get the height of any peak (i.e.eigenvalue).To compute ρ(λ) we split the spectrum into monotonic sequences of peaks as it has been discussed in the previous section.Sequences S 1 and S 2 shown in Fig. 5b are typical representatives.For them we have: Equations ( 42)-( 43) are tightly linked to the so-called visibility diagram (our notations are slightly different with the definition of the visibility diagram defined in [48]).Consider the square lattice of integer points (m, n) and add a weight p n−1 to each vertical row.Emit rays from the point (0, 0) at rational tangents of angles α m,n = arctan m n (m and n are coprimes, 1 ≤ m ≤ n − 1) within the wedge [π/4, π/2] in the positive direction as shown in Fig. 7 and sum up the weights, p n−1 , of all integer points along each emitted ray.
Let us show that the visibility diagram provides a straightforward way of calculating the spectral density ρ(λ).To see this, consider again the monotonic sequences of peaks S 1 and S 2 shown in Fig. 4 and in Fig. 6.For better visualization we keep colors throughout the 0 p 0 p 1 p 2 p 3 p 4 p 5 p 6 p 7 p 8 p 9 p 10 p 11 p 12 p 13 p 14 p text: the sequence S 1 is shown in red and S 2 -in blue.Summing rational points on the visibility diagram in Fig. 7 along a ray with the angle α m,n = arctan m n we get exactly the same results as given in ( 42)- (43).We demonstrate that on two examples: • Pick up the ray with tan α 2,3 = 2 3 , which corresponds to λ = −2 cos 2 3 π .Summing rational points (weighted with the corresponding power of p) along this ray we get the height: p 2 + p 5 + p 8 + p 11 + ... = ∞ s=1 p 3s−1 = p 2 1−p 3 , which coincides with the value for λ 2 (at n = 2) in (42).
The analysis of the visibility diagram allows to formulate the following general prescription for the evaluation of the spectral density ρ(λ) where n is the denominator of the fraction m n (m and n are coprimes).The angle α m,n = arctan m n uniquely defines the eigenvalue Comparing (44) with the definition of the Riemann-Thomae function g(x) in (11) we can immediately identify n with the denominator of the corresponding rational fraction x = m n .So, we have n = 1/g(x + iy) where the cutoff y defines the maximal denominator, n max .At irrational x the function g(x) is 0, which provides n = ∞, thus giving ρ(n) = 0 according to (44).From (45) we have λ = −2 cos(πx).Inverting this expression and taking into account the symmetry of the spectrum, we get x = 1 π arccos λ 2 which together with (44) leads to (32).

IV. RIEMANN-THOMAE FUNCTION AND PHYLLOTAXIS
Amazing connection of cell packing with Fibonacci sequences, known as phyllotaxis [49] was observed a long time ago in the works of naturalists and remains till now one of the most known manifestations of number theory in natural science.The generic description of growing plants based on symmetry arguments allowed researchers to uncover the role of Farey sequences in the plant's structure formation (see, for example, [50][51][52]), however, the question why the nature selects the Fibonacci sequence, among other possible Farey ones, was hidden until modern time.A tantalizing answer to this question has been given by L. Levitov in 1990 in [53], who proposed an "energetic" approach to the phyllotaxis, suggesting that the development of a plant is connected with an effective motion along the optimal path on a Riemann surface associated with the energetic relief of growing tissue.
The energetic mechanism suggested in [54] was applied later in [55] to the investigation of the geometry of flux lattices pinned by layered superconductors.It has been shown that under the variation of a magnetic field, the structure of the flux lattice can undergo a sequence of rearrangements encoded by the Farey numbers.However, lattices emerging in sequential rearrangements are characterized by the specific subsequence of the Farey set, namely, by the Fibonacci numbers.Very illuminating experiments have been provided in [56] for lattice formed by drops in rotating liquid, and in [57] for the equilibrium structure of a "magnetic cactus".
Here we consider, following L. Levitov, the model system of N strongly repulsive particles disposed and equilibrated on the surface of a cylinder of fixed diameter, D, and height, H and look at the rearrangement of these particles when the cylinder is compressed along its height under the condition that N and D remain unchanged -see Fig. 8a.At the continuous compression, for each height, particles form a triangular "Abrikosov" lattice with minimal energy [58].Various lattice topologies parametrized by the modular parameter, τ = D + iH, represent the valleys separated by energetic barriers on the manifold Γ with the non-archimedean ultrametric structure [59].
The notion of ultrametricity deals with the concept of hierarchical organization of energy landscapes [60].A complex system is assumed to have a large number of metastable states corresponding to local minima in a complex potential energy landscape.These minima are clustered in hierarchically nested basins: larger basins consist of smaller basins, each of these consists of even smaller ones, etc.The basins of local energy minima are separated by a hierarchically arranged set of barriers: large basins are separated by high barriers, and smaller basins within each larger one are separated by lower barriers.Ultrametric geometry fixes taxonomic (i.e.hierarchical) tree-like relationships between elements and, speaking figuratively, is much closer to Lobachevsky geometry, rather to the Euclidean one.
We provide an explicit construction of the energetic relief in a phase space of all possible patterns of compressed lattices for symmetric and asymmetric interaction potentials and demonstrate that the ground state is related to the deepest valley in Γ constructed via the Dedekind η-function.The lattice rearrangement caused by the compression of the cylinder along its axis is associated with the adiabatic flow along the geodesic in the energetic relief Γ, which can be understood as an RG flow.At each height, H, particles on the cylindric surface form a lattice with a minimal energy.Different lattice topologies, parameterized by the modular parameter τ = D + iH, have valleys separated by barriers on the manifold Γ with the non-archimedean ultrametric structure.For strongly compressed lattices (y ≪ 1) the energy U (x, y) as a function of x has a peak at every rational point, x = m n as it is shown in Fig. 8b.One sees that U (x, y) shares the hierarchical (ultrametric) behavior, which should be understood as follows: the transitions between two arbitrary local minima at x 1 and x 2 , are determined by the passage over the highest barrier U max (x 1 , x 2 ), separating the points x 1 and x 2 .
It is instructive to remind the definition of the ultrametric space.The ultrametric space, M is a set of elements supplied with the metric, i.e. the pairwise distance, d(x 1 , x 2 ) between elements x 1 and x 2 which meets three requirements: The ultrametric organization of the energy relief means that we identify the energy with the metric.Namely considering U (x, y = const) as a function of x, we may set d(x 1 , x 2 ) = U max (x 1 , x 2 ) as it is depicted in Fig. 8b.
This Section is organized as follows.We begin with the derivation of the symmetric potential U (x, y) separating valleys between different equilibrium configurations of particles on the cylinder when the cylinder is compressed along its height, H, under the condition that N and D remain unchanged -see Fig. 8.The corresponding analytic expression for potential barriers separating the valleys matches the generalized Riemann-Thomae function g(x) defined in (13) for h(n) = n −2 .We consider the RG flow of the minimum of the potential U (x, y) when y is tending to 0 (i.e. the lattice is strongly compressed) and propose the topological interpretation of the corresponding flow in terms of the diffusion of a particle in the triangular (equal-sided) lattice of obstacles tessellating the Euclidean plane.Finally, we propose the generalization of the model to non-symmetric potentials acting between particles on the cylinder and show that the corresponding RG flow might differ from the Fibonacci sequence which has a Golden ratio 1 2 √ 5 − 1 ≈ 0.618 as a fixed point.Specifically, we demonstrate that for some non-symmetric potentials which have different strengths along the cylinder axis and along its circumference, we find a set of "metallic fixed points" expressed in terms of the so-called "metallic ratios" among which the so-called "Silver ratio" [61] is one of the known representatives.The basin of attraction of a "Silver ratio" is lower than that of a Golden ratio, which could be a reason why the Golden ratio is distributed in nature much wider than the Silver ratio.
A. Construction of the potential and RG flow Any particle on the cylinder can be parameterized by a pair (z n , α n {mod 2π}), where n ∈ N, subject that all particles are arranged according to monotonic growth of z n .Projecting the cylindrical surface conformally onto the plane, we get new coordinates, r n,m (x, y), of particles on the planar lattice, where the connection between cylindrical and planar lattices is set by the following change of variables: Strong repulsive potential acting between particles can be approximated by the conformallyinvariant 1/r 2 potential.Consider two arbitrary particles one of which is located at the origin of the (x, y)-plane and the second -at some point (x m,n , y m,n ).Suppose that the potential U (r m,n ) has the following form: where q > 0 is some arbitrary parameter having sense of a charge.The energy of a whole lattice reads Substituting ( 47) into (48), we get: Comparing ( 50) to ( 20) and ( 23), we conclude that U (x, y) ≈ qE(x + iy, s → 1) → 4πq ln y 1/4 |η(x + iy)| + const (51) where E(x, s) is the non-holomorphic Eisenstein series (see (19)).Recall that here again we have exploited the 1st Kronecker limit formula, and dropped out the term divergent at s → 1 (which is independent on x, y).The self-similarity of the function U (x, y) is clearly seen in Fig. 9 where we have plotted a set of curves U (x|y) ≡ U (x, y) taken at different values of y.For better visualization the curves at different y (0.003 < y < 0.05) are shifted in the vertical direction as shown in Fig. 9. Having the function U (x|y) we may construct a trajectory that describes the continuous flow of the minimum of U (x|y) as a function of x (where x ∈ [0, 1]) when y is continuously changing from +∞ down to 0. The corresponding flow is depicted in Fig. 10a by a sequence of white dots for a family of plots U (x|y m ) where y m = y 0 − cm, and m = 0, 1, 2, ..., M .The parameters c > 0 and M are chosen such that y M > 0 (in our numeric computations y M = 0.005, y 0 = 1.00, c = 0.005, M = 199).
The flow of the minimum of the potential U (x) when y → 0, passes through the successive reflections of the fundamental domain of the free group Γ 2 is shown in Fig. 10b.The corresponding Cayley graph is a 3-branching Cayley tree.Recall that the 3-branching Cayley tree can also be viewed as the Cayley graph of the group Λ, which has the free product structure: Λ ∼ Z 2 ⊗ Z 2 ⊗ Z 2 , where Z 2 is the cyclic group of second order.The matrix representation of generators h 1 , h 2 , h 3 of the group Λ is well known: The optimal flow shown in Fig. 10a 2 i, we can find its image, z N , after N recursive applications of generators from the set {h 1 , h 2 , h 3 } according to the following formula: where z means complex conjugation of z and {a N , b N , c N , d N } are the coefficients of the matrix Using ( 52)-( 54) we reproduce the coordinates of the points A, B, C in Fig. 10b.The sequence which converges to the Golden ratio is as follows: where N = 3M , M = 1, 2, 3, ....The limiting value of x ∞ = Re z N →∞ is the Golden ratio: The sequence of "zigzag" reflections is encoded in the continued fraction expansion of the Golden ratio, ϕ: where interlacing odd and even "1" correspond to the left and right turns of a zigzag path.

B. Topological view at the potential U (x, y)
The construction of the optimal path in the relief constructed on the basis of the Dedekind η-function allows for geometric interpretation.Take the triangle ABC centered at w 0 in the complex plane w with the Euclidean metric.Tessellate the plane w by images of ABC obtained by reflections of this triangle with respect to its sides.Consider two patterns of the elementary cell: (a) the triangle with angles ( π 3 , π 3 , π 3 ), (b) the triangle with angles ( π 2 , π 4 , π 4 ).The corresponding tessellations are schematically shown in Fig. 11 where red points designate the images of the point w 0 under reflections.Suppose that we have made k reflections of the triangle ABC with respect to its sides (k = 7 in Fig. 11).Let us address the following question: which sequence of k successive reflections corresponds to the maximal Euclidean distance, d(w 0 , w k ), between the initial point w 0 and its image after k reflections, w k ?The answer seems straightforward: d(w 0 , w k ) is the "most aligned" sequence which is encoded in the following set of reflections lying in a gray strip in Fig. 11a,b: To see the connection of the longest Euclidean distance d(w 0 , w k ) in the plane w covered by the triangular lattice, with the optimal flow in the modular domain, let us make the conformal mapping of the equal-sided triangle lying in w = u + iv to the fundamental domain of the modular group -the zero-angled triangle bounded by arcs (see Fig. 10b).The corresponding mapping can be performed in two steps: (i) we map the triangle ABC on w with branching points at the corners onto the upper half plane of the complex plane ζ = ξ + iχ with the branching points at 0, 1, i∞, and then (ii) we map the upper half plane of ζ onto the fundamental domain of the modular group.Such a composite mapping has been described in detail in [25,62], so below we reproduce the final result only: (ii) Conformal mapping ζ → z is realized via the k 2 (z) modular function:  where θ i (0, q) (i = 1, ..., 4) are the Jacobi elliptic θ-functions.
In what follows we will need only the Jacobian J of the composite mapping w(ζ(z)) which can be easily computed: θ 2 (0, e iπz )θ 3 (0, e iπz )θ 4 (0, e iπz ) 8/3 = π 4/3 2 14/3 3Γ 6 (1/3) |η(z)| 8 (61) where the following relations between Jacobi theta-functions have been used: Consider now the diffusion-like problem in the plane w equipped with the triangular lattice of obstacles.The corresponding probability distribution of random paths, P (w, t) obeys the parabolic equation: where w = u + iv and w = u − iv.To classify topological states of trajectories from the point of view of their entanglements with obstacles, it is instructive to pass to the covering space as it has been explained in [25,28].Taking into account the conformal invariance of the Laplace operator and making use of the Laplace transform, P (w, λ) = ∞ 0 P (w, t)e −λt dt, we perform the conformal mapping w → z, and rewrite Eq.( 63) in a form of a stationary diffusion equation in z-plane in the effective "potential" W (z): where the potential W (z) is defined by the Jacobian J of the conformal mapping expressed again via the Dedekind η-function: Eq. ( 65) shows that the effective potential emerging in the topological problem of diffusion in the array of obstacles has exactly the same structure of minima and maxima as the potential acting between repulsive particles located at the surface of the cylinder in the phyllotaxis problem.The topological meaning of coordinates x and y is as follows: y describes the "complexity" of the entanglement, which in the polymer language is known as the "length of the primitive path" [24,27] and has a meaning of the geodesic length in the covering space [25,28], while x describes the local winding of the path around obstacles.
One can immediately see now that the "most straight" sequence (58) for the tessellation of the plane by the equal-sided triangle ( π 3 , π 3 , π 3 ) shown in Fig. 11a exactly matches the "zigzag" paths coded by the cyclic sequence (h .. of generators of the group Λ in Fig. 10.The extension to the asymmetric potential is considered in the next section.

C. Optimal flows in asymmetric potentials
Let us return to the problem of finding an optimal lattice of repulsive particles on the cylinder.Suppose now that the potential acting between particles is asymmetric: it is stronger along the circumference of the cylinder and weaker along its axis.In that case, the equilibrium configuration of particles will not be anymore the Abrikosov triangular lattice depicted in Fig. 12a, but rather the triangular lattice with the elementary placket in a form of an "isosceles triangle" as it is shown in Fig. 12b.
The question which we address here is as follows.If we squeeze the cylinder along its axis, should the stable lattice patterns for particles interacting with the asymmetric potential follow again the Fibbonacci sequence, or the corresponding optimal flow will choose another minimal energy valley with a fixed point distinct from the Golden ratio?As we shall see, particular asymmetries of the potential force the flow to select an optimal path in the phase space distinct from the Fibonacci sequence.In considered examples of lattices with elementary triangles ( π 2 , π 4 , π 4 ) and ( 2π 3 , π 6 , π 6 ) the flows follow so-called "metallic ratios" (in particular, "Silver ratio" for the triangle ( π 2 , π 4 , π 4 )).Taking into account the geometrical interpretation of the Jacobian of the conformal transform provided in the previous section, we conjecture that one can mimic the asymmetry in the interaction between particles on the cylinder by considering the isosceles triangle (instead of the equal-sided one) tessellating the plane w.The conformal mapping of the elementary cell in the form of an isosceles triangle with angles (απ, βπ, γπ) = (aπ, aπ, (1 − 2a)π) to the fundamental triangle of the modular group can be constructed by a straightforward generalization of the conformal transform described by equations ( 59)- (60).Namely, instead of (59) one has where α = a, β = a, γ = 1−2a.Equation ( 60) remains unchanged.Computing the Jacobian J(z, a) of the composite conformal mapping, one gets At a = 1 3 we return to the expression (61) for the Jacobian J(z, 1 3 ) ≡ J.The isosceles triangle completely (without gaps and overlays) tessellates the plane by reflections with respect to its sides for values a = 1 3 ; 1 4 ; 1   6   only.If we do not restrict ourselves by the condition to tessellate the plane by isosceles triangle, there is one with angles ( π 3 , π 6 , π 2 ) which tessellates the Euclidean plane completely, however this case is not considered in our work because it does not correspond any physical choice of interaction potential acting between particles.All other triangles lead to an incomplete tessellation of the plane or to the tessellation with overlays.
For a = 1 4 and a = 1 6 the potentials J(z, a = 1 4 ) and J(z, a = 1 6 ) have the following explicit expressions For triangles with a = 1 4 and a = 1 6 we construct the sets of sequential reflections of the fundamental domain of the modular group operating with generators {h 1 , h 2 , h 3 } as it has been done for the equal-sided triangle with a = 1 3 (see ( 52)-( 55)).The flows of optimal paths and corresponding reflections are shown in Fig. 13a,b.Sets of reflections for a = 1 4 and a = 1 6 have the following explicit expression: where interpreted as the re-distribution of the system of repulsive particles (equilibrated at the surface of the cylinder) when the cylinder is squeezed along its axis.The contour plot of U (x, y) for c = 1 in the region 0.01 < x < 0.99, 0.005 < y < 1 within the bounding box −0.27 < U (x, y) < −0.25 is shown in Fig. 14a.The (x s , y s ) coordinates of saddle points have the generic expression: where (m 1 , m 2 , n 1 , n 2 ) are some integers.In particular, white dots in Fig. 14a have the following coordinates: ( 1 2 , 1 2 ), ( 3 5 , 1 5 ), ( 813 , 1 13 ), ( 2134 , 1  34 ).From the topological point of view there is no difference between all these saddle points, however the orientation of saddles with respect to the x-axis is different and the geodesic (cyan line in Fig. 14a) passes through different saddle points at different angles.The coordinates (x s , y s ) of saddles constituting the Fibonacci series are: where To proceed, let us expand the potential U (x, y) in the vicinity of some saddle point (x s , y s ) and find the explicit form of the corresponding surface U (x, y) near (x s , y s ).The Taylor expansion of U (x, y) up to the second order reads: where the derivatives U xx , U xy = U yx , U yy are taken at the point (x s , y s ).The first derivatives in the Taylor expansion (75) are nullified at the point (x s , y s ), and the condition U xx U yy − U 2 xy < 0 ensures that the point (x s , y s ) is actually a saddle.Let us define the coefficients U xx = a 1 , U x,y = a 2 , U yy = a 3 .Note that the coefficients a 1 = a 1 (k), a 2 = a 2 (k) and a 3 = a 3 (k) depend on k, where k = 0, 1, ..., ∞ is the corresponding member in the Fibonacci sequence.In Fig. 14a we depict the relief U (x, y) where the white dots are saddle points with the following (x, y) coordinates: (0, 1), ( Let us consider the RG flow in the complex z = x + iy plane in vicinity of saddle points (x s , y s ) of the surface U (x, y).Introducing the coordinates u = x − x s and v = y − y s and separating real and imaginary parts, we may write down the system of nonlinear first-order differential equations describing the RG flow in complex plane w = u + iv in the vicinity of the point (x s , y s ): where µ is the RG time.
Equations (76) imply that the RG flow near the bifurcation points is fully determined by the topology of the Riemann surface U (x, y).It is worth mentioning that our construction is consistent with ideas expressed in works [63][64][65].In particular, in [63] the connection between the RG flows and the topological structure of β-function has been discussed in the context of CFT/ADS 2 duality, while in [64] and in [65] the equations for RG flows ideologically similar to (76) have been derived to describe the behavior of RG flows in the FQHE in the vicinity of critical points.The emergence of BKT fixed points in similar context has been also studied in [66] for layered high-T c superconductors.
Dividing the first equation of ( 76) by the second one we can convert the system (76) into the following single equation Introducing the new function h and writing u = hv, we arrive at the equation in which the variables can be separated: Solving (78) we get where G remains invariant along the RG flow (i.e.G does not depend on the scale µ).
Plugging the function h = u/v in (78) and denoting G a 1 /a 2 −2 by ∆, we have Substituting u(v) into the second equation in (76) and performing the integration, we obtain an non-explicit solution for v(µ) Despite ( 81) looks rather complicated, it is essentially simplified in the limit of small u.Computing explicitly the shape U (u, v) (recall that u = x − x s and y = y − y s ) in vicinity of the saddle point (x s (k), y s (k)), we see that with k → ∞ the coefficient a 1 = U uu tends to zero, while the coefficients a 2 = U v and a 3 = U vv remain finite.To demonstrate this, we have depicted in Fig. 14(b)-(e) the potentials U (x, y) in vicinity of four first terms of the Fibonacci series for k = 0, 1, 2, 3: One sees from ( 82) that with increasing k the coefficient a 1 (k) in front of the term u 2 relatively decreases.Substituting a 1 = 0 (corresponding to k → ∞) into (76) we get equations describing the RG flow in the XY -model in vicinity of the BKT transition.The critical scale (the correlation length) near the transition point is defined by the condition √ −2a 2 ∆ ln µ c ∼ 1 which implies the BKT dependence of the correlation length, µ c , on ∆: One can see from (82) that the coefficient a 2 in front of the uv term periodically changes the sign.So, one can expect the signature of the BKT-like transition (83) when a 2 < 0.
The physical meaning of encountered critical behavior could have the following interpretation.When the cylinder is squeezed along its principal axis, the corresponding lattice of repulsive particles experiences a set of successive rearrangements ("bifurcations").Each bifurcation is a collective effect that is accompanied by the melting of the lattice.Our analysis permits us to conjecture that some of these bifurcations in the strong compression limit have signatures of Berezinsky-Kosterliz-Thouless (BKT) transtions.

V. RIEMANN-THOMAE FUNCTION, DEVIL'S STAIRCASE AND LONG-RANGE 1D LATTICE MODELS ON A RING
Consider a one-dimensional system of n particles positioned on a ring of N sites.Particles interact via the repulsive pairwise long-range potential, V (|i − j|), which depends on the distance between particles along the ring.In addition, there is an external field, h, acting on all particles.The corresponding Hamiltonian, H, reads where n i is the indicator function of a particle at a site i, i.e. n i = 1 if the particle is present at the site i, and n i = 0 if the site i is empty.The Hamiltonian (84) corresponds to the discrete Coulomb gas on a ring considered in [67].
The exact structure of the ground state of such a system for a fixed number of particles has been discovered independently by Hubbard [68] and Pokrovsky and Uimin [69].Later, in [70,71] it has been demonstrated that ground states for a model with h ̸ = 0 form a complete Devil's staircase structure.Some mathematical aspects of the appearance of the Devil's staircase have been discussed in [72].The recursive algorithm for constructing the corresponding staircase is fairly simple and can be described as follows.
• First, we pick up a filling density ρ 0 = 1/N .This choice fixes the smallest step in density that we can detect.We also define, for completeness, h − (0) = 0.
• On the next step we recursively determine h + (ρ 1 ) by the following equations where the value ∆h is set by the sum ∆h = 2 • The algorithm runs T times (k = 1, 2, ..., T ) until the requested density ρ T ≡ ρ is reached.
Two remarks concerning the above algorithm should be made.First of all, the potential V should be convex.Secondly, we need to be careful to always write the running density ρ i as an irreducible fraction (the nominator and denominator should not have any common divisors).The denominator of this irreducible fraction is used to compute ∆h in (86).
The densities ρ(h) of systems of particles minimizing the energy described by the Hamiltonian (84) are computed by the described algorithm for two potentials, V (r) = 1/r and V (r) = 1/r 2 on the lattice of N = 7560 and N = 665280 sites.The corresponding plots of Devil's staircases are depicted in Fig. 15.It should be noted that the number of steps (plateaus) in staircases depends on the factorization of N .To see as many steps as possible, N should be a "superior highly composite number" [73].To the contrary, if N is prime, no Devil's staircase structure emerges.From a physical point of view, jumps between plateaus are the phase transitions between the states minimizing the system energy at a varying external field, h, and the density ρ(h) plays a role of an order parameter.
Comparing functions ρ(h) in Fig. 11 and G 3, it is eligible to ask the following questions: (i) Whether there is a coincidence of the function h(ρ) (horizontally reflected) for some potential V (r) = 1/r γ and the function G 2 (x) (subject to affine deformations); (ii) If such a coincidence actually exists for some values of γ, what could be the physics behind it?
In (13) we have defined the generalized Riemann-Thomae (gRT) function g α (x).Figures Fig. 2a,b provide sample plots of g α (x) for two values, α = 0.41 and α = 2.76 (for n = 100).Let us consider now the function g α,β (x) which extends the definition of g α (x) as follows: where A is some scaling factor and α and β are exponents such that α < β.
Since the Devil's staircases shown in Fig. 11a,b consists of two symmetric branches, in what follows, we will consider only one of them.Also, in view of further comparison with the Riemann-Thomae function g(x), it is convenient to work with the derivative of the staircase.So, we consider the derivative of the "inverted Devil's staircase" h(ρ).The question which we address is as follows: is it possible to find such values α and β in the definition of the function g α,β (x) (see (87)) that one can match ψ(ρ) for any potential V (r)?The answer is positive and below we demonstrate an excellent agreement of the generalized Riemann-Thomae function g α,β (x) with the function ψ(ρ) for algebraically decaying potentials, V (r) = 1/r γ , where γ = 1, 2, 1/2, as well as for V (r) = exp(−r) and V (r) = − ln r.Note that since x in the definition (11) lies on the segment [0, 1] and the density ρ changes, by definition, within the interval [0, 1], we can identify x with ρ.The plots for the potentials V (r) = {1/r 2 ; 1/r; exp(−r); − ln r} are shown in Fig. 12.The inserts demonstrate the magnification of the region near the horizontal axis.The parameters α, β and A for various algebraically descending potentials defined in (87) are presented in Table I (to save the space in Fig. 16 we have not shown the plot for the potential V (r) = 1/ √ r).From the Table I we can conjecture the following relation between α in the definition of the generalized Riemann-Thomae function g α,β and the exponent γ in the definition of the algebraic potential V (r) = 1/r γ : The exponent β is apparently non-universal and describes the finite-size corrections to the leading behavior of g α,β (ρ) at small densities (i.e. when ρ → 0).
A. Universality of the numerical parameters α, β To find the values of the parameters α and β in g α,β (x) that best describe the Hubbard model on a ring, we should check if α and β depend on the initial filling density ρ 0 = 1/N .Values presented in Table I correspond to ρ 0 = 1/7560.Here 7560 has 64 divisors and belongs to a group of highly composite numbers [73].They are defined as natural numbers that have more divisors than all smaller numbers.
Since 7560 belongs to a very specific group of numbers, we study the sensitivity of α and β to ρ 0 .To investigate the universality of these parameters, we numerically compute their values for different N .For the potential V (r) = 1/r and V (r) = 1/r 2 the results are presented in Fig. 17.Let us discuss in detail the case of V (r) = 1/r 2 .As expected, the exponent β that controls the finite size corrections is non-universal and varies essentially.Values of α stay mostly near α = 2, with the exception of anomalous points that gather around α = 1.These exceptional points occur when N is a prime number or has a small number of divisors.In this case, the devil's staircase has almost no jumps, and consequently, the derivative ψ(ρ) is constant or has only a few jumps what makes the comparison with the generalized Thomae function senseless.We have also checked that all values of N producing α = 1 and β = 1 in the interval N ∈ [853, 947] are prime numbers.Furthermore, these are also all prime numbers in this interval.We can also spot some points that are neither at α = 1 nor α = 2.We checked some of them, and found that they appear at N , that have a very low number of divisors, usually four.The results for V (r) = 1/r 2 are presented in Fig. 17.
In the same manner, we analyzed optimal values of α and β for other potentials, and found the same dependence of α on the divisibility of N .We can conclude that the parameter α is stable with respect to the choice of initial filling density ρ 0 (N ) = 1/N .Exceptions are "anomalous" points that correspond to N with a low number of divisors.The case of V (r) = − log(r) is slightly different since only one parameter suffices for describing the data.Nevertheless, we verified the universality of the obtained value α ≈ 1.Looking at Fig. 17, we can see that regardless of N , the value of α is close to one.

B. Fibonacci series in a Hubbard model on a ring
Let us formulate the rules which select the Fibonacci sequence in the Hubbard model on a ring.Look at the forest of barriers in Fig. 16 and pay more detailed attention to the potential V (r) = 1/r 2 .The corresponding plots is shown in Fig. 18 for the set of potential barriers between different ground states (Fig. 18a) and for the corresponding integrated function having a Devil's staircase structure (Fig. 18b).It is worth noting that the rule formulated below works for any long-ranged potential producing the structure of the generalized Riemann-Thomae function.
So, we want to define a natural rule for the Hubbard model which selects in Fig. 18a the sequence 1 − 2 − 4 − 6 − ... (corresponding to the Fibonacci "zigzag" sequence), but not 1 − 2 − 3 − 5..., or not 1 − 2 − 3 − 8 − ....The idea of the rule is as follows.First, build a graph as a 3-branching tree on vertices that are close by heights.In Fig. 18 the vertex 2 has neighbors 3 and 4, the vertex 3 has neighbors 7 and 8, the vertex 4 has neighbors 5 and 6, etc.Now, we would like to select an "optimal" path on the constructed tree that corresponds to the Fibonacci sequence 1 − 2 − 4 − 6 − ....Note that in (86) the expression under the sum sign has the structure of the 2nd difference relation (the discrete version of the 2nd derivative) and is defined on each triple of adjacent points on the constructed tree.The prescription of path selection is simple: we locally minimize a 2nd difference relation for each triple of neighboring vertices.For example, comparing the subsequences 1 − 2 − 3 and 1 − 2 − 4 in Figs.18c-f we see that the 2nd difference for the triple 1 − 2 − 4 is smaller than that of the triple 1 − 2 − 3. So, we select the subsequence 1 − 2 − 4. Coming to the point 4 we should choose between 2 − 4 − 5 and 2 − 4 − 6.The heights of barriers are organized in such a way that the subsequence 2 − 3 − 6 has smaller 2nd difference than that of the subsequence 2 − 4 − 5, etc.Interestingly, the Fibonacci sequence survives for any symmetric potentials.

VI. DISCUSSION
Here we comment on a few points related to our study, restricting ourselves by qualitative arguments.We postpone more detailed analysis for a separate publication.A. Jack polynomials on the Devil's staircase The existence of the Devil's staircase structure in the fractional quantum Hall effect in the thin torus limit was the subject of the studies [3,4].In [75] it has been shown that there exists a precise mapping of FQHE in this limit onto the Hubbard model on a ring discussed above.In the thin torus limit the Hamiltonian of the FQHE for interacting fermions reads as where k i are momenta of the particles forming the lattice on the circle.The Devil's staircase structure for FQHE has been reformulated in terms of Jack polynomials in [76] which looked rather surprisingly and the question "What Jacks are doing on the Devil's staircase?"seems eligible.Taking into account that the Jack polynomials are the wave functions of the Calogero many-body system, let us argue that the relations found in [75,76] for V (r) = 1 r can be understood via the chain of dualities known in the framework of the integrable many-body systems.
First, recall the duality between the Calogero-Moser (CM) / Ruijsenaars-Schneider (RS) family of long-range many-body systems and the family of inhomogeneous twisted spin chains [77][78][79][80][81].This duality, for instance, provides the relation between the quantum inhomogeneous twisted XXZ chain and the trigonometric RS model [79].Taking into account the realization of the FQHE on the torus in terms of the trigonometric RS model [82] and combining this relation with the duality mentioned above, we arrive at the qualitative understanding how the Devil's staircase emerges in terms of the particular counting problem at the Calogero-Moser side.
To be more specific, consider the quantum trigonometric RS model per se whose wave functions are McDonald polynomials, and both, coordinates and momenta, live respectively on the circles R q and R p .Geometrically the phase space of the trigonometric RS model with N particles corresponds to the moduli space of the SU (N ) flat connections of the punctured torus [83].The coupling constant corresponds to the operator inserted at one marked point.To get the system of particles, one performs the T -duality transformation for both torus cycles.Under such a transformation, the eigenvalues of holonomies over the cycles become the coordinates and momenta of particles.
On the other hand, the quantum trigonometric RS model is dual to the quantum inhomogeneous twisted XXZ spin chain [81,84] via the quantum-quantum duality.The trigonometric N -body RS model is described by the Hamiltonian where momenta live on the circle of radius R p = η −1 .The mapping between the trigonometric RS model and the spin chain goes as follows: (i) the coordinates, q i , in the RS model are the inhomogeneities in the XXZ spin chain, (ii) the momenta qi are the non-local Hamiltonians, H i , at the XXZ spin chain side, (iii) the eigenvalues of the Lax operator at the RS side are the twists at the XXZ spin chain side, and (iv) the coupling constant at RS side is the Planck constant at the spin chain side.This correspondence is described in detail in [84].
The radius of "coordinate circle", R q , is related to the anisotropy χ in XXZ chain and as the first step we can consider the limit χ → ∞ when the XXZ spin interaction term gets reduced to S z i S z j .Note that in this limit we have R q → 0 instead of more familiar limit R q → ∞, when the trigonometric RS model gets reduced to the rational RS model and the XXZ spin chain becomes the XXX chain.In our regime, to avoid fast oscillating behavior, particles tend to form the momentum lattice structure in the R q → 0 limit.To get the Hubbard-type Hamiltonian discussed above (see (84)), one has to make the second degeneration taking the semiclassical limit of the spin chain which transforms it to the Gaudin model.At the RS side this second degeneration corresponds to the R p → ∞ limit for the "momentum circle" when the relativistic RS model gets reduced to the nonrelativistic CM model in the peculiar limit.Note that the considered limit of the RS model [83] corresponds to the k → ∞ limit for the Kac-Moody level.
Taking the limits R q → 0 and R p → ∞ and transforming the trigonometric RS model to the Hubbard-like Hamiltonian, we have built the formal basis of our consideration.It is worth adding a more physical flavor to this construction.To this aim recall that the trigonometric RS model has been related to the FQHE on the torus [82] by extending the approach developed in [85,86].Now we can provide a physical interpretation of considered limits for the radii R q and R p .First, from the relation between the FQHE on the disc, the rational Calogero system supplemented with the oscillator potential ω q 2 i can be mapped it to the trigonometric Sutherland model on the circle of radius R q , where ω = 1 Rq , as it is shown in [87].Since ω is defined by the magnetic field [86], the limit R q → 0 corresponds to the strong magnetic field.Second, in the limit R p → ∞ we arrive at the specific version of the trigonometric CM model, which means that we indeed find ourselves with the Jack polynomials for the "momentum lattice", as expected.
Let us complete the discussion by comparing the Devil's staircase structures at different extremities of our mapping.In the Hubbard model we focus at the chemical potential and the density of particles.In the spin chain before taking limit, these objects are involved into twist term in the non-local Hamiltonian Tr GS with the twist matrix G.In our case we can identify the twist with the diagonal S z and eigenvalue h, hence the term Tr GS gets reduced to the term hn i in the Hubbard Hamiltonian.The Hubbard Hamiltonian in terms of the non-local Gaudin Hamiltonians H i is the sum of individual terms H = N i H i .In the Hubbard case we can recognize the Devil's staircase ρ(h) for the density as the function of twist h, or the inverse function h(ρ) -see [2].
What is the meaning of ρ(h) function and Hubbard Hamiltonian at CM side?As we discussed above, the twist h ↔ E cal corresponds to the eigenvalue of the Lax operator and therefore to the eigenvalues of CM Hamiltonians.Since we have a single eigenvalue of the twist matrix, we consider the highly degenerate state.The H i Hamiltonians correspond to the momenta at the CM side, hence the total Hubbard Hamiltonian is nothing but the total momentum of CM particles P = i p i .
In the iterative procedure which yields the Devil's staircase we add at each step one additional particle demanding the total Hubbard energy to be constant.Being translated to the CM side it means that we add particles keeping the total momentum fixed.It is worth mentioning the difference between the Devil's staircase structure in the CM model in the context of Fibonacci numbers and in Hubbard case.In the phyllotaxis problem we fix the energy E = 0 and count the weighted degeneracy of this level.In the context of the Hubbard model we focus on the multiplicity of the P = const state in a specific limit of RS system.Since the momentum of the RS particles is determined via the Bethe ansatz equations [77] one could a bit loosely say that the Devil's staircase structure emerges in the space of solutions to Bethe ansatz equations.

B. On Fibonacci universality class at out-of-equilibrium
It was suggested in [88] that critical exponents describing fluctuations in the nonequilibrium dynamics are of more generic nature than it is usually assumed.It was claimed that the Gauss and KPZ scalings are just two first representatives in the generic "Fibonacci hierarchy" with critical exponents z n = F n+1 Fn , where F n is the nth Finonacci number.This result has been obtained by analyzing the hydrodynamic equations with several conservation laws.As a toy example, the 3-species TASEP model has been discussed and the numerical simulations indeed exhibit the third critical exponent for the large-time asymptotics of the maximal value of the two-point correlator max ⟨ϕ(x, t)ϕ(0, 0)⟩ ∝ t −1/z 3 (93) The existence of Fibonacci universality class is still under discussion.For instance, it was demonstrated rigorously that for the systems with the non-abelian global symmetries including the integrable spin models the large-time asymptotics of two point function enjoys the KPZ critical exponent [89] generalizing the initial observation in [90].
It is eligible to ask a question of whether there is any relation between the Fibonacci hierarchy and our study.Let us speculate on the possible connection and assume that we got somehow the correlator of a particular operators in the modular domain, for instance the two-point correlator, ⟨ϕ(τ )ϕ(0)⟩ = f (τ ) (94) where τ is a modular complex parameter, for which we interpret Im τ as a time variable.Such a viewpoint is valid at least in two situations.First, in the phyllotaxis problem Eq.( 94) has a sense of the propagator in the modular domain with Im τ as an evolution parameter.Second, in the Whitham dynamics for the Seiberg-Witten solution, the Nekrasov partition function can be mapped via the AGT correspondence [36] onto the conformal block in the Liouville or Toda models where the modular parameter corresponds to the insertion point for some vertex operator ⟨V 1 (0)V 2 (1)V 3 (∞)V 4 (τ )⟩ and Im τ is the time for the Whitham dynamics.
Let us suppose a simple typical scaling behavior for a correlator (94) in the modular domain, say f (τ ) ∝ exp(τ ) and consider the parameter τ = θ + i log t from the perspective discussed at the phyllotaxis side.In this representation at small values of Im τ one gets the critical exponents dictated by asymptotic values of θ for which we know that they are given by truncated continued fraction expansion (57) for sequential quotients of Fibonacci numbers.This point of view supports the idea of Fibonacci universality raised in [88].Certainly, these arguments are very superficial and an accurate analysis is required, however it immediately poses the following challenging question: do we have non-equilibrium dynamical systems producing a "Silver ratio universality class"?

C. Dedekind in the proper place
The Riemann-Thomae function has been discussed recently from an interesting perspective in the theory of massless free bosons and fermions on the circle of radius R at finite temperature T = β −1 upon the Wick rotation [91,92].The twisted boundary conditions involving both circles are imposed and for the rational p/q the twist corresponds exactly to the T p,q torus knot for the closed space-time trajectory of a particle.The emergence of the RT function is not a surprise since the partition function of the 2D theory on the torus with the global U (1) × U (1) involves the Dedekind functions hence the free energy in the proper limit indeed has a RT structure leading to the Devil's staircase.
It was argued in [91,92] that for the torus knot boundary conditions the thermodynamic properties demonstrate interesting fractal behavior.It was suggested that the statistics of free particles with the torus knot boundary condition is level-dependent and the negative pressure regime can be found.Since the Dedekind enters the partition function of the massless scalar, the RT function defines the dependence of the pressure on the twist parameter in the limit L = R β → ∞.It would be interesting to match the fractal thermodynamic properties for system with torus knot boundary condition and our RG approach.Indeed for the torus knot trajectories we have specific knot invariants corresponding to particular multiplicities of states.Hence, their impact on the partition function is expected.It would be also interesting to recognize the possible fractal thermodynamics for the massless scalar on the mirror torus when the complex and Kahler structures get interchanged.

VII. CONCLUSION
This work provides a modest attempt to add some flavor of universality to the interplay between the modular group acting in the parameter space of physical systems (spectra of random operators, phyllotaxis, Hubbard model on a ring) and the RG flows in the peculiar region of the fundamental domain of the modular SL(2, R) group when the real part of the modular parameter tends to zero.We have argued that in this regime the systems possess the universality described by the generalized Riemann-Thomae (gRT) function, and the generalized Devil's staircase emerges.Using a natural regularization of the fractal gRT by the modular Dedekind η-function we were able to interpolate between "neighboring" fractal states and connect this interpolation with the RG flows on the modular group.Saying differently we were looking for additional arguments supporting the universality of commensurabilityincommensurability transitions.The problem can be reformulated as the derivation of the RG flow for the deformations of lattices of different nature via some disorder.We argued that the analysis of the lattice structure by studying the corresponding spectral properties of propagating probe could be very useful and in the "thin torus limit" the rearrangement (bifurcation) of highly squeezed lattice is a collective effect with a signature of the BKT transiton.
The limit y → 0 of the modular parameter τ = θ + iy we are focused at, has the clear physical interpretation.It corresponds to the situation when the disorder y associated with the imaginary part of the modular parameter in some frame tends to zero, while the θ-term associated with the real part of τ and serves as the chemical potential for some version of the topological charge, remains finite.The very notion of the disorder y is model-dependent and in some systems it can be identified with the diffusion constant as for example in the Anderson transition problem, or with the magnetic coupling constant as in the Yang-Mills theory.The notion of the "weak disorder" is frame-dependent since the weak coupling limit in the magnetic frame corresponds to the strong coupling regime in the electric frame.We have argued that in this regime the modular the non-perturbative "instanton" renormalization dominates since the disorder is weak.In particular, our analysis suggests that one could expect the Devil's staircase in some version of multiplicities of BPS states near Argyres-Douglas point.
There are several issues that certainly deserve additional study.It would be interesting to find the place to gRT function in the group-theoretic framework.There are finite algebras that involve several parameters, such as Sklyanin algebra which has the p-adic and quantum groups as the peculiar degenerations.The models which enjoy the devil staircase and generalized Devil's staircase have extended symmetries and it would be interesting to recognize the structures discussed at length of our paper in some limits of Sklyanin algebra.Another possible group-like structure concerns the algebra of BPS states which was identified as the hyperbolic Kac-Moody algebra [93].The flow between the lattices in this framework correspond to the interpolation between hyperbolic algebras which are related with the Fibonacchi numbers [94,95].
It would be also interesting to include into the framework of our study the structures associated with general T n,m torus knots and links related with the instanton counting via the instanton-torus knot duality.This should generalize the relation between invariants of T 2,n knots and Fibonacci numbers.The last note concerns the resurgence theory (see [96] for the review) providing the interplay between the non-perturbative and perturbative contributions to different objects including β-function.
Fig. 5a.The corresponding generic recursive algorithm constitutes the Farey sequence construction.Schematically, the construction goes as follows: take the segment [0, 1] and draw two circles O and O ′ both of radius r = 1 2 touching each other, and the segment at the terminal points 0 and 1.Now inscribe a new circle O 3 touching O 1 , O ′ and [0, 1].Where is the position of the new circle along the segment?The position of the newly generated circle projected to the segment [0, 1] is determined via the ⊕ operation (36).For example, the centers x O 2 of the circle O 2 , x O 3 of the circle O 3 and x O 4 of the circle O 4 are located correspondingly at the points:

2 Figure 5 .
Figure 5. (a) The Ford circles as illustration of the Farey sequence construction: each circle touches two neighbors (right and left circles) and the segment.The position of newly generated circle is determined via the ⊕ addition: p i q i ⊕ p j q j = p i +p j q i +q j ; (b) The spectral density ρ ε (λ) for the ensemble of tridiagonal matrices of size N = 10 3 at p = 0.5 and its relation to Farey numbers; (c) The same Farey sequence generated by recursive fractional-linear transformations of the fundamental domain of the modular group SL(2, Z).

Figure 6 .
Figure 6.Family of 20 curves f (k, n) = −2 cos πk n+1 .Each curve corresponds to a particular n = 1, ..., 20 (from left to right), points along each curve mark values k = 1, ..., n.Each horizontal dashed line corresponds to the multiplicity of the eigenvalue and contributes to the height of the peak in the spectral density.

Figure 7 .
Figure 7. Visibility diagram.Each point in the vertical row n carries a weight p n−1 .Integer points within the wedge [π/4, π/2] are designated by circles.The dashed rays are emitted at rational tangents m n , where m and n are coprimes.The weights corresponding to marked integer points, are summed up along the rays.

2 Figure 8 .
Figure 8.(a) Repulsive points distributed on the surface of the cylinder form a lattice, characterized by the parameter α, with a minimal energy.The lattice is rearranged when the cylinder is compressed along its vertical axis; (b) Dependence U (x, y = const) defined in (51) for the compressed lattice (y = 0.001 ≪ 1 and β = 1) as a function of x.

Figure 9 .
Figure 9. Set of plots of U (x|y) on x taken at different values of y in the region 0.003 < y < 0.05.For better view the curves for different y are shifted in the vertical direction.As smaller y as more generations of peaks are proliferated.

3 Figure 10 .
Figure 10.(a) Evolution of a minimum of a potential U (x) when y is continuously changing from y 0 = 1.00 towards 0. For better visualization each minimum at a given value of y is marked by a white dot; (b) The fundamental domain of the modular group and few sequential reflections highlighted by the arcs of geodesics:A 0 → A 1 → A 2 → A 3 → A 4 ..., i.e. (h 3 h 1 h 2 )(h 3 h 1 h 2 )(h 3 h 1 h 2 )....

Figure 12 .
Figure 12.(a) Particles on the cylinder interacting with the symmetric potential form the equalsided triangular "Abrikosov lattice"; (b) Particles on the cylinder having stronger interaction along the circumference of the cylinder and weaker along its axis form the lattice with the elementary placket in a form of isosceles triangle.

Figure 13 .
Figure 13.Evolution of a minimum of the potential J(x|y, a) when y is continuously changing from y 0 = 1.00 towards 0. For better visualization each minimum at a given value of y is marked by a white dot; (a) the flow for a = 1 4 ; (b) the flow for a = 1 6 .

Figure 18 .
Figure 18.(a) Generalized Riemann-Thomae function emerging for the Hubbard model on a ring for the potential V (r) = 1/r 2 .The Fibonacci sequence 1 − 2 − 4 − 6 − ... is "optimal" with respect to other sequences from the point of view of local minimization of the relation (86); (b) The associated Devil's staircase.Steps corresponding to the Fibonacci sequence are shown in red; (c)-(f) Sequential steps of the optimal path selection. )