Phase structure of NJL model with weak renormalization group

TL;DR

This paper employs the functional renormalization group to analyze the chiral phase structure of the NJL model at finite temperature and density, introducing a weak solution method to track RG flows beyond phase transition points.

Contribution

It introduces a weak solution approach to the RG equation, enabling analysis of RG flows after chiral symmetry breaking in the NJL model.

Findings

Successfully captures vacuum structures and critical phenomena.

Provides a detailed chiral phase diagram.

Demonstrates the effectiveness of weak solutions in RG analysis.

Abstract

We analyze the chiral phase structure of the Nambu--Jona-Lasinio model at finite temperature and density by using the functional renormalization group (FRG). The renormalization group (RG) equation for the fermionic effective potential $V(\sigma;t)$ is given as a partial differential equation, where $\sigma:=\bar \psi\psi$ and $t$ is a dimensionless RG scale. When the dynamical chiral symmetry breaking (D$\chi$SB) occurs at a certain scale $t_c$, $V(\sigma;t)$ has singularities originated from the phase transitions, and then one cannot follow RG flows after $t_c$. In this study, we introduce the weak solution method to the RG equation in order to follow the RG flows after the D$\chi$SB and to evaluate the dynamical mass and the chiral condensate in low energy scales. It is shown that the weak solution of the RG equation correctly captures vacuum structures and critical phenomena within the pure fermionic system. We show the chiral phase diagram on temperature, chemical potential and the four-Fermi coupling constant.

Figures

Tables

Table 1. Table 1: The (pseudo-)critical temperature at vanishing chemical potential, the CEP and the curvature of the chiral phase boundary ( 62 ) obtained from several models and methods in the two-flavor and three-color case. The abbreviation “QM” stands for the quark-meson model. “HQM” and “PQM” are the QM with higher order quark-mesonic scattering processes and the Polyakov-QM models, respectively. The numbers (I), (II) in the PQM denote different values of the parameter γ ^ ^ 𝛾 \hat{\gamma} governing the curvature of the critical Polyakov-loop temperature T 0 ​ ( μ ) subscript 𝑇 0 𝜇 T_{0}\!\left(\mu\right) ; see Schaefer et al. ( 2007 ); Herbst et al. ( 2011 ) .

Model	Method	pion mass	$T_c$($\mu =0$)	CEP $(T_c,{\mu }_c)$	${\kappa }_{\mu }$
		(or bare quark mass)	[MeV]	[MeV]
This work	FRG (LPA)	chiral limit	$182$	$(63,278)$	$1.900$
NJL (I) Buballa (2005)	MFA	chiral limit	$222$	$(112,286)$	$1.974$
NJL (II) Buballa (2005)	MFA	$135$ MeV	—	$(81,330)$	$1.974$
NJL (III)	MFA	$138$ MeV	$210$	(48, 328)	$1.947$
QM Schaefer and Wambach (2005)	FRG (LPA)	chiral limit	$143$	$(52,251)$	$1.135$
QM Braun et al. (2012)	FRG (LPA)	$138$ MeV	—	—	1.375
HQM Pawlowski and Rennecke (2014)	FRG	$138$ MeV	$166$	$(50,291)$	$1.397$
PQM (I) Schaefer et al. (2007)	MFA	$138$ MeV	$184$	$(150,168)$	—
PQM (II) Schaefer et al. (2007)	MFA	$138$ MeV	$184$	$(163,164)$	—
PQM (I) Herbst et al. (2011)	FRG	$138$ MeV	$190$	$(23,292)$	—
PQM (II) Herbst et al. (2011)	FRG	$138$ MeV	$190$	$(32,293)$	—
QCD de Forcrand and Philipsen (2002)	Lattice ($i\mu$)	$m_{q}a=0.025$	$174$	—	0.498
QCD Wu et al. (2007)	Lattice ($i\mu$)	$0.899\le \frac{m_{\pi }}{m_{\rho }}\le 0.943$	$171$	—	0.641
QCD Allton et al. (2003)	Lattice (Taylor)	$m_{q}a=0.1$	$170$	—	0.691

Equations

$$\displaystyle{\mathcal{L}}_{\rm NJL}=\bar{\psi}i{\ooalign{\hfil/\hfil\crcr$\partial$}}\psi-\frac{G_{0}}{2}\left[(\bar{\psi}\psi)^{2}+(\bar{\psi}i\gamma^{5}\psi)^{2}\right].$$

$$\displaystyle S_{\rm NJL}=\int\text{d}^{4}x\left[\bar{\psi}{\ooalign{\hfil/\hfil\crcr$\partial$}}\psi-\frac{G_{0}}{2}(\bar{\psi}\psi)^{2}\right],$$

$$\psi$$

$$Z=\int \mathcal{D}\varphi {\mathrm{e}}^{-S[\varphi ;{\Lambda }_0]}.$$

$$Z=\int \mathcal{D}{\varphi }_{<}\mathcal{D}{\varphi }_{>}{\mathrm{e}}^{-S[{\varphi }_{<}+{\varphi }_{>};{\Lambda }_0]}=\int \mathcal{D}{\varphi }_{<}{\mathrm{e}}^{-S_{\mathrm{eff}}[{\varphi }_{<};\Lambda ]},$$

$$\frac{\text{d}{S}_{\mathrm{eff}}}{\text{d}t}=-\frac{1}{2}{\int }_{\mathrm{shell}}\frac{{\text{d}}^{d}p}{(2\pi {)}^d}\left\{\frac{\delta S_{\mathrm{eff}}}{\delta {\varphi }_p}{\left(\frac{{\delta }^2{S}_{\mathrm{eff}}}{\delta {\varphi }_p\delta {\varphi }_{-p}}\right)}^{-1}\frac{\delta S_{\mathrm{eff}}}{\delta {\varphi }_{-p}}-\mathrm{str}\log \left(\frac{\overrightarrow{\delta }}{\delta {\varphi }_p}{S}_{\mathrm{eff}}\frac{\overleftarrow{\delta }}{\delta {\varphi }_{-p}}\right)\right\},$$

$$t:=\log \left({\Lambda }_0/\Lambda \right).$$

$$\displaystyle S_{\rm eff}[\psi,\bar{\psi};\Lambda]=\int\text{d}^{4}x\left[\bar{\psi}{\ooalign{\hfil/\hfil\crcr$\partial$}}\psi-V\!\left(\bar{\psi}\psi;\Lambda\right)\right].$$

$$V\left(\bar{\psi }\psi ;\Lambda \right)=\frac{G_{\Lambda }}{2}(\bar{\psi }\psi {)}^2+\frac{G_{4,\Lambda }}{4}(\bar{\psi }\psi {)}^4+\cdots .$$

$$V\left(\bar{\psi }\psi ;{\Lambda }_0\right)=\frac{G_0}{2}{\left(\bar{\psi }\psi \right)}^2,$$

$$\sigma \frac{\overleftarrow{\partial }}{\partial \psi }$$

$${\int }_{\mathrm{shell}}\frac{{\text{d}}^{4}p}{(2\pi {)}^4}=T\sum _{n=-\infty }^{\infty }\int \frac{{\text{d}}^{3}p}{(2\pi {)}^3}\Lambda \delta \left(|\vec{p}|-\Lambda \right).$$

$$\displaystyle\partial_{t}V\!\left(\sigma;t\right)=\frac{\Lambda^{3}}{\pi^{2}}\bigg{[}E+T\log\!\left(1+{\rm e}^{-\beta(E-\mu)}\right)+T\log\!\left(1+{\rm e}^{-\beta(E+\mu)}\right)\bigg{]},$$

