Density of States FFA analysis of SU(3) lattice gauge theory at a finite density of color sources

We present a Density of States calculation with the Functional Fit Approach (DoS FFA) in SU(3) lattice gauge theory with a finite density of static color sources. The DoS FFA uses a parameterized density of states and determines the parameters of the density by fitting data from restricted Monte Carlo simulations with an analytically known function. We discuss the implementation of DoS FFA and the results for a qualitative picture of the phase diagram in a model which is a further step towards implementing DoS FFA in full QCD. We determine the curvature $\kappa$ in the $\mu$-$T$ phase diagram and find a value close to the results published for full QCD.


Definition of the model and the density of states
We study SU(3) lattice gauge theory with static color charges. The dynamical degrees of freedom are SU(3)-valued gauge links U ν (n), ν = 1, 2, 3, 4, where n = ( n, n 4 ) denotes the sites of a N 3 s × N t lattice with periodic boundary conditions. The action is given by Re Tr U µ (n)U ν (n +μ)U † µ (n +ν)U † ν (n) − η n e µNt P( n) + e −µNt P( n) , (1) where β is the inverse gauge coupling, η the coupling strength of the static color sources and µ is the chemical potential. The static color sources are represented by Polyakov loops P( n) = 1 3 Tr Nt−1 n 4 =0 U 4 ( n, n 4 ). In the action (1) the chemical potential µ gives a different weight to charges P( n) and to anti-charges P( n) , such that the theory has a complex action problem which is equivalent to the one of QCD. For and vacuum expectation values of moments of X[U] can be computed as moments of x in the corresponding integrals over the density ρ(x). The expression (4) makes clear the emergence of the complex action problem in the DoS formulation: The density ρ(x) is integrated with the oscillating function cos( ξ µ x ), and from the definition of the coupling ξ µ in (2) it is obvious that the frequency of the oscillation increases exponentially with µ (and linearly with η), such that ρ(x) has to be computed with sufficient accuracy. The recent developments of DoS techniques [16] - [27] are based on new strategies for calculating ρ(x) with very high precision.

