pQCD running couplings finite and monotonic in the infrared: when do they reflect the holomorphic properties of spacelike observables?

We investigate a large class of perturbative QCD (pQCD) renormalization schemes whose beta functions $\beta(a)$ are meromorphic functions of the running coupling and give finite positive value of the coupling $a(Q^2)$ in the infrared regime (``freezing''), $a(Q^2) \to a_0$ for $Q^2 \to 0$. Such couplings automatically have no singularities on the positive axis of the squared momenta $Q^2$ ($ \equiv -q^2$). Explicit integration of the renormalization group equation (RGE) leads to the implicit (inverted) solution for the coupling, of the form $\ln (Q^2/Q^2_{\rm in}) = {\cal H}(a)$. An analysis of this solution leads us to an algebraic algorithm for the search of the Landau singularities of $a(Q^2)$ on the first Riemann sheet of the complex $Q^2$-plane, i.e., poles and branching points (with cuts) outside the negative semiaxis. We present specific representative examples of the use of such algorithm, and compare the found Landau singularities with those seen after the 2-dimensional numerical integration of the RGE in the entire first Riemann sheet, where the latter approach is numerically demanding and may not always be precise. The specific examples suggest that the presented algebraic approach is useful to find out whether the running pQCD coupling has Landau singularities and, if yes, where precisely these singularities are.


