Oscillating propagators in heavy-dense QCD

Using Monte Carlo simulations and extended mean field theory calculations we show that the 3-dimensional ℤ3 spin model with complex external fields has non-monotonic spatial correlators in some regions of its parameter space. This model serves as a proxy for heavy-dense QCD in (3 + 1) dimensions. Non-monotonic spatial correlators are intrinsically related to a complex mass spectrum and a liquid-like (or crystalline) behavior. A liquid phase could have implications for heavy-ion experiments, where it could leave detectable signals in the spatial correlations of baryons.


JHEP10(2016)055
although in the presence of a first order transition, it is of course easy to identify the liquid as the denser, less compressible state.
To understand better when to expect such non-monotonic behavior, Ogilvie et al. have in a series of papers [7][8][9] studied models which break charge conjugation C, but remain invariant under the combined action of C and complex conjugation K. QCD at nonzero chemical potential µ has this property, but also simpler models like the Polyakov-Nambu-Jona Lasinio (PNJL) model with nonzero µ, SU(3) (Polyakov loop) spin models with nonzero µ, and even the three-state Potts model with nonzero µ have the same property. So before tackling full QCD one can hope to learn the implications of this symmetry pattern from simpler models, which may even in some cases be mapped to limiting cases of QCD itself. It is well known that for QCD, the expectation value of the Polyakov loop differs from the expectation value of its Hermitian conjugate, Tr F L = Tr F L † when the chemical potential µ is nonzero. However, the free energies are real because of the CK symmetry. As a further consequence of the breaking of C, the transfer matrix T is not Hermitian, which means that the eigenvalues are not all necessarily real. Because of the invariance under CK, however, if λ is an eigenvalue of T , then so is λ * , i.e. the eigenvalues are either real or occur in complex conjugate pairs. This is interesting because it implies, in the case where complex eigenvalues occur, that the Polyakov loop correlator is non-monotonic.
The dependence of a correlation function on the eigenvalues of the transfer matrix is most easily demonstrated with a one-dimensional lattice model, but generalizes to any correlator with fixed momenta in the orthogonal directions. Let φ i , i = 1, . . . , N, be a field living on a circle with N sites. If T is the transfer matrix connecting neighboring sites, then the correlation function of φ is given by φ(x)φ(0) = Tr T N −x φT x φ Tr (T N ) = Tr Λ N −xφ P −1 φP Λ xφ P −1 φP where P −1 T P ij = Λ ij = λ i δ ij is the diagonalized transfer matrix with the eigenvalues λ of T sorted in magnitude such that Reλ 0 ≥ Reλ 1 ≥ · · · andφ ij is φ in the eigenbasis of T . For simplicity we assume a discrete spectrum in the description below. We also consider the N → ∞ limit. In general three scenarios are possible. 1 Firstly, all eigenvalues can be real and the correlator is a conventional, exponentially decaying function, Here we parametrize λ n /λ 0 = e −mn , with m n ≥ 0. Secondly, if the largest eigenvalue is real and the next two are a complex-conjugate pair, 2 then the correlator also decays 1 The eigenvalues of the transfer matrix are either real or come in complex conjugate pairs, since the model is invariant under the simultaneous action of charge and complex conjugation. 2 Strictly speaking, there can be more real eigenvalues above the complex-conjugate pair, with a consequently weaker oscillation in the correlator. We do not treat that case separately here.

