Chiral symmetry restoration at finite temperature within the Hamiltonian approach to QCD in Coulomb gauge

The chiral phase transition of the quark sector of QCD is investigated within the Hamiltonian approach in Coulomb gauge. Finite temperature T is introduced by compactifying one spatial dimension, which makes all thermodynamical quantities accessible from the ground state on the spatial manifold $\small\mathbb{R}^2 \times S^1(1/T)$. Neglecting the coupling between quarks and transversal gluons, the equations of motion of the quark sector are solved numerically and the chiral quark condensate is evaluated and compared to the results of the usual canonical approach to finite-temperature Hamiltonian QCD based on the density operator of the grand canonical ensemble. For zero bare quark masses, we find a second-order chiral phase transition with a critical temperature of about 92 MeV. If the Coulomb string tension is adjusted to reproduce the phenomenological value of the quark condensate, the critical temperature increases to 118 MeV.


I. INTRODUCTION
Understanding the phase diagram of quantum chromodynamics (QCD) is still one of the most challenging problems in particle physics [1,2]. Lattice calculations can shed some light on its structure for vanishing baryon density but still suffer from the so-called sign problem in the general case of finite densities [1,3]. To overcome this problem, in the past two decades several non-perturbative continuum approaches, which do not suffer from the sign problem, have been developed [4], one of them being the variational approach to Hamiltonian QCD in Coulomb gauge [5], see ref. [6] for a recent review.
In ref. [7], the dressed Polyakov loop, the order parameter for confinement, and the chiral quark condensate, the order parameter for the spontaneous breaking of chiral symmetry, have been evaluated within this approach for vanishing chemical potential (i.e. baryon density). Thereby, finite temperatures were introduced by compactifying one spatial dimension using the alternative formulation of finite-temperature Hamiltonian quantum field theory proposed in ref. [8]. While the pseudo-critical temperatures of the chiral and, respectively, deconfinement phase transition were in good agreement with lattice data, the width of the transition region and the order of the chiral phase transition turned out to be at odds with the lattice predictions. This was suspected to be correlated to the neglect of the temperature dependence of the quark and gluon propagator, which were replaced by their zero-temperature limits to avoid the numerically highly expensive solution of the finite-temperature equations of motion.
In the present paper, we solve the quark part of these equations numerically. Thereby, we ignore the coupling of the quarks to the (transversal) spatial gluons. This corresponds to a confining quark model -the so-called Adler-Davis model [9] -which was considered in refs. [10][11][12][13] in the standard canonical formulation of finite-temperature quantum field theories. From our solution, we calculate the chiral condensate and compare it with the result of previous work.
The organization of the rest of this paper is as follows: In section II, we briefly review the essential ingredients of the novel approach to finite-temperature Hamiltonian quantum field theory developed in ref. [8] and its application to QCD in Coulomb gauge given in ref. [7]. The numerical solution of the quark equations of motion is described in detail in section III. The results for the mass function and the chiral condensate are presented in section IV, and we conclude the manuscript with a brief summary, some comments and an outlook on future directions in section V.

