A comparative lattice analysis of $SU(2)$ dark gluebal

We study the mass and scattering cross section of $SU(2)$ glueballs as dark matter candidates using lattice simulations. We employ both naive and improved $SU(2)$ gauge actions in $3+1$ dimensions with several $\beta$ values, and adopt both the traditional Monte Carlo method and the flow-based model based on machine learning techniques to generate lattice configurations. The mass of dark scalar glueball with $J^{PC}=0^{++}$ and the NBS wave function are calculated. Using a coupling constant of $\beta=2.2$ as an illustration, we compare the dark glueball mass calculated from the configurations generated from the two methods. While consistent results can be achieved, the two methods demonstrate distinct advantages. Using the Runge-Kutta method, we extract the glueball interaction potential and two-body scattering cross section. From the observational constraints, we obtain the lower bound of the mass of scalar glueball as candidates of dark matter.


I. INTRODUCTION
Dark matter (DM) that gives rise to about a quarter of the universe's mass remains one of the most mysterious objects in particle physics.Primarily recognized through its gravitational influence on galaxies and cosmic structures, it might be also undetectable through electromagnetic means.The weakly interacting massive particles miracle (WIMP), the suggestive coincidence of the weak coupling DM with a proper density that would have thermally been generated, motivated experimental tests through their direct scattering, decays to visible cosmic rays, and DM productions at collider experiments, but no positive result has been reported so far (for an incomplete list on this topics, please see Refs.[1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]).
Aside from WIMP, several other interesting DM models exist.An interesting scenario suggests dark matter is likely to be glueballs, particles composed of strongly interacting dark gluons bound together without fermions, stemming from a confining dark SU (N ) gauge theory extrapolated from Quantum Chromodynamics (QCD).These glueballs in the dark sector are hypothesized to interact with standard model particle mainly through gravitational forces and may contribute significantly to the dark matter content [16][17][18][19][20][21].
Glueballs within QCD are bound states of gluons, characterized by the absence of quark constituents.These entities emerge from non-Abelian gauge theories, notably SU (N ) lattice gauge theories, where non-perturbative effects support the formation of such composite particles [22][23][24].This feature surprisingly bridges the gap between the domain of strong interactions and cosmological phenomena, suggesting a novel view where the dynamics of the strong force support the mysterious nature of dark matter.The similarities underscore the appealing hypothesis that glueballs may indeed constitute the dark matter component in these scenarios [20].Theoretical studies and simulations within this framework have been important in substantiating the glueball dark matter hypothesis, offering a robust statement that integrates the microcosmic interactions of particle physics with the large-scale structures and dynamics observed in the universe.
In nonperturbative studies of QCD glueball, lattice QCD calculations play a crucial role [25][26][27][28][29][30][31][32].This nonperturbative method offers detailed insights into key properties of glueballs, including the spectrum and decay properties.Consequently, it reveals their potential significance across a spectrum of physical phenomena.For the lattice investigation of self-interacting dark matter, the ratio of two-body scattering cross section to the mass of dark matter, denoted as σ/m is of significant interest.It is not just a crucial quantity for experimental physicists seeking to constrain the characteristics of dark matter [33][34][35] but also holds key importance in understanding cosmic phenomena, as highlighted in previous research [36].Recent advances in lattice calculations [37,38] have been also crucial in determining the scattering cross section of dark glueballs within lattice SU (2) Yang-Mills theory.
It should be noticed that traditional Monte Carlo (MC) methods for lattice simulations face issues such as critical slowing down [39] and topological freezing [40].Due to the chain-like nature of the Markov process, configurations generated by traditional methods are not strictly independent unless a very long Markov chain is achieved.Many improvements are developed in the past decades, but another potentially promising method that was recently inspired from machine learning concepts and techniques is the flow-based model [41][42][43][44] (for a review, see Refs. [45,46]).A flow-based model is a type of neural network that can learn a complicated probability distribution from data, and generate new samples from it, without using Markov chain Monte Carlo methods.Moreover, flow-based model can also learn the potential features and correlations of the data, and provide insights into the structure and dynamics of the system.Flow-based models impose an invertible transformation constructed by neural networks transforming a simple distribution to a complex one to perform the maximum likelihood estimation.
In the remnant of this work, we will use the traditional Monte Carlo method and machine learning approach, in a comparative manner, to generate lattice configurations for the SU (2) gauge theory with zero fermion flavors (N f = 0) and several coupling constants β.Subsequently, we compare the results for the dark glueball obtained from these distinct simulations.Furthermore, we calculate the interaction potential for dark glueball and determine the interaction coupling constants in effective Lagrangian.Finally, we obtain the relation between the dark glueball mass and scattering cross section.By confronting theoretical results with experimental data, we extract the lower bound on the mass of scalar glueball as candidates of dark matter.Some details are collected in the appendix.

