**1. Introduction**

One of the major open challenges for numerical lattice field theory is the treatment of QCD (Quantum Chromo Dynamics) at finite density. The central problem is the fact that at finite density, the fermion determinant is complex and cannot be used as a probability in Monte Carlo simulations. Density of states (DoS) techniques have been among the possible strategies for overcoming the complex action problem since the pioneering days of lattice QCD [1–6]. The key challenge for DoS techniques is accuracy, since for computing observables, the density needs to be integrated over with a highly oscillating factor. A simple sampling of the density with histogram techniques will allow one to access only very low densities.

An important step for the further development of DoS techniques was presented in [7] where, based on ideas from statistical mechanics [8], a suitable parameterization of the density combined with restricted vacuum expectation values was used to improve the accuracy for the determination of the density of states considerably. In a subsequent series of papers, this so-called LLR method was developed further and assessed for several test cases [9–16]. A related DoS technique, the so-called functional fit approach (FFA), was proposed in [17] and successfully tested in [18–21].

However, all these DoS techniques were formulated for bosonic systems, and no approach to finite density lattice QCD with modern DoS techniques had been presented. Finally, in [22], two possible formulations of DoS techniques for lattice field theories with fermions were suggested. One of the two formulations is the canonical DoS approach (CanDoS) where the density is computed as a function of the imaginary chemical potential *μ* ≡ *iθ*/*β*, where *β* is the inverse temperature. The canonical partition sum and observables are then obtained as Fourier moments of the density, and the FFA can be used to

#### *Particles* **2020**, *3*

obtain sufficient accuracy also for the highly oscillating integrals for the higher Fourier modes at large net particle numbers.

The second DoS approach presented in [22] is a direct grand canonical DoS formulation (GCDoS) based on rewriting the grand canonical partition sum of lattice QCD with a suitable pseudo-fermion representation and identifying the imaginary part of the action in this representation. Subsequently, FFA can be applied to evaluate the density as a function of the imaginary part, and again, suitable integrals over the density give rise to vacuum expectation values of observables.

In [22], the two new DoS approaches were presented for the formulation of lattice QCD with Wilson fermions, and the first tests were presented for free Wilson fermions at finite density. In this paper, we now discuss the CanDos formulation and the direct GCDoS approach for the formulation of lattice QCD with staggered fermions. For the CanDos approach, we also present some exploratory tests in the free case, which allows one to assess the accuracy of the method with exact results and to explore the parameters of the new techniques.

#### **2. The Canonical Density of States Approach**

In this section, we present the basic formulation of the canonical DoS approach (CanDos) for finite density lattice QCD. We stress, however, that the CanDoS approach can easily be implemented for other fermionic theories, e.g., theories with four Fermi interactions generated with auxiliary Hubbard–Stratonovich fields.

#### *2.1. Canonical Ensemble and Density of States*

We study lattice QCD in *d* dimensions with two degenerate flavors of quarks. The canonical partition sum at a fixed net quark number *N* is given by:

$$Z\_N = \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \int \mathcal{D}[\mathcal{U}] \, e^{-S\_G[\mathcal{U}]} \, \det D[\mathcal{U}, \mu] \,^2 \Big|\_{\mu = \text{i}\frac{\theta}{\mathcal{B}}} e^{-i\theta N} \,\tag{1}$$

where *SG*[*U*] is the Wilson gauge action (we dropped the constant additive term),

$$S\_G[\mathcal{U}] = -\frac{\beta\_G}{3} \sum\_{\mathbf{x}, \mathbf{v} < \rho} \text{Re } \text{Tr } \mathcal{U}\_\mathbf{v}(\mathbf{x}) \mathcal{U}\_\theta(\mathbf{x} + \boldsymbol{\vartheta}) \, \mathcal{U}\_\mathbf{v}(\mathbf{x} + \boldsymbol{\rho})^\dagger \, \mathcal{U}\_\theta(\mathbf{x})^\dagger. \tag{2}$$

*βG* is the inverse gauge coupling, and the path integral measure D[*U*] in (1) is the product of Haar measures for the link variables *<sup>U</sup>ν*(*x*) ∈ SU(3). We already integrated out the fermions and obtained the fermion determinants for the two flavors. *<sup>D</sup>*[*<sup>U</sup>*, *μ*] is the Dirac operator at finite chemical potential *μ*. In this study of the canonical DoS approach, we use the staggered Dirac operator, but stress that it is straightforward to implement the formalism also for different discretizations of the Dirac operator, e.g., for Wilson fermions (compare [22]). The staggered Dirac operator *<sup>D</sup>*[*<sup>U</sup>*, *μ*] is given by:

$$D[\mathcal{U}, \mu]\_{\mathbf{x}, \boldsymbol{y}} = \boldsymbol{m} \, \delta\_{\mathbf{x}, \boldsymbol{y}} \, \mathbb{1}\_3 + \frac{1}{2} \sum\_{\nu=1}^d \eta\_{\nu}(\mathbf{x}) \left[ e^{\mu \, \delta\_{\boldsymbol{\vartheta}, \boldsymbol{\vartheta}}} \mathcal{U}\_{\boldsymbol{\vartheta}}(\mathbf{x}) \, \delta\_{\mathbf{x} + \boldsymbol{\vartheta}, \boldsymbol{y}} - e^{-\mu \, \delta\_{\boldsymbol{\vartheta}, \boldsymbol{\vartheta}}} \mathcal{U}\_{\boldsymbol{\vartheta}}(\mathbf{x} - \boldsymbol{\vartheta})^\dagger \, \delta\_{\mathbf{x} - \boldsymbol{\vartheta}, \boldsymbol{y}} \right], \tag{3}$$

where *ην*(*x*)=(−<sup>1</sup>)*<sup>x</sup>*1<sup>+</sup> ... +*xν*−<sup>1</sup> are the staggered sign factors and 13 is the unit matrix in color space. We work on a *d*-dimensional lattice of size *Nd*−<sup>1</sup> *S* × *NT*, where the temporal (*ν* = *d*) extent *NT* gives the inverse temperature in lattices units, i.e., *β* = *NT*. All boundary conditions are periodic, except for the anti-periodic temporal (*ν* = *d*) boundary conditions for the fermions. *m* denotes the bare quark mass and *μ* the chemical potential.

In order to project the partition function *ZN* to fixed net quark number *N*, in (1), the chemical potential *μ* is set to *μ* = *iθ*/*β* = *iθ*/*NT* and subsequently integrated over the angle *θ* with a Fourier factor *e* −*iθN*. This Fourier transformation with respect to the imaginary chemical potential sets the

*Particles* **2020**, *3*

net quark number to *N* and thus generates *ZN*. The corresponding free energy density is defined as *fN* = − ln *ZN*/*V*, where *V* = *Nd*−<sup>1</sup> *S NT* denotes the *d*-dimensional volume.