II. THE QUARK SECTOR OF FINITE-TEMPERATURE QCD
Below, we briefly discuss the main ingredients of the Hamiltonian approach to the quark sector of QCD when finite temperatures are introduced by compactifying a spatial dimension, for which we choose w.l.o.g. the 3-axis. For a more detailed description and a discussion of full QCD, the interested reader may consult refs. [7,8,14].
Let H be the QCD Hamiltonian in Coulomb and Weyl gauge on the compactified spatial manifold R 2 × S 1 (β), where β = 1/T denotes the inverse temperature. One can then show [8] that the grand canonical partition function Typeset by REVT E X arXiv:1806.04493v1 [hep-lat] 18 May 2018 at finite temperature T and chemical potential µ is given by where → ∞ is the length of the uncompactified spatial dimensions and E 0 is the smallest eigenvalue of the pseudo-Hamiltonian Here, α i denotes the usual Dirac matrices and ψ is the quark field which has to fulfill the anti-periodic boundary condition on the compactified manifold, while for the bosonic fields A the periodic condition holds (m and a are color indices in the fundamental and adjoint, respectively, representation). Furthermore, we have introduced the short-hand notation for the spatial integration. Let us stress that the novel finite-temperature Hamiltonian approach proposed in ref. [8] and leading to eq. (1) is equivalent to the familiar finite-temperature (imaginary-time) approach for any O(4) (i.e. relativistic) invariant quantum field theory. It is, however, advantageous in a Hamiltonian formulation in the sense that it does not require to explicitly carry out the thermal expectation values with the grand canonical density operator (with H being the Hamiltonian on R 3 and N being the fermionic particle-number operator) over the whole Fock space. Rather the thermal quantities like the partition function (1) are obtained from the vacuum state on R 2 × S 1 (β) alone (see below). Thus, the novel approach avoids introducing additional approximations to the Hamiltonian in the grand canonical density operator exp −β[H − µN ] , which is certainly an advantage. The novel approach is not manifestly spatial O(3)-invariant in the same way as the standard Hamiltonian approach based on canonical quantization is not manifestly Lorentz-invariant. However, the O(3) and O(4) invariance is hidden and recovered when the approach is carried out exactly. Let us also mention that in the novel approach the T → 0 limit can be easily taken after Poisson resummation, see Refs. [8] for more details. This fact will be exploited in the discussion following eq. (27) below. The QCD Hamiltonian H entering eq. (2) is given by [15] where is the quark single-particle Dirac Hamiltonian with g being the strong coupling constant, m Q the bare quark mass, γ 0 the usual Dirac matrix and t a the color generator in the fundamental representation. The second term in eq. (7) is the gluonic Yang-Mills Hamiltonian where Π = −iδ/δA is the canonical momentum operator (which agrees with the color electric field), is the color magnetic field and denotes the Faddeev-Popov determinant witĥ being the covariant derivative in the adjoint representation. Finally, is the so-called color Coulomb interaction which contains, besides the color density of quarks and gluons, the non-Abelian Coulomb kernel From eq. (1) it follows that all thermodynamical quantities can be obtained from the ground state |φ of the pseudo-Hamiltonian H which fulfills the functional Schrödinger equation H|φ = E 0 |φ [8]. Solving the functional Schrödinger equation is, thus, the aim of the Hamiltonian approach. On the compactified manifold R 2 × S 1 (β), this has been first tackled in ref. [16] for the Yang-Mills sector and was recently extended to full QCD in ref. [7]. Thereby, the ground state was calculated in an approximative way by using the variational principle: Using Gaussian type ansätze for both the bosonic and fermionic 1 parts of the vacuum wave functional |φ , the expectation value φ| H|φ was calculated on two-loop level. From its minimization, a set of coupled integral equations for the variational kernels contained in the ansatz for the wave functional |φ was obtained. While the so-called gap equation for the Yang-Mills sector was solved numerically in ref. [16], the full coupled equations were left unsolved in ref. [7] due to the high numerical expense. Instead, the zero-temperature propagators obtained in ref. [17] were used to calculate the dressed Polyakov loop and the temperature dependence of the chiral quark condensate for µ = 0. Remarkably, within these approximations the inclusion of the coupling of the quarks to the transverse spatial gluons showed only a negligible effect on the pseudo-critical temperatures of the deconfinement and chiral phase transitions.
In the present paper, we will give the numerical solution of the finite-temperature variational equations of motion for the quark sector and calculate the chiral condensate from it. Since the numerical cost is substantially higher for solving the full coupled equations, we will thereby neglect the coupling between quarks and transversal gluons. 2 Although it is not clear whether the effect of the coupling of the quarks to the transversal spatial gluons is still subleading at finite-temperature, this will enable us to study the effects of the temperature-dependence of the solution on the order and width of the chiral phase transition. Furthermore, it also allows for comparison between the compactified theory and the usual grand canonical approach to finite temperatures in Hamiltonian QCD considered in refs. [10][11][12][13].
Neglecting the coupling between quarks and transverse gluons, the fermionic part of the QCD Hamiltonian 3 reduces to where H 0 D [Eq. (8)] is the free Dirac Hamiltonian and H QQ C follows from the Coulomb term (13) after substituting ρ → ρ Q [Eq. (14)]. Note that this implies the cancellation of the Faddeev-Popov determinant in eq. (13). Furthermore, on two-loop level the non-Abelian Coulomb kernel can be replaced by its (Yang-Mills) vacuum expectation value, which plays the role of a confining quark potential, V C (|x − y|) = σ C |x − y| at |x − y| → ∞, where σ C is the Coulomb string tension [18]. Neglecting the coupling between quarks and transversal gluons, the ansatz for the fermionic part of the vacuum wave functional from ref. [7] reduces to the BCS-type functional where S is a scalar variational kernel, ψ ± denotes the positive/negative spectral projection of the quark field ψ and |0 is the bare vacuum of the Dirac sea, fulfilling ψ + |0 = ψ † − |0 = 0. This type of ansatz together with the Hamiltonian (16) corresponds to the confining quark model (Adler-Davis model) considered e.g. in refs. [9,[19][20][21] at zero temperature and in refs. [10][11][12][13]22] in the usual canonical approach to finite temperatures and densities. For explicit calculations, it is convenient to switch to the momentum space representation using where p ⊥ = p 1ê1 + p 2ê2 is the planar momentum and are the fermionic Matsubara frequencies resulting from the Fourier transformation of the (compactified) spatial component x 3 . Furthermore, we have introduced the short-hand notation In the following, we focus on the limit of vanishing chemical potential (µ = 0) and chiral quarks (m Q = 0). From the variational principle φ|H Q |φ → min one finds then the following integral equation for the variational kernel S 2NC is the value of the quadratic Casimir of the color group SU (N C ) [7] and p ⊥ = |p ⊥ |. For the numerical solution it is, however, more convenient to rewrite the scalar kernel S in terms of the effective quark mass function which transforms the gap equation (22) to Assuming the linearly rising form 4 V C (x) = σ C |x| for the non-Abelian Coulomb potential (17), its Fourier transform in the gap equation is given by The Coulomb string tension σ C entering this expression sets the overall scale in the present model. Lattice and continuum calculations [24][25][26][27][28] favour values σ C /σ ≈ 2 . . . 4 in terms of the Wilson string tension σ, with the rather large uncertainties coming from the extrapolation of the lattice Coulomb potential in the deep infrared. With the standard value σ = (440MeV) 2 for the Wilson string tension, this puts √ σ C in the range 650 MeV . . . 880 MeV. In the present work, we will use a standard value of √ σ C ≈ 700 MeV corresponding to σ C /σ ≈ 2.5, but we should be aware that this stipulation easily has uncertainties of up to 20%.
For a numerical evaluation, eq. (24) is not directly useful, since the entire calculation is dominated by the pole of the Coulomb potential for the single frequency Ω n = Ω l which -in constrast to the T = 0 equation discussed below -is not lifted by the integration measure. We thus have to introduce a small mass parameter µ (not to be confused with the chemical potential) to regularize the potential, and the entire calculation becomes very sensitive to the actual value of µ and the number of Matsubara frequencies included in the numerical code. We will present more details on the Matsubara type of gap equation (24) in section IV, but follow, for the main part of this paper, we follow a different route which also brings the underlying physics to the fore, viz. we Poisson resum the Matsubara series. For fermions, this is based on the simple distributional identity valid for a suitable test function f . If we use this equation to replace the Matsubara sum we find, after combining terms, shifting the loop momentum p → q ≡ p − k, and moving the Poisson sum outermost: For any index m, the integral under the sum is bound by the m = 0 contribution, i.e. the T = 0 limit. The zero temperature equation is, however, known to be both ultraviolet and infrared finite [7,17] and the same must hence hold for each integral in the Poisson sum (27) separately. (We will corroborate this assertion further below.) The infrared singularity of the Matsubara formulation now re-appears as convergence issue of the Poisson sum but, as we will demonstrate below, this issue can be handled analytically.
To close this section, we note that the ansatz (18) leads to the following expression for the chiral quark condensate:

III. NUMERICAL METHOD
In this section, we sketch the numerical techniques necessary to solve eq. (27). To fix our notation and discuss some numerical optimiziation, we briefly revisit the T = 0 case.

A. The zero temperature case revisited
The zero-temperature gap equation is simply the m = 0 contribution from eq. (27). To study it numerically, we measure all dimensionfull quantities in units of the mass scale As explained earlier, this stipulation has rather large uncertainties from the lattice calculations of σ C so will have all absolute numbers quoted in the present work. In the discussion, we will also present results for the quark condensate and the critical temperature, when m 0 is adjusted to match the lattice findings for the condensate at T = 0.
Next, we introduce spherical coordinates and exploit the rotational symmetry of the T = 0 system to eliminate the azimuthal angle. This gives where and we have indicated that the scalar mass function can only depend on k = |k| due to spherical symmetry. The prefactors in the mass scale eq. (29) were chosen such that all clutter is removed from the Coulomb potential, which now simply reads It is easy to see that the momentum integral in eq. (30) is ultraviolet convergent as long as M (k) is bounded at k → ∞. In the infrared, the superficial 1/q pole in the integrand disappears after integration over ξ, and the equation is infrared finite as well. However, solving eq. (30) by iteration is very unstable and requires substantial underrelaxation for convergence. As a consequence, a huge number of iterations (up to 20,000) is necessary to find the solution with high accuracy. To better understand this behaviour, it is convenient to rewrite eq. (30) in quotient form by collecting all pieces that contain the mass as a function of the external momentum, At first sight, this seems like a very bad way to rewrite the equation, since both the numerator and denominator are now infrared divergent. If we regularize the divergence by a lower cutoff µ to the momentum integral, it is easy to see that the leading µ → 0 contributions in numerator and denominator differ only by a factor M (k). Thus, the r.h.s. of eq. (33) reads, at small regulators, Iterating eq. (33) therefore produces changes to the mass function which are of order O(µ), i.e. the small infrared regulator also limits the speed at which the iteration progresses. This is what gives eq. (33) its inherent stability: any form of the gap equation in which the integrals are infrared finite will produce iteration changes of order O(1) which are way too large and must hence be substantially underrelaxed. By contrast, eq. (33) can be overrelaxed without loosing stability. As often, stability does not automatically imply efficiency: since eq. (33) makes very small progress in each step, a large number of more than 7,000 iterations is still necessary to solve it. This can, however, be cured by using sequence accelerators for the iteration, which improve convergence speed without sacrifying stability. In the present case, we have tested both a variant of Aitken's ∆ 2 -process [29] and Anderson's higher degree secant method [30]. Both algorithms must be vectorized, i.e. their transformation must affect the entire solution at all momenta k uniformly, because the convergence speed would otherwise differ at different k and the solution M (k) would become progressively distorted by the acceleration.
The accelerators generally buffer a certain number of iteration elements, and predict an improved estimator for the next iteration based on its history. Anderson's method, in particular, comes with a level k ≥ 1 that describes the dimension of the sub-space in which the univariate secant method is applied. It requires to store (k + 1) previous iterations, both accelerated and un-accelerated, and must solve a k × k linear system for each iteration. Usually, k = 2 and k = 3 give the best results, while levels k ≥ 6 rarely show any improvement. The Adler-Davis equation is different, however: Since it converges extremely slowly, the iteration can benefit from a much larger level k, combined with a moderate overrelaxation. We found that k = 18 with an overrelaxation α = 1.5 give the best results. The outcome is pretty impressive: The unaccelerated iteration requires more than 7, 000 iterations to reduce the residual As can be seen from the double logarithmic plot, the convergence of the accelerated sequences is less smooth but much faster than the standard iteration, and Anderson's method is clearly superior.
The resulting solution to the T = 0 equation is shown in the right panel of Fig. 1. As mentioned in the introduction, this mass function originates from the instantaneous part of the quark propagator in Coulomb gauge. It cannot be compared directly to the constituent mass function in Landau gauge, and attempts to match the two definitions reveal that the infrared limit M (0) ≈ 133 MeV observed here could still be compatible with the standard findings in Landau gauge [31]. The mass function computed here gives rise to a chiral condensate [7] if the standard scale eq. (29) corresponding to σ C /σ = 2.5 is used.

