Calculating the EoS of the dense quark-gluon plasma using the Complex Langevin equation

TL;DR

This paper computes the equation of state of dense quark-gluon plasma at finite chemical potential using an advanced Complex Langevin method with improved actions, showing good agreement with traditional approaches at small chemical potentials.

Contribution

It introduces a generalized stout smearing procedure for the Complex Langevin equation on the SL(3,C) manifold, enabling more accurate calculations of the quark-gluon plasma properties.

Findings

Good agreement with Taylor expansion at small chemical potentials.

Successful implementation of improved actions in the Complex Langevin framework.

First application of generalized stout smearing for this purpose.

Abstract

The pressure and energy density of the quark-gluon plasma at finite baryon chemical potential are calculated using the Complex Langevin equation. The stout smearing procedure is generalized for the SL(3,$\mathcal{C}$) manifold allowing the usage of an improved action in the Complex Langevin setup. Four degenerate flavors of staggered quarks with $m_\pi=500-700$ MeV are used with a tree-level Symanzik improved gauge action on $16^3 \times 8 $ lattices. Results are compared to the Taylor expansion and good agreement is found for small chemical potentials.

Tables

Table 1. Table 1: The lattice spacing and the pion mass using the plaquette action with naive staggered fermions with N F = 4 , m = 0.01 , formulae-sequence subscript 𝑁 𝐹 4 𝑚 0.01 N_{F}=4,\ m=0.01, measured on a 24 3 × 48 superscript 24 3 48 24^{3}\times 48 lattice.

$\beta$	$a$(fm)	$m_{\pi }a$	$T$(MeV) for $N_T=8$	pion mass(MeV)
5	$0.2892\pm 0.0002$	$0.2595\pm 0.0002$	85.3	177
5.1	$0.1895\pm 0.0005$	$0.2881\pm 0.0002$	130	300
5.2	$0.1105\pm 0.0004$	$0.2965\pm 0.0004$	223	529
5.3	$0.0822\pm 0.0003$	$0.2727\pm 0.0005$	300	654
5.4	$0.0633\pm 0.0005$	$0.2496\pm 0.0016$	389	777
5.5	$0.0503\pm 0.0005$	$0.2253\pm 0.0016$	490	883
5.6	$0.0433\pm 0.0006$	$0.2229\pm 0.0015$	570	1020

Table 2. Table 2: The lattice spacing and the pion mass using the Symanzik gauge action with stout smeared staggered fermions with N F = 4 , m = 0.02 , n = 2 , ρ = 0.15 , formulae-sequence subscript 𝑁 𝐹 4 formulae-sequence 𝑚 0.02 formulae-sequence 𝑛 2 𝜌 0.15 N_{F}=4,\ m=0.02,\ n=2,\ \rho=0.15, measured on a 24 3 × 48 superscript 24 3 48 24^{3}\times 48 lattice.

$\beta$	$a$(fm)	$m_{\pi }a$	$T$(MeV) for $N_T=8$	pion mass(MeV)
3.5	$0.1474\pm 0.0004$	$0.3111\pm 0.0004$	167	417
3.6	$0.1159\pm 0.0003$	$0.2790\pm 0.0004$	213	475
3.7	$0.0946\pm 0.0005$	$0.2515\pm 0.0005$	261	525
3.8	$0.0769\pm 0.0004$	$0.2259\pm 0.0009$	321	579
3.9	$0.0644\pm 0.0004$	$0.2088\pm 0.0016$	383	640
4	$0.0535\pm 0.0004$	$0.1987\pm 0.0024$	461	733
4.1	$0.0415\pm 0.0006$	$0.2119\pm 0.0052$	594	1010

Table 3. Table 3: The coefficients of the Taylor expansion of the pressure calculated at μ = 0 𝜇 0 \mu=0 using eq. ( 8 ) (label “Taylor exp.”) and by fitting a polynomial to the results of the integration method (label “CLE”). The unimproved action with N F = 4 , m = 0.01 formulae-sequence subscript 𝑁 𝐹 4 𝑚 0.01 N_{F}=4,\ m=0.01 is used on a 16 3 × 8 superscript 16 3 8 16^{3}\times 8 lattice.

$\beta$	$c_2$ Taylor exp.	$c_4$ Taylor exp.	$c_2$ CLE	$c_4$ CLE
5.2	$2.102\pm 0.059$	$0.233\pm 0.17$	$2.21\pm 0.2$	$0.15\pm 0.05$
5.3	$2.277\pm 0.026$	$0.095\pm 0.06$	$2.24\pm 0.1$	$0.18\pm 0.05$
5.4	$2.333\pm 0.016$	$0.146\pm 0.095$	$2.39\pm 0.1$	$0.16\pm 0.02$
5.5	$2.376\pm 0.019$	$0.125\pm 0.019$	$2.35\pm 0.1$	$0.18\pm 0.02$