JHEP10(2016)055
exponentially but is modulated by a cosine, and the system behaves as a liquid, where λ 1 /λ 0 = e −m R +im I . Finally, if the eigenvalue with the largest real part is part of a conjugate pair λ 1 /λ 0 = e im I , then the correlator is a pure trigonometric function and a crystalline behavior is observed, of the form |φ 00 | 2 + A cos(m I x + θ), which reveals long-range order for arbitrarily large system size N . If m I is small in units of the lattice spacing then one oscillation spans several lattice spacings and has nothing to do with the underlying structure of the lattice. For a continuous spectrum, the above categorization is still valid. In case one, the eigenvalues are all distributed on the real line whereas in case two, the eigenvalues branch off into the complex plane somewhere below the largest eigenvalue. The transition between cases one and two occurs at a so called disorder line. Case three is obtained when the branching point reaches the largest real eigenvalue. All these three cases have been found by Ogilvie et al. [7] in 1-dimensional models, where the complete phase diagram can be obtained using transfer-matrix methods. Such a 1-dimensional model can for example serve as a dimensionally reduced effective models of 1 + 1-dimensional QCD at finite temperature. Recently [9] it has been proposed that also higher dimensional models show these characteristics, based on the fact that the 1dimensional solution can be seen as the first order in a character expansion. It has, however, to our knowledge, not been demonstrated with first-principles lattice simulations that this is actually the case.
As mentioned above, it has been suggested [3] that the conditions in the fireball after a heavy-ion collision might be such that the baryon-number correlations have an oscillatory character. This conjecture is based on an effective flux-tube model introduced in [1,2] which can be mapped into an XY -model with external magnetic fields which break charge symmetry, such that it falls in the same category of models discussed above. Another flux-tube model, which can be mapped into a three-state Potts model, is treated in [10]. In general, the Hamiltonian and partition function for such a flux-tube model are given by where l x,ν denote flux tubes with string tension σ living on the links, q x denote quarks with mass m and chemical potential µ living on the sites and j x denote junctions with vertex JHEP10(2016)055 energy v living on the sites. All occupation numbers are integer valued and, depending on their allowed range and on whether v is zero or nonzero, the model can be mapped to either an XY model (for v = 0) or a Z N spin model (for v = 0). The junctions j call for further explanation. In SU(N ), they are related to the invariant -tensor, i.e. N flux lines emanating from N (anti-)quarks join at a junction and form an SU(N ) singlet, and thus the (anti-)quarks together with the flux lines are identified with a (anti-)baryon.
In this report we study the Z 3 spin model with nonzero chemical potential µ in 1 and 3 dimensions and show that a complex mass spectrum can occur in both cases. The rest of the report is organized as follows. In section 2 we define the model and review the solution in 1 dimension using the transfer matrix as well as the approximate solution in any dimension using Extended Mean Field Theory (EMFT) [11]. In section 3 we present our Monte Carlo results and then we draw our conclusions in section 4.

Model
The model we will be studying is the three-states Potts model with nonzero chemical potential or, more accurately, the Z 3 spin model with complex external fields. 3 In d dimensions, it can be seen as the crudest approximation of (d + 1)-dimensional QCD in the static-dense limit. The action is given by and the Z 3 spins P ∈ {1, e i 2π 3 , e −i 2π 3 } at each site represent the center of the Polyakov loops Tr F L. The usual interpretation of the external fields is h R = e −M/T cosh(µ/T ), h I = e −M/T sinh(µ/T ), where M and µ are the mass and chemical potential of the quarks respectively [12], but as mentioned in the introduction, it is also possible to map (β, h R , h I ) of eq. (2.1) into (σ, m, µ) of eq. (1.5), as described in [10] (see eqs. (14)(15)(16)(17)(18)). We will primarily use the first mapping but will evaluate the results also in the light of the second one. Note, however, that the mapping between eq. (1.5) and eq. (2.1) is not possible for all parameter values (see figure 2).
In the formulation (2.1) the action is complex, and the model clearly suffers from a sign problem, but as long as h R , h I ∈ R and h R > |h I |, which corresponds to the physical case of M, µ ∈ R, there exists a sign-problem-free representation 4 that can be sampled by a worm algorithm. The model can however be interesting in its own right also in the unphysical region h I > h R , but it is a shortcoming that it does not have a continuum limit in three dimensions, which could make it harder to clearly separate the lattice spacing a from the correlation length ξ and in extension, the wavelength λ of the oscillations we are looking for. In one dimension the model can be solved for general external fields using a transfer-matrix method and we can use EMFT to obtain an approximate solution in any number of dimensions.

