QCD mesonic screening masses using Gribov quantization

The screening masses of mesons provide a gauge invariant and definite order parameter of chiral symmetry restoration. Different mesonic correlation lengths for flavor non-singlets, at least up to NLO, are well-defined gauge invariant physical quantities calculated earlier using the perturbative resummation techniques. The NLO perturbative results match the available non-perturbative lattice QCD results at the high-temperature regime. We have studied the spatial correlation lengths of various mesonic observables using the non-perturbative Gribov resummation, both for quenched QCD and (2 + 1) flavor QCD. The study follows the analogies with the NRQCD effective theory, a well-known theory for studying heavy quarkonia at zero temperature.


Introduction
Understanding experimental results at the relativistic heavy-ion colliders in high-energy nuclear and particle physics and the early cosmos' cosmic history depends heavily on quantum chromodynamics (QCD) at finite temperatures. When the temperature rises, the theory switches from a confined phase with hadronic degrees of freedom, where chiral symmetry is spontaneously broken, to a deconfined phase of quarks and gluons, where chiral symmetry is restored.
The temporal and spatial directions are generally disconnected because the heat bath destroys the Lorentz symmetry at finite temperatures. The temporal correlation functions in the Minkowski space can be used to define spectral functions for various operators in Fourier space. The spectral functions then indicate the plasma's essential "real-time" features, such as the particle production rates [1]. In investigations involving heavy ion collisions, these are then immediately observable. Nevertheless, the correlation functions in the spatial direction answer questions like At what length scales thermal fluctuations correlation occurs? External charges are screened on what length scales? These static observables are also physical quantities eminently suited for lattice experiments' measurements.
The usual perturbation theory has a severe infrared problem at the higher orders in coupling [2,3], but the infrared problems are only connected with the static, soft gluons. It has led to the development of the dimensional reduction technique [4,5], which reduces the non-perturbative infrared behavior to a more straightforward effective theory for the soft modes only, enabling the precise estimation of the observables in the first order of expansion in perturbation theory. However, the systematic expansion in the perturbation theory for arbitrary orders remains unknown. This way, the terms involving dynamical quarks can be computed perturbatively in four dimensions. At the same time, the computationally costly non-perturbative methods can be used for a more straightforward three-dimensional theory. Naturally, the QCD coupling is significant at moderate temperatures, and higher-order corrections cannot be disregarded.
The dimensional reduction has been extensively used in QCD and electroweak theory to study thermodynamics [5,6,7,8,9,10,11,12] and different gluonic correlators [13,14,15,16,17] in the bosonic sector. The fermionic modes, on the other hand, only affect the dimensionally reduced theory's parameters, which are generally integrated out from the effective theory. However, several intriguing observables built of quark fields are sensitive to infrared physics and, as a result, require some resummations, elegantly arranged as successive effective theories. The initial attempts to include this class of correlators in theory with reduced dimensions go back almost 30 years [18,19]. These works, however, did not consistently include every term in a fixed order. In ref. [20], it was demonstrated that the dimensionally reduced theory could be expressed in nonrelativistic quarks, and the accurate power counting of various operators was established. The other motivation for studying the observables built of quark fields, i.e., "mesonic and baryonic" observables, is related to their more closely signs with experimental signatures. In recent years, significant work has been done on one class of these observables, namely the various quark number susceptibilities (see ref. [21,22,23,24,25,26,27,28,29,30]).
The other class of observables, more closely related to IR physics and more sensitive toward it, is the correlation lengths of the mesonic operators, which are gauge-invariant quark bilinears constructed from the light quark flavors [31]. The quark bilinears' two-point spatial correlation functions can determine the correlation lengths. Generally, these correlation functions are dominated by the screening masses at large distances, defined as the inverse of the correlation lengths. These screening masses give information about the quark-gluon plasma (QGP) response when a meson is included in the system. A large number of studies have been done in the literature to study the meson correlation function using perturbative methods as well as using non-perturbative lattice computations. On the perturbative side, the mesonic correlation functions using hard thermal loop approximation (HTL) have been studied in refs. [32,33]. Additionally, the meson correlation function at finite momentum in QCD plasma has been studied in ref. [34]. Using the perturbative QCD, the leading order result (2 ) for the meson screening mass is obtained in ref. [35,36]. The next-to-leading (NLO) correction for meson and gluon screening masses has been calculated using an effective theory at zero chemical potential in ref. [37,38]. The connection between thermal screening masses and real-time rates has been explored in ref. [39] in the spectral representation. The NLO correction of meson screening masses at finite chemical potential has also been extended in ref. [40].
On the lattice side, the temporal and spatial hadronic correlation function in QCD plasma has been studied in ref. [41,42,43]. Also, screening masses in purely SU (2) and SU(3) gauge theories using lattice QCD is studied in ref. [44,45,46,47]. Recently, mesonic screening masses have been studied for the first time in a large range covering the temperature from ∼ 1 GeV to ∼ 160 GeV using lattice QCD in ref. [48]. Also, the recent lattice QCD results for the meson screening masses in (2 + 1) flavor QCD in the temperature range from 0.14 GeV ≤ ≤ 2.7 GeV is presented in ref. [49] and for the recent study of mesonic screening mass at finite chemical potential using lattice QCD see ref. [50,51]. In the high-temperature limit, the perturbative result of meson screening masses calculated in ref. [37] is comparable with the lattice results. However, no such analytic calculation in the literature can explain the lattice data for ≤ 2 GeV in the low-temperature regime. In this article, we have tried to overcome this gap by using a non-perturbative scheme offered by Gribov quantization.
Since QCD's infrared (IR) region is strongly coupled, the conventional resummed perturbative approach will not give appropriate results [52]. One of the effective ways to handle the IR region is to consider the Gribov-Zwanziger approach [53,54] (For recent reviews on Gribov-Zwanziger approach, see [55,56] and some recent works on extended Gribov-Zwanziger approach follow [57,58,59,60,61,62,63]) and the references therein. In these extended Gribov-Zwanziger models, first introduced in ref. [57], the effect of local composite operators, which consist of a mass term, has been explored, and it has been found that the inclusion of the mass term in the propagator gives results which match very well with the lattice simulations results in the infrared domain. Also in ref. [64] the authors have calculated the gluon and ghost propagators to study the IR physics at zero temperature by including the mass term in usual Faddev Popov action, which are in excellent agreement with lattice results. Although the infrared domain of QCD is nonperturbative, once the Gribov copies effects are introduced, QCD can be treated perturbatively (or semi-perturbatively) in the infrared. Also, it has been observed in Monte-Carlo simulations that the Yang-Mills coupling constant in the infrared is compatible with a perturbative expansion once the mass effects in the propagator have been included. For more details on this, see the recent review in ref. [65] In recent years the general Gribov approach, without any mass term in the propagator, has been utilized to improve the thermodynamic quantities in QCD by evaluating free energy [66]. Also, in kinetic theory, the transport coefficients have been explored using this scheme in ref. [67,68,69]. In this approach, the gluon propagator modifies its form, and a new (chromo)magnetic scale enters the theory through the mass parameter. Using the modified gluon propagator in Gribov action, quark dispersion relation [70], the dilepton production rate along with quark number susceptability [71] and electromagnetic debye mass [72] have been calculated, which sheds light on some of the new physical insights of the said observables. Recently, this approach has been applied to heavy quark phenomenology to study the heavy quark potential [72,73,74] and heavy quark diffusion coefficient [75].
The following is how the paper is set up. We shall establish necessary notations and review some of the mesonic correlators' well-known characteristics at high temperatures in section 2. Section 3 covers how a dimensionally reduced effective field theory [20] can describe high-temperature QCD in general while incorporating the mesonic correlators along the lines of work [37]. In section 4, we will do the required matching computation to fix the unknown parameter in the effective lagrangian by matching the dispersion relation of QCD and NRQCD 3 , using Gribov formalism. Section 5 evaluates the dynamics of the effective theory with Gribov's gluon propagator inclusion. Results of the screening masses for the quenched QCD case and = 3 are compared with the recent lattice data in section 6. We summarize the work in the last section, 7.