Table 4. Table 4: The coefficients of the Taylor expansion of the pressure calculated at μ = 0 𝜇 0 \mu=0 using eq. ( 8 ) (label “Taylor exp.”) and by fitting a polynomial to the results of the integration method (label “CLE”). The 2-stout improved action is used with N F = 4 , m = 0.02 , formulae-sequence subscript 𝑁 𝐹 4 𝑚 0.02 N_{F}=4,m=0.02, on a 16 3 × 8 superscript 16 3 8 16^{3}\times 8 lattice.

$\beta$	$c_2$ Taylor exp.	$c_4$ Taylor exp.	$c_6$ Taylor exp.	$c_2$ CLE	$c_4$ CLE	$c_6$ CLE
3.7	$2.206\pm 0.009$	$0.156\pm 0.016$	$0.016\pm 0.013$	$2.33\pm 0.1$	$0.13\pm 0.02$	$0.002\pm 0.001$
3.8	$2.293\pm 0.007$	$0.171\pm 0.017$	$-0.01\pm 0.01$	$2.32\pm 0.1$	$0.14\pm 0.02$	$0.002\pm 0.002$
3.9	$2.312\pm 0.007$	$0.150\pm 0.007$	$0.001\pm 0.005$	$2.36\pm 0.04$	$0.14\pm 0.01$	$0.002\pm 0.001$
4.0	$2.371\pm 0.012$	$0.124\pm 0.009$	$-0.001\pm 0.006$	$2.43\pm 0.02$	$0.13\pm 0.01$	$0.002\pm 0.001$

Equations

$$U_{x,\nu }(\tau +ϵ)=\text{exp}\left[i\sum _a{\lambda }_a(ϵK_{ax\nu }+\sqrt{ϵ}{\eta }_{ax\nu })\right]U_{x,\nu }(\tau ),$$

$$D_{ax\nu }f(U)={{\partial }_{\alpha }f(e^{i\alpha {\lambda }_a}{U}_{x,\nu })|}_{\alpha =0}.$$

$$S_{eff}=S_g-\text{ln det}M(\mu ),$$

$$\frac{1}{4\Omega }\sum _{x,\mu }\text{Tr}((U_{x,\mu }{U}_{x,\mu }^{+}-1{)}^2),$$

$$\frac{p}{T^4}=\frac{\ln \mathcal{Z}}{V{T}^3}$$

$$\Delta \left(\frac{p}{T^4}\right)={\frac{p}{T^4}|}_{T,\mu }-{\frac{p}{T^4}|}_{T,\mu =0}.$$

$$\Delta \left(\frac{p}{T^4}\right)=\sum _{n>0,even}{c}_n(T){\left(\frac{\mu }{T}\right)}^n,$$

$$c_n(T)=\frac{1}{n!}\frac{N_t^3}{N_s^3}{\frac{{\partial }^n\ln \mathcal{Z}}{\partial (\mu N_t{)}^n}|}_{\mu =0}.$$

$$\frac{{\partial }^2\ln \mathcal{Z}}{\partial {\mu }^2}=⟨\frac{N_F}{4}\frac{{\partial }^2\ln \det M}{\partial {\mu }^2}⟩+⟨{\left(\frac{N_F}{4}\frac{\partial \ln \det M}{\partial \mu }\right)}^2⟩,$$

$$\frac{\partial \ln \det M}{\partial \mu }$$

$$\frac{{\partial }^2\ln \det M}{\partial {\mu }^2}$$

$$\Delta \left(\frac{p}{T^4}\right)$$

$$\displaystyle n={1\over\Omega}{\partial\ln\mathcal{Z}\over\partial\mu}={N_{f}\over 4\Omega}\big{\langle}\textrm{Tr}(M^{-1}\partial_{\mu}M)\big{\rangle}$$

$$\frac{ϵ-3p}{T^4}=-\frac{1}{V{T}^3}a{\left(\frac{\partial \beta }{\partial a}\right)}_{\text{LCP}}\left[\frac{\partial \text{ln}\mathcal{Z}}{\partial \beta }+{\left(\frac{\partial m}{\partial \beta }\right)}_{\text{LCP}}\frac{\partial \text{ln}\mathcal{Z}}{\partial m}\right],$$

$$\frac{{\partial }^3\text{ln}\mathcal{Z}}{\partial \beta \partial {\mu }^2},\frac{{\partial }^3\text{ln}\mathcal{Z}}{\partial m\partial {\mu }^2},$$

$$\frac{\partial ⟨O⟩}{\partial \beta }=-⟨O\frac{\partial S_g}{\partial \beta }⟩+⟨O⟩⟨\frac{\partial S_g}{\partial \beta }⟩,$$