II. METHODOLOGY
In this section, we will exhibit the theoretical basis of our lattice simulations.This comprises the application of the SU (N ) gauge theory action within both the Monte Carlo method and machine learning approach.We elaborate on the specific machine learning methodology, especially the flow-based model.Additionally, this section also includes the strategy to investigate the dark glueball, displaying how to determine its mass, coupling constants of effective Lagrangian, and scattering cross section.

A. SU (N ) lattice gauge theory
The SU (N ) gauge theory was established for describing interactions of non-Abelian gauge boson, with SU (3) gauge theory being crucial for modeling strong interactions.Explicitly one can discretize the spacetime to simulate QCD in continuum with gauge links U µ (n) in path integral approach [49].Standard simulations of SU (N ) gauge theory make use of plaquettes U µν (n), defined as: with n indicating the lattice coordinates, and µ, ν representing gauge link directions.The discretized action is constructed as: where β is related to the coupling constant: β = 2N/g 2 .N signifies the color degree of freedom within the SU (N ) gauge field, which is set to 2 in this study.The action described previously is referred to as the naive action.
To reduce lattice spacing discretization effects, the Lüscher-Weisz (improved) action is frequently employed [50], in which both plaquettes U µν (n) and rectangles R µν (n) are taken into account: The corresponding action is constructed as: where the coefficients for plaquettes and rectangles are c 0 = 5/3 and c 1 = −1/12, respectively.
Our study aims to apply self-interacting gauge fields to explore dark matter, where we opt for simplicity by setting N = 2.We plan to employ both naive and improved actions for Monte Carlo simulations, and the naive action in combination with machine learning simulations, as detailed in Section III.

