Generalized Ising Model on a Scale-Free Network: An Interplay of Power Laws

We consider a recently introduced generalization of the Ising model in which individual spin strength can vary. The model is intended for analysis of ordering in systems comprising agents which, although matching in their binarity (i.e., maintaining the iconic Ising features of ‘+’ or ‘−’, ‘up’ or ‘down’, ‘yes’ or ‘no’), differ in their strength. To investigate the interplay between variable properties of nodes and interactions between them, we study the model on a complex network where both the spin strength and degree distributions are governed by power laws. We show that in the annealed network approximation, thermodynamic functions of the model are self-averaging and we obtain an exact solution for the partition function. This allows us derive the leading temperature and field dependencies of thermodynamic functions, their critical behavior, and logarithmic corrections at the interface of different phases. We find the delicate interplay of the two power laws leads to new universality classes.


Introduction
It is almost futile, and perhaps impossible, to comprehensively list the advances in understanding of various phenomena in physics and beyond that were achieved due to the Ising model. Excellent reviews of the one-hundred year history of the model [1][2][3][4][5][6] are supplemented by discussions in other papers of this Special Issue. This paper has been written for the Special Issue of Entropy 'Ising Model: Recent Developments and Exotic Applications'. We think it is therefore more beneficial to open our paper with two first-hand accounts that concern Ernst Ising, the person and the model. The first of these is of a historical nature and concerns another body of work by the present authors and their colleagues. The second, rather methodological account, will bring us closer to the subject of studies of new physics presented in this paper.
For a quarter of a century, the Ising lectures have facilitated the emergence of different initiatives, both spontaneously and by design, that both review and advance Ising modelrelated research [7]. This workshop started in Lviv (Ukraine) in 1997 with 'traditional' statistical physics and has recently broadened its scope to encompass a more general context of complex systems. The lectures became the subject of a review series [8][9][10][11][12][13] and gradually the workshop gave rise to various research projects centered around the Ising model and its history. Historical documents collected to date, and displayed publicly with permission of Ernst Ising's family, include his dissertation [14] and its shortened version which was published in Hamburg in 1924 [15]. They also include memoirs of Ernst's wife, Johanna (Jane) Ising [16], as well as a recent publication that includes memoirs of their son Thomas [17]. It was through this collaborative atmosphere of the workshop, and in the context of a broader L 4 Collaboration in Statistical Physics of Complex Systems [18], that the problem considered below emerged.
As mentioned, the second remark brings us closer to the scientific subject of this paper; it concerns a special feature which made the Ising model so popular for descriptions of collective behavior in multitudes of systems. In its original form, as presented in Ising's thesis, this feature is binarity-representation of the state of an agent as from a pair of binary oppositions. It is to a large extent due to this feature that the model has been (and we believe will continue to be) applied in almost all fields where binarity plays a core role [17,19,20]. Some generalizations of the Ising model lose this feature. An example is the q-state Potts model [21,22] which keeps the discrete symmetry of the Ising model, generalizing it from Z 2 to Z q . As a result, although each agent (spin) can take on only a finite number of states, the binarity is lost for any q = 2. Another popular generalization, the O(m)-symmetrical model [23,24], enables an infinite number of states for a single agent because the symmetry is continuous at m = 1.
Here, we address ordering phenomena in systems of agents that are not necessarily physical in nature with the special role that is played by spin models in complex networks in mind [20,25]. Recently, we have suggested another generalization of the Ising model that tackles such circumstances by keeping binarity of the Ising model but relaxing the condition of fixed spin length on each site [26]. Within the model, the length of each spin is considered as a quenched random variable with a given distribution function and hence the observables are calculated by the usual Gibbs averaging over the (up and down) spin configurations as well as over the random spin length distribution. The model is related to (but differs from) other spin models that are used to study the impact of structural disorder on collective behavior [27][28][29][30][31][32][33][34] and it may be useful in analysis of ordering in magnetic or ferroelectric systems of particles with polydisperse elementary moments [35,36]. Another obvious field of applicability of this model is understanding peculiarities of ordering processes in systems containing agents that, although being of binary character ('+' or '−', 'up' or 'down', 'yes' or 'no'), differ in strength of expression [37,38].
An example is illustrated in Figure 1. The structure of the network is used to model the underlying interactions in a system of interest, be they of specific chemical, biological, social, or economic origin. In a recent short communication [26], we reported on the peculiarities of the generalized Ising model when the random spin length is governed by a power-law decaying distribution function. We obtained an exact solution for this model on complete and Erdős-Rény graphs as well as commented on the phase diagram of this model on an annealed scale-free network. The analytic solution for this last case has never been displayed to date and is a subject of this paper. The rest of the paper is organized as follows. In Section 2, we formulate the model and demonstrate that the partition function of the model possesses an important feature: it is self-averaging. This fact essentially facilitates calculations of thermodynamic functions as displayed in Section 3. We apply the steepest descent method to get exact results on the thermodynamic limit. We also analyze the phase diagram and show how an interplay between two different power laws, one governing the network structure and another one governing spin properties, defines universal features of critical behavior. Conclusions and outlook are given in Section 4 and asymptotic estimates for the integrals that enter thermodynamic functions are derived in Appendix A. S Ising model with varying spin length (strength) as a model for a social phenomenon. Each individual is represented as a complex network node of a given degree k i (i.e., a number of persons connected to it via social links) and given strength S i . One may consider spreading of positive (spins up) and negative (spins down) emotions in a social network.