$$\frac{\partial ⟨O⟩}{\partial m}=⟨\frac{\partial O}{\partial m}⟩+\Omega ⟨O\chi ⟩-\Omega ⟨O⟩⟨\chi ⟩$$

$$\Delta \left(\frac{ϵ-3p}{T^4}\right)={\frac{ϵ-3p}{T^4}|}_{\mu }-{\frac{ϵ-3p}{T^4}|}_{\mu =0}$$

$$C_{\nu }(x)=\sum _{\sigma \ne \nu }{\rho }_{\nu \sigma }\left(U_{\sigma }(x)U_{\nu }(x+\hat{\sigma })U_{\sigma }^{-1}(x+\hat{\nu })+U_{\sigma }^{-1}(x-\hat{\sigma })U_{\nu }(x-\hat{\sigma })U_{\sigma }(x-\hat{\sigma }+\hat{\nu })\right),$$

$$Z_{\nu }(x)=\sum _{\sigma \ne \nu }{\rho }_{\nu \sigma }\left(U_{\sigma }(x+\hat{\nu })U_{\nu }^{-1}(x+\hat{\sigma })U_{\sigma }^{-1}(x)+U_{\sigma }^{-1}(x-\hat{\sigma }+\hat{\nu })U_{\nu }^{-1}(x-\hat{\sigma })U_{\sigma }(x-\hat{\sigma })\right).$$

$${\Omega }_{\nu }(x)$$

$${\Omega }_{\nu }^i(x)$$

$$X_{\nu }(x)=i{Q}_{\nu }(x)$$

$$U\to U^{(1)}\to ...\to U^{(n)}$$

$$D_{ab\nu \sigma }=-\frac{i}{2}\text{Tr}\left({\lambda }_b(D_{a\nu }{U}_{\sigma }^{\prime })U_{\sigma }^{\prime -1}\right),$$

$$U_{\sigma }^{\prime }(e^{i{\alpha }_a{\lambda }_a}{U}_{\nu })=e^{i{\alpha }_a{D}_{ab\nu \sigma }{\lambda }_b}{U}_{\sigma }^{\prime },$$

$$F_{\nu }=i{\lambda }_a{D}_{a\nu }S[U]=-\frac{i{\lambda }_a}{2}\text{Tr}\left(F_{\sigma }^{\prime }(D_{a\nu }{U}_{\sigma }^{\prime })U_{\sigma }^{\prime -1}\right)=-\frac{i{\lambda }_a}{2}\text{Tr}\left(e^{-X_{\sigma }}{F}_{\sigma }^{\prime }{D}_{a\nu }{e}^{X_{\sigma }}\right)+e^{-X_{\nu }}{F}_{\nu }^{\prime }{e}^{X_{\nu }}-\frac{1}{N}\text{Tr}\left(e^{-X_{\nu }}{F}_{\nu }^{\prime }{e}^{X_{\nu }}\right)$$

$$F_{\nu }(x)$$

$$\exp (tX)=V\text{diag}(\exp (t{\lambda }_0),....,\exp (t{\lambda }_n))V^{-1}$$

$$\frac{\partial e^{tX}}{\partial \theta }=V(G\times E)V^{-1},$$

$$E_{ij}$$

$$E_{ii}$$

$$f(X)=f_0+f_{1}X+f_2{X}^2,$$

Full text

I Introduction and overview

The strong interactions are described by Quantum Chromodynamics (QCD). The determination of the QCD phase diagram is one of the great challenges of the theoretical study of this theory. It is phenomenologically relevant in many areas, such as in the early Universe, in relativistic heavy-ion collision experiments as well as in astronomy describing neutron stars. The lattice discretisation of the theory allows for the precise calculation of the Equation of State (EoS), the exploration of the hadronic and the quark-gluon plasma phase at zero baryonic density and the phase transition between them Petreczky (2012); Philipsen (2013); Borsanyi (2017). The EoS of the QCD matter can be calculated on the lattice using various methods. Usually the partition function is constructed by measuring its derivatives and integrating from a starting point at vacuum Engels et al. (1990); Boyd et al. (1996), but other approaches are also available Meyer (2009); Giusti and Meyer (2011); Suzuki (2013); Caselle et al. (2016). However, the lattice formulation of QCD suffers from a problem at nonzero chemical potential: the partition sum of the theory is written in terms of a complex measure due to the fermion determinant, thus the standard importance sampling approaches are invalid. This is called the QCD sign problem. Several different strategies have been proposed to circumvent the sign problem (see reviews Fodor and Katz (2009); de Forcrand (2009); Aarts (2012); Sexty (2014a)). The EoS at $\mu>0$ is traditionally calculated via Taylor expansion or reweighting Fodor et al. (2003); Allton et al. (2002, 2003, 2005); Borsanyi et al. (2012a); Bazavov et al. (2017) and analytical continuation D’Elia and Sanfilippo (2009); Guenther et al. (2017). In this study the Complex Langevin equation (CLE) is used to carry out simulations directly at $\mu>0$.