B. Flow-based model
In the flow-based model [41][42][43][44], the process begins with a basic prior distribution r(z) that is easy to sample.Utilizing this, a collection of samples z i can be produced and subjected to a reversible one-to-one transformation f (z).Consequently, the transformed samples x i = f (z i ) will inherently conform to a new distribution q(x): The objective is to align the new distribution q(x) more closely with the target distribution p(x), which takes the form p(x) = e −βS(x) /Z for lattice gauge theory.To assess this alignment, one can make use of the Kullback-Leibler (KL) divergence [47].The KL divergence between two distributions q(x) and p(x), denoted as D KL (q|p), is defined as A low KL divergence suggests a close similarity between the two distributions, and D KL (q|p) = 0 only when q(x) = p(x).In this context, one can construct the transformation as a sequence of neural networks and employ the KL divergence between the resultant distribution q(x) and the target distribution p(x) as a loss function for network training.Throughout the training phase, the KL divergence can be approximated using the sample mean value: Within the flow-based model, the design of the invertible transformation f (z) demands careful consideration due to the Jacobian determinant of f (z) influencing the resulting distribution q(x).To optimize efficiency, the transformation is structured to yield an upper or lower triangular matrix for the Jacobian.More explicitly, one can firstly choose certain components of z to transform, such as z 1 , z 2 ...z k (where upper indices denote different components of the variable), and keep the remaining components z k+1 , z k+2 ...z n to be frozen, i.e. f (z j ) = z j , k + 1 ≤ j ≤ n.The overall transformation, which is composed of many single transformations with different frozen freedoms, can transform the samples completely, such that every component of z is transformed at least once.Next, one can define a set of invertible one-to-one mappings F (z; α i ) for 1 ≤ l ≤ k with the α i s being parameters.For example, if z l ∈ R, one can set the one-to-one maps as F (z l ; α i ) = α 1 z l +α 2 with α 1 > 0. Subsequently, one can construct neural networks that take the fixed components as inputs and output the parameters of the mappings, i.e., the networks map the fixed components to the parameters: Here F N N denotes the neural network transformation, which is typically non-invertible or difficult to invert, while the transformation applied to the samples remains invertible.The neural network parameters, denoted as ω λ , need to be optimized.Incorporating fixed components as inputs to the neural networks enhances the overall transformation's expressiveness and fosters correlations among the various components.Under the above transformation, the Jacobian is structured as follows: This results in an upper triangular matrix whose determinant can be efficiently computed.By following these steps, one can illustrate the implementation of the flow-based model in 4-dimensional SU (2) lattice gauge theory.For further insights, refer to [52].The prior distribution is established as a uniform distribution with respect to the Haar measure of the SU (2) group due to its numerical sampling convenience.
To start, the fixed components of the gauge fields should be chosen.For a 4-dimensional lattice with coordinates n i , i = 0, 1, 2, 3, n i ∈ Z, a translational phase is defined based on a translation vector t, where t i ∈ Z: This method divides the entire lattice into four segments based on their respective phase values.During a single transformation, one direction of the link variables within a segment is designated as the actively transformed component, while the others are kept fixed.By combining a minimum of 16 individual transformations, a comprehensive transformation is accomplished.In our numerical experiments, we set the translation vector to (0, 1, 2, 1) and designate U 0 (n) as the active transformed link direction.
Secondly, one-to-one mappings are designed.In gauge theory, one needs to preserve the gauge invariance of the theory.Therefore, as link variables are not gauge covariant, they are not directly transformed.Instead, the plaquettes containing the actively transformed links are the ones that undergo transformation.Describing the oneto-one mapping on plaquettes or in SU (2) group space directly in numerical terms is challenging.However, by considering the theory's requirement for gauge covariance, we can confine the mapping to be f (XP X −1 ) = Xf (P )X −1 , where P and X belong to SU (2).It can be proved that, given the diagonalization of P , P = SU S −1 the requirement is equivalent to Namely only the eigenvalues need to be transformed.The two eigenvalues are not independent of each other due to the constraint det P = 1, which leads to θ 1 + θ 2 = 0.One can design diagonalization procedure to make θ 1 ≥ 0, θ 2 ≤ 0, then The one-to-one mapping now can be simplified as The requirement now reduce to that f (θ) should be a one-to-one map on interval [0, π].In our numerical simulations, we opt for the rational-quadratic neural spline flow method with 5 nodes as the one-to-one mapping on the interval [0, π] due to its high flexibility.This choice is based on the work [53].
In constructing the neural network, maintaining gauge covariance necessitates that inputs consist solely of plaquettes containing frozen link variables.To ensure the translation invariance and cyclic boundary conditions inherent in lattice theory, the optimal architecture involves using convolutional neural networks (CNNs) with cyclic padding.In our simulations, a 4-dimensional CNN with a kernel size of 3 and 48 hidden channels is utilized, employing LeakyReLU as the activation function and the Adam optimizer for training.A mask, shaped to match the lattice plaquettes, is applied to select the frozen plaquettes before they are inputted.
It is necessary to point out that although the flowbased model has shown great success, it encountered the mode collapse problem when β was increased [54].To address this problem, a quenched method was proposed in a study in Ref. [55].The quenched method involves first training a model to generate samples that adhere to the distribution e −β1S(x) /Z 1 , where β 1 is small.Subsequently, this distribution with a small β 1 is utilized as the prior distribution to train a second model that generates samples following e −β2S(x) /Z 2 with a larger β 2 .This process is iterated until the value of β reaches the desired level.In this work we will start with β = 0 and incrementally increase β by 0.1 or 0.2 every 10 or 20 rounds of training until β reaches the target value.Through this way the total training time is reduced compared to the original method, and an advantage is that only one model needs to be saved.The overall training time increases approximately linearly with β and the lattice volume (the number of training parameters increases approximately linearly with the lattice volume), while the time required for generating configurations is only dependent on the lattice volume.At this point, it is important to note that the flow-based model is inefficient in generating configurations with a very large volume.For a comparative analysis in this work, we set β up to 2.2 and the lattice volume up to 10 4 .

