New DoS approaches to finite density lattice QCD

We present two new suggestions for density of states (DoS) approaches to finite density lattice QCD. Both proposals are based on the recently developed and successfully tested DoS FFA technique, which is a DoS approach for bosonic systems with a complex action problem. The two different implementations of DoS FFA we suggest for QCD make use of different representations of finite density lattice QCD in terms of suitable pseudo-fermion path integrals. The first proposal is based on a pseudo-fermion representation of the grand canonical QCD partition sum, while the second is a formulation for the canonical ensemble. We work out the details of the two proposals and discuss the results of exploratory 2-d test studies for free fermions at finite density, where exact reference data allow one to verify the final results and intermediate steps.


I. INTRODUCTION
Finding a suitable approach for Monte Carlo simulations of finite density QCD that extends the accessible part of the QCD phase diagram towards larger values of the chemical potential is currently one of the great challenges for lattice field theory. The problem is that at finite chemical potential the fermion determinant is complex and cannot be used as a probability in a Monte Carlo process. Among the approaches that have been explored to solve this so-called complex action problem are density of states (DoS) techniques, introduced to lattice field theory in [1,2]. The key challenge for a DoS approach is to compute the density with sufficiently high accuracy, such that it can be integrated over with fluctuating integrands that appear when evaluating observables at finite density. A naive determination with, e.g., simple histogram techniques turned out to be useful only for very low densities [3][4][5][6].
Inspired by the Wang-Landau [7] approach an interesting new development was presented by Langfeld, Lucini and Rago [8][9][10][11][12][13][14]. The idea is to use a parameterization of the density as the exponential ρ(x) = exp(−L(x)) of a piecewise linear and continuous function L(x). Vacuum expectation values restricted to the intervals where L(x) is linear are used to determine ρ(x) with very high precision. A variant of the Langfeld-Lucini-Rago method is the DoS Functional Fit Approach (FFA) [15][16][17][18][19] and in recent years both techniques were used to obtain interesting results for several bosonic lattice field theories at finite density (see, e.g., the review [20]).
However, no modern DoS formulation for systems with fermions was presented so far, and thus no clear path towards precise DoS calculations for finite density QCD has been outlined yet. The challenge is to formulate the DoS approach such that it is compatible with conventional pseudo-fermion Monte Carlo techniques that may be applied to a real and positive fermion determinant. * Member of NAWI Graz.
In this paper we discuss two proposals how to implement the DoS FFA for finite density lattice QCD. The first of the two is based on using a suitable pseudofermion representation of the QCD grand canonical partition sum. The imaginary part of the action is identified and the density is considered as a function of that imaginary part. DoS FFA is used to determine the corresponding density and observables are then obtained as integrals of the density.
The second proposal implements DoS FFA in a canonical setting. The canonical partition functions at fixed net quark number are written as the Fourier moments with respect to imaginary chemical potential µ = iθ/β. Considering θ as an additional degree of freedom in the path integral allows one to implement the DoS FFA and compute the density as a function of θ. Observables at fixed net quark number are then again obtained as integrals of the density.
In both formulations only Monte Carlo simulations without sign problem are needed, which furthermore can be implemented using well established techniques of standard lattice QCD simulations. We work out the details of the two new approaches and present the results of small exploratory 2-d studies of the free case where exact results can be used to assess the results and intermediate steps of the new proposals.

II. GENERAL FORMULATION OF THE DOS FFA APPROACH
Before we can discuss our two new DoS approaches to finite density QCD we first need to discuss the details of the DoS formulation we use, the Functional Fit Approach. The FFA [15][16][17] is here presented for a general bosonic theory with a complex action problem and we will later show that with suitably chosen pseudo-fermion representations finite density lattice QCD can be brought into the general form introduced in this section.

