Glueball scattering cross section in lattice SU(2) Yang-Mills theory

We calculate the scattering cross section between $0^{++}$ glueballs of the SU(2) Yang-Mills theory on lattice using the indirect (HAL QCD) method. We employ the cluster-decomposition error reduction technique and use all space-time symmetries to improve the signal. The relation between the interglueball cross section and the scale parameter $\Lambda$ is determined as $\sigma_{\phi \phi}$ = (3.5 - 8.0) $\Lambda^{-2}$ (stat.+sys.). From the observational constraints of galactic collisions, we obtain the lower bound $\Lambda$>60 MeV. We also discuss the naturalness of the Yang-Mills theory as the theory explaining dark matter.


I. INTRODUCTION
The study of dark matter (DM) [1][2][3][4] is one of the most fundamental subjects of physics. Its presence is providing us the most consistent explanation for many astrophysical and cosmological phenomena and problems. There are, however, no candidates of DM in the standard model (SM). In this context, an additional Yang-Mills theory (YMT) which does not or very weakly interact with the SM particles is an attractive candidate, since the lightest particle of the spectrum is a glueball, fulfilling the conditions required for being DM [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]. The glueballs, like the other hadrons, are generated by the nonperturbative physics of the nonabelian gauge theory, so the lattice gauge theory simulation is required to quantify its dynamics .
Here we challenge the quantification of the interglueball cross section. The DM self-interaction (scattering) may affect the structure of the galactic halos and their collisions [62]. The DM scattering cross section is actually constrained by observations [63]. We expect the lattice calculation of the interglueball scattering to yield quantitative relation between the cross section and the unknown scale parameter Λ of the dark YMT, which will be bounded by the observational data. In this work, we use the HAL QCD method [64,65], which is quite successful in the determination of the interhadron potential on lattice, to quantify the glueball scattering.
This paper gives the complete and detailed discussion of the letter [66], in which the DM cross section within the YMT was derived for the first time on lattice. In the next section, we will explain the naturalness of the * ynodoka@yukawa.kyoto-u.ac.jp YMT as a candidate of theory beyond the SM. We then describe in Section III our setup of lattice simulation and the calculation of the interglueball potential obtained by employing the HAL QCD method. In Section IV, we show the result of the glueball cross section in SU (2) lattice YMT and constrain the scale parameter from observation. The final section is devoted to the summary.