Dark Glueball Spectrum
In lattice simulations, the operator for dark glueball with J P C = 0 ++ is given: To account for the vacuum state sharing the glueball's quantum numbers J P C = 0 ++ , the vacuum expectation value is subtracted, reformulating the operator as Ô′ (t, ⃗ x) = Ô(t, ⃗ x) − ⟨ Ô(t, ⃗ x)⟩.We can use two-point correlation functions to determine the mass of dark glueball, When the time slice t is large enough, the two-point correlation function is reduced to The effective energy E 0 represents the energy of the lowest dark glueball state (ground state) and when the momentum is zero, it is reduced to the mass m g .In contrast to the naive action, the operator formulation for the dark glueball in the context of the improved action incorporates an additional rectangular term, as defined in Eq.( 3).Thus the operator becomes: For subsequent analyses, we employ this refined operator representation for the dark glueball within the framework of configurations derived from the improved action.

The two body scattering cross section
In the following, we will determine the two-body dark glueball scattering cross section in various ways.
The cross section of dark glueballs scattering can be inferred from the interaction glueball potential in a nonrelativistic limit.In an approach one can employ the Bonn approximation to determine the cross section from the interaction glueball potential.This potential is modeled in this work, for instance, as a spherically symmetric Yukawa and Gaussian potential: V Gaussian (r) = Z (1) e − (mg r) 2 8 where m g represents the mass of the dark glueball, and Y, Z (1) , Z (2) are parameters.For both Yukawa and Gaussian potentials, one can consider the differential cross section in nonreltivistic limit with k → 0 which is derived as follows: Here ⃗ K = ⃗ p 3 − ⃗ p 1 represents momentum transfer, and | ⃗ K| = 2k sin θ 2 with k and θ being the magnitude and scattering angle of the relative momentum ⃗ k between two interacting particles respectively.V (r) denotes either the Yukawa or the Gaussian potential.
Alternatively, using an effective theory for scalar particles, and considering the ϕ 3 and ϕ 4 interaction, one can construct the following Lagrangian for dark glueballs as: Corresponding to the contributions in the potential, the ϕ 4 interaction term gives a delta function, while it is difficult to be obtained numerically.Therefore, only ϕ 3 term is considered for the scattering.Perturbative theory yields the differential cross section for two-body elastic scattering ϕ(⃗ p 1 )ϕ(⃗ p 2 ) → ϕ(⃗ p 3 )ϕ(⃗ p 4 ).The Feynman diagram in Fig. 1 provides the iM matrix element: Setting | ⃗ K| = 0, we obtain the differential cross section: The ϕ 3 interaction generates Yukawa potential in the coordinate space in non-relativistic limit.If V Yukawa (r) is inserted to Eq.( 20), the differential cross section is expressed as: Comparing Eq.( 23) and Eq.( 24), one can find that the value of λ 3 is determined by the coefficient Y : The s-wave scattering cross section can also be directly extracted from the wave function.This method is grounded in the radial Schrödinger equation, in which the wave function asymptotically approaches the asymptotic form as r → ∞: By solving the radial Schrödinger equation, we can extract the wave function Ψ(r) and determine the phase shift δ l (k).This enables the calculation of both the differential and total cross sections for the s-wave scattering in the k → 0 limit:

Nambu-Bethe-Salpeter wave function
To calculate the coefficient Y in the Yukawa potential or Z (1) and Z (2) in the Gaussian potential, we make use of the Nambu-Bethe-Salpeter (NBS) wave function derived from lattice calculations.The NBS wave function was initially developed for hadrons [57,58], which is also suitable for analyzing glueballs.It is defined as: where Ô′ is a glueball operator, and |GG⟩ represents a hadron state comprising two glueballs.
The NBS wave function related to the Schrödinger equation with the interaction dark glueball potential.The spatial Schrödinger equation with a non-local potential is expressed as: where U (⃗ r ′ , ⃗ r) approximates a local, central potential The relation between V (r) and Ψ g (⃗ r) is: Replacing V (r) with Eq.( 18) or Eq.( 19) allows the determination of the coefficient Y or Z (1) , Z (2) from Ψ g (⃗ r).
In lattice simulations, correlation functions to extract the NBS wave function are constructed as: with Ŝ as the glueball operator source.Then, the wave function is given as In principle a two-body state could serve as the source Ŝ. However for the dark glueball with J P C = 0 ++ , |GG⟩ and |G⟩ share the same quantum numbers.Additionally, one cannot split one-body source into spatially separated correlators, unlike two-and three-body sources.Hence, we use Ŝ = Ô′ .

III. NUMERICAL SIMULATIONS AND RESULTS
In this section, we give our numerical results obtained from lattice simulations.This includes an in-depth overview of the configuration generation process, the results for dark glueball mass, as well as the determination of the coupling coefficient for the effective Lagrangian and the dark glueball scattering cross section.

A. Lattice setup
When generating lattice configurations, we have varied the coupling constant β in the set {2.0, 2.2, 2.4, 2.6} for naive actions in SU (2) simulations.In contrast, for the improved actions, we specifically employed β = 3.5 for comparative purposes.
Table .I delineates the setup of configurations employed in our study.During the generation process, we have used both heatbath and overelaxation updates for the SU (2) gauge field.To enhance signal clarity, we implemented APE (Array Processor with Emulator) smearing for the gauge links as described in [26], In our simulations, we select a smearing parameter of α = 0.1, and a total of 20 APE smearing iterations are uniformly applied across all configurations.Employing the established relation between the string tension σ and the scale parameter Λ as documented in [59], we determine the lattice spacing in terms of Λ −1 for each simulation setting [60].Detailed information is given in Appendix.A, and the results are listed in Table II.In contrast, for the machine learning approach, we restrict our simulations to the naive action at β = 2.2 with a lattice volume of 10 4 .While machine learning yields configurations that closely resemble SU (2) gauge configurations, they are not precise.Therefore, we use the heatbath method to update the gauge links from machine learning.Fig. 2 compares the plaquette values from MC and ML with heatbath updates.This comparison reveals that MC configurations align with ML ones after 25 times heatbath updates.Eventually, we performed 40 times.Furthermore, we have compared the static potential aV (n x ), derived from large Wilson loops as shown  in Eq.(A1), which is illustrated in Fig. 3.In this figure, the static potential displays a positive value of B for the B/r term, which is attributed to APE smearing.Further details are discussed in Appendix A. This comparison supports our conclusion that machine learning with heatbath updates can generate reasonable results as the traditional Monte Carlo method.
In this investigation, we also undertake a comparative analysis of the computational time consumption between the Monte Carlo method and machine learning approach.It should be warned that this comparison is qualitative since for machine learning approach due to the lack of GPU nodes we have used the CPU re- sources across 6 DCU nodes.When implementing the Monte Carlo calculation for 10 4 lattice configurations, the processing time is approximately 0.17 seconds per configuration.This implies that the initial warm-up process requires around 120 seconds, and the total time for generation, including both warm-up and production, is calculated as (Conf×0.17+120)seconds.In contrast, the machine learning method involves distinct training and production phases, consuming 21 hours for training and Conf×0.12seconds for production.Fig. 4 illustrates a cost comparison between the Monte Carlo and machine learning methods.Based on this analysis, it can be inferred that the Monte Carlo method is more efficient for small-scale computations, whereas machine learning demonstrates greater suitability for extensivescale simulations.Therefore, configurations derived from machine learning exhibit reduced correlation and demand less computational time in this example.From this perspective, machine learning may emerge as an efficient method for generation.On the other side, it is necessary to point out that the implementation of machine learning is still limited.A significant issue is the substantial memory requirement during the training process, which may exceed the capacity of limited-core systems.Additionally, the persistent problem of mode-collapse in machine learning necessitates the continued use of the Monte Carlo method for configuration filtering.Given these considerations, our primary simulations employ the Monte Carlo approach, while machine learning is only utilized for a comparative analysis at this stage.