High-temperature static correlators in QCD
The Euclidean Lagrangian of the quarks at finite temperature QCD is given by Here, the covariant derivative ( ) of the fermionic field is given by where are the generators defined in fundamental represenation of group ( ), which are hermitian in nature. The quark field is an -component vector in flavor space. For the sake of simplicity, we will assume M to be a diagonal and degenerate matrix, M = diag( , ⋯ ), and we will take the assumption of = 0. We may use quark fields to define bilinear objects with varying spin and flavor structures. We are interested in the operators having the form, where Γ can take the values {1, 5 , , 5 } for the different channels, namely scalar, pseudoscalar, vector, and pseudo vector , , and respectively. The traceless matrices along with the identity matrix provides the flavor basis, with , = 1, 2, ⋯ 2 − 1. For the operators defined in eq. (3), we will focus on the correlators having the structure as or, in position space, Using the rotational invariance, one can choose the correlation measurement direction as in eq. (5). With this choice, one can take the average over the 1 2 surface, giving a correlator that is easier to handle where is given in eq. (6). The spatially separated two-point correlation function mentioned in eq. (7) defines the screening masses as The above eq. (8) describes the correlation function's exponential falloff at large distances. Let us focus on the correlation function's behavior at high temperatures. Due to asymptotic freedom, the correlators with free fermions depicted in figure (1a) can be evaluated using perturbation theory. By using the dimensional regularisation technique, one would get the result for the diagram (1a) as Here, denotes number of colors, Γ is the dirac matrix coming from the opeartor form , ∕ ≡ and refers to the fermionic Matsubara modes having the form, = 2 ( + 1∕2). After doing the trivial Dirac algebra calculation, which gives some constant terms, the above correlation function in eq. (9) contains the function 3d ( 2 ) which comes out to be From the eq. (10), it is clear that the singularity of the function 3 (2 0 ) appears at the point 2 . Thus, the correlator dominates at large distances for the zeroth Matsubara frequencies defined as ± 0 = ± . Now, this frequency excitation occurs in the 3-dimensional theory. One can also consider this three-dimensional theory as (2 + 1)dimensional and consider the time direction as the one in which the correlation is assessed. In this case, screening mass is equivalent to determining the screening states of a pair of on-shell heavy quarks, each with a mass of " 0 ".

