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

TL;DR

This paper applies the Density of States Functional Fit Approach to SU(3) lattice gauge theory at finite density, providing insights into the phase diagram and estimating the curvature parameter near known QCD results.

Contribution

It introduces and implements the DoS FFA method for SU(3) lattice gauge theory with static color sources, advancing towards full QCD applications.

Findings

Determined the curvature $$ in the phase diagram close to full QCD results.

Demonstrated the feasibility of DoS FFA in a finite density SU(3) model.

Provided a qualitative phase diagram for the model.

Abstract

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.

Figures

Equations

$$S[\mathrm{U}]=-\frac{\beta}{3}\sum\limits_{n}\sum\limits_{\mu<\nu}\mathrm{Re}\Bigl{[}\mathrm{Tr}\,\mathrm{U}_{\mu}(n)\mathrm{U}_{\nu}(n+\hat{\mu})\mathrm{U}_{\mu}^{\,\dagger}(n+\hat{\nu})\mathrm{U}_{\nu}^{\,\dagger}(n)\Bigr{]}-\eta\!\sum\limits_{\vec{n}}\Bigl{[}\,e^{\mu N_{t}}\,\mathrm{P}(\vec{n})+\,e^{-\mu N_{t}}\,\mathrm{P}(\vec{n})^{\star}\Bigr{]},$$

$$\begin{split}S_{\rho}[\mathrm{U}]&\;=\;-\frac{\beta}{3}\sum\limits_{n}\sum\limits_{\mu<\nu}\mathrm{Re}\Bigl{[}\mathrm{Tr}\,\mathrm{U}_{\mu}(n)\mathrm{U}_{\nu}(n+\hat{\mu})\mathrm{U}_{\mu}^{\,\dagger}(n+\hat{\nu})\mathrm{U}_{\nu}^{\,\dagger}(n)\Bigr{]}-2\eta\cosh(\mu N_{t})\sum\limits_{\vec{n}}\mathrm{Re}[\mathrm{P}(\vec{n})]\;,\\ X[\mathrm{U}]&\;=\;\sum\limits_{\vec{n}}\mathrm{Im}[\mathrm{P}(\vec{n})]\quad\;\;\mbox{and}\;\;\quad\xi_{\mu}\;=\;2\eta\sinh(\mu N_{t})\;.\end{split}$$

$$\rho (x)=\int \mathcal{D}[\mathrm{U}]e^{-S_{\rho }[\mathrm{U}]}\delta (X[\mathrm{U}]-x),$$

$$Z=\int \mathcal{D}[\mathrm{U}]e^{-S[\mathrm{U}]}={\int }_{-x_{max}}^{x_{max}}dx\rho (x)e^{i{\xi }_{\mu }x}=2{\int }_0^{x_{max}}dx\rho (x)\cos ({\xi }_{\mu }x),$$

$$\rho (x)=A_n{e}^{-k_{n}x}\text{mbox}forx\in I_n\text{mbox}where{A}_n=e^{-\sum _{j=0}^{n-1}{\Delta }_j(k_j-k_n)}.$$

$$\begin{split}\langle\langle X\rangle\rangle_{n}(\lambda)=\frac{1}{Z_{n}(\lambda)}\int\!\mathcal{D}[\mathrm{U}]\,e^{-S_{\rho}[\mathrm{U}]+\lambda\,X[\mathrm{U}]}\;X[\mathrm{U}]\;\theta_{n}\bigl{(}X[\mathrm{U}]\bigr{)}\,,\;Z_{n}(\lambda)=\!\int\!\mathcal{D}[\mathrm{U}]\,e^{-S_{\rho}[\mathrm{U}]+\lambda\,X[\mathrm{U}]}\,\theta_{n}\bigl{(}X[\mathrm{U}]\bigr{)},\end{split}$$

