Chiral crossover in QCD at zero and non-zero chemical potentials

TL;DR

This study uses lattice QCD to determine the pseudo-critical temperatures of QCD chiral crossovers at various chemical potentials, providing precise parametrizations and insights relevant for heavy-ion collision conditions.

Contribution

It offers the first detailed lattice QCD determination of the curvature parameters of the QCD phase boundary at finite chemical potentials.

Findings

Precise value for $T_c(0) = 156.5 \,\mathrm{MeV}$.

Determined curvature coefficients $\,\kappa_2^X$ and $\,\kappa_4^X$ for different chemical potentials.

QCD phase boundary lines are close to chemical freeze-out conditions in heavy-ion collisions.

Abstract

We present results for pseudo-critical temperatures of QCD chiral crossovers at zero and non-zero values of baryon ($B$), strangeness ($S$), electric charge ($Q$), and isospin ($I$) chemical potentials $\mu_{X=B,Q,S,I}$. The results were obtained using lattice QCD calculations carried out with two degenerate up and down dynamical quarks and a dynamical strange quark, with quark masses corresponding to physical values of pion and kaon masses in the continuum limit. By parameterizing pseudo-critical temperatures as $ T_c(\mu_X) = T_c(0) \left[ 1 -\kappa_2^{X}(\mu_{X}/T_c(0))^2 -\kappa_4^{X}(\mu_{X}/T_c(0))^4 \right] $, we determined $\kappa_2^X$ and $\kappa_4^X$ from Taylor expansions of chiral observables in $\mu_X$. We obtained a precise result for $T_c(0)=(156.5\pm1.5)\;\mathrm{MeV}$. For analogous thermal conditions at the chemical freeze-out of relativistic heavy-ion collisions, i.e., $\mu_{S}(T,\mu_{B})$ and $\mu_{Q}(T,\mu_{B})$ fixed from strangeness-neutrality and isospin-imbalance, we found $\kappa_2^B=0.012(4)$ and $\kappa_4^B=0.000(4)$. For $\mu_{B}\lesssim300\;\mathrm{MeV}$, the chemical freeze-out takes place in the vicinity of the QCD phase boundary, which coincides with the lines of constant energy density of $0.42(6)\;\mathrm{GeV/fm}^3$ and constant entropy density of $3.7(5)\;\mathrm{fm}^{-3}$.

Figures

Tables

Table 1. Table 1: Continuum-extrapolated values of second- ( κ 2 X superscript subscript 𝜅 2 𝑋 \kappa_{2}^{X} ) and fourth-order ( κ 4 X superscript subscript 𝜅 4 𝑋 \kappa_{4}^{X} ) Taylor coefficients, defined in Eq. ( 8 ), of pseudo-critical temperature T c ​ ( μ X = B , Q , S , I ) subscript 𝑇 𝑐 subscript 𝜇 𝑋 𝐵 𝑄 𝑆 𝐼 T_{c}(\mu_{X=B,Q,S,I})  obtained from the chiral order parameter Σ ​ ( T , μ X ) Σ 𝑇 subscript 𝜇 𝑋 \Sigma(T,\mu_{X}) and the disconnected chiral susceptibility χ ​ ( T , μ X ) 𝜒 𝑇 subscript 𝜇 𝑋 \chi(T,\mu_{X}) . Also listed are the continuum-extrapolated values of κ 2 B , f superscript subscript 𝜅 2 𝐵 𝑓 \kappa_{2}^{B,f}  and κ 4 B , f superscript subscript 𝜅 4 𝐵 𝑓 \kappa_{4}^{B,f}  for thermal conditions resembling the freeze-out stage of relativistic heavy-ion collisions, i.e. , μ Q ​ ( T , μ B ) subscript 𝜇 𝑄 𝑇 subscript 𝜇 𝐵 \mu_{Q}(T,\mu_{B}) and μ S ​ ( T , μ B ) subscript 𝜇 𝑆 𝑇 subscript 𝜇 𝐵 \mu_{S}(T,\mu_{B}) fixed by strangeness-neutrality and isospin-imbalance of the colliding heavy-ions. The last row is obtained from unweighted average of the first two rows.

	${\kappa }_2^B$	${\kappa }_4^B$	${\kappa }_2^S$	${\kappa }_4^S$	${\kappa }_2^Q$	${\kappa }_4^Q$	${\kappa }_2^I$	${\kappa }_4^I$	${\kappa }_2^{B,f}$	${\kappa }_4^{B,f}$
$\Sigma$	0.015(4)	-0.001(3)	0.018(3)	0.001(3)	0.027(4)	0.004(5)	0.023(3)	0.004(4)	0.012(2)	0.000(2)
$\chi$	0.016(5)	0.002(6)	0.015(4)	0.007(5)	0.031(4)	0.011(9)	0.028(3)	0.006(6)	0.012(3)	0.000(4)
Average	0.016(6)	0.001(7)	0.017(5)	0.004(6)	0.029(6)	0.008(1)	0.026(4)	0.005(7)	0.012(4)	0.000(4)

Equations

$$\Sigma =\frac{1}{f_K^4}\left[m_s⟨\bar{u}u+\bar{d}d⟩-(m_u+m_d)⟨\bar{s}s⟩\right].$$