Effective Theory
We want to determine the next-to-leading order (NLO) correction to the mesons screening mass using the nonperturbative effects. To go beyond the leading order, we must consider the diagrams shown in fig. (1b) and (1c). Nevertheless, as we know, in finite temperature QCD, it will not be sufficient to calculate only these two diagrams. In principle, we can have a large number of diagrams that can behave in the same order as the two depicted in fig. (1b) and (1c) do. In other words, to find the bound state (mesons) poles, one needs to do the entire sum of ladder diagrams since we can not obtain bound states by considering a finite number of Feynman diagrams. In our approach, we will use the diagrams above to find the interaction potential, which will be further used in the Schrödinger equation to find the non-relativistic bound states. This approach is equivalent to summing all the ladder diagrams in the non-relativistic limit.
To calculate the first-order correction to the free quarks propagator in a consistent way, it is necessary to do the resummation of the diagrams with an arbitrary number of zero modes of the low-momentum exchanged gluon between quarks. The effective field theory approach provides a convenient way to perform such resummations. Only the diagram containing one-gluon exchanged would contribute to order 2 depicted in fig. (1c). As shown in fig. (1c) for non-zero gluon modes, the gluon exchange diagrams do not need to be computed because this diagram alone cannot alter the pole location or the screening mass. These diagrams would contribute to the overall normalizing factor of the correlator under consideration.
Conversely, in the bosonic section, the dynamics of the gluonic zero modes after the dimensional reduction technique describe the effective three-dimensional gauge theory [5], known as Electrostatic QCD(EQCD). The action takes the form as where the dots refers to higher dimensional operators [76] and = 1, 2, 3 , The temporal component of gauge field 0 acts as a scalar field with mass . However, we are interested in the processes which occur at ( 2 ); then one can integrate out the scalar field from the theory in eq. (11). After integrating the temporal part, we are left with the effective theory known as Magnetostatic QCD (MQCD), which is given by This three-dimensional theory has non-perturbative dynamics; thus, it must be handled completely in a non-perturbative way [2].