The complexification of the Langevin equation was proposed long ago Parisi (1983); Klauder (1983). The idea is to circumvent the sign problem by complexifing the field manifold of the theory and defining a stochastic process on this manifold using analiticity. After an initial excitement it was noticed that the Complex Langevin equation sometimes gives wrong results Ambjorn and Yang (1985); Ambjorn et al. (1986) and also practical problems (runaway trajectories) appeared. Recently, the method has enjoyed renewed attention and many of its problems have been solved. It has been proved that provided the action is holomorphic and the distributions of the variables decay fast enough, the complex Langevin equation will converge to the correct results Aarts et al. (2010a, 2011); Seiler (2018). In a recent study it has been proposed that using a certain observable the magnitude of the ’boundary term’ at infinity can be estimated Scherzer et al. (2019), with an observable that is cheap to calculate also for lattice systems Scherzer et al. . Gauge theories pose an additional problem: the complexification of the gauge degrees of freedom, which in turn leads to large fluctuations and the breakdown of the simulation. The procedure of gauge cooling Seiler et al. (2013); Aarts et al. (2013) was introduced to mitigate this problem. With the help of gauge cooling it became possible to simulate HDQCD (heavy dense QCD) where the quarks are kept static Seiler et al. (2013); Aarts et al. (2013), as well as full QCD using light quarks in the staggered Sexty (2014b); Schmalzbauer and Bloch (2016); Bloch and Schenk (2018); Nagata et al. (2018); Tsutsui et al. (2018); Kogut and Sinclair (2019) and the Wilson formulation Aarts et al. (2014a), as well as QCD with a theta term Bongiovanni et al. (2014).

Reaching is the continuum limit in lattice calculations is a non-trivial task. The usual strategy is to use improved actions which are more expensive numerically but they ensure a quicker convergence. Using an improved gauge action with CLE is straightforward, but fermionic improvements can be more involved. In this paper the stout smearing procedureMorningstar and Peardon (2004) is generalized to ensure applicability in the Complex Langevin setup.

In Section II a brief overview of the complex Langevin method is given. In Section III the strategy of the calculation of the pressure and energy density is discussed. Afterwards, the stout smearing procedure is generalized to the complexified manifold of the link variables in Section IV. In Section V the numerical results are presented. Finally, conclusions are offered in Section VI.

II The complex Langevin equation

The Langevin equation for $U_{x,\nu}\in$ SU($N$), the link variables of a gauge theory, in discretised form with Langevin timestep $\epsilon$ is written as Batrouni et al. (1985):

[TABLE]

with $\lambda_{a}$ the generators of the gauge group, i.e. the Gell-Mann matrices, the drift force $K_{ax\nu}=-(D_{ax\nu}\rho[U])/\rho[U]$ calculated from the measure $\rho$ using the left derivative

[TABLE]

For a complex measure the drift terms become complex with $K_{ax\nu}\in\mathcal{C}$. The manifold of the link variables is then complexified to SL(3,$\mathcal{C}$). In the case of lattice QCD with fermions the action is written as

[TABLE]

with the determinant of the fermionic Dirac matrix $M(\mu)$. The measure $e^{-S_{eff}}=e^{-S_{g}}\textrm{det}M(\mu)$ generally has zeroes on the complexified field manifold, resulting in meromorphic drift terms. Simulating such a theory with the CLE can potentially lead to incorrect results. It has been shown that in the case of QCD at large temperatures these zeroes are not reached by the process, thus the formal justification of the Complex Langevin method goes through Aarts et al. (2017).

The non-unitarity of the link variables can be monitored using the unitarity norm

[TABLE]

where $\Omega=N^{3}_{s}N_{t}$ is the space-time volume of the lattice. The uncontrolled growth of the unitarity norm observed in naive complex Langevin simulations can be countered using complexified gauge transformations after each update such that the unitarity norm is decreased, i.e. gauge cooling Seiler et al. (2013); Aarts et al. (2013) (see also Nagata et al. (2015) for the inclusion of gauge cooling into the formal proof of correctness). It has been observed that gauge cooling is effective as long as the $\beta$ parameter of the theory is not too small Aarts et al. (2014b). The minimal $\beta$ corresponds to a maximal lattice spacing, such that the continuum limit can be carried out in the safe region, allowing the mapping of the phase diagram of the HDQCD theory Aarts et al. (2016).

