Developing and testing the density of states FFA method in the SU(3) spin model

The Density of States Functional Fit Approach (DoS FFA) is a recently proposed modern density of states technique suitable for calculations in lattice field theories with a complex action problem. In this article we present an exploratory implementation of DoS FFA for the SU(3) spin system at finite chemical potential $\mu$ - an effective theory for the Polyakov loop. This model has a complex action problem similar to the one of QCD but also allows for a dual simulation in terms of worldlines where the complex action problem is solved. Thus we can compare the DoS FFA results to the reference data from the dual simulation and assess the performance of the new approach. We find that the method reproduces the observables from the dual simulation for a large range of $\mu$ values, including also phase transitions, illustrating that DoS FFA is an interesting approach for exploring phase diagrams of lattice field theories with a complex action problem.


Introduction
An ab initio lattice simulation of QCD at finite density is one of the great current challenges in lattice field theory. The reason is that for many field theories the action S becomes complex at finite chemical potential µ and consequently the Boltzmann factor e −S does not allow for a probabilistic interpretation which is necessary for a Monte Carlo Simulation. Over the years various methods to overcome this so-called "complex action problem" have been developed. Examples are complex Langevin, the Lefshetz thimble, Taylor expansion around µ = 0, analytical continuation to imaginary chemical potential or the dual approach, where the system is rewritten to new variables such that all weights are real and positive. In the end, probably only the combination of the results from different approaches will lead to a reliable understanding of finite density lattice QCD.
Another approach are Density of States (DoS) methods which appear in the literature early on [1,2] and were revisited regularly, see, e.g., [3,4,5,6]. For the application of DoS techniques to finite density lattice field theory the main challenge is to compute the density of states ρ(x) with very high precision, since for the evaluation of observables it is integrated over with a highly fluctuating function.
An important new technique for pushing up the precision for the density, related to a proposal by Wang and Landau [7], was introduced by Langfeld, Lucini and Rago in [8,9,10,11,12]. The key idea of this so-called LLR method is to divide the variable x of the density ρ(x) into small intervals and to determine the variation of ρ(x) in each interval using restricted vacuum expectation values. These new techniques make it plausible that DoS methods could become competitive also for the analysis of finite density lattice QCD [13].
Here we present an exploratory study of a variant of DoS techniques, the Density of States Functional Fit Approach (DoS FFA) [14,15,16]. Similar to the LLR approach we parameterize the density on small intervals and evaluate restricted vacuum expectation values. These depend on a free parameter λ, and when using the density ρ(x) one can derive a closed form for their functional dependence on λ. With a one-parameter fit of the Monte Carlo data one then can precisely determine ρ(x) on the corresponding interval. The Dos FFA has been introduced and tested for the Z 3 spin system and was found to reproduce the results of a reference simulation with dual variables for a relatively wide range of parameters [14,15,16].
In this article we develop the DoS FFA further and apply it to a system with continuous degrees of freedom, the SU(3) spin model with chemical potential (the Z 3 model has discrete degrees of freedom). The SU(3) spin model is an interesting candidate theory since it is an effective model for the Polyakov loop and inherits the complex action problem from QCD. Maybe even more important is the fact that the model has a dual representation in terms of worldlines [17], such that only real and positive weights appear. Thus one can run Monte Carlo simulations directly in terms of the dual variables [18,19] and in this way generate reference data for assessing the efficiency and application range of DoS FFA.
In our exploratory study we present the generalization of the DoS FFA to continuous variables, discuss the implementation for the SU(3) spin model and use the reference data from the dual simulation for evaluating the new DoS approach. In addition we discuss various techniques that allow for an optimal use of the numerical resources in DoS calculations.
2 Definition of the model and the density of states Using lattice QCD as a starting point one can apply strong coupling expansion for the Wilson gauge action and a hopping expansion of the fermion determinant of the Wilson-Dirac operator at finite chemical potential to obtain an effective action for the Polyakov loop. The resulting model is the so-called SU(3) spin model with action where the dynamical degrees of freedom P( n) are the traces of SU(3) matrices representing the Polyakov loop, which live on the sites n of a 3-dimensional lattice Λ with periodic boundary conditions. The nearest neighbor coupling τ is an increasing function of the temperature T in the underlying lattice QCD theory. κ is a decreasing function of the quark mass and is proportional to the number of (mass degenerate) flavors. µ is the chemical potential, and we have absorbed the inverse temperature of the underlying lattice QCD such that our µ corresponds to the dimensionless combination µ/T of QCD.
The partition function is obtained by integrating the Boltzmann factor e −S[P] over all configurations of the dynamical degrees of freedom P( n) where the path integral measure is the product over the Haar measures on all sites n, where Since the action depends only on the trace of the SU(3) matrices we need only two angles θ 1 ( n), θ 2 ( n) ∈ [−π, π] to parameterize our degrees of freedom, and we can work with the reduced Haar measure (4) It is obvious that for finite chemical potential µ = 0, the action (1) has a non-zero imaginary part such that the Boltzmann factor e −S[P] has a phase and the model in the conventional representation has a complex action problem. However, as already mentioned in the discussion, the SU(3) spin model has a dual representation in terms of worldlines [17] where the complex action problem is solved and Monte Carlo simulations become possible for arbitrary µ [18,19]. These dual results will be used as reference data for the assessment of the DoS FFA.
The first step of a density of states method is to define the density to be used. Here we work with a weighted density and for its definition the action is divided into a real and an imaginary part, i.e., S[P ] = S R [P ] + iS I [P ]. Only the second part of (1) contains µ and we can write, κ n e µ P( n) + e −µ P( n) = 2κ coshµ n Re P( n) + i2κ sinhµ n Im P( n) .
We now explore symmetries of the system to further simplify the partition sum and expectation values. Charge conjugation C corresponds to the transformation ∀ n : θ i ( n) → −θ i ( n), i = 1, 2. Real and imaginary part of the action, and the reduced Haar measure transform as, such that the partition sum transforms as Thus we can write the partition function as Based on this form of the partition sum we now define the weighted density ρ, where we use as weight the Boltzmann factor with the real part of the action, Here we have already used, that X[P] is bounded such that x is restricted to the interval where V denotes the number of sites of or lattice. Using again the transformation properties under charge conjugation one easily finds that the density is an even function of x, i.e., ρ(x) = ρ(−x). Therefore, the partition function expressed in terms of the density is given by Expectation values of observables O(X[P]) which are a function of X[P] are given by where we again used the fact that ρ( The partition function (12) and the expectation values (13) are obtained by integrating the density ρ(x) with the factors cos(2κ sinhµ x) and sin(2κ sinhµ x). While the density ρ(x) is strictly positive, these factor are oscillating with x and the frequency of the oscillation increases exponentially with the chemical potential µ and linearly with the parameter κ. So for larger values of µ and κ the density has to be computed with very high accuracy. This is how the complex action problem manifests itself in the density of states approach.
3 Computing the density of states with the FFA For the numerical evaluation we need to parameterize the density ρ(x) in a suitable form. For that purpose we divide the interval [0, x max ] into N intervals [x n , x n+1 ], n = 0, 1, ... N−1 with x 0 = 0 and x N = x max . We stress that the intervals can differ in their size ∆ n = x n+1 − x n , but clearly the sum rule x n = n−1 j=0 ∆ j must hold. It will turn out that working with variable interval sizes allows one to use finer resolutions for the density in regions of x where ρ(x) shows strong variations.
We now parameterize the density as where L(x) is a continuous function which is piecewise linear, i.e., a straight line on each interval [x n , x n+1 ]. We stress at this point that the parameterization in the form (14) with a piecewise linear function in the exponent is only an approximation, since for continuous degrees of freedom ρ(x) cannot be represented with a finite number of parameters, which here are given by the slopes in the intervals 1 . Thus an exact representation can be obtained only in the limit ∆ n → 0, and below we will discuss a strategy how to determine an optimal choice for suitable interval sizes ∆ n . It is an interesting question how the interval sizes ∆ n affect the vacuum expectation values of observables. Obviously the ∆ n enter vacuum expectation values in a highly non-linear way and a thorough theoretical analysis of the effect of a finite discretization is beyond the scope of this study and has to be postponed to future work. However, in Section 5 we demonstrate numerically that already with moderate ∆ n we can reproduce the reference data from the dual simulation for a rather large range of chemical potential values.
Denoting the slopes of the straight line on the interval [x n , x n+1 ] by k n , one can work out the explicit representation of ρ(x), where we have fixed the (irrelevant) normalization of the density by setting ρ(0) = 1, i.e., L(0) = 0. It is obvious, that in our parameterization the density depends only on the slopes k n of the piecewise linear function L(x). For later use we have introduced the notation ρ(x) = A n e −knx in the second equality of (15), i.e., we have collected all factors independent of x into the constant A n .
To determine the density as given in (15) we need to find the slopes k n . For the calculation of the k n we use so-called restricted vacuum expectation values O n (λ), n = 0, ... N −1, which depend on a free parameter λ ∈ R. They are defined as with The support function θ n (x) restricts the values of X[P] in Z n (λ) and O n (λ) to the interval [x n , x n+1 ]. Furthermore we have also introduced an additional Boltzmann factor e λ X[P] , which is real and positive, such that the generalized vacuum expectation values O n (λ) can be evaluated with standard Monte Carlo techniques. The additional Boltzmann factor plays a twofold role: By varying λ we can systematically explore the density in all of the interval [x n , x n+1 ]. In addition we can evaluate explicitly the dependence of Z n (λ) on the parameter λ and in this way determine the k n . This second role will be outlined now.
Writing Z n (λ) with the density and using the fact that on the interval [x n , x n+1 ] the density is given by ρ(x) = A n e −knx (see (15)) we find, A simple calculation then gives (use x n = n−1 j=0 ∆ j ), which we rearrange to new, rescaled observables Y n (λ) defined as, where in the last step we introduced the abbreviation h(s) Using Eq. (20) we can now determine the slopes k n in a simple way: On the left hand side we have the observable Y n (λ) which is related to the restricted vacuum expectation value X[P] n (λ) which can be evaluated with standard Monte Carlo techniques for several values of λ. These Monte Carlo results are described by h (λ−k n )∆ n where h(s) is a simple smooth function. Thus on the right hand side of Eq. (20) we have a single free parameter, namely the slope k n which we want to determine. Its value is then obtained by a simple one-parameter fit of h (λ − k n )∆ n to the numerical data. This procedure is the reason for referring to our method as the "Functional Fit Approach" (FFA). Since we determine the k n from a fit of the data for Y n (λ), the normal 1/ √ N dependence of the errors on the statistics N applies.
It is straightforward to show that h(s) is monotonically increasing with h(0) = 0, h (0) = 1/12 and lim s→±∞ h(s) = ±1/2. Thus the function h (λ − k n )∆ n used for the fit has a single zero at λ = k n where it has a slope of ∆ n /12. In principle one could determine the zero of h (λ−k n )∆ n to find k n 2 , but using all values of λ in a fit with h (λ−k n )∆ n makes use of all generated Monte Carlo data. A second important aspect is that the quality of the fit also allows for a consistency check of the method: If one finds that the Monte Carlo data are not well described by the function h (λ − k n )∆ n this is an indication that the interval size ∆ n was chosen too large. We remark again, that we here work with variable interval sizes ∆ n = x n+1 − x n and in regions of x where the density ρ(x) shows a large variation one can choose smaller interval sizes ∆ n . In Fig. 1 we show an example of the fit procedure on a small lattice. The circles are the results of the evaluation of Y n (λ) with the restricted Monte Carlo simulations. We show the results for several values of n between n = 0 and n = 180 (circles of different colors; the data sets move to the right with increasing n). The full curves represent the functions   We find that preconditioning describes well the overall shape but the zoom on the right hand side shows that the preconditioning density is too coarse for evaluating observables.
h (λ − k n )∆ n where the k n were determined from a fit of the Monte Carlo data. It is obvious that the Monte Carlo data are very well represented by the functions h (λ − k n )∆ n . Note that for the largest two n we show (n = 160 and n = 180) a smaller ∆ n was used, which leads to a smaller slope ∆ n /12 at the point where h (λ − k n )∆ n crosses zero. The reason for choosing a smaller ∆ n is the fact that for large x the density varies faster such that a finer resolution, i.e., finer intervals [x n , x n+1 ] are advisable. Once the k n are obtained from the fits one can determine the density ρ(x) from the explicit expression (15). Examples for the density are shown in Figs. 2, 3 and 7. Let us finally comment on the choice of the values λ where the simulations for determining Y n (λ) are done. In Fig. 1 we observe that for larger n the curves shift to the right, which is due to the fact that for larger values of x, which corresponds to intervals [x n , x n+1 ] with larger n, the density ρ(x) has larger slopes k n . To obtain optimal fit results one should have Monte Carlo data in a range of λ that properly brackets the value λ = k n where the fit function h (λ − k n )∆ n has its zero.
To avoid exploring a large range of values for λ we developed a strategy which we refer to as "pre-conditioning". The idea is to do a first run with only a small number N of intervals, which then of course have rather large sizes ∆ n and give rise to only a crude approximation of ρ(x). However, analyzing this crude approximation we already have an idea what size the slopes k n are. This information about the k n can then be used for: 1. Determining the optimal interval sizes ∆ n (small for regions with a large variation of  ρ(x) and larger for regions with little change of ρ(x)).
2. Identifying suitable ranges for the values of λ in the determination of Y n (λ) such that the corresponding value λ = k n is properly bracketed by the Monte Carlo data.
In Fig. 2 we illustrate preconditioning and compare the final results for the density to a coarse approximation with only few and large intervals. The two plots (the rhs. is a zoom into the small-x range) show that the preconditioning results already give a good first estimate of the final density which can be used to optimize the runs for the target resolution with fine intervals [x n , x n+1 ]. On the other hand it is clearly seen in the zoom on the right hand side that the coarse density obtained in the preconditioning step is not suitable for computing observables with the fluctuating integrals of (13).
We remark at this point that the preconditioning idea can be pushed further: One can include the coarse density obtained from a preconditioning calculation into the definition of the density such that the final density is obtained as a correction to the preconditioning density. It is relatively simple, although somewhat technical, to work out the corresponding equations along the lines of Eqs. (11) -(20). One finds that the final density has much smaller values of k n since they are only corrections relative to the slopes of the preconditioning density. Thus the corresponding relevant ranges of λ are all centered near λ ∼ 0 and the effort for the determination of the optimal range for different n is considerably reduced.
It is clear that the shape of the density ρ(x) will be different for different parameter values. For a first illustration in Fig. 3 we show ρ(x) for fixed values of µ and κ but different values of τ . For a comparison of ρ(x) at different values of µ see Fig. 7 below.

Details of the restricted Monte Carlo simulations
Before we come to a more detailed discussion of our results, in this section we provide a short overview over the technical aspects of our simulations used in this paper, designed for assessing the DoS FFA method and to further develop this approach.
The DoS FFA method is based on fitting the Monte Carlo data for Y n (λ) which is related to the restricted expectation value X[P] n (λ) as defined in (16). Since the restricted vacuum expectation value requires X[P] ∈ [x n , x n+1 ], we first need to generate an initial configuration of the angles θ 1 ( n) and θ 2 ( n) such that this constraint is obeyed. Such a configuration with constant spin values P( n) for all n ∈ Λ can easily be constructed by hand, but of course needs to be equilibrated before taking measurements. For this and the subsequent computation of observables the restricted vacuum expectation values require a slightly modified Monte Carlo update. It contains an additional step which rejects proposal configurations P that violate the condition X[P] ∈ [x n , x n+1 ]. In the simulation of our SU(3) spin system we did not observe problems due to this additional rejection step and the acceptance rate remained reasonably high for all parameter values we tested. As we will see, with the parameters chosen here we can cover a rather large range of chemical potential values, but is is also clear that for pushing to even higher values one would have to use considerably smaller intervals and the additional rejection step might slow down the algorithm.
In addition to running the simulations for all N intervals, in each interval [x n , x n+1 ], n = 0, 1 ... N − 1 we have to produce data for a suitably chosen set of λ values such that we can perform an optimal fit of Y n (λ) with h (λ − k n )∆ n . Thus the approach requires a total of N × N λ restricted Monte Carlo simulations, where N λ is the number of λ values used for the fit. The number N of intervals of a given size ∆ n scales with the volume V of the lattice since X[P] is an extensive quantity, while N λ is a fixed number independent of the system size.
In this exploratory study we show data for lattice volumes of 8 3 and 12 3 . We use between N = 200 and N = 600 intervals and N λ ≈ 60 − 100. We use O(10 4 ) configurations for evaluating each X[P] n (λ) and these configurations are separated by 10 sweeps for decorrelation. All errors we display are statistical errors. For their determination we computed the error for each X[P] n (λ) via the standard deviation, determined the error for the k n from the fit with h (λ − k n )∆ n and subsequently performed a Monte Carlo propagation of this error. This procedure allows to simultaneously take into account the influence of the errors of all Monte Carlo data on the statistical error we show for the final observables.
As outlined in the previous section we use a preconditioning phase where we perform a simulation with only few and large intervals and low statistics. Typically we use N = 50 to N = 100 intervals and N λ ≈ 20 − 30 with a statistics of only O(10 2 ) configurations. However, this additional low-cost simulation considerably speeds up the final simulation since we can use suitably chosen interval sizes ∆ n and optimal ranges for the values of λ.

Observables and their comparison to a dual simulation
So far we have only discussed techniques how to compute the density ρ(x) and have shown some of the corresponding results. However, the true test of a density of states approach is of course the assessment of observables, since in their evaluation according to (13) the density is integrated over with highly oscillating factors that probe the fine details of ρ(x). Thus only observables test if the determination of ρ(x) is sufficiently accurate and delimit up to which values of µ the approach is reliable since, as we have discussed in Section 2, the frequency of oscillation increases exponentially with µ. It would be very interesting and important to perform a systematic comparison of accuracy and computational cost of physical observables to the results from other modern DoS techniques, in particular the LLR method [8,9,10,11,12], where, however, the focus so far was on technical development and on the density itself.
As pointed out in the introduction, one of the main reasons for choosing the SU(3) spin model as a test case for developing and assessing the DoS FFA method is the fact that the model has a dual representation [17] that solves the complex action problem. Dual Monte Carlo simulations in terms of worldlines [18,19] are possible for arbitrary µ and provide reference data for assessing the new DoS FFA techniques.
For the comparison to the dual simulation we use the particle number n and the corresponding susceptibility χ n . They are obtained as derivatives of ln Z and can be computed easily in the dual approach. Their exact definition and the representation in terms of the density ρ(x) are given by In the dual formulation one can perform the derivatives of ln Z in the same way and obtains the observables as moments of dual variables [17,18,19]. We begin the discussion of our results with Fig. 4, where we show n and χ n as a function of µ on a 8 3 lattice at κ = 0.005 and τ = 0.066. Here we use a relatively coarse discretization with N = 450 intervals and a low statistics of 5000 measurements for each X[P] n (λ). We find that with these parameters the density n (lhs. plot) from DoS FFA agrees well with the dual results for chemical potentials up to µ ∼ 3.0, while for the corresponding susceptibility we find good agreement up to µ ∼ 2.5. Note that for the second moment also the dual results show quite large errors due to long autocorrelation times of the dual simulation in this region of the parameter space.
An interesting question is how one should choose the overall scale for the interval size ∆ n . We already discussed that in regions of x with a large variation of ρ(x) one chooses   We compare the DoS FFA results to the reference data from the dual formulation.
smaller ∆ n , but we have not yet addressed the question what is a suitable choice for ∆ n in general. For computing observables according to Eq. (13) we integrate the density with sin(2κ sinhµ x) and cos(2κ sinhµ x). Thus, for a given κ and µ one full oscillation corresponds to ∆x = 2π/2κ sinhµ, which is the scale at which the fluctuating factors probe ρ(x). Clearly the resolution of ρ(x) given by the interval size ∆ n must be considerably smaller than ∆x and we obtain as criterion ∆ n 2π 2κ sinh µ . (23) In Fig. 5 we show results where we used a much higher statistics of 75000 measurements and chose intervals that are slightly smaller compared to those used in Fig. 4. Furthermore we here switched to τ = 0.130 (keeping κ = 0.005), a value where the system undergoes a crossover near µ = 2.0 (see [18]). Determining the observables across this crossover is clearly a more challenging task than the calculation at τ = 0.066 in Fig. 4, where there is no crossover in the range of µ values considered [18]. Nevertheless, with the improved statistics and smaller ∆ n we find very good agreement with the dual formulation all the way up to µ = 4.0 (although the dual data for χ n again suffer from large errors due to autocorrelation at these parameter values). It is interesting to note that the agreement remains good across the crossover near µ ∼ 2, which is visible in a maximum of n. Careful inspection of the error bars for the DoS FFA results show that they even become smaller again for µ > 2, which is an unexpected behavior that needs some extra consideration which we present below. For completeness we also considered the observables on larger lattices, and for the same parameters as in Fig. 5, i.e., κ = 0.005 and τ = 0.130 we show our results for lattice size 12 3 in Fig. 6. The behavior is almost exactly the same as on the smaller lattices used in Fig. 5.
Let us come back to the interesting observation we made in Figs. 5 and 6: For µ > 2 the error bars of the DoS data become smaller again. This is unexpected because we have pointed out that when increasing µ the frequency of the oscillating factors in (13) increases exponentially and the density is probed on even finer scales. The solution of this riddle is the fact that also the density ρ(x) depends on the chemical potential, since µ enters the real part of the action which is used for defining ρ(x) in Eq. (11). Indeed it turns out that when increasing µ above the crossover value µ ∼ 2 the shape of ρ(x) changes considerably. This is illustrated in Fig. 7 where we show the density ρ(x) at different values of µ for the parameter values used in Fig. 5. It is obvious that the behavior changes considerable when increasing µ above 2, and the density decreases much faster as a function of x for larger values of µ. This fast decrease suppresses contributions at large x and thus leads to faster convergence in (13) and thus smaller errors. We conclude that the applicability, accuracy and convergence range in µ of DoS methods not only depends on the parameters N , ∆ n and the statistics used, but also on the underlying physics.
We continue our discussion of observables with showing the results for a vertical section in the µ-τ plane. In Fig. 8 we work at fixed κ = 0.005 and fixed µ = 1.0 and study the observables n and χ n as a function of τ . In this case one crosses a first order transition near τ = 0.13 as determined from the Polyakov loop susceptibility in [18]. Although the volume is rather small in our study, it is still remarkable that also here we find very good agreement between the DoS FFA results and the reference data from the dual simulation. This indicates that DoS methods could indeed be competitive with other approaches also in the vicinity of phase transitions. We finalize our discussion of observables in Fig. 9, where we show e iS I pq , the phase quenched expectation value of the complex phase, which is a measure for the severeness of the complex action problem. The lhs. plot shows e iS I pq as a function of µ at κ = 0.005, τ = 0.13 for 8 3 and 12 3 , while the rhs. plot shows e iS I pq as a function of τ at κ = 0.005, µ = 1.0 for 8 3 . In the lhs. plot one nicely sees that the transitory behavior which we also observe in the corresponding observables shown in Figs. 5 and 6 is also reflected in the behavior of e iS I pq , with the complex action problem being more severe in the vicinity of the transition and then again for large µ. We also find that, as expected, the complex action problem becomes more severe when the volume increases. When considered as a function of τ (rhs. plot) we find that e iS I pq shows a monotonous behavior, matching the monotonous behavior seen in the corresponding observables in Fig. 8.

Truncating the density
It is obvious from Figs. 2, 3 and 7 that the weighted density ρ(x) drops over many orders of magnitude when increasing x. Although the amount of decrease depends on the parameters, it is an interesting question whether it is necessary to evaluate the density ρ(x) for all values of density is computed and the error for observables from truncation is negligible compared to errors from other sources, in particular from the discretization of x and from finite statistics.
The basis for experimenting with truncation are Eqs. (21) and (22) which make explicit that regions of x where ρ(x) is very small will contribute little to observables. In this section we explore the idea of truncation and show that the numerical cost can be reduced considerably without significant changes of the results. Let us begin with an estimate of the truncation effect for the partition sum Z. By x * we denote the value x where we truncate the density, i.e., we set ρ(x) = 0 for x > x * . The partition sum changes by an amount δZ which is given by A (very conservative) bound for this integral is given by |δZ| ≤ 2 x max × max x≥x * ρ(x).
Considering that for reasonable lattice sizes x max = V 3 √ 3/2 is of order O(10 3 ) − O(10 4 ) we conclude that when the density is smaller than, e.g., e −50 ≈ 10 −22 for all x > x * , then the contribution of the values x > x * is negligible. A density smaller than e −50 corresponds to ln ρ(x) = −L(x) < −50. This is the quantity we show in Figs. 2, 3 and 7 and it is obvious that for some parameter sets the value of −50 is reached already at relatively small x.
After this first estimate for the partition sum we now study the effect of truncation for observables numerically. More specifically, we define for the particle density n the error which is a simple generalization of δZ. It is straightforward to give an upper bound for δn in a similar way as for δZ. However, instead of discussing this bound, in Fig. 10 we show the numerical results for the relative truncation error δn/n as a function of x * /x max . The data are for µ = 2.0, k = 0.005 and τ = 0.066, and we find that when cutting at x * /x max = 0.2, i.e., we take into account only 20 % of the full range of x, the relative error has already dropped to 10 −46 . This plot impressively demonstrates that it is possible to considerably truncate the range of x values where the density is evaluated and obtain an error which is completely negligible in comparison to other error sources. This saves a lot of numerical resources which can be invested in reducing these other errors, e.g., by using smaller discretization intervals ∆ n or by increasing the statistics. We stress again that the functional form of the density changes when varying the parameters (see Figs. 3 and 7) and truncation can be applied only after checking whether the density is suitable. However, for that task we can again invoke preconditioning, i.e., the observation that a reasonable estimate for the overall behavior of the density can already be obtained from a first numerically inexpensive coarse estimate of ρ(x) using few and large intervals ∆ n .

Summary and outlook
Density of states techniques are an interesting approach to the simulation of lattice field theories that have a complex action problem, e.g., theories with a chemical potential µ. As discussed, the main challenge is to calculate the density of states ρ(x) with sufficient accuracy. When evaluating observables the density is integrated over with a highly oscillating factor and the frequency of the oscillation increases exponentially with µ. Thus, when one tries to reach reasonably large values of µ, very high precision for ρ(x) is mandatory.
In this paper we present the results of an exploratory implementation of the so-called Density of States Functional Fit Approach (DoS FFA) for the SU(3) spin model. The method uses a parameterization of the density ρ(x) by an exponential of a piecewise linear function. The slopes of the linear pieces are computed with restricted vacuum expectation values on the respective intervals. These restricted vacuum expectation values depend on a free parameter λ and can be computed with standard Monte Carlo techniques. In each interval the functional form of the restricted vacuum expectation as function of λ is known and depends only on the slope determining ρ(x) on that interval. The slope then is obtained via a simple oneparameter fit to all Monte Carlo data and we show that the quality of the fit can be used as a self-consistency check of the method. From the slopes one determines ρ(x) and with this density the observables. We introduce a strategy which we refer to as "preconditioning": we show that already a coarse parameterization and low statistics are sufficient to get a good estimate for the overall behavior of the density. This information can then be used to optimize the interval size ∆ n for the piecewise linear parametrization of the exponent of ρ(x), to determine suitable ranges for the values of λ to be used, and to assess a possible truncation of the density, a step which we illustrate to have a great potential for saving computer time that can be better used to reduce the error from discretization effects of ρ(x) or from finite statistics.
We use the dual formulation of the SU(3) spin model for generating reference data for assessing the DoS FFA method, which is possible because the dual formulation is free of the complex action problem. In particular we study the particle number density and the corresponding susceptibility and find that the DoS FFA results match the dual simulation reference data for a surprisingly large range of µ. Part of the reason for this success is the fact that the shape of the density changes in a crossover transition that is crossed when increasing µ, such that ρ(x) drops faster with x, a feature that is beneficial for the accuracy of the DoS FFA.
Encouraged by the success of the DoS FFA in the SU(3) spin model and the Z 3 model, studied in an earlier paper, we are currently working on developing the method further towards QCD. The essential steps will be the change from SU(3) spins to SU(3) gauge fields, and the formulation and implementation for theories with fermions. For both these issues we make progress and expect that the DoS FFA can be formulated for non-abelian gauge fields interacting with fermions as a generalization of the implementation presented here.