B. The finite temperature case: Poisson resummation
At non-zero temperatures, the presence of the heat bath singles out a rest frame and the original spatial O(3) symmetry is broken to O(2). As explained earlier, we have put the heat bath in the spatial 3-direction and compactified this dimension. The use of polar coordinates k = (k z , k ⊥ , ϕ k ) is, however, discouraged since the mass function would then depend on two non-compact coordinates k z and k ⊥ , which complicates the UV and, in particular, the IR limit considerably. A better strategy is to keep the spherical coordinates (k, ϑ k , ϕ k ). The remaining axial O(2) symmetry of rotations about the 3-axis entails that we can place the external momentum into the 1, 3-plane and set the azimuthal angle ϕ k = 0. Also, the mass function must be invariant under the reflection k z → (−k z ), as a remainder of the original O(3) symmetry 6 and we can take ξ k ≡ cos ϑ k ≥ 0 without loss of generality. The mass function is hence 5 Here and in the following, we are measuring the distance of solutions by the normalized L 2 metric, where {k i } are the grid positions on which the solutions are defined, and N is the number of grid points. (This formula can be extended to 2D grids for mass functions M (k, ξ k ) at finite temperature in an obvious manner.) The L 2 norm is a good compromise which measures convergence on average while still giving each individual grid position enough weight so that outliers will not go by unnoticed. 6 This corresponds to the reflection n → −(n + 1) of the Matsubara indices which flips the sign of the Matsubara frequency.
For the loop integration, we also adopt spherical coordinates (q, ϑ, ϕ). The angles only enters through their cosine via the scalar product q · k = q z k z + q ⊥ · k ⊥ = q z k z + q ⊥ k ⊥ cos ϕ = qk cos ϑ cos ϑ k + cos ϕ 1 − cos 2 ϑ 1 − cos 2 ϑ k = qk ξξ k + η , where we have introduced the cosines and defined the useful abbreviation In these coordinates, the Poisson re-summed gap equation (27) in quotient form becomes with the shifted momentum If we compare this to the T = 0 version in eq. (30), it is evident that the m = 0 term in the Poisson sums reproduces the T → 0 limit, if we assume that the mass function does not depend on ξ k because of the restored O(3) symmetry. Also, the shifted momentum eq. (42) agrees with the T = 0 limit in the equation below (30) when ξ k = 1, i.e. when the external momentum points into the direction of the heat bath. This direction of the external momentum therefore gives the closest analogue of the T = 0 mass function, and we will therefore compare M (k, ξ k = 1) to the T = 0 limit in section IV below. Eq. (41) is not yet suited for numerical investigation. As we have explained in the previous section, a small regulator µ for the infrared divergence of the quotient form is required and provides for a stable iteration, At finite temperatures, however, the infrared divergence also leads to a poor behaviour of the Poisson series, whose terms typically decay very slowly when µ 1. It is therefore convenient to subtract an analytic helper function in the integrands of eq. (41), which will render the q-integral IR finite at q → 0. Of course, we have to add back in what we subtracted, and this extra term will now carry the infrared divergence and the poor Poisson sum. The advantage of this procedure is that the part of the calculation which depends on the (numercially expensive) mass function is finite and quickly converging, while the problematic terms can be handled analytically. The gap equation now takes the form where the integrands read and the subtractions are compensated by the inhomogeneities g and h. An obvious choice for the subtractions is the q = 0 limit of the integrands, ∆u(q, ξ, γ ; k, ξ k ) = M (k, ξ k ) ∆v(q, ξ, γ ; k, ξ k ) = 1 With this subtraction, the integrands u and v behave as O(q) at small momenta, but the leading O(q) term is linear in γ and thus γ-integrates to zero. The result of the γ-integration is hence O(q 2 ) which, together with the factor q 2 U (q) ∼ q −2 , yields a finite loop integral at q → 0, even in the absence of a cutoff. The infrared divergence now re-appears in the inhomogeneities where it can be handled analytically: after performing the integrations with a small infrared regulator µ > 0 in the potential U (q) as in eq. (43), we obtain .
Note that the full Poisson sum appears to be finite when the regulator is removed, while the m = 0 term is 1/(2µ) and thus really diverges in the infrared. As explained in the previous section, both numerator and denominator of the gap equation in quotient form should be infrared divergent at T = 0 and this property should also persist at T > 0. This apparent contradiction comes from an order of limits issue: since the relevant terms in eq. (48) have the argument βµ, the limit µ → 0 at finite β gives a finite result, while the limit β → ∞ at finite µ yields the expected T = 0 divergence 1/(2µ). This indicates that the correct formulation (the one which is continuously connected to the T = 0 case) must retain a small, but finite IR cutoff µ > 0. Taking µ → 0 too early will lead to a different (finite) formulation in which the T → 0 limit disagrees with the well-known T = 0 result. We will thus always keep a small but non-zero IR cutoff µ > 0 which ensures a smooth limit T → 0 and also stabilizes the iteration as explained in the previous section.
With the Fourier integrand going as O(q 0 ) at q → ∞, simple dimensional analysis suggests that the terms in the Poisson sum decay as (βm) −1 which, together with the alternating sign, amounts to a poorly converging series. 7 This type of slowly converging alternating series can, however, be handled quite efficiently using the -algorithm [32,33]. The alternative would be to attempt one more subtraction of the O(q 2 ) behaviour under the integrands u and v. This leads, however, to a rather formidable expression involving up to second order derivatives of the mass function. Since the latter is only known numerically on a rather coarse momentum grid, the second subtraction cannot be carried out with sufficient accuracy and we stick to eq. (46).
Eq. (44) is the final form of the gap equation which we solve iteratively: we start with an arbitrary function M 0 (k), either the T = 0 solution or a constant, and use Anderson's algorithm as a sequence accelerator as in the T=0 case, cf. fig. 1. Once the system has been iterated to convergence, we can extract the chiral condensate from which is the finite temperature extension of eq. (36).