A. Partition sum and density of states
The vacuum expectation values O for some observable O that we consider here can be written as bosonic path integrals, where Φ denotes an arbitrary set of general bosonic lattice fields that can be based on sites or links, and D [Φ] is the corresponding product measure. We have already separated the exponent of the Boltzmann factor into two terms, the real part S R [Φ] of the action and the imaginary part αX [Φ]. We have allowed for a real-valued coupling α ∈ R multiplying the imaginary part, which is useful in some of the applications we have in mind.  (2) where J [Φ] is an arbitrary real and positive functional of the fields. Below we will identify J [Φ] with some observable which in general can be decomposed into pieces that obey these requirements. Note that different choices of J [Φ] result in different densities ρ (J ) (x) and we use a superscript J to indicate which density we refer to.
With the densities ρ where in the second line we have explicitly listed also the particularly simple case where the observable is some The range of integration for the integrals dx in (3) depends on the properties of the imaginary part X [Φ]. If X[Φ] is bounded by some number x max , so is the integration range. We will see that this is the case for the canonical DoS formulation discussed in Section 4. Furthermore, usually one can identify symmetries to show that the densities ρ (J ) (x) are even or odd (depending on J ), such that the integration interval where we need to determine the densities is [0, x max ].
In case X[Φ] is unbounded the integration runs up to x = ∞, which is the case we will encounter in the direct DoS approach discussed in Section 3. Again symmetries can be used to show that ρ (J ) (x) is even (or odd) such that the actual integral is ∞ 0 dx. Furthermore we will see that the densities ρ (J ) (x) quickly decrease with x such that the range of integration can be truncated such that also in this second case we need to determine the densities in an interval x ∈ [0, x max ].
Having defined the densities ρ (J ) (x) and expressed observables as integrals over these densities we now have to address the problem of finding a suitable representation of the densities and how to determine the parameters used in the chosen representation.

B. Parametrization of the density
The densities ρ (J ) (x) are functions of the parameter x and we are interested in the densities in some finite interval [0, x max ]. For parameterizing the densities, we divide the interval [0, x max ] into N subintervals as follows: The densities ρ where the L (0) = 0 we can completely determine the constants a n as functions of the slopes k n . A simple calculation shows that L (J ) (x) can be written in the following closed form, from which we obtain the explicit form of the density ρ Thus our parameterized density ρ (J ) (x) depends only on the set of slopes k (J ) n , one for each of the intervals I n . We point out that the parametrization allows one to work with intervals I n that have different sizes ∆ n . In particular in regions where the density ρ (J ) (x) varies quickly one should use smaller intervals, while in regions of slow variation larger ∆ n can be used to reduce the computational cost. For a coarse scan of the density ρ (J ) (x) with the goal of determining the regions of quick variation, one can do a first numerically cheaper determination with large ∆ n which subsequently is refined with finer intervals. These techniques are referred to as preconditioning and are discussed in detail in [15][16][17].

C. Evaluation of the density parameters with FFA
To determine the density, we need to compute the slopes k with the corresponding restricted partition sums Z (J ) n (λ) given by where we have introduced the support functions Θ n (x) = 1 for x ∈ I n , 0 for x / ∈ I n .
In the restricted expectation values X (J ) n (λ) and the partition sums Z (J ) n (λ) we have introduced a free real parameter λ which couples to the imaginary part X [Φ] and enters in exponential form. Varying this parameter allows one to properly explore the x-dependence of the density in the whole interval I n . The expectation values X  In the first step we have rewritten the restricted partition sum as the integral of the density ρ (J ) (x) over the interval [x n , x n+1 ]. In the second step the parameterized form (9) was inserted for that particular interval, which gives rise to a simple integral of an exponential that can be evaluated in the closed form on the right-hand side.
Comparing (10) and (11) it is obvious that the restricted vacuum expectation value X (J ) n (λ) can be computed as the derivative X (J ) n (λ) = d ln Z (J ) n (λ)/dλ, such that we find the closed expression, After multiplicative and additive normalization we can express the result for X (J ) n (λ) (which in its normalized form we denote as V (J ) n (λ)) in terms of a function h(s), and has the properties The strategy for determining the slope k (J ) n for an interval I n now is as follows: Using a standard Monte Carlo simulation without sign problem we compute the restricted vacuum expectation value X

III. DIRECT DOS APPROACH FOR LATTICE QCD WITH A CHEMICAL POTENTIAL
In this section we discuss the first of our two implementations of the new DoS approach to finite density QCD.
Here we use a suitable pseudo-fermion representation of the grand canonical partition sum and separate the part with the complex action problem. For this factor we set up the DoS FFA formulation, discuss its properties and present results of a first exploratory test in the free case.