Transfer matrix
In 1d the partition function of a chain of N Z 3 spins with periodic boundary conditions is given by where T is the transfer matrix. It is easy to verify that the characteristic polynomial of T is a cubic polynomial with real coefficients so there are either three real roots or one real root and a pair of complex conjugate roots, as claimed above. For a given β, it is now straightforward to determine the phase diagram which contains the three phases described in the introduction. The phase diagram at fixed β = 0.08 can be seen in figure 1. The color coding and labels are as follows: I (blue) marks the region where all eigenvalues of T are real. In Ia they are all positive and the connected correlator is a pure sum of exponentials.
In Ib two eigenvalues are negative (the product of all three, i.e. the determinant of T , is always positive) and the connected correlator is in general a sum of two oscillating functions with wavelength 2, due to factors (−1) x . Depending on how the signs and magnitudes of the eigenvalues are distributed, this may or may not be detectable on a discrete lattice. II (green) denotes the region where the largest eigenvalue is real and the other two are a complex conjugate pair. The connected correlator is a cosine-modulated exponential, this is characteristic of a liquid. III (red) marks the region where the complex conjugate pair is larger in magnitude than the real eigenvalue and the connected correlator at long distance is a pure trigonometric function, this is the long-range order characteristic of a crystal. The two black lines bound the wedge where h R > |h I | and mark the region where the flux-variables representation is sign-problem free and the worm algorithm can be used. It is evident that the crystalline phase is out of reach of the worm algorithm but some parts of the liquid phase lie within the physical region h R > |h I |, so that the non-monotonic behavior of the connected correlator there can be reproduced by lattice simulations. Initially the transfer-matrix method is only defined for integer separations but it is straight forward to extend it to any real separation via the matrix power-function.
In the liquid phase, the connected correlator is given by eq. (1.3), which is made periodic (exp → cosh) at finite N to obtain where f (P ) is either P, ReP or ImP . The parameters a f and φ f can be calculated from the eigenvectors of T . These functions can be directly compared to the correlators obtained by the worm algorithm and will serve as a consistency check for the algorithm before going on to three dimensions where no exact results are available. A comment about this "phase diagram" is in order. Actually, the different phases are not separated by phase transitions in the strict sense; there is no singularity in the free energy anywhere in the (h R , h I )-plane, since the zeros of the characteristic polynomial of T are smooth functions over the whole plane. Instead, the boundary of the different phases are disorder lines, which mark a smooth change in the characteristic of the correlator, for example from a non-oscillatory exponential decay to an oscillatory exponential decay. In general, however, it is not necessarily so that the change from non-oscillatory to oscillatory behavior take place at a disorder line, it can also occur at a first order transition, as is evident from for example the water-vapor transition.

EMFT
In more than one dimension, and especially in the physically interesting case of three dimensions, the transfer-matrix method is not practical anymore. It is reasonable to assume that the structure of the phase diagram will remain [9] but one is totally at a loss when it comes to the exact location of the disorder lines. In the light of the one-dimensional results, it is unlikely that the crystalline phase can be probed by lattice simulations, but one may hope to find evidence of a liquid phase. In this case the three largest eigenvalues of the transfer matrix will be given by (up to a trivial overall multiplicative, real constant) where m R , m I > 0 are real numbers chosen to paramterize the eigenvalues. The decay of the spin-spin correlator will thus be governed by P (0)P † (r) ∼ e −m R r cos(m I r). It becomes clear that our prospects for detecting this characteristic behavior of the correlator depend rather sensitively on m R and m I ; we require a point in phase space where m R is not too large at the same time as m I is not too small, so that the first maximum in the correlator occurs before the signal is too damped.