$${\chi }^{\Sigma }=m_s\left(\frac{\partial }{\partial m_u}+\frac{\partial }{\partial m_d}\right)\Sigma .$$

$$\chi =\frac{m_s^2}{f_K^4}\left[⟨{\left(\bar{u}u+\bar{d}d\right)}^2⟩-{\left(⟨\bar{u}u⟩+⟨\bar{d}d⟩\right)}^2\right].$$

$$\begin{array}{rl}& \Sigma (T,{\mu }_X)=\sum _{n=0}^{\infty }\frac{C_{2n}^{\Sigma }(T)}{(2n)!}{\left(\frac{{\mu }_X}{T}\right)}^{2n},\mathrm{where}\\ & C_{2n}^{\Sigma }(T)={\frac{{\partial }^{2n}\Sigma }{\partial {\left({\mu }_X/T\right)}^{2n}}|}_{{\mu }_X=0}.\end{array}$$

$$\Sigma (T,{\mu }_B)\sim m^{1/\delta }{f}_G;\chi (T,{\mu }_B),{\chi }^{\Sigma }(T)\sim m^{(1-\delta )/\delta }{f}_{\chi }$$

$$\begin{array}{rl}{\partial }_T{\chi }^{\Sigma }(T),{\partial }_T{C}_0^{\chi }(T),& C_2^{\chi }(T)\sim m^{(\beta -\beta \delta -1)/\beta \delta }{f}_{\chi }^{\prime };\\ {\partial }_T{C}_0^{\Sigma }(T),& C_2^{\Sigma }(T)\sim m^{(\beta -1)/\beta \delta }{f}_G^{\prime }.\end{array}$$

$$\begin{array}{rl}& {\partial }_T^2{C}_0^{\Sigma }(T)=0,{\partial }_T{C}_2^{\Sigma }(T)=0,\\ & {\partial }_T{\chi }^{\Sigma }(T)=0,{\partial }_T{C}_0^{\chi }(T)=0,C_2^{\chi }(T)=0,\end{array}$$

$$T_c({\mu }_X)=T_c(0)\left[1-{\kappa }_2^X{\left(\frac{{\mu }_X}{T_c(0)}\right)}^2-{\kappa }_4^X{\left(\frac{{\mu }_X}{T_c(0)}\right)}^4\right].$$

$${{\partial }_T^2\Sigma (T,{\mu }_X)|}_{{\mu }_X}=0,{{\partial }_T\chi (T,{\mu }_X)|}_{{\mu }_X}=0.$$

$$\begin{array}{rl}{\kappa }_2^X& =\frac{1}{2T^2{\partial }_T^2{C}_0^{\chi }}\left[T{\partial }_T{C}_2^{\chi }-2C_2^{\chi }\right],\\ {\kappa }_4^X& =\frac{1}{24T^2{\partial }_T^2{C}_0^{\chi }}[-72{\kappa }_2^X{C}_2^{\chi }-4C_4^{\chi }+T{\partial }_T{C}_4^{\chi }\\ & +12{\kappa }_2^X\left(4T{\partial }_T{C}_2^{\chi }-T^2{\partial }_T^2{C}_2^{\chi }+{\kappa }_2^X{T}^3{\partial }_T^3{C}_0^{\chi }\right)],\end{array}$$

$$T_c(0)=(156.5\pm 1.5)\mathrm{MeV},$$

Full text

1 Introduction

The spontaneous breaking of the chiral symmetry in quantum chromodynamics (QCD) is a key ingredient for explaining the masses of hadrons that constitute almost the entire mass of our visible Universe. Lattice-regularized QCD calculations have demonstrated (near) restoration of the broken chiral symmetry in QCD at high temperature ($T$) through a smooth crossover [1]. The chiral crossover temperature of QCD marks the epoch at which massive hadrons were born during the evolution of the early Universe. The chiral crossover in the early Universe took place at vanishingly small baryon chemical potential $\mu_{B}$, although the electric charge chemical potential $\mu_{Q}$ at that stage might have been non-vanishing [2]. For $\mu_{B}>0$, i.e., when QCD-matter is doped with an excess of quarks over antiquarks, the chiral crossover in QCD might lead to a rich phase diagram in the $T$-$\mu_{B}$ plane [3]. The phase structure of QCD-matter in the $T$-$\mu_{B}$ plane can be probed in various ongoing and upcoming relativistic heavy-ion collision experiments [4]. The phase diagram of QCD can be explored in these experiments if the so-called chemical freeze-out takes place in the proximity of the chiral crossover phase boundary in the $T$-$\mu_{B}$ plane [5]. Since the colliding heavy-ions do not carry any net strangeness, the medium formed in the process is strangeness-neutral, i.e., characterized by $n_{S}=0$, $n_{S}$ being the net strangeness-density. Additionally, the proton-to-neutron ratio of the colliding nuclei determines the ratio of net charge-density ($n_{Q}$) to net baryon-density ($n_{B}$) of the produced medium. For the most common relativistic heavy-ion collisions with Au+Au and Pb+Pb this ratio turns out to be $n_{Q}/n_{B}=0.4$; consequently, the corresponding chemical freeze-out stages also respect the conditions $n_{S}=0$ and $n_{Q}=0.4n_{B}$.