A. Grand canonical partition sum and pseudo-fermion representation
We consider lattice QCD with N f mass-degenerate flavors of Wilson fermions. After integrating out the fermions the corresponding grand canonical partition sum with quark chemical potential µ is given by We consider the theory in d = 2 and d = 4 dimensions using lattices of size V = N d−1 s × N t . The SU(3)valued gauge variables U ν (x) live on the links (x, ν) of the lattice and obey periodic boundary conditions. Their path-integral measure is the product of Haar measures is the Wilson gauge action (we dropped the constant additive term), β g is the inverse gauge coupling and P [U ] the sum over the real parts of the traced plaquettes. By D[U, µ] we denote the Wilson Dirac operator with chemical potential µ in the background of a gauge field configuration U . We write the Dirac operator in the form, with the matrix elements of the hopping terms given by By γ ν we denote the Euclidean γ-matrices in d = 2 or d = 4 dimensions, and κ is the hopping parameter κ = 1/(2d + 2m) with m the bare quark mass. To be specific, we use a representation of the Euclidean γ-matrices where γ d is symmetric, which in d = 4 is, e.g., the chiral representation and in d = 2 the choice γ 1 = σ 2 , γ 2 = σ 1 with γ 5 = σ 3 . In (21) we use matrix/vector notation for the d Dirac indices of the γ-matrices and the 3 color indices of the link variables U ν (x). The chemical potential gives different weight for hopping in forward and backward temporal direction, i.e., the ν = d direction. The fermions obey periodic boundary conditions in the spatial direction(s) and anti-periodic boundary conditions in time, i.e., the terms in (21) that connect sites with x d = N t − 1 and x d = 0 have an additional minus sign.
In order to introduce a pseudo-fermion representation that is suitable for the DoS FFA we write the fermion determinant as In the second step we have written 1/ det D[U, µ] † as a bosonic integral over a complex-valued scalar field Φ(x) with 3d components for Dirac and color degrees of freedom, and in the exponent we use vector/matrix notation for all indices. The constant C is given by In the third step we have organized the exponent into real and imaginary parts such that the pseudo-fermion integral matches the general form introduced in (1), where here we set α = 1. The corresponding real and imaginary parts are given by where we also write the gauge field U as an argument in These two matrices are defined as A straightforward evaluation gives the matrix elements,  (23) are real. Thus the pseudo-fermion integral in (22) has the form that allows one to use DoS FFA to evaluate that integral. This will be discussed in more detail in the next section.
Let us add a few comments on the first factor in (22), i.e., the determinant det D[U, µ] † D[U, µ] . Using the well known generalized γ 5 -hermiticity property corresponds to the fermion determinant of two mass-degenerate quark flavors with an isospin chemical potential which is free of complex action problems. We stress, however, that this isospin determinant is of course only a part of the weight and its coupling to the pseudo-fermion factor in (22) generates the full dynamics (see also the comments below).
The matrix D[U, µ] † D[U, µ] is obviously hermitian and has real and non-negative spectrum, such that it is directly accessible with pseudo-fermion methods. Possible approaches are a direct pseudo-fermion representation (below χ and χ j denote bosonic complex-valued pseudofermion fields), or an order-n Chebychev multi-boson representation [21,22] of the form where u j = e i2πj/(n+1) are the coefficients for the Chebychev factorization and we have used (20) can be obtained as follows: Using the triangle inequality one finds Using the definition (21) (32) Thus we find that there is a finite range of µ where the factor det D[U, µ] † D[U, µ] which is free of the complex action problem can be treated with conventional pseudofermion techniques. We stress again, that the estimate (32) is only a crude non-exhaustive bound that essentially reflects the situation of the free case, where condensation sets in at µ = m. For a dynamical background gauge configuration U the spectrum of H[U, µ] is known to contract such that values of µ that exceed the bare quark mass parameter m become accessible. To precisely delimit the range where the pseudo-fermion treatment of det D[U, µ] † D[U, µ] is possible beyond the bound (32) obviously is a dynamical question that has to take into account the emerging finite density physics as well as possible numerical instabilities of the HMC algorithm that can only be assessed in a full QCD simulation, which clearly goes beyond the scope of this presentation. However, already with the simple bound (32) we have established an interesting minimal region where the direct DoS approach is applicable in principle.