In this study the naive plaquette and staggered action is used as well as the tree-level Symanzik improved gauge action with stout smeared staggered fermionsMorningstar and Peardon (2004). The implementation of the gauge actions is straightforward: one ensures the holomorphicity of the action by replacing the matrix adjungate with the matrix inverse for the plaquette and extended plaquette variables appearing in the action. The naive fermionic drift is calculated with the help of noise vectors Sexty (2014b), and the implementation of the stout smearing is detailed in Sec. IV.

III Thermodynamics at nonzero chemical potential

Using the grand canonical ensemble the pressure in units of $T^{4}$ is calculated from the grand partition function $\mathcal{Z}(T,\mu)$:

[TABLE]

For the purposes of this study we assume that the pressure calculation at zero chemical potential has been carried out by some method. Our primary interest here is the change of the pressure as the chemical potential is increased at a fixed temperature, since a direct calculation at $\mu>0$ is not possible with the usual importance sampling simulations

[TABLE]

Usually one Taylor-expands this difference at $\mu=0$ to allow calculations of the coefficients using Monte-Carlo simulations Allton et al. (2003, 2005)

[TABLE]

with

[TABLE]

From the symmetry of the partition function $Z(\mu)=Z(-\mu)$ we see that only the even coefficients $c_{2k}$ are nonzero. The derivatives in $c_{n}$ can be expressed as expectation values of traces of operators involving $M^{-1}$ and $\partial_{\mu}M$, measured at $\mu=0$. For example, $c_{2}$ is evaluated using

[TABLE]

where the derivatives of $\ln\det M$ are given by

[TABLE]

Higher derivatives involving more and more terms and higher powers of $M^{-1}$ and $\partial_{\mu}M$ can be found in e.g. Allton et al. (2005).

Using the Complex Langevin equation we can simulate at nonzero chemical potential so the pressure is accessible as:

[TABLE]

where we have defined the charge density (using the space-time volume $\Omega=N_{s}^{3}N_{t}$)

[TABLE]

This means we can calculate the pressure at high chemical potentials at the cost of measuring the density at several chemical potentials in between and performing the integral (11). The density is a cheap observable with relatively small fluctuations. In contrast, for the Taylor expansion one needs to measure the $c_{n}$ coefficients at $\mu=0$, however these are quite costly, as they involve many inversions as $n$ increases, and they tend to be very noisy with increasing $n$, such that state of the art calculations can measure coefficients up to $c_{6}$ with a great effort Bazavov et al. (2017), but also other approaches exist based on imaginary chemical potentials D’Elia et al. (2017); Borsanyi et al. (2018). The extrapolated results to $\mu>0$ have error bars increasing such that they quickly lose predictive power above $\mu/T\sim 1$.

Once the pressure is calculated, i.e. the grand partition function is reconstructed, other thermodynamical observables can be calculated from it using various derivatives. The $\mu$ dependence of the density and fluctuations of the density can be directly measured in a simulation at the $\mu$ value of interest. Below the calculation of the energy density is detailed, the calculation of further quantities such as entropy density, speed of sound, charge susceptibilities, etc. is beyond the scope of this study.

The energy density $\epsilon$ can be accessed from the grand partition function through the trace anomaly

[TABLE]

where $\beta$ and $m$ are the bare parameters of the action, and this formula is also valid at $\mu>0$ Allton et al. (2003). As indicated, the derivatives in the formula above are understood to be defined along the line of constant physics (LCP), where the pion mass is kept fixed in physical units. In the importance sampling formulation eq.(13) is Taylor expanded in $\mu$ similarly to the pressure. For the first nonzero coefficient of the Taylor expansion (at the second order) one then measures the observables

[TABLE]

using

[TABLE]

with the chiral condensate $\chi=(N_{F}/4\Omega)\partial\textrm{ln}\textrm{det}M/\partial m$. In the complex Langevin setup however the $\mu$ dependence

[TABLE]

can be directly calculated using two simulations at the $\mu$ value of interest and at $\mu=0$. The observables needed for this calculation are the gauge action average $\langle S_{g}\rangle=-\partial\textrm{ln}\mathcal{Z}/\partial{\beta}$ and the chiral condensate $\chi=(\partial\textrm{ln}\mathcal{Z}/\partial m)T/V$. The beta function $a(\partial\beta/\partial a)_{\textrm{LCP}}$ and the derivative $\left({\partial m/\partial\beta}\right)_{\textrm{LCP}}$ can be estimated by independent simulations at zero temperature and $\mu=0$.

IV Stout smearing

To use stout smearing in Complex Langevin simulations, we must generalize its domain of definition from SU(N) to SL(3,$\mathcal{C}$) matrices, using a holomorphic function. The weighted staple sum corresponding to a link variable $U_{\nu}(x)$ is given by