Model
Well-studied generalizations of the Ising model include the m-vector [23,24] and the Potts [21,22] model. Instead of a discreet scalar variable σ i = ±1, the former considers a classical vector variable σ i that can point in any direction in an m-dimensional space. The Potts model, on the other hand, maintains discrete variables, but relaxes the number of single-site spin states. Here, we consider another generalization of the Ising model. The new model preserves the binary character of the spin variables but allows them to change their absolute value in a continuous and random manner [26]. To achieve this, we endow the spins with 'strength' that can vary through a random variable S with a given probability distribution function q(S ). Below, we consider the case where this distribution function is characterized by a power-law decay: with the normalization constant c µ and µ > 2 to ensure finiteness of the mean strength S at S max → ∞. As mentioned in the Introduction, the model mimics inhomogeneities in many-particle (multi-agent) systems of different natures, that may range from polydisperse magnets or ferroelectrics [27][28][29][30][31][32][33][34][35][36] to various complex social or economical systems [37,38]. In turn, the choice of the distribution function in the form of a power law allows both to proceed with analytic calculations as well as to gain access to various regimes of polydispersity by tuning exponent µ.
Considering the critical behavior of a spin system on a complex network, special attention has been paid to scale-free networks, which are characterized by a power-law decay of a node degree distribution function: where p(K) is the probability that any given node has degree (number of links) K, c λ is a normalization constant, and λ > 2. It is well established by now that the Ising model on a scale-free network has a non-trivial critical behavior: depending on the value of λ, it is characterized by different critical exponents [39][40][41]. For example, when λ > 5, the critical exponents coincide with the mean-field ones observed for regular lattices. In the region 3 < λ < 5, the exponents become λ dependent. When λ = 5, logarithmic corrections to scaling appear. Below, we consider a generalized Ising model with varying spin strength on a scalefree network. Doing so, we analyze how an interplay of power laws (1) and (2)-the first governing network structure and the second governing agents' strengths-impacts critical behavior. To proceed, we first formulate the annealed network approximation we will be dealing with.