$$Y_{n}(\lambda)\equiv\frac{\langle\langle X\rangle\rangle_{n}(\lambda)-\!\!\sum^{n-1}_{j=0}\!\!\Delta_{j}}{\Delta_{n}}-\frac{1}{2}=\frac{1}{1\!-\!e^{-(\lambda-k_{n})\Delta_{n}}}-\frac{1}{(\lambda\!-\!k_{n})\Delta_{n}}-\frac{1}{2}\ =h\Big{(}\!(\lambda-k_{n})\Delta_{n}\!\Big{)},$$

$$⟨\mathrm{Im}\mathrm{P}⟩$$

$${\chi }_{\mathrm{Im}\mathrm{P}}$$

$$\frac{T_c(\mu )}{T_c(0)}=1-\kappa {\left(\frac{3\mu }{T_c(\mu )}\right)}^2+\mathcal{O}\left({\left(\frac{3\mu }{T_c(\mu )}\right)}^4\right).$$

Full text

1 Introductory remarks

The success of numerical calculations in lattice field theory relies on the availability of probabilistic polynomial algorithms for computing observables in a Monte Carlo (MC) simulation. The key point is the interpretation of the Boltzmann factor $e^{-S}$ as a probability. However, in some situations the action $S$ acquires an imaginary part that spoils the probabilistic interpretation necessary for a MC simulation, a problem usually referred to as ”complex action problem” or ”sign problem”.

An important class of systems with a sign problem are lattice field theories with non-zero chemical potential. In many cases the complex action problem is the main obstacle for an ab-initio determination of the full phase diagram at finite density. Different methods such as complex Langevin, Lefshetz thimbles, Taylor expansion, fugacity expansion, reweighting, and worldline formulations were applied to finite density lattice field theory (see, e.g., the reviews [1] – [9] at the annual lattice conferences).

Another important general approach are Density of States (DoS) techniques [1], [10] – [27]. Here we develop further the Density of States Functional Fit Approach (DoS FFA) and apply it to SU(3) lattice gauge theory with static color sources (SU(3) LGT-SCS). The DoS FFA was already presented in depth in [24] – [27] and we refer to these papers for a detailed discussion of the method. Here we review only the parts specific for the SU(3) LGT-SCS and the respective observables, which are related to the particle number and its susceptibility used to determine a qualitative picture of the phase diagram.

We stress at this point that the results presented here are not meant as a detailed systematic study of the phase diagram of SU(3) LGT-SCS, which would imply a controlled thermodynamical limit followed by extrapolating to vanishing lattice spacing. This paper serves to document the developments and tests of the DoS FFA in a model which is a further step towards a Density of States calculation in full lattice QCD at non-zero chemical potential.

2 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 $\mathrm{U}_{\nu}(n),\nu=1,2,3,4$, where $n=(\vec{n},n_{4})$ denotes the sites of a $N_{s}^{3}\times N_{t}$ lattice with periodic boundary conditions. The action is given by

[TABLE]

where $\beta$ is the inverse gauge coupling, $\eta$ the coupling strength of the static color sources and $\mu$ is the chemical potential. The static color sources are represented by Polyakov loops $\mathrm{P}(\vec{n})=\frac{1}{3}\mathrm{Tr}\,\prod^{N_{t}-1}_{n_{4}=0}\,\mathrm{U}_{4}(\vec{n},n_{4})$. In the action (1) the chemical potential $\mu$ gives a different weight to charges $\mathrm{P}(\vec{n})$ and to anti-charges $\mathrm{P}(\vec{n})^{\star}$, such that the theory has a complex action problem which is equivalent to the one of QCD.

For defining the density of states we decompose the action $S[\mathrm{U}]$ into real and imaginary parts and write it in the form $S[\mathrm{U}]=S_{\rho}[\mathrm{U}]-i\xi_{\mu}X[\mathrm{U}]$, where it is easy to see that

[TABLE]