B. Implementation of the DoS FFA
To set up the density of states approach and to define the corresponding densities as outlined in the general presentation in Section 2 we need to write the grand canonical definition with the pseudo-fermion representation. Since we consider the general case of N f flavors, we need N f copies of the pseudo-fermion fields, Φ j , j = 1, ... N f , where by {Φ} we denote the set of all these fields. Based on the discussion of the previous section we thus write the grand canonical partition sum of QCD in the form that matches Eq. (1) with α = 1 (irrelevant overall constants were dropped), i.e., where the real and imaginary parts of the pseudo-fermion action, as well as the path integral measure were gener-alized to N f flavors, As we have outlined in the previous section the term (22) can be treated with conventional pseudo-fermion techniques and we combined the corresponding factor for N f flavors together with the gauge field action Following the general DoS FFA strategy outlined in Section 2.1 we now define the densities as where we allow for general observables J [{Φ}, U ] that can be functionals of both, the set {Φ} of pseudo-fermion fields Φ j , as well as the gauge fields U .
In the general outline of the method in Section 2 we have already announced that symmetries can be used to establish that the densities ρ (J ) (x) are either even or odd functions, depending on the observables J , which we assume themselves to be even or odd (general J may be decomposed into even and odd pieces). As an example we briefly discuss the simplest case of J = 1 and show that ρ (1) (x) is even. The symmetry transformation we consider is charge conjugation that for the gauge links and the pseudo-fermion fields is implemented as where * denotes complex conjugation and T transposition. It is straightforward to show that Equally straightforward is to show that the gauge action S g [U ] defined in (19) is invariant under the charge conjugation transformation (37), i.e., S g [U ] = S g [U ]. The invariance of the factor det D[U, µ] † D[U, µ] can be shown using the representation (27) and charge conjugation: Denote by C the charge conjugation matrix that where in the first step we used (27) Thus we have established that ρ (1) (x) is an even function and in a similar way one may analyze the symmetry properties for the general densities ρ (J ) (x) that contain the insertion of some observable J .
As the final step for the implementation of the DoS FFA we need to identify the restricted expectation values as defined in the general description of the method in Section 2. Comparing with the general form (2), (10) we can read off from the densities (36) the necessary restricted expectation values for full QCD, where in the second step we have written the combina- We stress that the ensemble considered in the restricted vacuum expectation values (41) is not simply QCD with isospin chemical potential reweighted to quark chemical potential, where a serious overlap problem would emerge. Instead the exponent of the Boltzmann factor in (41) is given by of the pseudo-fermion terms, which contribute the dynamics of the quark chemical potential.