Ising Model on an Annealed Network
Following Refs. [42][43][44][45], we define an annealed network as an ensemble of networks of N nodes each, with a given degree arrangement {K} = (K 1 , K 2 , ..., K N ), maximally random under the constraint that their degree distribution is a given one. The linkage between nodes is taken to fluctuate for each fixed sequence {K}. Therefore, in the spirit of the concept of annealed disorder [46], the partition function is to be averaged with respect to these fluctuations. This is different from quenched disorder, when for each fixed sequence {K} network links are fixed too and therefore the free energy is averaged. In this latter case, the configurational model serves as a counterpart of the annealed network (see, e.g., [47]).
To construct an annealed network of N nodes, one assigns to each node i a random variable (label) k i taken from the distribution p(k) and the probability of a link between two nodes is defined as: with k = 1 N ∑ l k l . One can show that the value of the random variable k i indicates the expected value of the node degree: EK i = ∑ j p ij = k i whereas its distribution p(k) defines node degree distribution p(K).
In the presence of a homogeneous external magnetic field H, the Hamiltonian of the (usual) Ising model on an annealed network reads: where the second sum spans all N network nodes, the first is over all their pairs and J ij is an adjacency matrix with matrix elements equal to J if nodes are connected and 0 otherwise: For the fixed sequence of random variables {k} = (k 1 , ......, k N ), the partition function is obtained by averaging with respect to random annealed links {J}: where β = T −1 is the inverse temperature and the averaging over links reads, cf. Equation (5): In turn, obtained after averaging over random linking, the partition function Z N ({k}) depends on the particular choice of random variable (label) sequence {k}. Recall that this sequence was taken as a fixed one, i.e., quenched. Therefore, the observable free energy F N is to be obtained by averaging the sequence-dependent free energies F N ({k}) as: It is worth mentioning here another prominent feature of the annealed network: as we will explicitly show below, the partition function Z N ({k}) is self-averaging, i.e., it does not depend on a particular choice of {k}: Z N ({k}) ≡ Z N . This leads to an obvious relation: which means that the free energy is a self-averaged quantity too and avoids averaging of the logarithm of partition function, facilitating calculations on annealed networks.

Ising Model with Random Spin Length on an Annealed Network
The model we consider in this study [26] relaxes the restriction on the fixed spin length in the Hamiltonian (4). Similar to the Ising model, we preserve the binary character of spin variables keeping global Z 2 symmetry of the whole system, however, we allow each spin to change its absolute value in a continuous and random fashion. Namely, we endow the spins σ i with 'strengths' which vary from site to site through a random variable |σ i | ≡ S i . The Hamiltonian of the model reads: where all notations are as in Equation (4) and S i are independent identically distributed (i.i.d.) random variables with a given distribution function q(S ) each. The Hamiltonian (11) can be equivalently rewritten in terms of usual Ising spins of unit length, choosing variables We consider the case when the sequence {S} = (S min , ..., S max ) is maximally random under the constraint that their distribution is a given one. For the fixed sequence of random variables {k} (that define network linkage) and {S} (that define local spin strength), the partition function is obtained by averaging with respect to random annealed links {J}, cf. Equation (6): with the trace defined in (7). Generally speaking, after the trace over spins has been taken, the partition function also remains dependent on the (randomly distributed) spin strengths {S}, as explicitly denoted in Equation (13). However, in the next subsection, we show that in the case of annealed networks, the partition function Z N ({k}, {S}) is a self-averaging quantity both with respect to random variables k and S (Z N ({k}, {S}) = Z N ). Therefore, for the free energy, similar to (10), one obtains: Our task now is to proceed in deriving the partition function of the Ising model with varying spin length S on an annealed scale-free network when distributions of the random variables q(S ), p(k) follow power-law behavior (1), (2). In the course of derivation, we arrive at the conclusion about its self-averaging properties.