With the aid of state-of-the-art lattice-regularized QCD calculations this work aims at determining chiral pseudo-critical temperatures in QCD at zero and non-zero chemical potentials $\mu_{B,Q,S}$, as well as for the situation analogous to the chemical freeze-out stage of relativistic heavy-ion collision experiments. We will begin by providing the necessary backgrounds in Sec. 2, describe our methods in Sec 3, follow up with our results in Sec. 4, and end with comparisons of our results with extant lattice QCD results and a short summary in Sec. 5.

2 Observables and definitions

2.1 Chiral observables

To define the chiral order parameter we choose the combination

[TABLE]

Here, $\langle\bar{q}{q}\rangle=T(\partial{\ln Z}/\partial{m_{f}})/V$ denotes chiral condensates of the up ($u$), down ($d$), and strange ($s$) quarks; $m_{f}$ denotes the masses of the quarks; $Z$ is the partition function for $2+1$ flavor QCD, with $m_{u}=m_{d}=m_{s}/27$, volume $V$, temperature $T$, and $\langle\cdot\rangle$ denotes average over gauge configurations corresponding to $Z$. The susceptibility corresponding to the chiral order parameter is defined as

[TABLE]

$\chi^{\Sigma}$ contains both quark-line connected, as well as quark-line disconnected pieces. Since the singlet-axial $U_{A}(1)$ symmetry of QCD is expected to remain broken at all $T$, the quark-line connected piece is expected to remain finite even for $m_{u}=m_{d}\to 0$. Thus, we also separately consider the quark-line disconnected chiral susceptibility

[TABLE]

Note that, all chiral observables defined here are free of additive power divergences and renormalization group invariant, ensuring existence of a continuum limit up to small logarithmic corrections in $m_{f}$. Additionally, all chiral observables are defined to be dimensionless in units of the kaon decay constant $f_{K}=156.1/\sqrt{2}$ MeV, the quantity used to determine lattice spacing ($a$) [6].

2.2 Taylor expansions in chemical potentials

The chemical potentials $\mu_{u,d,s}$ of quarks in $Z$ can be traded with the chemical potentials $\mu_{X}$ corresponding to any 3 other linearly independent conserved charges, such as $\mu_{B}$, $\mu_{S}$, $\mu_{Q}$ or $\mu_{I}$. Here, we choose to work with 2 independent sets $\{B,Q,S\}$ and $\{B,I,S\}$. $\mu_{B,Q,S}$ are related to $\mu_{u,d,s}$ through $\mu_{u}=\mu_{B}/3+2\mu_{Q}/3$, $\mu_{d}=\mu_{B}/3-\mu_{Q}/3$, and $\mu_{s}=\mu_{B}/3-\mu_{Q}/3-\mu_{S}$. Similar relations for the $\{B,I,S\}$ set are $\mu_{u}=\mu_{B}/3+\mu_{I}/2$, $\mu_{d}=\mu_{B}/3-\mu_{I}/2$, and $\mu_{s}=\mu_{B}/3-\mu_{S}$.

The $\mu_{X}$ dependence of an observable, e.g., of $\Sigma$, can be obtained by following the well-established Taylor expansion method [7, 8, 9] given by

[TABLE]

For simplicity, here we have assumed all $\mu_{Y\neq X}=0$. Due to $\mathcal{CP}$-symmetry of $Z$, Taylor expansions of the chiral observables contain only even powers of $\mu_{X}$. Similar expansions can be written for $\chi(T,\mu_{X})$, with Taylor coefficients $C^{\chi}_{2n}(T)$. For brevity, we have introduced the notations $C^{\Sigma}_{0}(T)=\Sigma(T,0)$ and $C^{\chi}_{0}(T)=\chi(T,0)$. The detailed expressions for $C^{\Sigma}_{2n}$ and $C^{\chi}_{2n}$ in terms of the $u,d,s$ quark propagators can be found in Refs. [10, 11].

If $\mu_{u,d,s}$ in $Z$ are replaced by $\mu_{B,Q,S}$ and, subsequently, $\mu_{Q}=\mu_{S}=0$ are imposed, then $\mu_{B}$ will be given by the combination $\mu_{B}/3=\mu_{u}=\mu_{d}=\mu_{s}$. Exactly the same will happen for the $\mu_{B,I,S}$ basis if $\mu_{I}=\mu_{S}=0$ conditions are imposed. Similarly, for $\mu_{B}=\mu_{Q}=0$ or $\mu_{B}=\mu_{I}=0$ both bases will lead to $\mu_{S}=-\mu_{s}$, $\mu_{u}=\mu_{d}=0$. Thus, while computing the Taylor coefficients for $B$ and $S$ there is no need to distinguish between $\{B,Q,S\}$ and $\{B,I,S\}$ bases. However, for $\mu_{B}=\mu_{S}=0$, in contrast to the $\mu_{B,I,S}$ basis, Taylor coefficients with respect to $\mu_{Q}$ will receive additional contributions from the strange quark.

2.3 Definitions of pseudo-critical temperatures

