Fourier coefficients of the net-baryon number density and chiral criticality

We investigate the Fourier coefficients $b_k(T)$ of the net--baryon number density in strongly interacting matter at nonzero temperature and density. The asymptotic behavior of the coefficients at large $k$ is determined by the singularities of the partition function in the complex chemical potential plane. Within a QCD-like effective chiral model, we show that the chiral and deconfinement properties at nonzero baryon chemical potential are reflected in characteristic $k$-- and $T$-- dependences of the Fourier coefficients. We also discuss the influence of the Roberge-Weiss transition on these coefficients. Our results indicate that the Fourier expansion approach can provide interesting insights into the criticality of QCD matter.


I. INTRODUCTION
Exploring the phase diagram of quantum chromodynamics (QCD) at finite temperature and chemical potentials is a challenging problem in experimental and theoretical studies. The Beam Energy Scan (BES) at the Relativistic Heavy Ion Collider (RHIC) [1] has been dedicated to the search for the conjectured QCD critical point (CP) and the onset of deconfinement through the systematic studies of various observables such as the fluctuations of conserved charges. In particular, those of the net-proton number have been shown to exhibit a non-monotonic behavior with collision energies which may potentially be attributed to critical chiral dynamics [2][3][4][5][6][7][8].
While LQCD provides first principle insights into the thermodynamics of QCD at finite temperature T and at vanishing and moderate values of the baryon-chemical potential µ B < 3T , the well-known sign problem still inhibits systematic calculations at larger baryon densities. Several strategies to overcome this problem are being pursued [17,18].
In particular, LQCD calculations at imaginary baryon chemical potential offer a possibility to circumvent the sign problem at real µ B . On the one hand, these calculations, which do not suffer from the sign problem, can, in principle, be analytically continued to the real axis [19][20][21]. On the other, a Fourier expansion in imaginary µ B can be applied to the partition function in order to study the properties, in particular the phase structure, of QCD at finite baryon density in the canonical ensemble [22][23][24][25]. Moreover, the canonical partition function can be used to obtain the probability distribution function P (N ) of the net-baryon number N , which in turn yields the cumulants of the net baryon-number fluctuations [26][27][28][29][30][31][32]. As noted in Ref. [30], knowledge of the probability for large-amplitude fluctuations, i.e., fluctuations with large N , is required for a correct identification of the critical properties associated with the chiral phase transition.
Furthermore, a novel strategy for locating the confinement-deconfinement transition by exploring the complex phase of the Polyakov loop at imaginary chemical potential has been proposed in [33,34]. Thus, by studying QCD at imaginary chemical potential, one can gain insight into the critical properties of QCD related to deconfinement, chiral symmetry restoration and the Roberge-Weiss transition.
In this paper we explore the Fourier decomposition of the (dimensionless) net baryon-number density where p = p(T,μ B ) is the pressure,μ B = µ B /T is the reduced baryon chemical potential and T the temperature. The analytic continuation of χ B 1 to imaginary baryon chemical potential is an odd, periodic function of θ B = Imμ B , and can thus be expanded in a Fourier sine series The aim of our study is to identify the influence of chiral symmetry restoration and deconfinement on the Fourier coefficients b k (T ). We examine the effects of criticality by making use of QCD-like chiral effective models. In particular, we discuss the effect of singularities in the complex chemical potential plane associated with firstand second-order phase transitions, crossover transitions and the contribution of the so called thermal singularities. Model-independent results will be reported in a separate publication [35].
To be specific, we employ the two-flavor Polyakov-Quark-Meson (PQM) model, which emulates the characteristic properties of QCD both at real [36] and imaginary [37] baryon chemical potentials. At nonzero T and θ q = θ B /N c , the resulting thermodynamic potential exhibits chiral symmetry restoration and statistical confinement, as well as the Roberge-Weiss symmetry [38], which implies a periodicity of 2π/N c in θ q .
We show that the k-and T -dependence of the Fourier coefficients b k (T ) exhibits characteristic features, reflecting the chiral and Roberge-Weiss transitions at imaginary baryon chemical potential.
The results are discussed in light of the Fourier expansion coefficients that were recently obtained in LQCD simulations at imaginary chemical potential for a wide range of temperatures around and above the chiral and deconfinement transitions [39]. Furthermore, we examine the recently proposed cluster expansion model [40] for Fourier coefficients and discuss its predictive power and analytic properties in the context of chiral criticality.
The paper is organized as follows: in the next section, we analyze phenomenological models for the Fourier coefficients and their applicability to the description of criticality and interpretation of recent LQCD results. The properties of the Fourier expansion coefficients in a QCDlike effective chiral model are examined in Sec. III. Finally, section IV is devoted to a summary and conclusions.