Self-Averaging
Substituting into (12) the adjacency matrix (5) and averaging over spin configurations, we obtain: Taking into account that the spin product in (15) can attain only two values (σ i σ j = ±1), we can make use of the equality to obtain the partition function (15) Simplifying the expression for the partition function, one arrives at: Making use of the equality (16) to represent ln(a ij + b ij σ i σ j ) in (18), we obtain for the partition function: with The latter coefficients implicitly depend on p ij via (19). Substituting these dependencies into (21), one obtains: Substituting p ij into the expression for the partition function (20) and evaluating d ij (23) in the thermodynamic limit N → ∞ (i.e., in the limit of small p ij ), we get: Now the interaction term in (25) attains a separable form and one can apply Stratonovich-Hubbard transformation to take the trace over spins σ i exactly and to obtain the following expression for the partition function: In this and all other partition function integral representations, we omit the prefactors that are irrelevant for our analysis. As long as the functional dependence on the random variables S i , k i in (26) is of the unary type, it is convenient to pass from sums over nodes i to sums over the random variables k i , S i with a given distribution function p(k), q(S ).
Considering the random variables to be continuous, one arrives at: For an infinite system, we put k max = S max → ∞ and, without a loss of generality, we choose the lower bonds equal to k min = S min = 2 and J = 1. Note, that the peculiarities of the critical behavior we are interested in are caused by the behavior at k max , S max → ∞.
Although it is more natural to choose the lower integration bond equal to unity, scale-free networks with k min = 1 do not possess a spanning cluster for λ > λ c (with λ c = 3.48 for discrete node degree distribution and λ c = 4 for the continuous one) [48][49][50]. We avoid this restriction by choosing k min = 2. To have expressions symmetric in k, S, we choose S min = 2 too. Now it is straightforward to see that the partition function Z N ({S}, {k}) does not depend on random variables k and S and is self-averaging: As one can see from Equation (28), the self-averaging property is quite general and concerns any form of distributions p(k), q(S ). Below, we use this expression to analyze thermodynamics in the case when these distributions attain power-law forms (1), (2).

