{'fraction_non_alphanumeric': 0.03975530906861379, 'fraction_numerical': 0.01730672760715679, 'mean_word_length': 4.58395580454404, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Progress in vacuum susceptibilities and their applications to the chiral phase transition of QCD

The QCD vacuum condensates and various vacuum susceptibilities are all important parameters which characterize the nonperturbative properties of the QCD vacuum. In the QCD sum rules external field formula, various QCD vacuum susceptibilities play important roles in determining the properties of hadrons. In this paper, we review the recent progress in studies of vacuum susceptibilities together with their applications to the chiral phase transition of QCD. The results of the tensor, the vector, the axial-vector, the scalar, and the pseudo-scalar vacuum susceptibilities are shown in detail in the framework of Dyson-Schwinger equations.


Introduction
Nowadays, the commonly accepted theory that describes the strong interaction is Quantum Chromodynamics (QCD). Due to the asymptotic freedom nature of QCD, in the high energy and high momentum transfer region reliable calculations can be done using the perturbative theory. Numerous comparisons between the results from theoretical calculations and experimental measurements have shown that QCD is correct in the high energy field. However, in the low energy and low momentum transfer region the coupling constant α s of QCD becomes large and running, consequently things have to be treated non-perturbatively; and thanks to the non-perturbative nature, the QCD vacuum is not trivial. People believe that the non-trivial properties of QCD vacuum are closely linked to spontaneous chiral symmetry breaking and confinement, which are the two main features of low-energy hadron physics. In order to characterize the QCD vacuum, people often introduce various vacuum condensates (such as the two-quark condensate, the gluon condensate, the mixed quark gluon condensates, and four-quark condensate, etc.) and many vacuum susceptibilities (such as the tensor vacuum susceptibility, the vector and the axial-vector vacuum susceptibilities, the scalar and the pseudo-scalar vacuum susceptibilities) phenomenologically [1,2,3,4,5,6]. In this paper, we review the recent progress in studies of vacuum susceptibilities together with their applications to the chiral phase transition of QCD. The results of various vacuum susceptibilities are shown in detail in the framework of Dyson-Schwinger equations (DSEs).
Metaphorically speaking, we can regard the QCD system as a black box, one way to study the QCD vacuum is to cause a perturbation by adding an external field to the system, such that we can learn something about the QCD vacuum indirectly by studying its responses to the external field. 1 In the theory of QCD, the quark propagator is the simplest Green function (two-point Green function), and its linear responses to the external fields reflect the properties of QCD vacuum directly. QCD vacuum susceptibilities are just the parameters that are associated with the linear responses of the quark propagator to the external fields, which can then characterize the non-perturbative properties of QCD vacuum. In the QCD sum rules external field formula, QCD vacuum susceptibilities play important roles in determining the hadron properties [2,3,5,6]. For example, the strong and parity-violating pion-nucleon coupling depends crucially upon the vacuum susceptibilities of pion [7,8], and the tensor vacuum susceptibility is closely related to the tensor charge of nucleon [9,10,11,12,13,14,15,16,17]. The nonlinear susceptibilities are also found to be correlated with the cumulant of baryon-number fluctuations in experiments [18,19,20,21,22]. Due to the importance of QCD vacuum susceptibilities in low energy hadron physics field, many different methods and models to study various QCD vacuum susceptibilities have been utilized over the years. Take the scalar vacuum susceptibility as an example, there are Lattice QCD calculations [23,24,25,26], calculations using the DSEs method [27,28,29], the multi-flavor Schwinger model [30,31], the linear sigma model [32,33], the (Polyakov loop extended) Nambu-Jona-Lasinio ((P)NJL) model [34,35,36], the nonlocal chiral quark model [37], and the chiral perturbation theory [38,39], etc. However, owing to the lack of a model independent expression for calculating the vacuum susceptibilities, different models or approximations often give different results, some of which may even be distinct. In this case, in order to provide a good starting point for studies of vacuum susceptibilities, the authors of Ref. [40] proposed a model independent equation for calculating the vector vacuum susceptibility firstly, using the QCD sum rules external field formula. It should be pointed out that, their results are based on the linear response theory to an external field of fermions, so that it is universal, and can also be used in other quantum field theories (e.g., Quantum Electrodynamics in (2+1) dimensions, QED 3 ) [41,42,43]. Based on this, people have also studied the axial-vector and the tensor vacuum susceptibilities in the framework of DSEs and under the rainbow-ladder approximation [44,45], while the authors of Ref. [46] studied the scalar vacuum susceptibility beyond the rainbow-ladder approximation, as well as the pseudo-scalar vacuum susceptibility [47]. Here we want to stress that the authors of Ref. [47] derived a model independent result for the pseudo-scalar vacuum susceptibility using the isovector-pseudo-scalar vacuum polarization, whereas the authors of Ref. [48] obtained model independent results for the vector and axial-vector vacuum susceptibilities using the vector and axial-vector Ward-Takahashi identities (WTIs). In Ref. [49], the authors calculated the tensor vacuum susceptibility employing the Ball-Chiu (BC) vertex [50,51], and also studied the dressing effect of the quark-gluon vertex on the tensor vacuum susceptibility. Their results show that the tensor vacuum susceptibility obtained in the BC vertex approximation is reduced by about 10% compared to its rainbow-ladder approximation value, which means that the dressing effect of the quark-gluon vertex is not very large in the calculation of the tensor vacuum susceptibility within the DSEs framework. These works promoted our understandings of the five kinds of vacuum susceptibilities, and up to now, only the scalar and the tenser vacuum susceptibilities are still model-dependent, to choose a more reliable model then becomes very important.
It is commonly accepted that with increasing temperature and/or quark chemical potential, strongly interacting matter will undergo a phase transition from the hadronic phase to the quark-gluon plasma (QGP) phase, which is also expected to appear in the ultra-relativistic heavy ion collisions [52], and the chiral symmetry which has been broken dynamically is also supposed to restore partially [53]. This is a very important as well as a very hot topic. One of the main goals for RHIC (Relativistic Heavy Ion Collider) and LHC (Large Hadron Collider) is just to create and then study such a new state of matter. These studies also play crucial roles in researches on the evolution of the early universe. The properties of the QCD vacuum will change when the temperature and/or quark chemical potential vary, as a result the vacuum condensates and vacuum susceptibilities which reflect the nature of the QCD vacuum should also change with temperature and/or quark chemical potential. Theoretically, the QCD phase transitions via the calculations of temperature and quark chemical potential dependence of the two-quark condensate or QCD vacuum susceptibilities can then be studied. In the standard definition, an order parameter is a quantity which is zero at one side of a phase transition and nonzero at the other side, so that the two-quark condensate in the chiral limit is such a quantity for the chiral phase transition (therefore, in some cases it is also referred to as "chiral condensate") [54]. However, while away from the chiral limit the two-quark condensate is no longer an exact order parameter (since there is no exact phase transition), although it can be stretched to encompass first order transitions where the quantity has jumps from one nonzero value to another. Various susceptibilities, which are related to the derivatives of the two-quark condensate, are usually used to locate the transition, and hence are also important quantities to characterize the properties of the system. 2 For example, the scalar susceptibility is often used to describe the QCD chiral symmetry restoration at finite temperature and finite chemical potential [23,28,35,55,56,57,58,59,60], while studies of the temperature and chemical potential dependence of the quark number susceptibility can provide useful information to the vicinity of the critical end point (CEP), where the first order phase transition meets with the crossover and a second order phase transition takes place [61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,46,79,80,81,82].
The cases with finite temperature and chemical potential represents a big branch in the field of studies on QCD vacuum susceptibilities. Hence it is difficult to give a comprehensive presentation on its theories and applications to the chiral phase transition of QCD. In this paper, we mainly discuss the calculations of various QCD vacuum susceptibilities at zero temperature and zero chemical potential, and then present some results in studies of the chiral phase transition of QCD. We will review the QCD sum rules external field formula in a nutshell, then take the constant vector field as an example to give the general derivation of the vector vacuum susceptibility, and obtain a model independent expression at last. We will also discuss the tensor, the vector, the axial-vector, the scalar, and the pseudo-scalar vacuum susceptibilities in detail together with their applications. From the model independent form of the vacuum susceptibilities, we see that their calculations are closely linked to the quark propagator (two-point Green function) and the corresponding vertex function (three-point Green function). Lattice QCD, which is commonly accepted as a first principle calculation of QCD itself, can give good results for the two-point Green functions, but for three-and more points Green functions there is still a long way to go. Accordingly, in order to do the calculations we have to, at present, turn to some effective models or approaches that are based on QCD itself. In this paper we mainly use the framework of DSEs approach, which is a very useful and successful non-perturbative method to treat the non-perturbative strong interactions as well as the properties of hadrons [83,84,85,86,87,88,89,90,91,92,93,94,95,96]. It must be emphasized here that, DSEs is an infinite series of coupled non-perturbative integral equations that the Green functions of continuous space-time should satisfy, of which a given n-point Green function is related with the m-point Green function, where m > n. In order to solve practical problems, people have to take "appropriate" approximations to the DSEs. 3 So far, the most commonly used approximation scheme is the rainbow-ladder approximation.
The rainbow approximation means that in the DSE for quark propagator the full quark-gluon vertex is replaced by the bare vertex, while the ladder approximation is that in the Bethe-Salpeter equation (BSE) of the quark-meson vertex the quark-anti-quark scattering kernel is replaced by its lowest order perturbative contribution, details of which can be found in Refs. [83,85,88]. The rainbow-ladder approximation is self consistent [97]. 4 Besides, rainbow-ladder approximation also promises the DSEs to satisfy the axis vector 2 Moreover, the two-quark condensate is actually unmeasurable, while various susceptibilities are reachable in the experiments. 3 A good approximation should holds the original symmetries of QCD itself as many as possible, and meanwhile, easy to calculate. 4 In the global color symmetry model (GCM) [86,98], if adopting the saddle point approximation, it can be found that the WTI, which ensures that pion is the Goldstone boson, and its mass is zero in the chiral limit. Through many years of research, it has been found that the DSEs under rainbow-ladder approximation have achieved great success in describing the pseudo-scalar and the vector mesons, since in these cases the contributions from higher order correction terms is little [85,87]. However, the results for the scalar meson are not satisfactory, since in this case the higher order corrections contribute significantly [84]. Therefore, for the scalar vertex and the scalar vacuum susceptibility the rainbow-ladder is not a good approximation. In other words, the results obtained in this approximation are not reliable, and one must go beyond it [46]. Accordingly, in the following discussions the tensor vacuum susceptibility is within the rainbow-ladder approximation (the results for the BC vertex can be found in Ref. [49], with a decrease of about 10%), while for the scalar vacuum susceptibility we try to go beyond it.
2. Derivation of the vector vacuum susceptibility using the QCD sum rules external field formula In this section we derive a model independent expression of the vector vacuum susceptibility using the QCD sum rules external field formula. To be specific, we will work in Euclidean space throughout this paper.
The vacuum susceptibilities are closely related to the linear responses of the full quark propagator to the external fields, so we add to the Lagrangian a term that the quark current coupled to an external field, here J Γ (y) is the quark current, V Γ (y) is the external field, q(y) is the quark field, Γ is a matrix in the direct product space of Dirac, flavor, and color. The linear response term of the full quark propagator to the external field can be written as [40], where (cc ′ ) and (αβ) are the color and spinor indices of the quark, |0 is the exact vacuum; G Γ,N P q (x) and G Γ,P T q (x) denote the non-perturbative and perturbative effects of the full quark propagator G cc ′ Γ αβ,q (x) under an external field, respectively. In the paper we adopt such a concept that, G cc ′ Γ αβ,q (x) is calculated in the Nambu vacuum, of which dynamical chiral symmetry breaking (DCSB) and quark confinement are two of the most important features; while G Γ,P T q (x) is calculated in the Wigner vacuum, where the system holds chiral symmetry (in the chiral limit case) and the dressed quarks are not confined (more discussions can be found in Sec. 3.1 of this paper, and Refs. [40,83,88]). As we will show in the following part, G Γ,N P q (x) is related to the vacuum susceptibilities and hard to calculate directly, while G cc ′ Γ αβ,q (x) and G Γ,P T q (x) can be calculated, and then we can get G Γ,N P q (x) as well as the vacuum susceptibilities indirectly.
In the QCD sum rules theory, the vacuum susceptibility is defined as [2,3,5,8], equation which quark propagator satisfy is just the unrenormalized DSE under rainbow approximation, while the one which quark-meson vertex satisfy is just the unrenormalized BSE under ladder approximation.
|0 represents the perturbative vacuum, χ Γ is the vacuum susceptibility, while H(x) is a function that characterize the non-local property of the non-local two-quark condensate, and H(0) = 1. From Eq. (3) we can see that, in order to calculate χ Γ we should know G cc ′ Γ,N P αβ,q (x) first, which means we should know G Γ,P T q (x) and G Γ q (x) simultaneously. Therefore, how to solve G Γ,P T q (x) and G Γ q (x) self-consistently becomes the key to get a model independent result. Next we will take a constant vector field V µ as an example to show how to get a model independent expression for the vector vacuum susceptibility [40].
When there is an external field, a term which is the coupling of the vector current and the vector external field should be added to the action of the system, namely, ∆S ≡ d 4 xq(x)iγ µ q(x)V µ (x). In this case, the dressed quark propagator in the chiral limit is where and Here we do not show the gauge fixing term, as well as the ghost field and its integration measure. Considering the linear response of the quark propagator to the external field, we can obtain T τ represents the time-ordering operation on Euclidean time τ , the full quark propagator without an external field (V µ = 0), while is the linear response term of the quark propagator to the external field. Hence, now the main task is to know d 4 y 0 |T τ [q(x)q(0)q(y)γ µ q(y)]|0 , we can try to do Taylor expansion of the inverse of dressed quark propagator to V µ [40,44,45], and then only keep the first order term of the external field, Here the vector vertex Γ µ is defined as 7 The derivation of Eq. (9) is general, so is applicable to both full and perturbative quark propagators. In configuration space it can be written explicitly as Then the linear response term of the full quark propagator G Z (x) to an external field is, and similarly for the perturbative quark propagator G Z,P T (x), where G P T (x) ≡ 0|T τ [q(x)q(0)] |0 , and Γ P T µ (p; 0) is the corresponding perturbative vertex. Now let's turn to the calculation of vacuum susceptibilities, the definition of the vector vacuum suscep- where G abV,N P αβ,q is the local two-quark condensate. It can be seen that the two-quark condensate is an explicitly gauge-invariant quantity, and is indeed computed in a gauge-invariant way for example in Lattice QCD simulations. That it depends on the renormalization scheme is on the other hand obvious. 5 If we multiply δ ab γ βα ν on both sides of Eq. (14), then From Eqs. (12) and (13) we can find that substituting them into Eq. (15), we can get a QCD sum rules based, also model independent expression of the vector vacuum susceptibility as, Here we note that, from the viewpoint of Feynman diagram, tr DC [γ ν GΓ µ G] and tr DC [γ ν G P T Γ P T µ G P T ] are respectively the values of the full and the perturbative vector vacuum polarizations at zero momentum, so the calculation of the vector vacuum susceptibility comes down to the calculation of the vector vacuum polarization. Here we want to stress that, as pointed out in Ref. [40], the derivation in this section is general and model independent for a constant external field, and so it can also be used to derive other vacuum susceptibilities similarly, such as the tensor vacuum susceptibility. However, for the axial-vector vacuum susceptibility, we need to take a variable external field, details can be found in Ref. [44]. The significance of the works in Refs. [40,44] is that they give theoretically model independent results of various vacuum susceptibilities, which provides a good starting point for further model calculations. In the next section we will discuss the tensor, the vector, the axial-vector, the scalar, and the pseudo-scalar vacuum susceptibilities respectively.

Calculations of various vacuum susceptibilities
After the model independent expressions of QCD vacuum susceptibilities are obtained, the next step is to use these expressions to calculate the values of these vacuum susceptibilities. From Eq. (18) we know that, in order to calculate the vector vacuum susceptibility, we must first know the quark propagator G(P ) (two-point Green function) and the vertex function Γ µ (P, 0)(three-point Green function). This holds for the other vacuum susceptibilities. However, as we mentioned earlier, nowadays it is difficult or even impossible to know these Green functions from the first principles of QCD directly, and people have to resort to various non-perturbative approaches besides the Lattice QCD, such as the QCD sum rules, the chiral perturbation theory, the DSEs theory, etc. Generally speaking, each approach has its own advantages and disadvantages. Since the early ninety's of the last century, the DSEs has made great progress, and also gained a lot of results which are in good agreement with the experimental measurements and Lattice QCD calculations [83,88,89,90,91,92,93,94,95]. This approach, therefore, is getting more and more attentions in many fields. In the following part of this section, we will discuss the calculations of the quark propagator and the quark-meson vertex in the framework of DSEs, and then discuss the calculations of several vacuum susceptibilities.

Calculations of quark propagator and quark-meson vertex in the framework of the DSEs
Using the functional path integral techniques, we can get the DSEs for Green functions from the Lagrangian (details can be found in the textbooks of Quantum Field Theory). Since the derivation of functional path integral does not rely on the perturbative theory, DSEs is strict (in the weak coupling limit, DSEs give the standard Feynman diagram expansion of the Green functions). The DSE of the renormalized quark propagator is here g 2 s D ab µν (p) is the renormalized gluon propagator, λ a is the Gell-Mann matrices with a the color index, Γ a ν (q, p) is the renormalized quark-gluon vertex, m bm is the bare mass of current quark and is related to Λ, Λ q ≡ Λ d 4 q/(2π) 4 represents some regularization scheme that keeps translation invariance, where Λ is the mass scale, we should take the limit Λ → ∞ at last. Z 1 (µ 2 , Λ 2 ) and Z 2 (µ 2 , Λ 2 ) are the renormalization constants of the quark-gluon vertex and quark wave function, which depend on both the renormalization point µ and the renormalization mass scale Λ. According to the Lorentz structure analysis, the solution of the quark DSE Eq. (19) has the following form Eq.(19) must be solved under certain renormalization condition. Since the asymptotic freedom nature of QCD, people usually demand that at a big space-like momentum square there should be where m(µ) is the renormalized current quark mass at µ. When chiral symmetry is explicitly broken, the renormalized mass and the bare mass satisfy this relation, is the renormalization constant of mass; while there is no explicit chiral symmetry breaking, Z 2 m bm = 0, which is just the case of chiral limit. Most of the following discussions on the vacuum susceptibilities are for the chiral limit case.
It can be seen clearly from Eq. (19) that, the quark propagator, the gluon propagator, and the quarkgluon vertex are related to each other. The gluon propagator and the quark-gluon vertex also satisfy their own DSEs, which are coupled to the Green functions with more points. Therefore, in order to obtain a closed equation of the quark propagator, we must make some appropriate approximations to the gluon propagator and the quark-gluon vertex. On the gluon propagator side, the researches of its DSE have made great progress, and the gluon propagator in Landau gauge is for all phenomenological purposes sufficiently well-known [83,88,95,96,101,99,102,100,107,108,109]. The lack of knowledge refers much more to the quark-gluon vertex. It is generally believed that, in studies of the color confinement, the DCSB, and the properties of bound states, the non-perturbative quark-gluon vertex plays an important role [110,111,112,113,114,115,116]. One of the most complete investigations of the quark-gluon vertex can be found in Ref. [100], and further guides as well as some important recent works and studies beyond the rainbowladder truncation can be found in Refs. [117,118,119] (and references therein). Therefore, when studying the dynamical chiral symmetry breaking problems of the quark propagator, people usually choose an appropriate gluon propagator model as input. A gluon propagator model often contains some parameters, which can be fixed by fitting the low-energy hadron physics experiments (such as the mass and decay constant of pion).
For the quark-gluon vertex, the case is more complicated. The most simple approximation is just the rainbow approximation, namely, use the bare quark-gluon vertex γ µ to replace the full one Γ µ (q, p). This is a big approximation, since the momentum dependence of the quark-gluon vertex is completely neglected, and of course, the non-perturbative information is also lost completely. In studies of low-energy hadron physics, such as the calculation of hadron masses, decay constants, and other physical quantities, it is necessary to solve the BSE that the quark-meson vertex satisfy. As mentioned above, the rainbow approximation of the quark DSE and the ladder approximation of the quark-meson vertex BSE are consistent with each other, and over many years of studies it has been demonstrated that although the rainbow-ladder approximation is the most simple one, it can give successful descriptions to some meson properties. The reason why it can not describe the scalar meson well has also become known. There have been many attempts that try to study the scalar meson beyond the rainbow-ladder approximation. However, this is rather difficult work, since if we do some approximations to the quark-gluon vertex of the quark DSE, 6 we also need to find a corresponding quark-meson BSE that is consistent with these approximations, which is really a difficult work. As we mentioned earlier, for studies of scalar vacuum susceptibility people should go beyond the rainbow-ladder approximation. The authors of Ref. [46] noted that in the calculations of the scalar vacuum susceptibility, only the scalar vertex with zero momentum exchange Γ(p, 0) is used, rather than the full scalar vertex Γ(p, q). Using this feature and the WTI, they give a method that can describe the scalar vacuum susceptibility beyond the rainbow-ladder approximation. In the following we will introduce the equations of the quark propagator and the quark-meson vertex under rainbow-ladder approximation first, then present a method to go beyond it.
Under the rainbow approximation and taking the chiral limit (m bm = 0), and substituting Eq. (20) into Eq. (19), it is easy to know the equations that A(p 2 , µ 2 ) and B(p 2 , µ 2 ) satisfy where we take the Landau gauge for the gluon propagator D ab µν (p) = δ ab δ µν − p µ p ν /p 2 D(p 2 ). After inputting the model gluon propagator, Eq. (22) can be solved by numerical iteration. From Eq. (22) we can see clearly that it has two distinct solutions: the Nambu-Goldstone solution (or Nambu solution) that B(p 2 , µ 2 ) = 0 and the Wigner-Weyl solution (or Wigner solution) that B(p 2 , µ 2 ) ≡ 0. In the Nambu-Goldstone phase: 1) chiral symmetry breaking occurs dynamically, the originally massless current quarks obtain masses by DCSB; 2) quarks are confined, since the quark propagator does not have the Lehmann representation. While in the Wigner phase, there is no DCSB, and the quarks are not confined [83,86,88].
In the Wigner phase, Eq. (22) becomes, here A ′ (p 2 , µ 2 ) represents the vector part of the self-energy function. Therefore, the chiral limit case of the quark propagator in the Wigner phase is In the ladder approximation, quark-meson vertex (three-point Green function) also satisfies its own BSE.
For the tensor vertex, its BSE is where the superscript m denote "meson", and then Γ means quark-gluon vertex while Γ m represents quarkmeson vertex. Similar to the quark propagator, we choose such a the renormalization condition According to Lorentz structure analysis, the general form of Γ m µν can be written as the following Now substituting Eq. (27) into Eq. (25) and after some algebra, we get five coupled integral equations of the The BSE that the vector vertex satisfy is Here one way to solve this is to write down the general Lorentz structure of Γ m µ (p, 0) and then do the numerical iteration, another more direct way is to use the Ward identity (which is obtained via the Ward-Takahashi identity by taking the infrared limit of the photon) by calculating the partial derivative of the quark propagator to the momentum p, we can know the vector vertex directly. Similar treatment can also be used to the axial-vector vertex, for the details please see Ref. [44].
For the scalar vertex, the rainbow-ladder approximation is not a good one and we have to go beyond it.
Let us see its BSE first, here K(q, p) is the quark-anti-quark scattering kernel, which is a four-point Green function that contains many non-perturbative information. As mentioned above, when beyond the rainbow-ladder approximation it is very difficult to find the corresponding BSE of the vertex. But if take a closer look we can find that only the scalar vertex (it is known from the functional path integral theory that taking the derivative of the two-point Green function with respect to the scalar external field will give the scalar vertex) Γ m (p, 0) = ∂G(p) −1 ∂m | m=0 is needed to calculate the scalar vacuum susceptibility. The authors of Ref. [46] skillfully used this point.
By the DSE of quark propagator, Eq. (19), it is easy to know the equation of scalar vertex where the summation over color is carried out to be 4/3 and Λ ν (q, p) = ∂Γν (q,p) ∂m .
Using the vector Ward-Takahashi identity, 7 we can know another identity that Λ ν (q, p) and Γ m (p) satisfy, This is an important step towards the consistent solutions of quark propagator and the scalar vertex.
By the Lorentz structure analysis, the inverse of the quark propagator G(p) −1 and the scalar vertex Γ has such general forms, And according to the works of Ball and Chiu [50,51], we take the following forms of Γ ν (q, p) and Λ ν (q, p) and F = A, B, C, D. Here it should be pointed out that, there are also many other forms of vertex ansatz that can also satisfy the vector WTI. The main reasons that the authors of Ref. [46] chose the BC vertex is that: firstly, there is no kinematic singularity; secondly, the momenta p and q are symmetrical. The main purpose of their work is to discuss the differences of the results between the bare vertex and a corrected one. 7 Actually, the QCD vertex does not satisfy this identity, but rather the more complicated Slavnov-Taylor identity (STI), which can be derived from the Becchi-Rouet-Stora invariance of QCD and corresponds to the WTI of QED. The difficult challenge is to find an Ansatz for the unknown renormalised propagators, vertices, etc., of QCD which satisfy the DSEs and respect the STIs of the theory. Progress can be made by satisfying a subset of the DSEs and STIs, and then supplementing these with information gleaned from Lattice QCD, phenomenology, etc. However, without ghosts Eq. (33) is identical to the corresponding WTI of QED [83]. Now substituting Γ ν (q, p) into the DSE of the quark propagator Eq. (19), we can find the numerical solutions of the vector function A and the scalar function B of the quark propagator. Then by substituting Γ ν (q, p), Λ ν (q, p), A, and B into Eq. (32) we can solve the scalar vertex. This method can be generalized to the calculations of vacuum susceptibilities at finite temperature and finite quark chemical potential. In the following part we will show the details of the gluon propagators and the numerical results.