$${\partial }_{t}V\left(\sigma ;t\right)=\frac{{\Lambda }^3}{{\pi }^2}\sqrt{{\Lambda }^2+{\left({\partial }_{\sigma }V\right)}^2}.$$

$${\partial }_{t}G=\frac{G^2{\Lambda }^2}{{\pi }^2}.$$

$${\partial }_t\tilde{G}=-2\tilde{G}+2{\tilde{G}}^2.$$

$$\tilde{G}\left(t\right)=\frac{{\tilde{G}}_{\mathrm{c}}{\tilde{G}}_0}{{\tilde{G}}_0-({\tilde{G}}_0-{\tilde{G}}_{\mathrm{c}}){\mathrm{e}}^{2t}},$$

$$t_{\mathrm{c}}=\frac{1}{2}\log \left(\frac{{\tilde{G}}_0}{{\tilde{G}}_0-{\tilde{G}}_{\mathrm{c}}}\right).$$

$${\partial }_{t}V\left(\sigma ;t\right)=-F\left(M\left(\sigma ;t\right),\sigma ;t\right),$$

$${\partial }_{t}M\left(\sigma ;t\right)+{\partial }_{\sigma }F\left(M\left(\sigma ;t\right),\sigma ;t\right)=0.$$

$${\int }_0^{\infty }\text{d}t{\int }_{-\infty }^{\infty }\text{d}\sigma \left({\partial }_{t}M+{\partial }_{\sigma }F\right)\phi \left(\sigma ;t\right)=0.$$

$${\int }_0^{\infty }\text{d}t{\int }_{-\infty }^{\infty }\text{d}\sigma \left(M\frac{\partial \phi }{\partial t}+F\frac{\partial \phi }{\partial \sigma }\right)+{\int }_{-\infty }^{\infty }\text{d}\sigma (M\phi ){|}_{t=0}=0.$$

$$\left(\begin{array}{ccc}{\partial }_{M}F& 1& 0\end{array}\right)\left(\begin{array}{c}{\partial }_{\sigma }M\\ {\partial }_{t}M\\ -1\end{array}\right)=0,$$

$$\left(\begin{array}{ccc}\text{d}\sigma & \text{d}t& \text{d}M\end{array}\right)\left(\begin{array}{c}{\partial }_{\sigma }M\\ {\partial }_{t}M\\ -1\end{array}\right)=0,$$

$$\frac{\text{d}\sigma \left(s\right)}{\text{d}s}=\frac{\partial F}{\partial M},$$

$$\frac{\text{d}t\left(s\right)}{\text{d}s}=1,$$

$$\frac{\text{d}M\left(\sigma ,s\right)}{\text{d}s}=0,$$

$$\frac{\text{d}V\left(s\right)}{\text{d}s}$$

$$=M\frac{\partial F}{\partial M}-F,$$

$$\frac{\text{d}\sigma \left(t\right)}{\text{d}t}=\frac{\partial F}{\partial M},$$

$$M\left(\sigma \left(t\right);t\right)=M\left({\sigma }_0;0\right)=:M_0,$$

$$\sigma \left(t\right)={\sigma }_0+{\int }_0^t\text{d}\tau {\partial }_{M}F\left(M_0;\tau \right),$$

Full text

††preprint: KANAZAWA-17-04

I Introduction

Solving Quantum Chromodynamics (QCD) is a challenging problem in elementary particle physics. Due to the strong dynamics of QCD in macroscopic scales, the phenomena such as the quark confinement and the dynamical chiral symmetry breaking (D$\chi$SB) are predicted. The exploration of the phase structure and the equation of state in hot and dense QCD are quite important for understanding the evolution of quark matters in the universe.

Numerous studies have been performed and have shown that QCD at finite temperature and density has various phase structures, e.g., the quark-gluon plasma and the color superconductivity. In particular, the study by lattice Monte Carlo simulation has substantially contributed to the understanding of nature of QCD. At present, lattice QCD, however, suffers from the sign problem. Since the Boltzmann weight at finite density has imaginary parts, the generated configurations oscillate and thus the statistical errors become large.

In contrast, the analytical methods, e.g., the mean-field approximation (MFA) and the Schwinger–Dyson equation (SDE) with the ladder approximation, have been applied to QCD and have succeeded to understand its dynamics and phase structures both qualitatively and quantitatively. Furthermore, effective theory approaches have also helped us to understand QCD. The Nambu–Jona-Lasinio (NJL) model Nambu and Jona-Lasinio (1961a, b) and its extended models Fukushima (2004, 2008), especially, have also played crucial roles as effective models describing the D$\chi$SB and the confinement in QCD; see also Hatsuda and Kunihiro (1994); Du et al. (2015); Rozynek and Wilk (2016). The Lagrangian of the NJL model is given by the following four-Fermi interaction structure:

[TABLE]

This Lagrangian is invariant under the chiral $\text{U}\!\left(1\right)$ transformation: $\psi\to{\rm e}^{i\gamma^{5}\theta}\psi$ and $\bar{\psi}\to\bar{\psi}{\rm e}^{i\gamma^{5}\theta}$. Such a four-Fermi interaction would be generated via the interaction between gluon and quark fields within QCD; see e.g., Mitter et al. (2015). For a larger coupling constant $G_{0}$ than a certain critical coupling constant $G_{c}$, the non-trivial vacuum, i.e., the chiral condensate $\langle\bar{\psi}\psi\rangle$, is generated.

Let us here consider the analysis of the D$\chi$SB in the NJL model at finite temperature and density with the functional renormalization group (FRG) method Aoki et al. (2000a, 1999); Kodama and Sumi (2000); Aoki et al. (2000b); Braun (2012); Braun and Janot (2011); Braun and Herbst (2012); Jakovac and Patkos (2013); Aoki et al. (2014a); Aoki and Yamada (2015); Braun et al. (2017); see also Aoki (2000); Bagnuls and Bervillier (2001); Berges et al. (2002); Polonyi (2003); Pawlowski (2007); Gies (2012); Rosten (2012) as the reviews of the FRG. The advantage of the FRG is that one can systematically improve approximations for the FRG analysis. However, as will be explained in the next section, the renormalization group (RG) equation of the NJL model with the strong interaction has singularities originated from the divergence of four-Fermi coupling constant. Therefore, it is difficult to follow the solution of the RG equation after arising the singularities. The method, which have been employed in order to overcome the singularities, is the auxiliary field method based on the Hubbard–Stratonovich transformation Aoki et al. (1999).111 Even if the auxiliary field is introduced at the initial scale, the four-Fermi interactions are generated in the middle of scale by the box diagrams made of the Yukawa interaction between the quark and the auxiliary field. To avoid it, the dynamical hadronization method (or rebosonization) has been developed and used to investigate the vacuum structure of QCD Gies and Wetterich (2002, 2004); Jaeckel and Wetterich (2003); Pawlowski (2007); Gies (2012); Braun (2009); Floerchinger and Wetterich (2009); Mitter et al. (2015); Braun et al. (2016); Rennecke (2015); Fu and Pawlowski (2015).

In this method, the four-Fermi vertex is replaced by exchanges of the bosonic auxiliary field $\phi$ via the Yukawa interactions $y\phi\bar{\psi}\psi$, namely, $G_{0}\sim y^{2}/m^{2}$.222 Note that it may seem that the number of free parameters increases from one to two (from $G_{0}$ to $y$ and $m^{2}$). However, one of them should be redundant. It has been shown in Aoki et al. (1999); Mitter et al. (2015); Braun et al. (2016) that the low energy physics actually does not depend on the initial value of the Yukawa coupling constant.