II. NATURALNESS OF DARK YANG-MILLS THEORY
In this section, we shall show that the glueball DM is very natural among the composite DM scenarios and dark gauge theories. The assumption of the existence of hidden gauge sectors is indeed a reasonable answer to the problem of the ad hoc gauge group of the SM, since the dark gauge theories may arise naturally in many contexts such as in string or grand unification theories [5,6,14,[67][68][69][70][71][72][73][74][75][76][77][78]. Here we do not discuss their origin, but rather the constraint on the dark gauge sector that could be imposed assuming the naturalness.
The DM model has to be conceived respecting many phenomenological constraints such as the electric neutrality, the nonrelativisticity, the consistency with the bigbang nucleosynthesis, experimental data from direct and indirect detections, etc. The simplest possibility is the massive dark photon. To give it a finite mass, the Higgs potential has to be present with a scalar quadratic term. The latter is known to introduce a severe hierarchical problem due to the quadratic divergence, unless the Higgs potential is just an effective one below an energy scale close to the dark photon mass. The same reasoning applies for non-abelian gauge theories with scalar field(s), whether it (they) induces the Higgs mechanism or not.
The second class of theories is the gauge theories with-out spontaneous breakdown by the Higgs mechanism. In this case, we can first conceive an asymptotically free gauge theory. However, we have to note that the fermions must be massless (or very light compared to the current temperature of the Universe), which will conflict with the nonrelativistic property and the bigbang nucleosynthesis. In this class of gauge theories, the colored particles should, therefore, be confined, and the mass of hadrons is generated dynamically. The chiral fermions are not favored for several reasons as shown below. If the fermions are chiral i.e. massless, the general case will be the spontaneous break down of chiral symmetry. Massless Nambu-Goldstone modes are then generated, which are also not allowed phenomenologically, due to the constraint on the number of relativistic particles at the bigbang nucleosynthesis. We may equally think of a case where the chiral symmetry is not spontaneously broken, requiring the generation of massless composite fermions with the same global symmetry as the original theory, as required by the 't Hooft anomaly matching. This scenario is again forbidden for the same reason as the previous case. If the elemantary fermions become massive due to another spontaneous breakdown of the chiral symmetry through the Yukawa interaction (like the SM), this means that there is at least one additional scalar field which forms the Higgs potential, thus generating again the hierarchical problem. We finally arrive at the conclusion that the most natural scenarios are nonabelian gauge theories with vectorlike fermions. In QCD-like theories, the vectorlike fermions may have arbitrary masses. If the vectorlike fermions are lighter than the confining scale (Λ Nc ) and have gauge charges of the SM sector, the dark baryon number asymmetry may be generated through the sphaleron process. On the contrary, if the masses of vectorlike fermions are heavier than Λ Nc , the glueball of this gauge sector will become the DM. If we wish to explain the DM within the nonabelian gauge theory with vectorlike fermions, we therefore have two choices, the DM composed of baryons made of vectorlike fermions, or lightest glueballs for which the YMT is the relevant theory. Needless to say, we may also conceive a pure YMT without any other fields. In this case, the number of input parameters is minimal, so that the YMT is the most natural theory explaining the DM. Another feature which has to be emphasized is that, in grand unification scenarios, Λ Nc is controlled by the integer number of colors N c which runs the coupling logarithmically over the energy scale, and it may generate a variety of energy scales.

A. Simulation setup
Let us now present the detail of the simulation of the YMT on lattice. In this work, we use one lattice spacing corresponding to β = 2.5, which is sufficient in the HAL QCD method, because the formalism does not depend on the renormalization scale [64,65]. We generated about 2M and 100k configurations with the volumes 16 3 × 24 and 32 3 × 24, respectively, using the pseudo-heat bath method. The first configuration was generated after 5000 (10000) sweeps with cold start, followed then by taking 150 (250) sweeps of separations with the 16 3 × 24 (32 3 × 24) lattice.
The relation between the lattice spacing and Λ, the scale parameter of the SU (2) YMT, was calculated for the general N c [79,80], giving By using the result of the calculation of the string tension a √ σ = 0.1834(26) calculated on 16 4 lattice [45], we have a = 0.107(8)Λ −1 for β = 2.5. The scale of the lattice is now expressed in the unit of the scale parameter Λ. We note that Λ is an unknown parameter, since the property of the DM particle is totally unknown. The aim of our work is precisely to calculate the dependence of the interglueball scattering cross section on Λ, so that the constraint from observational data will give a bound on it.
Let us now define the 0 ++ glueball operator: Here P ij are the plaquette operator in the i − j direction, with a e 1,2,3 the unit vector. The 0 ++ glueball has the same quantum number as the vacuum, so it has an expectation value which corresponds to a divergence in the continuum limit. To extract physical information, we have to subtract it from the glueball field operator φ.
The glueball correlators are then expressed in terms of In order to improve the glueball operator, we also apply the APE smearing [39,50,51]. The smeared link operator U (n) i is constructed by maximizing where Here we choose α = 2.0 and n = 28 to obtain an optimal improvement of the glueball effective mass. We compare in Fig. 1 the effective mass plots with the smeared and unsmeared glueball operators. We see that the smeared operator requires much less imaginary time to form plateau, and the statistical error is much smaller. With the above setup, we found the glueball mass m φ = 0.6857 (28) in lattice unit (value at t = 3).
In the physical unit, we have m φ = 6.40(48)Λ, with the largest uncertainty coming from the relation (1).

B. Nambu-Bethe-Salpeter amplitude
To extract the scattering cross section between two hadrons, we have to calculate the Nambu-Bethe-Salpeter (NBS) amplitude, defined as follows: (5) Here J is the source operator which has the same quantum number as the two-glueball state, and it may be improved using the APE smearing. However, the sink operators should not be smeared, since the nonlocality of the operator will affect the definition of the interglueball distance. For the 0 ++ glueball, all n-body operator (φ n ) (n ∈ N, n = 0) may be chosen. We also note that J also has expectation value, so we have to subtract it. We can easily show that the removal at one side, either at the source or at the sink, is sufficient, since we have For the computational convenience, we choose to remove the expectation value of the source J . In Fig. 2, we show the NBS amplitude with the 1-, 2-, and 3-body sources. We see that for the 1-body source calculation, the NBS amplitude becomes zero in the long range (large r), while it is finite for the cases of 2-and 3-body sources. This is because the NBS amplitude with 1-body source cannot be splitted into two spatially separated correlators, while that with 2-and 3-body sources can. We note that this fact does not occur if the expectation values of the source (and sink) operator are not correctly subtracted, and this proves that the cluster decomposition principle is correctly working. Since the lattice calculation becomes noisier as the mass dimension of the correlator increases, we will mainly discuss the NBS amplitude with the 1body source. The correlator (5) is purely gluonic, and the statistical error is significant in the lattice calculation. To improve the accuracy, we use all space-time symmetries (space-time translation and cubic rotation) to effectively increase the statistics.

C. Finite volume effect
Let us compare the lattice calculations with two different volumes 16 3 × 24 and 32 3 × 24. We plot in Fig. 3 the results of the NBS amplitude (1-body wall source) calculated with approximately the same statistics, which are considered to be proportional to the volume thanks to the use of translational symmetries. We see that the results are in agreement, although the statistical error is large for the case of 32 3 × 24 lattice. It is actually known that the noise enhances when distant gluonic correlations (disconnected insertions) contribute to the correlator, since such contribution increases with the volume, as we chose the wall source. In the next section, we precisely use this fact to reduce the statistical uncertainty, by removing the uncorrelated contribution.

D. The cluster-decomposition error reduction technique
As we saw in Section III C, the contribution to the correlator from gluonic operators, such as the glueball or disconnected quark loops, that are four-dimensionally well separated yields statistical fluctuation which accumulates with the increase of the volume. This is because these operators have expectation values and they are fluctuating even when they are isolated. The idea then came to remove these meaningless fluctuations originating from the distant positions of the operators, and just keep the contribution from the true correlation with closely located interpolating fields. Since the YMT has a mass gap, the correlation is exponentially suppressed in the Euclidean space, so that we can set a four-dimensional cutoff which removes the contribution from the uncorrelated region [81]. This is also an important technical application of the cluster decomposition principle. We plot in Fig. 4 the example of the calculation of the glueball two-point correlator using the above mentioned method, the clusterdecomposition error reduction technique (CDERT). We found that, with the cut ρ = 8 (lattice unit), the systematic error is less than the statistical one. In this work, we apply the CDERT used in Ref. [81] to the calculation of the NBS amplitude. We set cutoffs to the relative four-dimensional distances between the source operator and the sink ones. The NBS amplitude with the 1-body source is then calculated as where C(t, v) is the 3-dimensional projection of the 4dimensional hypersphere with the center (t, v) and with the radius ρ onto the t = 0 (3-dimensional) hyperplane. The cutoffs are simultaneously applied to both glueball operators of the sink, which means that the wall source is changed to a source with a restricted region where two spheres of the cutoffs overlap. We note that this simultaneous application of the cutoff to the above two pairs of operators upsets the possibility to reduce the computational cost using the Fourier transform, which was very efficient in the case of two-point functions [81].
In the case of the SU (2) YMT, the calculation of the NBS amplitude with the CDERT is computationally the most costly step of the whole project. We also apply the same cutoff to the distance between the two interpolating fields of the sink, but this manipulation is just equivalent to not considering the radial plot of the NBS amplitude for the radius beyond the cutoff. In Fig. 5, we plot the result of the application of the CDERT to the NBS amplitude with the cutoff ρ = 8 (lattice unit). We see that the CDERT could successfully reduce the statistical error by more than a factor of two, keeping consistency with the wall source calculation.

E. Problems with the direct method
We now extract the scattering phase shift. It may be calculated by Fourier transforming the NBS amplitude and by inspecting the momentum modulation of the energy of the two-glueball system (the so-called Lücher's method) [82]. This approach is expected to be applicable when there are no states with smaller energy or when the quantum number of the two-body system forbids transition to lighter states. However, in the case of the two-glueball state, it mixes with the single glueball state. Indeed, the NBS amplitude forms a plateau with the effective mass of a single glueball (see Fig. 6). It is, of course, possible to remove the contribution from the one-glueball state by hand or by diagonalizing the source operator, but we have to keep in mind that the glueball spectrum has other states and resonances with energy close to the two-body threshold. The extraction of the two-glueball scattering in the direct method is therefore very challenging.

F. HAL QCD method
An alternative approach to calculate the scattering phase shift is to indirectly calculate it via the potential obtained from the HAL QCD method [65]. This method Effective mass of NBS (unit: t/a BS (p=0,10000Conf) Teper (16 4 ,1998) FIG. 6. Effective mass plot of the glueball NBS amplitude in the momentum space, with 1-body source (J =φ). We see that the energy of the system saturates at the single glueball mass.
has the crucial advantage that we do not need the ground state saturation for extracting the potential [83]. In particular, the glueball correlators are in general very noisy, so the use of this method is almost mandatory if one wants to keep good statistical accuracy. In addition, the potential handled in the HAL QCD method does not depend on the renormalization scale [64,65]. We, however, have to keep in mind that the potential is not an observable, and it may depend on the choice of the operators. We now describe the formalism of the HAL QCD method. The object is to extract the nonlocal potential U ( r, r ′ ) according to the following time-dependent Schrödinger-like equation [83] 1 4m φ Here we have to choose t so that 1/t is less than the inelastic threshold E T = 3m φ −2m φ = m φ , if one wants to study the elastic scattering. In our setup, we have m φ = 0.6857 (28), which implies that we have to take the data for t ≥ 2 (t = 0 corresponds to the contact term). The elasticity is an important feature of the HAL QCD method, since it makes the physics essentially nonrelativistic. The potential should then be local and central, U ( r, r ′ ) ≈ V φφ ( r)δ( r − r ′ ), to a good approximation. We then have We neglect higher order terms of the derivative expansion which contain the nonlocal and angular momentum dependent terms. In Fig. 7, we plot the results of the calculation of the interglueball potential with 1-and 2-body wall sources. We see that the results are consistent at short distance. We also remark that the error bars are very large in the region r ≥ 0.5 Λ −1 . This is partly due to the fact that the potential is calculated by dividing by the NBS amplitude [see Eq. (9)], which is zero in the long range for the case of the 1-body source, and finite but covered by the statistical fluctuation for the 2-body one (see Fig. 2). We fitted the region r ≥ 0.5 Λ −1 with a constant, and obtained negative values for both 1-and 2-body source calculations. These nonzero values are maybe suggesting the attraction at long distance, but it is hard to set conclusive statement due to the statistical error. We do not discuss this topic further in this paper.
Let us now fit the potential calculated with the CDERT. We plot in Fig. 8 the result of our calculation. First, by noting that the glueball operator has a spread of one lattice unit, we remove the data points at r = 0 and r = 1 (in lattice unit) since they are contact terms. We see from our result that the interglueball potential is repulsive in the short range region r ∼ 0.2 Λ −1 . We also see there a kind of oscillation, which should be due to the discretization effect of partial waves [84]. Also, as pointed above, the potential becomes very noisy at long distance.
Let us try to analyze our result in terms of a simple scalar effective field theory. The general renormalizable effective theory of the 0 ++ glueball is given by the following trilinear + φ 4 Lagrangian The above trilinear interaction (term with A) generates an attractive Yukawa potential which has the longest range. The φ 4 interaction may be repulsive depending on the sign of the coupling λ, but this one is rather responsible for a shorter range potential. This naive effective picture is apparently in contradiction with our result depicted in Fig. 8. The trilinear coupling A of the glueballs in SU (3) YMT has actually been extracted in Refs. [38,85] using Lüscher's finite volume method [82,86], and A ∼ O(100) Λ was predicted. To check the consistency with our result, we will have to refine the data points in the long range region, since the one-glueball exchange process with A induces an attractive force of the Yukawa type with the longest tail. As an aside, we note that the glueball mass m 2 φ , the trilinear coupling A, and the φ 4 coupling λ scale as O(N 0 c ), O(N −1 c ), and O(N −2 c ), in the large N c limit, respectively. This results in the scaling of the interglueball cross section σ φφ = O(N −4 c ) in the same limit.
We choose to fit the potential with two fitting forms, i.e. the Yukawa function and the Yukawa+Gaussian After the fit, we find The results are plotted in Fig. 8. The fit with the single Yukawa function has a large χ 2 , whereas the one with the Yukawa+Gaussian function works reasonably well. This is obviously because the latter one has two variables. We, however, also consider the fit with the single Yukawa function to estimate the systematic error.

B. Dark matter cross section and constraint on Λ
We now solve the s-wave Schrödinger equation with the fitted potentials: and extract the scattering phase shift δ(k) at vanishing momentum k → 0, which is precisely the relevant case for the DM scattering. Asymptotically, the solution of Eq. (15) behaves as φ(r) ∝ sin[kr + δ(k)]. The quantity we have to calculate is then We finally obtain the following scattering cross sections for the two fitting forms we used: σ φφ = (7.5 − 8.0)Λ −2 (Yukawa + Gaussian), (18) where the ranges correspond to the statistical error bar. The systematic error is just estimated by taking the difference between the two fits. Finally, the interglueball scattering cross section in the SU (2) YMT is σ φφ = (3.5 − 8.0)Λ −2 (stat. + sys.).
We can then derive the constraint on Λ from observational data. By equating the upper bound on the DM cross section obtained from the galactic collision [63] σ φφ /m φ < 0.47 cm 2 /g, we finally have In this work, we do not discuss the lower bound on Λ which might be set by the inspection of the astrophysics at the scale below kpc [62].
Using the large N c argument, we can qualitatively extend the discussion to all N c . As we saw in Section IV A, the cross section scales as 1/N 4 c , while the mass of the 0 ++ glueball is constant at large N c . The lower limit of the scale parameter is then extended to N c ≥ 3 as MeV.
We note that the contribution of nonplanar diagrams to σ φφ is of O(N −6 c ) (see Fig. 9), so higher order corrections in the 1/N c expansion are not small. To control the systematics down to the percent level, we have to calculate the interglueball cross section up to N c = 10. . We note that all planar diagrams are of the same order. The next-to-leading order correction is given by nonplanar diagrams, such as (c), which are of O(N −4 c ).

V. CONCLUSION
In this paper, we calculated the interglueball scattering cross section using the HAL QCD method, and derived the constraint on the scale parameter of the SU (2) YMT from the observational data of galactic collision. In our work, the use of the HAL QCD method was almost mandatory, since this allowed us to take data with small imaginary time, which was crucial for extracting physical information from the very noisy glueball correlators. We also used the CDERT which was shown to be efficient in the improvement of the signal.
Theoretically, the glueball DM is a natural conception because the theory does not depend on any massive parameters except the scale parameter. It is also part of the nonabelian gauge theory with heavy vectorlike fermions. The other dark gauge theories have more or less hierarchical problems, thus proving the attractiveness of the YMT.
The interglueball potential showed us several important characteristics. At short distance (r ∼ 0.2 Λ), the interglueball potential is repulsive. This feature is difficult to explain in a simple renormalizable scalar field theory in which the glueball potential is attractive at the distance r ∼ 1/m φ due to the one-glueball exchange process. A more sophisticated effective field theory based on the conformal Ward identity has to be investigated in the future. We will also have to confirm the presence of the attractive long range potential due to the one-glueball exchange which is currently covered by statistical noise, although our constant fit is somewhat suggesting its existence.
In our work, we assumed the local potential, but the systematics due to the nonlocality has to be checked in the future, since we defined the glueball operator with plaquettes which have spatial extent. The investigation of the nonlocality can be rephrased as the inspection of the operator dependence, which was already discussed for the mesonic and baryonic systems [87][88][89].
Moreover, we have only calculated the DM cross section in the SU (2) YMT, and the cases of N c ≥ 3 were just qualitative extrapolations using the 1/N c expansion in the present paper. Since the interglueball cross section, being an O(N −4 c ) quantity, receives a correction of O(N −6 c ), we expect to complete the analysis of the 0 ++ glueball of the SU (N c ) YMT as the DM candidate for all N c with the accuracy of O(1%) by accomplishing the calculations until N c = 10. This project is in our view not impossibly challenging.