Density of states method for the Z(3) spin model

We apply the density of states approach to the Z(3) spin model with a chemical potential mu. For determining the density of states we use restricted Monte Carlo simulations on small intervals of the variable for the density. In each interval we probe the response of the system to the variation of a free parameter in the Boltzmann factor. This response is a known function which we fit to the Monte Carlo data and the parameters of the density are obtained from that fit (functional fit approch; FFA). We evaluate observables related to the particle number and the particle number susceptibility, as well as the free energy. We find that for a surprisingly large range of mu the results from the FFA agree very well with the results from a reference simulation in the dual formulation of the Z(3) spin model which is free of the complex action problem.


Introductory comments:
It is well known that at finite density Monte Carlo simulations of many lattice field theories are plagued by the complex action problem (or sign problem). At non-zero chemical potential µ the action S acquires an imaginary part and the Boltzmann factor e −S cannot be used as a probability weight. Several different strategies to overcome the complex action problem have been proposed over the years, among them the density of states (DoS) method, which is the topic of this letter.
Originally the DoS method was introduced to lattice field theories in [1] and for some recent applications and a critical review see, e.g., [2,3]. The key problem of DoS techniques is that the density of states ρ varies over many orders of magnitude. When applying DoS to finite-µ problems the density ρ is integrated over with a highly oscillating factor and the frequency of the oscillation increases exponentially with µ. Thus the numerical challenge is twofold: The density ρ has to be determined over many orders of magnitude and this determination has to be very precise since ρ is probed with a highly oscillating factor.
An interesting new approach to the accuracy problem in the DoS for lattice field theories was presented in [4]. The idea is related to a proposal by Wang and Landau [5] and divides the variable of the density ρ into small intervals where restricted Monte Carlo simulations are performed to determine ρ in that interval. This leads to exponential error suppression [4] and was shown to allow for a precise evaluation of observables in pure gauge theory and in SU(2) gauge theory with a Polyakov loop source [4]. Furthermore, in [6] the density and the complex phase were studied in the Z 3 spin model with chemical potential, which is a simple effective theory for the center degrees of freedom of QCD [7]. This model is particularly interesting since it can be mapped to a dual representation without complex action problem where precise Monte Carlo simulations at finite µ are possible [8].
In this letter we develop the DoS for the Z 3 spin model further (functional fit approach; FFA), using again the dual simulation [8] as our reference data. We evaluate observables related to the quark number density, the quark number susceptibility as well as the free energy and show that the corresponding dual results can be reproduced for a surprisingly large range of the chemical potential. First results with the FFA were presented in [9].

Definition of the model and the density of states:
The Z 3 spin model in an external field with strength κ, a chemical potential µβ (in units of the underlying QCD inverse temperature β) and an effective temperature parameter τ is described by the action where the dynamical degrees of freedom are the spins P x ∈ Z 3 = {1, e i2π/3 , e −i2π/3 }, living on the sites x of a 3-dimensional lattice with periodic boundary conditions. The action is normalized such that S[P ] = 0 if P x = 1 ∀x. The partition function is obtained by summing the Boltzmann factor over all configurations P , i.e., Z = {P } e S[P ] . The first and second derivatives of ln Z with respect to µβ have the interpretation of the quark number density and the quark number susceptibility. It is obvious that for finite chemical potential, µβ = 0, the action (1) has a non-zero imaginary part and the model has a complex action problem. For a convenient notation we introduce abbreviations for the total numbers of spins pointing in each of the three possible directions as N 0 [P ] = x δ (P x , 1) and N ± [P ] = x δ P x , e ±i2π/3 . Obviously N 0 [P ] + N + [P ] + N − [P ] = V , where V denotes the lattice volume, i.e., the total number of sites of the lattice. Using these we split the action into real and imaginary parts, i.e., where We now define a weighted density of states (δ here denotes a Kronecker-delta), Exploring again the symmetry properties of the action one trivially finds that ρ(d) is an even function of d, i.e., ρ(−d) = ρ(d). Using the density of states, the partition sum (3) can be written as Expectation values of observables O(∆N ) which are a function of ∆N are given by where O E and O O denote the even and odd parts of O(∆N ). The partition sum (5) and the expectation values (6) are obtained by reweighting the density ρ(d) with the factors cos κ √ 3 sinh(µβ) d and sin κ √ 3 sinh(µβ) d . While the density ρ(d) is strictly positive, these factor are oscillating with d and the frequency of oscillation increases exponentially with the chemical potential µβ and linearly with the strength parameter κ. Thus for larger values of µβ (or κ) the density ρ(d) has to be computed very accurately. This is how the complex action problem manifests itself in the density of states approach.