Tree-level NRQCD 3
In the fermionic sector, we will consider only the lowest fermionic modes 0 ≡ because the other modes will not contribute much to the correlator at large distances as seen in eq. (10). The Euclidean quark lagrangian for the field ( ) having Matsubara mode reads as where = 1, 2 and 0 is the zero gluonic mode interacting with ( ). We considered the Euclidean Dirac matrices a different form than the standard one. In this representation, we have As the quarks are very heavy, the quark fields can be considered static fields, and the Dirac spinors can be rewritten in terms of and where , are the two-component spinor objects that led to the lagrangian where is antisymmetric two rank tensor and 12 = +1. The quarks in the considered correlators are almost "onshell," and for the free theory, the on-shell point is given by 2 0 + 2 = 0, i.e., 3 = ± 0 . At the tree level, quarks with a fixed Matsubara frequency interact with zero gluonic modes only; thus, it is expected that quarks' off-shellness is related to the gluonic momentum scale as | 3 ± 0 | ≲ . Now, after expanding around a state consisting of a free pair of quark-antiquark, one of the components is heavy, and another is light compared to the dynamical scale, viz. , 2 , one can solve the equation of motion for the heavy component and the light mode. After doing an expansion in powers of 1∕ 0 , the effective action for the fermionic mode can be written as

Power Counting Arguments
We must consider the quantum corrections subjected to the parameters in eq. (17) to ascertain the radiative corrections. One must take into account all additional operators that symmetries might permit. Setting up a power counting for the likely operators is crucial.
In the NRQCD 3 side, quarks have momentum scale as |p ⟂ | ≲ and their off-shellness Δ 3 = 3 + 0 is of the same order. By requiring the action mentioned in eq. (17) to be of the order of unity, we have On-shell relativistic gluons have the same energy and momentum order viz., 3 ∼ p ⟂ . On the other hand, as the considered quarks are nonrelativistic, their kinetic energy is directly proportional to momentum squared in (2+1) dimensional theory. For a nearly on-shell quark with transverse momentum |p ⟂ | ≲ in NRQCD 3 , the off-shellness in longitudinal momentum becomes Δ 3 ∼ 2 ⟂ ∕ 0 ∼ 2 and this gives 3 ∼ 2 that acts on quarks.
Similar power counting arguments demonstrate that ∼ 3∕2 , which means that this term will be of higher order compared to the derivative term present in the transversal covariant derivative. As 2 ⟂ is already of ( 2 ), one can left out the trasnverse gluons. Thus the only parameter which needs to be matched beyond tree-level is the zero point energy 0 ≡ . Thus to find out the ( 2 ) corrections to meson correlation lengths, the following lagrangian is sufficient to use where  3 = 3 − 3 . Note that, here 0 and 3 play an important dynamical scale role. On the other hand, transverse gluons 1 , 2 can be ignored as long as one is interested in the energy shift of ( 2 ). This means that gluons will not transfer transverse momentum to quarks in this order. To be consistent at ( 2 ), we should replace, energy i.e., 0 of tree-level effective Lagrangian by a matching coefficient = 0 + ( 2 ), which we will determine in the next section.