Various vacuum susceptibilities
In this part we discuss the calculations and results of various vacuum susceptibilities in detail, namely, the tensor, the vector, the axial-vector, the scalar, and the pseudo-scalar vacuum susceptibilities.

The tensor vacuum susceptibility
Studies of tensor vacuum susceptibility are closely related to the tensor charge of nucleon [10]. for the non-local condensates sum rule and −0.0053 ± 0.0021 GeV 2 for standard sum rule [16]; while W.
Broniowski et al (1998) got -(0.0083 ∼ 0.0104) GeV 2 using the chiral constituent model [14], and H. T. Yang et al (2003) got -(0.0014 ∼ 0.0020) GeV 2 within global color symmetry model [17]. Clearly, in order to get a reliable theoretical prediction of the tensor charge, we should determine the tensor vacuum susceptibility as precise as possible. The work of Ref. [45] has provided a general treatment of the issue, and in Ref. [49] the calculations have been performed by employing the Ball-Chiu (BC) vertex [50,51]. The results of Ref. [49] show that the tensor vacuum susceptibility obtained in the BC vertex approximation is reduced by about 10% compared to its rainbow-ladder approximation value, which means that the dressing effect of the quark-gluon vertex is not very large in the calculation of the tensor vacuum susceptibility within the DSEs framework.
Through the previous general discussions of the vacuum susceptibilities (Sec. 2 of this paper) we know that, if we introduce a term that the quark tensor current coupled to a constant tensor external field q(x)σ µν q(x)Z µν (Z µν is the tensor external field), using the QCD sum rules approach we can obtain the general expression of the tensor vacuum susceptibility (here is the chiral limit case, but it is easy to generalize to the case with nonzero current quark mass) The expression in the earlier literatures is [14,15,16,17] Comparing these two equations it is easy to see that, the difference between χ Z and χ ′Z comes from the fact that in the previous calculations instead of using the exact dressed vertex Γ µν and the perturbative one Γ P T µν , they used the bare vertex σ µν , in the course of which the non-perturbative and perturbative dressing effects on the vertex were ignored. Here we want to stress that, both the full and the perturbative tensor vertexes have complicated momentum dependence. Especially for the full tensor vertex, there may exist mass singularities of various meson bound states, so compared with the bare vertex it contains quite a lot of sophisticated non-perturbative effects of QCD.
The discussions in Ref. [45] are under the rainbow-ladder approximation, and the renormalization procedure is necessary. The renormalized version of the tensor vacuum susceptibility is (for the details please see Ref. [45]) where G R , Γ R , G P T R and Γ P T R represent the renormalized versions of the full quark propagator, the full tensor vertex, the perturbative quark propagator and the perturbative tensor vertex, respectively. 0 | :q R (0)q R (0) : |0 is the renormalized two-quark condensate.
To show the numerical difference between χ Z and χ ′Z , we need to specify a model gluon propagator. The one Ref. [45] selected is the popular as well as famous Maris-Tandy model [87], which reads 234 GeV, and the renormalization point is chosen as µ = 19 GeV. As Ref. [87] pointed out, the parameters ω and m t are not independent, to chose m t = 0.5 GeV and ω = 0.3 GeV is to ensure that in the region 2 GeV 2 ≤ k 2 there is g 2 s D(k 2 ) ≃ 4πα s (k 2 ), where α s (k 2 ) is the "running" coupling constant of QCD. In Ref. [87] they chose D = 1.25 GeV 2 , and there studies show that this set of parameters can give very good results for the properties of pion and kaon.
Once the model gluon propagator is chosen, we can then calculate the quark propagator and the vertex through numerical iteration, as shown in Fig. 1 and Fig. 2 respectively. The final numerical results of the tensor vacuum susceptibilities are We see that, χ Z and χ ′Z are indeed different, but the difference is not very large. Let us now analyze the reasons for this. As stated above, the difference between them comes from the choices of vertexes. We know that in the perturbative region the bare vertex σ µν can be regarded as the zeroth approximation of the dressed vertex Γ µν , and therefore it is physically interesting to analyze the effects caused by their differences in the non-perturbative region. It can be seen clearly from Fig. 2 that, for large enough momentum square s = p 2 , Λ 1 (s, µ 2 ) tends to 1 while Λ 2 (s, µ 2 ) tends to 0, which means that in the large momentum region Γ µν tends to σ µν ; however, for small s, Λ 1 (s, µ 2 ) differs from 1 and Λ 2 (s, µ 2 ) differs from 0, although the differences are not significant. Therefore, the dressed vertex Γ µν (p, 0) obtained in the rainbow-ladder approximation of the DSEs approach is similar to that of the bare vertex in the small momentum region, and it is just this  small difference that accounts for the small numerical difference between the numerical values of χ Z and χ ′Z .
The small magnitude of the numerical difference between χ Z and χ ′Z suggests that the non-perturbative dressing effects on the dressed tensor vertex are not important, and therefore in this case to use σ µν as a replacement of Γ µν (p, 0) is a rather good approximation. This is quite different from the case of vector vacuum susceptibility [40].