Bulk observables and their moments can be obtained as derivatives of *fN* with respect to couplings of the theory. A simple example, which we also will consider in our numerical tests below, is the chiral condensate *ψ*(*x*)*ψ*(*x*)*N* = *∂ fN*/*∂<sup>m</sup>*,

$$\langle \sqrt{\boldsymbol{\Psi}}(\mathbf{x})\boldsymbol{\Psi}(\mathbf{x})\rangle\_{\mathcal{N}} = -\frac{2}{V} \frac{1}{Z\_N} \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \int \mathcal{D}[\boldsymbol{\mathcal{U}}] \, e^{-S\_G[\boldsymbol{\mathcal{U}}]} \, \det D[\boldsymbol{\mathcal{U}}, \boldsymbol{\mu}]^2 \, \text{Tr} \, D^{-1}[\boldsymbol{\mathcal{U}}, \boldsymbol{\mu}] \bigg|\_{\boldsymbol{\mu} = \boldsymbol{i}\frac{\theta}{\theta}} e^{-i\theta \boldsymbol{N}}.\tag{4}$$

The mass derivative leads to the insertion of Tr *<sup>D</sup>*−<sup>1</sup>[*<sup>U</sup>*, *μ*] in the path integral. Similarly, general vacuum expectation values of some observable O at fixed net quark number *N* have the form:

$$\langle \mathcal{O} \rangle\_N = \frac{1}{Z\_N} \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \int \mathcal{D}[\mathcal{U}] \, e^{-S\_C[\mathcal{U}]} \, \det \, D[\mathcal{U}, \mu]^2 \, \mathcal{O}[\mathcal{U}, \mu] \Big|\_{\mu = i\frac{\theta}{\tilde{\mathcal{P}}}} \, e^{-i\theta N} \,. \tag{5}$$

The partition sum (1) and the expressions for the vacuum expectation values (5) can be written with suitable densities *<sup>ρ</sup>*(*J*)(*θ*), which we define as:

$$\rho^{(l)}(\theta) = \int \mathcal{D}[\mathcal{U}] \, e^{-\mathcal{S}\_{\mathbb{C}}[\mathcal{U}]} \, \det \mathcal{D}[\mathcal{U}, \mu]^2 \left. f[\mathcal{U}, \mu] \right|\_{\mu = \frac{\theta}{\mathcal{B}}} \tag{6}$$

where *J*[*<sup>U</sup>*, *μ*] is set to *J*[*<sup>U</sup>*, *μ*] = 1 for the partition sum and to *J*[*<sup>U</sup>*, *μ*] = O[*<sup>U</sup>*, *μ*] for the vacuum expectation values of observables. With the densities *<sup>ρ</sup>*(*J*)(*θ*), we may express O*N* and *ZN* as:

$$\langle \mathcal{O} \rangle\_N = \frac{1}{Z\_N} \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \rho^{(\mathcal{O})}(\theta) \, e^{-i\theta N}, \qquad Z\_N = \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \rho^{(1)}(\theta) \, e^{-i\theta N}. \tag{7}$$

Note that charge conjugation symmetry can be used to show that *ρ*(1)(*θ*) is an even function such that *ρ*(1)(*θ*) needs to be determined only in the range *θ* ∈ [0, *<sup>π</sup>*], which cuts the numerical cost in half (see, e.g., [22]). A general observable O[*<sup>U</sup>*, *μ*] can be decomposed into even and odd parts under charge conjugation such that also here, the corresponding densities *ρ*(*J*)(*θ*) need to be evaluated only for *θ* ∈ [0, *<sup>π</sup>*].

Having defined the densities *ρ*(*J*)(*θ*) and expressed observables in the canonical ensemble as integrals over the densities, we now have to address the problem of finding a suitable representation of the density and how to determine the parameters used in the chosen representation.

#### *2.2. Parametrization of the Density*

We need to determine the densities *ρ*(*J*)(*θ*) for different operator insertions *J* as discussed in the previous section. For notational convenience, in this section, where we now discuss the parameterization of the densities, we denote all densities as *ρ*(*θ*), but stress that we need to determine the parameters of the different *ρ*(*θ*) independently for every choice of *J*.

The densities *ρ*(*θ*) are general functions of *θ* in the interval [0, *<sup>π</sup>*], which for a numerical determination, we need to describe with only a finite number of parameters. To obtain a suitable parameterization, we divide the interval [0, *π*] into *M* subintervals as,

$$[0, \pi] \quad = \bigcup\_{n=0}^{M-1} I\_{\mathbb{H}} \quad \text{with} \quad I\_{\mathbb{H}} = [\theta\_{\mathbb{H}}, \theta\_{n+1}]\_{\prime} \tag{8}$$

where *θ*0 = 0 and *θM* = *π*. Introducing Δ*n* = *θn*+<sup>1</sup> − *θn* for the length of the intervals *In*, we find *θn* = ∑*<sup>n</sup>*−<sup>1</sup> *j*=0 Δ*j* for *n* = 0, 1, ... *M*. For the densities *ρ*(*θ*), we now make the ansatz:

*ρ*(*θ*) = *e* −*<sup>L</sup>*(*θ*), (9)

where the *L*(*θ*) are continuous functions that are piecewise linear on the intervals *In*. We use the normalization *L*(0) = 0, which in turn implies *ρ*(0) = 1. Introducing a constant *an* and a slope *kn* for the linear function in every interval *In*, we may write *L*(*θ*) in the form:

$$L(\theta) = a\_n + k\_n \left[ \theta - \theta\_n \right], \text{ for } \theta \in I\_n = \left[ \theta\_{n\prime} \theta\_{n+1} \right]. \tag{10}$$

Since the functions *L*(*θ*) are normalized to *L*(0) = 0 and are required to be continuous, we can uniquely determine the constants *an* as functions of the slopes *kn* and write *L*(*θ*) in the following closed form:

$$L(\theta) = d\_n + \theta \, k\_{n\prime} \quad \theta \in I\_{n\prime} \quad d\_n = \sum\_{j=0}^{n-1} \left[ k\_j - k\_n \right] \Delta\_j \text{ for } n = 0, \dots \, M\_{\prime} \tag{11}$$

and express the densities *ρ*(*θ*) as:

$$\rho(\theta) = A\_n e^{-\theta k\_n}, \quad \theta \in I\_n, \quad A\_n = e^{-d\_n}. \tag{12}$$

Obviously, the parameterized density *ρ*(*θ*) depends only on the *kn*, i.e., the set of slopes of the linear pieces in the intervals *In*. We point out that our parametrization allows one to work with intervals *In* of different sizes Δ*n* such that in regions where the density *ρ*(*θ*) varies quickly, one may choose small Δ*<sup>n</sup>*, while in regions of slow variation, one may save computer time by working with larger Δ*<sup>n</sup>*.

#### *2.3. Evaluation of the Parameters of the Density*

To compute the slopes *kn* that determine the densities, we introduce so-called restricted expectation values *θ n*(*λ*) that are defined as:

$$\mathcal{I}(\boldsymbol{\theta})\_{\boldsymbol{n}}(\boldsymbol{\lambda}) \equiv \frac{1}{Z\_{\boldsymbol{n}}(\boldsymbol{\lambda})} \int\_{\boldsymbol{\theta}}^{\boldsymbol{\theta}\_{\boldsymbol{n}+1}} \int\_{\boldsymbol{\theta}} d\boldsymbol{\theta} \int \mathcal{D}[\boldsymbol{\mathcal{U}}] \, e^{-S\_{G}[\boldsymbol{\mathcal{U}}]} \, \boldsymbol{\theta} \, e^{\boldsymbol{\theta}\boldsymbol{\lambda}} \, \det \boldsymbol{D}[\boldsymbol{\mathcal{U}}, \boldsymbol{\mu}]^{2} \, \boldsymbol{J}[\boldsymbol{\mathcal{U}}, \boldsymbol{\mu}] \Big|\_{\boldsymbol{\mu} = \boldsymbol{i}\frac{\boldsymbol{\theta}}{\boldsymbol{\theta}}} \tag{13}$$

where again either *J*[*<sup>U</sup>*, *μ*] = 1 or *J*[*<sup>U</sup>*, *μ*] = O[*<sup>U</sup>*, *μ*] is chosen, depending on whether the slopes of the density for the partition sum *ZN* or the vacuum expectation O*N* are being computed. The corresponding restricted partition sums *Zn*(*λ*) we use in (13) are given by:

$$Z\_{\pi}(\lambda) \equiv \int\_{\theta\_{\pi}}^{\theta\_{\pi+1}} d\theta \int \mathcal{D}[\mathcal{U}] \ e^{-S\_{\mathbb{G}}[\mathcal{U}]} \ e^{\theta \lambda} \, \det D[\mathcal{U}, \mu]^2 \left. f[\mathcal{U}, \mu] \right|\_{\mu = i\frac{\theta}{\theta}}.\tag{14}$$

In the restricted expectation values *θ n*(*λ*) and the partition sum *Zn*(*λ*), the phase angle *θ* is integrated only over the interval *In*. We have also introduced a free real parameter *λ*, which couples to *θ* and enters in exponential form. Varying this parameter allows one to explore the *θ*-dependence of the density in the whole interval *In* fully. Since for imaginary chemical potential *μ* = *iθ*/*β*, the fermion determinant is real and after squaring also positive, the expectation values *θ n*(*λ*) can be evaluated without complex action problem in a Monte Carlo simulation as long as the insertions *J* are real and positive (for general insertions, *J* needs to be decomposed into pieces that obey positivity). This is