C. The finite temperature case: Matsubara formulation
The poisson resummation technique described in the last section is convenient at low temperatures were only a few terms are required. In addition, the T = 0 limit is recovered from the lowest term of the series. As the temperature increases, more and more terms of the Poisson series have to be included. At very high temperatures, we may thus reach a point were the original Matusbara formulation is more convenient, as the relevant sums as saturated by the first few Matsubara frequencies. In particular, only the lowest Matsubara frequency survives in the high-temperature limit T → ∞.
Though our numerical procedure mainly relies on the Poisson technique, we have also solved quark gap equation (24) in the Matsubara representation and compared it with the results of the Poisson formulation. This will provide an independent test for the accuracy of our numerics.
For the Matsubara formulation, we employ the residual O(2) symmetry of rotations about the 3-axis to let the external momentum component in the plane perpendicular to the heat bath point into 1-direction, k ⊥ = k ⊥ e 1 . For the loop integration, we use polar coordinates q ⊥ and ξ = cos (q ⊥ , e 1 ). We can then express eq. (24) in these coordinates, scale all dimensionfull quantities in the units of eq. (29) and finally go over to the more stable quotient form. This gives where Q ⊥ ≡ k 2 ⊥ + q 2 ⊥ + 2k ⊥ q ⊥ ξ. Because of the symmetry M (k ⊥ , Ω ) = M (k ⊥ , −Ω ) = M (k ⊥ , Ω − −1 ), we can restrict the external Matsubara index to ≥ 0. On the rhs of eq. (52), we have also combined the terms with index n and −(n + 1), since they only differ in the sign of the frequency Ω n . This ensures that only mass functions with Matsubara index n ≥ 0 appear on both sides of eq. (52) and the coupled integral equation system closes. The chiral condensate can be expressed in the Matsubara formulation as For a numerical evaluation, the potential U (q) must be infrared regularized as in eq. (43), and the number of Matsubara frequencies included in the system must be restricted to n < N . The system (52) then resembles the T = 0 equation, however with an N -component solution M (k ⊥ ) ≡ M (k ⊥ , Ω ) and a different integration measure. It is the latter property which makes the Matsubara formulation less convenient: for small regulators µ 1, the potential U has a strong singularity at q ⊥ = 0 if the external and loop frequency match, Ω = Ω n . This singularity is only partially cancelled by the integration measure and the n = term dominates the entire Matsubara sum by a relative factor 1/µ 2 . As before, this singular factor is canceled between the numerator and denominator, and the remaining O(µ 2 ) contributions from the other terms n = in the Matsubara sum carry the actual corrections to the mass functions. This means that the iteration progresses much slower than at T = 0 and, more problematic, the n = terms from the Matsubara sums must be computed to a very high accuracy. In addition to the high accuracy demand, the Matsubara formulation has the property that all components M (k ⊥ , Ω ) are coupled by the system (52), so that the index cutoff N must be fixed once and for all and cannot be adjusted dynamically. 8 If the external index approaches the cutoff N , there are only a few frequencies larger than the dominating contribution n = , i.e. the Matsubara sum is truncated unsymmetrically and the higher frequencies ≈ N thus have a systematic bias. The only solution is to include a very large number of frequencies N , so that the inaccurate modes near the cutoff give such a small contribution that their combined error does not matter. Unfortunately, the computational effort of the system (52) scales strictly as O(N 2 ), so that the inclusion of higher frequencies is limited by practical considerations. In our studies, we were able to push the Matsubara mode count to N = 100 for the lower temperatures, which means that we have 10, 000 times the numerical effort of the T = 0 solution per iteration, plus a much larger iteration count due to the slow convergence of eq. (52) and an increased accuracy demand. Even with this considerable effort, the frequency count was just enough to map the transition region, but temperatures below T = 50 MeV give incorrect results and require a different (massively parallel) strategy, cf. Fig. 7. By contrast, the Poisson formulation -though much more costly per iteration -scales better with increasing temperature as it can adjust its mode count dynamcially and hence will always give a correct result, albeit with a (moderately) increased effort at higher temperatures.
Finally, it should also be mentioned that the convergence of the Matsubara system (52) is non-uniform in the frequency index, i.e. the lowest frequencies are stable after a relatively small iteration count, while the highest frequencies (which contribute the least) are the slowest to converge. In practice, we have to stop the iteration at some point where the highest frequencies may not have fully converged. Since we relax from above, this means that the highest frequency mass functions are systematically too large, and, although each frequency contributes very little, their combined effect may lead to overestimate the condensate eq. (53). We will see this effect in Fig. 7 below, where the Matsubara values for the condensate in the transition region are all slightly larger than the corresponding Poisson results.