JHEP10(2016)055
Much time and effort can be saved by quickly, albeit approximately, solving the model for extended regions of parameter space. Mean field theory is one candidate which falls short since it does not give access to the mass spectrum. EMFT [13] on the other hand does exactly that and is thus an apt choice.
It will be useful to consider the real part ReP and the imaginary part ImP of the Potts spin P as independent variables here. Since the imaginary part of the action (2.1) is odd in ImP , the expectation value of iImP will be real and we have P = ReP + iImP = ReP − iImP = P † . The Z 3 spin P is then decomposed into its mean value and fluctuations around the mean, We now formally integrate out all fields except the one at the origin and assume that this amounts to the introduction of effective couplings for the bilinears δReP δReP, δImP δImP and δReP δImP [11]. The effective EMFT action can then be written So far we have not assumed anything about the variables P , so the effective action above is generally valid for any action of the form (2.1). For P ∈ Z 3 the action can be simplified slightly by using (ImP ) 2 = 1 − (ReP ) 2 and ReP = − 1 2 whenever ImP = 0. We then obtain It is obvious how to self-consistently determine the linear expectation values, whereas the bilinears may need some more explanation. The details of their determination will reveal how a complex spectrum can arise. As usual in EMFT [11], we fix the effective quadratic couplings ∆ i by matching the bilinear expectation values to an approximation to the pointto-point correlator of the full model, (2.10)

JHEP10(2016)055
This is a matrix equation where G 0,c (k) is the connected Green's function of the free theory. It is not immediately clear what the free theory of a spin model is, but the approximation above is in fact valid for any choice. A good choice will be close to the model we want to study and at the same time allow for an efficient numerical treatment.
We have chosen the free model to have the same action as the original Z 3 model, eq (2.1), but with the variables P ranging freely over the complex plane. With this choice the free connected Green's function is given by G −1 0,c = −2βId 2 ν cos k ν . The self-energy Σ(k) in eq. (2.10) then arises due to the restriction of the field to take values in Z 3 . The EMFT self-energy Σ EMFT is likewise identified as the difference between the variance of eq. (2.5) with (ReP, ImP ) ∈ R 2 and the variance when P ∈ Z 3 and is given by Hence, the final self-consistency equation becomes It is clear that G −1 EMFT,c + ∆ − 2βId 2 ≡ βM plays the role of a mass matrix and we should diagonalize it to obtain the mass spectrum. It will also be vastly more efficient to integrate over the momenta when M is diagonal. The mass matrix can be parametrized as where a, b, c ∈ R. The eigenvalues are then given by m ± = a ± √ b 2 − c 2 , such that if |c| > |b| the spectrum will consist of a pair of complex conjugated masses m R ± im I with m R = a and m I = √ c 2 − b 2 . This implies cosine-modulated exponential fall-off in the correlators in the (ReP, ImP ) basis, as expected. By solving the model in the (h R , h I )plane, a phase diagram analogous to what was obtained in one dimension with the transfer matrix, figure 1, can be constructed by studying the behavior of the masses. In figure 2 we show the results for fixed β = 0.08, with the most interesting features being the disorder lines (in red), where the masses are degenerate, and the blue dashed lines where the real part of the complex masses vanishes. Beyond these lines the momentum integral in the self-consistency equation no longer converges, since the integrand is no longer decaying at large distances. One may guess that with purely imaginary masses, the system would enter a crystalline phase with a purely trigonometric correlator but there is no way to verify that using EMFT. This phase diagram can then be compared both to the mapping (h R , h I ) → (e −M/T cosh(µ/T ), e −M/T sinh(µ/T )) and to the alternative mapping in [10]. It is found that the second case covers a subspace of M, µ ∈ R and there are indeed regions in parameter space where the mapping is valid and where one expects a complex spectrum. However, in that region h R is substantially larger than h I which means in the M, µ variables that both M and µ are rather small, which is presumably far away from the region where the model is expected to be a valid approximation of QCD. Now that an approximate phase diagram has been obtained, we can select points in the liquid phase which are favorable in terms of m R and m I where full Monte Carlo simulations using the worm algorithm will be performed.