I. INTRODUCTION
According to general principles of Quantum Field Theories, the physical amplitudes D(Q 2 ) (such as the quark current correlators) and even unphysical amplitudes (such as the dressing functions of quark and transverse gluon propagators in QCD) are holomorphic (analytic) functions in the complex Q 2 -plane (where Q 2 ≡ −q 2 = −(q 0 ) 2 + q 2 ) except on the negative Q 2 semiaxis [1,2]. On the other hand, QCD running coupling a(Q 2 ) ≡ α s (Q 2 )/π can be defined as a product of the gluon dressing function and the square of the ghost dressing function in the Landau gauge [3]. Further, the leading-twist part of the spacelike physical QCD amplitudes D(Q 2 ) is a function of the running coupling, D(Q 2 ) (l.t.) = F(a(κQ 2 ); κ) (where κ ∼ 1 is a positive renormalization scale parameter). Therefore, a natural consequence of the holomorphic behaviour of QCD amplitudes D(Q 2 ) is that QCD running coupling a(Q 2 ) should reflect these properties, i.e., that a(Q 2 ) should be a holomorphic function in the complex Q 2 -plane with the exception of a negative semiaxis, Q 2 ∈ C\(−∞, −M 2 thr ], where M thr is a threshold mass, 0 ≤ M 2 thr 0.1 GeV 2 . The perturbative QCD (pQCD) frameworks usually used in the literature are the MS-type mass independent renormalization schemes (such as MS, 't Hooft, Lambert schemes), which give the running coupling a(Q 2 ) which is not holomorphic in the mentioned sense, but has a (Landau) branching point Q 2 * > 0 (∼ 0.1-1 GeV 2 ) for the cut, i.e., the cut reaches beyond the negative semiaxis to the positive IR regime, i.e., there is a Landau ghost cut (0, Q 2 * ). Further, the coupling often diverges at the branching point, a(Q 2 * ) = ∞ (Landau pole). These properties are mathematical consequences of the form of the beta function β(a) which appears in the RGE determining the flow of a(Q 2 ) with the squared momentum Q 2 . These properties contradict the earlier mentioned holomorphic properties for a(Q 2 ) which are motivated physically.
In our work, the considered QCD coupling a(Q 2 ) ≡ α s (Q 2 )/π will be such that it has so called freezing in the infrared regime, i.e., a(0) = a 0 is finite positive. This behaviour is suggested by the scaling solutions for the gluon and ghost propagators in the Landau gauge in the Dyson-Schwinger equations (DSE) approach [3,4], in the functional renormalization group (FRG) approach [5], stochastic quantization [6], and by . Further, 0 < a(0) ≡ a 0 < +∞ is also obtained in various physically motivated models for the running QCD coupling, among them the minimal analyticity dispersive approach [8][9][10][11] and its modifications or extensions [12], 1 and the AdS/CFT correspondence modified by a dilaton backgound [15]. For reviews, we refer to [16,17]. Such a behaviour has also been suggested in [18], where the running coupling definition involves explicitly the dynamical gluon mass and thus gives positive (nonzero) a(0) even in the case of so called decoupling solution of DSE [19] for gluon and ghost propagator in the Landau gauge. 2 All these approaches lead to nonperturbative (NP) running coupling A(Q 2 ), which in general differs from the underlying perturbative coupling a(Q 2 ) (i.e., the pQCD coupling in the same renormalization scheme) by power terms ∼ 1/(Q 2 ) N , i.e., terms of the type exp(−C/a) which cannot be Taylor-expanded around the pQCD point a = 0.
In Sec. II we define the class of considered pQCD beta functions β(a) (i.e., renormalization schemes), which are meromorphic functions leading to a finite positive a(0), and present the implicit solution of the RGE in the complex Q 2 -plane. In Sec. III we then present a practical algebraic procedure which allows us to find for a chosen beta function (in the considered class) the Landau singularities in the complex Q 2 -plane, i.e., the points where the behaviour of the running coupling a(Q 2 ) does not reflect the holomorphic properties of the spacelike Green functions D(Q 2 ) as required by the general principles of Quantum Field Theories. In Sec. IV we present some practical examples, and check with (2-dimensional) numerical integration of the RGE in the complex Q 2 -plane that the algebraic procedure gives us the correct answer. In Sec. V we summarize our results.

II. IMPLICIT SOLUTION OF THE RENORMALIZATION GROUP EQUATION
The renormalization group equation (RGE) for the coupling parameter F (z) ≡ a(Q 2 ) ≡ α s (Q 2 )/π, where z ≡ ln(Q 2 /Q 2 in ) is in general complex (and the initial scale is Q 2 in > 0), can be written in the following way: where the beta function β(F ) characterizes a mass independent renormalization scheme in perturbative QCD (pQCD), i.e., it has a well defined expansion around F = 0 Here, β 0 and β 1 are universal constants, β 0 = (11 − 2n f /3)/4 and β 1 = (102 − 38n f /3)/16, where n f is the number of active quark flavours. In the low-momentum regime (|Q 2 | 10 1 GeV 2 ), this number is usually taken to be n f = 3, corresponding to the three lightest, almost massless, active quarks u, d and s. The coefficients β j (j ≥ 2) characterize the pQCD renormalization scheme [28]. As mentioned in the Introduction, there exist several theoretical arguments which suggest that the running coupling F (z) ≡ a(Q 2 ) is a finite function for all positive Q 2 and that it possibly acquires a finite positive value in the infared limit, 0 < a(0) ≡ a 0 < +∞. In this case, it turns out that β(F ) for positive couplings F ≡ a has a root at F = a 0 [and double root at F = 0 according to Eq. (2)], and for 0 < F < a 0 it has no roots. In view of this, we will consider the following class of beta functions: where T M (Y ) and U N (Y ) are polynomials of degree M and N , respectively, both normalized in such a way that T M (0) = 1 = U N (0). Specifically, we denote as 1/t j the roots of T M (Y ), and 1/u j the roots of U N (Y ) We will restrict ourselves, for physical reasons, to such beta functions of the form (3) in which those t j and u k which are real and positive are all below unity: 0 < t j < 1 and 0 < u k < 1. This means that: • a = a 0 is the smallest positive root of the beta function; • and that all those poles of the beta function which are positive are larger than a 0 .
If the latter conditions were not fulfilled, the running coupling a(Q 2 ) would obviously have (Landau) singularities on the positive Q 2 -axis, contradicting the theoretical arguments mentioned in the Introduction. The former condition only means that we define a 0 as the smallest positive root of the beta function, and demand that at least one such positive root exist. An important practical consequence of these restrictions will be highlighted in Sec. III (the first paragraph). The first universal coefficient β 0 in the expansion of the beta function (2) is reproduced automatically by our construction. The second universal coefficient β 1 in Eq. (2) imposes the following restriction on the polynomials T M (Y ) and U N (Y ): where are the coefficients at Y 1 of the two polynomials T M (Y ) and U N (Y ). In addition, we will restrict the considered class of meromorphic beta functions to M + 1 ≥ N . In such a case, it turns out that the RGE (1) can be integrated algebraically and leads to an implicit solution of the form z = G(F ) [for F ≡ F (z)]. Namely, the integration of the RGE (1) gives and if we introduce a new integration variable t ≡ a 0 /F , this can be written as where a in = a(Q 2 in ) = F (z = 0) has a real positive value, 0 < a in < a 0 . When M + 1 ≥ N , the integrand can be written as a sum of simple partial fractions 1/(t − t j ), where t 0 = 1 and t j (j = 1, . . . , M ) are the roots of the (M -degree) polynomial t M T M (1/t) Namely, we have where the M + 1 constants B j are As a special case, we see that which is a real number. Using this, and the expression (3), we also obtain the following relation: Incidentally, in the limit of large t the relations (10) imply the following sum rule: where the last equality is obtained by using the relation (5). Using the form (10b) for the integrand in Eq. (8) leads us immediately to the implicit solution of the RGE Each logarithm has an ambiguity (winding number) because where ln (pb) is the principal branch: −π < Im ln (pb) A = Arg(A) ≤ +π; further, n A is the winding number representing the ambiguity. When A is positive, we consider that ln A is automatically the principal branch. This would then suggest that the right-hand side of Eq. (15) has M + 1 independent winding numbers n j , correspondig to each logarithm there. The physically acceptable winding numbers of the logarithms on the right-hand side of Eq. (15) are such that they give for the expression (15) a number z = ln(Q 2 /Q 2 in ) corresponding to the squared impulse Q 2 on the first Riemann sheet, i.e., |Imz| ≤ π (cf. also the discussion in the beginning of Sec. III A).
However, in general some of the roots t j of the polynomial T M (Y ), Eq. (4a), are not real, but form complex conjugate pairs. For example, if the first complex conjugate pair is (t 1 , t 2 = t * 1 ), then it is straightforward to check that the corresponding coefficients B 1 and B 2 are mutually complex conjugate, and the corresponding two terms in the sum (10b) are and the coresponding contribution to the integral (8) is Therefore, the expression on the right-hand side of Eq. (15) can be rewritten more explicitly, for the case when t j (j = 1, . . . , 2P ) are P complex conjugate pairs and t j (j = 2P + 1, . . . , M ) are real We note that among the P z-dependent ArcTan terms, which are in general complex, each has one winding number because for A = |A| exp(iθ) (|θ| ≤ π) 3 where we regard as the principal branch ArcTan (pb) (A) the one which fulfills the inequality −π/2 < ReArcTan (pb) (A) ≤ +π/2. When A is real, we consider that ArcTan is automatically the principal branch. Further, each of the M − P + 1 z-dependent logarithms appearing in Eq. (19) has a winding number according to the relation (16). This means that we have in general in total M + 1 winding numbers. This realization will play a role in the next Section III.
We recall that the considered β(a) functions are such that a(Q 2 ) is a holomorphic function in and around any positive point Q 2 > 0. However, at Q 2 = 0, where a = a 0 < ∞, the function a(Q 2 ) could be nonholomorphic (nonanalytic), i.e., certain (high enough) derivative (d/dQ 2 ) n a(Q 2 ) at Q 2 = 0 could be infinite. In our considered cases we have for the Taylor expansion around This implies The use of relation (13) then gives the power index κ in terms of the parameters contained in the considered beta function Eq. (3) In general, κ is noninteger, and consequently the coupling is in general not analytic at Q 2 = 0 (z = −∞).
We wish to point out that the class of the β-functions considered here, Eq. (3), in addition to having a Padé form P [M + 3/N ](a), have specific restrictions which result in a finite positive and monotonically decreasing running coupling a(Q 2 ) < a 0 on the entire nonnegative . This is reflected in the formal requirement that those of the parameters t j and u k of Eqs. (4) which are not complex and are positive must fulfill the restrictions 0 < u k < 1 and 0 < t j < 1.

III. SINGULARITIES (LANDAU) OUTSIDE THE REAL Q 2 -AXIS
We note that, by restrictions on the beta function mentioned in the previous Section, the running coupling a(Q 2 ) ≡ F (z) has no singularities on the real positive Q 2 -axis, i.e., there are no positive-Q 2 Landau singularities. This is so because a(Q 2 ) is constrained there to run between the value a(0) = a 0 (> 0) and a(+∞) = 0, the latter equality being valid by the asymptotic freedom of QCD reflected in the form (2) of beta function when F → 0. Namely, by our restrictions on the class of considered β(a) functions, when a(Q 2 ) is RGE-running with increasing positive Q 2 , beta function β(a(Q 2 )) will be negative finite all the time, since no new roots or poles of the beta function are encountered. Therefore, by the mentioned restrictions on the roots and poles of the beta function (3) we ensured in advance that the positive-Q 2 Landau singularities (poles and/or cuts) do not exist.

A. Landau poles
We will now construct an algebraic algorithm which allows us to verify whether in the (first Riemann sheet of the) complex Q 2 -plane the solution (15) has poles outside the real Q 2 -axis (Landau poles). We assume that only the first sheet of the complex Q 2 -plane has physical meaning, 7 i.e., Q 2 = |Q 2 | exp(iφ) where −π ≤ φ < +π. This corresponds for the z ≡ ln(Q 2 /Q 2 in ) variable to be a band in the complex z-plane with −π ≤ Imz < +π, cf. Figs. 1 (a), (b). As argued in the Introduction, if a(Q 2 ) is to reflect the holomorphic properties of spacelike Green functions and observables, such as current correlators or structure functions, then a(Q 2 ) can have singularities (cut) only along the negative Q 2 -axis: −∞ < Q 2 ≤ −M 2 thr , where the threshold mass M 2 thr is either positive (∼ 0.1 GeV 2 ) or zero. This cut corresponds in the z-stripe to the cut along the Imz = −π border line.
As explained, the possible Landau singularities in the considered pQCD renormalization schemes are within the Q 2 -complex plane outside the real Q 2 -axis. In the z-plane this corresponds to the possible Landau singularities within the interior of the z-stripe, −π < z < +π, but not along the real axis, z ∈ R.
Landau pole z is usually a branching point of a cut singularity of F (z), such that F (z * ) = ∞, and it is situated on the first Riemann sheet outside the timelike semiaxis (|Imz * | < π). Let us denote a * ≡ F (x * ) (0 < a * < a 0 ). We then apply the implicit solution Eq. (19) to the points z 1 = x * and z 2 = x * + iy * , and subtract the two equations; this then gives us the equation where Here, we accounted for the nonuniqueness of the (z-dependent) logarithms Eq. (16) and ArcTan Eq. (20) lim y→y * where (k = 0, . . . , P − 1), and we denoted the M + 1 winding numbers where n j , N k , N k = 0, ±1, ±2, . . .. We note that the terms ln (pb) (−t j ) may have t j either negative or positive, and therefore As a special case, we used in Eq. (25): ln(−t 0 ) = ln(−1) = iπ. We note that, since the real roots t j fulfill the inequality t j ≤ 1 [our initial physical restrictions on β-function, cf. comments after Eqs. (4)], the logarithm ln(a 0 /a * − t j ) in Eq. (25) is a real number because it has positive argument. For the same reason, also the P logarithms of the trinomials in (a 0 /a * ) in Eq. (25) are real. We point out that that winding numbers n appear when integrating the RGE (1) from z 1 = x * to z 2 = x * + iy * .
If, for example, all B j coefficients are real, then Im G n (a) = Im G 0 (a); if in such a case Im G 0 (a * ) has no zero in the positive interval 0 < a * < a 0 , then one necessary condition for the existence of Landau poles is not fulfilled, i.e., there are no Landau poles.

B. Landau branching points
In the previous Subsection we presented an algorithm which allows us to find, inside the complex z-stripe, the (Landau) poles where the coupling is infinite F (z * ) = ∞. However, the complex function F (z) can have also another type of Landau singularities, namely a cut with a finite-valued branching point z * .
One illustrative mathematical example is F (z) = (z − z * ) 1/2 , where z * = x * + iy * is such a branching point, F (z * ) = 0 and F (z) = ∞. The cut in this case is usually defined along the semiaxis to the left of z * : x + iy * (x ≤ x * ).
However, we may worry at first that other, even more "finite," type of Landau branching points z * ∈ R may appear, such as F (z) = (z − z * ) 3/2 , for which F (z * ) < ∞ and F (z * ) = ∞. We show that this does not occur for the considered class of meromorphic beta functions (3)-(4). Namely, The poles of the right-hand side are at the same values F = a 0 /u s (s = 1, . . . , N ) as in the beta function β(F ) itself, cf. Eqs. (3)- (4). This means that, if F (z * ) = ∞, then F (z * ) = ∞. We can continue this argumentation, by applying further derivatives (d/dz) n to Eq. (31). E.g., if F (3) (z * ) = ∞, then F (z * ) = ∞. Therefore, the only relevant situation of finite-valued Landau branching points z * for the considered beta functions is: F (z * ) = ∞ and F (z * ) < ∞. Since F (z * ) = β(F (z * )), such a branching point is one of the poles of the beta function, z where a The winding numbers are generated in a limiting process analogous to that in Eqs. (26) lim when integrating the RGE (1) in the z-plane along the vertical direction, from z 1 = x (s) * toward z 2 = x (s) * + iy (s) * . We note that one of the physically motivated restrictions on the β-function, from the outset, was that those roots u s which are real satisfy u s < 1 [cf. the comments after Eqs. (4)]. This means that for such u s , the branching point z (s) * where F (z * ) = a 0 /u s (> a 0 ) cannot be achieved at real z (s) * = x (s) * , i.e., also in such cases z (s) * must have y (s) * = 0, and thus we can have also in such a case nonzero winding numbers n = 0.
Here, the procedure described in the previous Sec. III A for Im G n (a * ) and Re G n (a * ) [for a * ≡ F (x * ) in the interval 0 < a * < a 0 , and for n j , N k , N k = 0, ±1, . . . ,], is now performed for Im K n (a (s) * ; u s ) and Re K n (a where ImK n (a and ReK n (a We notice that ImK n (a (s) * ; u s ) depends only on the winding numbers N (s) k , cf. Eq. (33). The existence of a Landau branching point means that Eqs. (36) have a solution, for an s, and a (s) * and y (s) * such that: 0 < a (s) * < a 0 and |y (s) * | < π. The procedures described in this Section III for Im G n (a * ) and Re G n (a * ), and for Im K n (a (s) * ; u s ) and Re K n (a (s) * ; u s ), represent a relatively simple algebraic instrument for practical verification of whether the pQCD scheme with a given beta function of the form (3) described in Sec. II has Landau singularities or has no such singularities, and where these singularities are.

A. Polynomial β with real roots
Here we consider the case of (M = 2, N = 0) where t 1 and t 2 are real [and t 1 , t 2 < 1 by physical requirements, cf. the text after Eqs. (4)]. We choose the input values The conditions (5)-(6) then give where the numerical value is obtained by using in the universal β-coefficients β 0 and β 1 the number of active quark flavours N f = 3. This then immediately gives for the κ coefficients of Eqs. (21)-(23) the value Since t 1 and t 2 are real (hence: M = 2; P = 0), the only winding numbers (27) are n = (n 0 , n 1 , n 2 ), and thus ImG n (a * ) is independent of n. The first condition of Eq. (29) then immediately gives for a * (we recall: 0 < a * < a 0 ) ImG n (a * ) = 0 ⇒ a * ≈ 0.315656, (43) and the corresponding x * is 8 as can be easily checked by the implicit solution (19) when using there for F (z) the value of a * Eq. (43). When we now numerically integrate the RGE (1) along the line Re(z) = x * in the z-plane, 9 we obtain for the real and imaginary 8 We use throughout the reference value αs(M 2 Z ; MS) = 0.1179 [45]. This corresponds to the N f = 3 regime at Q 2 in = (2mc) 2 = 2.54 2 GeV 2 to a(Q 2 in ; MS) (≡ αs(Q 2 in ; MS)/π) = 0.0834921. We will use this reference value throughout (although, by using a different reference value is equivalent to changing the value of Q 2 in which does not affect our conclusions). The RGE-running from M 2 Z down to (2mc) 2 in MS is performed by using the five-loop RGE [46] with four-loop quark threshold conditions at µ 2 thr. = (2mq) 2 [47], where the MS quark mass values formq ≡mq(m 2 q ) was takenm b = 4.20 GeV andmc = 1.27 GeV. The transition from the (five-loop) MS scheme to the scheme of the considered β-function was performed at the scale Q 2 = (2mc) 2 and N f = 3, according to the approach as explained, e.g., in Ref. [26] [Eq. (13) there]. This gives, in the considered scheme of the β-function (39), the value a(Q 2 in ) = 0.0737597.
The obtained points z (j) * ,± are the Landau poles. We can cross-check that these points are really the Landau poles by evaluating the algebraic expression G n (a * ), Eq. (25), for various winding numbers n = (n 0 , n 1 , n 2 ), and we find that G n (a * ) = +y On the other hand, without the algebraic seminumeric approach described above, it would be difficult to find the (four) Landau poles on the first Riemann sheet of Q 2 . In Figs. 3(a),(b) we present |β(F (z))| for the considered couplings, obtained by the 2-dimensional numerical integration of the RGE in the complex z-stripe. In Fig. 3(a) it is difficult to see that there are two Landau poles close to each other, at positive (and negative) values of Im(z) = y; only the strongly "zoomed" Fig. 3(b) suggests that there are two poles near to each other, as clearly obtained in Eqs. (46) by our algebraic seminumeric analysis.
The conditions (5)-(6) then give where, as in Sec. IV A, the numerical value is obtained by using in the universal β-coefficients β 0 and β 1 with N f = 3. This then gives for the κ coefficient of Eqs.  3: (a) The numerical values of |β(F (z))| in the the physical z-stripe (−π ≤ y ≡ Imz < π), corresponding to the first Riemann sheet of the complex momenta Q 2 . The numerical results indicate only one Landau pole in this region, and its complex conjugate. (b) "Zoomed" numerical calculation indicates two mutually close Landau poles in this region (and their complex conjugates). The calculation was performed using Mathematica software [48].
Since t 1 and t 2 are complex nonreal (hence: M = 2; P = 1), the winding numbers (27) are n = (n 0 , N 0 , N 0 ), and thus ImG n (a * ) depends on N 0 and ReG n (a * ) depends on n 0 and N 0 . The first condition of Eq. (29) then gives for a * the acceptable solution (i.e., in the interval 0 < a * < a 0 ) only when As in Sec. IV A, we conclude that the obtained points z (j) * ,± are the Landau poles. We cross-check that these points are really the Landau poles by evaluating the algebraic expression G n (a * ), Eq. (25), for various values of the winding numbers n = (n 0 , N 0 , N 0 ), and we find G n (a * ) = +y On the other hand, the fully numerical (2-dimensional) integration of the RGE (1) in the first Riemann sheet of the complex squared momenta Q 2 (i.e., in the complex z-stripe with |Imz| ≤ π) gives us the results in Figs. 5(a),(b) where we present |β(F (z))| for the considered couplings. In Fig. 5(a) it is hard to see two of the four mentioned Landau poles, namely those with Im(z) = ±2.02315. Only the strongly "zoomed" Fig. 5(b) suggests that there are Landau poles also at z = x * ± i 2.02315.

C. Padé β with a real pole
Here we consider a numerical example for the types of β-function of Sec. III B where a finite-valued Landau branching point is realized. We will take the simplest case M = 1 and N = 1 (and P = 0) where β-function Eq. (3) has a Padé form with one real pole Here, both u 1 and t 1 are real and related via the relation (5). The β-function has a pole at the coupling value F (z * ) = a 0 /u 1 . We choose the input values The conditions (5)-(6) then give where, as in the previous examples, we use the values of β 0 and β 1 with N f = 3. This then gives for the κ coefficient the value Since in the considered case we have M = 1 and P = 0 (and N = 1), the only winding numbers are n = {n (1) 0 , n 1 }. The (real) value of the coupling F (x (1) * ) = a (1) * (0 < a (1) * < a 0 ) is then obtained by the condition ImK n (a (1) * , u 1 ) = 0, cf. Eq. (36), where ImK n (a (1) * , u 1 ) is independent of the winding numbers n ≡ {n (1) 0 , n (1) 1 }. This then immediately gives us The corresponding value of x (1) * is [we use α s (M 2 Z ; MS) = 0.1179 as in Sec. IV A] is obtained from the implicit solution (19) Now performing the simple (1-dimensional) numerical integration of the RGE (1) On the other hand, the second condition in Eq. (36) should give us in this case the same values ±y (1) * = ±1.14448. Indeed, the evaluation of the algebraic expression K n (a (1) * , u 1 ), Eq. (33), gives This is consistent with the results (61), and clearly shows that in the considered case the Landau branching point is achieved in the first Riemann sheet only at the two complex conjugate points z 1 } equal to {0, 0, } and {−1, 0}, respectively. Further, Figs. (6) indicate that the coupling F (z (1) * ) at this point achieves the (real) value 0.60 which coincides with the value a 0 /u 1 , i.e., the value where β-function diverges (but not the coupling).
The fully numerical (two-dimensional) integration of the RGE (1) in the first Riemann sheet of the complex squared momenta gives us the results in Figs. 7. Figure 7(a) shows |β(F (z))| and indicates the Landau singularities at z (1) * ,± = x (1) * ± iy (1) * . Figure 7(b) shows ImF (z) and indicates that the previously mentioned singularities are indeed branching points, with the cut in the complex-z stripe extending from z    On the other hand, the algebraic seminumeric analysis above, Eqs. (59)-(62) and Figs. 6, shows that the Landau singularities z (1) * ,± = x (1) * ± iy (1) * are indeed branching points (with cuts) and correspond to specific winding numbers, and that no other branching points exist in the first Riemann sheet.

V. SUMMARY
In this work we presented an algebraic algorithm for finding possible Landau singularities of the pQCD running coupling a(Q 2 ) in the complex plane of the squared momenta Q 2 (first Riemann sheet). We considered a large class of β-functions, representative of the scenarios where the running coupling a(Q 2 ) is a monotonic function of Q 2 at positive Q 2 and "freezes" in the IR sector, a(Q 2 ) → a 0 for Q 2 → 0, where the IR freezing value a 0 is considered positive finite. The consideration of the running coupling a(Q 2 ) ≡ F (z) was performed on the corresponding complex z-stripe, −π ≤ Im(z) < π, where z = ln(Q 2 /Q 2 in ) and Q 2 in > 0 was an initial scale for the integration of the RGE. The analysis was performed by explicit integration of the RGE which led to the implicit (inverted) solution of the form z = H(F ). An analysis of this implicit solution than led us to an algebraic procedure for the search of the Landau singularities of F (z) on the z-stripe. We considered two types of such singularities, the poles F (z) = ∞ and the branching points (for cuts) β(F (z)) = ∞. For illustration, we then presented the mentioned algebraic (and seminumeric) analysis for three specific representative cases of the β-function, and compared the found Landau singularities with those seen directly by the numerical 2-dimensional integration of the RGE in the entire complex z-stripe, the latter approach being numerically demanding. The presented specific cases suggest that our algebraic seminumeric approach is reliable and has high precision in finding the Landau singularities, while the 2-dimensional integration of the RGE gives these singularities with less precision and sometimes we may miss some of the singular points with this purely numerical method, especially if the numerical scanning over the entire z-stripe is made with limited density. Therefore, the presented algebraic seminumeric formalism appears to be useful when we want to find out whether the pQCD running coupling has Landau singularities, and if there are any, to find the location of these singularities with high precision.
sheet in the squared momentum plane Q 2 (= −q 2 = Q 2 in exp(z))] was always performed first along the entire real z axis, and then at each fixed real value of z = x the RGE was integrated along the imaginary (y) direction of z = x + iy (−π ≤ y < +π).