Using this method one can not only avoid the singularities but also include the mesonic quantum effects. However, analyses of the system becomes complicated due to the additional bosonic degrees of freedom and ambiguities of approximations within the FRG method. Although the method of bare mass is also used to system Aoki and Miyashita (2009), the operator expansion for the effective potential (15) does not converge in the small bare mass limit, and then it is difficult to obtain physical values in the chiral limit Aoki and Sato (2013). Moreover, the most serious problem in the method of bare mass is that one cannot address the chiral phase transition of first-order.

Recently, to overcome the singularities, a novel method called the weak solution method was suggested Aoki et al. (2014b). The solution of the RG equation with the singularities is mathematically defined as the “weak” solution. Hence, we call the RG equation applied the weak solution method “weak renormalization group equation”. It has been shown in Aoki et al. (2014b) that the weak solution method actually can address the chiral phase transitions of both second- and first-order and evaluate the physical values such as the dynamical mass and the chiral condensate.

In this paper, we analyze the D$\chi$SB in the NJL model at finite temperature and density using the FRG with the weak solution method. The chiral phase diagram of the NJL model will be investigated. It should be emphasized that we will obtain the chiral phase diagram of the NJL model using the FRG without the auxiliary field method. Previously, only the auxiliary field method within the FRG has been able to capture the first-order phase transition. The main purpose of this paper is to show that one can obtain the chiral phase diagram with both second- and first-order phase transitions by the weak solution method without the help of the auxiliary field method.

This paper is organized as follows: In the section II, we introduce the model and its effective model. The FRG method is briefly explained in order to apply it to the model. In the section III, we introduce the weak solution method and give the weak RG equation. The numerical results are exhibited in Sec. IV. We summarize and discuss our approach to the system and the results in Sec. V. In the appendix A, we show how to derive the RG equation. In the appendix B, the convexity and concavity of the beta function is discussed.

II Nambu–Jona-Lasinio model and weak renormalization group

Our purpose is to investigate the D$\chi$SB described by the NJL model in $d=4$ and the chiral phase transition at finite temperature and density. In this section, we give the action of the NJL model and the RG equation for the effective potential. We see that the effective potential in the NJL model becomes singular within the RG evolution.

II.1 NJL model

The Lagrangian of the NJL model with one flavor is given in Eq. (3). Since in this paper we would like to investigate the fundamental behavior of the D$\chi$SB by using the weak RG, we consider the following simplified NJL model with one flavor in four dimension Euclidean space:

[TABLE]

where we assume that the quark fields have the degree of freedom of color, $N_{\text{c}}$. Since the pseudo-scalar operator $(\bar{\psi}i\gamma^{5}\psi)^{2}$ is not introduced, this action has no continuous chiral symmetry. Instead, it is invariant under the discrete type transformation,

[TABLE]

which forbids the operators with negative power of $\bar{\psi}\psi$. Hence, the mass term of fermion $m\bar{\psi}\psi$ does not appear in the action. This model describes that the chiral symmetry breaks down due to the strong four-Fermi interaction and the non-trivial vacuum, i.e., the non-zero chiral condensate $\langle\bar{\psi}\psi\rangle\neq 0$, is realized. In this vacuum, the fermion obtains the dynamical mass which is proportional to the chiral condensate. It should be noted that we do not really work for the discrete chiral symmetric model with only scalar four-fermi interactions. Starting with the continuous chiral symmetric model with scalar and pseudo-scalar four-fermi interactions together, we define the renormalized fermionic interactions expressed in terms of scalar bilinear $\bar{\psi}\psi$ and pseudo-scalar bilinear $\bar{\psi}i\gamma_{5}\psi$ operators. The renormalized interactions have the same continuous chiral symmetry, since the RG beta functions are set to respect the symmetry.

After integrating the RG equation to the infrared end, we investigate the free energy (vacuum amplitude) as a function of scalar and pseudo-scalar source fields. The free energy is readily obtained making use of the translational invariance feature of our beta functions. By the Legendre transformation in a usual way, we define the effective potential as a function of $\sigma=\langle\bar{\psi}\psi\rangle$ and $\pi=\langle\bar{\psi}i\gamma_{5}\psi\rangle$, which has the same chiral symmetric form, i.e., it is a function of $\sigma^{2}+\pi^{2}$.

Therefore only a section of the effective potential where $\pi=0$ is sufficient to get the radial profile of it. It is shown straightforwardly that this section does not depend on the coupling constant of elementary pseudo-scalar interactions in the bare Lagrangian. Then we set the vanishing pseudo-scalar coupling constant. This is our strategy of model set up, which gives the shortest way of evaluating the effective potential profile.333 We comment on the finite number flavor case. In a general flavor case, the chiral symmetry is given as $\text{SU}\!\left(N_{\text{f}}\right)_{V}\times\text{SU}\!\left(N_{\text{f}}\right)_{A}\times\text{U}\!\left(1\right)_{V}\times\text{U}\!\left(1\right)_{A}$. $\text{U}\!\left(1\right)_{A}$ symmetry is broken by the chiral anomaly. This effective is described by the ’t Hooft–Kobayashi–Maskawa term Kobayashi and Maskawa (1970); ’t Hooft (1976, 1986) within the effective model approach. Since this term with $N_{\text{f}}=2$ is given as the four-Fermi interaction, we would expect no crucial difference from the present setup. In contrast, in the $N_{\text{f}}=3$ (or $2+1$) case, ’t Hooft–Kobayashi–Maskawa term is given by the six-Fermi interaction. Although the anomaly effect on the chiral phase transition is investigated in several literatures, e.g. Lenaghan et al. (2000); Roder et al. (2003); Mitter and Schaefer (2014), it might depend on chiral effective models.

We should, however, take account of the pseudo-scalar operator to calculate, for example, critical exponents, which are sensitive to the difference of the discrete and continuous symmetries.

II.2 FRG equation

We introduce the FRG to investigate the D$\chi$SB in the NJL model. Now, we briefly explain its idea and derivation. Let us consider the bare action $S[\phi;\Lambda_{0}]$ and the path integral,

[TABLE]

This integration can be interpreted as the summation of the quantum fluctuations with the momentum $0<|p|<\Lambda_{0}$. In the FRG, the path integral is evaluated by integrating out the quantum fluctuation with the higher momentum $\Lambda<|p|<\Lambda_{0}$, namely,

[TABLE]

where we divided the field $\phi\!\left(p\right)$ into the higher momentum mode $\phi_{>}\!\left(p\right)$ with $\Lambda<|p|<\Lambda_{0}$ and the lower momentum mode $\phi_{<}\!\left(p\right)$ with $0<|p|<\Lambda$. The action $S_{\rm eff}[\phi_{<};\Lambda]$ is called the Wilsonian effective action and is defined by the lower momentum mode. The FRG describes the change from the bare action $S[\phi;\Lambda_{0}]$ to the effective action $S_{\rm eff}[\phi_{<};0]$ through integrating the shell momentum $\Lambda-\delta\Lambda<|p|<\Lambda$, which is the so-called coarse graining. This step is formulated as a functional differential equation, and solving this equation is equivalent to evaluating the path integral (8). In this paper, we employ the Wegner-Houghton (WH) equation Wegner and Houghton (1973),

[TABLE]

where $\phi_{p}:=\phi\!\left(p\right)$ and we introduced the dimensionless scale parameter,

[TABLE]

