Gluon bound state and asymptotic freedom derived from the Bethe–Salpeter equation

... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . In this paper we study the two-body bound states for gluons and ghosts in a massive Yang– Mills theory which is obtained by generalizing the ordinary massless Yang–Mills theory in a manifestly Lorentz-covariant gauge. First, we give a systematic derivation of the coupled Bethe– Salpeter equations for gluons and ghosts by using the Cornwall–Jackiw–Tomboulis effective action of the composite operators within the framework of the path integral quantization.Then, we obtain the numerical solutions for the Bethe–Salpeter amplitude representing the simultaneous bound states of gluons and ghosts by solving the homogeneous Bethe–Salpeter equation in the ladder approximation. We study how the inclusion of ghosts affects the two-gluon bound states. Moreover, we show explicitly that the approximate solutions obtained for the gluon– gluon amplitude are consistent with the ultraviolet asymptotic freedom signaled by the negative β function. .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Subject Index B69


Introduction
As is well known, the basic ingredients relevant to the strong interactions are quarks and gluons as described by quantum chromodynamics (QCD). However, we believe that they are confined and never observed in their isolated form, and that only the color-singlet composite particles such as mesons, baryons, and glueballs consisted of quarks and gluons can be observed in experiments. In particular, glueballs were predicted in the middle of the 1970s by Fritzsch and Minkowski [1], in contrast to mesons and baryons being explained by the quark model proposed by Gell-Mann in the middle of the 1960s [2]. Since then, the properties of glueballs have been investigated in various theoretical frameworks and also by experiments.
The equation for describing the bound states in the framework of the relativistic quantum field theory is known as the Nambu-Bethe-Salpeter equation or the Bethe-Salpeter (BS) equation [3,4]. The BS equation is a self-consistent equation written in terms of the BS amplitude as the wave function of the bound state. See, e.g., Refs. [5,6] for reviews. The BS equation is specified for a given set of propagators and the interaction vertex functions describing the constituent particles. In the non-perturbative framework, propagators and vertex functions can be obtained by solving the Schwinger-Dyson (SD) equation, which is a self-consistent equation written in terms of propagators and vertex functions for the fundamental fields in a given field theory. In order to describe the bound state in the non-perturbative framework, therefore, we must take into account the SD equation and the BS equation simultaneously: The SD equation must be solved to give propagators and vertex functions as inputs for specifying the BS equation. Then the BS equation is to be solved to the same PTEP 2017, 053B05 H. Fukamachi et al. level of approximation as those used in solving the SD equation, as far as they cannot be solved exactly.
The purpose of this paper is to understand glueball formation from the underlying dynamics of the Yang-Mills theory in a way consistent with confinement and ultraviolet asymptotic freedom. In order to avoid complications arising from the mixing of glueballs with mesons in QCD, we restrict our consideration to the pure glue dynamics, namely, pure Yang-Mills theory. Therefore, a glueball is supposed to be a bound state of gluons. In this paper, we study a glueball as a two-gluon bound state to simplify the treatment.
In the late 1970s, Fukuda [7] studied a two-body bound state of gluons and ghosts as a solution of the BS equation in the ladder approximation in the framework of the Yang-Mills theory with the covariant gauge fixing. He showed the existence of the tachyon bound state due to the dominating attractive force between two gauge (gluon) particles in the color-singlet and spin-singlet channels, for whatever small coupling constant. The gauge particles interact among themselves through a three-gluon interaction and a contact four-gluon interaction, and interact with ghosts through the gluon-ghost interaction due to the covariant gauge fixing. The three-gluon interaction gives an attractive force between two gauge particles in the color-singlet channel by exchanging a gauge particle (this is the analog of the attraction of opposite charge in the Abelian gauge theory). The contact four-gluon interaction produces a repulsive force. However, if the attraction due to the former dominates the repulsion from the latter, the two gauge particles can form a bound state by the net attractive force. This is nothing but a pairing phenomenon familiar in the BCS theory of superconductivity. Since gauge particles are massless, however, the resulting bound state should necessarily correspond to a tachyon pole. However, Fukuda's observations play the very important role of obtaining glueballs consisting of massive gluons in this paper. It was also shown [8] that such a tachyonic bound state continues to exist even in the presence of n flavors of quarks, and gluons acquire a mass due to vacuum rearrangement, when n is smaller than or equal to a critical value n c , n ≤ n c , while all particles remain massless for n > n c .
At the classical level, indeed, the Yang-Mills theory in four-dimensional spacetime is a massless theory with no scale. At the quantum level, however, the Yang-Mills theory should have a mass scale as suggested from the dimensional transmutation, and could have the mass gap. Recently, a solution corresponding to a confining gluon of the massive type was discovered as a solution of the coupled SD equation for the gluon and ghost propagators in the SU (N ) Yang-Mills theory with the Landau gauge fixing, which is called the decoupling solution [9][10][11][12][13], in sharp contrast to the scaling solution [14][15][16][17]. The existence of the decoupling solution is supported, and its confining properties are investigated by various analytical and numerical frameworks. However, it is not yet clear how to understand the decoupling solution as a solution to be consistent with gluon confinement and more general color confinement; see, e.g., Ref. [18].
The Brout-Englert-Higgs mechanism is accepted as a unique way for providing the mass for gauge bosons in quantum field theories, because it is the only established method that enables one to maintain both renormalizability and physical unitarity. In this paper, we consider a specific model, including a mass term for the Yang-Mills field, which is a special case of the Curci-Ferrari model [19]. We call this model the massive Yang-Mills theory for later convenience. We adopt this model to perform an analytical investigation for gluons of massive type, instead of using the decoupling solution obtained by solving the SD equation, since the decoupling solution has so far been obtained 2/32 Downloaded from https://academic.oup.com/ptep/article-abstract/2017/5/053B05/3835935/Gluon-bound-state-and-asymptotic-freedom-derived by CERN -European Organization for Nuclear Research user on 03 October 2017 only in a numerical way, and the analytical expression for the decoupling solution is not yet known explicitly.
We started our analysis on massive Yang-Mills theory in the previous papers [20,21] without introducing the Higgs field, which causes the Higgs mechanism. Nevertheless, we are not plagued by the unitarity violation in the high-energy region of the massive Yang-Mills theory without the Higgs scalar. This is because we regard our massiveYang-Mills theory just as a low-energy effective theory (valid below a certain cutoff ) of the quantum Yang-Mills theory with mass gap and confinement. To mimic precisely the decoupling solution, we must introduce the momentum-dependent gluon mass m( p) which vanishes in the high-energy region and remains non-zero in the low-energy region. In this paper, we have adopted a constant mass m valid up to a momentum cutoff to simplify the analyses without assuming the specific function of the momentum p. Our standpoint mentioned above has already been given more explicitly in Ref. [21], and is not repeated here.
It is also known that this model of the massive Yang-Mills theory is still renormalizable but loses the physical unitarity, at least in the perturbation theory [20,21]. We leave the issue of how to realize it as a sound model in a quantum field theoretical framework. We have a hope that this issue will be resolved in a non-perturbative framework. One of the motivations of this paper is to obtain some information toward this direction to consider bound states in a massive Yang-Mills theory without the Higgs field as a low-energy effective theory of QCD.
The issue of a tachyon bound state will be avoided by starting from the massive Yang-Mills particles instead of the massless gauge particles. Consequently, we will obtain the massive bound states consisting of two massive gluons described by the massive Yang-Mills theory. In obtaining such a solution, we pay special attention to how our results for the bound state are consistent with the ultraviolet asymptotic freedom [22]. This is a feature that has not been discussed in the preceding works [23,24]. In their works, the precise determinations of the glueball spectrum have been the main purpose of the investigations. Such an issue is not the aim of the present paper, although it is a nice goal for future investigations to be achieved by developing our method.
We give a systematic derivation of the coupled BS equation for gluon and ghost by starting from the the Cornwall-Jackiw-Toumboulis (CJT) effective action for the composite operator with the bilocal source term. We show that the bound state obtained as a solution of the homogeneous BS equation for the gluon-gluon BS amplitude is consistent with the ultraviolet asymptotic freedom. This paper is organized as follows. In Sect. 2, we introduce the massive Yang-Mills theory from which we start the study of the bound state. In Sect. 3, we give the CJT effective action for the bilocal composite operator, which gives a systematic derivation of the BS equation. In Sect. 4, we derive the BS equation by differentiating the CJT effective action with respect to the full Green functions. In Sect. 5, we write down the coupled BS equation for gluon and ghost at vanishing total momentum P = 0 in the massless case, and obtain the numerical solutions. This is a preliminary study of obtaining the bound state. In Sect. 6, we write down the homogeneous BS equation at the general total momentum P in the massive case, but restricted to the gluon-gluon BS amplitude alone. In Sect. 7, we solve the BS equation for the gluon-gluon BS amplitude at the general total momentum P in the massive case. The final section is devoted to conclusion and discussion. Some technical materials are collected in the appendices. Lorenz gauge [25] plus the "mass term" L m formulated in a manifestly Lorentz-covariant way. The total Lagrangian density L tot mYM is written in terms of the Yang-Mills field A μ , the Faddeev-Popov ghost field C , the antighost fieldC , and the Nakanishi-Lautrup auxiliary field N , if we use the terminology adopted in the usual massless Yang-Mills theory: where D μ is the covariant derivative defined by andN is defined byN Here we have adopted the notation: where f abc are the structure constants of the Lie algebra of the SU (N ) group. Here, α and β are parameters that are identified with the gauge-fixing parameters in the m → 0 limit. The α = 0 case is the Curci-Ferrari model [19] with the coupling constant g, the mass parameter m, and the parameter β. In the Abelian limit with vanishing structure constants f abc = 0, the Faddeev-Popov ghosts decouple and the Curci-Ferrari model reduces to the Nakanishi model [26].
In what follows, we put β = 0 to simulate the decoupling solution representing massive gluon and massless ghost. By eliminating the Nakanishi-Lautrup field N , L tot mYM reads By starting from the massive Yang-Mills theory given by the Lagrangian density (5) using the CJT effective action for the composite operators, as explicitly shown in Ref. [22]. See Refs. [20,21] for more details on this model.