a technical issue that may be solved also in other ways, e.g., for a bounded observable, the addition of a positive constant is a simple option.

The important observation now is that for the parameterization (12) we have chosen for the densities, *θ n*(*λ*) and *Zn*(*λ*) can be computed also in closed form. Writing the partition sum with the density and then inserting the form (12), one obtains:

$$Z\_n(\lambda) = \int\_{\theta\_n}^{\theta\_{n+1}} d\theta \,\rho(\theta) \, e^{\theta \lambda} = e^{-d\_n} \int\_{\theta\_n}^{\theta\_{n+1}} d\theta \, e^{-\theta \, k\_n} e^{\theta \lambda} = e^{-d\_n} \frac{e^{\theta\_n[\lambda - k\_n]}}{\lambda - k\_n} \left( e^{\Lambda\_\mathbb{Z}[\lambda - k\_n]} - 1 \right). \tag{15}$$

From a comparison of (13) and (14), one finds that the restricted vacuum expectation value *θ n*(*λ*) can be computed as the derivative *θ n*(*λ*) = *d* ln *Zn*(*λ*)/*dλ*, such that also *θ n*(*λ*) can be found in closed form:

$$\langle \theta \rangle\_n(\lambda) \equiv \frac{d \ln Z\_n(\lambda)}{d \lambda} = \theta\_n + \frac{\Delta\_n}{1 - e^{-\Lambda\_n[\lambda - k\_n]}} - \frac{1}{\lambda - k\_n}.\tag{16}$$

Using a multiplicative and an additive normalization, we bring *θ n*(*λ*) into a standard form *Vn*(*λ*) where the result is expressed in terms of a simple function *h*(*s*),

$$V\_n(\lambda) \equiv \frac{\langle \theta \rangle\_n(\lambda) - \theta\_\mathbb{R}}{\Delta\_\mathbb{R}} - \frac{1}{2} = h\left(\Delta\_\mathbb{R}[\lambda - k\_\mathbb{R}]\right) \quad \text{with} \quad h(s) \equiv \frac{1}{1 - e^{-s}} - \frac{1}{s} - \frac{1}{2}.\tag{17}$$

The function *h*(*s*) obeys *h*(0) = 0, *h*(0) = 1/12, and lim*s*→±∞ *h*(*s*) = ±1/2.

The determination of the slope *kn* for the interval *In* now consists of the following steps: For several values of *λ*, one computes the corresponding restricted vacuum expectation values *θ n*(*λ*) defined in (14) and brings them into the normalized form *Vn*(*λ*) defined in Equation (17). Fitting the corresponding data with *<sup>h</sup>*<sup>Δ</sup>*n*[*<sup>λ</sup>* − *kn*] allows one to determine the *kn* from a simple stable one-parameter fit. From the sets of the slopes *kn*, we can determine the densities *ρ*(*θ*) using (11) and (12) and finally compute the observables via the integrals (7).

#### **3. An Exploratory Test of the Canonical DoS Approach in the Free Case**

As a first assessment of the new canonical density of states approach, we tested the new method for the case of free fermions at finite density in two dimensions. This served to verify the method and the program and allowed for exploring the parameters of the method, such as the number of intervals *In* and suitable choices for the values of *λ*. In addition, for the free case, all steps of the CanDoS approach could be cross-checked with exact results obtained from Fourier transformation.

#### *3.1. Setting and Reference Results from Fourier Transformation*

For this first test, we used the chiral condensate at fixed particle number *ψ*(*x*)*ψ*(*x*)*N* = *∂ fN*/*∂m* as our main observable. For the free case, the corresponding expression (4) reduces to:

$$<\langle \overline{\psi}(\mathbf{x})\psi(\mathbf{x})\rangle\_{N} = -\frac{2}{V} \frac{1}{Z\_{N}} \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \det D[\mu]^{2} \operatorname{Tr} D^{-1}[\mu] \Big|\_{\mu=i\frac{\theta}{\overline{\theta}}} e^{-i\theta N} \Big. \tag{18}$$

where all links in the Dirac operator (3) were set to *<sup>U</sup>ν*(*x*) = 1. For implementing the CanDoS approach for the condensate, we need the two densities,

$$\rho^{(1)}(\theta) = \det D[\mu]^2 \Big|\_{\mu = i\frac{\theta}{\tilde{\mathbb{B}}}} \quad \text{and} \quad \rho^{(\text{Tr}\,D^{-1})}(\theta) = \left. \det D[\mu]^2 \text{Tr} \, D^{-1}[\mu] \right|\_{\mu = i\frac{\theta}{\tilde{\mathbb{B}}}}.\tag{19}$$

For determining the slopes *kn* of these two densities, we thus have to compute the restricted expectation values (13) for *J* = 1 and *J* = Tr *D*−1. Normalizing the corresponding Monte Carlo data *Particles* **2020**, *3*

according to (17) and fitting them with *<sup>h</sup>*<sup>Δ</sup>*n*[*<sup>λ</sup>* − *kn*] gives rise to the slopes *kn*. From the respective sets of slopes, we find the densities *ρ*(1)(*θ*) and *ρ*(Tr *<sup>D</sup>*−<sup>1</sup>)(*θ*) using (11) and (12), and finally, the vacuum expectation value *ψ*(*x*)*ψ*(*x*)*N* is obtained as:

$$\langle \langle \overline{\Psi}(\mathbf{x}) \Psi(\mathbf{x}) \rangle\_{\mathcal{N}} = -\frac{2}{V} \frac{1}{Z\_N} \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \rho^{(\text{Tr}\, \mathcal{D}^{-1})}(\theta) \, e^{-i\theta \mathcal{N}}, \qquad Z\_{\mathcal{N}} = \int\_{-\pi}^{\pi} \frac{d\theta}{2\pi} \rho^{(1)}(\theta) \, e^{-i\theta \mathcal{N}}.\tag{20}$$