The nature of the QCD chiral transition for $m_{u}=m_{d}\to 0$ and $m_{s}>0$ remains an open issue. Nevertheless, increasing numbers of sophisticated lattice QCD calculations are now showing that, in this limit, the QCD chiral transition is most likely a genuine second order phase transition that belongs to the 3D, $O(4)$ universality class [12, 13, 14, 15, 16]. On the other hand, for physical values of the quark masses and vanishing chemical potentials, it is well established that chiral symmetry restoration takes place via a smooth crossover [17, 6, 18]. The present work solely focuses on physical $m_{u,d,s}$. To ascribe precise meaning to chiral crossover temperatures we resort to the well-defined notion of pseudo-critical temperatures $T_{c}(\mu_{X})$.

In the vicinity of the second order chiral phase transition, behaviors of chiral observables are governed by scaling properties of the 3D, $O(4)$ universality class [16, 19]:

[TABLE]

and

[TABLE]

Here, $C^{\Sigma}_{2n}$ and $C^{\chi}_{2n}$ are the coefficients of the Taylor series for $\mu_{B}>0$ and $\mu_{Q}=\mu_{S}=0$. The two relevant scaling functions of the 3D, $O(4)$ universality class, $f_{G}(z)$ and $f_{\chi}(z)$ [20, 21], are functions of the so-called scaling variable $z=t/m^{1/\beta\delta}$, where $m\sim m_{u,d}/m_{s}$, $t\sim(T-T_{c}^{0})/T_{c}^{0}+K(\mu_{B}/T)^{2}$, $T_{c}^{0}$ is $T_{c}(0)$ in the chiral limit $m\to 0$, $\beta$ and $\delta$ are the critical exponents, and $K$ is a non-universal constant.

The chiral critical temperature $T_{c}^{0}$ is defined as the temperature at which $\partial_{T}\Sigma$ and $\chi^{\Sigma}$ diverge in the limit $V\to\infty$ and $m\to 0$. For any $m>0$, residing within the scaling regime, universality dictates that $\partial_{T}\Sigma$ and $\chi^{\Sigma}$, scaled with appropriate (non-integer) powers of $m$, will have maxima located exactly at the maxima of the corresponding scaling functions $f_{G}^{\prime}(z)$ and $f_{\chi}(z)$. Thus, for $m>0$ the locations of the maxima of $f_{G}^{\prime}(z)$ and $f_{\chi}(z)$, denoted by $z^{G}_{p}$ and $z^{\chi}_{p}$, respectively, define two pseudo-critical temperatures $T_{c}^{G,\chi}(0)$. As $m\to 0$, $\partial_{T}\Sigma$ and $\chi$ diverge, and $T_{c}^{G,\chi}(0)$ reduce to $T_{c}^{0}$ according to the scaling relation $T_{c}^{G,\chi}(0)=T_{c}^{0}+Az^{G,\chi}_{p}m^{1/\beta\delta}$, with a non-universal constant $A$.

Physical values of $m_{u,d}$ might not reside within the scaling regime of the second order chiral phase transition; consequently, chiral observables may also contain additional non-singular, polynomial in $m$, corrections. Thus, for physical values of $m_{u,d}$ we define $T_{c}(0)$ using the following criteria

[TABLE]

where $C^{\Sigma}_{2n}$ and $C^{\chi}_{2n}$ are the coefficients of the Taylor series for $\mu_{B}>0,~{}\mu_{Q}=\mu_{S}=0$. Each of these 5 criteria may lead to 5 different values of $T_{c}(0)$, all of which will reduce to the unique $T_{c}^{0}$ as $m\to 0$. If the physical values of $m_{u,d}$ happen to be in the scaling regime, then all these 5 criteria will lead to only two values of pseudo-critical temperatures $T_{c}^{G,\chi}(0)$. The above definitions exhaust all second order fluctuations of the chiral order parameter through which locations of the maxima of $f_{G}^{\prime}$ and $f_{\chi}$ can be determined.

Following the spirit of Taylor expansions, $\mu_{X}$ dependence of pseudo-critical temperatures, up to $\mathcal{O}(\mu_{X}^{4})$, can be written as

[TABLE]

As we will see later in Sec. 4.1, for $\mu_{X}=0$, the $T_{c}(0)$ defined through all 5 criteria listed in Eq. (7) actually lead to the same result, within our errors, in the continuum limit. Thus, for $\mu_{X}>0$ it is sufficient to define $T_{c}(\mu_{X})$ by the 2 criteria

[TABLE]

Expressions for $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$ can be obtained by: (i) Expanding $\Sigma(T,\mu_{X})$, $\chi(T,\mu_{X})$ in $\mu_{X}$ using Eq. (4); (ii) Taylor expanding $C^{\Sigma}_{2n}$, $C^{\chi}_{2n}$ in powers of $(T_{c}(\mu_{X})-T_{c}(0))$; (iii) Expanding $(T_{c}(\mu_{X})-T_{c}(0))$ using Eq. (8), keeping terms up to $\mathcal{O}(\mu_{X}^{4})$; (iv) Taking $\partial_{T}$ at fixed $\mu_{X}$ of the fully-expanded expression up to $\mathcal{O}(\mu_{X}^{4})$, and imposing Eq. (9) order-by-order in $\mu_{B}$. Since all quantities are assumed to be analytic in $\mu_{X}$ around $\mu_{X}=0$, all expansions in $\mu_{X}$ and taking $\partial_{T}$ can be carried out in any order, as long as all terms contributing up to $\mathcal{O}(\mu_{X}^{4})$ are systematically included at each step. E.g., for $\chi$ we obtained [10]