C. First tests for the free case
For a first exploratory study of the new DoS approach we analyze the free case in two dimensions. Note that due to the restricted expectation values that need to be evaluated, this analysis already requires Monte Carlo simulations also for the free case and indeed provides a nontrivial test of the method. Insight about suitable sizes ∆ n for the intervals, the numerical cost, the accuracy that is needed for the density et cetera, can be obtained. Furthermore, the free case allows for a systematical comparison of the final results and the intermediate steps against analytical results that may be computed with Fourier transformation.
For the free case the density ρ (1) (x) defined in (36) with the help of the pseudo-fermion representation simplifies to (we consider the case of N f = 1 flavor) The gauge field integration has been dropped for the free case and also the Boltzmann factor (35) for the effective action since it is independent of x, such that it only would affect the overall normalization of the density which is set by requiring ρ (1) (0) = 1. Following the steps of the implementation of DoS FFA in the previous section, for determining the parameters k (1) j we need to evaluate the restricted expectation values defined in (41) which for the free case reduce to (44) The imaginary part X[Φ] can be obtained from (23) and (25) (drop the link variables U ν (x) there) as, and the kernel M[λ] in the Boltzmann factor of (44) is given by (see (42)), This matrix is obviously hermitian such that its eigenvalues are real. However, for the existence of the path integral needed for the evaluation of the restricted vacuum expectation values (44), the eigenvalues also have to be positive, and we now address this issue that has been neglected previously. It is easy to see that eigenvalues can become negative for large values of the parameter λ. The result for such an analysis is shown in Fig. 1, where we plot the value λ max as a function of µ/m and compare our results for different masses m and lattice sizes L × L. The value of µ/m where λ max becomes zero signals the breakdown of the method. We remark that the polygon-like behavior of the curves for the smaller volumes reflects the fact that for small volumes the momenta populate the interval [−π, π] with only a few values such that also a "sparse" spectrum emerges, and the different sections of the "polygon" correspond to a different eigenvalue becoming negative.
The data in Fig. 1 is organized in groups where in each group we consider a sequence of values L → ∞ and m → 0 at a fixed value of mL, i.e., we study the fixed volume continuum limit of the free theory. The dotted curves are for mL = 0.64, the dashed curves for mL = 1.28 and the full curves correspond to mL = 2.56. Note that the curves for different mL cluster according to the respective values of m. The figure shows that with increasing µ/m the values for λ max decrease and at a critical value of µ/m the boundary λ max becomes zero, signaling the breakdown of the method. We observe that for all three values of mL we study, the critical value of µ/m converges from above to a critical value of µ/m = 1, which is the value of the chemical potential where condensation sets in. Thus we expect that we can use DoS FFA all the way to the condensation point.
For the dynamical case one expects a similar behavior: For non-trivial gauge links the spectrum of the Dirac operator is known to contract, giving rise to an additive renormalization of the mass and a critical κ that is larger than the free value κ = 1/2d. Qualitatively one finds that for a larger critical κ a larger value of µ is accessible, and one may expect that also for the full case the critical value of µ coincides with the point where condensation sets in. We stress, however, that obviously this is only a very qualitative discussion of the situation in the fully dynamical case and future explicit Monte Carlo calculations will be necessary for a detailed analysis.
Having identified non-vanishing windows of λ where we can safely evaluate the restricted vacuum expectation values X (1) n (λ) defined in (44), we show some of these results for illustration in Fig. 2. We plot the restricted vacuum expectations X (1) n (λ) normalized to the form V (1) n (λ) defined in (15) as a function of λ. The symbols represent the data that we determined in a small Monte Carlo simulation on a 16 × 16 lattice using m = 0.1 and µ = 0.05. For x we use intervals of length ∆ n = 1 ∀n, such that the intervals are given by I n = [n, n + 1]. The symbols shown in Fig. 2 are the data for the intervals I n with n = 0, n = 10, n = 20, n = 50, n = 80 and n = 120. The lines are the fits of V (1) n (λ) with h(∆ n [λ−k (1) n ]) where h(s) is defined in (17).
From the fits of the restricted vacuum expectation value data with h(∆ n [λ−k (1) n ]) we can determine all slopes k (1) n , and from those compute the density ρ (1) (x) using the closed expressions (9). In Fig. 3 we show our results for ln ρ (1) (x) as a function of x, again using the DoS FFA data for V = 16 × 16, m = 0.1 and µ = 0.05. Note that this now is a quantity where for the free case we can compute analytical reference results. These are also shown in Fig. 3 and we find excellent agreement between the DoS FFA data and the analytic results, and stress at this point that we use a logarithmic scale on the vertical axis in Fig. 3.
We point out that further smoothening of the density with suitable fits will be part of a final strategy for DoS techniques -see, e.g., the recent systematic comparison of such techniques in [14].
We conclude this section with commenting on how the analytic reference results shown in Fig. 3 were obtained: Starting from the definition (43) of the density ρ (1) (x) we may use the integral representation of the Dirac delta and find

IV. DOS FFA FOR THE CANONICAL FORMULATION OF LATTICE QCD
In this section we present the second new DoS approach to finite density lattice QCD, now working with the canonical ensemble. The canonical partition sums at different net-quark numbers are expressed as Fourier moments of the grand canonical partition sum at imaginary chemical potential µ = iθ/β and then θ is considered as an additional degree of freedom in the path integral. In this form we may implement the DoS FFA and compute the density ρ(θ).