A. Models for Fourier coefficients
A model for the Fourier coefficients of the net baryon density, the cluster expansion model (CEM), was recently proposed in Ref. [40]. Based on the popular hadron-resonance gas equation-of-state with excludedvolume corrections, the authors suggested a simple prescription for computing the higher order Fourier coefficients of the net baryon density b k (T ) in Eq. (2) in terms of the first two, b 1 (T ) and b 2 (T ), are the Fourier coefficients of the density of a noninteracting gas of massless quarks with N c = 3 colors and N f = 3 flavors. The coefficients (3) are constructed such that they all approach the corresponding Stefan-Boltzmann (SB) values (4), when b 1 (T ) and b 2 (T ) approach b SB 1 and b SB 2 , respectively.
The first four Fourier coefficients, b 1 (T ) − b 4 (T ), have been obtained in LQCD calculations [39]. It was shown in Ref. [40] that Eq. (3), using the lattice results as input for b 1 (T ) and b 2 (T ), provides a description of b 3 (T ) and b 4 (T ) consistent with LQCD. Furthermore, the model predictions for the temperature dependence of the sixthorder baryon-number susceptibility are qualitatively consistent with LQCD findings.
Thus, the CEM describes the basic features of the available lattice results on baryon-number fluctuations. Nevertheless, one may ask whether the critical behavior of QCD associated with the restoration of chiral symmetry, which is reflected in the asymptotic behavior of the Fourier coefficients, can be captured by the model. More generally, it is of interest to assess to which extent the modeling of Fourier coefficients is unique, when only the first two coefficients are provided as input. In order to illustrate these points, we introduce an alternative model for b k (T ), which we dub the Rational Fraction Model 1 Also in this model, lattice results for b 1 (T ) and b 2 (T ) are used as input. As shown in the left panel of Fig. 1, the RFM model provides a description of b 3 (T ) and b 4 (T ), which is comparable in quality to that of the CEM.
There is, however, an essential difference between CEM in Eq. (3) and RFM in Eq. (5), illustrated in the right panel of Fig. 1. The asymptotic behavior of the Fourier coefficients b k (T ) is power law in the RFM and exponential in the CEM. Consequently, observables that depend on high-order Fourier coefficients are expected to behave very differently in these models. In particular, this is the case for higher cumulants of the net baryon-number density. The difference in asymptotic behavior of the Fourier coefficients is a consequence of the difference in analytic structure in the complex chemical potential plane.
For real values of the baryon chemical potential, the net-baryon number density (1) is given by 1 Note that the prefactor in (5) is an even function of k, and thus preserves the symmetry of the sine Fourier coefficients (4) [40] and the rational function model (RFM), defined in Eq. (5), compared to LQCD data [39]. Right: the Fourier expansion coefficients obtained in these models at T = 230 MeV.
whereμ B = µ B /T and b k (T ) are the Fourier coefficients of the density at imaginary chemical potential (2), Higher order fluctuations of the net-baryon number are obtained by taking derivatives of Eq. (8), In Fig. 2 we show Im[χ B 1 (T, i θ B )], obtained with Eq. (2) using the Fourier coefficients (3) and (5) at several temperatures. Since the first four coefficients coincide in CEM and RFM, the difference between the two models for the baryon number density at imaginary chemical potential is relatively small and discernible only at θ B ≃ π. There, by construction, both functions drop rapidly to zero and the convergence of the Fourier sum is slow.
Differences between the higher order Fourier coefficients imply very different predictions for the baryon number susceptibilities χ B n , in particular at large n. In Fig. 3 we show the fourth and tenth order cumulants normalized by χ B 2 in CEM and RFM. Each point represents the result of a model calculation, where the Fourier coefficients were obtained using (3)(4) and (5-7), respectively, while b 1 (T ) and b 2 (T ) are given by LQCD data. The lines are obtained by interpolating the LQCD values for b 1 (T ) and b 2 (T ) as functions of temperature in the range T ∈ [165 MeV, 220 MeV]. We thus obtain the following fit to the LQCD data log b 1 = −41.537 + 83.059x − 56.803x 2 + 13.073x 3 , with x = T /155 MeV. As shown in the left panel of Fig. 3, the two models yield very similar results for the ratio χ B 4 /χ B 2 . We conclude that the cumulants χ B 2 and χ B 4 are well constrained by the leading Fourier coefficients. However, higher-order cumulants are qualitatively different in the two models. This is expected due to qualitative differences in the higher-order Fourier coefficients. In the right panel of Fig. 3 we show that χ B 10 /χ B 2 in the two models is indeed very different in the crossover region. This clearly demonstrates that the characteristic features of higherorder fluctuations, which are sensitive to criticality, are not uniquely determined by requiring that the first four Fourier coefficients are reproduced. Thus, we find that models of this type are not suited for exploring chiral criticality.  [40] and the rational function model (RFM) of Eq. (5). The points are computed directly using the LQCD datapoints [39] for b1 and b2 and the two models for the higher-order Fourier coefficients. The solid curves were obtained by fitting b1 and b2 to the LQCD data, as described in the text.