It is rewritten as $\Lambda=\Lambda_{0}{\rm e}^{-t}$, and then the increasing $t$ corresponds to the decreasing the cutoff $\Lambda$. Since the two point function $\frac{\overrightarrow{\delta}}{\delta\phi_{p}}S_{\rm eff}\frac{\overleftarrow{\delta}}{\delta\phi_{-p}}$ is generally given as the supermatrix in field space, “str” denotes taking the trace for the supermatrix; see Clark et al. (1993); Oda and Yamada (2016) for the treatments of the supermatrix and its manipulations. For details of derivation of the WH equation, see e.g. Clark et al. (1993).444 The WH equation is also obtained from the sharp cutoff limit of the Plochinski equation Polchinski (1984); Sumi et al. (2000)

The Wilsonian effective action is generally spanned by an infinite number of effective operators generated by the interactions in the bare action. Note that symmetries which the bare action has are not broken through the shell momentum integral, hence, the generated effective operators respect their symmetries.

Although the WH equation (10) itself is exact, we cannot solve it without some approximations. That is, we have to restrict the theory space of the Wilsonian effective action to its subspace. We give the following action as the approximated Wilsonian effective action for the NJL model (6):

[TABLE]

If the effective potential is expanded into the polynomials of $\bar{\psi}\psi$, we have

[TABLE]

We call this potential the “fermionic” effective potential. As mentioned above, the odd powers of $\bar{\psi}\psi$ is not generated due to the (discrete) chiral symmetry. The potential at the initial scale $\Lambda_{0}$ is given by

[TABLE]

that is, $G_{\Lambda_{0}}\equiv G_{0}$ and other effective coupling constants vanish. Hereafter, we will use the notations $\sigma:=\bar{\psi}\psi$. Note that the variable $\sigma$ is not an auxiliary field, accordingly its derivatives with respect to $\psi$ and $\bar{\psi}$ are given as

[TABLE]

We will discuss the chiral phase transition by the effects of temperature and density in the next section. For this purpose, the density operator $\int\text{d}^{4}x\,\mu\bar{\psi}\gamma_{0}\psi$ should be introduced in the bare action (6) and the Wilsonian effective action (14).555 We assume that the chemical potential $\mu$ is not corrected by quantum fluctuations.

Moreover, the momentum in the time direction $p_{0}$ is replaced with the Matsubara frequency $\omega_{n}$ and its integration becomes the Matsubara summation. It is complicated to evaluate the shell momentum integral for four dimension momentum. To avoid it, we insert the cutoff $\Lambda$ into the momentum in the spatial direction. Then the shell momentum integral becomes

[TABLE]

The RG equation for the effective potential is given by

[TABLE]

where $\beta=1/T$, $E:=\sqrt{\Lambda^{2}+M^{2}}$ and $M:=\partial_{\sigma}V$. We show the explicit derivation of the RG equation in the appendix A. Note that the large-$N$ approximation ($N_{\text{c}}\to\infty$) and the uniform $\bar{\psi}\psi$ approximation were employed in order to obtain Eq. (19). The function $M$ corresponds to the dynamical mass of quark. We here call it the mass function.

Here, we revisit the NJL model at vanishing temperature and density. The RG equation (19) at zero temperature and zero density becomes

[TABLE]

Using the expansion of the effective potential given in Eq. (15), we obtain the RG equations of each coupling constant. In particular, the RG equation of the four-Fermi coupling constant is given by

[TABLE]

If we use the dimensionless rescaled coupling constant $\tilde{G}:=G\Lambda^{2}/2\pi^{2}$, it becomes

[TABLE]

Obviously, this equation has the fixed point $\tilde{G}_{\rm c}=1$ and yields the blow-up solution,

[TABLE]

for the initial value $\tilde{G}_{0}>\tilde{G}_{\rm c}$. We see the critical scale at which the denominator becomes zero,

[TABLE]

The divergence of $\tilde{G}$ itself is physically relevant since it corresponds to the chiral fluctuation $\tilde{G}\sim\langle(\bar{\psi}\psi)^{2}\rangle$. Thus, the divergence is the signal of the second order phase transition. However, the divergence disturbs solving the RG equation up to the infrared (IR) scale $\Lambda\to 0$. Since the four-Fermi coupling constant is lowest order of the effective potential (15), the effective potential becomes non-differentiable at its origin $V\!\left(\sigma=0\right)$ due to the singularity of the four-Fermi coupling constant. In other words, the derivative $\partial_{\sigma}V$ cannot be defined after the critical scale $t_{\rm c}$.

We see in the next section that the weak solution method Aoki et al. (2014b) allows us to analyze the D$\chi$SB without introducing the auxiliary field and to investigate the chiral first-order phase transition.

III Weak solution method

In this section, the weak solution of the RG equation (19) is defined. To construct it numerically, we introduce the method of characteristics which are given as the system of ordinary differential equations made from the partial differential equation. Moreover, to obtain a single-valued solution from a multi-valued solution, we use the Rankine–Hugoniot (RH) condition which are introduced by the definition of the weak solution. It is demonstrated how the solutions of characteristics with the RH condition choose a unique vacuum. Finally, we comment on the RG equation for the four-Fermi coupling constant from the viewpoint of the weak solution method.

III.1 Definition of weak solution

As discussed in the previous section, the RG equation for the effective potential becomes singular when the D$\chi$SB occurs. In Aoki et al. (2014b), it was shown that the weak solution method is valid for overcoming the singularity of this system. Here, we briefly review basic notions of the weak solution method and show the weak RG equation for Eq. (19).

The RG equation (19) is given as the partial differential equation,

[TABLE]

where $F\!\left(M,\sigma;t\right)$ is the beta function with a negative sign. Differentiating the both sides of RG equation with respect to $\sigma$, we obtain

[TABLE]

Let us now call this equation the “strong” RG equation in order to compare with the “weak” RG equation introduced below.

Here, we introduce a test function $\varphi\!\left(\sigma;t\right)$ which is smooth and converges to zero for both infinite $t$ and $\sigma$: $\varphi\!\left(\pm\infty;t\right)=\varphi\!\left(\sigma;\infty\right)=0$. Integrating Eq. (26) multiplied by this function, we have

[TABLE]

By utilizing the integration by parts, it becomes

[TABLE]

We call the RG equation (28) “weak RG equation”, and its solution is called “weak solution”. The derivatives with respect to $t$ and $\sigma$ move from the mass function and the beta function to the test function. Hence, $M\!\left(\sigma;t\right)$ given as a solution of the weak RG equation (28) can have an arbitrary number of singularities at any point on the $\sigma$–$t$ plane.

III.2 Method of characteristics

The RG equation of the mass function (26) is given as a partial differential equation. Here, we introduce the method of characteristics which reduces from the partial differential equation to the system of ordinary differential equations. The method of characteristics allows us to numerically solve the RG equation using e.g., the Runge-Kutta method and the finite-difference methods.

Let us consider a curved surface $M\!\left(\sigma;t\right)$ in the $\sigma$–$t$–$z$ space and rewrite the RG equation (26) as

[TABLE]

where we used $\partial_{\sigma}F=(\partial_{M}F)(\partial_{\sigma}M)$. The infinitesimal displacement vector $(\text{d}\sigma,\text{d}t,\text{d}M)$ is a tangent vector for a point $(\sigma,t,M)$ on the curved surface $M\!\left(\sigma;t\right)$, and satisfies

[TABLE]

since the total derivative of $M$ is given by $\text{d}M=\partial_{\sigma}M\,\text{d}\sigma+\partial_{t}M\,\text{d}t$. Thus, the vector $(\partial_{\sigma}M~{}\partial_{t}M\,-1)^{\rm T}$ is orthogonal to $(\text{d}\sigma~{}\text{d}t~{}\text{d}M)$, and we find that the vector $(\partial_{M}F~{}1~{}0)$ is also the tangent vector for a point on the curved surface $M\!\left(\sigma;t\right)$. Taking $(\text{d}\sigma~{}\text{d}t~{}\text{d}M)$ as the vector being proportional to $(\partial_{M}F~{}1~{}0)$, we have