Matching Conditions from QCD to NRQCD 3 , with Gribov
We will determine the one-loop correction to the variable by matching the Green function, calculated on the QCD side and NRQCD 3 side, using Gribov formalism. Here, the matching will be done by finding the finite temperature Euclidean dispersion relation so that we do not need to worry about the overall normalization factors arising from the fields. This computation produces gauge-invariant results in both sectors, and we will do the matching order by order in 1∕ 0 on the NRQCD 3 side so that dimensional regularisation copes up with the power counting inputs as done in ref. [77].
In the conventional quantization of QCD, the gauge condition is not ideal, as proposed by Faddeev and Popov. Thus, there remains a residual gauge discrepancy in the infrared region, which is pointed out by Gribov and known as the Gribov copy problem [53]. Gribov proposed that the functional integral should be constrained to the first Gribov region, which collects gauge fields without zero modes for the Faddeev-Popov operator. The expression for the Gluon propagator with Gribov quantization in the general covariant gauge is written as where = 0, 1 corresponds to Landau and Feynman gauges. On the QCD side, using the Gribov propagator mentioned in eq. (20), the quark propagator's inverse, defined as −1 ( ), is given by, where the gluonic four-momentum in Euclidean space-time reads , with 0 = 2 , and represents the number of colors. This calculation uses a dimensional regularization technique with MS renormalization scheme. In this technique, the sum-integral is given as with = 4 − 2 as the space-time dimensions. Also, (⋯) , (⋯) refers to fermionic and bosonic Matsubara modes respectively, and = ( 2 −1)∕2 . As we are calculating the ( 2 ) corrections to Σ( ), one can utilize the free theory constraints namely 2 = 0 and ∕ ( ) = 0 to handle the terms which are in proportion to 2 and ∕ . Thus the longitudinal gauge-dependent part then vanishes, leaving only the transverse terms of Gribov's gluon propagator to contribute to the calculation. The outcome is irrespective of the gauge component as expected.
Because of the rotational invariance of the remaining terms in eq. (21), we can do the matching of the Green's function at some particular momentum and set ⟂ = 0. After the multiplication of 0 from the left to eq. (21), the expression becomes block-diagonal. Thus we can focus on the particular term, say, [ 0 −1 ( )] 11 , which is given by In eq. (23), we have the contribution coming from the transverse gluons, i.e., 1 and 2 component only. In contrast, the gauge field components 0 and 3 vanish altogether at the free theory pole position 0 = − 3 as they occur with opposite signs. Thus the only integral which we need to evaluate to find the ( 2 ) correction is Now, the first integral 1 becomes Here, ± represents the Bose-Einstein (B.E) distribution function with ± = √ 2 ± 2 . Additionally, the second integral 2 takes the form as where is given by wherẽis the Fermi-Dirac (F.D) distribution function with energy = − 3 (1 + cos ). Thus, the Euclidean dispersion relation on the QCD side using the Gribov propagator becomes On the NRQCD 3 side, one can extract the Feynman rules from the lagrangian in eq. (19). The one-loop contribution to the quark self-energy vanishes on NRQCD 3 side because the two gauge field components, namely 0 and 3 , are the same and come with opposite signs as is the case in QCD side. Thus the inverse propagator in NRQCD 3 becomes So, the pole location is simply 3 = on NRQCD 3 side. Now, after doing the matching, we will get The integrals 1 and 2 need to be evaluated numerically to find the matching parameter . This matching accounts only for the hard gluons in the loop integral of figure (1b).

Dynamics of effective theory
After considering the hard gluons for figure (1b), we still need to take care of the soft gluons involved in figure (1b) and as well as fig. (1c), i.e., we need to solve the dynamics of the effective theory at ( 2 ). The dynamics will be studied by computing the correlators, as done in section 2, by using the effective theory mentioned in eq. (12), (19). The modified form of the correlation function of eq. (7) is One may assume that the exchange of gluons between quarks and antiquarks occurs instantly in the nonrelativistic domain and may find out the static potential for the quark-antiquark pair by integrating the gauge fields. After obtaining the potential, we may use the usual Schrödinger equation to solve for the screening states and arrive at a solution.
Thus we need to establish Green's function for the quantity of the form ∫ d 2 , where Γ is a constant term which consists of 2 × 2 matrix which does not affect the overall calculation if we drop out that term. To find out the static potential, one can use the trick of introducing a point-splitting in the correlator, which means putting the quark pairs at a particular distance apart, which is finite, and then determining the Schrödinger equation obeyed by the correlator and set → 0 afterward. Thus the updated correlator is For the Gribov's gluon propagator ( ), we have where the Hamiltonian can be read from the eq. (19). For the one-loop diagrams shown in fig. (1b) and (1c), the equation of motion is given by where the kernel  is a dimensionless quantity that can be expanded to its first two arguments to find out the one-loop static potential, Now the eqs. (34), (35), (36) can be clubbed together to get the final equation of motion upto one-loop order as For the Gribov propagator, the expectation values of temporal gauge fields and spatial gauge fields are the same in the computation of the static potential. After solving the above integral in eq. (36), we get the final form of one-loop static potential as where 0 is the modified Bessel function whose asymptotic behavior is given by 0 ( ) = − ln 2 − + ( ). Now, this potential will determine the correlation length −1 = , through where Ψ 0 represents the ground state wave function. We will numerically solve this eq. (39) to find the screening mass.
For the easiness of the problem, we can do the re-scaling as In the polar coordinates, the Schrödinger equation (39) modifies as . To find out the numerical solution of eq. (41), we need to find out the wave function 0 ( ′ ) around the origin. Thus one finds that [37], Ψ 0 (0) is some finite number to make the wave function bounded. To find out the 0 , one can integrate the eq. (41) for a large value of ′ and by requiring the square integrability condition. Contrary to the perturbative case [37], we find temperature dependent 0 that saturates to a constant value at high temperature.