The vector and axial-vector vacuum susceptibilities
For the vector vacuum susceptibility χ V and the axial-vector vacuum susceptibility χ A , in order to distinguish them from the later vertex wave function χ ν , here we use different symbols κ V and κ A instead.
According to the previous derivations, κ V can also be written as where Π V µν (P ) is the vector vacuum polarization while for κ A , with Π A µν (P ) the axial-vector vacuum polarization in which the trace is to be taken in color and Dirac space. V µ (z) and V 5µ (z) represent the variable vector and axial-vector external fields, χ ν , χ 5ν are the exact vector and axial-vector vertex wave functions, which can be written as GΓG using the no external field versions of the full quark propagator G and the corresponding vertex functions Γ. Z 2 is the renormalization constant of the quark wave function. According to the WTI, the vertex renormalization constant is also equal to Z 2 . qq is he chiral quark condensate.
As can be seen from Eqs. (43) and (45), in order to obtain the non-perturbative vector and axial-vector vacuum susceptibilities, we should subtract Π V,P T µν (P ) and Π A,P T µν (P ) respectively, which arises from the perturbative effects. In other words, when calculating the vacuum susceptibilities, the mean values of the perturbative vacuum should be subtracted. This is very important and is deserving of additional attention.
It is well known that the separation of the perturbative and the non-perturbative contributions from the mean values of the vacuum is somewhat arbitrary. Usually, this arbitrariness is avoided by introducing some normalization point [123]. In such a formula, the condensates will depend on the choice of the normalization point. There are also other methods besides this one, for example, in studies of the mixed quark-gluon condensate, the authors of Ref. [124] identified the perturbative vacuum with the Wigner vacuum, since both of them are trivial in the sense that there are no chiral symmetry breaking and confinement, in contrast to Nambu-Goldstone vacuum (the non-trivial vacuum) which corresponds to DCSB (more details can be found in Ref. [124]). In Refs. [40,44,45], the authors adopted the viewpoint of Ref. [124] to calculate the vector,  [126]. In the following, we will show how the authors of Ref. [48] obtain model independent results for the vector and axial-vector vacuum susceptibilities using the vector and axial-vector WTIs.
Let us discuss the vector vacuum susceptibility κ V first. The conservation of the vector current ensures that the vector vacuum polarization is purely transverse, namely, If we take the constant external field limit V µ (z) ≡ V µ , and use the following integration equation we can get It can be seen that the vector vacuum susceptibility is closely related to the vector vacuum polarization at zero total momentum. Now contracting both sides of Eq. (44) with δ µν , we then have The case for Π V,P T T (P 2 = 0) is similar.
In the chiral limit case, the wave function of the vector vertex satisfy the WTI where q ± = q ± P/2. Expanding the right side of Eq. (51) to P µ and taking the limit P µ → 0 leads to the Ward identity Substituting Eq. (52) into the right side of Eq. (50), and adopting the following parametrization of the quark propagator G(q) we obtain Since for large q 2 , σ V ∼ q −2 , this integral is quadratically divergent. However, this divergence is not genuine. Note that the integrand is a total divergence, so the above integral vanishes if a translation invariant regularization is adopted, 8 and hence For similar reasons we can get Π V,P T T (P 2 = 0) = 0 too. Thus we may draw the conclusion that the vector vacuum susceptibility is zero so long as the Ward identity is satisfied. In other words, the vanishing of the vector vacuum susceptibility is a direct consequence of gauge invariance. The authors of Ref. [40] also made use of the vector Ward identity in constructing the vertex function, but at that time they did not realize that the vector vacuum susceptibility can be expressed with such a simple formula, but their numerical result (of the order 10 −5 ) is very close to this model independent one.
The case for the axial-vector vacuum polarization is a little bit more complicated, since the axial-vector vertex contains a massless bound state pole. Nevertheless, the analysis is still straightforward, and the general form of the axial-vector vacuum polarization can be written as which contains both a transverse part and a longitudinal part. Following the same procedure of deducing the vector vacuum susceptibility, we can get From Eqs. (56) and (46), we know that Now we will determine the axial-vector vacuum polarization at P 2 = 0. In the chiral limit, the WTI for the axial-vector vertex can be expressed as where we also express the identity in terms of the dressed axial-vector vertex wave function χ 5µ instead of the dressed vertex Γ 5µ . Since χ 5µ possesses a longitudinal massless bound state pole, its generally expression is then [54] χ 5µ (q; P ) = γ 5 χ R µ (q; P ) +χ 5µ (q; P ) + f π P µ P 2 χ π (q; P ) .
In the limit P 2 → 0, the first term of Eq. (66) is (4 − ζ)4N c Z 2 Λ q (F R (q; 0) + 1 4 q 2 G R (q; 0)), the second term does not contribute, while the third term is (ζ − 1)f 2 π . Now using Eq. (65) we then get Note that the result is independent of ζ. can be expressed using the vector part of the quark propagator and the BS wave function of pion (details can be found in Ref. [54]) Substituting Eq. (68) into Eq. (67), and using Eq. (63), we get To obtain this result we have used the previous discussion that when using a translation invariant regularization, the integral to q in Eq. (69) is 0.
The axial-vector vacuum susceptibility is then From Eqs. (55) and (70), we found that in the chiral limit case which is just the result from Weinberg sum rules [127].
In summary, in this part we have derived model independent results of the vector and the axial-vector vacuum susceptibilities, which will play an important role in the related calculations of the QCD sum rules external field method. Using the vector and the axial-vector WTIs, Ref. [48] has demonstrated that in the chiral limit the vector vacuum susceptibility is 0, while the axial-vector vacuum susceptibility equals 3/4 of the square of the pion decay constant, which is the same as the result from Weinberg sum rules.