[TABLE]

where $\rho_{\nu\sigma}$ are some real weights. Since we cannot use adjungation, we also need the sum of the inverse paths:

[TABLE]

We than define

[TABLE]

and finally the smeared field is given by $U_{\nu}^{\prime}(x)=e^{iQ_{\nu}(x)}U_{\nu}(x)$. This definition coincides with the usual stout smearing if the gauge fields are in SU(N), and the matrix $Q_{\nu}(x)$ is Hermitian in this case. On SL(3,$\mathcal{C}$) $Q_{\nu}(x)$ is no longer Hermitian but it is still traceless, so $U^{\prime}_{\nu}(x)$ is also an element of SL(N,$\mathcal{C}$). Typically one takes multiple smearing steps with

[TABLE]

and the measure becomes $e^{-S(U)}=e^{-S_{g}(U)}\textrm{det}\left(M(U^{(n)})\right)$, where $S_{g}(U)$ is the gauge action and $M(U)$ is the Dirac matrix describing the fermionic degrees of freedom.

For the calculation of the drift terms we need to evaluate $D_{a\nu x}S(U)$. Since the gauge part does not involve smearing we write $S(U)=S_{g}(U)+S_{f}(U)$ and we only consider the fermionic drift $D_{a\nu x}S_{f}(U)$ below.

Let’s consider one smearing step first where $U^{\prime}_{\nu}=e^{X_{\nu}}U_{\nu}$, and $i\lambda_{a}D^{\prime}_{a\sigma}S_{f}(U^{\prime})=F^{\prime}_{\sigma}$ is the standard force for unimproved fermions (with $D^{\prime}_{a\sigma}$ the left derivative with respect to variable $U^{\prime}_{\sigma}$), and the space-time indices are suppressed. Our aim is to calculate $F_{\nu}=i\lambda_{a}D_{a\nu}S_{f}(U)$, the force of the unsmeared field. For multiple smearing steps the procedure detailed below is repeated iteratively. For the drift term we will need to evaluate the left derivative $D_{a\nu}U^{\prime}_{\sigma}$, which can be represented as

[TABLE]

such that to first order in $\alpha_{a}$

[TABLE]

and the chain rule is written as $D_{a\nu}S[U]=D_{ab\nu\sigma}D_{b\sigma}^{\prime}S[U^{\prime}]$. The drift term is then written (also using the product rule $D_{a}(e^{X}U)=(D_{a}e^{X})U+e^{X}D_{a}U$ and the identity $-(i\lambda_{a}/2)\textrm{Tr}(i\lambda_{a}W)=W-(1/N)\textrm{Tr}W$)

[TABLE]

We write $\textrm{Tr}(e^{-X_{\sigma}}F^{\prime}_{\sigma}D_{a\nu}e^{X_{\sigma}})=\textrm{Tr}(L_{\sigma}(X,F^{\prime})D_{a\nu}X_{\sigma})$ where one can take $L_{\sigma}(X,F^{\prime})$ to be traceless (and anti-Hermitian for unitary link variables). Using the definition of $X$ we obtain (for isotropic smearing with $\rho_{\nu\sigma}=\rho$):

[TABLE]

Finally, to calculate the matrix $L_{\sigma}(X,F^{\prime})$ one can proceed using the following theorem Jennrich and Bright (1976): For a matrix $X$ of size $N\times N$, we write $\exp(tX)$ as

[TABLE]

Where $\lambda_{i}$ are the eigenvalues and $V$ is the matrix whose $j$th column is the eigenvector of $\lambda_{j}$. We than have

[TABLE]

where the cross-product is defined as $(G\times E)_{ij}=G_{ij}E_{ij}$ (no summation), and $G=V^{-1}(\partial X/\partial\theta)V$. The matrix $E$ is defined as

[TABLE]

This leads to $\textrm{Tr}(RD_{a}e^{X})=\textrm{Tr}(V((V^{-1}RV)\times E)V^{-1}D_{a}X)$. Alternatively, using the Cayley-Hamilton theorem any analytical function of a traceless $3\times 3$ matrix can be written as:

[TABLE]

where $f_{i}$ depends on the invariants of the matrix, $c_{0}=\textrm{det}X=\textrm{Tr}(X^{3})/3,\ c_{1}=\textrm{Tr}(X^{2})/2$ (recall that $\textrm{Tr}X=0$). Consequently the derivative is written as:

[TABLE]

Calculating the coefficients $f_{i}$ and their derivatives for the exponential function needed here proceeds by using a polynomial approximation to a fixed order ensuring correct results up to machine precision. Finally we write $\textrm{Tr}(RD_{a}e^{X})=\textrm{Tr}(BD_{a}X)$, with

[TABLE]