[TABLE]

where we introduced an infinitesimal parameter $\text{d}s$ as a proportionality constant. The solutions $\left(\sigma\!\left(s\right),t\!\left(s\right),M\!\left(s\right)\right)$ of the system of ordinary differential equations with each initial conditions are called characteristics. Since the solution of Eq. (32) is simply $t=s$, we can replace from $s$ to $t$ in Eqs. (31) and (33).

The fermionic effective potential $V$ can be evaluated within the method of characteristics. Considering $V$ on the characteristics, we obtain

[TABLE]

where we used Eqs. (26), (31) and $M=\partial_{\sigma}V$.

To summarize, we have the following coupled equations for the present system:

[TABLE]

Actually, these equations are independent of each other, therefore, their solutions are given by

[TABLE]

where $\sigma_{0}:=\sigma\!\left(0\right)$. Note here that $F\!\left(M,\sigma;t\right)$ given in the right-hand side of Eq. (19) (with negative sign) does not explicitly depend on $\sigma$; $F\!\left(M,\sigma;t\right)\equiv F\!\left(M;t\right)$. We also note that this simplification happens in the NJL model with a certain approximation. Although for more general cases such as QCD, we have to solve the system of ordinary differential equations for a given system, solving it is easier than the case of the partial differential equation.

A schematic figure for the evolution of the mass function described by Eqs. (36)–(38) is shown in Fig. 1. Note that for a fixed $M_{0}>0$, $\sigma$ described by Eq. (37) moves towards the left-hand side since $\partial^{2}F/\partial M^{2}<0$ for Eq. (19) is satisfied; see the appendix B. The system of ordinary differential equations (36)–(38) does not yield any singularity. Instead, as a schematic figure shown in Fig. 2, after the critical scale $t_{c}$, the mass function becomes a multi-valued function around a value $\sigma^{*}$, at which the fermionic effective potential becomes the swallowtail form (the right-hand side of Fig. 2). However, the solution of the fermionic potential and the mass function after the D$\chi$SB have to be a single-valued function of $\sigma$ since there should exist only one physical vacuum. To uniquely determine the physical vacuum, we give the Rankine–Hugoniot condition in the next subsection.

III.3 Rankine–Hugoniot condition

To uniquely determine the weak solution from Eq. (28) (or equivalently Eq. (35)), we give a condition called the Rankine–Hugoniot (RH) condition. This condition is derived from the definition of the weak solution and reads

[TABLE]

where $\sigma^{*}\!\left(t\right)$ is a discontinuity point, and $\sigma_{+}^{*}$ and $\sigma_{-}^{*}$ denote closing to the point $\sigma^{*}$ from the right-hand side ($+$) and the left-hand side ($-$), respectively (see Fig. 2). This condition describes the evolution of $\sigma^{*}\!\left(t\right)$ and implies geometrically the differential-area invariance, namely,

[TABLE]

In the present case, since at the critical scale $t=t_{\text{c}}$, $S_{A}\!\left(t_{\text{c}}\right)=S_{B}\!\left(t_{\text{c}}\right)=0$ is satisfied, the constant value in Eq. (40) vanishes, and then $S_{A}\!\left(t\right)=S_{B}\!\left(t\right)$. In Fig. 2, the RH condition means that the vertical line (gray line) is defined such that the equal area rule $S_{A}\!\left(t\right)=S_{B}\!\left(t\right)$ is satisfied and the dotted line is ruled out as the solution. The weak solution being the solution of Eq. (28) with the RH condition is given by the solid line with a discontinuity at $\sigma^{*}$. The discontinuity determined by the RH condition is called “shock”.666 The terminology comes from the “shock wave” in gas dynamics.

As an example, the evolution of the mass function at vanishing temperature and density is shown in Fig. 3, where the dimensionless mass function $\tilde{M}=M/\Lambda_{0}$ and the dimensionless field $\tilde{\sigma}=\sigma/\Lambda_{0}$ are defined; see Eqs. (55) and (56) in the section IV. At the critical scale $t_{c}$, the slope of the mass function at $\tilde{\sigma}=0$ becomes infinity, which corresponds to the divergence of the four-fermi coupling constant. After $t_{c}$, the solution of the mass function becomes the multi-valued function. Applying the RH condition, the solution is uniquely determined; see the right-hand panel of Fig. 3. Next, let us consider the case at zero temperature and finite density. The evolution of the mass function in this case is shown in Fig. 4. We see that two singularities arise at the non-zero field values, $\tilde{\sigma}\neq 0$. They move towards the origin $\tilde{\sigma}=0$ and eventually merge with lowering the scale $t$. For such a behavior of the mass function, we observe the first-order phase transition as will be seen in the next section.

The mass function given as the weak solution satisfies

[TABLE]

since it reflects the chiral transformation $\sigma\to-\sigma$. Then, the dynamical mass is given by

[TABLE]

To summarize so far, the multi-valued mass function as a solution of the characteristics (36)–(38) is uniquely determined for any field value of $\sigma$ by the RH condition. Instead, the mass function with the RH condition has the discontinuities. Such a solution is defined mathematically as the weak solution of Eq. (28) in which the mass function and its beta function are not differentiated by $t$ and $\sigma$. The physical dynamical mass is given by Eq. (42).

Finally, we comment on the Legendre effective potential and its convexity and concavity. Let us start with defining the generating functional for the connected Green function,

[TABLE]

Here, the initial condition for the fermionic potential is given by

[TABLE]

where the last term is the source term. From the Legendre transformation of Eq. (43), we can define the Legendre effective potential,

[TABLE]

This potential is a function of the chiral condensate given by

[TABLE]

Note that the source $j$ is regarded as a function of $\phi$. It is known that the Legendre effective potential $V_{\rm L}$ becomes the convex form Fukuda (1976); O’Raifeartaigh et al. (1986); Tetradis and Wetterich (1992); Litim et al. (2006). The weak solution method automatically gives us the convex Legendre effective potential as the convex envelope of the nonconvex function obtained from the solution of characteristic equations (36)–(38). See discussions in Aoki et al. (2014b) on its proof.

III.4 RG flow of four-Fermi coupling constant within weak solution method

We discuss the RG flow of the four-Fermi coupling constant in view of the weak solution. Since the mass function $M\!\left(\sigma;t\right)$ is given by

[TABLE]

the behavior of the origin of $M\!\left(\sigma;t\right)$ is described by the four-Fermi coupling constant. At the critical scale $t_{\rm c}$ at which the four-Fermi coupling constant diverges, the slope of the mass function at $\sigma=0$ becomes infinity; see the orange line in Fig. 3. After the critical scale, it becomes negative whereas the RG flow of the four-Fermi coupling constant given by Eq. (22) cannot be evaluated due to the divergence. We have seen in the previous section that in the weak solution method with the RH condition, the mass function with the discontinuities is defined mathematically and uniquely determined for any field value of $\sigma$. Note that in the case where the first-order phase transition takes place, the slope of the mass function at the origin does not become infinity; see Fig. 4. That is, the RG flow of the four-Fermi coupling constant does not diverge at any scale.

Next, we discuss the relation between the inverse four-Fermi coupling constant $g:=1/G_{\Lambda}$ and the Legendre effective potential $V_{\rm L}$. We start with the well-known relation,

[TABLE]

where $W\!\left(j;t\right)$ is given in Eq. (43). The second derivative of $W$ with respect to $j$ corresponds to the susceptibility:

[TABLE]

Together with Eq. (48), we obtain

[TABLE]