[TABLE]

where $C^{\chi}_{2n}$  are the expansion coefficients of $\chi$ with respect to $\mu_{X}$, and the expressions are to be evaluated at $T=T_{c}(0)$. Similar expressions can be obtained for $\Sigma$ [10].

The expression for our $\kappa_{2}^{B}$ corresponding to the order parameter is different from that used in Ref. [22], where $T_{c}(\mu_{B})$ was defined through temperature derivatives at constant $\mu_{B}/T$, rather than at constant $\mu_{B}$. We have checked that the numerical results using both definitions are same within our errors.

3 Computational details

All computations presented in this study were carried out with the lattice actions previously used by the HotQCD collaboration [6, 23, 24], viz., the $2+1$ flavor highly improved staggered quarks (HISQ) [25] and the tree-level improved Symanzik gauge action. The bare parameters of the lattice actions, $m_{u}=m_{d}$, $m_{s}$, and the bare gauge coupling, are fixed by the line of constant physics determined by the HotQCD collaboration [6, 23, 24]. The temperature is given by $T=1/(aN_{\tau})$, where $N_{\tau}$ is the extent of the lattices along the Euclidean temporal direction. The extents of the lattices along all 3 spatial directions were always chosen to be 4$N_{\tau}$, and the temporal extents were varied from $N_{\tau}$=6, 8, 12, and 16, going towards progressively finer lattice spacing at a fixed $T$. Bare quark masses were chosen to reproduce, within a few percent, the physical value of the kaon mass and a pseudo-Goldstone pion mass of 138 MeV in the continuum limit at vanishing temperature and chemical potentials.

The fermionic operators needed to construct $C^{\Sigma}_{4}$ and $C^{\chi}_{2,4}$ were obtained using the so-called linear-$\mu$ formalism [26, 27, 24], but the traditional exponential-$\mu$ formalism [28] was used for $C^{\Sigma}_{2}$. On dimensional grounds, within the linear-$\mu$ formalism no additive ultraviolet divergence (or constant) is expected in $C^{\chi}_{2n}$ for all $n$, and in $C^{\Sigma}_{2n}$ for $n>1$ [10]. To confirm these theoretical expectations we computed $C^{\chi}_{2}$ by employing both linear- and exponential-$\mu$ formalism and found identical results for both cases [10]. The $n^{\mathrm{th}}$ order Taylor coefficients of the chiral observables contain up to $n+1$ quark propagators, compared to $n$ quark propagators for that in case of the pressure ($TV^{-1}\ln Z$). Hence, the computational cost of $C^{\Sigma}_{2n}$ and $C^{\chi}_{2n}$ increases accordingly.

All fermionic operators needed to construct the chiral observables and their Taylor coefficients were measured on about 100K, 500K, 100K and 4K gauge field configurations for $N_{\tau}$=6, 8, 12 and 16 lattices, respectively. In each case, the gauge field configurations were separated by 10 rational hybrid Monte-Carlo trajectories of unit length. The fermionic operators were calculated using the standard stochastic estimator technique; more details about these computations can be found in Ref. [10].

As discussed in Sec. 2.3, determinations of $T_{c}(0)$, $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$ involve computing derivatives of the basic chiral observables and their Taylor coefficients with respect to the temperature. To compute these derivatives, we interpolated the basic observables in $T$ between the computed data via the following procedure. For each observable several $[m,n]$ Padé approximants were used for $N$ ($>m+n$) computed data, and $N$ was varied by leaving out data away from the crossover region. Statistical error of each Padé approximant was estimated using the bootstrap method; the bootstrap samples for each computed data were drawn from a Gaussian distribution centered around the mean value of the data and with a standard deviation equal to the $1\sigma$ statistical error of that data. The final $T$-interpolation for each observable was obtained by weighted averaging over all the Padé approximants where the weight for an approximant was determined using the Akaike information criterion [29, 30]. This procedure gave reliable results for all the required $T$-derivatives, especially for $T$ in the vicinity of the chiral-crossover [10].

We assumed that for all observables the leading discretization errors are of the type $a^{2}\propto 1/N_{\tau}^{2}$. Extrapolations to the continuum limit $a\to 0$ were carried out by fitting data at different $N_{\tau}$ to a function linear in $1/N_{\tau}^{2}$ and extrapolating it to $N_{\tau}\to\infty$ limit. The error on each continuum-extrapolated result was obtained using the above described bootstrap method. For all observables we found that $1/N_{\tau}^{2}$-fits were satisfactory. To check the systematics of our continuum extrapolations, we used fits including higher order $1/N_{\tau}^{4}$ corrections, as well as carried out the extrapolation procedure using an alternative $T$-scale determined using the Sommer parameter $r_{1}$; all results were found to be consistent within our errors [10].

4 Results

4.1 Zero chemical potential: $T_{c}(0)$