Results
For our lattice simulations we used the flux-variables representation described in [14], with quark occupation number n x ∈ {−1, 0, 1} on each site and flux occupation number l x,ν ∈ {−1, 0, 1} on each link. Gauss' law requires that the flux at each site is a multiple of three. Allowed configurations consist of flux-tube networks with or without attached quarks. If there are no quarks attached the flux-network can be thought of as a glueball. There are also neutral networks with any number of quarks and an equal number of antiquarks attached, for example networks connecting one quark with an anti-quark can be thought of as mesons. The third possibility is to have a surplus of 3n (anti-)quarks. This is equivalent to having the junctions of the network to sum up to n, we say that the network has junction charge n. These charged networks are associated with baryons.
The worm algorithm generates a Markov chain of allowed configurations by temporarily violating the constraint, something which can be exploited to obtain improved estimators for spin-spin correlation functions. In addition to the usual P (0)P † (x) we use a  modification introduced in [15,16], which allows us access to improved estimators of also ReP (0)ReP (x) and ImP (0)ImP (x) . This is crucial because the best signal-to-noise ratio will be found in the correlator of the imaginary parts of the spins, since it has the smallest constant background. We first reproduced the results obtained by the transfer matrix method in one dimension in order to verify that the algorithm was properly implemented. A typical correlator in the liquid phase is shown in figure 3 and there is perfect agreement with the analytic result for all three propagators. It should be noted that the real part of the mass m R is in general always large when the imaginary part m I is of order one or larger, this makes it very difficult to resolve the first local maximum of the correlator.
We also measured the junction-junction correlator on the configurations generated by the worm algorithm. The junction j x takes the value n if 3n, n ∈ Z units of flux flow into the site x. With the flux variables described above there are in general four types of junctions, depicted in figure 4 but in one dimension only junction A with one quark and two in-going fluxes (or its reverse) attached to the site is possible. In figure 5 we show the correlation between positive j + and negative j − junctions for two different parameter values. Here the oscillation is even clearer due to a less noisy observable, although we do not have an improved estimator for this correlator. The dashed line is obtained by fitting the amplitude and phase in eq. (2.3) while keeping the masses fixed at the exact values obtained by the transfer matrix. The mass is the same as for the spin-spin correlator since the junction is a local object and there is only one (complex) mass in the one-dimensional case. It should be noted that these parameter values have been selected to give a maximally clear first maximum in the oscillation. For general parameter values in the liquid phase it is only possible to see the first minimum, while the first maximum is drowned in noise. This will be especially true in three dimensions where the real part of the mass is larger than in the one-dimensional case. Figure 4. The different junctions allowed in the flux-variable representation of the Z 3 model described in the text. The red crosses represent quarks and the lines represent the directed fluxtubes. The junction is located in the center of each network where the flux sums up to three. Note that the quarks bounding the network may also be replaced by arbitrary larger networks of charge one. In one dimension only junction A is possible. The three-dimensional junction D is only present in dimension three or higher.  We then move to the physically interesting case of three dimensions, and guided by the phase diagram calculated by EMFT we select a few points assumed to be in the liquid phase and look for the corresponding signals in the correlators. Also here, however, the damping of the correlator is always strong, as is illustrated in figure 6. In the left panel m R and m I , obtained by EMFT, are plotted as a function of tanh µ/T for fixed M/T and β and in the right panel m I is plotted as a function of m R to emphasize the approximately linear growth relation. Both masses increase with µ/T . As a consequence, it is typically only possible to resolve the first minimum of the oscillating correlator. In figure 7 we show correlators of ImP , as a function of the Euclidean distance r = x 2 + y 2 + z 2 , obtained by our worm simulations for three different values of the chemical potential µ/T ∈ {2.0, 2.5, 3.3} at fixed β = 0.08 and e −M/T = 0.05. There is a clear staggered component in the correlators, which makes it very hard to fit the data to a simple ansatz. This short-distance effect, whose sign is µ-dependent (figure 7, blue vs red), stems from the size and shape of the junctions, shown in figure 4: in d = 1, only A is possible. For all µ there is a clear minimum whose position moves toward zero and whose width decreases as the chemical potential increases. This suggests that the imaginary part of the mass increases with µ, as it does in one dimension and as EMFT predicts. Also, in neither of these correlators is it possible to see a maximum. This is not very surprising, but a discernible maximum would be indisputable evidence of a complex spectrum. The depletion in the left correlator and the enhancement in the right correlator around distance 2 support the proposition that the system behaves as a liquid. However, the strong staggered character still leaves some doubt. In the right panel, the data point at √ 2 and 2 are far below 0.99 and are omitted such that it can be clearly seen that the points at √ 3 and √ 5 are above 1.