Cornwall-Jackiw-Tomboulis effective action
For the Yang-Mills theory, we introduce the local sources J μ a (x),η a (x), and η a (x) for the gluon field A a μ (x), the ghost field C a (x), and the antighost fieldC a (x), respectively. Moreover, we introduce the bilocal sources I μν ab (x, y) and θ ab (x, y) for the bilocal composite operators A a μ (x)A b ν (y) andC a (x)C b (y), respectively, where the bilocal source I is symmetric, I ab μν (x, y) = I ba νμ (y, x), and θ ab (x, y) is a general matrix. Then, the generating functional W [J ,η, η, I , θ] is defined by where we have introduced the abbreviated notation x := d D x. For the local source, we find that the (left) derivatives of W with respect to the local sources are given by where D and˜ are the full gluon and ghost-antighost propagators. For the bilocal source, we find that the derivatives of W with respect to the bilocal sources are given by We define the CJT effective action [27] of the Yang-Mills theory by the Legendre transform: Then we find the derivatives of with respect to the field expectation, and the derivatives of with respect to the full propagators, It is shown that the CJT effective action is rewritten as where 2 is the two-particle irreducible (2PI) contributions, D and˜ are respectively the full gluon and ghost propagators, and D AA andD CC are respectively the tree-level gluon and ghost propagators defined by The original Yang-Mills theory is characterized by the bare gluon and ghost propagators, and the bare one-particle irreducible (1PI) vertex functions: three-gluon, four-gluon, and gluon-ghostantighost vertices given in Fig. 1. The full gluon and ghost propagators and the full 1PI vertex functions, three-gluon, four-gluon, and gluon-ghost-antighost vertices, are given in Fig. 2. It is shown that the 2PI vacuum diagrams in the Yang-Mills theory are written as Fig. 3.