In Figs. 1 and 3, we show all observables used for the determination of pseudo-critical temperatures as defined in Eq. (7) for lattices with $N_{\tau}$=6, 8, 12, and 16. The results of the temperature interpolations, obtained following the procedure described in Sec. 3, are shown by the corresponding solid bands. Using the interpolated results, and applying the definitions in Eq. (7), we obtained 5 values of $T_{c}(0)$ for $N_{\tau}$=6, 8, and 12. These results are shown in Fig. 2. Since we have not computed $C^{\Sigma}_{2}$ and $C^{\chi}_{2}$ for $N_{\tau}$=16, we only show results for the other 3 definitions of $T_{c}(0)$. On coarser lattices, different definitions resulted in different values of $T_{c}(0)$. These differences progressively reduce with increasingly finer lattice spacing. Results of $T_{c}(0)$ for each of the definitions were separately extrapolated to the continuum (see Sec. 3 for details). The continuum-extrapolated results for all 5 definitions of $T_{c}(0)$ were all consistent with each other within errors. We took an unweighted average of all the 5 continuum results, and added the statistical errors of each continuum-extrapolation in quadrature to quote our final result for the chiral crossover temperature at zero chemical potentials $T_{c}(0)=(156.5\pm 1.5)$ MeV. It is an interesting fact that continuum results for different pseudo-critical temperatures coincide within a couple of MeV. However, if the value of $T_{c}^{0}$ [12] is significantly different from $T_{c}(0)$, then, based on the scaling properties of $T_{c}^{G,\chi}(0)$, it is natural to expect more dispersion among the values of $T_{c}(0)$. Coincidence of different pseudo-critical temperatures for physical quark masses may accidentally arise due to the presence of non-singular and/or sub-leading corrections to scaling. Further work will be needed to clarify this issue.

4.2 Non-zero chemical potentials: $\kappa_{2}^{B,Q,S,I}$ and $\kappa_{4}^{B,Q,S,I}$

Now, we present continuum-extrapolated results for the expansion coefficients $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$, defined by Eq. (8), of $T_{c}(\mu_{X})$ for all conserved charges $X=B,S,Q,I$. In all cases, extrapolations to the continuum were carried out using results for $N_{\tau}$=6, 8, and 12. We discuss an example in detail, viz., $\kappa_{2}^{B}$ and $\kappa_{4}^{B}$ at $\mu_{Q}$=$\mu_{S}$=0. When $T_{c}(\mu_{B})$ is defined as the temperature where $\chi(T,\mu_{B})$ peaks at a given $\mu_{B}$, the corresponding $\kappa_{2}^{B}$ and $\kappa_{4}^{B}$ can be obtained using Eq. (10). The zeroth-, $C^{\chi}_{0}(T)$, second-, $C^{\chi}_{2}(T)$, and the fourth-order, $C^{\chi}_{4}(T)$, expansion coefficients of $\chi(T,\mu_{B})$ in $\mu_{B}/T$ (with $\mu_{Q}$=$\mu_{S}$=0) are shown in Fig. 1 (middle), Fig. 3 (top-left) and Fig. 3 (top-middle), respectively. The interpolations in $T$ are shown by the corresponding solid bands. Having determined $T_{c}(0)$, $\kappa_{2}^{B}$ and, subsequently, $\kappa_{4}^{B}$ were obtained using the $T$-interpolations of $C^{\chi}_{0,2,4}$. Similarly, $\kappa_{2}^{B}$ and $\kappa_{4}^{B}$ were computed also from the inflection point of $\Sigma(T,\mu_{B})$ in $T$, for a given $\mu_{B}$, using the expansion coefficients $C^{\Sigma}_{0,2,4}$, which are shown in Fig. 1 (left), Fig. 3 (bottom-left) and Fig. 3 (bottom-middle), respectively. Fig. 3 (bottom-right) exemplifies the very mild dependence of $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$ on lattice spacing.

We also carried out similar computations to determine continuum-extrapolated $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$ corresponding to (i) $T_{c}(\mu_{S})$ at $\mu_{B}$=$\mu_{Q}$=0, (ii) $T_{c}(\mu_{Q})$ at $\mu_{B}$=$\mu_{S}$=0, and (iii) $T_{c}(\mu_{I})$ at $\mu_{B}$=$\mu_{S}$=0; the values are listed in Tab. 1. In all the cases, for both $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$, the results obtained using two different definitions of $T_{c}(\mu_{X})$, given in Eq. (9), gave the same result within our errors. In each case, we took unweighted averages of continuum-extrapolated results corresponding to both definitions for $T_{c}(\mu_{X})$, and added the respective statistical errors in quadrature to arrive at the final values for $\kappa_{2}^{X}$ and $\kappa_{4}^{X}$; these final results also are listed in the third row of Tab. 1. In all cases, $\kappa_{4}^{X}$ were found to be zero within errors, with central values about an order of magnitude smaller than the corresponding $\kappa_{2}^{X}$. Also, $\kappa_{2}^{Q,I}$ were found to be about a factor 2 larger compared to $\kappa_{2}^{B,S}$.

4.3 Heavy-ion collisions: $\kappa_{2,4}^{B,f}$ for $n_{S}=0$, $n_{Q}=0.4n_{B}$

In this case, i.e., for the thermal condition resembling the chemical freeze-out stage of heavy-ion collision experiments, we introduce the notations $\kappa_{n}^{B,f}$ as the Taylor coefficients of the corresponding pseudo-critical temperature $T_{c}^{f}(\mu_{B})$.