B. Dark glueball mass
The mass of the dark glueball in the framework of SU (2) gauge theory is ascertainable through the analysis of the two-point correlation function C 2 (t), as delineated in Eq. (15).Utilizing the dark glueball operator provided in Eq.( 17), we select µ and ν along the spatial directions to characterize the dark glueball.Considering the ground state contribution, C 2 (t) exponentially decays with respect to time slice t.In light of the periodic boundary conditions applied during the configuration generation, C 2 (t) is modeled as: To clearly demonstrate the value of the dark glueball mass, we presented the effective mass plots in Fig. 5.
Here, we calculated the effective mass m eff , using the following formula: Despite considerable uncertainties in each case, the trend suggests that as t increases, m eff tends to stabilize roughly towards a plateau.To enhance the stability and clarity of our fitting process, we performed a constant fit for m eff .The results are illustrated in Fig. 5 as red bands, whose lengths represent the ranges of used data in fitting.The results are also listed in Table III, which indicates that as β increases, the dark glueball mass exhibits a gradual increase.Moreover, we extrapolate the dark glueball mass to the continuum limit by using the equation: and determine the results m g (a → 0) = 6.35(51)Λ.Since only the naive action is employed in the extrapolation, we retain the O(a) term, which should be neglected with improved action.In parallel, our simulations using  III.The red bands in these plots represent the fitting results, with their lengths corresponding to the fitting ranges.
the improved action suggest an effective glueball mass of m g = 5.72(31)Λ.Comparing this result with those obtained using the naive action, we find that the improved action configurations at such a large lattice spacing a = 0.1900(10)Λ −1 has smaller discrete effect than naive action.We have put together all these results in Tab.III, which also shows very similar values for the naive action at β = 2.2, achieved through both machine learning and Monte Carlo methods.This comparison confirms that traditional Monte Carlo methods and machine learning techniques are consistent with each other.These results from different methods underline the usefulness of machine learning in particle physics, especially for complicated tasks like calculating the dark glueball mass.