A. Canonical ensemble and density of states
The setting is as in the previous section, i.e., we study lattice QCD in d dimensions with N f degenerate flavors of quarks, and the grand canonical partition sum Z(µ) is defined in (18) - (21). The canonical partition sums Z N at a fixed net quark number N can be obtained as Fourier integrals over an imaginary chemical potential µ = i θ /β, where β is the inverse temperature in lattice units, i.e., β = N d , with N d being the number of lattice points in time direction (= d-direction), The corresponding free energy density at fixed N is de- Simple bulk observables can be obtained as derivatives of f N with respect to couplings of the theory. An example is the vacuum expectation value of the scalar fermion bilinear, The derivative generates the insertion of Tr D −1 [U, µ], i.e., the traced inverse Dirac operator (quark propagator) as an additional factor in the path integral. Note that also in the quark propagator the chemical potential µ appears and is set to the complex value µ = iθ /β, used for projecting to fixed net quark number N . General vacuum expectation values at fixed N have the form The expressions for the observables at fixed net quark number N can be rewritten with the help of densities ρ (51) J [U, µ] is an arbitrary functional of the gauge fields, which, if it contains the quark propagator, may also depend on the chemical potential µ. Note that again different choices of J [U, µ] result in different densities ρ (J ) (θ) and as before we use a superscript J to make clear which density we refer to.
With the densities ρ (J ) (θ) vacuum expectation values O N at fixed net quark number can be expressed as It is important to note that the densities ρ (J ) (θ) have symmetries that should be identified, because this allows one to reduce the range of θ that one needs to integrate over. Thus also the ρ (J ) (θ) need to be determined only in the reduced range of θ which lowers the numerical cost. In the previous section we have used charge conjugation symmetry to show that the density ρ (1) (x) for the imaginary part x ≡ X[U, Φ] is even in x. Also here it is straightforward to establish that ρ (1) (θ) is even. As before In a similar way as in (53) one can show that also the general densities ρ (J ) (θ) are either even or odd functions, depending on the symmetry of the insertion J [U, µ] (after decomposition into C-even and C-odd parts if necessary). Thus the integrals (52) for evaluating observables only run from 0 to π and exploring charge conjugation symmetry cuts the numerical cost in half.
We conclude this subsection with discussing another interesting symmetry property of the density, which not necessarily can be used to reduce the numerical cost, but reflects an important aspect of the underlying physics: If QCD is in a purely hadronic phase this is equivalent to ρ (1) (θ) being 2π/3 periodic. This property corresponds to the Roberge-Weiss symmetry and can be seen as follows: The statement that QCD is in a purely hadronic phase means that Z N = 0 for all net quark numbers N that are not multiples of 3. We first assume that ρ (1) (θ) is 2π/3-periodic. Then we find which shows that a 2π/3-periodic density ρ (1) (θ) implies that only Z N where N is a multiple of 3 are nonvanishing.
For the inverse statement we can use the completeness and orthogonality of the Fourier factors e iθN and sum over N the Z N in the form (52) with factors e iθN and find where in the second step we used that ρ (1) (θ) is even which in turn leads to Z N = Z −N . The relation (55) implies that if the Z N vanish for values of N which are not multiples of 3 the density ρ (1) (θ) is 2π/3-periodic. We remark, that the representation (55) of course holds in both, the hadronic and a possible non-hadronic phase, and in our small numerical test below we will use the form (55) to determine the canonical partition sums Z N from a fit of the density according to (55).