When the second-order phase transition occurs, the susceptibility diverges: $\chi\to\infty$. That is, the left-hand side $\frac{\partial^{2}V_{\rm L}\!\left(\phi;t\right)}{\partial\phi^{2}}$ vanishes at the critical scale $t=t_{\rm c}$. More explicitly, we see

[TABLE]

At the origin of $V_{\rm L}$, the mass-squared $m^{2}$, i.e. the curvature of $V_{\rm L}$, corresponds to the inverse susceptibility $1/\chi$. The divergence of $\chi$ means the vanishing mass-squared $m^{2}=0$. As mentioned in the section II, the four-Fermi coupling constant corresponds to the chiral susceptibility:

[TABLE]

Then, the divergence of $G_{\Lambda}$ is the signal of the second-order phase transition, and we see that the inverse four-Fermi coupling constant corresponds to the curvature at the origin of the Legendre effective potential:

[TABLE]

The RG flow of the dimensionless rescaled inverse four-Fermi coupling constant $\tilde{\kappa}=\tilde{G}^{-1}=2\pi^{2}/(G_{\Lambda}\Lambda^{2})$ is given by

[TABLE]

where $\tilde{\kappa}_{0}$ is the initial value of $\tilde{\kappa}$ at the initial scale $t=0$. For $\tilde{\kappa}_{0}<1$, $\kappa\!\left(t\right)$ becomes negative after $t=t_{\rm c}$ at which $\tilde{\kappa}\!\left(t_{\rm c}\right)=0$. Thus, from the relation (53), one can see that the negative $\tilde{\kappa}$ corresponds to the negative curvature at the origin of the Legendre effective action. The negative curvature at the origin means that the non-vanishing vacuum $\langle\phi\rangle=\langle\bar{\psi}\psi\rangle\neq 0$ exists. However, one cannot evaluate the value of the chiral condensate by only observing the divergence of the four-Fermi coupling constant in the second-order phase transition. Besides, the first-order phase transition is not captured by only looking for the behavior of the four-Fermi coupling constant.

In the next section, we analyze the weak renormalization group at finite temperature and density and show the chiral phase structures of the NJL model.

IV Numerical analysis

In this section, we show the chiral phase transitions within the NJL model described by Eq. (19). To this end, we first give the RG equation written by the dimensionless variables. Then, the evolutions of the mass function with second- and first-order phase transitions are shown for a benchmark value of the four-Fermi coupling constant. The chiral phase diagram of the NJL model is shown. In the last subsection, we show our results at a benchmark point in comparison with results of other works.

IV.1 Dimensionless RG equations and initial conditions

We make the dimensionless beta function of $F\!\left(\sigma;t\right)$ given in Eq. (19) at finite temperature and density:

[TABLE]

where we defined the dimensionless variables,

[TABLE]

One obtains evolutions of the dimensionless mass function and the dimensionless effective potential by solving the dimensionless characteristic equations,

[TABLE]

with the initial conditions at $t=0$ ($\Lambda=\Lambda_{0}$),

[TABLE]

where $g_{0}:=G_{0}\Lambda_{0}^{2}/2\pi^{2}$. Note that the critical value of $g_{0}$ at vanishing temperature and density is $g_{0}{}_{\text{c}}=1$. Hence, for $g_{0}>1$, the D$\chi$SB is observed.

IV.2 Phase transitions

Here setting $g_{0}=2$ as a benchmark point, we show the dependence of the mass function on temperature and density. First, let us consider the zero density case $\mu=0$. We show the dependence of the mass function in the IR limit $t\to\infty$ ($\Lambda\to 0$) on temperature in Fig. 5. With increasing $\tilde{T}$, the value of the mass function at $\tilde{\sigma}\to 0$ smoothly becomes smaller. This behavior corresponds to the second-order phase transition. Note that the dynamical mass is given in Eq. (42). The mass function at $\tilde{\sigma}\to 0$ vanishes at ${\tilde{T}}_{\text{c}}\simeq 0.589$ which is the critical temperature. Next, we consider the finite chemical potential and vanishing temperature case. In this case, the evolution of the mass function in the IR limit is shown in Fig. 6. We see that the slope of the origin of the mass function changes from negative to positive at $\tilde{\mu}\simeq 0.707$; see the left-hand side of Fig. 6. This means that the four-Fermi coupling constant becomes positive via the divergence, and thus the Legendre effective potential has a local minimum at its origin. Here we call this behavior the “second-order” although it does not always correspond to the phase transition. Note that there is still a non-trivial global minimum at $\tilde{\mu}\simeq 0.707$. Above $\tilde{\mu}_{\text{c}}\simeq 0.861$, the discontinuities of the mass function in the IR limit are located at non-vanishing $\sigma$. Then, when taking the limit $\sigma\to 0$, the physical dynamical mass suddenly vanishes with increasing $\tilde{\mu}$, i.e., $M\!\left(\sigma\to 0;t\to\infty;T=0,\mu<\mu_{\text{c}}\right)\neq 0\to M\!\left(\sigma\to 0;t\to\infty;T=0,\mu>\mu_{\text{c}}\right)=0$. This phase transition is of first-order. We see that the weak solution method can capture the phase transitions with both second- and first-order. In the next subsection, we show the phase diagram on the $g_{0}$–$\tilde{\mu}$–$\tilde{T}$ plane.

IV.3 Phase diagram

In Fig. 7, the chiral phase diagram on the $g_{0}$–$\tilde{\mu}$–$\tilde{T}$ plane are shown. The red line stands for the second-order phase boundary. As discussed in the section II, this boundary can be obtained by only the behavior of the four-Fermi coupling constant. The blue line corresponds to the first-order phase boundary. The biggest and smallest phase boundaries correspond to the cases of $g_{0}=2$ and $g_{0}=1.0101$, respectively. We show the phase diagram on the $\tilde{\mu}$–$\tilde{T}$ plane in Fig. 8, which is obtained by the projection from the three dimensional phase diagram of Fig. 7. A schematic figure of a phase boundary with a fixed four-Fermi coupling constant is shown in Fig. 9. Note that the region, which is surrounded by the second- and first-order phase boundaries towards the high density side, is the broken phase. As discussed in the previous subsection, although the four-Fermi coupling constant diverges at the second-order phase boundary, that is, the curvature of the origin of the effective potential becomes positive, the global minimum is still located at the non-vanishing chiral condensate. In Fig. 10, the behaviors of the critical end-point (CEP) and the second-order phase boundary at $T=0$ with varying the four-Fermi coupling constant are shown. At a finite value of the four-Fermi coupling constant $g_{0}\simeq 1.2$, they merge and the CEP vanishes.777 We see in Fig. 8 that the first-order phase boundary abruptly disappears at a certain four-Fermi coupling constant. However, the CEP should completely merge with the second-order phase boundary at $T=0$. In order to see this, we need very precise numerical calculations. Here, extrapolating the CEP towards the horizontal axis, we obtain the value of the four-Fermi coupling constant at which the critical end-point merges with the second-order phase boundary at $T=0$.

The first-order phase boundary vanishes for the smaller four-Fermi coupling constant than $g_{0}\simeq 1.2$. We again emphasize that the chiral phase diagrams shown in this work are obtained by the RG equation without introducing the auxiliary field.

Here, we note the case where the finite bare mass is taken into account. In this case, we could observe the crossover for smaller chemical potential rather than the second-order phase transition. This behavior within the weak renormalization group is shown in Aoki et al. (2014b).

IV.4 Comparisons with other models and methods