Thermodynamic Functions
It is convenient to pass in Equation (28) to integration over positive values of x and to present the partition function as Being interested in the leading asymptotics of the partition function at N → ∞ and keeping the first leading term in H, we present the expression (29) in the following form: with where and we have substituted distributions q(S ), p(k) in power-law forms (1) and (2). The lower integration bound ε = 2 x N tends to zero, when N → ∞. The asymptotic expansions of the integral (32) at small ε (large N) are evaluated in the Appendix. Substituting these expansions at different values of parameters λ, µ into Equation (30), we arrive at corresponding expressions for the partition function that is evaluated at large N by the steepest descent method. The final expression for the partition function reads: where and the linear term in H originates from the large N asymptotics of the hyperbolic cosine in Equations (30) and (31). Now it is straightforward to write for the Helmholtz free energy F N (T, H) per node: with m being the coordinate of function Φ µ,λ (x) minimum: The resulting free energy is symmetric upon an interchange of indices µ ↔ λ. Therefore, below, we give the corresponding expressions for two cases: µ > λ and µ = λ. For the first case, µ > λ, an asymptotic of the free energy at small m is governed by the lower value of the exponents, i.e., by λ. Keeping the leading terms, we arrive at: with where S 2 = ∞ 2 S 2 q(S )dS, k 2 = ∞ 2 k 2 p(k)dk, the distribution functions q(S ), p(k) are given by Equations (1) and (2), and we have taken into account that S min = k min = 2 (see explanation below Equation (27)). The coefficients i µ are listed in the Appendix and c µ , c λ are normalizing factors of the distribution functions (1), (2).
For the case λ = µ, the leading behavior at small m reads: with the notations explained above. The signs of the coefficients i µ,λ do not matter in our analysis.
The estimates obtained above for the free energy asymptotics (37), (39) give one access to the thermodynamic properties of the system of interest. As we will see below, parameters µ and λ play a crucial role in governing the onset of ordering and define the universality class of the generalized Ising model on a scale-free network. Before proceeding in analyzing these expressions, it is instructive to recall the main peculiarities of the critical behavior of two models, where each of these parameters has been considered separately: these are the Ising model on a scale-free network with a node-degree distribution (2) [39,40,51] and the generalized Ising model with a power-law spin strength distribution (1) on a complete graph [26]. As is well established by now, the Ising model on a scale-free network remains ordered at any finite temperature at low values of the node-degree distribution exponent 2 < λ ≤ 3. The order parameter decays with temperature as a power law m ∼ T 1/(λ−3) at 2 < µ < 3. The decay is exponential for λ = 3: m ∼ e −bT . With a further increase in λ, a second order phase transition occurs for λ > 3 at finite T = T 0 and H = 0: m = 0 at the high-temperature phase, whereas the order parameter emerges as m ∼ τ 1/(λ−3) in the vicinity of the transition point at H = 0 with τ = |T − T 0 |/T 0 . The power-law temperature behavior of the order parameter attains its usual mean-field value only when λ exceeds five: m ∼ τ 1/2 , λ > 5. Logarithmic correction to scaling appears at marginal λ = 5: m ∼ τ 1/2 | ln τ| −1/2 . The phase diagram described above is sketched in Figure 2a. A similar picture is observed when one analyzes the generalized Ising model with a power-law spin strength distribution on a complete graph, i.e., when, in the spirit of the Kac model [52][53][54][55][56][57][58], each graph node is connected to all other nodes. As has been demonstrated in Ref. [26], the role of the global parameter is played in this case by the spin strength distribution exponent µ. In turn, we summarize the behavior of the order parameter m for different values of µ in Figure 2c.  Table 1. Now, with the free energy asymptotics for the generalized Ising model on a scale-free network (37), (39) at hand, we are in a position to analyze the interplay of two parameters: the first one governing individual spin strength (µ) and the second one governing its connectivity (λ), on the emergent critical behavior. Temperature behavior of the order parameter and the phase diagram that originate from this analysis are shown in Table 1 and in Figure 2b. The behavior is controlled by the parameter (λ or µ) with the smaller value. When at least one of the parameters (λ or µ) is less than three, the system remains ordered at any finite temperature and the order parameter decays as a power-law function of T: When either λ or µ equals three, and the other one is larger than three, m decays exponentially. A second order phase transition occurs when both λ, µ > 3. Depending on the values of λ, µ, the order parameter is characterized by different asymptotics. In the region 3 < µ < 5 (µ < λ), the critical exponents are µ dependent, and in region 3 < λ < 5 (µ > λ), they are λ dependent and logarithmic corrections appear in these regions at λ = µ: Logarithmic corrections to scaling, however, of different values, also appear when λ = 5 or µ = 5. We discuss these corrections in more detail later. Table 1. Temperature behavior of the order parameter m at different values of µ and λ. The asymptotic is governed by the smaller parameter from the pair (µ, λ).
The phase diagram in Figure 2b visualizes the behavior discussed above. There, we show different regions in the λ − µ plane that are characterized by different critical behaviors. The last is governed by the distribution with a 'fatter' tail (smaller value from the pair λ, µ). It is instructive to compare this diagram with those of Figure 2a,c. Indeed, when one of the exponents in Figure 2b is larger than five (very fast decay of one of the distributions (1) or (2)), the resulting diagram does not depend on this exponent any more. One may speak about degeneracy of the critical behavior with respect to this exponent and about reduction of the phase diagram Figure 2b to one of its corresponding counterparts, as shown in Figure 2a,c. Interesting new phenomena emerge along the lines of the diagram in Figure 2b, that separate regions with different asymptotics of the order parameter. Usually, changes in the power law asymptotics of thermodynamic observables are accompanied by logarithmic correction-to-scaling exponents (see, e.g., [59] and references therein). For d-dimensional lattices, such corrections appear at upper critical dimensions, and for the scale-free networks they are known to accompany the leading asymptotics at λ = 5. In our analysis, we complete the picture by observing the lines in the λ − µ plane, where such corrections appear. Furthermore, new scaling laws are observed at the intersection of these lines, as further outlined below.
To proceed with the analysis of critical behavior, we obtain expressions for the other thermodynamic functions in the vicinity of the second order phase transition that occurs for µ, λ > 3 at T = T 0 , H = 0. In particular, besides the order parameter, we evaluate the leading critical exponents for the isothermal susceptibility χ T , specific heat c H , and magnetocaloric coefficient m T (the magnetocaloric coefficient is defined by the mixed derivative of the free energy over magnetic field and temperature, m T = −T(∂m/∂T) H ): We also find the logarithmic terms that appear at marginal values of λ, µ and define the logarithmic correction exponents for each of the above quantities: where A is one of the thermodynamic functions (43), Θ is the critical exponent, andΘ is a corresponding logarithmic correction exponent. Values of the leading critical exponents for thermodynamic functions (42) and (43) are summarized in Table 2. The corresponding logarithmic corrections to scaling exponents are collected in Table 3. Table 2. Critical indices of the generalized model with power-law distributed spin strength on an annealed scale-free network in different regions of the phase diagram Figure 2b. Line 4: 3 < (λ, µ) < 5, λ = µ; region III: 3 < µ < 5, µ < λ; region IV: 3 < λ < 5, λ < µ; region V: λ, µ ≥ 5.
Region V, Lines 5-6, B 0 0 Similar to the case of scale-free networks, the logarithmic corrections to scaling appear at λ = 5, µ > 5, and µ = 5, λ > 5, along Lines 5 and 6 in Figure 2b. The values of the logarithmic correction exponents coincide with those for the usual Ising model on a scale-free network [39][40][41]. However, two new types of logarithmic corrections emerge in the model under consideration: in region 3 < (λ = µ) < 5 (line 4 in Figure 2b ) as well as at λ = µ = 5 (point B). For λ = µ = 5, all logarithmic correction exponents are twice as large in comparison with those for the Ising model on a scale-free network at λ = 5. In the region 3 < (λ = µ) < 5, all logarithmic correction exponents are λ dependent. All of them obey the scaling relations for logarithmic corrections [60][61][62]. Table 3. Logarithmic correction exponents of the generalized model with power-law distributed spin strength on an annealed scale-free network in different regions. Exponents for lines 5-6 coincide with those found previously [39][40][41]. Here, we find two new sets of exponents that govern logarithmic corrections along line 4 and in point B.