IV. RESULTS
We split our result section in two parts: first we consider numerical details on the individual parts of our calculation to demonstrate that the fairly complicated process actually works as intended. In the second part, we discuss the final results for the mass function and the chiral condensate at different temperatures. 8 It is possible to increase the number of frequencies included in the sum by assuming that O(3) invariance is restored for large frequencies.
This allows to approximate M (k ⊥ , Ω ) ≈ M (k ⊥ , Ω N −1 ) for large frequencies ≥ N , wherek ⊥ = k 2 ⊥ + Ω 2 − Ω 2 N −1 > k ⊥ . For small N , this approximation is not applicable, and for large N it is unnecessary since the sum will converge without it. The extrapolation is hence most useful in the intermediate region N ≈ 40 where it can save a factor of 2-4 in computation time.

A. Details on the numerical method
In the following, we present some typical results of intermediate steps in the calculation. Numerical issues appear predominantly in the earlier steps of the iteration and the eventual mass function has a similar shape to the initial zero-temperature solution, cf. below. For simplicity, we can therefore assume the zero-temperature mass function for the qualitative arguments in this subsection. Furthermore, we fix the external momentum to a typical value k = 200 MeV and ξ k = 0.5, which is in the region where the mass function changes most quickly.
It should also be noted that we generally combine terms with both signs ±ξ in all internal calculations, i.e. we symmetrize the ξ-integrand To keep the formulas simple, we will not always indicate this symmetrization, which is implicitly understood. We begin with the integrand of the momentum integral omitting the Fourier cosine factor for clarity and combining terms from ξ and (−ξ), The integration here can be done very efficiently using Gauss-Chebychev integration, which automatically takes care of the square root factor in the denominator. In Fig. 2, we have plotted eq. (55) for two values of ξ close to the boundary, for both the numerator (left panel) and the denominator (right panel) of the gap equation (44). The shape of the functions is quite similar in all cases: for small momenta q, the functions approach a constant, which is maintained for about two to three orders of magnitude, before the influence of the regulator µ = 0.1 MeV sets in and the functions quickly vanish for q < µ. For the numerator, the functions with the larger ξ lies above the one with the smaller ξ, while this order if reversed in the denominator. The general shape of all these functions is in agreement with our discussion of the subtraction above. The actual integrand of the momentum integral still has the Fourier cosine factor, which leads to the oscillating functions depicted in Fig. 3.
Next, we consider the integrand of the ξ-integral after performing the Fourier momentum integration, Note that the argument ξ in this function can be restricted to ξ ∈ [0, 1] due to the symmetrization of ±ξ. Besides the external momentum (which we have fixed to the same standard value as in the previous figures), this function now depends on both the temperature and the Poisson summation index. In the left panel of Fig. 4, we have plotted eq. (56) for a fixed temperature T = 50 MeV and two Poisson indices m = 1 and m = 10. As can be seen, the function eq. (56) is non-oscillating and rather smooth, except for a steep drop in the vicinity of ξ = 0. This is the region where the Fourier momentum integral has a low frequency and thus the cancellations due to rapid oscillations are absent. Numerically, an accurate integration of eq. (56) requires a large number of sampling points if a uniform ξ-sampling is adopted. Alternatively, it is more efficient to spread out the low ξ-behaviour by a change of variables ξ = t n with n > 1, In the right panel of Fig. 4, we show the n = 4 transformation eq. (57) of the plots in the left panel. The detailed structure at low ξ is now spread out and the resulting function can be accurately integrated using Gauss-Legendre with about 30 sampling points. (By contrast, more than 1200 uniformly distributed sampling points were used for the left panel of Fig. 4). The absolute value of the ξ-integrals in the left panel of Fig. 4  Finally, we check the convergence of the Poisson sums in the gap equation (44). We fix the external momentum again at our preferred value k = 200 MeV and ξ k = 0.5, and plot the partial Poisson sums in eq. (44), including the prefactor 1/π 2 . The sums are even in m, i.e. we can combine terms with ±m, . These oscillations are much more pronounced for the higher temperature in the right panel of Fig. 5. Numerically, we have found that up to 10, 000 terms would have to be summed in this case, in order to suppress the oscillations and predict the value of the infinite series to a relative accuracy of 10 −4 . By contrast, the -algorithm is able to reach the same accuracy from only the first 25 terms in the series.