Let us compare our results with other models and methods. We first note that as shown in Refs. Aoki et al. (2000a); Terao (2001), the FRG analysis in the large-$N$ limit is equivalent to the SDE with the ladder approximation. The ladder approximated SDE is also equivalent to the MFA. Therefore, these methods yield quantitatively the same results.

We compare the present work with the other ones in the two-flavor and three-color case ($N_{\text{f}}=2$, $N_{\text{c}}=3$).888 In the present setup, there are no significant differences between one-flavor and two-flavor cases. See the comments in the footnote 3.

For the mean-field analysis of the NJL model, we employ the following Lagrangian,

[TABLE]

Note that the critical four-Fermi coupling is $g_{0}{}_{\text{c}}=1/N_{\text{f}}N_{\text{c}}=1/6$. In Ref. Buballa (2005), the UV cutoff, the four-Fermi coupling constant and the bare quark mass are set to $\Lambda_{0}=587.9$, $g_{0}=0.247$ and $m_{q}=5.6$ MeV, respectively. These parameters yield $M=400$ MeV, $f_{\pi}=92.4$ MeV and $m_{\pi}=135$ MeV at vanishing temperature and density. We have reexamined the chiral phase structure with the following initial values: $\Lambda_{0}=631$ MeV, $g_{0}=0.2217$ and $m_{q}=5.5$ MeV Hatsuda and Kunihiro (1994). In this case, we obtain $M=336$ MeV, $f_{\pi}=93$ MeV and $m_{\pi}=138$ MeV at zero temperature and density. For the present FRG analysis, we used the following initial values as a benchmark point: $\Lambda_{0}=610\,\text{MeV}$, $g_{0}=1.3158$. These values yield the dynamical mass $M=301$ MeV at $T=\mu=0$. In the other works with the analytic methods Schaefer and Wambach (2005); Braun et al. (2012); Pawlowski and Rennecke (2014); Schaefer et al. (2007); Herbst et al. (2011), the free parameters are set such that one obtains the pion decay constant $f_{\pi}=93$ MeV and the pion mass $m_{\pi}=138$ MeV in the IR regime.999 The work in Schaefer and Wambach (2005) obtains $f_{\pi}=87$ MeV and $M\sim 300$ MeV in the chiral limit. When the effects of the finite and realistic pion mass are taken into account, the decay constant becomes $f_{\pi}=93$ MeV.

Refs. de Forcrand and Philipsen (2002); Wu et al. (2007) performed the lattice QCD simulation with imaginary chemical potential $\mu_{I}=-i\mu$.101010 Ref. de Forcrand and Philipsen (2002) uses two-flavor staggered quarks, whereas two-flavor Wilson quarks is used in Ref. Wu et al. (2007).

By using analytic continuation to real (quark) chemical potential $\mu$, the dependence of the critical temperature on small chemical potential is investigated. In Ref. Allton et al. (2003), the Taylor-series expansion method is employed.

In Table. 1, the (pseudo-)critical temperature $T_{\text{c}}\!\left(\mu=0\right)$ and the CEP obtained from various models and methods are collected. We also show the value of the curvature of the chiral phase boundary $\kappa_{\mu}$ which is defined as

[TABLE]

Note that $\kappa_{\mu}$ does not depend on the cutoff scale $\Lambda_{0}$. We evaluate the curvature in a range $0\leq\mu/(\pi T_{\text{c}}\!\left(0\right))\lesssim 0.1$ for the present analysis and the mean-field analysis of the NJL model (III).

We see from Table. 1 that there are no drastic differences in the critical temperature at vanishing chemical potential and the CEP. In contrast, the curvatures obtained from our work and the mean-field analysis of the NJL model are somewhat larger than that of the others. This could be because the mesonic fluctuations are not taken into account in the present FRG and the MFA computations of the NJL model. Fig. 11 exhibits the dependence of the curvature $\kappa_{\mu}$ on the dimensionless four-Fermi coupling $g_{0}$. We see that $\kappa_{\mu}$ tends to be larger with increasing $g_{0}$. Note that the value of the critical four-Fermi coupling constant is $g_{0}{}_{\text{c}}=1$, and as we have seen in the previous subsection, the CEP vanishes below $g_{0}\simeq 1.2$.

V Summary and discussion

In this paper, we have studied the chiral phase structure of the NJL model at finite temperature and density using the FRG. We have seen that when the D$\chi$SB takes place, the solution of the RG equation has to become singular and cannot be evaluated after a certain critical scale. To overcome this situation, we have introduced the weak solution method to the RG equation. It has been shown that the weak solution method can appropriately choose the vacuum and describe the chiral phase transition with both second- and first-order. Previously, in the FRG method, the physical values such as the chiral condensate and the dynamical mass have been evaluated by only introducing the auxiliary field method. In particular, it was difficult to obtain the chiral phase diagram with the first-order phase transition. We have shown that the weak solution method allows us to evaluate the chiral dynamics within the pure fermionic system.

Using the FRG with the weak solution method, the D$\chi$SB in QCD at zero temperature and zero density was studied in Aoki and Sato (2013), where effects beyond the ladder approximation (non-ladder effects) on the dynamical mass and the chiral condensate are involved.111111The crossed ladder diagrams and the anomalous dimension of the quark field are included. These effects cannot be taken into account by the SDE with the ladder approximation and the MFA.

It was shown that the gauge dependence of the physical values is suppressed thanks to including the non-ladder effects. As a next step, we should investigate the D$\chi$SB in QCD at finite temperature and density. The weak solution method also can be applied to systems with the color superconductivity and ones in an external magnetic field.121212 See e.g., Fukushima and Pawlowski (2012); Mueller and Pawlowski (2015); Aoki et al. (2016) for the analyses of the system in an external magnetic field using the FRG.

In future works, we will apply the FRG with the weak solution method to such systems.

Here, we discuss several issues of the weak solution method. First, we comment on the cutoff scheme and the convexity and concavity of the beta function. Through this paper, we have analyzed the RG equation obtained by using the WH equation which is formulated with the sharp cutoff regularization. When the smooth cutoff scheme, e.g., the optimized cutoff regularization Litim (2001), is employed, we obtain the different form of the beta function for the effective potential; see Eq. (83). As shown in the appendix B, for the sharp cutoff regularization, the beta function for the fermionic effective potential (19) is concave, i.e., $\partial^{2}F/\partial M^{2}<0$ for an arbitrary RG scale, temperature and density. Therefore, the motion of $\sigma\!\left(t\right)$ described by Eq. (37) behaves as shown in Fig. 1. In contrast, for the case where the optimized cutoff regularization is used, the concavity of the beta function is not guaranteed. Hence, there are solutions in which the motion of $\sigma\!\left(t\right)$ comes back to the right-hand side as shown in Fig. 12.131313Such a behavior in terms of the four-Fermi coupling constant was reported in Aoki et al. (2014a); Aoki and Yamada (2015). In this case, however, we need some additional information about the solution of the mass function removed once by the RH condition, which corresponds to the dotted line in Fig. 12, in order to define the mass function as a function of $\sigma$. At present, we do not know how to mathematically define the weak solution for such a “come-back” solution.141414 It should be stressed here that the “come-back” solution (physically correct equal area law solution) is a weak solution because it satisfies the RH condition, but it is not a unique solution since there always appear “rarefaction” solutions in this type of situation. We do not even know if we can define some additional condition to assure uniqueness of weak solution in this situation, nor if such condition, if any, may pick up our equal area law solution.

