Chiral phase transition of QCD with $N_f=2+1$ flavors from holography

Chiral phase transition for three-flavor $N_f=2+1$ QCD with $m_u=m_d\neq m_s$ is investigated in a modified soft-wall holographic QCD model. Solving temperature dependent chiral condensates from equations of motion of the modified soft-wall model, we extract the quark mass dependence of the order of chiral phase transition in the case of $N_f=2+1$, and the result is in agreement with the"Colombia Plot", which is summarized from lattice simulations and other non-perturbative methods. First order phase transition is observed around the three flavor chiral limit $m_{u/d}=0, m_{s}=0$, while at sufficient large quark masses it turns to be a crossover phase transition. The first order and crossover regions are separated by a second order phase transition line. The second order line is divided into two parts by the $m_{u/d}=m_s$ line, and the $m_s$ dependence of the transition temperature in these two parts are totally contrast, which might indicate that the two parts are governed by different universality classes.


Introduction
Quantum Chromodynamics(QCD) is widely accepted as the fundamental theory of strong interactions, and QCD vacuum is characterized by spontaneous chiral symmetry breaking together with color charge confinement. The dynamically generated chiral condensate ψ ψ in the vacuum serves as quarks' dynamical mass and spontaneously breaks the chiral symmetry, which is an exact symmetry of QCD lagrangian when quarks are massless. It is believed that at sufficient high temperature and/or density, quark condensate might be destroyed completely and the spontaneous breaking symmetry would be restored. Understanding the property of chiral phase transition has been an important topic in both non-perturbative QCD and cosmology for decades [1]. The order of chiral phase transition depends sensitively on the inertial degrees of freedom of the system, such as the number of flavors(N f ) and the mass of quarks(m u , m d and m s ). Based on theoretical consideration and lattice QCD simulations [2][3][4], the expected three flavor phase diagram in the quark mass plane, describing quark mass dependence of the order of QCD phase transitions, is summarized in the sketch plot (it is also called "Columbia Plot") shown in Fig.1(a). In this sketch plot, the whole m u/d − m s plane are divided into three simply connected region: two first order zones in the bottom left corner and the upper right corner as well as the crossover region in the middle. The upper right corner is near the infinite quark mass limit m u = m d = m s = ∞, where the breaking and restoration of Z 3 centre symmetry, related to confinement/deconfinement phase transition, are well defined. The bottom left corner is near the chiral limit m u = m d = m s = 0, where the breaking and restoration of chiral symmetry, related to chiral phase transition, are well defined. In between these two regions, there are no known exact symmetries and the phase transitions are expected to be a continual and rapid transition between different phases, usually named "crossover". The boundary of these three regions are second order lines, where the phase transitions are expected to be of second order. Furthermore, we noted that the second order line between the bottom left first order region and the crossover region is divided by the SU (3) diagonal line with m u = m d = m s into two parts, the upper one of which is expected to be governed by the O(4) universality class [6] while the lower one of which is expected to be governed by the Z(2) universality class. Theoretically, the dynamics of QCD near phase transition is non-perturbative, and normal perturbative methods of quantum field theory become invalid here. Lattice QCD has been considered as the most reliable non-perturbative method to study non-perturbative properties of QCD. However, lattice QCD simulations still require further improvements in several aspects, especially the difficulty called sign problem at finite chemical potential, which is waited to be solved in order to extract full understanding on QCD phase diagram. Hence, it is quite necessary to develop other non-perturbative methods. The recent progress of anti-de Sitter/conformal field theory (AdS/CFT) correspondence and the conjecture of the gravity/gauge duality [7][8][9] does provide such a new powerful tool to tackle the strong coupling problem of gauge theory like QCD.
In the framework of holography, by breaking the conformal symmetry of the original AdS/CFT correspondence in different ways, efforts towards realistic holographical description of non-perturbative physics of QCD, such as hadron physics  and hot/dense QCD matter [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50], have been made both in top-down approaches and in bottom-up approaches(see Refs. [51][52][53][54][55] for reviews). Different from top-down approach, bottom-up holographic QCD starts from QCD phenomena and try to build up more realistic holographic models. In this approach, chiral phase transition has been studied in several different models [5,[56][57][58][59][60][61][62][63][64][65]. Most of these studies considered the case with equal quark masses for all quarks. As can be seen from Fig.1(a), it is also interesting to consider the cases when m u/d = m s , where the physical point locates. In our previous studies [5,63], based on soft-wall AdS/QCD model [11], the mass dependence of chiral phase transition is extracted, as shown in Fig.1(b). The qualitative results for SU (2) and SU (3) cases are in good agreement with the current understanding from Fig.1(a): for two flavor case it starts from a second order phase transition and turns to be crossover at any finite quark mass while for three flavor case it starts from a first order phase transition and only at sufficient large quark masses it turns to be crossover. Furthermore, in [65], we extend these studies to finite magnetic field, and find that it can provide good description on inverse magnetic catalysis effect, which was discovered in lattice QCD [66,67] and studied in other methods [68][69][70][71][72][73][74][75][76][77][78] recently. However, the cases when m u/d = m s have not been examined in this model. Therefore, in this work, we will extend our studies in [5,63] to N f = 2 + 1 when m l ≡ m u = m d = m s and study the property of chiral phase transition.
The paper is organized as follows. In Sec.2, we give a short introduction on the model and the numerical method we used, especially on how to introduce quark masses and chiral condensates in the model. Then in Sec.3 we show numerical results from our model study, especially the quark mass dependence of the order of the chiral phase transition like Fig.1. Finally in Sec.4 a brief discussion will be given.
2 Soft-wall model in N f = 2 + 1 case

Background
In the original paper of soft-wall model [11], the 4D global chiral symmetry SU (N f ) L × SU (N f ) R is promoted to 5D and becomes local gauge symmetry of the following action with Φ the dilaton field, X a complex scalar field, V X the scalar potential, F mn the field strength defined as F ,n ] in terms of the left/right hand gauge potential A L/R , g 5 the gauge coupling, g the determinant of metric g mn , and the covariant derivative D m defined as D m X = ∂ m X −iA L m X +iXA R m . The scalar potential V X only takes the mass term and has the form of From the AdS/CFT prescription M 2 5 = (∆ − p)(∆ + p − 4) [9], the mass of the complex scalar field X M 2 5 can be determined as M 2 5 = −3 L 2 (we will take the AdS radius L = 1 in this work) by taking ∆ = 3, p = 0. The dilaton field is taken to be a simple quadratic form Φ(z) = µ 2 z 2 . In this way, the meson spectra are shown to be linear with respect to the radial excitation quantum number n at large n, which gives good description of the linear behavior of meson spectra. However, in the original soft-wall model, there is no spontaneous chiral symmetry breaking in QCD vacuum and also no restoration at sufficient high temperature.
As pointed out in [5,63], further modifications of the dilaton field and the scalar potential are necessary in order to describe the spontaneous chiral symmetry breaking in QCD vacuum and its restoration. The specific profile of the dilaton field are proposed [5,63] and it takes the following form The positive quadratic behavior of Φ(z) is responsible for the linear spectra, which is well known in the soft-wall AdS/QCD. The scalar potential takes the first several leading powers of V X and has the form of In [5,63], we have shown that the negative part as well as the quartic term λ|X| 4 in the scalar potential are essential for the spontaneous chiral symmetry breaking in the vacuum as well as its restoration at sufficient high temperature. In SU (2) or two-flavor case, the t'Hooft determinant term γRe[det(X)] is taken to be zero, and we find a second order phase transition in the chiral limit and a crossover transition at any finite quark mass case. In SU (3) or three-flavor case, we consider the t'Hooft limit, and we find a first order chiral phase transition in the chiral limit, while only at sufficient quark mass the phase transition turns to be a crossover one. Furthermore, in our recent study [65], we show that after introducing magnetic field through the Einstein-Maxwell sector, the above SU (2) model can describe inverse magnetic catalysis in the soft-wall model quite well. All the qualitative results are in agreement with the current understanding from lattice simulations [2-4, 66, 67] and other non-perturbative studies. Nevertheless, in our previous study, we have only considered equal mass in both the SU (2) and SU (3) cases, i.e., m u = m d and m u = m d = m s . Therefore, we can only obtain the top line and the diagonal line in Fig.1. As shown in Fig.1, it is also interesting to consider the case N f = 2 + 1 with m l ≡ m u = m d = m s case. It would be a quite natural extension of our previous study in SU (3) case to the N f = 2 + 1 case. In the following, A L , A R will be set to be zero, since only the scalar field X is relevant for the chiral phase transition. Since m l ≡ m u = m d = m s , the complex scalar field should take the following form instead of a simple χI 3 with I 3 the 3 × 3 matrix. Here, we assume that χ l (z), χ s (z) are functions of the 5D coordinate z only and the factor 1 √ 2 is just a normalization constant 1 . When m l = m s , we should have χ l = χ s , since the boundary values of χ l , χ s are related to the quark masses. Under this ansatz, we can get the effective form of the action Eq.2.1 as and L the AdS radius, which will not affect the final results and will be taken to be 1 later.
As in [5,63], we will neglect the back-reaction of χ l , χ s to the background metric, and take the simple AdS-Schwarzchild(AdS-SW) black hole solutions where z h is the black hole horizon defined at f (z h ) = 0 and could be related to the temperature T by the following relation Under this ansatz, the equations of motion for χ l , χ s could be derived as the following form Please notice that in the SU (3) case with equal quark mass m l (≡ m u = m d ) = m s , we can have χ l = χ s ≡ χ, and the above two equations will be reduced to the same one which is exactly the same as the one in [5,63]. In [5,63], we have taken v 3 = −3, v 4 = 8 to show the qualitative behavior of chiral phase transition in this model. Therefore, in this work we will continue to use this group of parameters. The dilaton profile Eq.2.3 has been shown to give well description of both chiral symmetry breaking and linear confinement. Thus in this work we will stick to this profile and extend the model to m l = m s case. The parameters µ 0 , µ 1 , µ 2 in the dilaton profile Eq.(2.3) will be taken as the same value (2.14) as in [5,63]. In the next section, we will show how to solve the order parameter of chiral phase transition from this model.