The scalar vacuum susceptibility
Studies on the QCD phase transitions play important roles in the studies of the evolution of the early universe and the high energy heavy ion collisions. The results of the Lattice QCD show that the QCD "phase transition" in the early universe is not a real phase transition, but a smooth crossover [58]. In Ref.
[58] they adopted the scalar vacuum susceptibility χ as an indicator, and studied the temperature dependence of χ.
The scalar vacuum susceptibility is very important in studies of QCD phase transitions, however, it is known that there is a square divergence in the χ(T ) calculation. In order to eliminate it, a usual treatment is to subtract χ(T = 0) from χ(T ), since the divergence only exists in χ(T = 0), and the temperature dependent part is not divergent [57,58]. Obviously, such method is not applicable in the zero temperature case, but the scalar vacuum susceptibility is meaningful at zero temperature. We will give a new treatment to the square divergence in χ(T ) later.
In QCD, the two-quark condensate associated with f flavor of quarks is defined as [128] qq where m f (µ) is the renormalized current quark mass, with µ the renormalization scale; Z 4 is the renormalization constant of the mass term in the Lagrangian, which depends implicitly on the gauge parameter.
In the chiral limit case the two-quark condensate can be used as an order parameter to describe the chiral phase transition, the scalar vacuum susceptibility (or chiral susceptibility) just characterize its response to the change of quark mass [31,129] herem is the auxiliary current quark mass that can be regarded as renormalization invariant. The order of integration and differentiation can be interchanged because the theory is properly regularized. Now use the Ward identity ∂ ∂m G(q; m; µ) = −G(q; m; µ)Γ(q, 0; µ)G(q; m; µ) then where Γ is the renormalized fully-dressed scalar vertex, which satisfies the inhomogeneous BSE Γ(k, P ; µ) = Z 4 I D + Λ q [S(q + )Γ 0 (q, P )S(q − )] sr K rs tu (q, k; P ). (78) with k the relative and P the total momenta of the quark-anti-quark pair; q ± = q ± P/2; r, s, t, u represent color and Dirac indices; and K is the fully-amputated quark-anti-quark scattering matrix.
The quantity m 2 Π 0 (P ) is renormalization point invariant in QCD, where Π 0 (P ) is the scalar vacuum To compare with Eq. (77), we can get the general result Hitherto we have not specified a regularization procedure for the scalar vacuum susceptibility. In this connection it is noteworthy that if a hard cutoff is used, then in the chiral limit of a non-interacting theory Ref. [128] pointed out that, this result can be traced to the dependence on current quark mass in Eq. (74).
On the other hand, the Pauli-Villars regularization would yield zero, which is then the procedure that the authors of Ref. [46] recommend and employ in models that preserve the one-loop renormalization-group behavior of QCD, and in this case the scalar vacuum susceptibility is defined self-consistent.
For the calculation of χ(µ) we only need the scalar vertex at P = 0, at which total momentum it has the general form Owing to the Ward identity, Eq. (76), it is not necessary to solve the inhomogeneous BSE since Hence a solution of the gap equation suffices completely to fix Γ 0 (k, 0).
A simplified form of the renormalization group improved effective interaction is employed in Ref. [46]; viz., they retain only that piece which expresses the long-range behavior [85,87] with t µν (q) the transverse momentum projection operator. Eq. (84) delivers an ultraviolet finite model gap equation. Hence, the regularization mass scale Λ can be moved to infinity and the renormalization constants equal to one.
The active parameters in Eq. (84) are D and ω, but they are not independent. In reconsidering a renormalization group improved rainbow-ladder fit to a selection of ground state observables [87], Ref. [93] noted that a change in D is compensated by an alteration of ω. This feature has further been elucidated and exploited in Refs. [130,131].
Herein, we employ ω = 0.5 GeV, and therefore D = m 3 g /ω = 1.0 GeV 2 corresponds to what might be called the real-world reference value for the bare vertex. Here we should note that, if we choose Eq. (84) as the model gluon propagator, but for the quark-gluon vertex using the BC ansatz Eq. (37), there is no reason to suppose that ω and D still satisfy the relation in Eq. (85), since the BC vertex contains considerable unknown non-perturbative information of interaction. Therefore, we need to refit the parameter D. After some tries, it has been found that if one assumes D = 0.5 GeV 2 , we can get the value of the two-quark condensate shown in the fourth column of Table 1. Hence we postulate that D = 0.5 GeV 2 is the real world reference value under the BC vertex.  The vector vacuum polarization is given by here Γ ν is the vector quark-anti-quark vertex that satisfies an inhomogeneous BSE similar to Eq. (78). This polarization couples to the photon and hence must be transverse while contain no mass term. Consider therefore [48] where the vector Ward identity was used in the second line and we have written the dressed quark propagator in the following form If we consider a free field and use a hard cutoff, then Eq. (87) becomes namely, the same sort of divergence is encountered as in the scalar vacuum polarization, Eq. (81). Naturally, any regularization scheme that preserves the vector Ward-Takahashi would yield zero as the result [132]. So, considering Eq. (87), we can define the regularized scalar vacuum susceptibility as Here, as long as each term in Eq. (90) is regularized independently in a valid fashion, then this is a "nugatory" transformation: since the photon is massless and hence in any valid scheme there must be Π V µµ (0; 0, Λ) = 0. In other words, with this definition we can move on to calculate the susceptibility without specifying a particular regularization scheme. The last two contributions in the integrand act as the Pauli-Villars terms, and Z 4 = 1 in the ultraviolet-finite model. for the BC vertex, as listed in Table 1.  It can be seen clearly from Fig. 3 that, the BC vertex ansatz has a quantitative impact on the magnitude and pointwise evolution of the gap equations solution. This can be anticipated from Ref. [133]. Moreover, the pattern of behavior can be understood from Ref. [112]: the feedback arising through the ∆ B term in the BC vertex always acts to alter the domain upon which A(p 2 ) and B(p 2 ) differ significantly in magnitude from their respective free-particle values, especially in the intermediate momentum region. The behaviors of  C(p 2 ) and D(p 2 ) shown in Fig. 4 can be understood easily from those of A(p 2 ) and B(p 2 ), so that there is no need to discuss it further.
The integrand in Eq. (90) is depicted in Fig. 5 for each vertex ansatz at the associated real world interaction strength. The resulting scalar vacuum susceptibilities are presented in the fifth column of Table 1.
It can be found from this figure that the regularization has served to eliminate the far-ultraviolet tail of the integrand, thereby ensuring convergence of the integral. We have varied the detailed form of the regularizing subtraction, namely, using free-field propagators and vertices instead of gap and BSE solutions, and the behavior of the integrand changes a little. Moreover, at real world values of the interaction strength for both ansatzs the integrands have negative support in the infrared and positive support for p 2 > (0.5GeV ) 2 .
These results tell us that when studying the scalar vacuum susceptibility, to go beyond the rainbow-ladder truncation is important and necessary, which then deserves further investigation.
In Fig. 6 we depict the evolution of the scalar vacuum susceptibility with increasing interaction strength, I = D/ω 2 . For I = 0 there is no interaction and then the "vacuum" is unperturbed by a small change in the current quark mass, hence the susceptibility keeps zero. With increasing I the susceptibility will grow since the interaction in theqq channel is attractive and therefore magnifies the associated pairing. This is equivalent to stating that the scalar vertex is enhanced above its free field value. Then the growth continues and accelerates until at some critical value I c the susceptibility becomes infinite, namely, a divergence appear.
The critical values are I c = 1.93 for the rainbow-ladder approximation, and I c = 1.41 for the BC vertex ansatz.
This divergence can be understood in this way: the models we have defined contain a dimensionless parameter I, which characterizes the interaction strength, and a current quark mass, which is an explicit source of chiral symmetry breaking. In the general theory of phase transitions, the latter is analogous to an external magnetic field while 1/I is kindred to a temperature. Consider the free energy for such theories  is a homogeneous function of its arguments in the neighborhood of t = 0 = m. From this it follows that the theory's magnetization exhibits the following behavior (e.g., see the Appendix of Ref. [27] for the details): and the associated magnetic susceptibility evolves in this way For a mean-field theory, β = 1/2, and δ = 3.
Now the nature of the critical interaction strength is easy to understand. In the class of theories we are considering, the quark condensate is analogous to the magnetization, and it is attended by the scalar vacuum susceptibility. For I < I c , the interaction is insufficient to generate a nonzero scalar term in the dressed-quark self-energy in the absence of a current quark mass; namely, DCSB is impossible and the model realizes chiral symmetry in the Wigner mode. But the situation changes at I c , and for I > I c a B = 0 solution is always possible. Moreover, the behaviors of the susceptibility show that each model undergoes a second-order phase transition and realizes chiral symmetry in the Nambu-Goldstone mode for interaction strengths above their respective values of I c . These observations emphasize the usefulness of the scalar vacuum susceptibility. With the bare vertex people can construct the pressure explicitly, 9 and thereby show 9 The pressure is defined as the negative of the effective action, and then the ground state of a system is that configuration for which the pressure is a global maximum. For the rainbow-ladder approximation, the corresponding pressure is just constructed from the famous CJT effective action [134]. that for I > I c the DCSB solution is dynamically favored because it corresponds to the configuration of the maximum pressure. Nevertheless, the diagrammatic content of the BC vertex ansatz is not knowable and hence an expression for the pressure cannot be derived. In this case we may rely on the behavior of some susceptibilities to conclude that DCSB is favored, as illustrated within many model studies. Once past the critical point, the susceptibility decreases as I increases. This is because the magnitude of the condensate order parameter grows in tandem with I and, therefore, the influence of any perturbation associated with a current quark mass must steadily diminish. We will do some related discussions in Sec. 4 of this paper. For the theoretical studies of the scalar and the following pseudo-scalar vacuum susceptibilities using the QCD sum rules, please refer to Ref. [15].