To check that the implementation is correct I have benchmarked the CLE results with results from the usual Hybrid Monte Carlo (HMC) implementation at $\mu=0$, see in Fig. 1. The comparison used the Symanzik gauge action and $n=2$ stout smeared staggered fermions with $N_{F}=4$ and isotropic smearing with $\rho=0.15$. The simulations were started from an SU(3) configuration, the observables are averaged between Langevin times $10<\tau<20$. Agreement within statistical errorbars is observed as long as the $\beta$ parameter is chosen large enough. For smaller $\beta$ values the gauge cooling becomes less effective, the unitarity norm rises quickly and the simulations become instable, just as it was observed for the naive action Fodor et al. (2015).

In Fig. 2 the effect of the smearing on a typical configuration from a CLE simulation is shown. The plaquette average nears 1.0 as it does also in the usual smearing of an SU(3) configuration. The unitarity norm of the configuration also increases slightly during the smearing procedure. If the initial unitarity norm of the configuration is higher $\gtrsim 0.1$, the smearing might cause a numerical overflow on the computer, especially if the $\rho$ parameter of the smearing is not small. This is similar to the ’runaway’ behavior known to occur in complex Langevin simulations. For the simulations in this study parameters are chosen such that this breakdown is very unlikely to occur.

V Results

Two actions used in this study, this gives a very rough estimate of the cutoff effects, and it allows for the testing of the stout staggered fermionic action with the Complex Langevin equation. First I use the plaquette gauge action with the naive staggered formulation using $N_{F}=4$ and the mass parameter $m=0.01$. Second the Symanzik gauge action is used with a stout smeared staggered action using $n=2,\ \rho=0.15,\ N_{F}=4,\ m=0.02$. The lattice spacing (measured with the $w_{0}$ parameter Borsanyi et al. (2012b)) and the mass of the lightest pion taste is shown in Tables 1,2.

In Fig. 3 the pressure difference (6) is shown for the naive ensemble for two different lattice spacings. To estimate the Taylor coefficients, $\approx 1000$ configurations were generated using a HMC simulation and on each configuration the $c_{n}$ were estimated using 128 noise vectors. The temperature of the system is above the deconfinement transition for both lattice spacings. The Taylor coefficients are listed in Table 3. Note that while in the continuum limit the Stefan-Boltzmann(SB) limit of $c_{2}$ is $N_{F}/2=2$, in the $N_{t}=8$ discretisation the SB limit is expected to be $\approx 2.8$ Allton et al. (2003). To apply the integration method (11) the integral is discretised with the stepsize $a\Delta\mu=0.025$, and CLE simulations are carried out at each chemical potential. The simulations used a partially second order update scheme Fukugita et al. (1987) with adaptive control of the Langevin stepsize Aarts et al. (2010b), using control parameters such that the timestep was typically in the range $(0.5-1)10^{-4}$. The simulations are started from a configuration where the link variables are initialized with white noise in SU(3) directions only. The thermalization of physical quantities such as the plaquette average, Polyakov loop average, density, etc. follows the expected exponential relaxation $\sim e^{-\tau/\tau_{0}}$ with $\tau_{0}<1$ for all parameters. 3 runs are used to collect averages for Langevin times $10<\tau<20$. The pressure is then reconstructed numerically and statistical errors are estimated using the jackknife method by splitting the stream of measurements to 10 pieces. Since the density as a function of the chemical potential is reasonably smooth at the high temperatures employed here, the systematic error coming from the discretisation of this integral is small (smaller than the statistical errors in this case), as can be estimated by employing different $\Delta\mu$ stepsizes. Quark chemical potentials up to $\mu=4T$ are used, this corresponds to $\mu a=0.5$. The Complex Langevin setup can be used for calculations at even higher chemical potentials, however cutoff effects will become important there. One observes good agreement of the Taylor expansion and the integration method. Note that while the errorbars of the pressure calculated from the integration method are small, the estimation of the coefficients of the Taylor expansion includes the systematic error corresponding to the choice of the fitting range.

In Fig. 4 the pressure difference is shown for the improved ensemble for two lattice spacings. The parameters were chosen such that the setup roughly corresponds to the same physical lattice spacings and pion masses as the setup using the unimproved action. To measure the Taylor coefficients $\approx 2000$ configurations from a HMC simulation were used with 64 noise vectors each. The Taylor coefficients are listed in Table 4. Using the improved action the importance sampling calculation of the $c_{n}$ coefficients is slightly less noisy such that the $c_{4}$ is also calculated with relatively small errors. The calculation of the $c_{6}$ coefficient can also be attempted, however since it is quite small only an upper limit on its magnitude is obtained. One observes good agreement of the Taylor expansion and the integration method, with relatively small discrepancy of the CLE and 4th order Taylor expansion results also at large $\mu/T$, suggesting that 6th and higher order terms have very small coefficients.