Boundary condition and numerical solutions
Under the background Eqs.(2.3,2.7,2.8,2.9), the equations of motion for χ l , χ s could be solved numerically. Before that, we should specify the boundary condition. Firstly, near the Ultraviolet (UV) boundary z = 0, one can extract the perturbative expansion solution of χ l , χ s as with c l , d l , c s , d s four integral constants of the two coupled second order ordinary derivative equations Eqs.(2.11,2.12), which could be related to the current quark masses m l , m s and chiral condensates σ l ≡ ūu = d d , σ s ≡ ss by the following equations [11,79] c l = m l ζ, (2.17) and ζ = √ 3 2π . As mentioned above, in this work we will try to study chiral phase transition under different values of m l , m s , so we will tune m l , m s in the later calculation. At a first sight, the other two integral constants σ l , σ s could be chosen independently on m l , m s . However, if one checks the two equations of motion Eqs.(2.11,2.12), there are terms like with f (z) in the denominator. At the horizon z = z h , we have f (z h ) = 0. Thus, z h is an apparent singularity of Eqs.(2.11,2.12), where χ l , χ s might be divergent. To avoid this divergence, physical solutions of χ l , χ s should also satisfy the infrared (IR) boundary conditions to cancel the singularity at the horizon where f (z h ) = 0. It is easy to understand that with these two additional conditions from the requirement of the regularity of χ l , χ s , the other two UV coefficients σ l , σ s cannot be considered as free integral constants with given m l , m s . Instead, they should be solved from the equations of motion. Therefore, to find a physical solution with given m l , m s , one has to solve the boundary value problem At low temperature, χ l and χ s are almost the same while at higher temperature they would be separated.
One can use the "Shooting Method" to solve this boundary value problem. After it was solved, one can extract the chiral condensates σ l and σ s . As an example, we take m l = 100MeV, m s = 0 , T = 50MeV, z h = 1 πT ≈ 6.37GeV −1 and T = 190MeV, z h = 1 πT ≈ 1.67GeV −1 . For T = 50MeV, z h = 1 πT ≃ 6.37GeV −1 , we get σ l = 0.10GeV 3 ≈ (470MeV) 3 , σ s = 0.12GeV 3 ≈ (495MeV) 3 using "Shooting Method" and plot the corresponding regular χ l , χ s solutions in Fig.2(a). The non-vanishing values of σ l , σ s at low temperature are signal of chiral symmetry breaking of the vacuum. From the figure, we could see that, at small z, χ l , χ s are slightly separated from each other, since the leading terms in this region are m l ζz and m s ζz, which are different when m l = m s . In the IR region when z is large, χ l , χ s is almost overlap and approach a constant value χ h l ≡ χ l (z h ) ≈ χ h s ≡ χ s (z h ) ≈ 0.47. Like in m l = m s case, the UV region of the solutions are governed by the trivial vacuum χ l = χ s = 0 while the IR region are governed by the non-trivial vacuum χ l = 0, χ s = 0.
Then at temperature up to T = 190MeV, we get σ l = 0.06GeV 3 ≈ (389MeV) 3 , σ s = 0.07GeV 3 ≈ (412MeV) 3 , which are smaller than the corresponding values at T = 50MeV, showing that chiral condensate are partly destroyed by temperature. In Fig.2(b) we plot the solutions of χ l , χ s at temperature T = 190MeV. From the figure, we could see that at high temperature, χ l , χ s are still interpolation of the trivial vacuum and non-zero horizon values χ h l , χ h s . The separation of χ l , χ s become larger than that at low temperature. As a short summary, we have shown that the regular condition of χ l , χ s would require the condensates σ l , σ s as functions of quark masses m l , m s and temperature T . Since chiral condensates are the order parameters of chiral phase transition, we will work out the quark mass and temperature dependence of σ l , σ s to extract the information of chiral phase transition in next section.