The pseudo-scalar vacuum susceptibility
The color-singlet current-current correlators (or equivalently, the associated vacuum polarizations) are directly related to physical observables, and hence play important roles in QCD. For example, the vector vacuum polarization is coupled to real and virtual photons, so that it is basic to the analysis and understanding of the process e + e − →hadrons [135,136]. The standard Lattice QCD methodology related to the correlators can be found in Ref. [137]; on the other hand, correlators are also amenable to analysis via the operator product expansion and are therefore fundamental in the application of QCD sum rules. In the latter connection, the pseudo-scalar vacuum susceptibility (or the pion susceptibility) plays a role in the QCD sum rules estimate of numerous meson-hadron couplings, for example, the strong and parity-violating pion-nucleon couplings, g πN N and f πN N , respectively [7,138,139]. Furthermore, the pseudo-scalar vacuum susceptibility is also a probe of QCD vacuum structure as the scalar susceptibility which is discussed in Sec. 3.2.3 of this paper [46], while its analysis is more subtle, with conflicts and misconceptions being common in different model calculations [7,8,33,139,140].
In this part we derive a model independent result for the pseudo-scalar vacuum susceptibility using the isovector-pseudo-scalar vacuum polarization [47], which can be written as where ζ represents the renormalization scale, the trace is over flavor and spinor indices, τ j is the Pauli matrices of the isospin, P is the total momentum of the quark-anti-quark pair; G(q ± ) represents the dressed-quark propagator, and Γ 5 is the fully dressed pseudo-scalar vertex, both of which depend on the renormalization point. Physical quantities obtained from Eq. (93) are then gauge invariant. The propagator can be obtained from the gap equation, Eq. (19).
The pseudo-scalar vertex is solved from an inhomogeneous BSE, namely here k is the relative momentum, r, s, t, and u represent color, flavor, and spinor indices, and with k ± = k±P/2, and K(q, k; P ) is the fully amputated two-particle irreducible quark-anti-quark scattering kernel. Now, let us consider the case with a space-time independent pseudo-scalar source, s 5 = 0, associated with the term in the action, we can define a pseudo-scalar vacuum condensate, the gauge-invariant as well as properly Then following the previous discussions, we can know that χ µν 5 (ζ) = −2Π µν 5 (P = 0; s 5 = 0,m, ζ).
Hitherto we have not specified a regularization procedure for the susceptibility, actually it can rigorously be defined via a Pauli-Villars procedure, as discussed in Ref. [46]. Now let us discuss the value of χ µν 5 (ζ) in the neighborhood of the chiral limit, therein we may write [54] iΓ µ 5 (k; 0) = where E R 5 (k; P ) is a part of the inhomogeneous pseudo-scalar vertex that is regular as P 2 + m 2 π → 0, Γ π µ(k; P ) is the pion bound state's canonically normalized Bethe-Salpeter amplitude, and r π (ζ) determined by iδ µν r π (ζ) = 0|q is the residue of this bound state in the inhomogeneous pseudo-scalar vertex. Using this, the pion's leptonic decay constant can be expressed as E R 5 (k; P ) can be determine via the axial-vector WTI with Γ µ 5ν the inhomogeneous axial-vector vertex. At P = 0 withm = 0 there is no pole contribution on the left side and hence in other words, this regular piece of the pseudo-scalar vertex is completely determined by the scalar part of them = 0 quark self-energy. Using the systematic, non-perturbative, and symmetry-preserving DSE truncation scheme introduced in Refs. [84,141], we can verify this equation order-by-order via the gap and Bethe-Salpeter equations. Now substituting Eq. (101) into Eq. (100), we then get so that, in the neighbourhood ofm = 0, the pseudo-scalar susceptibility splits into a sum of two terms: the first one expresses the contribution of the pion pole, O(m −1 ), and can be expressed in the closed form of which the last equality is proved in Ref. [54]; the second term O(m 0 ) can be determined via where qq (m; ζ, Λ) is the vacuum quark condensate defined in Eq. (74), and this entails where χ(ζ) is the scalar vacuum susceptibility defined in Eq. (75). Hence we arrive at a model independent consequence of chiral symmetry and the pattern by which it is broken in QCD, namely, : Pseudo-scalar vacuum susceptibility and related quantities computed using the two kernels of the BSE described in connection with Eqs. (37) and (84), also similar to those of Table 1. Dimensioned quantities are listed in GeV, and all are taken from Ref. [46] Vertex Now we show the numerical results of Ref. [47] in Table 2, which were computed from two models for the gap equations kernel.
Eq. (111) is a meaningful result. Recall that, in the absence of a current-quark mass, the two-flavor action has a SU L (2) SU R (2) symmetry and, moreover, that ascribing scalar-isoscalar quantum numbers to the QCD vacuum is a convention contingent upon the form of the current-quark mass term. It follows that, the massless action cannot distinguish between the continuum of sources specified by Therefore, the regular part of the vacuum susceptibility must be identical when measured as the response to any one of these sources, so that there is χ R = χ for all choices of θ. This is the content of the so-called "Mexican hat" potential, which is used in building effective models for QCD. The magnitude of χ depends on whether the chiral symmetry is dynamically broken and the strength of the interaction as measured with respect to the critical value required for DCSB, as discussed in Sec. 3.2.3 of this paper [46]. When the symmetry is dynamically broken, then the Goldstone modes appear, by convention, in the pseudo-scalarisovector channel, and thus the pole contributions appear in χ 5 but not in the chiral susceptibility. It is valid to draw an analogy with the Weinberg sum rules [48,127]. Eq. (111) also provides a novel, model independent perspective on a mismatch between the evaluation of the pseudo-scalar vacuum susceptibility using either a two-point or a three-point sum rules. The two-point study in Ref. [140] produces the pion pole contribution, χ π 5 , which is also the piece emphasized in Ref. [33], whereas a three-point method in Ref. [139] isolates the regular piece, χ R 5 , since a vacuum saturation ansatz is implemented in the derivation. Thus their analyses are not essentially in conflict, but emphasize different and independent pieces of the susceptibility, which can be distinguished. However, only the regular piece should be retained in a sum rules estimate of the pion-nucleon coupling constants [138]. We note in closing that other vacuum susceptibilities can be analyzed similarly.