Derivation of the Schwinger-Dyson equation
By taking the first derivative of the CJT effective action with respect to the full gluon and ghost propagators D, we obtain the SD equations for the gluon propagator D and ghost propagator δ˜ respectively (after the sources are set to be zero, I , θ → 0):

The equation for the scattering Green function
The BS equation is obtained by taking the first derivative of the SD equation (i.e., the second derivative of the CJT effective action) with respect to the full propagator. For the general field φ, the scattering Green function D 4 with four external lines is defined by the second derivative of W with respect to the bilocal source I : while the inverse D −1 4 of the scattering Green function D 4 is given by the second derivative of with respect to the full propagator D: Therefore, the first derivative of the SD equation yields the equation for the inverse D −1 4 of the scattering Green function D 4 : where D 4 is the tree part of the scattering Green function D 4 and K denotes the kernel. The diagrammatic representation is given as   We make use of the identities: ½ and ´¼µ ´¼µ ½ (18) to rewrite the equation in terms of D 4 , not the inverse D −1 4 . We multiply D 4 from the right to obtain and then multiply D Thus we obtain the self-consistent equation for the scattering Green function D 4 after multiplying by i: For the Yang-Mills theory, we obtain a set of the BS equations for the four BS amplitudes: gluon-gluon, gluon-ghost, ghost-gluon, and ghost-ghost scattering Green functions as given in