B. Analytic properties and criticality
Given the Fourier coefficients of the net-baryon density, one can compute susceptibilities and consequently obtain information on the location of singularities in the complex µ B plane, e.g., by computing the radius of convergence of the Taylor expansion of the pressure. In order to assess whether the CEM exhibits any singularities that could provide insight into the critical behavior of QCD, we examine the analytical properties of the model. Of particular interest in this context are the questions concerning the location of singularities associated with the QCD chiral transition.
We first consider the Fourier expansion of the density in the CEM in the Stefan-Boltzmann (SB) limit. In this case, the baryon density is computed using the Fourier coefficients (4) in (8), where Li n (z) denotes the polylogarithm of order n defined by We note that Eq. (12) is valid also for complex values of the chemical potential. Moreover, for |Imμ B | < π, Eq. (12) reduces to while for |Imμ B | > π, the function is consistent with the periodicity in the imaginaryμ B direction, which is implemented in the model. Now, the density of an ideal gas of massless quarks and antiquarks is given by the well-known expression wherep = p/T . For |Imμ B | < 3π, Eq. (15) reduces to the polynomial form (14). Consequently, χ B,SB 1 and χ B,IG 1 coincide within a band in the complexμ B plane defined by |Imμ B | < π (as well as in bands obtained by shifting Imμ B by multiples of 6 π). Outside these bands, the two densities differ, due to the different periodicities: χ B,SB 1 is invariant under translations of the baryon chemical potential by multiples of 2 π i T , i.e.,μ B →μ B + 2 π i N , while χ B,IG 1 is invariant under shifts by multiples of 6 π i T , i.e.,μ B → µ B + 6 π i N , where N is an arbitrary integer. The periodicity of χ B,SB Hence, the singularities of χ B,SB 1 closest toμ B = 0 are located on the imaginary axis, atμ B = ±i π. On the other hand, the closest singularities of the ideal quark gas, χ B,IG 1 , are found atμ B = ±3 i π. We note that the latter is generated by the pole of the Fermi-Dirac function atp = 0. For nonzero quark mass m, these thermal branch points are shifted away from the imaginary axis toμ B = ±m ± 3 i π, wherem = m/T . Thus, the singularity structure of the net baryon density is, in the Stefan-Boltzmann limit, completely determined by the analytic properties of the polylogarithm. In the following, we show that this is the case also for the CEM.
We first note that the Fourier coefficients of the baryon density in Eq. (3) can be expressed in the following form, Also in this case, a closed-form expression for the density can be obtained by resumming the Fourier series, Now, just as in the Stefan-Boltzmann limit (Eq. (12)), the singularity structure of χ B,CEM 1 can be readily deduced using the analytic properties of the polylogarithms. The two expressions differ only in the prefactor c(T ), and the fugacity parameter λ(T ). The latter, being less than unity, shifts the location of the branch cuts away from the imaginaryμ B axis, into the complexμ B plane. Thus, the singularities nearest toμ B = 0 are located atμ B = ± log λ ± iπ. This closely resembles the effect of a nonzero quark mass on the thermal branch points.
The radius of convergence of a Taylor expansion about µ B = 0 is given by the distance to the closest singularity, In Fig. 4 we show R(T ) computed using the Fourier coefficients obtained in LQCD as input. The results obtained with Eqs. (20) and (18) for each lattice point coincide with those obtained in Ref. [40] by summing the Fourier series numerically and using the Mercer-Roberts estimator [41] of the radius of convergence. It is thus clear,  [40] and the distance of the nearest singularity of the polylogarithm to the origin (20). The solid curve was obtained using the fitted values (11) for b1 and b2, as explained in the text.
that the CEM does not contain information on singularities in the chemical potential plane with |Imμ| < π. In other words, the CEM exhibits only the singularities of Eq. (19), which are associated with the periodicity in the imaginary part of the baryon chemical potential. Consequently, the model cannot provide any insight on the location of chiral crossover transition at non-vanishing net baryon densities nor on the possible existence of a critical point in QCD matter.
As we discuss in the following section, criticality is reflected in characteristic properties of high-order Fourier coefficients, which are not reproduced by the CEM.

III. FOURIER COEFFICIENTS IN A QCD-LIKE EFFECTIVE MODEL
In the previous section, we have demonstrated that the Fourier coefficients b k cannot be fixed uniquely by knowledge of the first few. In this section we explicitly examine the effects of critical behavior on the Fourier coefficients within a QCD-like chiral effective model. For this purpose we employ the two-flavor Polyakov-loop augmented quark-meson (PQM) model [36], which captures characteristics of QCD at imaginary baryon chemical potential [37].