Computing the density of states :
For the numerical computation we parameterize the density of states ρ(d) as with real parameters a j . Note that this parameterization is exact in the sense that it contains V + 1 parameters, precisely the number of independent degrees of freedom ρ(d) has (remember that ρ(d) is an even function). We also remark, that an overall normalization of ρ(d) can be chosen freely, since it cancels in the expectation values (6). Here we choose the normalization ρ(0) = 1, which corresponds to setting a 0 = 0. For the calculation of the coefficients a j we define restricted expectation values O n (λ), n = 0, 1, ... V − 1, which depend on a free parameter λ, [P ] . (8) Here we have defined Varying the parameter λ in (8) probes the density in the interval set by θ n (d) and the response of the system to changing λ can be used to determine the parameters of ρ(d) in that interval. In the restricted expectation values (8) only real and positive weight factors appear, such that they can be evaluated with a restricted Monte Carlo strategy which we will discuss below.
In particular we are here interested in the observable O = ∆N , and now use the density of states ρ(d) in the form of (7) to evaluate the restricted expectation values ∆N n (λ). A straightforward calculation gives The right hand sides are simple functions of λ: They are monotonically increasing (the derivatives with respect to λ are easily shown to be positive) and for n ≥ 1 have a single zero (lim λ→−∞ is negative, lim λ→+∞ is positive). Examples for different n are shown in Fig. 1.
Using Monte Carlo simulations we can evaluate ∆N n (λ)−n for different values λ i , i = 1, 2, ... N λ (typically N λ = O(10)) and fit the results according to the right hand sides of (10). The one-parameter fit for ∆N 0 (λ) determines the first non-trivial coefficient a 1 (remember that we chose the normalization a 0 = 0). The fit value for a 1 can then be inserted in the right hand side of (10) for n = 1 such that with another one-parameter fit of ∆N 1 (λ) − 1 we can determine a 2 , which in turn is then inserted in the fit function for ∆N 2 (λ) − 2, which gives a 3 from a one-parameter fit, et cetera. Using this sequence of fits we can determine all coefficients a j from fits of the Monte Carlo data with simple functions that depend on only a single parameter (compare Fig. 1). Thus we refer to our approach sketched in this letter as "functional fit approach" (FFA).
We expect that this method has smaller statistical errors for the same numerical effort than iteration methods such as the LLR [4,6] or some general root finding procedure for the ∆N n (λ)−n (which also provides iterative equations for the a j ). The advantage of the FFA comes from the fact that all Monte Carlo data, i.e., the results for ∆N n (λ i ) − n at all values λ i of the parameter λ, are used to determine the coefficients a j . Furthermore, it can be seen that possible instabilities from the statistical errors of the Monte Carlo data are minimized here.

Restricted Monte Carlo:
The FFA variant of the DoS method described here is based on fitting the Monte Carlo data for the restricted expectation values ∆N n (λ) as defined in (8). For this purpose we first need to generate an initial configuration P of the spin variables, such that the constraint ∆N [P ] ∈ {n − 1, n, n + 1} is obeyed 1 . Such a configuration can easily be constructed by hand, but of course needs to be equilibrated before taking measurements. For this and the subsequent computation of observables a slightly modified Monte Carlo update can be used. It contains an additional restriction which rejects trial configurations that violate the constraint ∆N [P ] ∈ {n − 1, n, n + 1}. The acceptance rate is very good throughout and only for n very close to the maximum value of n = V (i.e., the cases n = V − 2, V − 1, V ) we observe a drop in the acceptance rate. In principle it is easy to compute ρ(d) for these largest values of d exactly with a low temperature expansion. However, since for the values of d where the quality of the Monte Carlo data decreases ρ(d) is already very small, we simply use the data as we obtain them from the simulation. In this letter we show results for lattice volumes of 10 3 and 16 3 and in both cases used 10 6 equilibration sweeps and a statistics of 10 6 measurements separated by 100 sweeps for decorrelation, and all errors we display are statistical errors. In Fig. 1 we show the Monte Carlo results (symbols) for ∆N n (λ) − n with n = 0, 2, 10, 50, 100, 300, 600, 900, 950, 990, 995 and 998 for several values of λ in the interval [−7, 7]. The data were generated on 10 3 lattices for the parameter values τ = 0.16, κ = 0.01 and µβ = 1.0. The figure demonstrates that the Monte Carlo data show the expected simple behavior as a function of λ and can easily be fit (we use a standard χ 2 procedure) with the functions given in the right hand sides of (10). In Fig. 1 the results of the fits are shown as full curves and obviously describe the numerical data very well.
Once the coefficients a j are determined from the fits, we can build up the density of states ρ(d) as given in (7). Results for the density ρ(d) at different values of µβ are presented in Fig. 2. The data we show in the lhs. plot are for 16 3 lattices with τ = 0.16, κ = 0.01, while in the rhs. plot τ = 0.178, κ = 0.001 were used. It is remarkable that the range of the values for ρ(d) strongly depends on the parameters, including also µβ. This is due to the fact that the weighted density we use here also includes the Boltzmann factor e S R .

