Strong Coupling Lattice QCD in the Continuous Time Limit

We present results for lattice QCD with staggered fermions in the limit of infinite gauge coupling, obtained from a worm-type Monte Carlo algorithm on a discrete spatial lattice but with continuous Euclidean time. This is obtained by sending both the anisotropy parameter $\xi=a_\sigma/a_\tau$ and the number of time-slices $N_\tau$ to infinity, keeping the ratio $aT=\xi/N\tau$ fixed. The obvious gain is that no continuum extrapolation $N_\tau \rightarrow \infty$ has to be carried out. Moreover, the algorithm is faster and the sign problem disappears. We derive the continuous time partition function and the corresponding Hamiltonian formulation. We compare our computations with those on discrete lattices and study both zero and finite temperature properties of lattice QCD in this regime.

The determination of the QCD phase diagram, in particular the location of the critical end point (CEP) is an important, long standing problem, requiring nonperturbative methods. In lattice QCD, several approaches have been developed to investigate the phase transition from hadronic matter to the quark gluon plasma, but either they are limited to rather small µ B /T , with µ B the baryon chemical potential [1][2][3], or they cannot yet address full QCD [4][5][6] or study only low dimensional QCD-like toy models [7][8][9].
The reason for this is the notorious sign problem [10], which arises because the fermion determinant for finite µ B becomes complex, and importance sampling is no longer applicable. In lattice QCD, the finite density sign problem is severe. There is however a limit where the sign problem can be made mild: this is the strong coupling limit, where a so-called dual representation in terms of color singlets is possible. In the strong coupling limit of lattice QCD (SC-LQCD) the sign problem is mild enough such that the full (µ B , T ) phase diagram can be measured via Monte Carlo methods based on the dual variables. The method of dual variables has been successfully used in models with Abelian gauge symmetry [11,12], and there have been attempts to dualize non-Abelian gauge theories [13,14], but it has not yet been possible to overcome the finite density sign problem. Our own approach discussed in [15][16][17] is based on the strong coupling expansion, i.e. an expansion in the inverse gauge coupling β = 2Nc g 2 . It is in principle possible to sample partition functions that include all orders via Monte Carlo, in a spirit of [18,19]. In practice, the sign problem is reintroduced for large β.
In this paper, we will restrict to the strong coupling limit, since the focus is on deriving the Euclidean continuous time limit and apply the new formulation to Monte Carlo studies of QCD thermodynamics. Despite the fact that the strong coupling limit is the converse of the continuum limit, i.e. the lattice is maximally coarse and it is not possible to set the scale, it nevertheless shares important features with lattice QCD on finer lattices: chiral symmetry breaking and its restauration at finite temperature as well as the nuclear liquid gas transition are also present in this model. We will extend the existing studies on SC-LQCD that are either based on mean field theory in the 1/d expansion [20][21][22][23][24][25][26] or on Monte Carlo [27][28][29]. In the past either the spectrum or the phase diagram and the nuclear properties [29] have been studied. We investigate these phenomena in the continuous time limit, where the continuum limit of the temporal lattice spacing a τ → 0 is taken while leaving the spatial lattice spacing a σ finite. First simulations of SC-LQCD in continuous time have been performed by one of us in [30]. Here, we improve upon the continuous time formulation and give many more results at zero and non-zero temperature. The main advantage of the continuous time limit (CT) is that ambiguities arising from the anisotropy parameter γ are circumvented. Also, the sign problem is absent, Quantum Monte Carlo (QMC) can be applied, and temporal correlation functions can be obtained with high resolution. This paper is structured as follows: in Sec. II we will derive the Quantum Hamiltonian formulation of strong coupling QCD and its generalization to an arbitrary number of colors. In Sec. III we will describe the worm algorithm operating in continuous time in detail and show that it indeed reproduces results consistent with the continuum extrapolation of simulations at finite N τ . In Sec. IV we apply SC-LQCD in the CT-limit to determine zero temperature observables. In Sec. V we investigate finite temperature properties, such as the grand-canonical phase diagram in the µ B − T plane as well as the canonical phase diagram in the n B − T plane, with n B the baryon number density. In Sec. VII A we discuss temporal correlation functions and how to extract pole masses. We provide both results at finite temperature and density. In Sec. VI A we show that the pressure at finite baryon density can also be reconstructed from Taylor coefficients. We conclude with remarks on the radius of convergence. In the appendix, supplementary material for the various crosschecks of continuous time Monte Carlo and possible extensions such as for finite quark mass, more flavors and isospin chemical potential are discussed.

A. Staggered Action of Strong Coupling QCD and its Dual Representation
In SC-LQCD, based on the Euclidean lattice action, the gauge coupling is sent to infinity and thus the coefficient of the plaquette term β = 2N c /g 2 is sent to zero. Hence the Yang Mills part F µν F µν is absent. Then, the gauge fields in the covariant derivative can be integrated out analytically. In fact, the order of integration is reversed compared to the standard representation of lattice QCD in terms of the fermion determinant: the gauge links U µ (x) are integrated out before the Grassmann fields χ,χ. Thus the final degrees of freedom of the partition function are color singlets composed of fermions: mesons and baryon. However, as a consequence of the strong coupling limit, the lattice becomes maximally coarse and there is no way to set the scale: the lattice spacing a cannot be specified in physical units. We will see however that specific dimensionless ratios can still be compared to continuum physics.
We shortly outline the procedure to obtain the dual representation for staggered fermions in the strong coupling limit where the action is only given by the fermionic part: S[U, χ,χ] = x γ η 0 (x) χ(x)e aτ µq U 0 (x)χ(x +0) −χ(x +0)e −aτ µq U † 0 (x)χ(x) Here, am q is the quark mass and µ q = 1 3 µ B the quark chemical potential. The bare anisotropy parameter γ in the temporal Dirac coupling is introduced to vary the temperature continuously.
Following the procedure discussed in detail in [28], the gauge link integration over the Haar measure of SU(N c ) can be performed analytically, as the integration factorizes in Eq. (1), i. e. the partition function can be written as a product of one-link integrals z µ (x): M (x) =χ(x)χ(x), The new degrees of freedom after link integration on the right-hand side are the mesons M (x) and the baryons B(x). The weight of the one-link integral is a sum over the so-called dimer number k µ (x) = 0, . . . N c which corresponds to the number of (non-oriented) meson hoppings on that link, and on ρ(x, y) which is the weight for a baryon hoppingB(x)B(y). The final partition function for the discrete system on a N σ 3 × N τ lattice, after performing the Grassmann integrals analytically, is an exact rewriting from Eq. (1) and is given by: b=(x,μ)∈ ημ(x) The sum over all configurations {k, n, } is restricted to those that fulfill on each site x the so-called Grassmann constraint (GC): which expresses the fact that every Grassmann variables χ i (x), χ i (x) (i = 1 . . . N c ) appear exactly once in the path integral. After this exact rewriting of the strong coupling partition function the system can be described by confined, colorless, discrete degrees of freedom: • Mesonic degrees of freedom: kμ(x) ∈ {0, . . . N c } (non-oriented meson hoppings called dimers) and n(x) ∈ {0, . . . N c } (mesonic sites called monomers).
• Baryonic degrees of freedom: they form oriented baryon loops and may wind ω( ) times in temporal direction, which results on its dependence on the chemical potential µ q . The sign σ( ) = ±1 of the loop depends on the loop geometry.
• The baryonic loops are self-avoiding and do not touch the mesonic degrees of freedom, which follows from the Grassmann constraint Eq. (6): for a given configuration, this gives rise to a decomposition of the lattice volume into mesonic sites and baryonic sites: It should be mentioned that this representation corresponds to unrooted staggered fermions. Due to the fermion doubling, one flavor of a staggered fermion comes in the multiplicity of four so-called tastes. However, in the strong coupling limit, the fermions are spinless and the taste breaking is maximal. Hence it is indeed a one-flavor theory with only one pseudoscalar meson as the Goldstone boson.
To be more precise, in the chiral limit the action is invariant under the symmetry group U B (1) × U 55 (1): χ(x) → e iθ B +i (x)θ55 χ(x), (x) = (−1) µ xµ (8) which is due to the even-odd decomposition of the bipartite lattice for staggered fermions, i. e. even and odd sites can be transformed independently. The symmetry e iθ B ∈ U B (1) corresponds to baryon conservation and e iθ55 ∈ U (1) 55 is a subgroup of the full SU L (4) L ×SU R (4) chiral symmetry for unrooted staggered fermions. In the spin-taste basis this corresponds to the channel γ 5 ⊗ ξ 5 . At finite quark mass U (1) 55 is explicitly broken, and in the dual representation this is due to the presence of monomers: the number of monomers on even sites equals its number on odd sites. In the chiral limit we expect O(2) critical exponents for the chiral phase transition. This is also the case away from the strong coupling limit, as long as the lattice spacing is finite. In this work we will restrict to the chiral limit, m q = 0, where monomers are absent: n x = 0. We discuss the prospects of the continuous time formulation at finite quark mass in the appendix IX E.