Conclusions and Outlook
The effects of structural disorder on the onset of magnetic ordering in regular (lattice) systems is of mainstream interest in the modern theory of phase transitions and critical phenomena [8][9][10][11][12][13]. It is well established by now that even a weak dilution by non-magnetic components may lead to crucial changes in the behavior of magnetically ordered systems. If such a dilution is implemented in a quenched fashion, changes in the universality class of the Ising model [34] are governed by the Harris criterion [63]. Annealed dilution, on the other hand, causes changes in the Ising model critical exponents via Fisher renormalization [64,65]. Another textbook example of structural disorder is given by frustrations that may be implemented in the lattice Ising model by (quenched) competing ferro-and anti-ferromagnetic interactions and they are known to cause the spin-glass phase [32,33].
The generalized Ising model we consider here relaxes the usual condition of a fixed spin length (spin strength) and considers it as a quenched random variable with a given probability distribution. In the particular case where this random variable is 1 with probability p and 0 with probability 1 − p, one arrives at the familiar quenched diluted Ising model. In this study we consider, however, another, richer case, whereby the random spin strength obeys a power-law distribution (1) governed by the exponent µ. The model mimics polydispersity in magnetic moments of elementary interacting spins. Being interested in possible applications of such a model in the broad area of complex system science, we have analyzed its behavior on an annealed scale-free network. In doing so, we make use of two advantages: the annealed network approximation leads to self-averaging properties of thermodynamic functions and the scale-free behavior of the node-degree distribution (2) allows us to study competition of power laws (1), (2) in defining critical behavior.
As appeared in the course of our study, the model under consideration possesses a number of interesting unexpected features. Some of them are summarized in Figure 2b and Tables 1-3. The phase diagram of Figure 2b is accompanied by two others, Figure 2a,c, that correspond to the usual Iisng model on a scale-free network (a) and to the generalized Ising model with the power-law distributed spin strength in the complete graph (c). As one can see from this sketch, the diagram is symmetric under µ ↔ λ interchange. This means that both factors (i.e., node connectivity and individual spin strength) influence criticality in a similar fashion. Moreover, the corresponding asymptotics are governed by the smaller of the pair of parameters (µ, λ): the 'fatter' tail of the distribution function wins the competition in defining universality class! For very low values 2 < (µ, λ) ≤ 3, the system remains ordered at any finite temperature. In turn, the second order phase transition regime (µ, λ > 3) is characterized by three different sets of critical exponents (see Table 2).
Peculiar phenomena emerge in the regions with µ = λ, where the changes in critical exponent µ or λ dependencies occur. As one observes from Table 1, such changes are accompanied by an emergence of logarithmic corrections in the form of Equation (44). The values of the logarithmic correction exponents are summarized in Table 3. It is instructive to compare this phenomenon with what happens to the critical behavior in d-dimensional Euclidean space. There, a special role is played by a concept of an upper critical dimension d u . By definition, this is the space dimension above which the universality class is trivially defined by the mean-field behavior [66]. A special type of logarithmic corrections to scaling appears at the upper critical dimension (see [59]). For the scale-free networks, the logarithmic corrections were known to appear at λ = 5, where leading exponents attain their mean-field values [39][40][41]. Similar corrections also emerge for the generalized Ising model with the power-law distributed spin strength on the complete graph at µ = 5 [26]. For the model considered here, these corrections (observed before at single points in Figure 2a,c) are now observed throughout along lines 5, 6 in Figure 2b. The crossing point of these lines, point B in Figure 2, is characterized by a new values of logarithmic corrections. Moreover, another new set of logarithmic corrections appears at 3 < µ = λ < 5.
We are deeply indebted to Thomas Ising for conveying to us many insights into Ernst and his story. We are grateful to Sigismund Kobe for further historical insights. We also thank Reinhard Folk for our common work on historical detail of the Ising model and its development over the past century. This work was supported in part by the National Research Foundation of Ukraine, project 2020.01/0338 (M.K.) and by the National Academy of Sciences of Ukraine, project KPKBK6541230 (Y.H).

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