Using DoS FFA for computing the density of states
For a DoS calculation the density ρ(x) has to be parameterized on the interval [0, x max ] in a suitable way. For our parameterization we divide the interval [0, x max ] into N sub-intervals I n ≡ [x n , x n+1 ], n = 0, 1, ... N − 1 with x 0 = 0 and x N = x max . The density is then written in the form ρ(x) = e −l(x) , where l(x) is a continuous function that is piecewise linear on the intervals I n . We normalize the density using the condition ρ(0) = 1, which corresponds to l(0) = 0. Together with the continuity of l(x) this implies that only the slopes k n , n = 0, 1, ... N −1, which determine l(x) in the intervals I n appear as the parameters of ρ(x). A short calculation [24] - [27] gives the explicit form of ρ(x) as function of the k n , Here ∆ j ≡ x j+1 − x j denotes the size of the j-th interval. We stress that the intervals can be chosen with different sizes such that in regions of x where ρ(x) varies quickly a finer discretization can be used. For computing the slopes k n in the DoS FFA we use restricted vacuum expectation values X n (λ), n = 0, ... N −1, which depend on a free parameter λ ∈ R. They are defined as where θ n (x) = 1 for x ∈ I n , θ n (x) = 0 for x ∈ I n . The free parameter λ, which the restricted vacuum expectation values X n (λ) depend on, can be used to probe the system. We stress that the restricted vacuum expectation values are free of the complex action problem and can be evaluated with standard Monte Carlo calculations.
The key observation of the DoS FFA is that with the parameterization (5) the restricted partition sum Z n (λ) and thus the restricted vacuum expectation value X n (λ) = ∂ ln Z n (λ)/∂λ can be computed in closed form. It is convenient to shift and rescale the X n (λ) into a new form Y n (λ) for which one finds the explicit expression [24] - [27]: where in the last step we introduced the function h(s) = 1/(1 − e −s ) − 1/s − 1/2. We find that the shifted and rescaled expectation value Y n (λ) is given by h((λ − k n )∆ n ), i.e., it depends only on one of the parameters of ρ(x), the slope k n in the respective interval I n . Thus we can compute Y n (λ) for several values of λ and determine k n from a fit of the data for Y n (λ) with h((λ − k n )∆ n ). The function h(s) approaches ±1/2 for s → ±∞, is increasing monotonically and obeys h(0) = 0. Thus it is gives rise to simple stable 1-parameter fits and the k n can be determined reliably [24] - [27]. Analyzing the quality of the fit is an important self consistency check and poor quality of the fit indicates that the size ∆ n of the corresponding interval should be chosen smaller [24] - [27]. Once the slopes k n are computed, the density ρ(x) can be determined with (5) and from ρ(x) we can evaluate the observables. Before we come to presenting our results for observables in the next section, we have a look at the density and how it changes when varying the parameters. In Fig. 1 we show ln ρ(x) as a function of x for different values of the inverse coupling β at fixed η = 0.04, µ = 0.15. When comparing the different curves for ln ρ(x) over the full range of x in the lhs. plot, the different couplings seem to give rise to essentially the same density. However, the zoom into the small-x region (rhs. plot) reveals that the curve for the largest β shows a quite different behavior. We will see below that this change corresponds to a phase transition between β = 5.60 and 5.70. We conclude that inspecting the qualitative behavior of the density ρ(x) already reveals physical properties of the system. However, we stress again that only the evaluation of physical observables is the true benchmark for a DoS calculation, since the density still has to be integrated over with the highly oscillating factors, which tests if the evaluation of ρ(x) is sufficiently accurate. We conclude this section with a short comment on the piecewise linear parameterization of the exponent of ρ(x). It is clear that the exact result is obtained only in the limit where one sends the number of intervals N to infinity and their sizes ∆ n to 0. In our studies for this paper, as well as in [27], where we analyzed a related model where we could systematically compare the DoS FFA results to reference data from a dual simulation free of the complex action problem, it was found that the accuracy of the Monte Carlo results has a considerably larger effect on the final results than the size of the discretization intervals we use here. More specifically, our discretization intervals ∆ n were chosen such that an optimal use of the data for Y n (λ) in Eq. (7) is obtained: One can show [27] that the slope of the fit function h((λ − k n )∆ n ) at λ = k n is given by ∆ n /12, and that ∆ n /12 ∼ 0.1 − 0.5 gives rise to an optimal fit of the data for Y n (λ). This choice from [27] was implemented in our study here.