A. Criticality in the PQM model
In the PQM model, the thermodynamic potential at non-vanishing temperature T and imaginary quark chemical potential θ q = Im µ q /T exhibits the Roberge-Weiss periodicity [38], a residual of the Z(N c ) symmetry, which is an exact symmetry of QCD in the limit of infinitely heavy quarks. The Polyakov-Nambu-Jona-Lasinio model The temperature is normalized by the (pseudo-)critical temperature Tpc(mπ), for each value of the pion mass. In the chiral limit, Tpc(mπ = 0) ≡ Tc. [42,43] is expected to yield compatible results, since the phase structures, including the Roberge-Weiss transition, are very similar in the two models [44,45].
In the following we adopt the mean-field approximation, including renormalized quark vacuum fluctuations [46]. The thermodynamic potential per unit volume Ω in the PQM model is obtained by extremizing the functional with respect to the thermal expectation values of the Polyakov loop Φ, its conjugate Φ * , and the sigma field σ.
Here U(Φ, Φ * ) and U (σ) are the Polyakov loop potential and the purely mesonic potential for the O(4) multiplet (σ, π), while Ω qq denotes the quark contribution to the thermodynamic potential. We employ the polynomial form for U [43], with The parameters in the potential are adjusted so as to reproduce the equation of state of the pure gluonic matter with u 3 = 0.75, u 4 = 7.5, a 0 = 6.75, a 1 = −1.95, a 2 = 2.625, a 3 = −7.44, and T 0 = 270 MeV. The mesonic potential is given by where the pion field is suppressed, since we do not consider pion condensation. The explicit chiral symmetry breaking parameter equals h = f π m 2 π .
The quark contribution consists of a vacuum fluctuation part and a purely thermal part where M is an arbitrary renormalization scale and Here the quark mass and energy are given by m q = gσ and E q = p 2 + m 2 , respectively. The Polyakov loop variables Φ and Φ * take real values for real µ q , such that one can pick L = |Φ| andL = |Φ * | as the two independent variables. On the other hand, for imaginary values of the chemical potential they are complex conjugates of each other. Thus, the two independent variables are conveniently chosen as the modulus L and the phase φ, Φ = Le iφ and Φ * = Le −iφ . The expectation values are then determined by the stationarity condition and all thermodynamic quantities can be obtained from the thermodynamic potential Ω(T, µ q ). With N c = 3 and N f = 2, the vacuum parameters are fixed to be f π = 93 MeV, m phys π = 138 MeV, and m σ = 600 MeV, while the Yukawa coupling is set to g = 3.35.
In general, the pseudocritical temperature of a crossover transition is not uniquely determined. By maximizing the chiral (χ σ ) and Polyakov loop (χ L ) susceptibilities, we find, at µ q = 0, the crossover temperatures 231 MeV and 213 MeV for the chiral and deconfinement transitions, respectively. An alternative determination, obtained by maximizing the temperature derivatives of the order parameters, dσ/dT and dL/dT , yields 226 MeV and 223 MeV respectively. In the chiral limit the two procedures yield a unique chiral critical temperature.
In the following, we control the strength of the explicit symmetry breaking by varying the pion mass in vacuum from m π = 0 (chiral limit) to m π = 10 m phys π , in order to assess the effect of criticality on the Fourier coefficients. Thereby, we keep the same vacuum values for f π and m σ by re-adjusting the parameters in U σ . The resulting chiral order parameter at µ B = 0 is shown in Fig. 5. The chiral critical temperature in the chiral limit is found to be 228.256 MeV. For a finite, but small, pion mass (0.1 m phys π ) the pseudocritical temperatures obtained using χ σ and dσ/dT are 229 and 227 MeV, respectively. Hereafter, for simplicity, we denote by T pc the pseudocritical temperature corresponding to a maximum of χ σ . However, in the the heavy pion mass case, the chiral pseudocritical temperature 2 T pc = 228 MeV is determined by the maximum of dσ/dT . The behavior of order parameters and the phase structure in the chiral effective models have been investigated in Refs. [37,45,47]. For the present study we summarize a few relevant features. In Fig. 6, we show contour maps for the temperature derivative of the chiral order parameter d(σ/f π )/dT in the imaginary chemical potential-temperature plane in the chiral limit and for two nonzero values of the pion mass.
In the left panel of Fig. 6, one can unambiguously identify the chiral critical line in the chiral limit, which extends to higher temperature at large θ B , and merges with the Roberge-Weiss transition line at θ B = π and the critical temperature T c (θ B = π) =311 MeV. This firstorder transition line extends from T = ∞ down to the so-called Roberge-Weiss endpoint, which exhibits a second order phase transition at (T = T RW , θ B = π). We find T RW = 308.1 MeV = 1.35 T c which is slightly below T c (θ B = π). At T > T RW , the thermodynamic quantities resemble qualitatively those of a Stefan-Boltzmann gas.
In the case of a non-vanishing pion mass (the middleand the right-panel in Fig. 6), the rapid change of the order parameter found for a small pion mass is smoothed, such that the region around the maximum indicates the location of the crossover transition. Moreover, we find that the curvature of the pseudo-critical line is reduced at small θ B as the pion mass is increased.
In LQCD calculations, it is found that the nature of the Roberge-Weiss endpoint depends on the quark mass [48]. Thus, for small quark masses it is a triple point The imaginary part of the baryon density at imaginary chemical potential in the chiral limit. Solid lines are obtained by differentiating the pressure with respect toμB. The dashed lines represent the density reconstructed from the Fourier series and the dash-dotted lines denote the density reconstructed from the Fourier coefficients obtained using the CEM ansatz.
point. Finally, for large quark masses, the RW endpoint is again a triple point. In the model considered here, the RW endpoint is located at T RW = 311 MeV ≃ 1.35 T pc and 327 MeV ≃ 1.43 T pc for the physical and heavy values of the pion mass, respectively, and is second order in both cases 3 . Given the phase structure in the imaginary chemical potential-temperature plane, shown in Fig. 6, one expects that the Fourier coefficients b k at intermediate tem- peratures, T pc ≤ T ≤ T RW , reflect the existence of the transition line.