B. Results
In the left panel of Fig. 6, we show the mass function M (k, ξ k ) for the two extreme directions ξ k of the external momentum (longitudinal and perpendicular to the heat bath). As explained earlier, the closest analog of the T = 0 mass function M 0 (k) in our finite temperature formulation is M (k, ξ k = 1), when the external momentum points in the direction of the heat bath. This curve has indeed a very similar form to the T = 0 solution plotted for comparison, while the ξ = 0 curve shows some deviations. Also, the mass function is already considerably smaller than at T = 0, even though the temperature T = 80 MeV in this plot is still in the confined phase. In order to recover the T → 0 limit, we would thus have to go to very small temperatures. This is shown in the right panel of Fig. 6, where we compare the mass function for small temperatures with the T = 0 limit. As can be seen, the mass function starts to drop already at temperatures below 10 MeV. From eq. (51), this does not directly translate into a drop of the condensate, since the mass function appears in the numerator and denominator of the integrand, which is thus less sensitive to the mass function at small momenta. (The temperature also appears through the Fourier sum which gives the main temperature sensitivity.) For these reasons, the intercept M (0, 1) of the mass function is not a good indicator for the phase transition, in particular since it is also gauge dependent. This is also the reason why it cannot be directly compared to the constitutent mass in conventional covariant gauges; instead, it represents a gauge-dependent mass parameter from the instantaneous part of the quark propagator in Coulomb gauge. In Ref. [31] it was demonstrated that such a mass parameter in Coulomb gauge could be as low as M (0) ≈ 100 MeV and still be compatible with the (much larger) constituent masses in Landau gauge.
The gauge-invariant order parameter for the phase transition is the chiral condensate plotted in Fig. 7. This shows the expected behaviour, i.e. it is roughly constant for small temperatures, and drops quickly to very small values at a characteristic transition temperature T * . For higher temperatures, it is practically zero within our numerical precision. 9 The exact location of the phase transition temperature T * depends on the order of the transition: If it is a strong cross-over, one usually defines T * from the inflexion point, while a second order transition would have a jump discontinuity in the derivative at the temperature where a non-vanishing condensate appears for the first time as we cool from the deconfined phase. This value of the transition temperature is always larger than the inflexion point. Our data indeed indicates at a (weak) second order transition with a critical temperature of Both the order and the critical temperature differ from the result of our previous investigation in ref. [7], where a broad crossover phase transition with pseudo-critical temperature 165 MeV was obtained in the limit of a vanishing quark-gluon coupling. It should, however, be mentioned that the solution of the zero temperature gap equation was used at all temperatures in Ref. [7], which is certainly a quite crude approximation. In Fig. 7, we have also included data from a direct summation of the original Matsubara series (24) with the same cutoff µ = 0.1 MeV to take care of the infrared-singularity and the transformation to quotient form for better stability. This system scales quadratically with the number of included Matsubara frequencies, and up to 100 frequencies (i.e. 10,000 times the effort as compared to T = 0) were necessary to even reach the transition region. 10 By contrast, the Poisson formulation, though much more expensive per iteration, has no infrared singularity and provides the correct results even at larger temperatures, where the increase in computational effort as compared to lower temperatures is still moderate. Combined with the physical transparency of the method, this warrants the high numerical effort of the Poisson resummation.
Returning to the results for the chiral phase transition, it should be emphasized again that the absolute value of T * as well as the size of the condensate depends on the overall scale set by the Coulomb string tension, cf. eq. (29). The cited values are for the preferred value σ C = 2.5 σ. In view of the uncertainties about the fundamental scale, it is better to cite our findings for the critical temperature as For instance, a somewhat larger value σ C /σ = 4.1 would reproduce the lattice results for the chiral condensate and push the transition temperature to around T * ≈ 118 MeV. A second-order chiral phase transition is the expected result for a system of two chiral quark flavors, as can be seen from the so-called Columbia diagram [1,2]. This might be surprising since we consider only one single quark flavor within our model. One should, however, notice that as long as the variational kernel S is flavor-diagonal, our results would not change even if we would take more flavors into account. Since the neglect of the bare quark masses is the best approximation for two quark flavors (i.e. up and down), a second-order phase transition is definitely meaningful for the considered model. Especially, it also agrees with the result found when finite temperatures are introduced within the grand canonical ensemble, see e.g. refs. [12,13]. The introduction of a flavor-dependent variational kernel might become necessary as soon as finite current quark masses and unquenching effects are incorporated into the model. This could yield a dependence of the order of the phase transition on the number of quark flavors, as one would also expect from the Columbia diagram.
While the order of the chiral phase transition is the same in both the canonical and the present approach to finite temperatures, the critical temperature T * can ≈ 64 MeV found 11 in the numerical calculations of ref. [13] for vanishing chemical potential µ = 0 is significantly smaller than our result. There are several possible reasons for this discrepancy: On the one hand, the evaluation of the partition function [see eq. (6)] where H is the QCD Hamiltonian on R 3 , necessitates some approximations in the canonical approach. This concerns especially the treatment of the density operator [Eq. (6)] of the grand canonical ensemble, where a quasi-particle approximation is required for the Hamiltonian H. 12,13 Note that such an approximation is not required in the present approach and no equations of motion for the quasi-particle energies hence emerge. On the other hand, however, the present approach is based on the finite temperature formalism of Ref. [8], which in turn relies on the O(4)-invariance of the Lagrangian -this is not entirely fulfilled in the Adler-Davis model since it contains the fermionic parts of the Coulomb interaction (13) -which is related to the A 0 -A 0 correlator [34] -but no contribution from the spatial gluons. Nevertheless, the value for the critical temperature eq. (60) is too small as compared to the one found in lattice simulations for physical quark masses, T * lat ≈ 155 MeV [35]. Note that also the gauge invariant chiral quark condensate (Fig. 7) falls behind the value expected from phenomenology. The fact that the Adler-Davis model predicts significantly too small results, e.g. for the pion decay constant, is well-known and also obtained in the canonical approach [36]. In this respect, we should stress that an increase of the critical temperature, and other quantities as the condensate, is expected when the coupling of the quarks to the transverse spatial gluons is included, as the numerical calculation performed in Ref. [7] shows. Such a more complete system would also allow for a more reliable comparison between the results obtained within the novel and the canonical approach to finite-temperature Hamiltonian QCD.
Unfortunately, the inclusion of the coupling to the transversal spatial gluons will drastically increase the numerical costs, which is already considerable in the present study. Although the coupling effects turned out to be small when the solution of the zero-temperature gap equation is used [7], it is not clear if this is still true in the full temperature-dependent calculations. In contrast, the inclusion of finite bare quark masses should only smear out the phase transition to a crossover without having major effects on the value of the critical temperature.
Finally, one should emphasize again that the Coulomb string tension could also be adjusted to σ C /σ ≈ 4.1, which reproduces the phenomenological value of the quark condensate, ψ ψ = − 236 MeV 3 , but is somewhat larger than the lattice prediction. With this arrangement, one would find a critical temperature of T * = 118 MeV.