The functional $X[\mathrm{U}]$ in the imaginary part is bounded, i.e., $X[\mathrm{U}]\in[-x_{max},x_{max}]$, with $x_{max}=\frac{\sqrt{3}}{2}N_{s}^{3}$. For $x\in[-x_{max},x_{max}]$ we introduce the weighted density of states $\rho(x)$ as

[TABLE]

where $\mathcal{D}[\mathrm{U}]$ is the product of Haar measures for all link variables. Exploiting the transformation $\mathrm{U}_{\nu}(n)\rightarrow\mathrm{U}_{\nu}^{\,\star}(n)$ one finds that $\rho(x)$ is an even function. Thus the partition sum $Z$ in terms of the density reads

[TABLE]

and vacuum expectation values of moments of $X[\mathrm{U}]$ can be computed as moments of $x$ in the corresponding integrals over the density $\rho(x)$. The expression (4) makes clear the emergence of the complex action problem in the DoS formulation: The density $\rho(x)$ is integrated with the oscillating function $\cos(\,\xi_{\mu}\,x\,)$, and from the definition of the coupling $\xi_{\mu}$ in (2) it is obvious that the frequency of the oscillation increases exponentially with $\mu$ (and linearly with $\eta$), such that $\rho(x)$ has to be computed with sufficient accuracy. The recent developments of DoS techniques [16] – [27] are based on new strategies for calculating $\rho(x)$ with very high precision.

3 Using DoS FFA for computing the density of states