Chiral condensate and phase diagram in mass plane
As mentioned in the introduction, it is also interesting to consider the property of chiral phase transition when m l = m s . In Sec.2.2 we have shown that in the extended N f = 2 + 1 model one can solve chiral condensates σ l , σ s from the equations of motion Eqs. (2.11,2.12) with given quark masses m l , m s and the temperature T . Hence, in this section we will try to extract the quark masses and temperature dependence of chiral condensates, which contains the information of chiral phase transition. Firstly, we take m s = 0 and m l = 5MeV. Using "Shooting Method", we solve σ l , σ s from Eqs. (2.11,2.12). The results are shown in Fig.3(a). From the figure, we could see that since m l = 5MeV ≃ m s = 0, the differences of σ l and σ s are not very large. At low temperature, below 100MeV, both σ l and σ s are almost constants 0.1GeV 3 , showing the breaking of chiral symmetry in the vacuum. At high temperature, above 185MeV, σ l and σ s decrease to number smaller than 0.001GeV 3 . As we showed in [5], the non-zero value of the high temperature tail comes from the non-zero quark masses other than spontaneous chiral symmetry breaking, since in the chiral limit, this tail will tend to be zero. Therefore, in fact the high temperature tail stands for the symmetry restoration phase. So from the  numerical results in Fig.3(a), we see a low temperature symmetry breaking phase and a high temperature symmetry restoration phase, which indicates a phase transition. Then, we look into the intermediate temperature region and found that within 163MeV < T < 178MeV there are three branches of solutions at the same temperature. As discussed in [5], this kind of behavior is a characteristic signal of first order phase transition. The exact transition temperature would located inside the temperature region 163MeV < T < 178MeV and can be worked out from the free energy. However, here we will focus on the order of the transition other than the critical temperature, so we would not try to extract the exact transition temperature for the first order transition. Furthermore, we also plot the temperature dependence of the horizon value χ h l ≡ χ s (z h ), χ h s ≡ χ s (z h ) in Fig.4(a). There, we can see the same behavior as Fig.3(a). At small temperature, χ l , χ s are dominant by the non-trivial vacuum of scalar potential, while at high temperature they are dominant by the trivial χ l = 0, χ s = 0 vacuum.
Then, we increase m l to m l = 100MeV while keeping m s = 0. After solving the equations of motion, the results of chiral condensates σ l , σ s and χ h l , χ h s are given in Fig.3(b) and 4(b), respectively. There we could see that at low temperature, both σ l and σ s increase with the increasing of m l . Below T = 120MeV, σ l and σ s are almost constants σ l ≈ 0.11GeV 3 , σ s ≈ 0.12GeV 3 . It is also easy to see that σ s increase faster than σ l . As a result, the separation of σ l and σ s becomes larger than that at m l = 5MeV, m s = 0. At high temperature, again we could see that σ l , σ s decrease to a very small value, showing the restoration of the spontaneous breaking symmetry(though explicit breaking is always there due to the non-zero quark masses). However, different from m l = 5MeV case, in this case σ l and σ s decrease monotonically from the vacuum expectation values to zero without the triple branches region. Furthermore, at around T = 200MeV, σ l and σ s decrease very fast from the value at symmetry breaking phase to the value at symmetry restoration phase. This kind of behavior shows a characteristic crossover phase transition. As for χ h l , χ h s , we see that they also decrease monotonically from a larger value to zero. But differently, at low temperature, χ h l , χ h s do not change much with the increasing of m l . The low temperature value of these two quantities are still around 0.47. At high temperature, they monotonically decrease to zero at different rate. From the figure, we could see that χ h l decreases faster than χ h s .
( a ) ( b ) From the above discussion, it seems that when m s = 0, at small m l the system undergoes first order phase transition while at large m l the phase transition turns to be crossover. To be more rigorous, we fix m s = 0 and scan m l . From Fig.5, we find that when m l is smaller than 59MeV, both σ l , σ s are non-monotonic, indicating a first order phase transition. Then the non-monotonic region shrinks as the increasing of m l . At the critical value m l = 59MeV, the triple branches region disappears. At this value, we found that both the derivatives of σ l (T ), σ s (T ) with respective to temperature T diverge at the same temperature T = 189.3MeV, which reveals a second order phase transition. Then, above m l = 59MeV, we find that σ l and σ s decreases monotonically from nonzero value to zero, showing a crossover phase transition.
Then we increase m s to m s = 200MeV and scan m l . We plot the results of σ l , σ s in Fig.6. We see that qualitatively the results are similar to those when m s = 0. When temperature increases, condensates would decrease from nonzero value to zero. When Furthermore, we tune m s from 0 to 200MeV. We find that for each value of m s , there is a critical m l , at which chiral phase transition becomes second order. Here we note that both at large m s and small m s , the triple solutions region of σ l , σ s disappears at the same critical mass m c l and dσ dT diverges at the same critical temperature for both σ l and σ s . Below this critical m c l , the transition is of first order while above it the phase transition is crossover. We plot the critical line in Fig.7(a). From the figure, we see that the critical line divides the whole plane into two parts: the bottom left part is first order phase transition region while the upper right part is crossover transition region. Qualitatively, this result is in agreement with the "Colombia Plot" in Fig.1(a), which is summarized from lattice simulations and other effective methods. Moreover, from the critical line, one can extract the exact transition temperature, where dσ dT diverges for both σ l and σ s . The results are given in Panel.(b). There one can read that the transition temperature of the second order phase transition decreases when it approaches m l = m s ≈ 0.037GeV from left, while it increases when from right. This might indicate that the two branches of the critical line separated by the m l = m s point might be governed by different universality classes, though the exact correspondence is out of the scope of this work. Roughly speaking, this is also consistent with Fig.1(a), where the upper part of the critical line is governed by O(4) classes and the bottom part is governed by Z(2) classes.