B. The Fourier coefficients in the chiral limit
In order to identify the influence of criticality on the Fourier coefficients, we first consider the PQM model formulated in the chiral limit. The expansion coefficients b k are computed by Fourier transforming the density as a function of the imaginary chemical potential, as in Eq. (9). The numerical procedure applied to compute b k was tested on the Stefan-Boltzmann gas, and the results are found to be accurate for Fourier coefficients b k on the order of ∼ 10 −9 or larger. We note that the applicability of the procedure is limited by the magnitude of the coefficients rather than their order k.
In Fig. 7 we show the normalized baryon density Imχ B 1 in the PQM model calculated in the chiral limit at imaginary chemical potential and at several temperatures near T c . The solid lines represent the density obtained directly from the PQM thermodynamic potential by differentiating with respect to the chemical potential, χ B 1 = ∂(p/T 4 )/∂μ B . Also shown are the densities reconstructed from the Fourier coefficients obtained in the PQM model using Eq. (2). Finally, we also show the reconstructed densities obtained by using the CEM ansatz (3) applied to the PQM model.
At T = T c , the densities coincide, as is expected from the smooth behavior of the density. However, at T /T c = 1.2 (blue solid curve), the chiral phase boundary, shown in Fig. 6-left, appears as a kink at θ B /π ≃ 0.88. While the kink is reproduced by the density reconstructed from the Fourier coefficients, the CEM ansatz yields a smooth, analytic function. This is expected, since the CEM Fourier coefficients do not capture the power-law behavior of the higher-order coefficients, which is needed to reproduce the singularity.
At T /T c = 1.4, the density does not vanish at θ B = π. This, together with the fact that the density is an odd function of θ B that exhibits the Roberge-Weiss periodicity, implies that it must be discontinuous at this point. The discontinuity is a manifestation of the first-order Roberge-Weiss transition [38].
Since sin(k θ B ) vanishes at θ B = π for any integer k, the density reconstructed from a finite number of terms in the Fourier series (2) must also vanish at θ B = π. We note that the density reconstructed from the CEM coefficients shows a stronger deviation near θ B = π. Moreover, even the closed form expression (19) for the resummed CEM Fourier series, does not reproduce the discontinuity in the density for λ < 1, i.e., T < ∞. Thus, we expect significant differences between high-order b k and b CEM k at temperatures beyond T RW . In the following, we explore the structure of the Fourier coefficients b k in more detail.
At temperatures significantly below T c , the PQM model exhibits statistical confinement. Consequently there is an oscillatory dependence of the thermodynamic potential on the imaginary baryon chemical potential, with ∆Ω(T, θ) ∼ cos(3θ q ) = cos θ B [45], as in the hadron resonance gas. Therefore the density is, to a good approximation, that of a classical Boltzmann gas, Imχ B 1 ∼ sin θ B with small corrections from quantum statistics and interactions. This, in turn, implies that the higher order Fourier coefficients are strongly suppressed, as in the CEM. Our results confirm this expectation up to T /T c ≤ 0.9 for Fourier coefficients of order k < 9. Owing to the numerical limitations mentioned above, coefficients of order higher order, i.e., k 10, cannot be computed reliably in this temperature range, within the present scheme.
Given the results on the signature of the Roberge-Weiss transition in the Fourier series, one may expect that the chiral critical line at real values of the baryon chemical potential [49][50][51] could also contribute significantly to the high-order Fourier coefficients. However, contributions to the Fourier coefficients from singularities located at real µ B are exponentially suppressed [35], i.e., ∆b k ∼ e −k |ReμB | . This leads to an exponential suppression also of the contribution from the thermal branch points discussed in Sect. II B for nonzero fermion masses, with ∆b k ∼ e −k m .
From general scaling considerations one finds that at T = T c and at large k, mean-field criticality leads to an asymptotic dependence of b k ∼ 1/k 4 [35]. However, since in this case the integration contour, except for the point θ B = 0, lies in the chirally broken phase, a major contribution to the Fourier coefficients is due to massive fermion degrees of freedom. Consequently, the initial k-dependence of the Fourier coefficients is exponential, while the contribution of the critical point θ B is fairly small, implying that the power law is visible only at very large k. Figure 8 shows b k around T = T c = 228.256 MeV. Here we have multiplied b k by k 4 in order to highlight the large k behavior of the Fourier coefficients. One observes a rapid suppression of b k with increasing order below and at T c , in line with a substantial contribution of massive degrees of freedom to the baryon density. Thus, at temperatures T ≤ T c , the density is saturated by the first few Fourier coefficients. Consequently, the CEM ansatz is able to reproduce the baryon density at T = T c , as shown in Fig. 7, although the model does not yield the correct asymptotic behavior of the Fourier coefficients. Moreover, owing to the dominant contribution of massive degrees of freedom, it is extremely difficult, if not impossible, to reliably extract the asymptotic 1/k 4 behavior of the coefficients in numerical calculations.
At a temperature slightly above T c , we find that b k exhibits oscillations, with a 1/k 2 dependence superimposed, as shown in the inset in Fig. 8 for T = 229 MeV. This characteristic dependence on k is expected at temperatures between the critical temperature at µ B = 0, T c , and the Roberge-Weiss temperature T RW . In this temperature range, the critical point is located on the imaginary µ axis (see Fig. 10). It follows that the Fourier coefficients take the asymptotic form A sin(kθ c )/k 2 [35], where θ c = Im µ c /T is the critical value of the baryon chemical potential. Indeed, a fit of the oscillatory dependence on k shown in Fig. 8 yields θ c = 0.421, in perfect agreement with the location of the critical point at T = 229 MeV.
The expectation that the oscillation frequency depends on the location of the singularity is confirmed by the behavior of b k at higher temperatures. In Figure 9 we show the Fourier coefficients at T /T c = 1.1, 1.2, and 1.3. Here we have multiplied b k by (−1) k−1 . The additional phase, changes the frequency of the oscillation from θ c to π − θ c , and thus yields oscillations that are more easily discernible. This is because it is difficult to identify a frequency larger than π/2 on the discrete k grid. The lines show fits with the functional form A sin[k(π − θ c )]. Also shown in Fig. 9 are fits, where a contribution of the regular part of the baryon density, outside the critical region, is included C K 2 (a k)/k + A sin[k(π − θ c )]. Here K 2 (x) denotes the modified Bessel function of the second kind. The values obtained for θ c = 2.15, 2.77, and 3.063 are in agreement with the location of the chiral critical point for the corresponding temperature. For instance, θ c = 2.77 at T /T c = 1.2 yields θ c /π = 0.88, which is the position of the kink in the density in Fig. 7.
The increase of θ c with temperature is a consequence of the shape of the phase boundary (see Fig. 6), which shifts to larger values of θ with increasing temperature. We conclude that the high-order Fourier coefficients carry information on the chiral phase transition at imaginary values of the baryon chemical potential.