The homogeneous BS equation for the (amputated) BS amplitude
We can use the completeness relation: where the sum c runs over all the bound states except for the continuum spectrum. The scattering Green function D 4 for particles described by the field φ a is defined by This is decomposed into the product of the BS amplitude χ: where the BS amplitude χ is defined by The momentum-space representation of D 4 has the pole: the productχχ of the BS amplitude corresponds to the residue at the pole where the squared momentum P 2 coincides with the bound state mass M squared: From the equation for the scattering Green function D 4 , we can extract the productχχ of the BS amplitude by identifying it with the pole residue at the pole corresponding to the bound state, since D (0) 4 does not have the pole: By multiplying D we finally obtain the homogeneous BS equation for the amputated BS amplitude D (0)−1χ : For the Yang-Mills theory, the scattering Green function D 4 for the gluon field is defined by This is decomposed into the product of the BS amplitude χ : In the momentum-space representation, we have We introduce the BS amplitudes χ ab P;μν and χ ab P for two-gluons and ghost-antighost respectively as Thus, the amputated BS amplitude obeys the homogeneous BS equation: where we have defined the momenta p ± by and K by The diagram of the homogeneous BS equation is given in Fig. 7. Here, the kernels of the BS equation are given as Fig. 8. The explicit form of the kernel is given in Appendix A. This BS equation can be compared with that obtained by Meyers and Swanson [23].

BS amplitude
In what follows, we consider the color-singlet and scalar (i.e., spin-zero) bound state. Fukuda considered the decomposition for the amplitude: where A 1 , . . . , A 5 are functions of k 2 , P 2 , and (k · P) 2 . Here, A 1 and A 2 are respectively the transverse and longitudinal parts in terms of k, since we can define the transverse projector t μν and longitudinal projector μν by with the properties:   In this paper, however, we regard the amplitude A μν as the function A μν (k + ; k − ) of k + := k + P/2 and k − := k −P/2, rather than the function A μν (k; P) of k and P. For this purpose, this decomposition in terms of k and P is not suited in the presence of the non-vanishing total momentum P μ . We construct another decomposition of the amplitude in terms of k + and k − . In particular, we require that the gluon BS amplitude A μν (k + , k − ) is transverse in the sense that In view of these, we can define the modified transverse projector in the presence of P by which is subject to Indeed, the modified transverse projector T μν reduces to the usual transverse projector (41) with the property (42d) in the limit P → 0.
In the presence of P, we can construct another type of transverse projector. For this purpose, we can introduce the vectors k +⊥ and k −⊥ which are orthogonal to k + and k − respectively: Indeed, such vectors k +⊥ and k −⊥ are constructed from k + and k − as 1 We can also define another transverse projector with the property Notice that k +⊥ μ and k −⊥ μ vanish in the limit P μ → 0. In the limit P → 0, therefore, L μν becomes illdefined and does not reduce to μν , which is longitudinal to k. They satisfy the following properties: which lead to Now we proceed to construct the transverse amplitude. First, k + μ A μν = 0 is satisfied by a tensor of the form: Second, the condition A μν k − ν = 0 yields 1 This is obtained as follows. We define the unit vectors e + and e − of k + and k − : and we can use the ortho-normalization method due to Gram-Schmidt to construct the vectors e +⊥ and e −⊥ which are respectively orthogonal to e + and e − :