For a DoS calculation the density $\rho(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}\equiv[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 $\rho(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 $\rho(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 $\rho(x)$. A short calculation [24] – [27] gives the explicit form of $\rho(x)$ as function of the $k_{n}$,

[TABLE]

Here $\Delta_{j}\equiv 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 $\rho(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 $\langle\langle X\rangle\rangle_{n}(\lambda)$, $n=0,...\,N\!-\!1$, which depend on a free parameter $\lambda\in\mathds{R}$. They are defined as

[TABLE]

where $\theta_{n}(x)=1\text{ for }x\in I_{n}$, $\theta_{n}(x)=0\text{ for }x\not\in I_{n}$. The free parameter $\lambda$, which the restricted vacuum expectation values $\langle\langle X\rangle\rangle_{n}(\lambda)$ 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}(\lambda)$ and thus the restricted vacuum expectation value $\langle\langle X\rangle\rangle_{n}(\lambda)=\partial\ln Z_{n}(\lambda)/\partial\lambda$ can be computed in closed form. It is convenient to shift and rescale the $\langle\langle X\rangle\rangle_{n}(\lambda)$ into a new form $Y_{n}(\lambda)$ for which one finds the explicit expression [24] – [27]:

[TABLE]

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}(\lambda)$ is given by $h((\lambda-k_{n})\Delta_{n})$, i.e., it depends only on one of the parameters of $\rho(x)$, the slope $k_{n}$ in the respective interval $I_{n}$. Thus we can compute $Y_{n}(\lambda)$ for several values of $\lambda$ and determine $k_{n}$ from a fit of the data for $Y_{n}(\lambda)$ with $h((\lambda-k_{n})\Delta_{n})$. The function $h(s)$ approaches $\pm 1/2$ for $s\rightarrow\pm\infty$, 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 $\Delta_{n}$ of the corresponding interval should be chosen smaller [24] – [27]. Once the slopes $k_{n}$ are computed, the density $\rho(x)$ can be determined with (5) and from $\rho(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\rho(x)$ as a function of $x$ for different values of the inverse coupling $\beta$ at fixed $\eta=0.04$, $\mu=0.15$. When comparing the different curves for $\ln\rho(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 $\beta$ shows a quite different behavior. We will see below that this change corresponds to a phase transition between $\beta=5.60$ and $5.70$. We conclude that inspecting the qualitative behavior of the density $\rho(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 $\rho(x)$ is sufficiently accurate.

We conclude this section with a short comment on the piecewise linear parameterization of the exponent of $\rho(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 $\Delta_{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 $\Delta_{n}$ were chosen such that an optimal use of the data for $Y_{n}(\lambda)$ in Eq. (7) is obtained: One can show [27] that the slope of the fit function $h((\lambda-k_{n})\Delta_{n})$ at $\lambda=k_{n}$ is given by $\Delta_{n}/12$, and that $\Delta_{n}/12\sim 0.1-0.5$ gives rise to an optimal fit of the data for $Y_{n}(\lambda)$. This choice from [27] was implemented in our study here.

4 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

[TABLE]

Note that in leading order we have $\xi_{\mu}=2\eta\;\sinh(\mu N_{t})\propto 2\eta\;\mu N_{t}$, such that in this order $\partial\ln Z/\partial\xi_{\mu}\propto 1/2\eta\;\partial\ln Z/\partial\mu N_{t}$, indicating that $\langle\mathrm{Im}\,\mathrm{P}\rangle$ is closely related to the particle number density, and ${{{\chi}}}_{\mathrm{Im}\,\mathrm{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 $\rho(x)$ for $x\in[0,x_{max}]$ (for our $8^{3}\times 4$ lattices – for smaller test volumes see the discussion of the volume dependence in the next paragraph). For each interval we used 51 values of $\lambda$ and a statistics of 4800 configurations generated with a local restricted Metropolis algorithm. The statistical errors for the raw data were computed with a jackknife blocking analysis. For some parameter values in the vicinity of the crossover we increased the number of intervals to $N=384$ and used statistics up to 30000 configurations.

In order to analyze the dependence of the cost on the volume we implemented a small case study at $\beta=5.45,\mu=0.25$: In addition to $V_{4}=8^{3}\times 4$, we also did runs at $V_{4}=6^{3}\times 4$ and $V_{4}=4^{4}$, and adjusted the number of intervals $N$ such that the interval size $\Delta_{n}$ and thus the discretization of $\rho(x)$ remained constant. Since $x_{max}=\sqrt{3}N_{s}^{3}/2$ is proportional to the 3-volume $V_{3}\equiv N_{s}^{3}$, keeping the discretization of $\rho(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 $\lambda$ 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}\times 4$ and $V_{4}=8^{3}\times 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}\times(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 $\mu$-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 $\mu$ used here. In Fig. 2 we show $\rho(x)\,\sin(\xi_{\mu}\,x)\;x$ (lhs. plot of Fig. 2) and $\rho(x)\,\cos(\xi_{\mu}\,x)\;x^{2}$ (rhs.) as a function of $x$ for different values of $\mu$. While the integrand of (8) remains predominantly positive in the range of $\mu$ 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 $\mu$ we access here. When increasing $\mu$ further, both integrands quickly develop a highly oscillating behavior.

From (8) it is obvious that $\langle\mathrm{Im}\,\mathrm{P}\rangle\equiv 0$ since $\xi_{\mu}=2\eta\,\sinh(\mu N_{t})$ vanishes at $\mu=0$, while ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ is non-zero also at vanishing $\mu$. Thus we can use ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ for a first assessment of the DoS FFA implementation at $\mu=0$ where we can compare ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ to the outcome of a conventional simulation at $\mu=0$. The lhs. plot of Fig. 3 shows the DoS FFA results as well as the results from a standard simulation at $\mu=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 $\mu\neq 0$ is more demanding since there the density is integrated with the oscillating factors.

For the runs at chemical potential $\mu\neq 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 $\mu=0$ conventional simulation at different values of $\eta$ to locate a possible transition as a function of $\beta$. Note that with increasing $\beta$ we decrease the lattice spacing $a(\beta)$, such that increasing $\beta$ corresponds to increasing the temperature $T=1/(a(\beta)N_{t})$.

In the rhs. plot of Fig. 3 we show the expectation value of the absolute value of the Polyakov loop $\langle|\mathrm{P}|\rangle=\langle|\!\sum_{\vec{n}}\mathrm{P}(\vec{n})|\rangle/N_{s}^{3}$. For all $\eta$ we observe a fast increase of $\langle|\mathrm{P}|\rangle$ for values of $\beta$ in the range between $\beta=5.2$ and 5.8. For $\eta=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 $\mathds{Z}_{3}$ center symmetry. A non-zero $\eta$ breaks this symmetry explicitly and thus one expects that above some critical value of $\eta$ the phase transition turns into a crossover. We also observe that the position of the transition is shifting towards smaller $\beta$ (i.e., towards smaller temperature) when increasing $\eta$. We illustrate the physical picture in the schematic plot on the lhs. of Fig. 5: The full red curve in the $\mu=0$ plane indicates the line of first order transitions, which bends towards smaller $\beta$ when $\eta$ is increased. Based on the argument with the explicit breaking of center symmetry discussed above, beyond some critical $\eta$ 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. 5 we also display the $\mu$ axis. When considered in the space of all three parameters $\beta,\eta$ and $\mu$ we have a (pseudo-) critical surface that separates the confined and deconfined phases. This surface runs through the (pseudo-) critical line in the $\beta$-$\eta$ plane. An interesting question is how this surface bends for $\mu>0$, and in the figure we show in blue the trajectories (straight lines parallel to the $\beta$-axis at finite $\eta$ and different values of $\mu$) along which we compute our observables to explore the bending of the (pseudo-) critical surface.

The results for $\langle\mathrm{Im}\,\mathrm{P}\rangle$ (lhs. plot) and ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ (rhs.) as a function of $\beta$ at different values of $\mu$ are presented in Fig. 5. Both observables show a maximum at the transition between confinement and deconfinement (an exception is $\langle\mathrm{Im}\,\mathrm{P}\rangle$ which vanishes at $\mu=0$ as discussed above). We observe that for increasing chemical potential the maxima and thus the transition shift towards smaller $\beta$, i.e., towards smaller temperature. In order to determine the critical line in the $\mu$-$\beta$ plane we fit the data points near the maxima of ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ with a cubic polynomial and so determine the peak positions of ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$ as a function of $\beta$ for different values of $\mu$, i.e., we determine $\beta_{c}(\mu)$. The results of this determination are used for the plot of the phase diagram in the rhs. plot of Fig. 5.

For the presentation of the results for the (pseudo-) critical line in the rhs. plot of Fig. 5 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 $\mu_{B}=3\mu$ in units of the critical temperature at that $\mu$, i.e., the combination $3\mu/T_{c}(\mu)$. The results of our determination of $T_{c}(\mu)$ from the maxima of the susceptibility are shown as asterisks in the rhs. plot of Fig. 5. With increasing $\mu$ 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 $\kappa$ defined via the relation (again we use $\mu_{B}=3\mu$)

[TABLE]

The fit of our data with this quadratic polynomial is shown as the full curve in the rhs. plot of Fig. 5. The small-$\mu$ data are described reasonably well and we obtain a value of $\kappa=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., $\kappa=0.0135(2)$ [29], $\kappa=0.0149(21)$ [30] and $\kappa=0.020(4)$ [31].

5 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 $\rho(x)$ with sufficient precision, such that one can reliably determine physical observables by integrating $\rho(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 $\langle{\mathrm{Im}\,\mathrm{P}}\rangle$ and ${{{\chi}}}_{\mathrm{Im}\,\mathrm{P}}$. Comparing the DoS results to a conventional simulation at $\mu=0$ shows good agreement and for $\mu>0$ the observables and the critical line could be determined up to moderately large values of $\mu$. A determination of the curvature $\kappa$ in the $\mu$-$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 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.

Acknowledgments: This work was supported by the Austrian Science Fund, FWF, DK Hadrons in Vacuum, Nuclei, and Stars (FWF DK W1203-N16) and we thank Kurt Langfeld, Biagio Lucini and Pascal Törek for discussions.