C. Fourier coefficients at nonzero pion masses
For nonzero pion masses, the chiral symmetry is explicitly broken, implying that the chiral transition is of the crossover type and that the analytic properties of the  baryon density in the complex chemical potential plane are modified. Specifically, there is a shift of the chiral critical point from the real or imaginary axis into the complex chemical potential plane [50,51], as illustrated in Fig. 10. Owing to the real part of the chemical potential at the singularity, the high-order Fourier coefficients exhibit exponential damping in addition to the powerlaw scaling and oscillatory behavior found in the chiral limit for temperatures above T pc [35]. The resulting Fourier coefficients for temperatures near the crossover temperature T pc are shown in Fig. 11. The upper panel shows the results for a small pion mass, m π = 0.1 m phys π . At T = T pc , we find an oscillatory behavior of k 4 b k , accompanied with a slow damping of the amplitude. A comparison with temperatures slightly above and below T pc , shows that the oscillation frequency increases with temperature, while the damping of the am-   plitude is reduced with temperature. This observation is in qualitative agreement with the temperature dependence of the location of the singularity in the complex chemical potential plane, shown schematically in Fig. 10.
12. Same as Fig. 9, but for the light pion mass.
We note that, given the limited range in k available (k 50), one cannot distinguish between a fit with a function proportional to e −ak /k 2 with large a from one proportional to e −ak /k 4 and small a. Nevertheless, these results clearly demonstrate the sensitivity of the Fourier coefficients to small changes in temperature and pion mass and thus to chiral criticality.
In the lower panel of Fig. 11, we show the Fourier coefficients for a somewhat heavier pion mass, equal to 0.5 m phys π . Here the temperature dependence of the oscillation frequency is reduced. This is expected, since the relevant energy scale, the pion mass, is now large compared to the few MeV change in the temperature. Nevertheless, the qualitative behavior of the oscillations as a function of temperature is still consistent with the expected motion of the singularity in the complex chemical potential plane. Thus, the frequency increases and the damping of the amplitude weakens with increasing temperature.
Unlike the critical temperature in the chiral limit, the pseudocritical temperature T pc is not clearly distinguished by a sudden change in the properties of the Fourier components b k . This is a rather natural consequence of the shift of the singularity into the complex chemical potential plane, indicated in Fig. 10. As the temperature is increased and passes T pc , the motion of the singularity in the complex plane leads to a slow increase of the oscillation frequency and a smooth reduction of the exponential damping of the amplitude.
As in Fig. 9, we show in Fig. 12 the Fourier coefficients multiplied by the factor (−1) k−1 k 2 at T /T pc = 1.1, 1.2, and 1.3. In order to fit the large k dependence of the coefficients b k , we introduce an exponential damping factor of the form e −a k , as well as a regular part. The temperature dependence of the oscillation frequency and the damping factor is again consistent with the expectations deduced from Fig. 10.
A comparison of the results for small pion masses with those obtained in the chiral limit, shows that the higherorder Fourier coefficients are exponentially damped by the explicit breaking of chiral symmetry. Therefore, it is a quantitative question whether chiral criticality can be identified in the Fourier coefficients also for larger pion masses. In order to address this issue, we present the Fourier coefficients obtained in the PQM model for the physical and for a heavy pion mass in Fig. 13.
As illustrated in Fig. 5, in the mean-field approximation the PQM model exhibits a rather rapid chiral crossover transition for a pion mass equal to the physical one, while LQCD yields a smoother transition at the physical point [52,53]. Therefore, a mean-field calculation with a larger pion mass, somewhere between the physical and the heavy pion mass employed here, may be more relevant for a comparison with lattice results.
In the left panel of Fig. 13 we show k 4 b k for physical (squares) and heavy (circles) pion masses at T = T pc . For a physical pion mass, the behavior is qualitatively similar to that shown in the lower panel of Fig. 11.
We note that the sign structure of the Fourier coefficients differs from the alternating signs of the LQCD results (see Fig. 1 and Ref. [39]). As indicated above, the oscillation frequency depends on the location of the contributing singularities in the complex chemical potential plane, in particular on their imaginary parts, as well as on the strength of the individual contributions. Thus, the sign structure of the Fourier components is determined by an interplay between the chiral singularity and, e.g., the thermal branch point. Consequently, the staggered sign structure seen in the LQCD results at T = T pc may be due to non-critical physics, not captured by the PQM model.
The oscillations lead to deviations from a pure exponential damping of the magnitude of the Fourier coefficients, as shown in the inset of the left panel of Fig. 13, while for the heavy pion mass, the oscillations are almost completely washed out by the strong damping.
In the middle and center panels of Fig. 13 we show (−1) k−1 k 2 b k at T /T pc = 1.1 and 1.2. The additional phase factor compared to the left panel, removes the dominant oscillation to a large extent and hence yields clearer plots. At 1.1 T pc , almost all points are at positive values, with a few exceptions for the physical pion mass. Moreover, at 1.2 T pc the sign changes are completely removed by the phase factor. This implies that the staggered sign structure of the LQCD results on the first four coefficients Fourier coefficients [39] is reproduced at temperatures somewhat above the pseudocritical temperature T pc .
As shown in the insets in the middle and right panels, we indeed find an exponential damping of the magnitude of b k for the heavy pion mass for k larger than 3 or 4. At 14. Fourier coefficients near T = TRW for the physical pion mass. In the inset the coefficients multiplied by k 4/3 are shown.
the analysis less clear-cut.