Susceptibilities and the chiral phase transition of QCD
The chiral phase transition of QCD and a possible CEP on the QCD phase diagram have been drawing much attention on both experimental and theoretical sides; for example, beam energy scan (BES) program at the RHIC [142,143,144,145], and some theoretical studies [23,57,58,60,146,147,148,149,150,151,152,153,154,155,156,157,158,159]. To study the chiral phase transition of QCD, in the chiral limit the two-quark condensate can serve as an order parameter. Nevertheless, the current quark mass is small but not zero in the real world, so when studying the partial restoration of chiral symmetry at finite temperature and quark chemical potential, people also need to resort to some alternative order parameters.
When there is a coexistence of different phases, one of the candidates is the pressure difference (the bag constant) [160,161,162,163]. As discussed above, various susceptibilities are the parameters that related to the linear responses of the system to the external fields, which can then characterize the properties of the system, as well as to serve as indicators of corresponding phase transitions. In this section, we focus on the applications of some QCD susceptibilities to the chiral phase transition of QCD.
In quantum field theory, the dynamic properties of a system are fully characterized by the generating functional, which corresponds to the partition function in thermodynamics. When the system is in a certain phase, the generating functional is usually analytic for some choice of parameters, such as the current masses of the fermions, the chemical potential, and the temperature. The generating functional often exhibits a nonanalytic character while the phase transitions occur. So the location and characteristics of the chiral phase transition in the system can be determined by the behaviors of this quantity with respect to the corresponding parameters. In this case, a phase transition, in which the first order derivative of the generating functional with some of the parameters is discontinuous, is referred to as first order or discontinuous phase transition.
Second order or continuous phase transition exhibits continuity in first order derivative and a discontinuity or infinity in second order derivative. Nevertheless, thanks to the non-perturbative properties of QCD in the low energy region, such discussions are difficult and complicated. Accordingly, researchers often resort to some effective models, for example, see Refs. [43,163,164,165,166,167,168,169,170,171,172]. In the following, we take the NJL model, which is successful as well as popular in chiral phase transition studies [173,174,175,176,177,178,179,161,180,181,182], as an example, and give some discussions on the applications of some susceptibilities [36], similar discussions with in the framework of QED 3 can be found in Ref. [183]. For the details of the NJL model, please see Refs. [173,184,185,186].
The Lagrangian of the NJL model is (here we still work in the Euclidean space, and take the number of flavors N f = 2 while the number of colors N c = 3) where g is a coupling constant with the dimension of mass −2 , and the flavor and color indices are suppressed.
Now we introduce the definitions of four kinds of susceptibilities: the chiral susceptibility χ s , the quark number susceptibility χ q , the vector-scalar susceptibility χ vs , and another auxiliary susceptibility χ m . For mathematical convenience, we first introduce these susceptibilities in the free quark gas case, where the interaction term in the Lagrangian is zero, L int = 0. Denoting them with the superscript (0) , their definitions and expressions are then where Λ is a three-momentum cutoff which is introduced to avoid the ultraviolet divergence in the NJL the subscript f represents the free quark gas systems. It should be noted that χ vs have the same analytical expression, which is reasonable from the viewpoint of statistical mechanics, with Z f the QCD partition function in the free quark gas case.
In the interacting case, these susceptibilities are coupled with each other, Using the iterative method, we can obtain the numerical results of these susceptibilities. For example, the vector-scalar susceptibility is the response of the effective quark mass (M = m − 2g qq ) to the chemical potential µ, and its results are shown in Fig. 7. The parameters used in this section are taken from Ref. [185], namely, m = 5.5 MeV, g = 5.074 GeV −2 , and Λ = 631 MeV.
It can be seen from Fig. 7 that, when T is smaller than a critical value T c = 35 MeV, there always exists a convergent discontinuity of χ vs , corresponding to a first order phase transition; when T = T c , χ vs displays a sharp and narrow divergent peak, which implies a second-order phase transition, or in other words, here is a CEP on the phase diagram; when T > T c , the discontinuity disappears and a rather broad peak of finite height is shown, corresponding to the crossover region. Comparing this with the results of the chiral susceptibility χ s and the quark number susceptibility χ q , shown in Fig. 8, we can conclude that at low temperature, the first order phase transition occurs at almost the same chemical potential; while in the crossover region, if we pick the peak of these susceptibilities as the artificial critical points, we see that they tend to occur at different chemical potentials as the temperature increases.
The calculated results of χ vs , χ s , and χ q in the crossover region are shown in Figs. 9 and 10, respectively.
We can find that they exhibit different behaviors: the chiral susceptibility χ s exhibits an obvious band, so it is convincing to define the peak of χ s as the artificial critical points; in the high T and/or low µ region, the vector-scalar susceptibility χ vs tends to vanish; while the global shape of the quark number susceptibility χ q is just similar to the ones of χ s and χ vs , but it is non-vanishing in the high T and/or high µ region whose behavior is closely linked to the quark number density. Therefore, χ q can not describe the crossover property well in the high T and/or high µ region.
It is very interesting and meaningful to compare the results of these susceptibilities with that of the thermal susceptibility χ T = ∂ qq /∂T . For mathematical convenience, we also define χ n = ∂ q † q /∂T . By the same process employed above, we obtain a set of coupled equations for χ T and χ n as with E = M 2 + p 2 . The behavior of χ T in the crossover region is shown in Fig. 10, which is very similar to that of χ s .   11. The corresponding lines in the crossover region are determined by the peaks χ vs , χ T , and χ s . Here it should be stressed that, actually there is a very large number of studies investigating the phase diagram together with the CEP of QCD, and in various models [52,161,162,171,187,188,189,190,191,192,193] as well as functional and Lattice QCD studies [23,194,195,196,197,198,199,200] the results vary quite widely, most of them find a higher value for T c and a lower value for µ c of the CEP we plotted here, as summarised in Fig. 12 (here µ B = 3µ), which is taken from Fig. 4 of Ref. [52]. The points in Fig. 12 represent effective model and Lattice QCD predictions, the two dashed lines are parabolas with slopes corresponding to Lattice QCD calculations, and the red circles are locations of the freeze-out points for heavy ion collisions at corresponding center of mass energies per nucleon (indicated by labels in GeV). It is shown clearly that, different models or parameter sets give quite different predictions (please refer to Ref. [52] for more details of these points and lines). The parameters we adopted here (NJL89a in Fig. 12) are just for the convenience of showing the qualitative results of different susceptibilities more explicitly. For a recent parameter-independent attempt of such studies, please see Ref. [201].
We see from Fig. 11 that these susceptibilities split when the system tends to lower quark chemical potential and/or higher temperature. As we discussed above, up to now there is no exact order parameter for studies of QCD chiral phase transition when beyond the chiral limit, and different susceptibilities would characterize different physical properties of the system. Hence, related topics deserve further discussions.