C. Scattering cross section
In this subsection, we address the determination of interaction dark glueball potentials using the NBS wave function [57,58], as outlined in Eq. (31).Tackling the numerical second partial derivative in this context is challenging.In the case that the potential is spherically symmetric, we can approach it by solving a radial partial differential equation: Here, R(r) is the radial part of Ψ g (⃗ r).The solution's angular component is expressed as spherical harmonics with quantum numbers l, m, leading to Ψ g (⃗ r) = R(r)Y l,m (θ, ϕ).Averaging Ψ g (⃗ r) over a spherical shell simplifies Y l,m (θ, ϕ) to 1, effectively reducing the NBS wave function Ψ g (⃗ r) to R(r) times a constant.
To proceed, we calculate the two-body correlation function as stated in Eq.( 32) and average it over spherical shells, approximating it as R(r).However, it is derived under the condition t → ∞ which cannot be achieved in lattice simulations.Therefore, we examine R(r) at fixed time separation, denoted as t fix , which are displayed in Fig. 6.We notice that R(r) exhibits modest variation with increasing t, prompting us to utilize the averaged results from several t fix s.Our findings for R(r) from different lattice setups are presented in Fig. 7.
Considering Yukawa and Gaussian forms of potential energy, we set boundary conditions and adjust coefficients Y or Z (1) , Z (2) to find solutions for R(r), noted as R(r, Y ) or R(r, Z (1) , Z (2) ) using the fourth-order Runge-Kutta method.We then measure the differences between the R(r, Y ) or R(r, Z (1) , Z (2) ) curves and the data of R(r), selecting the curve that most closely matches our data to determine the coefficients.Fig. 8 compares the solutions with Yukawa and Gaussian potentials to our lattice data.This shows that solving the differential equation Eq. (39) with the Runge-Kutta method aligns well with our simulations, indicating reliable results.We thus obtained the values for Y in the Yukawa potential, and Z (1) , Z (2) in the Gaussian potential.In addition, we also investigated the χ 2 /d.o.f values of both Yukawa and Gaussian fits, finding them to be 2.915 and 0.764, respectively, over the fitting range of r ∈ [0.4,0.9]Λ −1 .This indicates that the Yukawa form only partially captures the data.Consequently, we decided to increase the uncertainty of the fit parameters by a factor of χ 2 /d.o.f.For clarity, we reformulate these equations as follows: V Gaussian (r) Our current analysis shows that the Gaussian potential better matches the original data, as we see in Fig. 8.We use a combination of bands from several lattice sets for each potential, helping us to find the final coefficients in these potentials.Through the relation between the Yukawa potential's coefficient Y and the coupling coefficient λ 3 in the scalar glueball Lagrangian (see Eq.( 25)), we find λ 3 = 2m g √ Y = 53.3(5.8)Λ.
Because of these findings, we decide to use both potentials in our cross section results.We treat the difference between these two types of potential as a way to estimate the systematic uncertainties.
Having established the values of Y , Z (1) , and Z (2) , we are now equipped to formulate the Schrödinger equation for the radial wave function Ψ(r), which is related to rΨ g (r), in the presence of a non-zero relative momentum k.This equation is expressed as: When r → ∞, V (r) → 0. It suggests that Ψ(r) approaches an asymptotic form Ψ(r) ∝ sin(kr + δ(k)) at large distances.Utilizing the Runge-Kutta method, we can accurately solve for Ψ(r).Applying the parameters derived from the Yukawa and Gaussian potentials, we calculate the s-wave total cross sections separately.The results with only statistic uncertainties are as follows: Our final result for the cross section is given as: which contains both statistic and systematic uncertainties.
Given our finding that both σ and m g vary with Λ, we are able to use m g instead of Λ when representing σ/m.This allows us to match our σ/m result with experimental data and assign an appropriate value to m. Fig. 9 presents our findings alongside data from experimental studies [61][62][63][64].It is evident from our analysis that the ratio σ/m g tends to decrease as m g increases.Additionally, according to current constraint, it is estimated that the value of m g is likely around 0.3 GeV.