B. SC-LQCD at Finite Temperature and the Continuous Time Limit
In the staggered action Eq. (1) we have introduced a bare anisotropy γ in order to vary the temperature continuously. Hence also in the dual representation the weights for temporal meson or baryon hoppings in Eq. (5) contain the anisotropy parameter γ. We will now explain why this is necessary and why it is also a key step to derive the continuous time limit.
The main objective of SC-LQCD is to study thermodynamic properties. Since β = 0, we cannot vary the temperature T = 1/(N τ a(β)) continuously via the lattice spacing, but only with the lattice extent N τ . The chiral transition is however at temperatures much higher than 1/2, such that for temperatures 1/N τ we are always in the chirally broken phase. The solution is to introduce an anisotropy in the Dirac Operator to favor fermion propagation in temporal direction. In contrast to the chemical potential, the bare anisotropy does not distinguish between forward and backward temporal direction. The temperature on an anisotropic lattice is given by the inverse of the lattice extend in temporal direction but the functional dependence ξ(γ) of the ratio of the spatial and temporal lattice spacings on the bare anisotropy is not known a priori. Hence also the dependence of T on γ is unknown. The main motivation for this study is to overcome this difficulty. The weak coupling analysis of Eq. (1) suggests that ξ(γ) = γ, but this does not carry over to strong coupling, where quarks are confined on links to color singlets. In the mean-field approximation of SC-LQCD [24] based on 1/d-expansion (with d the spatial dimension) the critical temperature is given by suggesting that a σ T c ∝ γ 2 c Nτ is the sensible N τindependent identification in leading and next-to leading order in d.
It is however possible to determine the function ξ(γ) non-perturbatively on anisotropic lattices with by a bare anisotropy calibration γ 0 (ξ) via conserved currents in both spatial and temporal direction [31]. For large N τ (implying large ξ and γ), it turns out numerically that ξ diverges as The precise value of κ can only be determined nonperturbatively and a posteriori, based on the values γ 0 (ξ) measured via the anisotropy calibration (see also Sec. IV A) and has been extrapolated for SU(3) from the set ξ = {0.5, 1, 1.5, 2, 3, 4, 5, 6, 8} to ξ → ∞. The function ξ/γ 2 , which we call the anisotropy correction factor, can be either parameterized by ξ or γ and is well described by the Ansätze Clearly, the extrapolated value for κ based on ξ → ∞ will depend on the Ansatz, as shown in Fig. 1. The Taylor expansion in 1/ξ 2 ∼ a τ 2 in Eq. (13) is limited to the fit range ξ ≥ 2 and results in the value κ = 0.7824(1), which is consistent with the already determined value in [31]. The second Ansatz Eq. (14) has only κ as a free parameter, and interpolates the data for all ξ surprising well, although there are deviations. By construction, ξ/γ 2 = 1 for γ = 1. In order to improve on this one-parameter fit, the third fit Ansatz Eq. (15) introduces 3 additional independent fit parameters to connect the regime ξ > 1 with the opposite regime ξ < 1; with Q 1 > 0 and Q 2 < 0: The fit parameters A and B are thus not independent. This fit results in a non-monotonic behavior, which reflects the fact that ξ γ 2 ξ=8 = 0.7834(2) is larger than ξ γ 2 ξ=6 = 0.7828 (2). Also, it has the smallest reduced chi-squared. Thus we think that the extrapolated result κ = 0.8017(2) is more trustworthy. The error is purely statistical, and the systematic error due to the choice of the fit ansatz is unknown. In Sec. IV A) we will overcome the ambiguities of the extrapolation a τ → 0 by measuring κ directly in the continuous time limit.
We will see in Sec. III E that many observables and the phase diagram have a strong N τ -dependence, which can even be non-monotonic. This requires large N τ to have control over the extrapolation. Hence we want to eliminate γ and N τ all together from the partition function Eq. (5) and replace them by the temperature aT . The continuous time definition of the temperature in lattice units is where we have dropped the subscript, a ≡ a σ . The limit N τ → ∞, γ → ∞ is a joint limit and the second condition implies that γ diverges as γ = √ aT N τ for N τ → ∞. Likewise we can define unambiguously the continuous time chemical potential to replace the chemical potential a τ µ q in Eq. (5): which is also consistent with the γ-dependence of the mean-field critical chemical potential µ c (T = 0) obtained via 1/d-expansion [23], similar to Eq. (10): Now all discretization errors from finite a τ are removed.
The new partition function in continuous Euclidean time will be derived in the next section. We will then have to check numerically that the above limits are well defined for the typical observables. We will present a worm-type Monte Carlo algorithm which samples the partition function efficiently. We denote aT , aµ B as the bare temperature and bare chemical potential which are then renormalized by κ. We will see that we can determine κ nonperturbatively directly by Monte Carlo simulations in the continuous time limit.

C. Continuous Time Partition Function
We will now explain in detail how to derive the continuous time partition function from the discrete time partition function Eq. (5) by tracing the γ-dependence and neglecting subleading terms that vanish in the limit N τ → ∞. The first step to obtain these results is to factorize Eq. (5) into the temporal and spatial part: where the factor N c ! from the site weights for zero monomer number cancels the 1/N c ! in the temporal gauge link, and a prefactor γ NcΛ was pulled out such that spatial links are now suppressed by 1/γ 2 for mesons and 1/γ Nc for baryons. Also we have put the Grassmann constraint Eq. (6) into the above equation via a Kronecker delta and the decomposition Eq. (7). We will now neglect the sub-leading terms, i.e. we will only keep terms that survive in the limit Eq. (17). For any temperature, the average contribution per time location is 1/γ 2 . This will have drastic consequences, as spatial baryons for N c ≥ 3 and spatial dimer occupation numbers k i > 1 will vanish. We will later see how to interpret this outcome and also show numerically that this is well justified. For now we note that the average dimer density will depend on the temperature, and (anti-) baryons are static for N c ≤ 3 for all temperatures and chemical potential. For large γ, N τ the partition function becomes where we have dropped the overall prefactor γ NcΛ and we have used We have resummed static baryons and anti-baryons ω = ±1 in the second line, with |Λ B σ | the number of spatial sites occupied by (anti-) baryons with The sum over configurations contains all possible partitions of the spatial lattice The vertex weights introduced in the second line v ( x,τ ) depend on the dimers k − 0 = k 0 ( x, τ −1) and k + 0 = k 0 ( x, τ ) k − 0 = k 0 ( x, τ − 1) and simplify due to the Grassmann constraint, They come in pairs of adjacent spatial sites ( x, τ ), ( x +î, τ ) at both ends of a spatial dimer on a bond b = ( x, τ ; i). There are N c types of vertices since k − 0 can take the values from 0 to N c − 1. For N c = 3, there are only three types of vertices, The important observation is that the mesonic part of the partition function only depends on the number of vertices, and not on the precise temporal position. The temporal intervals between vertices attached to spatial dimers have a trivial weight: due to the Grassmann constraint on every site where no spatial dimer is attached, k − 0 + k + 0 = N c implies that the dimer numbers form alternating chains as shown in Fig. 2 and cancel in weight: Only the relative order of the vertices v is important, but not the length of the intervals between them. The par-tition function of SC-LQCD with N c = 3 can be written in terms of these vertices as follows: The temporal dimers in the first time slice k 0 ( x, τ = 0) are now dynamic variables in the partition sum. The -vertices and -vertices at sites , and the order in the high temperature expansion is given by the number spatial dimers: In the partition sum, not all temporal positions of the vertices are possible due to the Grassmann constraint (GC). We still need to replace γ by the temperature aT , which requires bookkeeping of possible locations for spatial dimers. We will provide the details in the appendix Sec. IX A. A simplified argument that allows to understand the temperature dependence is that for the first spatial dimer there are up to N τ possible locations between two adjacent spatial sites x, y , but due to the even-odd decomposition there are only N τ /2 possible locations for the second spatial dimer, and likewise for all other dimers, as long as N τ is large. Hence every spatial dimer, after summing over possible locations, has weight where Γ n = {n ( x, τ ), n ( x, τ )} is the set of all valid configurations on the mesonic sublattice Λ M σ with n ≡ N Ds spatial dimers and N ≤ 2n is the total number of -vertices, integrated over the compact temporal direction. Since v = 1, we do not need to include them into the weight. The prefactor 1/n! is due to time-ordering. In the next section we will simplify this result further by a Hamiltonian formulation, where we obtain a meaningful expression for Γ n We now want to discuss the interpretation of the final partition function: as illustrated in Fig. 5, as the temporal lattice spacing a τ a/ξ(γ) → 0, multiple spatial dimers become resolved into single dimers. The overall number of spatial dimers remains finite in the CT-limit, as the sum over O(γ 2 ) sites compensates the 1/γ 2 from spatial dimers. Its number is a function of the temperature and will signal spontaneous chiral symmetry breaking, see Sec. V A. As shown in Fig. 4 it takes large N τ such that double dimers vanish, but it does not require large N τ to make baryons static. The sign problem has completely vanished as σ( ) = 1 for static baryons loops . The set of all baryonic sites coincides then with the fermion bags that have been discussed in [32]. The ex-pansion in n is an all order high temperature expansion. It will also hold at very low temperatures and we will be able to address zero-temperature phenomena.

D. Hamiltonian Formulation
In order to rewrite the partition function further, we make use of a diagrammatic expansion. These methods, giving rise to Quantum Monte Carlo, are nowadays widely used in condensed matter [33,34]. The general idea is to decompose the Hamiltonian H = H 0 + H i and express the partition function in terms of an expansion parameter n which keeps track of the number of interactions described by H i . After summing over all configurations of a given order in n, one integrates over all possible times at which interaction events may take place We will take a step back and reformulate Eq. (21) in new degrees of freedom: the temporal dimers k 0 (x) are replaced by an occupation number m(x) by the following assignment: with (x) = ±1 the parity of a site introduced in Eq. (8).
As a consequence, the alternating dimer chains will be replaced by meson occupation numbers m(x) which is constant on the interval between attached spatial dimers (see Fig. 2), and the dimer-based vertices v(0|2), v(1|1), v(0|2) in Eq. (25) are replaced by occupation numberbased verticesṽ(m|m ), which change the meson state by one unit: m(x) → m (x) = m(x) ± 1: In fact there is a conservation law: if a quantum number m(x) is raised/lowered by a spatial dimer, then at the site connected by the spatial dimer the quantum number is lowered/raised. This is a direct consequence of its def-inition Eq. (30): the parity of the two sites connected by a spatial dimer is opposite. We therefore can replace the vertices by raising and lowering operators: This result is valid for N c = 3, N f = 1. A corresponding result for N f = 2 is given in the appendix Sec. IX D.
In the second line we have included the baryonic sites into the trace, and introduced the mesonic raising and lowering operatorsĴ + ,Ĵ − (which contain the vertices), and the baryon number operatorN . The block-diagonal structure expresses the fact that the Hilbert space of hadrons is a direct sum of mesonic states and baryonic states, |h = |m ⊕ |b , which results in the vanishing commutator The fact that mesons and baryons are mutually exclusive (leading to the factorization into mesonic and baryonic subvolumes) results in one of the two blocks being zero in both operatorsĤ andN . The meson states |m count pseudoscalars, and we will denote them as pions π (despite the fact they are flavorless for N f = 1 and they cannot be distinguished from the η or η mesons). The pion current is conserved, but only in the chiral limit. Monomers would generate a mass to the pion. Since Pauli saturation holds on the level of the quarks and pions have a fermionic substructure, we cannot have more than N c pions per spatial site. Due to the conservation of the pion current, if we start on each site with N c pions, or with no pions at all, there cannot be any spatial dimer that transfers a meson to an adjacent site: either all sites are already saturated with mesons, or there is no meson to be transferred. If we omitted the additive constant N c /2 from Eq. (30), particle-hole symmetry becomes evident. To see this, consider the anti-commutator of the mesonic operators (restricted on mesonic states): The corresponding algebra has the structure of a spin, and it generalizes via Eq. (24) to arbitrary N c : is expressed in terms of the Pauli matrices, and the continuous time partition function becomes that of the quantum XY model. Although the algebra resembles that of a particle with spin, it has nothing to do with the spin of mesons or quarks. The alternating chains are simply expressing the fact that for staggered fermions, the lattice spacing is 2a τ rather than a τ . By shifting the pion occupation numbers by its average value, we can identify the quantum state corresponding to this algebra: This remarkable result is due to the fact that pion occupation numbers on the lattice are not just bounded from below but also from above. We conclude this section by providing a physical interpretation of the dynamics on the hadronic states: the pion dynamics encoded in the Hamiltonian is that of relativistic pion gas [35]. In contrast, the fact that baryon becomes static is due to its non-relativistic nature. Its restmass is large but finite (see Sec. IV C).