B. Implementation of DoS FFA
Having discussed the densities ρ (J ) (θ) and their symmetries we can now start the implementation of DoS FFA. For convenience we introduce the notation det D[U, θ ] ≡ det D[U, µ] | µ = iθ /β . For imaginary chemical potential γ 5 -hermiticity guarantees that det D[U, θ ] is real, such that the factor det D[U, θ ] N f is real and positive for even N f (or sufficiently large mass in case N f is odd), and we may write [24]. Using (56) we may write the canonical partition sum as It is interesting to note that the gauge fields U ν (x) and the phase variable θ enter the path integral in the same way, i.e., both are integrated over in the path integral and appear in the exponent of the Boltzmann factor. Thus one may view θ as one more degree of freedom in the path integral and compare (57) with the generic form (1) used in the general discussion of the DoS FFA in Section 2. The exponent in the integral is the action S[U, θ ] for all dofs. and we have already identified the real part of the action as S R [U, θ ]. The imaginary part, which only depends on θ , may be identified as X[θ ] = θ , and the parameter α in (1) is identified with the negative of the net quark number, i.e., α = −N .
Having found a form of the problem that matches the generic form discussed in Section 2, we may identify the restricted vacuum expectation values needed for the determination of the densities ρ (J ) (θ). They are given by where, as mentioned before, the fermion determinant may be represented using pseudo-fermions. The restricted vacuum expectation values (58) do not have a complex action problem and may be computed with standard Monte Carlo techniques. Note that the imaginary chemical potential θ is an additional degree of freedom that is restricted to the interval I n = [θ n , θ n+1 ] and needs to be updated as well.
After evaluating the restricted vacuum expectation values θ (J ) n (λ) they need to be brought into the normalized form V (θ) are computed using (8), (9), and finally observables in the canonical picture at fixed net quark number N are obtained from the densities via (52).