The formalism for Taylor expanding an observable in $\mu_{B}/T$, with the constraints $n_{S}=0$ and $n_{Q}=0.4n_{B}$, was introduced in Ref. [31] and has been applied to various cases [32, 24, 33]. With these constraints, $\mu_{S}$ and $\mu_{Q}$ are no longer arbitrary, but become functions of $T$ and $\mu_{B}$. Following Ref. [31], $\mu_{S}(T,\mu_{B})/T=s_{1}(T)\mu_{B}/T+s_{3}(\mu_{B}/T)^{3}$ and $\mu_{Q}(T,\mu_{B})/T=q_{1}(T)\mu_{B}/T+q_{3}(\mu_{B}/T)^{3}$ were Taylor-expanded in $\mu_{B}/T$. Expanding $n_{B,Q,S}$ in powers of $\mu_{B}^{i}\mu_{Q}^{j}\mu_{S}^{k}$ ($i+j+k\leq 3$), substituting expansions for $\mu_{Q,S}(T,\mu_{B})$ in expansions of $n_{B,Q,S}$, and imposing the constraints $n_{S}=0$ and $n_{Q}=0.4n_{B}$ order-by-order in $\mu_{B}$, expressions for $s_{1,3}(T)$ and $q_{1,3}(T)$ were obtained in terms of the Taylor coefficients of the pressure. Explicit expressions for $s_{1,3}(T)$ and $q_{1,3}(T)$ can be found in Ref. [24]. By Taylor expanding $\Sigma(T,\mu_{B},\mu_{Q},\mu_{S})$ ($\chi(T,\mu_{B},\mu_{Q},\mu_{S})$) in powers of $\mu_{B}^{i}\mu_{Q}^{j}\mu_{S}^{k}$ ($i+j+k\leq 4$) and by using the expansions for $\mu_{Q,S}(T,\mu_{B})$, we obtained the expansions for $\Sigma(T,\mu_{B})$ ($\chi(T,\mu_{B})$) up to $\mathcal{O}(\mu_{B}^{4})$. As before, by invoking Eq. (9), expressions were obtained for $\kappa_{2,4}^{B,f}$.

Continuum-extrapolated results for $\kappa_{2}^{B,f}$ and $\kappa_{4}^{B,f}$ are given in Tab. 1. $\kappa_{2}^{B,f}$ came out to be same as $\kappa_{2}^{B}$ and $\kappa_{2}^{S}$ within errors, and $\kappa_{4}^{B,f}$ was found to be consistent with zero. On our $N_{\tau}$=8 lattices, where we analyzed half a million gauge configurations at all $T$, we also computed $\mu_{B}^{6}$ corrections to the chiral observables. The order-by-order $\mu_{B}$ corrections to $\chi$ are shown in Fig. 3 (top-right) at $\mu_{B}$=300 MeV and for $n_{S}=0$, $n_{Q}=0.4n_{B}$. In the vicinity of $T_{c}^{f}(\mu_{B})$, difference between $\mu_{B}^{4}$ and $\mu_{B}^{2}$ corrections are clearly significant; but $\mu_{B}^{6}$ and $\mu_{B}^{4}$ corrections are consistent within our errors. This shows that up to $\mu_{B}^{4}$ the expansion of $T_{c}^{f}(\mu_{B})$ is controlled till $\mu_{B}\lesssim 2T_{c}(0)$. The phase boundary of QCD for $n_{S}=0$, $n_{Q}=0.4n_{B}$ is shown in Fig. 4; also shown are the chemical freeze-out points extracted from heavy-ion collision experiments at various collision energies [5, 34], the line of constant energy density $\epsilon(T,\mu_{B})=\epsilon(T_{c}(0),0)=0.42(6)~{}\mathrm{GeV/fm}^{3}$ [24], and the line of constant entropy density $s(T,\mu_{B})=s(T_{c}(0),0)=3.7(5)~{}\mathrm{fm}^{-3}$ [24].

5 Discussions and summary

The value of $T_{c}(0)$ reported in this work compares quite well with the previous results from the HotQCD collaborations [6, 17], but the present result is about 6 times more accurate than the previous continuum-extrapolated result [6]. Compared to that of Ref. [6], use of 100-500 times more gauge configurations for $N_{\tau}=6,8,12$ in the present study resulted in the 6 times more accurate determination of the continuum-extrapolated $T_{c}(0)$. Our present value of $T_{c}(0)$  also is compatible with the chiral pseudo-critical temperatures reported by other groups [35, 36]. It is pertinent to note that all our calculations were carried out within a finite-size box of about 5 fm3 in the vicinity of $T_{c}(0)$; finite-size corrections might increase the value of $T_{c}(0)$ by an amount commensurate to our present error on that quantity [37]. $\kappa_{2}^{B}$ determined in the present work is about a factor 2 larger than that reported previously in Ref. [38].