V. SUMMARY AND CONCLUSIONS
In the present paper, we have revisited the alternative Hamiltonian approach to finite-temperature QCD of Ref. [8] and solved the temperature-dependent equations of motion of the fermion sector numerically. In a first study, we have ignored the coupling of quarks and transverse spatial gluons. The apparent infrared singularity of the resulting model could be resolved using Poisson resummation of the original Matsubara series, and the resulting integral equation system was solved using refined numerical techniques combined with analytical computations. Even though the iteration is stable and amenable to series accelerations, the loss of O(3) symmetry requires a Poisson sum and three nested integrations per temperature, of which the Fourier quadrature, in particular, is fairly expensive. As a consequence, the final solution amounts to 100+ CPU hours per temperature.
For zero bare quark masses, the results for the chiral quark condensate show the expected weak second-order phase transition with a critical temperature of T * ≈ 92 MeV. While this value is larger than the one found within the usual (canonical) approach to finite temperature Hamiltonian QCD [13], our findings for the critical temperature are definitely smaller than the result of lattice calculations using dynamical quarks, T * lat ≈ 155 MeV [35]. We suspect that the mismatch between our findings and the lattice results is related to the neglect of the quark-gluon coupling in the variational ansatz, which also leads to a value of the quark condensate which is too small. However, if the scale is adjusted to reproduce the physical value of the quark condensate instead of fixing the scale from the rather poorly determined Coulomb string tension, the critical temperature increases to about 118 MeV, which is much closer to the lattice data (while σ C /σ is still in the range supported by lattice calculations.) In forthcoming investigations, we therefore intend to study how this coupling affects the chiral phase transition. The solution of the full coupled equations of motion will also allow for the evaluation of the dressed Polyakov loop as order parameter for confinement. Furthermore, we plan to extend our calculations to the general case of a finite chemical potential, µ = 0, in order to obtain a description of the whole QCD phase diagram within the Hamiltonian approach.