13/32
Downloaded from https://academic.oup.com/ptep/article-abstract/2017/5/053B05/3835935/Gluon-bound-state-and-asymptotic-freedom-derived by CERN -European Organization for Nuclear Research user on 03 October 2017 g ν k − ν = 0 follows from g ν = gk −⊥ ν . Thus, A μν has the form: Consequently, we define the gluon BS amplitude A μν by Multiplying both sides of A μν (k + , k − ) with g νμ and L νμ from the left, we obtain, using Eq. (51), which has the matrix form: This is solved with respect to A 1 , A 2 : Thus, A 1 and A 2 are extracted from A μν by operating the projection operators P A 1 νμ and P A 2 νμ respectively: Consequently, the homogeneous BS equation for the gluon amplitude is rewritten using the abbreviated form: The homogeneous BS equation for the gluon and ghost amplitudes is obtained: where the matrix element is given as

The coupled BS equation in the massless case
In this section, we study the massless case m = 0 to reproduce Fukuda's results by correcting some errors [7] in our framework, before studying the massive case m = 0 in the subsequent sections. We obtain the numerical solutions of the amputated gluon and ghost BS amplitudes for the coupled homogeneous BS equations in the massless case m = 0. To write down manageable integral equations, we perform the Wick rotation to the relative momenta p and q to obtain the Euclidean momenta p E and q E [28]. For the total momentum P := (E, P), however, we restrict our analysis to the vanishing total momentum P μ = 0. Although this is unrealistic for the bound state, this investigation is a preliminary step for studying the true bound sates dictated by P μ = 0 in the subsequent sections. Even the special result obtained at P μ = 0 is expected to give some helpful information on the mass spectrum P 2 , provided that the Yang-Mills theory has some similarity to the Wick-Cutkosky model [29,30]. In what follows, we adopt the Landau gauge in the Lorenz gauge fixing.
We replace the full propagators and the full vertices by the bare propagators and bare vertices in the BS equation. Then, the homogeneous BS equation reduces to an integral equation of the Fredholm type [31].
First of all, we assume that the gauge coupling constant g does not run, i.e., is standing. After the angular integration in the integration measure, the homogeneous BS equation for the SU (N ) Yang-Mills theory in the Landau gauge at P μ = 0 is cast into a coupled integral equation of the form: with the kernel written in the matrix with the elements: gluon-gluon: ghost-gluon: ghost-ghost: Notice that the integration measure is dq 2 /q 2 , modified from the naive dq 2 , which affects specifying the kernel. It turns out that this choice is efficient for the bound state, as shown in the next section.  the kernel elements are slightly different from those given in Ref. [7]. See Appendix B.2 for more details. First, we take into account only the element a 11 by putting a 12 = 0 to solve the integral equation for the gluon amplitude A 1 alone. The result for the gluon amplitude A 1 is given in the left panel of Fig. 10. Then, the ghost amplitude B is obtained by taking into account a 21 by putting a 22 = 0 from the solution for A 1 . The result for the ghost amplitude B obtained in this way is given in the right panel of Fig. 10.
The BS equation is regarded as the eigenvalue equation, and a BS amplitude as a solution of the BS equation is identified with an eigenfunction associated to each eigenvalue of 1/C. To perform the numerical calculations, we need to divide the integration region [ 2 I , 2 U ] into a sufficient number of small steps. Then the integral equation has the same number of eigenvalues (and the associated eigenfunctions) as the number of partitions. Here we have adopted the infrared and ultraviolet cutoff to be 2 I = 1 and 2 U = 101 respectively, and tried to divide the interval [ 2 I , 2 U ] into 10, 100, and 200 pieces to see whether the result converges or not as the number of partitions increases. The results are shown in Table 1. Figure 10 is a plot of the eigenfunction corresponding to an eigenvalue. Here we present only the graphs obtained from the result of 200 partitions, which gives the convergent result. For this choice of partitions, indeed, we observe that the largest eigenvalue of 1/C and the associated eigenfunction is convergent as confirmed according to Table 1-see the red curves in Fig. 10. The largest eigenvalue of 1/C corresponds to the nodeless eigenfunction, and the node of the eigenfunction monotonically increases as the eigenvalue decreases. Here we define the coupling constant for SU (2) and SU (3) as follows: Second, we take into account all the matrix elements. The numerical solutions are given in Fig. 11. See also Table 2. Figure 11 shows no sizable difference from Fig. 10. This confirms that the above observation is also correct from the quantitative point of view: It is enough to include the gluon-gluon contribution a 11 to study the gluon BS amplitude A 1 .
The solution can also be analyzed by converting the integral equation to the equivalent differential equation, which is simplified by taking into account the dominant element a 11 alone. The resulting equation becomes a fourth-order differential equation with the boundary conditions imposed at the infrared lower bound and the ultraviolet upper bound. In fact, the solution is analyzed in Ref. [22], but the details are omitted here. 16