Results for physical observables:
Having determined the density ρ(d) we can finalize the calculation and evaluate the expectation values of observables using (6). We consider M − M * and Thus for the evaluation of M − M * only an odd part appears in (6) and thus this expectation value is real.
In Fig. 3 we show the results from the FFA for M − M * and χ M −M * on a 16 3 lattice as a function of µβ (red circles). We display results for two sets of parameters, τ = 0.16, κ = 0.01 on the lhs., and τ = 0.178, κ = 0.001 (rhs.). The data are for 16 3 lattices and the statistics is the same as used for the density discussed above. As reference data in Fig. 3 we also show the results from a dual simulation. For the τ = 0.178, κ = 0.001 data (rhs.) the FFA results agree well with the dual results for all values of µβ we show. However, for τ = 0.16, κ = 0.01 (lhs.) the dual and the plain FFA data disagree for µβ larger than 2. This discrepancy can be attributed to the fact, that for this parameter set the sign problem is much harder than for τ = 0.178, κ = 0.001, as can, e.g., be seen in plots for the expectation value of the phase of the action [10].
The numerical problems can be related to small statistical fluctuations of ρ(d) around its exact values [6]. However, it is straightforward to smoothen these local fluctuations by using a fit of ρ(d) and then computing the observables (6) with the fitted ρ(d). For the fit we use a polynomial in d 2 (note that ρ(d) is an even function and ln ρ(0) = 0) for the logarithm of ρ(d), i.e., ln ρ(d) = N n=1 c n d 2n . In Fig. 3 we also display the results obtained from the fitted ρ(d) with N = 15 (triangles) and find that we obtain a much larger range in µβ where FFA and dual simulation agree also for the τ = 0.16, κ = 0.01 data.
We stress that using a fit to smoothen ρ(d) is of course not a fundamental ingredient of the FFA, and increasing the statistics when determining ρ(d) will also improve the results. However, the source of the error is very clear: For large µβ the density ρ(d) is probed by the rapidly oscillating factors in (5) and (6) and the small fluctuations in ρ(d) become dominant. On the other hand ρ(d) is a smooth function, such that a fit is a much more cost efficient method than a drastic increase of the statistics.
As a further observable we consider the free energy F (µβ) defined as F (µβ) = ln(Z(µβ)/Z(0)). In Fig. 4 we show the results for the free energy at two values of the couplings and again compare the FFA method with data obtained from a dual simulation. As for the other observables we find excellent agreement of the FFA results with the data from the reference simulation in the dual approach. It is remarkable that for the free energy this excellent agreement is achieved already without fitting the density, which was necessary for the bulk observables. lattices as a function of µβ. We use two sets of parameters, τ = 0.16, κ = 0.01 on the lhs., and τ = 0.178 and κ = 0.001 on the rhs. We show results from the FFA algorithm, the FFA algorithm combined with a fit of ρ(d) and for comparison also the results from a simulation in the dual representation.

Summary:
The key issue in the application of the density of states method is to determine the density ρ(d) as precise as possible, since in the expectation values (6) ρ(d) is summed over with a highly oscillating function. Thus it is necessary to optimize the strategy for the computation of ρ(d) in every possible way. In these notes we test a strategy (the functional fit approach (FFA)) where in restricted Monte Carlo simulations on small intervals on d the response of the system to a free parameter in the Boltzmann factor is evaluated. The response is given by a known function which we fit to the Monte Carlo data to determine the parameters of the density ρ(d).
The results from the DoS calculation for M − M * , χ M −M * and the free energy F (µβ) are compared to reference data obtained in a dual simulation. We show that when using a fit for the density of states ρ(d) the observables agree very well with the dual results on a surprisingly large range of µβ values up to µβ ≈ 4.