D. The effect of the Roberge-Weiss transition on the Fourier coefficients
Irrespective of the value of the pion mass, the Roberge-Weiss transition appears at T ≥ T RW . It is a firstorder transition, which ends at T = T RW . Assuming that the RW endpoint is a second-order critical point, one finds that the leading contribution to the highorder Fourier coefficients at T = T RW is of the form b k ∼ (−1) k−1 k −(1+1/δ) [35], where δ is the critical exponent associated with the external field strength in the Z(2) universality class. In Fig. 14 we show the Fourier coefficients at T = T RW for the physical value of the pion mass. We find that the coefficients b k indeed follow a power law decay with the exponent 4/3, consistent with δ = 3 in the mean-field approximation. In order to illustrate the sensitivity, we also display b k at T slightly below and slightly above T RW . The inset shows that the deviation from the k −4/3 behavior at the slightly different temperatures is substantial for k > 10. We also note that, in the chiral limit, we find a deviation from the k −4/3 scaling. We attribute this to the contribution from the chiral phase boundary, which is located close to the RW endpoint.
At temperatures above T RW , the density is discontinuous at θ B = π, implying that the Fourier coefficients take the asymptotic form [35] b k ∼ (−1) k−1 /k. This is confirmed by Fig. 15, where we plot |k b k | for T /T pc = 1.5. The asymptotic value of k b k is given by (2/π)Imχ B 1 (T, iπ), and thus directly connected with the discontinuity in the density at the Roberge-Weiss transition [35].
For comparison we also show the Fourier coefficients obtained using the CEM scheme in Fig. 15. Clearly, the 10 100 15. The Fourier coefficients |b k | for a temperature above TRW calculated in the PQM model and in CEM. The dash-dotted horizontal lines indicate the expected asymptotic value, which is determined by the discontinuity in the density.
CEM ansatz does not reproduce the asymptotic behavior of the coefficients related with the RW transition, as anticipated in the discussion of Fig. 7, in section III B. Thus, the CEM ansatz cannot account for the Roberge-Weiss transition, irrespective of the value of the pion mass.
Although these results are obtained in a particular model, the qualitative features are of more general validity, since they are directly linked to the Roberge-Weiss transition.

E. Temperature dependence of low-order coefficients
A crossover or a true phase transition may also be signalled by the temperature dependence of the Fourier coefficients. In Fig. 16 we show the b k , up to k = 5, as functions of temperature for physical and heavy pion masses. For comparison we also show the coefficients obtained using the CEM ansatz.
For the physical pion mass, the chiral crossover transition is reflected already in the lowest order coefficients through the oscillations of b k for temperatures near T pc , as shown in the upper-right panel. On the other hand, for a heavy pion mass this oscillatory behavior is suppressed, owing to the stronger smoothening. The coefficients obtained using the CEM ansatz do not exhibit any oscillatory behavior and shows substantial deviations from the PQM results with increasing temperature.
These results indicate that already at order k ≥ 2, the Fourier expansion coefficients may exhibit a non-trivial temperature dependence in the crossover region. We note that the Fourier coefficients obtained in LQCD do not seem to oscillate as functions of temperature [39]. However, as discussed in Section III C, the interplay between  critical and noncritical physics, and consequently the sign structure of the Fourier coefficients, depends on the location and relative strength of the singularities, which may be different in the model and in QCD. Consequently, we expect that oscillations of the Fourier coefficients in the crossover region will appear also in LQCD calculations, when the location of the chiral singularity approaches the origin of the complex µ-plane, i.e., for pion masses smaller than the physical value (cf. Fig. 10).