Second, let us consider improvements of the approximations. In this work, we have applied the large-$N$ approximation. For beyond this approximation, the beta function for the effective potential generally involves its second-order derivative terms with respect to the field $\sigma$, namely $F\!\left(M,M^{\prime},\sigma;t\right)$ where the prime on the mass function denotes the $\sigma$-derivative and $M=\partial_{\sigma}V$; see Eqs. (66) and (67). In this case, we cannot naively transfer the $t$- and $\sigma$-derivatives from the mass function to the test function by performing the integration by parts as shown in Eq. (28) and then cannot define the weak solution. The notion of the weak solution that we have introduced in this work corresponds to the so-called “entropy solution” based on distributions. On the other hand, it is known that there is also “viscosity solution” as a different definition of the weak solution.151515 See e.g., Evans (2010); Crandall et al. (1992) for descriptions on their mathematical definitions.

The viscosity solution can be defined for even some second-order nonlinear partial differential equations to which the definition of Eq. (28) cannot be applied. The viscosity solution could be a breakthrough for improving the approximation. Therefore, it is important to study applications of the weak solution method in order to analyze critical phenomena in elementary particle physics and condensed matter physics.

Acknowledgements

We thank Yasuhiro Fujii and Daisuke Sato for valuable discussions. The work of K-I. A. is supported by JSPS Grant-in-Aid for Challenging Exploratory Research (No. 16K13848). The work of M. Y. is supported by the DFG Collaborative Research Centre SFB 1225 (ISOQUANT).

Appendix A Derivation of beta function

In this appendix, we derive the RG equation for the effective potential $V\!\left(\bar{\psi}\psi,\Lambda\right)$ in both cases of the sharp cutoff and the smooth cutoff schemes. We derive the RG equation using the WH equation with the sharp cutoff and the Wetterich equation with the smooth cutoff.

A.1 Sharp cutoff scheme

For the given effective action (14), the WH equation reads

[TABLE]

where $\phi^{\rm T}_{p}=\left(\psi^{\rm T}\!\left(-p\right),\bar{\psi}\!\left(p\right)\right)$. To derive the RG equation, we introduce the mean-field defined as

[TABLE]

We apply them after evaluating the functional derivatives of $S_{\rm eff}$ with respect to $\phi\!\left(p\right)$. The inverse two point function is given by

[TABLE]

where we can calculate

[TABLE]

Here we employ the following approximation:

[TABLE]

which corresponds to the so-called large-$N$ leading approximation ($N_{\text{c}}\to\infty$).

The matrix (65) is reduced to the simple form,

[TABLE]

where we defined {\ooalign{\hfil/\hfil\crcrp}}^{\pm}:={\ooalign{\hfil/\hfil\crcrp}}\pm i\mu\gamma_{0}. After using the formula of matrix, ${\rm tr}\,\log\!\left(A\right)=\log{\rm det}\,\!\left(A\right)$, and {\rm tr}\log\!\left({\ooalign{\hfil/\hfil\crcrp}}^{\pm}\mp i\partial_{\sigma}V\right)=2\log\!\left(\left(p^{\pm}\right)^{2}+(\partial_{\sigma}V)^{2}\right), the shell momentum integral for Eq. (63) is evaluated as follows:

[TABLE]

where $E:=\sqrt{\Lambda^{2}+(\partial_{\sigma}V)^{2}}$ and $\Omega:=(2\pi)^{4}\delta^{(4)}\!\left(0\right)$ is a volume of space-time. Here $C_{n}:=2\log\!\left(T^{2}+p_{0}^{2}\right)$ does not depend on the effective potential and hereafter it is ignored.

Next, we evaluate the Matsubara summation using the formula,

[TABLE]

The summation for Eq. (70) with $p_{0}=(2n+1)\pi T$ becomes

[TABLE]

where $C:=-2T+2T\log\!\left(1+{\rm e}^{-1}\right)$ is a constant value and can be ignored. To summarize, using $\partial_{t}V=-\partial_{t}S_{\rm eff}/\Omega$, we obtain the RG equation for the effective potential as follows:

[TABLE]

A.2 Smooth cutoff scheme

To obtain the RG equation in case of the smooth cutoff scheme, we use the Wetterich equation Wetterich (1993); Morris (1994),

[TABLE]

where the $\Gamma_{\Lambda}$ is the Legendre effective action and $\mathcal{R}_{\Lambda}$ is the cutoff profile function. See e.g. Gies (2012); Oda and Yamada (2016) for its explicit derivation.

For the level of approximation at present, the Legendre effective action is equivalent to the Wilsonian effective action Morris (1994), thus it is given by

[TABLE]

Therefore the two point function $\frac{\overrightarrow{\delta}}{\delta\phi_{p}}\Gamma_{\Lambda}\frac{\overleftarrow{\delta}}{\delta\phi_{-p}}$ becomes the same form as Eq. (69). As mentioned above Eq. (18), since we insert the cutoff scale $\Lambda$ into the momentum in the spatial direction, the cutoff profile function is given by

[TABLE]

We see

[TABLE]

where we defined p^{\pm}:=\left(p_{0}\pm i\mu,\,\vec{{\ooalign{\hfil/\hfil\crcrp}}}(1+r_{\Lambda})\right), and used the large-$N$ approximation. Here we employ the optimized cutoff function Litim (2001); Litim and Pawlowski (2006) as the smooth cutoff scheme, namely,

[TABLE]

Let us now compute the beta function for the effective potential. For the system (77) at finite temperature and density, the Wetterich equation (74) reads

[TABLE]

where $\phi^{\rm T}\!\left(p\right)=\left(\psi^{\rm T}\!\left(-p\right),\bar{\psi}\!\left(p\right)\right)$, and we introduced the derivative with respect to $t$, which acts on only the cutoff dependence of $\mathcal{R}_{\Lambda}$, thus, it is $\tilde{\partial}_{t}:=\frac{\text{d}\mathcal{R}_{\Lambda}}{\text{d}t}\cdot\frac{\partial}{\partial\mathcal{R}_{\Lambda}}=\frac{\partial r_{\Lambda}}{\partial t}\cdot\frac{\partial}{\partial r_{\Lambda}}$. To compare this form with the WH equation (63), we notice that both equations are almost the same structure and that $\tilde{\partial}_{t}$ plays the role of taking the shell momentum mode. Thereby, we can use the same calculations as the case of the sharp cutoff scheme, and thus we evaluate the RHS of Eq. (81) as follows:

[TABLE]

where we used the second line of Eq. (71) for the Matsubara summation. Then we obtain the RG equation,

[TABLE]

Appendix B Concavity and convexity of beta function

We see that the beta function for the fermionic effective potential is concavity. To this end, we first consider the second-order derivative of the beta function for the sharp cutoff regularization case (Eq. (73)) with respect to $M$. Here we define $N_{-}:={\rm e}^{\beta(E-\mu)}$ and $N_{+}:={\rm e}^{\beta(E+\mu)}$. Then the beta function is

[TABLE]

We have

[TABLE]

Since

[TABLE]

we find $f\!\left(M;t\right)>0$ and then $\partial^{2}F\!\left(M;t\right)/\partial M^{2}=-\Lambda^{3}f\!\left(M;t\right)/\pi^{2}<0$. Hence, the beta function at a mass function for arbitrary temperature and density is concave as a function of $M$.

Next, the case for the optimized cutoff regularization is considered. The beta function is

[TABLE]

Its second-order derivative with respect to $M$ is

[TABLE]

Since one cannot analytically see its convexity and concavity, we plot it as a function of $M$ in Fig. 13, where the dimensionless $\partial^{2}\bar{H}\!\left({\bar{M}};t\right)/\partial{\bar{M}}^{2}=(\partial^{2}H\!\left({M}/\Lambda;t\right)/\partial M^{2})/\Lambda^{2}$ is plotted. We see that the concavity of the beta function in the optimized cutoff regularization case is not guaranteed since Eq. (88) can have positive values.