In the free case, the reference results can be obtained with the help of Fourier transformation. Furthermore, for the case of two flavors in two dimensions, which we are using for our test, we can explore the relation det *<sup>D</sup>*[*μ*]<sup>2</sup> = det *Dnaive*[*μ*] between the determinant of the staggered Dirac operator *<sup>D</sup>*[*μ*] and the determinant of the naive Dirac operator *Dnaive*[*μ*], which in two dimensions is given by:

$$D\_{\rm naive}[\mu]\_{\rm x,y} = m \, \delta\_{\rm x,y} \, \mathbb{1}\_2 \times \mathbb{1}\_3 \, + \, \frac{1}{2} \sum\_{\nu=1}^2 \sigma\_\nu \times \mathbb{1}\_3 \left[ e^{\mu \delta\_{\rm v,2}} \delta\_{\rm x+\partial,y} - e^{-\mu \delta\_{\rm v,2}} \delta\_{\rm x-\partial,y} \right],\tag{21}$$

where *σ*1 and *σ*2 are the first two Pauli matrices acting on the Dirac indices of the two-component spinors used in the naive discretization and 12 is the corresponding unit matrix. All link variables were set to their trivial values, i.e., they were replaced by the 3 × 3 unit matrix 13. The determinant of the naive Dirac operator can be computed by first diagonalizing *Dnaive*[*μ*] in space-time with the help of Fourier transformation and then taking the product of the corresponding momentum space Dirac operator determinants over all momenta.