F. Reconstructing susceptibilities from b k
The characteristic behavior of the Fourier coefficients b k is reflected also in the baryon number fluctuations. Using Eq. (10), the cumulants at µ B = 0 can be obtained by using the sum Formally the sum should be extended to infinite order but, provided b k decreases sufficiently fast with k, the series converges and the summation may be truncated at some moderate value k max .
We find that for n ≤ 8 the baryon number cumulants χ B 2n computed using Eq. (29) reproduces those obtained directly by taking derivatives of the pressure χ B n = ∂ n (p/T 4 )/∂μ n B for temperatures up to T /T pc ≤ 1.1 − 1.3. For higher temperatures it is numerically difficult to obtain reliable results for the coefficients b k at sufficiently high order to reach agreement. We note that there is an intrinsic problem in calculating higher order cumulants using the Fourier expansion at temperatures above T RW , owing to the slow decrease of the b k with k, which is related to the discontinuity in density at the RW transition.
We first compare the cumulant ratios obtained using the CEM scheme to those computed directly from derivatives of the pressure. In the CEM scheme the summation in Eq. (29) converges, due to the exponential damping    of higher order terms, . Figure 17 displays the χ B 4 /χ B 2 , χ B 6 /χ B 2 , and χ B 8 /χ B 2 , cumulant ratios. As expected from the comparison of Fourier coefficients, these cumulant ratios are rather well reproduced in the CEM scheme for a heavy pion mass, where criticality is strongly suppressed. However, for a physical value of the pion mass, the CEM scheme does not reproduce the χ B 6 /χ B 2 and χ B 8 /χ B 2 cumulant ratios, 4 which are sensitive probes of the chiral crossover transition [14]. Although the CEM scheme roughly reproduces the sign structure of the cumulant ratios, it clearly fails to capture the location of sign changes and the magnitude of the fluctuations.
The reason for these differences between the cumulants of the PQM model and those obtained in the CEM scheme can be traced back to the asymptotic behavior of the Fourier coefficients. In Fig. 18 we show the higher order cumulants computed at T = T pc . Since each derivative with respect to the chemical potential brings one power of k in Eq. (29), the higher order cumulants are more sensitive to the higher order coefficients 5 .
In Fig. 18 we show the χ B n for n = 4, 6 and 8 as functions of the upper limit of the summation k max in Eq. (29). For a very heavy pion mass, there is almost no difference between χ B n computed using the Fourier coefficients of the PQM model and those obtained using the CEM scheme. This is due to the rapid decrease of the higher order coefficients b k : substantial differences between the CEM and the true Fourier coefficients appear only at large k. Since these are small, they do not contribute significantly to the cumulants, as indicated, e.g., by the rapid convergence of χ B 8 , shown in Fig. 18. There are, however, large differences between cumulants of the PQM model and those of the CEM scheme for a physical pion mass. This is due to the higher order coefficients, which in this case contribute substantially to the χ B n . While the CEM cumulants converge rapidly with k, the convergence of the PQM model results is much slower. In the latter, the asymptotic values of χ B 6 and χ B 8 are, at T = T pc , reached for k max = 10 and 15, respectively. The partial sums for the sixth and eighth order cumulants shown in Fig. 18, clearly demonstrate that, for a physical m π , the negative values of these cumulants at T ≃ T pc are by and large due to the large negative values of b 3 and b 4 (see the left panel of Fig. 13). With increasing temperature, the frequency of the oscillation increases, leading to the staggered sign structure and strong cancellations in the sum of Fourier coefficients in Eq. (29). Consequently, at temperatures above T pc , the cumulants χ B 6 and χ B 8 rapidly approach zero, as seen in Fig. 18.

IV. CONCLUDING REMARKS
We have discussed the behavior of the Fourier coefficients of the net-baryon density and their relation to the singularities in the complex chemical potential plane. We found that the presence of singularities is reflected in the asymptotic behavior of the Fourier coefficients. This implies, e.g., that the existence or non-existence of a phase transition cannot be settled by examining only a few loworder terms in the Fourier series.
In order to explore the influence of singularities associated with the chiral phase transition on the Fourier coefficients, we employed the Polyakov-Quark-Meson (PQM) model as a low-energy effective approximation to QCD.
We then showed that the large order behavior of the Fourier coefficients exhibits characteristic oscillations and exponential damping, due to the imaginary and real parts of the chiral singularity, respectively. Consequently, in the chiral limit, the Fourier coefficients at temperatures T ≤ T c show only an exponential decay, while at T > T c they exhibit a power-law decay superimposed on oscillations, whose frequency reflects the imaginary baryon chemical potential at the location of the critical point. As the pion mass is increased, the critical behavior is weakened because the singularity is removed from the imaginary axis and thus gains a larger real part, resulting in a stronger damping of the Fourier coefficients. Moreover, we have shown that both the Roberge-Weiss transition and the corresponding endpoint give rise to characteristic power law decays of the Fourier coefficients.
Based on the characteristics of the Fourier coefficients, we discussed the implications for the higher-order netbaryon number cumulants. We pointed out that the signature of the chiral crossover transition in the net-baryon number cumulants, [14] the sign change around T pc , is related to the oscillatory behavior of Fourier coefficients induced by the chiral branch-point singularity in the complex chemical potential plane.
Our results indicate that the Fourier coefficients provide valuable information on the QCD phase transitions. Thus, lattice QCD calculations of these coefficients will improve our the understanding of the phase structure of QCD. General results on the effect of critical singularities on the Fourier coefficients of the net baryon number will be reported in a separate publication [35].