C. Tests of the canonical DoS FFA in the free case
Again we use 2-d free fermions at finite density for a first test also in the canonical formulation of the DoS FFA. For the free case the density (51) for the choice J = 1 reduces to the particularly simple expression (we where D[µ] denotes the Wilson Dirac operator (20), (21) in d = 2 with all link variables set to U µ (x) = 1. It is straightforward to evaluate this quantity using Fourier transformation and the reference data used in Fig. 5 below for verification were computed in this way. The restricted vacuum expectation values θ (1) n (λ) defined in (58) reduce to θ (1)  Although we could use the symmetry of the density ρ (1) (θ) and restrict the determination of ρ (1) (θ) to the interval θ ∈ [0, π], we here determine the density for the full range θ ∈ [−π, π]. The symmetry of ρ (1) (θ) should emerge and serves as a consistency check for the calculation. The interval [−π, π] was divided into 100 equal size intervals I n of length ∆ n = 2π/100 ∀n. The Monte Carlo simulation for sampling the restricted θ-integral in each interval I n uses a statistics of 10 6 sweeps of local Metropolis updates separated by 20 sweeps for decorelation and 10 5 sweeps for initial equilibration. The determinant in the acceptance step was computed with Fourier transforma-   (1) (θ) as a function of θ. We compare the DoS FFA result (thin blue curve) with the exact result (thick magenta curve). Note that we did not use the fact that the density is known to be an even function and for evaluation purposes numerically determined ρ (1) (θ) in the full range θ ∈ [−π, π]. tion und we typically use 10 values of λ for the evaluation of the restricted vacuum expectation values θ (1) n (λ). In Fig. 4 we show the results for the restricted vacuum expectation values θ (1) n (λ) already in their normalized form V (1) n (λ) according to (15). The symbols represent the data from the Monte Carlo simulation and the full curves are the fits with h(∆ n [λ − k (1) n ]) according to (17). The values of λ where the curves cross 0 are the slopes k (1) n . These crossing points start near 0 for the smallest n (i.e., intervals I n near −π) become negative then, revert back to 0, move to positive values and finally revert again back to 0 for intervals I n near +π. This full oscillation of the corresponding slopes k (1) n reflects the 2π-periodicity of the density ρ (1) (θ) (compare Fig. 5). From the slopes k (1) n obtained with the fits of the restricted vacuum expectation values we determined the density ρ (1) (θ) using (8) and (9). In Fig. 5 we compare the density determined in this way with the analytic result from Fourier transformation. The analytic result is represented by the thick magenta curve on top of which we plot the DoS result (thin blue curve). We stress again that the density ρ (1) (θ) was determined with DoS FFA for the full range θ ∈ [−π, π] and the fact that ρ (1) (θ) indeed comes out as an even function is a consistency check of the method. Anyway, the much more stringent test is the comparison with the analytic result where the plot shows that the DoS FFA curve perfectly falls on top of the exact curve determined as discussed above.
We complete our first test of the canonical DoS formulation with FFA by evaluating the canonical partition sums Z N from the density ρ (1) (θ) via the integrals (52) and comparing these Monte Carlo based results to the exact calculation based on a direct evaluation of (48) with Fourier transformation techniques.
In Fig. 6 we show the corresponding results for Z N normalized with Z 0 as a function of N . The blue diamonds represent the exact results and the red dots the DoS FFA data obtained with the integrals (52). The distribution resembles a Gaussian, rapidly decreasing with increasing |N | (which is of course a volume dependent statement). We find that the DoS FFA data based on (52) match the exact results very well.
We have already pointed out in the discussion of the direct DoS FFA approach that fitting the density with a suitable function will be an important part of future DoS strategies. Usually a large polynomial would be used for such a fit (see, e.g., [14,16,17] for related discussions), but for the canonical DoS approach the representation (55) of the density suggests another option for a fit, namely using a superposition of cosines (sines for odd densities). For the particular case of the density ρ (1) (θ) the fit parameters are the canonical partition sum Z N . In order to test this possibility, we determined the Z N also from a fit of ρ (1) (θ) with (55). The corresponding results are shown as black circles in Fig. 6 and again we find a very good agreement with the analytical results. This demonstrates that smoothening techniques based on periodic representations of the type (55) should be an interesting option to be explored in future development of canonical DoS techniques.
Also for the CanDos we would like to stress that the tests presented here constitute merely a very first assessment of the new approach and only the implementation in a full QCD simulation will show how well the numerical challenges can be brought under control in a calculation that includes the full gauge field dynamics.

V. SUMMARY, DISCUSSION AND OUTLOOK
In this article we have discussed two proposals for a modern DoS approach to finite density lattice QCD based on representations of the theory with pseudo-fermions. In the direct grand canonical approach the fermion determinant is represented with pseudo-fermions and subsequently their effective action is separated into real and imaginary parts such that the latter can then directly be treated with DoS FFA. We worked out the details of the formulation and provided some bounds on the involved kernels of the pseudo-fermion bilinears, showing that the method is applicable in an interesting range of values of the chemical potential µ. We presented very preliminary tests in the free case where a comparison to exact results allows one to assess the new approach. The direct DoS formulation in the grand canonical picture is rather straightforward, but has the disadvantage that also the densities depend on the chemical potential µ. As a consequence the densities have to be re-calcuated when changing µ. Whether this approach can beat our second suggestion, the canonical version of DoS FFA, has to be seen in future more detailed tests.
In the canonical formulation observables at a fixed net quark number N are obtained as the Fourier moments of the partition sum at imaginary chemical potential µ = iθ/β. In this setting we promote the angle θ to a new dynamical variable and interpret the exponent of the Fourier factors e −iθN as the imaginary part of the action. Again we treat this imaginary part with the DoS FFA approach and compute the density ρ(θ) as a function of θ. Observables at different net particle numbers N are then obtained by integrating the same density ρ(θ) with different Fourier factors e −iθN . Obviously here the re-sulting density ρ(θ) can be used for different net particle numbers N , but of course the accuracy of the determination of ρ(θ) has to be higher for larger N . Also here further tests that go beyond the first numerical checks we have presented here will be necessary to assess whether this formulation will be able to compete with other approaches to finite density QCD.
Both formulations we have suggested here, for the first time implement DoS techniques directly in a pseudofermion representation. This has the advantage that these well established techniques can be used in the framework of a modern DoS setting (here the DoS FFA is used but it is also straightforward to implement the ideas proposed here in the LLR framework). Obviously the simple exploratory numerical tests we have presented in this paper only serve to check the plausibility of the two new formulations and a much more detailed assessment will be necessary to explore their potential. Such further numerical tests are currently in preparation.
We conclude with remarking that the techniques developed here go beyond applications to finite density QCD.
The two approaches are general and can be applied to any lattice field theory with fermions where the interaction can be written with the help of a bosonic field such that the fermion action has a bilinear form and a fermion determinant emerges when integrating out the fermions. The bosonic fields do not have to be gauge fields, but also auxiliary fields of a Hubbard-Stratonovich transformation of quartic fermion interactions are a suitable option. We have begun to explore also these possible applications of the newly proposed DoS FFA techniques. Finally we remark that very recently [23] we presented a first test of the new approaches, now for the case of lattice QCD formulated with staggered fermions.