A. Poisson Process
Before we address the algorithm that samples the partition function Eqs. (29,32), we want to emphasize an important property: spatial dimers are distributed uniformly in time. The interval length (interpreted as the inter-arrival time between spatial dimers) are then exponentially distributed and the number of spatial dimers in a fixed time interval is Poisson distributed. Hence they can be generated via a Poisson process: with λ the "decay constant" for spatial dimer emissions. Due to the presence of baryons, λ is space dependent: The Poisson process of emitting pions from ( x, t) to an adjacent site ( y, t) with probability λ gives rise to a decomposition of vertices into emission sites ( x, t) ∈ E and absorption sites ( y, t) ∈ A. Spatial dimers can be oriented consistently due to the underlying even/odd decomposition of lattice sites, but is also evident in the Hamiltonian representation, where J − is an emission and J + is an absorption event. The emission sites E are simply those that reduce the pion occupation number m in ✏(2) = +1 < l a t e x i t s h a 1 _ b a s e 6 4 = " L Z l 6 l + U X L k Q m o p E K v d r Y r w 2 Y e r g = " > A A A C 1 3 i c h V F L S 8 N A E J 7 G V 1 t f V Y 9 e g k W o C C W J g l 6 E g g + 8 C B X s Q 9 o q m 3 S t S / M i S Q u 1 i D f x 6 s 2 r / i v 9 L R 7 8 s q a C F u m G z c x + 8 8 2 3 M z u m b 4 s w 0 r T 3 l D I 1 P T M 7 l 8 5 k 5 x c W l 5 Z z K 6 v V 0 O s F F q 9 Y n u 0 F d Z O F 3 B Y u r 0 Q i s n n d D z h z T J v X z O 5 h H K / 1 e R A K z 7 2 I B j 5 v O a z j i h t h s Q j Q V Z P 7 o b A 9 t 2 B s H W z r 1 7 m 8 V t T k U s c d P X H y l K y y l / u g J r X J I 4 t 6 5 B A n l y L 4 N j E K 8 T V I J 4 1 8 Y C 0 a A g v g C R n n d E 9 Z 5 P b A 4 m A w o F 3 8 O z g 1 E t T F O d Y M Z b a F W 2 z s A J k q b W K f S E U T 7 P h W D j + E / c S + k 1 j n 3 x u G U j m u c A B r Q j E j F c + A R 3 Q L x q R M J 2 G O a p m c G X c V 0 Q 3 t y 2 4 E 6 v M l E v d p / e g c I R I A 6 8 q I S s e S 2 Y G G K c 9 9 v I A L W 0 E F 8 S u P F F T Z c R u W S c u l i p s o M u g F s P H r o x 6 M W f 8 7 1 H G n a h T 1 n a J x v p s v G c n A 0 7 R O G 1 T A V P e o R K d U R h 0 W l F / o l d 6 U S + V B e V S e v q l K K s l Z o 1 9 L e f 4 C o Y C T e Q = = < / l a t e x i t > ✏(0) = +1 < l a t e x i t s h a 1 _ b a s e 6 4 = " 7 z H O 6 N w u Q c D 8 f z B a d A T z 9 P 7 L X r o = " > A A A C 1 3 i c h V F L S 8 N A E J 7 G V 1 t f V Y 9 e g k W o C C W p g l 6 E g g + 8 C B X s Q 9 o q m 3 S N S / M i S Q u 1 i D f x 6 s 2 r / i v 9 L R 7 8 s q a C F u m G z c x + 8 8 2 3 M z u G b 4 s w 0 r T 3 l D I 1 P T M 7 l 8 5 k 5 x c W l 5 Z z K 6 u 1 0 O s F J q + a n u 0 F D Y O F 3 B Y u r 0 Y i s n n D D z h z D J v X j e 5 h H K / 3 e R A K z 7 2 I B j 5 v O 8 x y x Y 0 w W Q T o q s X 9 U N i e W 9 C 2 D r b 1 6 1 x e K 2 p y q e O O n j h 5 S l b F y 3 1 Q i z r k k U k 9 c o i T S x F 8 m x i F + J q k k 0 Y + s D Y N g Q X w h I x z u q c s c n t g c T A Y 0 C 7 + F k 7 N B H V x j j V D m W 3 i F h s 7 Q K Z K m 9 g n U t E A O 7 6 V w w 9 h P 7 H v J G b 9 e 8 N Q K s c V D m A N K G a k 4 h n w i G 7 B m J T p J M x R L Z M z 4 6 4 i u q F 9 2 Y 1 A f b 5 E 4 j 7 N H 5 0 j R A J g X R l R 6 V g y L W g Y 8 t z H C 7 i w V V Q Q v / J I Q Z U d d 2 C Z t F y q u I k i g 1 4 A G 7 8 + 6 s G Y 9 b 9 D H X d q p a K + U y y d 7 + b L p W T g a V q n D S p g q n t U p l O q o A 4 T y i / 0 S m / K p f K g P C p P 3 1 Q l l e S s 0 a + l P H 8 B n K y T d w = = < / l a t e x i t > ✏(1) = 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " N f G 3 a E X O k j U U n o W u a E Q o C + f / x I E = " > A A A C 1 3 i c h V F L S 8 N A E J 7 G V 1 t f V Y 9 e g k W o B 0 t S B b 0 I B R 9 4 E S r Y h 7 R V N u k a l + Z F k h Z q E W / i 1 Z t X / V f 6 W z z 4 Z U 0 F L d I N m 5 n 9 5 p t v Z 3 Y M 3 x Z h p G n v K W V q e m Z 2 L p 3 J z i 8 s L i 3 n V l Z r o d c L T F 4 1 P d s L G g Y L u S 1 c X o 1 E Z P O G H 3 D m G D a v G 9 3 D O F 7 v 8 y A U n n s R D X z e d p j l i h t h s g j Q V Y v 7 o b A 9 t 6 B v H W z r 1 7 m 8 V t T k U s c d P X H y l K y K l / u g F n X I I 5 N 6 5 B A n l y L 4 N j E K 8 T V J J 4 1 8 Y G 0 a A g v g C R n n d E 9

Z 5 P b A 4 m A w o F 3 8 L Z y a C e r i H G u G M t v E L T Z 2 g E y V N r F P p K I B d n w r h x / C f m L f S c z 6 9 4 a h V I 4 r H M A a U M x I x T P g E d 2 C M S n T S Z i j W i Z n x l 1 F d E P 7 s h u B + n y J x H 2 a P z p H i A T A u j K i 0 r F k W t A w 5 L m P F 3 B h q 6 g g f u W R g i o 7 7 s A y a b l U c R N F B r 0
t F y q u I k i g 1 4 A G 7 8 + 6 s G Y 9 b 9 D H X d q p a K + U y y d 7 + b L p W T g a V q n D S p g q n t U p l O q o A 4 T y i / 0 S m / K p f K g P C p P 3 1 Q l l e S s 0 a + l P H 8 B n K y T d w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " K f s c F z m 0 W + T Q z N 0 c W Q 6 p X 5 e H 5 5 Q = " > A A A C 0 X i c h V F L S 8 N A E J 7 G V 1 t f V Y 9 e g k X w V J I q 6 L H g A w 8 K F e 0 D 2 l I 2 6 R q X 5 k W S F m o R x K s 3 r / r H 9 L d 4 8 N s 1 F b R I N 2 x m 9 p t v v p 3 Z s U J X x I l h v G e 0 u f m F x a V s L r + 8 s r q 2 X t j Y r M f B I L J 5 z Q 7 c I G p a L O a u 8 e u r i E 6 n i u l A w 1 L n I V 7 A h 6 2 h A v n K E w V d d d y D Z c p y p e K n i g x 6 E a x 8 f d S D M Z t / h z r t 1 M s l c 7 9 U v j o o V s r p w L O 0 T T u 0 h 6 k e U o X O q Y o 6 5 D R f 6 J X e t G t t p D 1 q T 9 9 U L Z P m b N G v p T 1 / A W i a k e M = < / l a t e x i t >

FIG. 2.
Correspondence between discrete time configurations in terms of dimer coverings and baryon world lines (top) and in terms of hadron occupation numbers in continuous time (bottom). Multiple spatial dimers become resolved in single spatial dimers (which can be oriented consistently from emission sites E to absorption sites A, indicated by the arrow), baryons become static and only vertices of -shape or -shape survive as aτ → 0.
Euclidean time by one unit, the absorption sites A are those that increase m. Every spatial dimer corresponds to a pion hopping to an adjacent site and connects an E-site with an A-site. The number of E-sites equals the number of A-sites, and due to the periodic boundary conditions in time this even holds for every spatial site x: The continuous time worm algorithm (CT-WA) needs to fulfill detailed balance, such that the emission process is counterbalanced by an absorption process to obtain the equilibrium distribution of spatial dimers according to temperature and chemical potential.