Appendix A
In the Appendix, we evaluate integrals that enter formulas (31) and (34). In particular, we are interested in the behavior at small ε of the following integrals: We will consider the region where λ, µ > 2.

Integral I µ (ε)
Let us first consider integral (A1). At 2 < µ < 3, it does not diverge for ε → 0, therefore its leading behavior in this limit can be evaluated by numerical integration: Numerical values of this and further constants i µ are plotted as a function of µ in Figure A1. With a further increase in µ, first, the logarithmic singularity appears at µ = 3. It can be singled out, leading to: where i 3 = 0.64525. For µ > 3, to single out the leading singularities of the function under integration at small x, we integrate twice by parts, resulting in: Further analysis depends on the value of µ. In the region 3 < µ < 5, the integral on the right-hand side of Equation (A7) converges at ε → 0 and its leading asymptotics can be evaluated numerically. So, keeping the leading behavior of the first three terms in (A6) results in: with Logarithmic singularity appears in (A7) at µ = 5, leading to: with i 5 = −0.11309. For higher values of µ, analysis can be performed in a similar fashion. In particular, for 5 < µ < 7, again integrating twice by parts, one extracts a power-law singularity from the integral (A7): with Summarizing the above derived expressions for the leading behavior of the integral (A1) at small ε, we obtain the following useful formula: 3 < µ < 5 , ε −2 /4 − (ln ε)/12 + O(ε), µ = 5 , ε 3−µ /(2(µ − 3)) − ε 5−µ /(12(µ − 5)) + O(ε), 5 < µ < 7 .
(A13) Constants i µ for different µ can be evaluated numerically using formulas (A4), (A9) and (A12). Their dependence on µ is shown in Figure A1. They can be also checked against analogous constants evaluated in Ref. [45] using different integral representations.