Our present value of $\kappa_{2}^{B}$ also is about a factor 2 larger than the $\kappa_{2}^{B}$ estimated using the curvature of the chiral critical temperature along the light quark chemical potential directions [19], but is consistent, within errors, with the same reported in Ref. [39]. In contrast to Ref. [19], Ref. [39] used the much improved HISQ discretization. This clearly suggests that the discrepancy between the present result and that estimated from Ref. [19] arises mostly due to the use of improved HISQ discretization in the present study. On the other hand, $\kappa_{2}^{B}$ reported in this work is, within errors, compatible with those obtained in more recent works of Refs. [40, 41, 36, 42], obtained from analytic continuations from purely imaginary $\mu_{B}$. It is also similar with that obtained in Ref. [22] from Taylor expansion of chiral order parameter for $\mu_{B}>0$, $\mu_{Q}$=0 and $\mu_{S}$=$\mu_{B}$/3, in contrast to our choice of $\mu_{B}>0$ and $\mu_{Q}$=$\mu_{S}$=0. Our value of $\kappa_{2}^{B,f}$ is quite similar to that reported in Ref. [43], determined from analytic continuations from purely imaginary $\mu$. Moreover, the phase boundary in the $T$-$\mu_{I}$ plane that can be obtained using our $\kappa_{2,4}^{I}$ is quite similar to that determined in Ref. [44] from lattice QCD computations performed directly at $\mu_{I}>0$, $\mu_{B}$=$\mu_{S}$=0.

In summary, using state-of-the-art lattice QCD computations we have determined pseudo-critical temperatures, $T_{c}(\mu_{X})=T_{c}(0)[1-\kappa_{2}^{X}(\mu_{X}/T_{c}(0))^{2}-\kappa_{4}^{X}(\mu_{X}/T_{c}(0))^{4}]$, of QCD chiral crossover for 6 different scenarios: (i) $T_{c}(0)$ for $\mu_{B}=\mu_{Q}=\mu_{S}=0$; (ii) $\kappa_{2,4}^{B}$ for $\mu_{B}>0$, $\mu_{Q}$=$\mu_{S}$=0; (iii) $\kappa_{2,4}^{S}$ for $\mu_{S}>0$, $\mu_{B}$=$\mu_{Q}$=0; (iv) $\kappa_{2,4}^{Q}$ for $\mu_{Q}>0$, $\mu_{B}$=$\mu_{S}$=0; (v) $\kappa_{2,4}^{I}$ for $\mu_{I}>0$, $\mu_{B}$=$\mu_{S}$=0; (vi) $\kappa_{2,4}^{B,f}$ for thermal conditions resembling that at the chemical freeze-out of relativistic heavy-ion collision experiments, viz, for $\mu_{B}>0$, $n_{S}=0$, $n_{Q}=0.4n_{B}$. We have found

[TABLE]

and the values of $\kappa_{2,4}^{X}$ are listed in Tab. 1. The QCD phase boundary relevant for relativistic heavy-ion collision experiments have been summarized in Fig. 4. For $\mu_{B}\lesssim 300$ MeV, the chemical freeze-out takes place close to the QCD chiral crossover, which, in turn, seems to happen along lines of constant energy density of $0.42(6)~{}\mathrm{GeV/fm}^{3}$ and a constant entropy density of $3.7(5)~{}\mathrm{fm}^{-3}$. At vanishing baryon chemical potential $\mu_{B}$, the ALICE result [5] for the chemical freeze-out temperature is in agreement with $T_{c}(0)$. For $\mu_{B}\lesssim 300$ MeV, all STAR results [34], except the highest collision-energy, agree with $T_{c}(\mu_{B})$ within their 1-sigma errors. The STAR result for the chemical freeze-out temperature at the highest collision-energy agrees with $T_{c}(\mu_{B})$ within 1.5-sigma error. Thus, there is no discrepancy between $T_{c}(\mu_{B})$ and chemical freeze-out temperatures extracted using statistical model based fits to the experimentally measured hadron yields. However, it may pose a challenge to the statistical hadronization based chemical freeze-out scenario if future improved experiments determine freeze-out temperatures with statistical significance above $T_{c}(\mu_{B})$.

Acknowledgments

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics: (i) Through the Contract No. DE-SC0012704; (ii) Within the framework of the Beam Energy Scan Theory (BEST) Topical Collaboration; (iii) Through the Scientific Discovery through Advance Computing (ScIDAC) award Computing the Properties of Matter with Leadership Computing Resources.

This research also was funded by— (i) The Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project number 315477589-TRR 211; (ii) The grant 05P18PBCA1 of the German Bundesministerium für Bildung und Forschung; (iii) The National Natural Science Foundation of China under grant numbers 11775096 and 11535012 (HTD); (iv) The Early Career Research Award of the Science and Engineering Research Board of the Government of India (PH); (v) Ramanujan Fellowship of the Department of Science and Technology, Government of India (SS).

This research used awards of computer time provided by the INCITE and ALCC programs at: (i) Oak Ridge Leadership Computing Facility, a DOE Office of Science User Facility operated under Contract No. DE-AC05-00OR22725; (ii) National Energy Research Scientific Computing Center, a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231; (iii) Argonne Leadership Computing Facility, a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-06CH11357.

This research also used computing resources made available through: (i) The USQCD resources at BNL, FNAL and JLAB; (ii) The PRACE grants at CSCS, Switzerland, and CINECA, Italy; (iii) The Gauss Center at NIC-Jülich, Germany; (iv) Nuclear Science Computing Center at Central China Normal University.