Gluon-gluon BS equation in the massive case
In what follows, we take into account the gluon-gluon contribution alone in the BS equation. Then the homogeneous BS equation for the amputated BS amplitude D and K 3 and K 4 are the kernels defined by Here, is the ordinary massless gluon propagator defined by It should be remarked that we adopt the massless gluon (with the propagator D   possible to replace the massless gluon mediating the strong interactions by the massive gluon with the same mass m as the constituent gluon, which, however, leads to a very complicated expression for the BS kernel. This issue will be tackled in a subsequent paper (H. Fukamachi, K.-I. Kondo, S. Nishino, and T. Shinohara, work in progress).
In what follows, we restrict our consideration to the BS amplitudeχ ab μν ( p; P) of the transverse type:χ ab μν ( p; P) = δ ab g μν − Indeed, this amplitude is transverse in the sense that Therefore, the left-hand side of Eq. (66) does not depend on the parameter α and reads (0)ac (72) By taking the contraction on the Lorentz indices, the left-hand side reads

Gluon-gluon BS amplitude
We take the approximation in which the angular dependence between the total momentum P and the relative momenta p, q, namely, the terms P · p and P · q, are neglected. Then the left-hand side of Eq. (73) reduces to In the same approximation as that for the left-hand side of Eq. (73), the right-hand side of Eq. (74) reduces to Notice here that g μν K 4 ab;cd μν;ρσ ( p, q; P)χ cd ρσ (q; P) does not depend on P and agrees with the original expression (74).
We perform the Wick rotation to the Euclidean region to perform the momentum integration explicitly. Notice that we apply the Wick rotation to p and q to convert them to the Euclidean momenta p E and q E , while we leave the total momentum P as the Minkowski momentum. After performing the integration over the angular variables in the integration measure d 4 q E for the Euclidean momentum q E , we obtain the integral equation with respect to the magnitude √ q E · q E of the Euclidean momentum: (77b) 7.1. P μ = 0 case First, we consider the limit P μ = 0 of Eq. (77a). This is a first step to examine the existence of the solution of the BS equation, although this limit is obviously unphysical for the bound state problem. This helps us to check whether or not our method works in this problem. The BS equation (77a) reduces to In order to solve the integral equation, we leave only the term in the integrand that is expected to give the most dominant contribution to the integral, and put the momentum cutoff as the upper limit of integration: we can rewrite Eq. (79) into In Eq. (81), we have IR value at p 2 E = 0: while we have UV value at p 2 E = 2 : Therefore, we find that the ratio is given by By differentiating both sides of Eq. (81) with respect to p 2 E , we obtain a differential equation that is equivalent to the original integral equation (81): By solving this equation, we obtain the general solution and the BS amplitude as For the solution (86) of the differential equation (85) to be the solution of the original integral equation, the boundary condition (84) must be satisfied: which leads to the relation When is quite large compared with m, this relation reduces to ln m 2 m 2 + 2 = −γ , γ : Then the mass m obeys the relation Provided that C becomes dependent and that C → 0 as → ∞, γ is approximated as γ ∼ 2 C ln 4 3 and the gluon mass m obeys the scaling law for large , The coupling C must go to zero as → ∞ so that m converges to a finite and non-zero value. This is consistent with the (ultraviolet) asymptotic freedom.

P μ = 0 case
In what follows, we consider the P μ = 0 case.

P μ = 0 case: analytical study
First, we give an analytical treatment. For this purpose, we incorporate the leading terms and neglect the subleading terms among the P 2 -dependent and p 2 -dependent terms in the numerator of the kernel K on the right-hand side of Eq. (77b), where the denominator of the kernel K is regarded as P 2 +4q 2 E . Then the BS equation reads The neglected terms are incorporated in the numerical treatment afterwards.
In the case of P μ = 0, we introduce the reduced variables x, y, r in units of P 2 as For the momentum cutoff of the Euclidean momentum p E , we define the reduced variable λ by In order to consider the two-gluon bound state, we restrict the region of r to Then Eq. (92) is rewritten as Moreover, we introduce the scaled amplitude A by Then the integral equation (96) is rewritten in terms of A(x) as In order to obtain the differential equation equivalent to the integral equation (98), we differentiate Eq. (98) with respect to x. Then the equivalent differential equation is given by The differential equation is solved: Thus, we obtain the BS amplitude: . (101) For a solution of the differential equation to become the solution of the original integral equation, it must satisfy the boundary conditions. The boundary conditions are obtained as follows: • infrared boundary condition by putting x → 0 in Eq. (98): • ultraviolet boundary condition by putting x → λ in Eq. (98): By substituting the solution (101) into Eqs. (102) or (103), we can in principle obtain the relation for C, r, and λ, which we call the scaling relation. However, the resulting expression in this case will be too complicated to extract the physical results we need. Therefore, we move to the numerical investigations.

P μ = 0 case: numerical study
From now on, we proceed to perform the numerical study of the BS equation with the original kernel (77a). We encounter new features once we take into account the p 2 dependence of the kernel (77a). In the numerical solution of the BS equation the eigenvalue of C is not uniquely determined, even if the values of λ and r are fixed, namely, the ultraviolet cutoff and the constituent gluon mass m are given in units of P 2 . In other words, P 2 can have a number of possible values for a given coupling constant C, even if the ultraviolet cutoff and the constituent gluon mass m are fixed. These values converge to a set of eigenvalues if the number of partitions for approximating the integral are increased. Figure 12 shows the real part of the BS amplitude as a function of p 2 /P 2 (momentum in units of P) for the cutoff λ = 10 4 and a small value r = 1. The plots are given for the gluon-gluon BS amplitudes A associated to four different eigenvalues. The resulting BS amplitudes associated to the different eigenvalues have different numbers of nodes. In general, the eigenvalue and the associated BS amplitude can be complex valued. We have confirmed that the numerical solutions for the BS amplitude are real valued when the largest eigenvalue is real valued. In fact, we have checked that the imaginary parts vanish for the numerical solutions. The numerical solution of the nodeless BS  amplitude with the largest eigenvalue plotted in Fig. 12 is consistent with the continuum solution (101). Figure 13 is a plot of the coupling C as a function of the ratio 1/r := P 2 /m 2 for various but fixed values of the cutoff λ, which represents the scaling relation. Figure 14 is a plot of the coupling constant C as a function of the momentum scale λ for various but fixed values of 1/r := P 2 /m 2 , which is required for the scaling relation represented by Fig. 13 to hold. This is the renormalization group flow which gives the same physics, P/m = const., namely, the same bound state mass for a given gluon mass m. This result will be valid for relatively large λ due to the approximations adopted in this paper. As the cutoff λ increases, the coupling C decreases like 1/ ln λ and finally vanishes as λ → ∞: This behavior is consistent with the (ultraviolet) asymptotic freedom expected to hold in the Yang-Mills theory. 24 Figure 15 shows a plot of the β function β(C) := λ dC dλ calculated from the running of the coupling C with respect to λ (Fig. 14) for various values of P/m. We find that all the β functions are negative and coincide for small C or large λ independently of P/m. The resulting β function can be compared with the beta function β(C) = − 11 3 C 2 obtained by one-loop calculation in perturbation theory, which is indicated by the solid curve in Fig. 15. We find that all the points are well fitted to the single curve, although they are located slightly above the solid curve. We observe that the β function obtained from the numerical solution of the BS equation has the same sign and order as the one-loop β function. It should be remarked that the exact agreement of our result with the one-loop prediction was not anticipated, since we have already neglected some terms in obtaining the kernel (77a). Our results indicate that the bound state solution of the BS equation is consistent with the asymptotic freedom in the Yang-Mills theory. This is one of the main results of this paper.
Finally, we comment on how the derivative of C with respect to λ is calculated numerically. The derivative was obtained by taking the finite difference between the smallest eigenvalues of C at λ and 1.1λ for any fixed value of r −1 = P 2 /m 2 . This approximation for the derivative becomes worse for 25/32 Downloaded from https://academic.oup.com/ptep/article-abstract/2017/5/053B05/3835935/Gluon-bound-state-and-asymptotic-freedom-derived by CERN -European Organization for Nuclear Research user on 03 October 2017 larger λ. In view of this, we have used the data for the cutoff λ = 10 5 , which is the largest possible value we can take due to the limitation of the computer memory available to us.

Conclusion and discussion
We summarize the results obtained in this paper. First, we have given a systematic derivation of the BS equation for the gluon and ghost BS amplitude by starting from the CJT effective action for the massive Yang-Mills theory. The solutions of the BS equation represent the simultaneous bound states for gluons and ghosts.
Second, we have obtained the numerical solutions of the derived BS equation in the ladder approximation with the standing gauge coupling, and the improved ladder approximation with the running gauge coupling. As a warm-up problem, we restricted ourselves to the vanishing total momentum P μ = 0. The we have found that (i) the gluon-gluon contribution is dominant in the solution, (ii) the ghost BS amplitude is much smaller than (at most 10% of) the gluon BS amplitude, and (iii) the infrared behavior is influenced by the effect of the running coupling constant.
Third, armed with these preliminary investigations, we proceeded to obtain an approximate solution for the gluon-gluon BS amplitude alone, which corresponds to a two-gluon bound state as a candidate for a color-singlet scalar glueball. We have calculated the beta functions through the running of the coupling constant by using the numerical solutions of the BS equation with an approximated kernel for different choices of the parameters, the ultraviolet cutoff, and the constituent gluon mass. We have shown that the numerical results converge to a single function that has the same sign and order of magnitude as the one-loop beta function, although exact agreement with one-loop perturbation theory is not reproduced due to the approximation we adopted in obtaining the kernel of the BS equation. Thus, we have shown that the solution for the gluon bound state obtained in our approximation for the BS equation is consistent with the (ultraviolet) asymptotic freedom of the Yang-Mills theory.
It is obvious that we have a number of unsolved problems to be tackled: improvement of the kernel for the BS equation by incorporating more terms, examination of the high-energy behavior for performing the normalization of the BS amplitude, renormalization of the BS equation, improvement of the approximation to keep gauge invariance, e.g., through the Ward-Takahashi relations, confirmation of physical unitarity of the massive Yang-Mills theory, and so on. We hope to report some of the results in the near future.