Summary
In the Dyson-Schwinger equations (DSEs) framework of QCD, we summarize studies of various vacuum susceptibilities together with some of their applications to the chiral phase transition of QCD. In Section 1, we give a brief introduction of the vacuum susceptibilities together with their related applications. Then a general derivation of the vacuum susceptibility using the QCD sum rules external field formula is explicated in the following Section 2, as well as a model independent expression, which is expressed with the quark propagator (two-point Green function) and the corresponding vertex function (three-point Green function).
In the next Section 3, we review the calculations of the quark propagator and the vertex function under the framework of DSEs-BSEs first, and then use these theoretical methods and results to discuss the calculations of the tensor, the vector, the axis-vector, the scalar, and the pseudo-scalar vacuum susceptibilities further. For the tensor vacuum susceptibility, we introduce the calculation results of the Maris-Tandy gluon Here it also needs to be pointed out that, in the QCD sum rules external field approach, all kinds of vacuum condensates and vacuum susceptibilities are introduced as independent phenomenological parameters, and their values are determined by fitting the theoretical predictions from QCD sum rules to the experimental results. The authors of Ref. [202] proposed for the first time that the four-quark condensate may relate to the corresponding vacuum susceptibility. If this conclusion is correct, the QCD sum rules external field formula can then introduce less independent parameters, which will certainly increase its predictiveness.
Last but not least, it should be noted that up to now people still do not know how to rigorously define some of the vacuum condensates (including but not limited to the four-quark condensate) from the first principle of QCD. There is no doubt that, the issues related to the four-quark condensate and the relations between the four-quark condensate and the corresponding vacuum susceptibility are worth further studying. And of course, so are many discussions in this paper.  taken from Ref. [52] (please refer to this paper for more details).