IV. SUMMARY AND OUTLOOK
In this work, we have conducted an analysis of 0 ++ glueballs within the SU (2) gauge theory using both FIG. 9.This shows the ratio σ/mg against the dark glueball mass mg.It incorporates the constraints from several research groups that have reported on the cross section of selfinteracting dark matter.The band in the graph represents the range encompassed by our calculations.For context and comparison, the graph also includes results from established literature: σ/m = 0.7cm 2 /g as in Ref. [61], σ/m = 0.1cm 2 /g from Ref. [62], σ/m = 0.47cm 2 /g from Ref. [63], and σ/m = 2.0cm 2 /g from Ref. [64].
Monte Carlo simulations and machine learning techniques.We have used these two approaches to generate various lattice configurations, and accordingly extracted the dark glueball mass using these configurations based on naive and improved actions.Using a coupling constant of β = 2.2 as an illustration, we have compared the dark glueball mass calculated from the configurations generated from the two methods.While consistent results can be achieved, the two methods demonstrate distinct advantages.Given these considerations, our primary simulations have employed the Monte Carlo approach, while machine learning is only utilized for a comparative analysis.
Then the glueball interaction potential, crucial for extracting the interaction coupling constant in effective quantum field theory and determining the glueball scattering cross section, is extracted.We then established a connection between the scattering cross section and the dark glueball mass, a vital aspect for understanding dark glueball behavior in various physical scenarios.Our findings show that the ratio σ/m g decreases with increasing m g , aligning with values determined in experimental studies.This correlation, along with estimated dark glueball mass m g ≈ 0.3 GeV from experimental data, has significant implications for dark matter research.
The application of Lattice simulations to dark glueball not only showcases the potential of computational advancements in theoretical physics but also paves the way for more refined and efficient research methodologies in the future.Our results can contribute to the broader understanding of dark glueball dynamics and their interactions in the universe [65], particularly in relation to dark matter.This study encourages further exploration into complex particle systems and their interactions, which could be pivotal in understanding the mysteries of dark matter and its fundamental forces.to use a test coefficient in V (r) and solve the differential equation via the Runge-Kutta method, starting from R(r ini ) and R ′ (r ini ).The initial condition is derived from the NBS wave function of lattice data, replacing differentiation with difference.
Applying the fourth-order Runge-Kutta method, we rewrite the above equation as: R ′′ (r) = f (r, R, R ′ ) = 2 r R ′ (r) − m g R(r)V (r).(B2) Then we can obtain the solution R(r) by iterating the following equations with step length h: where r 1 , R 1 , and R ′ 1 are the iteratively updated values at r + h.Using test coefficients in V (r), we approximate the corresponding solution, denoted as R test (r).We then assess the distance between this test solution and our lattice data.Then by varying the coefficients within a certain range and finding the minimal distance, we eventually identify the corresponding coefficients for V (r).

FIG. 1 .
FIG. 1.The Feynman diagram of t channel for elastic scattering of scalar field particles.

FIG. 2 .
FIG. 2. The comparison of plaquette based on configurationsfrom Monte Carlo method and machine learning refined with the heatbath method.

2 FIG. 3 .
FIG.3.The comparison of static potential from MC and ML with 40 times heatbath updates.

FIG. 4 .
FIG. 4. Comparison of computer costs from Monte Carlo and machine learning methods.

2 FIG. 5 .
FIG. 5. Effective mass plots are displayed for all cases.The corresponding fitting results are shown in TableIII.The red bands in these plots represent the fitting results, with their lengths corresponding to the fitting ranges.

radial
FIG. 6. Radial NBS wave function R(r) at different time separations.

2 FIG. 7 . 5 FIG. 8 .
FIG. 7.This figure shows the radial part of the NBS wave function derived from various lattice setups.We normalize all data at r = 0.3727Λ, and show the chosen t range in the middle of each panel.As noted in the main text, the results labeled and 'improved' are based on configurations with naive and improved gauge action, and represents results obtained with machine learning methods.

aV 1 FIG. 12 .
FIG.12.Static potentials from configurations with or without APE smearing at β = 2.4.The number in the label after "smear " represents times of smearing.

TABLE I .
Setup for SU (2) Lattice simulations: The initial configurations are random SU (2) matrices, followed by warmup process with heatbath updates, and then production process employing combinations of one heatbath update and multiple overrelaxation updates.

TABLE III .
The results of effective masses of 0 ++ glueball are shown in lattice unit and in Λ.