B. Details of the Continuous Time Worm Algorithm
Worm algorithm at discrete time (DT-WA) are well established also for strong coupling lattice QCD [29,36,37]. Designing an algorithm that operates directly in the Euclidean continuous time limit has several advantages: (1) the ambiguities arising from the functional dependence of observables on the anisotropy ξ(γ) -in particular nonmonotonic behavior -will be circumvented and (2) we do not need to perform the continuum extrapolation N τ → ∞. This will allow us (3) to measure the phase boundaries unambiguously, as the baryonic part of the partition function simplifies such that the sign problem is completely absent, and (4) at all temperatures of interest, the CT algorithm is considerably faster than its discrete version, in particular since the baryon update does not require a worm algorithm but can be replaced by a heat bath update.
In Fig. 3 we clearly see that the CT-worm algorithm outperforms the DT-worm algorithm at temperatures in the vicinity of the transition temperature and above. The lower the temperature, the more spatial dimers are sampled, which makes the average worm update longer. At first glance it seems that the CT-worm becomes more expensive, but one needs to keep in mind that lower temperatures require larger N τ to get valid estimates for observables. On a lattice with time extent N τ , temperatures below 1/N τ (which have γ < 1) will have more spatial dimers than temporal dimers and suffer from saturation effects: The density of spatial dimers is limited to N c N τ /2, whereas it is unlimited at continuous time.
In Fig. 4 we show the N τ -dependence of various observables: they have a well-defined CT-limit. Also, this figure illustrates that the approximations which lead to Z CT in Eq. (32) are well justified. The extrapolation from discrete time to continuous time is difficult: large N τ require more statistics, and due to the sign problem, most observables get noisy due to sign reweighting. The first approximation is to make baryons static, which eliminates the sign problem. This step makes the extrapolation much more controlled, and even for N τ = 4, the static baryon approximation is not bad. Next we prohibit sites which have more than 3 spatial dimers, which has only a mild effect at the temperatures considered here. If we also prohibit sites with more than 2 spatial dimers, the deviation at finite N τ is drastic, but also this approximation extrapolates to the same CT-limit for the observable. The point at 1/N τ = 0 in Fig. 4 is the outcome of the CT-WA, which has much smaller error bars and better performance with the same number of worm updates.
Continuous time (CT) algorithms as for Quantum Monte Carlo are now widely used in condensed matter (see e.g. [33,38]), whereas CT methods in quantum field theories is rather new [30,[39][40][41] The basic idea of a worm algorithm introduced in [42] is to sample an enlarged configuration space with defects on the lattice known as worm tail x T and worm head x H . Every worm algorithm consists of two kinds of updates: (1) move updates, which move the head x H and tail x T to a new site x 0 , and (2) shift updates, which move the head x H through the lattice until the worm recombines with the tail. Worm algorithms are highly efficient: after recombination, the configuration has been globally updated, similar to cluster algorithms. Moreover, during the shift update, 2-point correlation functions can be measured. In order to apply a worm algorithm, the partition function needs to be written in terms of bond variables. Those representations are typically available in spin models from the high-temperature expansion. In the case of lattice QCD, a dual representation based on the strong coupling expansion also admits the applicability of worm algorithms.  Our Continuous Time Worm Algorithm (CT-WA) can be derived from the Discrete Time Directed Path Worm Algorithm (DT-WA) that has been developed for U(N c ) gauge group in the strong coupling limit [36], which does not include baryons. This worm algorithm is based on an even-odd decomposition of weights: if the parity of the head (x H ) is the same as that of the tail (x T ), then the head has an active site location, if the parities differ the head is an a passive site. The active sites correspond to the absorption sites A, and the passive sites correspond to the emission sites E as discussed above.
For SU(N c ) gauge group, it is required to have two separate worms, one in the mesonic sector and one in the baryonic sector [37]. The mesonic worm for the SU(N c ) group differs from the directed path Worm for U(N c ) in one important aspect: In the directed path version backtracking is prohibited to evolve faster through configuration space (if the update shifts the worm head from x to the adjacent site y, then in the next shift update the worm is not allowed to go back). With the simple baryon loop geometries in the CT-limit, we can supplement the continuous time version of the directed path worm algorithm by an additional heat bath update: after the mesonic worm has recombined, we propose for all sites x where no spatial dimers are attached (the so-called static sites) a new hadronic state with the probabilities The consequence is that if the worm head propagates in positive/negative temporal direction, it will continue to do so until it will either emit or absorb a pion, i. e. it will either add or delete a spatial dimer. It will not change the direction and diffuse: the CT-WA can be regarded as a Poisson process. The updating rules are outlined in Fig. 5. The probabilities for the various cases (approaching/leaving an absorption site A or emission site E) depend on the involved states m: (1) if an A-site is approached from the temporal direction, the spatial dimer is removed with a heat bath probability determined by J − , (2) if an A-site is approached from a spatial direction, the new temporal direction is also determined by J − , (3) if an E-site is approached from temporal direction, the emission probability to insert a spatial dimer is 1 − e −λ∆τ and the probability to continue in temporal direction is e −λ∆τ , (4) if it is approached from spatial direction, the new temporal direction is chosen equally likely. At high temperatures, λ = λ(aT ) 1 according to Eq. (38) and the worm head will very likely continue in temporal direction by some time ∆t with probability p τ 1 − λ∆t and emit a spatial dimer with probability p σ λ∆τ . The higher the temperature, the longer the worm propagates in temporal direction, possibly looping through the periodic boundary back to where it started.
In the discrete time algorithm, during worm evolution, whenever the worm head is on a site with opposite parity compared to the worm tail, (x H ) = − (x T ), both worm head and tail can be interpreted as monomers (if (x H ) = (x T ), the head is a sink rather a source for monomers). Even in the chiral limit, the monomer 2-point function can be accumulated in a histogram (due to translation invariance, only the relative lattice vector z = x 1 − x 2 is needed): with d(x) defined in Eq. 38. An equivalent definition holds in the CT-limit:  Nτ -dependence of the chiral susceptibility (top) and the energy (center ) and the baryon susceptibility (bottom). We compare the full discrete simulations and various approximations according to the steps in deriving the continuous time limit (static baryon approximation, exclusion of spatial triple dimers, , exclusion of spatial double dimers). We have fixed the bare temperature to aT = 1.2 < aTc and aT = 1.5 > aTc All observables extrapolate well into the continuum limit, with its Monte Carlo result at 1/Nτ = 0 having much smaller error. Top: an absorption site can be approached either from the temporal direction (left: a spatial dimer may be removed) or from the spatial direction (right: a dimer was emitted in the previous step). Bottom: an emission site can be approached either from the temporal direction (left: a spatial dimer may be emitted) or from the spatial direction (right: a dimer was removed in the previous step).

C. Observables
Almost all observables that can be measured via the DT-WA version can also be measured via CT-WA. This is obviously the case for all observables that can be obtained as derivatives of log Z CT . The discrete time observables in terms of the dual variables are discussed in [43]. The corresponding dimensionless thermodynamic observables in the CT-limit simplify because which should be compared to the isotropic case based on Eq. (14): Also, in the CT-limit we have no longer temporal dimers but only spatial dimers, and have to consider the chiral limit: We are now able to define the continuous time observables in terms of dual variables, which are always in dimensionless units with a = a σ and V = N σ 3 a 3 . Important observables are (1) the baryon density: which is given by the average winding number; (2) the energy density where the irrelevant additive constant C = 1 2 N c Λ can be neglected compared to discrete time as we dropped the prefactor γ NcΛ in Eq. (21) which contained both the contribution from static mesons and static baryons; (3) the pressure which in the strong coupling limit and chiral limit is just proportional to the energy density such that the interaction measure − 3p vanishes. At finite quark mass, the interaction measure is proportional to the chiral condensate, which here is zero in a finite volume as χχ ∝ n M (but see Sec. IV); (4) the chiral susceptibility with only has the connected contribution non-zero in the chiral limit and G 2 ( is the translation invariant monomer 2-point function that is measured during worm evolution, see Eq. (42); (5) the entropy density The chiral condensate vanishes in the chiral limit in a finite volume. This is also evident from the absence of monomers in the dual representation. It is possible to obtain the chiral condensate from a 1/V expansion via chiral perturbation theory in a finite box, as explained in Sec. IV. Note that the pressure defined in Eq. (49) is not equal to because on the lattice the system is not homogeneous. The identity p = p only strictly holds in the continuum.

D. Polymer Formulation and Wang Landau Method
So far we have treated the mesonic and baryonic sector separately, and there is no need for the resummation known as the Karsch-Mütter trick [28] for real chemical potential as there is no sign problem in the CT-limit. However, a resummation of static mesons and baryons proves to be advantageous in several respects: (1) it allows to extend simulations to imaginary chemical potential beyond the value of a τ µ q = iπT /6, where the baryon density becomes zero (discussed in Sec. V D), (2) we are able to adapt the Wang-Landau method [44] for determining the first order transition at low temperatures very accurately, and obtain also the canonical phase diagram from the density of states at high precision, see Sec. V.
Apart from the usual (anti-) baryons denoted by B, will discuss here two kinds of resummations of quantum states: the superposition of baryons and anti-baryons (P-Polymers), and including static mesons (Q-Polymers): where for a given configuration C, on each spatial site, the baryon and polymer numbers B ≤ P ≤ Q are related via (in the following V = N σ with m x = 1 iff the site is mesonic and static. The corresponding single site weights are: These weights will be used for the following binomial/trinomial distributions: with B ± = P ±B 2 the number of (anti-) baryon sites and Q − P is the number of static mesons. For some observables we need higher moments of the baryon number. We then only keep track of the histogram for Q-polymers, H Q V,aT ,aµB (Q) (normalized accordingly to be a probability distribution), and get the histogram in the baryon number H B V,aT ,aµB (B) from the above distributions: which improves drastically over the usual measurement of higher moments. In Fig. 6 we show histograms H Q V,aT ,aµB for various temperatures and µ B = 0. The temperature dependence gives insight into the number of static vs. dynamic sites: at high temperatures, almost all sites are static, and at low temperatures almost all sites are dynamic, e.g. they interact via pion exchange with adjacent sites. The critical temperature is characterized by a broad distribution. Another important application of histogram techniques is the Wang-Landau method, which computes the density of states g(aT , B). It will allow us to obtain the canonical phase diagram, see Sec. V. We use that the grand-canonical partition sum is related to the canonical partition sum via the Laplace transformation One method to determine the canonical partition sum Z C (aT , B) in the context of QCD is to obtain the Z GC for imaginary chemical potential and reweighting for the resulting Fourier coefficient [45]. In the dual representation, Z C (aT , B) can be determined directly by the Wang-Landau method, since it is in fact the density of states with respect to the canonical conjugate to aµ B and it is approximated by g(aT , B) up to the target precision. Then, observables in the GC-ensemble are immediately obtained: The accuracy even improves when the density of states using the polymer resummation g(aT , P ) is determined via Wang-Landau, and the canonical partition sum is recovered by the binomial transformation Eq. (56): The Wang Landau method applied to g(aT , P ) consists of the following steps: (1) A CT-worm update is run (which makes g(aT , P ) temperature dependent).
(2) We loop through all spatial sites x and check whether the site is static (has no spatial dimers attached).
(2b) If not, the configuration is unchanged and ∆P = 0.
(3) The new configuration is accepted with a metropolis acceptance step: (3a) If accepted, P = P + ∆P is the new polymer number, (3b) If rejected P = P .
(4) In any case, even if the site was non-static and P = P (option (2b)) the histogram and density of states are updated: H(P ) → H(P ) + 1, log(g(P )) → log(g(P )) + log(f ) with f the modification factor.
We loop through (1-4) until the histogram H(P ) is flat enough: withH the histogram average and δ defining the flatness condition. This step, which refines g(P ), is repeated until the final precision is reached, log(f ) ≤ log(f final ). Then g(P ) approximates the true density of states with that precision. In Sec. (V) we will show the density of states and the canonical phase diagram for various temperatures.
We perform simulations at a set of fixed temperatures and weight the obtained density of states to the critical aµ c , which is characterized by equal probability of the low and high density phase. In practice, we determine aµ c at which both peaks in the first order region have the same height (see Fig. (20)).

E. Crosschecks
To check the correctness of our CT-WA implementation, we have made extensive crosschecks. A comparison of the CT-algorithm on volumes with an analytic result extrapolated from 2 × N τ lattice for gauge group U(1) is discussed in the appendix Sec. IX B. Since there does not seem a simple analytic expressions for N c > 1, we are left with comparing continuous time simulations with the extrapolation of discrete time simulations. We already discussed the suppression mechanism that lead to the continuous time results for various observables in Fig. 4. In Fig. 7 we show a comparison of the discrete time extrapolation and the continuous time simulations for the chiral susceptibility as a function of the temperature, which agree within errors for all temperatures. Another aspect is to verify that the distribution of spatial dimers is indeed Poissonian, due to the fact that the weight of a configuration does not depend on the interval lengths between subsequent spatial dimers. This is illustrated in Fig. 8. The Poisson distribution P (N (∆τ ) = n) = (λτ ) n n! e −λ∆τ (66) has been fitted to histograms from Monte Carlo via CT-WA. The comparison with the expected values of λ (with λ = 3 4T for the distribution of spatial dimers per bond and λ = 6d 4T for the distribution of vertices per sites, with d = 3) is very good for small intervals ∆τ < 1. The deviations to the expected λ for large intervals ∆τ 1 is due to the periodic boundary conditions, where the Poisson distributions start to overlap.

A. Determination of κ and Pion Decay Constant
The first task that is also relevant to define the temperature and chemical potential non-perturbatively (Eqs. (17,18)) is to determine the anisotropy correction factor κ, see Eq. (12). The procedure of anisotropy calibration is discussed for anisotropic lattices at strong coupling in discrete time in [31,46,47] in detail. The coefficient κ is the strong coupling analogue of the Karsch coefficients at weak coupling that have been analyzed in [28,48] and numerically studied at fixed physical scale in [49]. Anisotropic lattices are also relevant when determining mesonic correlators, e. g. in the FASTSUM collaboration [50].
Our strategy to obtain κ is based on the variance of the pion current. In the chiral limit, the pion current for discrete time is a conserved current: Likewise, also the corresponding currents in the CT-limit see Eq. (30) -where we have dropped the baryonic contributions and the constant, as they do not contribute at continuous time -are conserved and now directly linked to the meson occupation numbers: for all τ 2 > τ 1 , and the temporal/spatial charges are which have the expectation values The variances are however temperature dependent. If the spatial and temporal variances are equal, that corresponds to equal physical extent in space and time: This allows us to measure κ: given the lattice extent N σ , we scan the bare temperature aT to determine its value aT 0 that corresponds to a physically isotropic lattice: This calibration is shown in Fig. 9, the results for κ for various volumes are shown in Table I and its extrapolation in Fig. 10 (left). The finite size effects are very small. Note that in contrast to the previous study [31], there is no reason to distinguish κ for gauge group U(3) and SU(3): the thermodynamic extrapolation N σ → ∞ coincides with the zero temperature extrapolation, and since the calibration is performed at aµ B = 0, static baryons are virtually absent (see also Fig. 6). This is not the case at finite ξ (finite a τ ). As discussed in Sec. II B, the determination of κ in [31] suffers from systematic uncertainties as the extrapolation in ξ is based on rather small ξ ≤ 8. Our final continuous time result κ = 0.797(1) is consistent with the extrapolations, favoring Ansatz 3. In Fig. 10 (right) we show the thermodynamic extrapolation of the helicity modulus, which yields the square of the pion decay constant: resulting in aF π = 0.7797 (1). This compares well with the extrapolation of discrete time [31] which yields aF π = 0.78171(4), taking into account that the extrapolation of a 2 F 2 π has similar uncertainties as κ, which are overcome by the continuous time simulations. Right: Thermodynamic extrapolation of the helicity modulus a 2 Υ, from which we extract pion decay constant at zero temperature.
The method of anisotropy calibration has also been extended by us to finite quark mass [43] and recently also to finite β. These results are a clear indication that it is possible to define the continuous time limit unambiguously for finite m q and finite β in the strong coupling regime, with κ = κ(m q , β).

B. Chiral Condensate and Chiral Susceptibility
Despite the fact that in the chiral limit, the chiral condensate is zero in a finite volume -in the dual representation this is due to the absence of monomers -it is nevertheless possible to extract the chiral condensate from the chiral susceptibility χ σ (which is non-zero in a finite volume). The corresponding chiral perturbation theory in a finite box -the so-called -regime -is an expansion in the inverse volume [51], and for the O(2) model in d = 4: where β 1 = 0.140461 and β 2 = −0.020305 are shape coefficients of a finite 4-dim. box. Note that the value Σ that can be extracted from this equation corresponds to the chiral condensate in the thermodynamic limit. In Fig. 11 we show the fit according to Ansatz Eq. (78) to obtain the chiral condensate from the Monte Carlo data of the chiral susceptibility for various volumes, all in the CT limit. Apart from Σ, we also treat α as a fit parameter as we do not know the values of the renormalization group invariant scales Λ Σ and Λ M , but it turns out that α is consistent with zero. The value of aF π determined in the previous section is used. The thermodynamic extrapolation N σ → ∞ coincides with the zero temperature extrapolation as the bare temperature is set to aT = (κN σ ) −1 to always obtain a physically isotropic lattice. Our result from continuous time simulations yields a 3 Σ = 1.305(3) and agrees well with the extrapolation of the Monte Carlo data at discrete time as discussed in [31].

C. Energy and Baryon Mass
The baryon mass m B is an important quantity to understand the nature of nuclear interaction, and its value in lattice units am B is also a good choice to scale other quantities to dimensionless ratios, such as T /m B , µ B /m B . At zero temperature, where the free energy F = E − T S coincides with the internal energy E, the static baryon mass in the strong coupling limit is given by the probability of a baryon to propagate in temporal direction. This can be immediately expressed by the probability of having a static baryon in the ensemble: The extrapolation of the static baryon mass towards continuous time has been discussed in [31] with the result am B = ξa τ m B = 3.556(6) = κ am MF B , am MF B = 4.553 (7), which is about 20% larger than the isotropic value am B = 2.877(2). Since p B 1, the mass is evaluated via the so-called snake algorithm at discrete time: The ratio Z k+2 Z k is the probability to extend a static baryon segment of length k by two segments, and the sum results in a static baryon of length N τ . The method unfortunately does not extent straight forwardly to continuous time: the ratios Z k+2 Z k cannot be measured, since at the end of a static baryon segment there is a finite probability that two spatial dimers are attached at the same location, in contrast to other observables discussed above (Fig. 4). However, we are able to determine the baryon mass from the energy difference based on Eq. (48): The energy density at zero temperature in the CT-limit, if one does not take the irrelevant constant C in Eq. (48) into account (rendering it negative), can be measured at very high accuracy: where the value for gauge group U(3) (which does not have baryons) coincides with the value for gauge group SU(3) (where baryons become suppressed with decreasing temperature). The fact that a 4 0 = − lim T →0 aT n Ds is finite implies that the number of spatial dimers diverges as ∝ 1/T . Note that a previous determination of 0 at discrete time [35] includes the diverging constant: a 4 0 = 0.66(2)ξ. We measured the energy density without ( 0 ) and with a static baryon ( B ), both on discrete and continuous time lattices. The discrete time measurements of a∆E are extrapolated via a polynomial Ansatz in 1/ξ, as shown in Fig. 12. The fit results are summarized in Tab. II, and are compared with the continuous time results. Indeed, we find very good agreement of all extrapolated estimates of the baryon mass with its continuous time result within errors. It should be pointed out that at γ = 1, where k 0 = Nc 2d = 3 8 , the static baryon mass from ∆F (via the snake algorithm) differs substantially from the baryon mass obtained from ∆E. But towards the CT limit, both definitions agree. The extrapolation of the discrete time data (obtained from ∆E or ∆F ) is in 1/ξ rather than 1/ξ 2 : it is more suitable as the extrapolation appears to be almost linear in 1/ξ, but clearly there are additional uncertainties related to the derivative dξ/dγ that are bypassed by simulations directly in the CT limit.
We distinguish between U(3) and SU(3) results for the baryon mass: in U(3) gauge theory, there is only the valence baryon and no µ B -dependence of the partition function, whereas SU (3) Fig. 12. ∆E has been evaluated at various temperatures and extrapolated to zero temperature. The result for the snake algorithm valid for SU(3) differs slightly from the value amB = 3.556(6) given in [31] due to the improved extrapolation used here.
fluctuations. At zero temperature, those baryon fluctuations are largely suppressed. Even though U(3) gauge theory has no baryons, there is no obstacle in measuring the baryon mass in U(3) via the response of a valence baryon to the pion bath, resulting in less statistical noise. Our best estimate of the baryon mass is thus the U(3) result in the CT-limit, as it does not suffer from any ambiguities due to extrapolation: This baryon mass receives contributions from a pion cloud surrounding the static point-like baryon.

A. Chiral Transition
In Sec. IV B we have determined the chiral condensate in the chiral limit at zero temperature. In principle this can be extended to finite temperature, and the chiral transition could be determined by the vanishing of the chiral condensate. It suffices in practice to determine the chiral transition from the chiral susceptibility, which is obtained from the worm algorithm to high precision. Also, this method readily extends to finite density: the chiral transition can be easily obtained from finite size scaling of the chiral susceptibility up to the chiral tricritical point (aµ TCP B , aT TCP ). The finite size scaling of the susceptibility in the -regime is illustrated in Fig. 13 for volumes up to 64 3 × CT at µ B = 0. We expect critical behavior in the O(2) universality class in 3 dimensions, resulting the scaling law [52] lim L→∞ χ(L, T c ) ∝ L γ/ν , γ = 1.3177 (5), The result for the transition temperature is We find that also the specific heat is sensitive to the chiral transition: Fig. 14 shows that a weak cusp develops in the vicinity of T c . Although the strong coupling limit is far away from the continuum for realistic quarks, we can nevertheless compare dimensionless ratios T /m B with continuum extrapolated ratios. With m B 938GeV and the pseudo-critical crossover temperature T pc 154MeV [53] we find that the ratio at strong coupling and in the chiral limit is more than twice as large: (9). (87) The comparison improves when a finite quark mass is considered at strong coupling, as the pseudo-critical transition temperature drops rapidly with the mass while the baryon mass is quite insensitive [54]. We note that the continuous time transition temperature for U(3) gauge group and its comparison with the N τ → ∞ extrapolation have been discussed in [39], with aT The determination of aT c at finite chemical potential is straight forward up to the tricritical point. Fig. 15 illustrates the chiral susceptibility χ σ in the full µ B −T plane. The second order chiral phase transition turns into first order for µ B > µ tric B , and the chiral susceptibility -which is ∝ (ψψ) 2 in the chiral limit -behaves as an order parameter and develops a gap. There is no back-bending of the first order transition, in contrast to discrete time (due to saturation of spatial dimers, N Ds ≤ N c Ω/2), which has been discussed in [31]. Similarly, the energy density (T ) − 0 can be measured in the full µ B − T plane, as shown in Fig. 16. For small chemical potential and for temperatures below T c it behaves according to the Stefan-Boltzmann law [55]: which corresponds to an ideal pion gas and has already been discussed at zero chemical potential for discrete time [35]. At zero temperature, the energy density jumps at the first order transition to the finite value − 0 given in Eq. (83), which is the maximal value corresponding to the absence of spatial dimers.  FIG. 14.
The specific heat, which is proportional to the susceptibility of spatial dimers (see Eq. 48). The typical λshape is apparent in the transition region The chiral susceptibility diverges in the chirally broken phase, but is much smaller in the chirally restored phase. Along the first order transition which is strong already for aT < 0.7 and hinders reliable results below aT < 0.3, artificial wiggles appear due to hysteresis of the overlapping low and high density branches.

B. The Nuclear Transition
Strong coupling lattice QCD exhibits not only spontaneous chiral symmetry breaking and its restoration along a second and first order boundary, but also a nuclear liquid gas transition signaled by the baryon density. In order to determine the first order transition line in the phase diagram, we measure the baryon density and its susceptibility, both by direct simulations at finite chemical potential, and by the Wang-Landau method explained in Sec. III D. The baryon density in the µ B − T plane is shown in Fig. 17. The volumes considered are 4 3 × CT, 6 3 × CT and 8 3 × CT at low temperatures and additionally 12 3 × CT, 16 3 × CT in the vicinity of the chiral tricritical point. Simulations at low temperatures across the first-order transition are challenging: for µ B < µ 1st B , the phase is described as an ideal pion gas, for µ B > µ 1st B the phase is that of a baryon crystal (liquid), resulting in a large latent heat. In a Monte Carlo simulation, tunneling between the phases is exponentially suppressed by the volume and hysteresis between the low and high density phase shows up. This difficulty is overcome by the Wang-Landau method: in Fig. 18 we show the logarithmic density of states for P-polymer and baryon number, and in Fig. 19 the density of states are applied to recover the baryon density via Eq. (58). We find that the full first order nuclear transition coincides with the chiral first order transition. The determination of µ 1st B and the boundaries of the mixed phase is illustrated in Fig. 20 for various volumes. The result of the thermodynamic extrapolation according to which is not very different from the discrete time determination aµ 1st B = 1.78(1) valid for isotropic lattices, γ = 1 [29]. Nuclear matter at strong coupling is in fact a quark saturated phase: the baryon density at zero temperature jumps from n B = 0 to the maximal value n B = 1, where every lattice site is occupied by a static baryon. It is no coincidence that chiral symmetry is restored in the nuclear phase: mesons cannot occupy baryonic sites, leaving no room for spontaneous chiral symmetry breaking. Away from the strong coupling limit, where baryons are no longer pointlike and become spread over several lattice spacings, the nuclear phase may have a non-vanishing chiral condensate.
We want to conclude this section by quantifying the interaction strength between baryons. In the CT-limit we find which should be compared to the discrete time (γ = 1) ratio [29] Hence the nuclear interactions are enhanced in the CTlimit.  It also shows a strong first order behavior and at low temperatures becomes insensitive to the chemical potential below µ 1st B , which is known as Silver-Blaze property. The first order line terminates in a critical end-point which coincides with the chiral tricritical point. The baryon density is not sensitive to the second order chiral transition.

C. The SC-LQCD Phase Diagram
We now want to summarize the previous results on the chiral and nuclear transitions and establish the phase boundaries both of the grand-canonical and canonical phase diagram, shown in Fig. 21. In the grand-canonical phase diagram, one can clearly see that the chiral first order phase boundary and the nuclear transition (obtained from the Wang Landau method, see Tab. III) are on top. In the grand-canonical phase diagram, a mixed phase of both nuclear gas and liquid persists. The low density boundary a 3 n (1) B tends to zero, whereas the high density boundary a 3 n (2) B tends to 1. A meaningful density of nuclear matter cannot be assigned at strong coupling.
There are various strategies to locate the chiral tricritical point, which is characterized as the end point of a triple first order line where the existence of three phases cease to coexist (the nuclear phase and two chirally bro-   ken phases for positive and negative quark mass). According to the Gibbs' phase rule, the upper critical dimension is 3, such that the tricritical exponents are analytic: To distinguish tricritical second order behavior from O(2) critical behavior, Eq. (85), large volumes are required. There is a better strategy, based on the fact that the tricritical point coincides with the nuclear critical endpoint FIG. 20. The probability density, obtained from reweighting the density of states to aµ 1st B (Nσ) such that the two maxima are of the same height, for various volumes and at a fixed temperature aT = 0.5. The first maximum denotes the baryon density a 3 n (which can be made plausible via a percolation analysis, see Sec. IX C). This is clearly only expected in the strong coupling limit, but also holds for small values of β at finite N τ [56]. The nuclear end point is characterized by the vanishing of the mixed phase, resulting in n If one does not take into account the rescaling with κ, then aT TCP = 0.98(3) and aµ B TCP = 1.92(6) compares quite well with its determination on a disrete lattice: aT TCP Nτ =4 = 0.94 (7), aµ B

TCP
Nτ =4 = 1.92(9) [37], indicating that the N τ corrections are small up to the critical point and become only large at lower temperatures [31]. We also note that the mean field tricritial point deviates substantially: aT TCP MF = 0.866, aµ B TCP MF = 1.731 [57]. As soon as a small finite mass is introduced, the chiral tricritical point turns into a chiral critical end point of Z(2) universality class. Close to the chiral limit, we estimate which may in principle be within reach with conventional hybrid Monte Carlo, based on the fermion determinant such as Taylor expansion [58]. But with increasing quark mass also the ratio µ CEP B /T CEP increases rapidly (aµ CEP B increases whereas aT CEP decreases), as has been studied for discrete time in [54]. The critical endpoint is quickly The grand canonical phase diagram in the a 3 nB − aT plane. Note that due to Pauli saturation, at zero temperature the mixed phase at zero temperature extends to the full range in a 3 nB out of reach for methods of circumventing the sign problem via HMC methods. In the appendix Sec. IX E we elaborate further on the prospects of finite quark masses in the continuous time limit.
Our new results eliminate systematic uncertainties in previous findings in Monte Carlo for fixed N τ [29].

D. Extension to Imaginary Chemical Potential
Lattice QCD at imaginary chemical potential is usually considered because in contrast to non-zero real chemical potential, the fermion determinant is sign problem-free and it allows to analytically continue to real chemical potential [2]. It is also interesting in its own right due to the Roberge-Weiss periodicity [59] and the Roberge-Weiss transition [60].
In the dual representation of SC-LQCD at discrete time, it is not straightforward to simulate at imaginary chemical potential.
However, at continuous time where baryons are static, we can use cosh(iµ im B /T ) = cos(µ im B /T ), and with the Pand Q polymer resummation (see Sec. III D): The second equation enables us to measure the chiral transition for arbitrary imaginary chemical potential. Our result is shown in Fig. 22. At the Roberge-Weiss point µ im B /T = π we do not find a cusp, in contrast what would be expected at weak coupling. We also cannot observe a first order transition in the chiral observables, which is expected as the partition function becomes analytic in the high temperature limit. By integrating out the gauge links, the center sectors are no longer distinct. Gauge observables such as the Polyakov loop should be able to signal a first order transition between the center sectors at high temperatures, which requires to include gauge correction. We also want to note that the point at µ im B /T = π/2 is special as it corresponds to the U(3) transition temperature aT = 1.8843(1) (as discussed in [39]) as P-polymers have weight w p = 0 according to Eq. (53).

A. Taylor Expansion
The dual representation of SC-LQCD is a great laboratory to benchmark other methods to circumvent the sign problem. One of the prominent methods in the context of lattice QCD is the Taylor expansion [3], which might allow to estimate the location of a possible chiral critical endpoint based on estimates for the radius of convergence of the Taylor series. The standard thermodynamic observable that is Taylor expanded for that purpose is the pressure. This requires high orders of the Taylor series, but the current state of the art is limited to 6. order (improved action) [58] and 8. order (unimproved action) [61]. It turns out that due to the continuous time limit and by taking into account both the polymer resummations and histogram method presented in Sec. III D, we are able to determine higher orders of Taylor coefficients, both for the pressure and the baryon susceptibility. The Taylor expansion of the pressure Eq. (52) at fixed temperature and about µ B = 0, where only even orders contribute, is given by where the cumulants κ n are defined in terms of the moments of the winding number ω via a cumulantgenerating function K(t): We can measure all Taylor coefficients from the baryon density fluctuations, as a 3 n B = ω according to Eq. (47). We also obtain immediately from the Taylor coefficients of the pressure c 2n those of the baryon susceptibility: A comparison of discrete and continuous time evaluations of the first cumulants as shown in Fig. 23 demonstrates the cumulants are less noisy in the CT-limit. But it further requires the polymer resummations and histogram method to determine the higher order cumulants up to κ 12 , shown in Fig. 24. From a thermodynamic extrapolation of the inflection points, we obtain an estimate for T c consistent with its determination in Sec. V A.
A comment on the definition of the pressure used in this section is in order: we have previously discussed that Eq. (52) is only valid in homogeneous systems, as is expected for the continuum limit of lattice QCD. In the strong coupling limit this is not the case. We can however only measure the pressure defined by a volume derivative according to Eq. (49) in terms of dual variables, and it is of course possible to Taylor expand the spatial dimer density n Ds as well. But this definition is proportional to the energy density and shows a gap along the first order transition. In contrast, 97 is well behaved as it proportional to the thermodynamic potential F = −T log Z, which is continuous along any transition.

B. Estimates for the Radius of Convergence
We are now in a position to estimate the radius of convergence [62] from these Taylor coefficients: The corresponding results for the various n are given in Figs. 25, 26, where the radii for are plotted within the phase diagram. Above aT c , the radius becomes imaginary (indicated in gray colors). Note that we are still in the chiral limit where the whole phase boundary is either second or first order. Hence we expect that the radius of convergence drops to zero at aT c . Below aT c , the first singularity is given by the phase boundary, and we find indeed that the higher orders converge to the phase boundary. This is in particular observed for r χ B n , where the first order line is well approximated for n = 10.

A. Staggered Euclidean Time Correlators
We have explained in Sec. III that the monomer 2point correlation function is sampled during worm evolution. We are mainly interested in temporal correlation functions, from which we can extract the ground state energy corresponding to the meson pole mass. In this section we will explain how to extract them and discuss their dependence on temperature and baryon chemical potential.
The basic definition of the temporal correlators at zero momentum p = 0 for staggered fermionsχ, χ, based on the local single-time-slice operators [63] is: where the spin S of the meson is given by the kernel operators Γ S in terms of phase factors g S x,τ ∈ {±1}. We will only consider operators that are diagonal in spin-taste space: Γ S ⊗ Γ T with Γ T = Γ S * . We will not consider any flavor structure as N f = 1 (but see App. IX D for N f = 2). In every mesonic correlator specified by Γ S , there is a non-oscillating part and oscillating part with additional phase factor (−1) τ , which is due to the evenodd decomposition for staggered fermions. This parity partner has opposite spin, parity and taste content. Thus the non-oscillating and oscillating part correspond to different physical states, see Tab. IV. Of particular interest is the pion π P S which is the Goldstone boson for the residual chiral symmetry, Eq. (8). Throughout the worm evolution, monomer two-point correlation functions are accumulated whenever head and tail are at opposite parities: (102) with Z the number of worm updates. Such worm esti- mators are incremented as with f (γ) given in Eq. (41) and f (T ) given in Eq. (42). Summing over the correlators yields immediately the corresponding discrete/continuous time susceptibilities: The non-oscillating and oscillating parts of the correlators for discrete time are shown in Fig. 27. It is advantageous to consider the linear combinations and fit the even/odd correlators instead: (1) the fit is more stable (2) it generalizes to the continuous time limit, where we can distinguish even and odd τ via emission and absorption events, see Sec. II C. We can reconstruct the physical states by inverting Eq. (107). The discrete time correlators for the pion are shown in Fig. 28. We observe that the correlators for increasing N τ become more continuous and their range extends to N τ /2. In Fig. 29, the continuous time correlators for the pion π PS is reconstructed from with τ ∈ [0, 1/2] and spatial kernel g π x = (−1) x+y+z , and likewise for other mesons. This requires book-keeping on which events contribute to C Odd or C Even depends on whether the worm head is located at an absorption event x H ∈ A or an emission event x H ∈ E. Even in the CT-limit, it is necessary to discretize the temporal correlators, due to memory limitations and finite statistics. The histograms will depend on the bin size with N the number of bins. The finer ∆τ , the less events are placed in each bin, which makes the determination of the correlator more difficult. On the other hand, the coarser ∆τ , the less data are available to reconstruct the correlator. In principle one could measure the continuous time correlators without introducing a binning [64], but in practice this seems not necessary as our measurements for N = 100, 200, 400 lead to almost identical results.

B. Temperature and Density Dependence of Meson Pole Masses
Since temporal correlators are measured at zero spatial momentum, the extracted meson masses are pole masses: E 0 ( p = 0) = M . We extract the ground state mass M as dimensionless quantity M/T by multi-state fits (including excited states) and by varying the fit range  [τ min /aT , 1/(2aT )]. To obtain good balance between the required number of states and the error estimation, we apply the Aikaike Information Criterion [65]. We adjust τ min to be most sensitive to mass plateau. To compare discrete time (where we extract a τ M to continuous time, we convert via as shown in Fig. 30. Making use of the same fitting scheme, the error bars for the extracted pole masses from CT-correlators are much smaller than the corresponding DT-correlators. Moreover, the uncertainties when extrapolating DT-correlators of about 3% are circumvented. We have measured the temperature-dependence of the pole masses and find that in particular the pion becomes heavy at the chiral transition, see Fig. ( 32). For N f = 1 we find a mass degeneracy for the pairs of states: which corresponds to a multiplication by the parity (x), compare Tab. IV. This is due to the strong coupling and the chiral limit (i.e. we are in the -regime): e.g. the pion π P S is mass degenerated with the sigma meson σ S . This degeneracy is lifted as soon as am q > 0, see IX E. The pion becomes indeed massless below T c in the thermodynamic limit, as seen in Fig. 31. But the pion and all other mesons do not acquire a thermal mass, as shown in Fig. (32). Rather, they all tend to the same high temperature value aM = 0.411 (1). We suspect that this is an artifact of the strong coupling limit: even at high temperatures, in the chirally restored phase, the quarks are still confined into mesons. Hence, they do not experience the anti-periodic boundary conditions [66] and will not receive contributions from the lowest Matsubara frequencies πT above T c .
The extension to finite chemical potential aµ B is straight forward, the results on the temperature dependence of the pole masses for various chemical potentials below aµ B TCP is shown in Fig. 33. The pole masses change most at the transition temperature for the respective chemical potential. Their hign-temperature limits become independent of the chemical potential.

VIII. CONCLUSION
We have demonstrated the power of continuous time simulations of lattice QCD in the strong coupling limit, which make extrapolations for N τ → ∞ obsolete. All ambiguities arising from such an extrapolation are removed. The Hamiltonian formulation gives further insight into the world-line formulation of strong coupling lattice QCD. We discussed in detail the continuous time worm algorithm in terms of a Poisson process, the dual observables, and resummation and histogram techniques to determine the phase diagram both in the µ B −T plane and n B − T plane via the Wang-Landau method. The phase boundary can be compared with estimates from the radius of convergence from Taylor coefficients which we can determine via baryon fluctuations at zero density up to c 12 . We have also investigated temporal correlation functions, which we can measure with high resolution and higher statistics compared to discrete time, and from which we could determine the temperature dependence of the meson pole masses, both at zero and nonzero density. Whether the continuous time correlation functions can also be extended on the Schwinger-Keldysh contour to extract transport coefficients is under investigation. Real time simulations in the dual formulation of SC-LQCD are not completely sign-problem free, but much less severe compared to the standard formulation based on the fermion determinant. Some first steps to extend our Hamiltonian formulation to more flavors and finite quark mass are presented in the appendix. We plan to include the gauge corrections from the Wilson gauge action in continuous time in a similar way as we have already successfully implemented in discrete time [15,17,67]. As the continuous time limit is well defined also at finite lattice gauge coupling β, we may improve on the phase diagram by reducing the spatial lattice spacing directly in the continuous time limit via quantum Monte Carlo simulations.

ACKNOWLEDGMENTS
W. Unger is grateful to Philippe de Forcrand for providing the initial idea to consider the continuous time limit, and is thankful for the many discussions on the continuous time formulation. We would like to thank our colleagues Olaf Kaczmarek and Sören Schlichting for discussions of some aspects related to Euclidean correlators, Christian Schmidt for discussions on Taylor expansion, Owe Philipsen for discussions on the canonical phase diagram and Jangho Kim for his contributions to our code for discrete time. For the extraction of pole masses from temporal correlators, we are thankful to Hauke Sandmeyer for providing numerical tools. We acknowledge contributions of the students Aaron von Kamen (on the Wang Landau method) and Ferdinand Jünnemann (on percolation) to this project.
Numerical simulations were performed on the OCu-LUS cluster at PC2 (Universität Paderborn In this section we want to explain how to derive Eq. (29) from Eq. (27) We start from the discrete partition function for gauge group U(N c ), neglecting the baryonic part for a moment. We have to investigate what sequence of vertices is admissible on each site and at the same time conserves the pion current.
We will use the vertices in the meson occupation num-bers and introduce the shorthand notation v(k|l|m) ≡v(k|l)v(l|m).
We classify admissible sequences via the length of the interval: whether it is even or odd. This is determined by the sequence of emission sites E or absorption sites A. The discussion applies to N c = 3 but generalizes straightforwardly to odd N c . For even N c , the meson state m = N c /2 needs a special treatment which will not be address here. We distinguish via even-odd parity: 1. Odd intervals are those where an A-site is followed by an E-site, or an E-site is followed by an A-site: which is exactly the case when we have two subsequent -vertices or two subsequent -vertices.
2. Even intervals are those where an A-site is followed by an A-site site, or an E-site is followed by an E- wich is exactly the case when a -vertex is followed by a -vertex or vice versa.
Since N τ is even, the number of odd intervals must be an even number. Also, on each site, the number of E-sites equals the number of A-sites. Any CT-configuration G with N c odd is completely determined by specifying the location and kind of the vertices and whether an interval is even or odd: an interval between two vertices of the same type is always of odd length, between two different vertices it is of even length.
If we now consider a spatial bond given by the nearest neighbor pair b = x, y such that there is at least one spatial dimer on b, then the sequence of E-sites and Asites is exactly opposite if we ignore vertices which do not belong to dimers on b (see Fig. 2). This implies that it is completely determined by the type of vertex whether we have an even or odd interval. The first dimer on b can be put on any of the N τ temporal locations, but the second dimer can only be put on N τ /2 locations, and all subsequent spatial dimers (N τ − k)/2 temporal locations. Given that the maximal number of spatial dimers is given by the order in O(γ −n b ), in the limit N τ → ∞ the probability of two spatial dimers on b to be at the same location (effectively forming a double-dimer) is zero, and we can disregard the finite N τ corrections N τ − k. In this limit, we have N τ (N τ /2) k b possible temporal locations. We have however not yet considered symmetry factors as in the above argument, the spatial dimers added to the bond are not time ordered. Time ordering is however a global aspect that cannot be considered in isolation of a single bond. If we force the whole set of spatial dimers with n = b n b to be time ordered, we have to divide by n! as only one of the permutations are time-ordered sequence. Another way to see how the symmetry factor arises from time ordering for N τ → ∞ with t → τ /N τ , τ ∈ [0, 1[ is to replace the sums by integrals: where t i is the temporal location of the i-th spatial bond. This holds because the weight w(t 1 , t 2 , . . . t k ) does not depend on the locations but just on the number of vertices that appear. We conclude that the total weight of a U(3) configuration is: where Γ n is the set of topologically inequivalent configurations (which differ in the distribution of -andvertices over sites. In summary, a CT-configuration G is completely determined by specifying whether the intervals between vertices are even or odd, up to translation by a τ which corresponds to time reversal, with equal weight: w(G) = w(G T ). With N τ /γ 2 = 1/aT we arrive at the mesonic part of the continuous time partition function Eq. (29). It remains to discuss the baryonic part of the partition function. Since baryons form self-avoiding loops, it suffices to note that spatial baryon hoppings are suppressed by γ −Nc . Hence, for N c ≥ 3 spatial hoppings are essentially absent as γ → ∞ and baryons become static. This does not happen for N c = 1 (electrons) and N c = 2 (di-quarks). In physical terms, only for N c ≥ 3 the baryon is heavy and non-relativistic. Hence a baryon-antibaryon pair cannot be created from the vacuum.

B. Analytic Result for U(1)
We have derived an analytic expression for strong coupling U(1) on 2 × N τ lattices for arbitrary values of N τ and γ 2 , enabling us to obtain the continuous time result. A generalization to SU(N c ) is not straight forward. The continuous time assumption that spatial dimers with spatial multiplicity k i > 1 are suppressed is (trivially) exact in U(1).
The partition function in the chiral limit for U(1) can be derived from considering all spatial dimers, making use of the fact that the interval length in units of a τ between subsequent spatial dimers must be odd: Note that n always has to be even as spatial hoppings have to come in pairs to be consistent with the boundary conditions in time. The factor 2 n is due to the fact that each spatial dimers can hop either in forward or backward direction due to the periodic boundary conditions in space. Note that we have not approximated α n by 2 n! ( Nτ 2 ) n/2 , as we did in the steps leading to Eq. (116). We also want to consider the contribution to the partition sum with a total number of 2 monomers. For lattices with spatial extend N σ = 2, the situation is considerably simple because it is not possible to separate the monomers by spatial dimers. If we decompose configura-tions into a piece with the monomers located but no spatial hoppings, and a piece with no monomers, where the first piece has length D and the second piece has length N τ − D, we can factorize the possible configurations by considering those on the 2×D lattice and the 2×(N τ −D) lattice where it is required not to make use of periodic boundary conditions. This restricts the possible configurations further (no temporal dimers connecting the first and the last site allowed). D may be odd or even, depending on whether the two monomers are on the same spatial site (D even) or on different spatial sites (D odd). The corresponding result in the 2-monomer sector is: Both Z 0 and Z 2 are divergent series in γ, but their ratio is not, and gives the chiral susceptibility in the chiral limit: where we have used the definition of aT , Eq. (17). This result is explicitly temperature dependent, with N τ quantifying the cut-off dependence. In the limit N τ → ∞, arcsch(x) 1/x, and the chiral susceptibility has a welldefined continuous time limit: .
In Fig. 34 the agreement of Monte Carlo data with the exact result is shown, both at finite N τ and continuous time.

C. Mean field and Percolation Analysis
The mean-field analysis for SC-LQCD based on a 1/dexpansion has been studied for many decades [20][21][22][23][24][25][26]57]. Also the continuous time partition function derived here can be used as a starting point for a mean-field analysis. Our mean-field analysis assumes that a single site only couples to a mean-field bath of spatial dimers, where the location of E-sites and A-sites on its 2d nearest neighbors does not matter. This is well justified at high temperatures, where bonds with spatial dimers are isolated, but may also hold approximately at lower temperatures. Our partition sum has only one dynamical site, and all other sites have a fixed number of vertices determined by a self-consistency relation. Neglecting the baryon sector, the resulting partition function in d spatial dimension is: This implies for the energy density: resulting in a 4 M F = 3 2 + √ 3 which should be compared to the discrete time value a 4 = 3 4 at γ = 1.
A qualitative understanding of the phase diagram can also be obtained via a percolation analysis on the spatial volume. We consider mixed percolation, both on bonds and sites [68]: (1) In the chirally broken phase, the pion correlation length diverges, thus the phase is characterized by bond percolation, where a bond is activated whenever there is at least one spatial dimer at that bond (for some time location). It can be related to the average bond occupation probability θ p bond with (2) In the nuclear phase, where every site is activated if it is occupied by a baryon or anti-baryon (P-Polymer), it can be related to the average site occupation probability, n P p site . Our criterion for percolation is that in the statistical average, the probability that a cluster spans around the periodic lattice in at least one spatial direction is close to 1. The percolation threshold is characterized by a step function in the thermodynamic limit. Clearly, at low temperatures, the vacuum phase is characterized by bond percolation and the nuclear phase is characterized by site percolation. At higher temperatures, we may have a phase where bond and site percolation coexist. It turns out that this mixed phase exists along the first order transition. The percolation threshold for each bond and site percolation is reached close to the tricritical point, see Fig. 35. It is not surprising that the identification of the bond occupation probability with θ works extremely well, as the pions form a free relativistic gas [35]. In contrast, the identification of p site with the baryons density is not as good, as they do not form a free gas but are subject to strong nuclear interactions. The two-flavor formulation admits more than one baryon per site and the Grassmann constraint allows for pion exchange between them, modifying nuclear interactions substantially. It also compares better to the strong coupling limit with Wilson fermions in a world-line formulation, as discussed in the context of Polyakov effective theory [69] which integrates out the spatial, but not the temporal gauge links. SC-LQCD with N f = 2 in the dual formulation has been discussed on discrete lattices for U(1) gauge group in [70]. For U(3) gauge group, the link integrals have been addressed in [37]. Here, we report on first steps towards a Hamiltonian formulation. The suppression of spatial bonds γ −k , k > 2 also applies here. Let us first consider the static lines. We want to establish the basis of quantum states that generalize the N f = 1 states |m and |b . To arrive at this basis, we consider the SU(3) one-link integrals [71]: The sum over n i (i = 0, . . . 3) terminates due to the Grassmann integration. The corresponding invariants x i can be evaluated for N f = 2: x 0 = B uuu + B uud + B udd + B ddd +B uuu +B uud +B udd +B ddd x 1 = Tr[M x M y ] = k U + k D + k π + + k π − π + π − ,U D + k U D,π + π − − k U k D − k π + k π − .
Note that the baryonic fluxes are spinless and spin arises only when measuring baryonic correlators with the corresponding staggered kernels. We still have to integrate out the Grassmann variables to obtain the quantum states in the occupation number basis, and the corresponding Hamiltonian, where we consider the chiral limit only. The Grassmann constraint then dictates that all quarks u,d and anti-quarksū,d are within mesons or baryons. The Grassmann integral in the chiral limit is [du α dū α dd α dd α ](ūu) k U (dd) k D (ūd) k π − (du) k π + = (−1) k π + +k π − 2 (N c !) 2 1 (k π + + k π − )/2 mod N c = 0 1 Nc otherwise (129) which simplifies due to flux conservation: Just as for N f = 1, we can define vertices in the same

E. Finite Quark Mass
The chiral limit is the most interesting regime when studying the chiral transition, but we need to extend the derivation of the continuous time partition function to finite quark mass to address the quark mass dependence of zero and finite temperature observables. Only then it is possible to study the p-regime where the pion correlation function fits on the lattice. Whereas in the chiral limit, the chiral condensate is strictly zero (in a finite volume), already a small quark mass will result in a nonzero chiral condensate. Likewise the sigma meson becomes much heavier compared to the pion. This can be best understood in the dual representation: The number of monomers on even sites always equals the number of monomers on odd sites. In the pion correlator, the contributions from monomers at even sites have the opposite sign from those at odd sites, resulting in a light pion mass. In the sigma correlator, the contributions from monomers at even sites have the same sign as those at odd sites, resulting in a heavy sigma meson.
When attempting to derive the continuous time partition function at finite quark mass in a naive way, i.e. at fixed quark mass am q , the monomer number will diverge in the limit N τ → ∞. We have illustrated in [43] that the continuous time limit is well defined also at finite quark mass, but it turns out that the constant κ is now quark mass dependent. This function κ(m q ) has been determined non-perturbatively with a condition for keeping the quark mass constant in the limit a τ → ∞. With this knowledge, the continuous time limit can also be derived at finite quark mass, but it is not the bare quark mass am q , but rather the ratio M π /T which is the input parameter of the continuous time partition function. We are working on an extension of the CT worm algorithm such that the Poisson process incorporates a finite quark mass.