The density *ρ*(1)(*θ*) then was simply obtained via numerically evaluating det *Dnaive*[*μ*] for *μ* = *iθ*/*β*. For the density *ρ*(Tr *<sup>D</sup>*−<sup>1</sup>)(*θ*), one may use Jakobi's formula (*d* det *M*/*dx* = det *M* Tr[*<sup>M</sup>*−<sup>1</sup> *dM*/*dx*]) for the derivative of a determinant and the fact that *dD*/*dm* = 1 to obtain:

$$\rho^{\left(\text{Tr}\,D^{-1}\right)}(\theta) = \det D[\mu]^2 \,\text{Tr}\,D^{-1}[\mu]\Big|\_{\mu = i\frac{\theta}{\beta}} = \frac{1}{2} \left. \frac{d}{dm} \det D[\mu]^2 \right|\_{\mu = i\frac{\theta}{\beta}} = \frac{1}{2} \left. \frac{d}{dm} \det D\_{\text{naive}}[\mu] \right|\_{\mu = i\frac{\theta}{\beta}}.\tag{22}$$

The vacuum expectation value *ψ*(*x*)*ψ*(*x*)*N* can be obtained from (20) by numerically integrating over *θ*. For the reference data in the plots below, we implemented this integration with Mathematica.

#### *3.2. Numerical Results for CanDos in the Free Case*

Having discussed the observables and the corresponding densities for the free case, as well as the evaluation of reference data with the help of Fourier transformation, we now come to a brief exploratory numerical test for the free case in *d* = 2 dimensions. The results in the plots below were computed on 16 × 16 lattices at a mass parameter of *m* = 0.1. We used 50 intervals *In* of equal size to parameterize the density in the range [0, *<sup>π</sup>*]. For each interval, we computed the restricted expectation values (16) for 20 different values of *λ* using Monte Carlo simulations based on 10<sup>6</sup> measurements, where in the simulation, the fermion determinant was evaluated exactly with Fourier transformation. The restricted expectation values were then normalized to the form (17) and the slopes *kn* determined from the corresponding fits with *<sup>h</sup>*<sup>Δ</sup>*n*[*<sup>λ</sup>* − *kn*]. From the slopes, the densities were computed using (11) and (12).

In Figure 1, we show the results for the densities *ρ*(1)(*θ*) (lhs plot) and *ρ*( Tr *<sup>D</sup>*−<sup>1</sup>)(*θ*) (rhs). The thin blue curves are the results from the CanDos determination and the thick magenta curves the reference data computed with Fourier transformation as discussed in the previous subsection. Obviously, the CanDos densities matched the reference data very well.

*Particles* **2020**, *3*

**Figure 1.** The densities *ρ*(1)(*θ*) (lhs) and *ρ*( Tr *<sup>D</sup>*−<sup>1</sup>)(*θ*) (rhs figure; denoted as *ρ*(cond)(*θ*) in the plot). We compare the data from the canonical DoS (CanDoS) determination (thin blue curves) to the analytic results obtained with Fourier transformation (thick dashed magenta curves). The data are for 16 × 16 lattices with *m* = 0.1 and densities are normalized to *ρ*(0) = 1.

Having determined the densities, we can compute the canonical partitions sums *ZN* and vacuum expectation values at fixed net fermion number using (7). In the lhs plot of Figure 2, we show our results for the canonical partition sums *ZN* normalized by *Z*0 as a function of *N*. The results from the CanDos determination are shown as red dots, the reference data from Fourier transformation as black diamonds. Here as well, we observed essentially perfect agreemen<sup>t</sup> for all values of the net fermion number *N* we considered. A more physical quantity is the corresponding free energy density *fN* = − ln *ZN*/*V* (here normalized to *f*0 = 0), which in the rhs plot of Figure 2, we show as a function of *N*. Again, we compared the CanDos results (red dots) to the corresponding reference data (black diamonds) and found very good agreement, and only for the largest net particle number *N* = 10 shown in the plot, we observed a slight deviation, indicating that for net quark numbers *N* > 10, the accuracy of the determination of the density would have to be improved, e.g., by using more and finer intervals *In*.

**Figure 2.** The canonical partition sums *ZN*/*Z*0 (lhs) and the corresponding free energy densities *fN* = − ln(*ZN*/*Z*0)/*<sup>V</sup>* (rhs) as a function of the net fermion number *N*. The parameters are *V* = 16 × 16 with *m* = 0.1, and we compare the results from the CanDoS determination (red dots) to the analytic results obtained with Fourier transformation (black diamonds).

We conclude our exploratory study with discussing the vacuum expectation value of an observable, i.e., a case where the ratio of two integrals over two different densities needs to be computed. The quantity we considered was the chiral condensate, and the two corresponding densities *ρ*(1)(*θ*) and *ρ*( Tr *<sup>D</sup>*−<sup>1</sup>)(*θ*) were the ones already shown in Figure 1. For both of them, we found very good agreemen<sup>t</sup> with the reference data, and the crucial question now was if this translated also into the corresponding physical observable matching the reference data well. In Figure 3, we show the CanDos results (red dots) for the condensate *ψ*(*x*)*ψ*(*x*)*N* as a function of the net quark number *N*. Indeed, we found a very satisfactory agreemen<sup>t</sup> with the results from Fourier transformation (black diamonds) up to *N* = 7 where the first deviations became visible. Again, for higher values of *N*, a more precise determination of the involved densities will be necessary.

**Figure 3.** The chiral condensate *ψ*(*x*)*ψ*(*x*)*N* (in the plot denoted as *condN* and normalized by the *N* = 0 value) as a function of the net fermion number *N*. The parameters are *V* = 16 × 16 with *m* = 0.1, and we compare the results from the CanDoS determination (red dots) to the analytic results obtained with Fourier transformation (black diamonds).

We close the discussion of our numerical test by stressing once more that the results presented here could only be considered a very preliminary assessment of the new CanDos approach. The tests were done in two dimensions, and only the free case was considered (although this already constituted a non-trivial test of the method). Currently, we are extending the assessment of CanDos by implementing a study in 2-dQCD, but also started to explore lattice field theories with four Fermi interactions.

#### **4. Direct Grand Canonical DoS Approach**

In this section, we now briefly discuss our second DoS approach, which is based on a suitable pseudo-fermion representation of the grand canonical QCD partition sum (GCDoS approach). We will determine the imaginary part of the pseudo-fermion action and set up the FFA to compute the density as a function of the imaginary part.

#### *4.1. Pseudo-Fermion Representation and Introduction of Densities*

The starting point was the grand canonical partition sum of QCD. We again considered two flavors of staggered fermions such that the grand canonical partition sum at chemical potential *μ* is given by:

$$Z\_{\mu} = \int \mathcal{D}[\mathcal{U}] \, e^{-S\_G[\mathcal{U}]} \, \det D[\mathcal{U}, \mu] \, ^2 \,, \tag{23}$$

where *SG*[*U*] is again the Wilson gauge action (2), and the staggered Dirac operator *<sup>D</sup>*[*<sup>U</sup>*, *μ*] is specified in (3).

*Particles* **2020**, *3*

We first identically rewrite the fermion determinant and subsequently express the part with the complex action problem in terms of pseudo-fermions,

$$\det D[\mathcal{U}, \mu] = \det(D[\mathcal{U}, \mu]D[\mathcal{U}, \mu]^\dagger) \frac{1}{\det D[\mathcal{U}, \mu]^\dagger} = \mathcal{C} \det(D[\mathcal{U}, \mu]D[\mathcal{U}, \mu]^\dagger) \int \mathcal{D}[\phi] e^{-\phi^\dagger D[\mathcal{U}, \mu]^\dagger \phi}, \tag{24}$$

where *C* is an irrelevant numerical constant and *φ*(*x*) are complex-valued pseudo-fermion fields that have three color components. The measure D[*φ*] simply is a product measure where at every site of the lattice, each component is integrated over the complex plane. The overall factor det(*D*[*<sup>U</sup>*, *μ*]*D*[*<sup>U</sup>*, *μ*] †) is obviously real and positive and can be treated with standard techniques [23,24]. The exponent of the pseudo-fermion integral on the other hand has a non-vanishing imaginary part and thus requires a strategy for dealing with the corresponding complex action problem.

To set up the direct DoS approach in the grand canonical formulation, we divided the exponent of the pseudo-fermion path integral into real and imaginary parts,

$$\boldsymbol{\phi}^{\dagger} \boldsymbol{D} [\boldsymbol{l} \boldsymbol{L} \boldsymbol{\mu}]^{\dagger} \boldsymbol{\phi} = \boldsymbol{S}\_{\mathcal{R}} [\boldsymbol{\phi}, \boldsymbol{l} \boldsymbol{L} \boldsymbol{\mu}] - \boldsymbol{i} \boldsymbol{X} [\boldsymbol{\phi}, \boldsymbol{l} \boldsymbol{L} \boldsymbol{\mu}],\\\boldsymbol{\mathcal{S}}\_{\mathcal{R}} [\boldsymbol{\phi}, \boldsymbol{l} \boldsymbol{L} \boldsymbol{\mu}] = \boldsymbol{\phi}^{\dagger} \boldsymbol{A} [\boldsymbol{l} \boldsymbol{l}, \boldsymbol{\mu}] \boldsymbol{\phi},\\\boldsymbol{X} [\boldsymbol{\phi}, \boldsymbol{l} \boldsymbol{L} \boldsymbol{\mu}] = \boldsymbol{\phi}^{\dagger} \boldsymbol{B} [\boldsymbol{l} \boldsymbol{l}, \boldsymbol{\mu}] \boldsymbol{\phi},\quad(25)$$

where we defined two matrices for the kernels of the real and imaginary parts of the pseudo-fermion action,

$$A[\![\![\!L\/\mu\]\!]]=\frac{D[\![\![\!L\/\mu\]\!]+D[\![\![\!L\/\mu\]\!]\!]}{2},\\B[\![\![\!L\/\mu\]\!]]=\frac{D[\![\![\![\!L\/\mu\]\!]\!]-D[\![\![\!L\/\mu\]\!]\!]}{2i}.\tag{26}$$

It is straightforward to evaluate *<sup>A</sup>*[*<sup>U</sup>*, *μ*] and *<sup>B</sup>*[*<sup>U</sup>*, *μ*] explicitly,

$$A[\![\mathcal{U},\mu]\_{\mathbf{x},\mathcal{Y}}\!] = \ \ m\delta\_{\mathbf{x},\mathcal{Y}}\! + \frac{1}{2}\sum\_{\nu=1}^{d}\eta\_{\nu}(\mathbf{x})\sinh(\mu\delta\_{\mathbf{v},\mathcal{U}})\left[\![\mathcal{U}\_{\nu}(\mathbf{x})\: \delta\_{\mathbf{x}+\boldsymbol{\uptheta},\mathcal{Y}}\!) + \![\mathcal{U}\_{\nu}^{\dagger}(\mathbf{x}-\boldsymbol{\uptheta})\: \delta\_{\mathbf{x}-\boldsymbol{\uptheta},\mathcal{Y}}\!)\right],$$

$$B[\![\mathcal{U},\mu]\_{\mathbf{x},\mathcal{Y}}\!] = -\frac{i}{2}\sum\_{\nu=1}^{d}\eta\_{\nu}(\mathbf{x})\cosh(\mu\delta\_{\boldsymbol{\upvartheta},\mathcal{U}})\left[\![\mathcal{U}\_{\boldsymbol{\upvartheta}}(\mathbf{x})\: \delta\_{\mathbf{x}+\boldsymbol{\uptheta},\mathcal{Y}}\!) - \![\mathcal{U}\_{\boldsymbol{\upvartheta}}(\mathbf{x}-\boldsymbol{\uptheta})\: \delta\_{\mathbf{x}-\boldsymbol{\uptheta},\mathcal{Y}}\!)\right].\tag{27}$$

The fermion determinant thus assumes the form:

$$\det D[\mathcal{U}, \mu] = \mathcal{C} \det(D[\mathcal{U}, \mu]D[\mathcal{U}, \mu]^\dagger) \int \mathcal{D}[\phi] \, e^{-S\_\mathcal{R}[\phi, \mathcal{U}] + iX[\phi, \mathcal{U}]}.\tag{28}$$

We already remarked that the real and positive overall factor det(*D*[*<sup>U</sup>*, *μ*]*D*[*<sup>U</sup>*, *μ*] †) could be treated with conventional simulation techniques [23,24], which we will not address in detail here (see [22] for a discussion of this term in the Wilson fermion formulation). Together with the Boltzmann factor for the gauge field action, we combined this term into a new effective action Boltzmann factor defined as:

$$\varepsilon e^{-S\_{eff}[\mathcal{U}, \mu]} = \varepsilon^{-S\_G[\mathcal{U}]} \det(D[\mathcal{U}, \mu]D[\mathcal{U}, \mu]^\dagger). \tag{29}$$

The grand-canonical partition sum thus can be written as:

$$Z\_{\mu} = \int \mathcal{D}[\mathcal{U}] \int \mathcal{D}[\phi] \,\, e^{-S\_{eff}[\mathcal{U},\mu]} \,\, e^{-S\_R[\phi,\mathcal{U},\mu]} \,\, e^{i\,\mathcal{X}[\phi,\mathcal{U},\mu]} \,. \tag{30}$$

The next step is to introduce suitable densities for the imaginary part:

$$\boldsymbol{\rho}^{(l)}(\mathbf{x}) = \int \mathcal{D}[\boldsymbol{\mathcal{U}}] \int \mathcal{D}[\boldsymbol{\phi}] \, e^{-S\_{eff}[\boldsymbol{\mathcal{U}}, \boldsymbol{\mu}]} \, e^{-S\_{R}[\boldsymbol{\phi}, \boldsymbol{\mathcal{U}}, \boldsymbol{\mu}]} \, f[\boldsymbol{\phi}, \boldsymbol{\mathcal{U}}, \boldsymbol{\mu}] \, \delta\left(\mathbf{x} - \boldsymbol{X}[\boldsymbol{\phi}, \boldsymbol{\mathcal{U}}, \boldsymbol{\mu}]\right), \tag{31}$$

where we again allow for the insertion of functionals *J*[*φ*, *U*, *μ*] in order to take into account different observables. As for the CanDos approach, one may use charge conjugation symmetry to show that the

densities are either even or odd functions of *x*, depending on the insertion *J*[*φ*, *U*, *μ*] (see [22]). Thus, it is sufficient to compute the densities only for positive *x*.

With the help of the densities vacuum, the expectation values of observables in the grand canonical picture at chemical potential *μ* can be written as:

$$\langle \mathcal{O} \rangle\_{\mu} = \frac{1}{Z\_{\mu}} \int\_{0}^{\infty} d\mathbf{x} \,\rho^{(\mathcal{O})}(\mathbf{x}) \, e^{i\mathbf{x}} \,, \, Z\_{\mu} = \int\_{0}^{\infty} d\mathbf{x} \,\rho^{(\mathcal{I})}(\mathbf{x}) \, e^{i\mathbf{x}} \,. \tag{32}$$

#### *4.2. Evaluation of the Densities with FFA*

Having defined the densities and expressed grand canonical vacuum expectation values as suitable integrals over the densities, we now can set up the FFA approach for evaluating the densities.

First, we remark that the densities *ρ*(*J*)(*x*) are expected to be fast decreasing functions of *x*, and in [22], this was indeed verified in test cases. Thus, we may cut off the integration range in (32) to a finite interval [0, *xmax*] and determine the density only for this range. As for the canonical case, we divided the interval [0, *xmax*] into *M* intervals *In* = [*xn*, *xn*+<sup>1</sup>], *n* = 0, 1, ... *M* − 1, with *x*0 = 0 and *xM* = *xmax*. As for the CanDos formulation, the densities were parameterized by the negative exponential of a function *<sup>L</sup>*(*x*) that was continuous and piecewise linear on the intervals *In*. Again, we assumed the normalization *L*(0) = 0, and the density thus was entirely determined by the slopes *kn*.

In the FFA approach, the slope *kn* in each interval *In* is determined from suitable restricted vacuum values, which we here define as:

$$\mathcal{I}\_{\boldsymbol{\Lambda}}(\boldsymbol{X})\_{\boldsymbol{n}}(\boldsymbol{\lambda}) = \frac{1}{Z\_{\boldsymbol{n}}(\boldsymbol{\lambda})} \int \mathcal{D}[\boldsymbol{l}\boldsymbol{l}] \int \mathcal{D}[\boldsymbol{\phi}] e^{-S\_{\boldsymbol{n}f\boldsymbol{\phi}}[\boldsymbol{l}\boldsymbol{l},\boldsymbol{\mu}]} e^{-S\_{\boldsymbol{n}}[\boldsymbol{\phi},\boldsymbol{l}\boldsymbol{l},\boldsymbol{\mu}]} e^{\boldsymbol{\lambda}\cdot\boldsymbol{X}[\boldsymbol{\phi},\boldsymbol{l}\boldsymbol{l},\boldsymbol{\mu}]} I[\boldsymbol{\phi},\boldsymbol{l}\boldsymbol{l},\boldsymbol{\mu}] \,\boldsymbol{\Theta}\_{\boldsymbol{n}}\left(\boldsymbol{X}[\boldsymbol{\phi},\boldsymbol{l}\boldsymbol{l},\boldsymbol{\mu}]\right), \tag{33}$$

where we have defined the support function <sup>Θ</sup>*n*(*x*):

$$\Theta\_n(\mathbf{x}) \, \, = \begin{cases} \quad \mathbf{1} \text{ for } \mathbf{x} \in I\_{n\prime} \\ \quad \, \, 0 \text{ else.} \end{cases} \tag{34}$$

As in the canonical case, also the generalized expectation values (33) can be expressed in terms of the parameterized density and computed in closed form, along the lines discussed above. After normalizing them to the form (17), the generalized expectation values are again described by the functions *<sup>h</sup>*<sup>Δ</sup>*n*[*<sup>λ</sup>* − *kn*] such that the slopes *kn* can be determined from one parameter fits. Subsequently, the densities are constructed from the slopes using (11) and (12), with *θ* replaced by *x*. Finally we can compute observables from the densities using (32).

The direct, grand canonical density of states approach discussed in this section for staggered fermions was discussed for Wilson fermions in [22]. There, also first exploratory numerical results were presented, and for free fermions it was shown that the density obtained with the FFA approach matched exact reference data from Fourier transformation very well.

#### **5. Summary and Outlook**

In this paper, we extended our previous work [22], where we presented two new DoS techniques for finite density lattice QCD with Wilson fermions, to the formulation of QCD with staggered fermions. The first formulation was based on the canonical formulation where the canonical partition sum and vacuum expectation values of observables at fixed net quark number were obtained as Fourier moments with respect to imaginary chemical potential. The functional fit approach (FFA) could then be used to compute the density with sufficient accuracy for reliably determining observables for reasonable net quark numbers. We presented exploratory tests of the canonical DoS approach for the case of free fermions in 2-dand found that observables such as the chiral condensate at finite net quark numbers reliably matched reference data obtained from a direct calculation with Fourier transformation that was possible in the free case.

Our second approach was set up directly in the grand canonical ensemble. The QCD partition sum was rewritten in terms of a suitable pseudo-fermion representation, and the imaginary part of the pseudo-fermion action was identified. Using FFA, the density was then computed as a function of the imaginary part, and grand canonical vacuum expectation values were again obtained as the corresponding oscillating integrals. The tests of the new approaches presented here were done for the staggered fermion formulation, but we would like to point out again that also the Wilson formulation could be used and refer to our paper [22] for the discussion of the corresponding results.

Two comments are in order here: Although the first tests were encouraging, the numerical results presented here clearly constituted only a very preliminary and exploratory assessment of the new techniques. We are currently extending these tests towards QCD in two dimensions as the next test case before approaching the full 4-dtheory. We furthermore stress that the techniques we presented here were not restricted to QCD or other gauge theories with fermions. Furthermore, theories with four Fermi interactions could be accessed after the introduction of suitable Hubbard–Stratonovich fields, and also for this direction of possible further development we have started exploratory calculations.

**Author Contributions:** Conceptualization and analytical work: C.G., M.M., and P.T.; software and numerical work: M.M. and P.T.; validation: C.G., M.M., and P.T.; paper writing: C.G. All authors read and agreed to the published version of the manuscript.

**Funding:** This work is supported by the Austrian Science Fund FWF, Grant I 2886-N27, and partly also by the FWF DK 1203 "Hadrons in Vacuum Nuclei and stars".

**Acknowledgments:** We thank Mario Giuliani, Peter Kratochwill, Kurt Langfeld, and Biagio Lucini for interesting discussions.

**Conflicts of Interest:** The authors declare no conflict of interest.