Results and Discussion
This section will outline our findings for the screening masses and contrast them with the most recent lattice data obtained in ref. [49]. The screening mass can be obtained from the eqs. (30) and (40) as In the calculation of screening mass, we have used the lattice-fitted running coupling obtained in ref. [66], which is given by where = 1.43 for the infrared (IR) case and = 2.97 for ultraviolet (UV) case respectively. To obtain the one loop coupling in eq. (44), the authors fitted the lattice data of running coupling extracted from the IR and UV behavior of heavy-quark free energy from ref. [78]. In our calculation, since we are looking at the non-perturbative region, we use the infrared coupling case. The integrals 1 and 2 consist of the Gribov parameter , which is often determined utilizing the one-loop or two-loop gap equation (for details, see [79]). Figure 2 shows the temperature variation of the scaled Gribov mass parameter ∕ has been obtained after doing the matching with the lattice thermodynamics data, as obtained in ref. [69]. For the values of ∕ shown in figure 2 the integrals 1 and 2 have been evaluated numerically. The integral 2 also contains the imaginary part, which eventually comes from the pole condition used to evaluate integral shown in eq. (27). However, since screening mass depends on the real value of parameter defined in eq. (30), by using eq. (43), we plot the scaled screening mass ∕ with temperature in figure 3 for quenched QCD ( = 0) case and for = 3 case. We found a good agreement with the lattice data reported in ref. [49] and a good improvement over the perturbative results obtained in ref. [37] in the lowtemperature domain. The main outcome of fig. 3 shows that screening mass significantly decreases from the free theory results in the low-temperature region. At the same time, in the high-temperature realm, it approaches the screening mass result obtained in ref. [37].

Summary
In this work, we used the non-perturbative resummation using the Gribov quantization approach to study the mesonic screening masses. We began by examining the fundamental characteristics of the static meson correlators at high temperatures and saw that the zeroth fermionic Matsubara mode dominates the large distance correlator. Then we discussed the framework of the effective theory, namely NRQCD 3 , in which we need to determine the parameters of the effective  lagrangian to find out the non-perturbative NLO correction to mesonic screening mass. This matching is done by finding the Euclidean dispersion relation on the NRQCD 3 side and QCD side, using the Gribov approach. After this, we calculated the static quark-antiquark potential using Gribov's gluon propagator in the context of effective theory to understand the dynamics of the theory. This static potential determines the coefficient of exponential falloff by numerically solving the usual Schrödinger equation. We plotted the screening mass in the temperature ranges from 300 MeV to 2700 MeV for = 0 and = 3, respectively. We compare our results with the recently obtained lattice data [49] and found a good agreement. Since our calculation does not respect the different channels like a scalar, vector, etc., an immediate extension can be done by adopting some strategy to include the different channels and different mesons in this framework. The other natural extension can be done by having the finite chemical potential in theory for a more general case study of mesonic screening masses which can give insights into the nature of the QCD phase diagram, which we will report soon.