Conclusion and discussion
To study QCD phase transitions in different situations are of great importance and it is interesting to consider chiral phase transition at different quark masses. In [5,63], we proposed a modified soft-wall AdS/QCD model and study chiral phase transition in N f = 2 and N f = 3 case. After extracting temperature dependent chiral condensate, it is found that chiral phase transition is of second order in two flavor chiral limit and turns to be crossover at any finite quark masses, while in three flavor chiral limit it becomes first order and turns to crossover only at sufficient large quark masses. Then, in [65], we have shown that this model can describe inverse magnetic catalysis after considering the Einstein-Maxwell sector. Therefore, in this work, we try to extend these studies to N f = 2 + 1 case when m u = m d = m s .
The extension of previous study is quite natural and simple. The main different is that the expectation value of the scalar field X in soft-wall model should be taken as 3 × 3 diagonal matrix diag{χ l , χ l , χ s } other than 2 × 2 diag{χ, χ}. If m l = m s , then it is expected that χ l = χ s and one has to deal with the two coupled second order derivative equations. The UV boundary condition of χ l , χ s can be related to quark masses m l , m s and chiral condensates σ l , σ s . The black hole horizon would come up with another boundary condition, which will require condensates as functions of temperature and quark masses, i.e. of the form σ l (m l , m s , T ), σ s (m l , m s , T ).
Fixing m s and solving the equations of motion, it is found that at both small m l and large m l , chiral condensate would decrease from finite value at low temperature to zero at high temperature, indicating a phase transition between symmetry breaking phase at low temperature and symmetry restoration phase at high temperature. Moreover, it is found that at small m l , σ l , σ s are triple valued in certain temperature range, giving the signal of first order phase transition. The triple valued temperature range would decrease with the increasing of m l , and at certain critical value it disappears and the phase transition become a second order one. If one continues to increase m l , then σ l , σ s will decrease monotonically and the phase transition becomes crossover. Varying m s , the qualitative behavior is similar. Thus, the whole m l − m s plane is divided into two regions: first order region and crossover region as shown in Fig.7(a). The boundary of these two regions is the second order line, which could be extracted by solving critical values of m l at different m s . The second order line are divided by the m l = m s line into two parts, and it is shown that the m s dependence of the transition temperature in these two parts are totally contrast, which might indicate that the two parts are governed by different universality classes. Qualitatively, these model results for chiral phase transition in Fig.7 is in agreement with the "Colombia Plot" in Fig.1(a) summarized from lattice simulations and other non-perturbative analysis. It confirms that the soft-wall AdS/QCD framework can provide good holographic description on chiral dynamics.