Observables and results
The observables we study are the expectation value of the imaginary part of the Polyakov loop and the corresponding susceptibility. Their definitions and expressions in terms of density integrals are given by Note that in leading order we have ξ µ = 2η sinh(µN t ) ∝ 2η µN t , such that in this order ∂ ln Z/∂ξ µ ∝ 1/2η ∂ ln Z/∂µN t , indicating that Im P is closely related to the particle number density, and χ Im P to the particle number susceptibility, which makes them suitable observables for assessing the phase diagram.
In our numerical simulations we use N = 256 intervals for the discretization of ρ(x) for x ∈ [0, x max ] (for our 8 3 × 4 lattices -for smaller test volumes see the discussion of the volume dependence in the next In order to analyze the dependence of the cost on the volume we implemented a small case study at β = 5.45, µ = 0.25: In addition to V 4 = 8 3 × 4, we also did runs at V 4 = 6 3 × 4 and V 4 = 4 4 , and adjusted the number of intervals N such that the interval size ∆ n and thus the discretization of ρ(x) remained constant. Since x max = √ 3N 3 s /2 is proportional to the 3-volume V 3 ≡ N 3 s , keeping the discretization of ρ(x) constant gives rise to a cost factor proportional V 3 , which is the cost factor that is specific for the DoS FFA. However, as for any other Monte Carlo method, we also need to take into account the longer correlation times and the increasing cost of an individual sweep when the volume is increased: Keeping the number of values for λ fixed at 51, we adjusted the number of Monte Carlo sweeps such that the errors for the density are the same on all volumes. This resulted in a statistics of 350, 1500 and 4800 for V 4 = 4 4 , 6 3 × 4 and V 4 = 8 3 × 4, which roughly scales like (V 4 ) 1.25 . Finally, the cost for one local Monte Carlo sweep is proportional to V 4 , such that we expect that the cost of FFA roughly scales like V 3 × (V 4 ) 2.25 , where only the factor V 3 is specific for DoS FFA, while the other factor is more general and will also depend on the couplings. It is clear that this brief study provides only a very rough assessment of the cost and a dedicated analysis is necessary for conclusive cost estimates. More complicated is the analysis of the µ-dependence of the cost, since this is expected to very strongly depend on the other parameters. Here we can only refer to our study [27], where this question could partly be assessed through a comparison with dual simulation data in a closely related model.
Before we come to analyzing the physical observables it is interesting to have a look at the integrands of (8) and (9) for the different values of µ used here. In Fig. 2 we show ρ(x) sin(ξ µ x) x (lhs. plot of Fig. 2) and ρ(x) cos(ξ µ x) x 2 (rhs.) as a function of x for different values of µ. While the integrand of (8) remains predominantly positive in the range of µ values we consider here, the integrand of the connected part of (9) develops essential negative regions illustrating that the complex action problem (sign problem) is challenging for that observable already at the values of µ we access here. When increasing µ further, both integrands quickly develop a highly oscillating behavior. From (8) it is obvious that Im P ≡ 0 since ξ µ = 2η sinh(µN t ) vanishes at µ = 0, while χ Im P is non-zero also at vanishing µ. Thus we can use χ Im P for a first assessment of the DoS FFA implementation at µ = 0 where we can compare χ Im P to the outcome of a conventional simulation at µ = 0. The lhs.
plot of Fig. 3 shows the DoS FFA results as well as the results from a standard simulation at µ = 0. We find very good agreement of the DoS FFA data with the conventional simulation, which reassures us about the correctness of the implementation of DoS FFA and its accuracy, but stress again that the situation at µ = 0 is more demanding since there the density is integrated with the oscillating factors.
For the runs at chemical potential µ = 0 we first have to determine the parameters for the simulation, i.e., we need to identify the region with transitory behavior. For this purpose we ran a µ = 0 conventional simulation at different values of η to locate a possible transition as a function of β. Note that with increasing β we decrease the lattice spacing a(β), such that increasing β corresponds to increasing the temperature T = 1/(a(β)N t ).
In the rhs. plot of Fig. 3 we show the expectation value of the absolute value of the Polyakov loop |P| = | n P( n)| /N 3 s . For all η we observe a fast increase of |P| for values of β in the range between β = 5.2 and 5.8. For η = 0 this increase corresponds to the first order phase transition (in the thermodynamical limit) that leads from the confined to the deconfined phase. This transition can be understood as spontaneous breaking of the Z 3 center symmetry. A non-zero η breaks this symmetry explicitly and thus one expects that above some critical value of η the phase transition turns into a crossover. We also observe that the position of the transition is shifting towards smaller β (i.e., towards smaller temperature) when increasing η. We illustrate the physical picture in the schematic plot on the lhs. of Fig. 4: The full red curve in the µ = 0 plane indicates the line of first order transitions, which bends towards smaller β when η is increased. Based on the argument with the explicit breaking of center symmetry discussed above, beyond some critical η we expect only crossover type of behavior which we indicate by using a dashed instead of a full curve for that pseudo-critical line.
In the diagram in the lhs. plot of Fig. 4 we also display the µ axis. When considered in the space of all three parameters β, η and µ we have a (pseudo-) critical surface that separates the confined and deconfined phases. This surface runs through the (pseudo-) critical line in the β-η plane. An interesting question is how this surface bends for µ > 0, and in the figure we show in blue the trajectories (straight lines parallel to the β-axis at finite η and different values of µ) along which we compute our observables to explore the bending of the (pseudo-) critical surface.
The results for Im P (lhs. plot) and χ Im P (rhs.) as a function of β at different values of µ are presented in Fig. 5. Both observables show a maximum at the transition between confinement and deconfinement (an exception is Im P which vanishes at µ = 0 as discussed above). We observe that for increasing chemical potential the maxima and thus the transition shift towards smaller β, i.e., towards smaller temperature. In order to determine the critical line in the µ-β plane we fit the data points near the maxima of χ Im P with a cubic polynomial and so determine the peak positions of χ Im P as a function of β for different values of µ, i.e., we determine β c (µ). The results of this determination are used for the plot of the phase diagram in the rhs. plot of Fig. 4.
For the presentation of the results for the (pseudo-) critical line in the rhs. plot of Fig. 4 we converted lattice units to physical units using the scale determined for the Wilson gauge action in [28]. On the vertical axis we use the temperature T in units of the critical temperature at vanishing chemical potential, i.e., we plot T /T c (0). On the horizontal axis we plot the baryon chemical potential µ B = 3µ in units of the critical temperature at that µ, i.e., the combination 3µ/T c (µ). The results of our determination of T c (µ) from the maxima of the susceptibility are shown as asterisks in the rhs. plot of Fig. 4. With increasing µ we observe the bending of the (pseudo-) critical line towards lower temperature values as expected also in full QCD. This bending can be be quantified with the curvature κ defined via the relation (again we use µ B = 3µ) The fit of our data with this quadratic polynomial is shown as the full curve in the rhs. plot of Fig. 4. The small-µ data are described reasonably well and we obtain a value of κ = 0.012(3) for the curvature. We stress that this result is of course a very crude estimate, since we work with only a single lattice spacing and do not attempt a thermodynamic limit. Nevertheless it is interesting to note that our result is in the vicinity of the curvature values published for full QCD in different settings, e.g., κ = 0.0135(2) [29], κ = 0.0149(21) [30] and κ = 0.020(4) [31].

Concluding remarks
In this letter we have presented an exploratory study where the Density of States Functional Fit Approach DoS FFA was implemented for SU(3) lattice gauge theory with static color sources. The purpose is to further develop the DoS FFA method towards its use in a full lattice QCD simulation at finite density.
The key challenge of DoS calculations is the determination of the density ρ(x) with sufficient precision, such that one can reliably determine physical observables by integrating ρ(x) with the oscillating factor, where the frequency increases exponentially with the chemical potential. In our test we demonstrate that for the system of SU(3) lattice gauge theory with static color sources the DoS FFA method can be implemented and the accuracy is sufficient for an evaluation of Im P and χ Im P . Comparing the DoS results to a conventional simulation at µ = 0 shows good agreement and for µ > 0 the observables and the critical line could be determined up to moderately large values of µ. A determination of the curvature κ in the µ-T phase diagram gives a value which is surprisingly close to the results published for full QCD.
We stress again, that the results presented here should not be considered as a final determination of the phase structure of SU(3) lattice gauge theory with static color sources, since infinite volume illustrates the phase boundary between the confined (small β) and the deconfined region. We use a dashed line to indicate that above some critical η one expects only a crossover type of behavior, while at small η the deconfinement transition is of first order (full curve). The blue lines at different non-zero values of µ illustrate the lines in parameter space along which we evaluate observables to probe the curvature of the critical surface. Rhs. plot: Results for the critical line in the µ-T plane at fixed η = 0.04. The data points on the (pseudo-) critical line were determined as the maxima of a cubic fit of the data points for χ Im P . In the µ-T plane we fit the data with a quadratic polynomial to determine the curvature κ as explained in the text (the result of this fit is shown as full curve).  extrapolation and continuum limit were not attempted here. The purpose of the paper is to document the further development of DoS FFA and to explore the steps necessary towards getting the method ready for a full QCD calculation.