JHEP10(2016)055
If the mass spectrum is complex and the system behaves like a liquid then the junctionjunction correlator should also show the characteristic, oscillatory behavior seen in the 1d model. Since the junction-junction correlator is less noisy than the spin-spin one, one may even hope that a maximum of the oscillating correlator can be resolved, thus establishing the complex spectrum without doubt. In figure 8 we show two junction-junction correlators for β = 0.08, e −M/T and µ/T = 3.3. In the left panel we show the correlator of the absolute values of the junctions whereas in the right panel the sign of the junctions is also taken into account. The difference of scales of the two correlators comes from the fact that they are both normalized to one at large distances and that |j| ∼ 3 j . To emphasize the staggered component of the correlators we plot the correlators on the two different sublattices with different colors. Inspecting first the charge-insensitive correlator (left panel ) we see that there is indeed a depletion in the density of junctions of any type within distance [1, 2.5] of a junction but it is not possible to tell if this minimum is followed by a maximum. In the charge-sensitive correlator (right panel ) there is a clear maximum in the correlator in roughly the same interval, but only in one of the sub-lattices. This strong staggered dependence is of course a lattice artifact. It should however be noted that for smaller values of µ/T (and thus longer wavelength oscillations, cf. figure 7) there is a clear, broader, minimum in both correlators and on both sub-lattices, which indicates that the effect is not merely a staggered effect, although the maximum which is predicted to follow is completely damped away. All in all, the behavior of the different correlators strongly suggests that there is a complex mass spectrum at the investigated parameter values, and that the prediction in [9] that the phase structure observed in one-dimension has an analogue also in three dimensions is correct.
Finally, we measured some statistics of the flux-tube networks and the junctions. Using the labeling of figure 4 we find that the ratio of C to D junctions is very close to 3/2 and

Conclusions
Using unbiased Monte Carlo simulations, we have shown that the Z 3 spin model in three dimensions with charge-symmetry breaking external fields has non-monotonic correlators, both in the original spin variables and in the flux variables, for some regions in parameter space. This strongly suggests that the spectrum in these regions is complex and this claim is also backed up by EMFT calculations of the three-dimensional model. Of special interest is the oscillatory nature of the junction-junction correlator, which is the analogue of the baryon-baryon correlator in heavy-dense QCD, for some regions of parameter space. The possibility that models with complex saddle points may have a complex mass spectrum and thus non-monotonic correlators has previously been established analytically in the one-dimensional case and has been argued to also hold true in higher dimensions. We have shown that the worm algorithm is capable of reproducing these results, even though the original spin model suffers from a strong sign problem. The phase diagram in one dimension contains regions where the system behaves like a liquid, with exponentially damped oscillations, and like a crystal with a purely oscillatory correlator. In general, it JHEP10(2016)055 is expected [9] that these features carry over also to higher dimensions. For those regions of parameter space where the model has a sign-problem free representation we have only found evidence of the liquid phase with exponentially damped but oscillating correlations between spins, as well as junctions. We have found no evidence of a crystalline phase and it is probable that it lies beyond the reach of the worm algorithm in three dimensions, as it does in one dimension. It should also be noted that even the liquid phase may lie in an unphysical region of parameter space. At least according to EMFT it lies outside of the region of parameter space which can be mapped to the more physical flux-tube model of Condella and Detar [10]. A complex mass spectrum can only be found in a parameter region where the mass M of the underlying heavy quark satisfies M T , whereas the validity of the effective description of QCD by a Potts model requires M T . This situation may change if the junctions are given a nonzero weight as in [3], but this possibility has not been investigated here. However, whatever the values of the other parameters, introducing a junction weight will further damp the signal we want to measure, making the search yet more difficult.
Our findings supports the claim that in general, it is plausible that models without charge-conjugation symmetry, but invariance under the combined action of charge conjugation and complex conjugation, will have regions with a complex mass spectrum in their phase diagram. However, more work is needed before precise statements can be made about whether or not this is a phenomenon which occurs under physical conditions. This first proof of principle should encourage the study of more realistic models, and the search for experimental signals in heavy-ion collisions as advocated in [3].
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.