In Fig. 5 the quantity

[TABLE]

is plotted as a function of $(\mu/T)^{2}$. This allows for an intuitive way of judging the performance of the Taylor expansion. In this quantity, the second order term has a constant contribution, the fourth order term gives a linear behavior while the sixth order term adds a curvature. Note that at small $\mu$ the magnitude of the density is small, therefore the relative errors of ${n/(T^{2}\mu)}$ are larger.

To calculate the energy density, the estimation of the LCP and the beta function is necessary. Using zero temperature simulations (with HMC) at $\mu=0$ with slightly shifted $\beta$ values on $24^{3}\times 48$ lattices the needed mass parameters to keep the physical pion mass fixed are found by bracketing and using a chiral perturbation theory inspired ansatz for the fitting of the pion mass dependence on the quark mass. Using finite differences we get the following results:

[TABLE]

The $\mu$ dependence of the trace anomaly is given by a linear combination of the $\mu$ dependence of the gauge action average and the chiral condensate (see in eq. (13)). These quantities can be directly measured at nonzero $\mu$ using CLE simulations. Alternatively their behavior can be extrapolated using Taylor expansion from configurations at $\mu=0$. In Fig. 6 the $\mu$ dependence of the average gauge term and the chiral condensate term from eq. (13) is shown (omitting the extra factors of the beta function and the mass derivative) for the naive action at $\beta=5.3$ as well as for the improved action at $\beta=3.8$. One observes that the chiral condensate term is much better behaved with less fluctuations for the CLE as well as the Taylor extrapolation, and one observes good agreement for small chemical potential. The gauge action term however is much more noisy in both the CLE and Taylor expansion approach. With the amount of configurations at hand even the sign of the second order term is not clear. In Fig. 7 the final result for the $\mu$ dependence of $(\epsilon-3p)/T^{4}$ including the beta function and the mass derivative is presented. Since the mass derivative is rather small the results are dominated by the gauge action term, and thus the fluctuations are large. With the statistics at hand we can see that the dependence of the anomaly term on $\mu$ is significantly weaker than that of the pressure, but further conclusions are hard to gather from the data.

VI Conclusions

In this paper the thermodynamical properties of QCD at non-zero baryon density are studied. The aim of the study is to establish new methods offered by the availability of the CLE simulations at $\mu>0$ where naive importance sampling calculations are invalidated by the sign problem. The results are compared to the approach relying on the Taylor extrapolation of the results from the $\mu=0$ axis. As the Complex Langevin equation has potential problems at smaller temperatures related to the zeroes of the measure, here only the high temperature phase, namely the quark-gluon plasma state is investigated.

The pressure difference of the plasma between zero density and finite density states is estimated using an integration method where simulations are needed at intermediate points. The chemical potential dependence of the trace anomaly $\epsilon-3p$ is also calculated, this quantity is directly accessible in a CLE simulation at $\mu>0$.

Two lattice actions were investigated, the Wilson plaquette action with 4 flavors of naive staggered fermions and an improved action with the Symanzik improved gauge action and stout smeared staggered fermions. To this end the stout smearing procedure is generalized to the complexified SL(3,$\mathcal{C}$) manifold of the link variables. To reduce the cost of the simulations relatively heavy pion masses $\sim 500-700$ MeV are used.

Comparing with the usual Taylor expansion approach, good agreement is found in the small chemical potential region where the expansion is valid. The CLE approach can be used to calculate at higher chemical potentials as well, with relatively small errors. The results suggest that the 4th order expansion formula describes the dependence of the pressure on the chemical potential relatively closely, while the trace anomaly $\epsilon-3p$ remains approximately independent of the chemical potential for baryon chemical potentials up to $\sim 9T$.

The findings in this study show that the complex Langevin equation is a useful tool to access thermodynamic quantities, and allows calculations at high chemical potentials with small errorbars. To allow for direct applicability for the physical world some more work is needed: the continuum limit and infinite volume extrapolations still have to be carried out at the physical quark mass values.

Acknowledgements.

I would like to thank Manuel Scherzer, Erhard Seiler, Ion-Olimpiu Stamatescu for many discussions and collaboration on related topics, and Szabolcs Borsányi for discussions. I gratefully acknowledge funding by the DFG grant Heisenberg Programme (SE 2466/1-2), as well as the Gauss Centre for Supercomputing (GCS) for providing computer time on the supercomputers JURECA/BOOSTER and JUWELS at the Jülich Supercomputing Centre (JSC) under the GCS/NIC project ID HWU32. The research was partially supported by the BMBF grant no. 05P18PXFCA. Some parts of the numerical calculations were done on the GPU cluster at the University of Wuppertal.