RG scaling relations at chiral phase transition in two-flavor QCD

TL;DR

This study examines the chiral phase transition in two-flavor QCD, using lattice simulations and renormalization group analysis, finding evidence supporting a second order transition through universal scaling behavior of meson masses.

Contribution

It derives and tests a scaling relation for meson masses at the chiral transition, providing numerical evidence for a second order phase transition in two-flavor QCD.

Findings

Effective meson masses follow universal scaling curves.

Numerical data aligns with second order phase transition predictions.

Scaling behavior is consistent across different lattice sizes.

Abstract

We investigate the nature of the chiral phase transition in the massless two-flavor QCD using the renormalization group improved gauge action and the Wilson quark action on $32^3\times 16$, $24^3\times 12$, and $16^3\times 8$ lattices. Based on the renormalization group equation, we derive the scaling relation for the effective masses of mesons at the chiral phase transition point. If the chiral phase transition is second order, the effective masses as a function of the rescaled time/space do not depend on the lattice size and show the universal behavior. We find that our numerical simulations on the three sizes of lattices are excellently on the scaling curves, which is consistent with the second order phase transition.

Figures

Tables

Table 1. Table 1: Simulation parameters at massless quark points K c subscript 𝐾 𝑐 K_{c} on the 32 3 × 16 superscript 32 3 16 32^{3}\times 16 , 24 3 × 12 superscript 24 3 12 24^{3}\times 12 and 16 3 × 8 superscript 16 3 8 16^{3}\times 8 lattices and numerical results for m q subscript 𝑚 𝑞 m_{q} , m P ​ S subscript 𝑚 𝑃 𝑆 m_{PS} , m V subscript 𝑚 𝑉 m_{V} and the chiral condensate ⟨ ψ ¯ ​ ψ ⟩ delimited-⟨⟩ ¯ 𝜓 𝜓 \langle\bar{\psi}\psi\rangle . ”traj.” represents the number of trajectories (including thermalization) and ”acc.”represents the acceptance ratio of the HMC Metropolis test. The ”NAN” represents that the simulation cannot be performed due to the zero eigenvalues. See the main text in section IV.

size	$\beta$	$K$	traj.	acc.	$m_q$	$m_{PS}$	$m_V$	$⟨\overline{\psi }\psi ⟩$
16x32	3.0	.1435	2000	0.84(5)	.0087(2)	.374(2)	.408(3)	0.0087(15)
16x32	2.9	.1445	2000	0.76(4)	.0021(3)	.379(3)	.427(7)	0.0619(7)
16x32	2.8	.1455	3000	0.69(2)	.0052(3)	.365(3)	.429(3)	0.0169(14)
16x32	2.7	NAN
12x24	2.7	.1475	2000	0.87(3)	-.015(1)	.494(10)	.575(2)	-0.060(1)
12x24	2.6	.148	3000	0.83(1)	.0091(4)	.481(50	.668(6)	0.032(10)
12x24	2.5	NAN
8x16	2.5	.150	2500	0.94(1)	.003(1)	.735(11)	.821(15)	-0.0033(44)
8x16	2.4	.152	3000	0.94(2)	-.001(1)	.718(6)	.822(8)	-0.0033(44)
8x16	2.3	.1547	6000	0.88(1)	-.009(6)	.680(5)	.824(6)	-0.0393(20)
8x16	2.2	NAN

Table 2. Table 2: Numerical results for m q subscript 𝑚 𝑞 m_{q} , m P ​ S subscript 𝑚 𝑃 𝑆 m_{PS} , and m V subscript 𝑚 𝑉 m_{V} on the 16 3 × 64 superscript 16 3 64 16^{3}\times 64 lattice at β = 2.3 𝛽 2.3 \beta=2.3 for determination of the lattice spacing

$K$	$m_q$	$m_{PS}$	$m_V$
.140	.372(1)	1.113(1)	1.148(1)
.145	.234(1)	0.879(1)	0.931(2)
.146	.209(1)	0.826(1)	0.879(2)
.147	.184(1)	0.769(1)	0.828(2)
.148	.160(1)	0.716(1)	0.781(1)
.149	.134(1)	0.656(2)	0.729(4)
.150	.111(1)	0.591(1)	0.676(3)
.151	.089(1)	0.523 (2)	0.609(5)
.152	.067(1)	0.448(3)	0.556(12)

Equations

$$m_q$$

$$=\underset{t\to N_t/2}{\lim }\frac{{\sum }_x⟨{\nabla }_4{A}_4(x,t)P(0)⟩}{2{\sum }_x⟨P(x,t)P(0)⟩}$$

$$G_t(t)=\sum _{x,y,z}⟨\bar{\psi }{\gamma }_H\psi (x,y,z,t)\bar{\psi }{\gamma }_H\psi (0,0,0,0)$$

$$G_s(x)=\sum _{t,y,z}⟨\bar{\psi }{\gamma }_H\psi (x,y,z,t)\bar{\psi }{\gamma }_H\psi (0,0,0,0)$$

$$\frac{\cosh (m(t)(t-N_t/2))}{\cosh (m(t)(t+1-N_t/2))}=\frac{G(t)}{G(t+1)}.$$

$$m(t)=\ln \frac{G(t)}{G(t+1)}.$$

$$G_t(n_t;g,m_q,N,\mu {)}_H$$

$${\left(\frac{N^{\prime }}{N}\right)}^{-2\gamma }$$

$$G_s(n_s;g,m_q,N,\mu {)}_H$$

$${\left(\frac{N^{\prime }}{N}\right)}^{-2\gamma }$$

$$\frac{dg}{d\log s}=\mathcal{B}$$

$$\frac{d{m}_q}{d\log s}=\gamma m_q.$$

$$\mathfrak{m}(n_t;g,m_q,N)=N\ln \frac{G(n_t;g,m_q,N)}{G(n_t+1;g,m_q,N)},$$

$$\mathfrak{m}(\tau ;g,m_q,N)=-{\partial }_{\tau }\ln G(\tau ;g,m_q,N).$$

$$\mathfrak{m}(\tau ;g,m_q,N)=\mathfrak{m}(\tau ;g^{\prime },m_q^{\prime },N^{\prime })$$

$$\tilde{G}(\tau ,N)=G(\tau ;g^{\ast },m_q=0,N).$$

$$\tilde{G}(\tau ;N)={\left(\frac{N{}^{\prime }}{N}\right)}^{-2{\gamma }^{\ast }}\tilde{G}(\tau ;N{}^{\prime }).$$

$$\mathfrak{m}(\tau ,N)=\mathfrak{m}(\tau ,N{}^{\prime })$$

$$\check{G}(\tau ;g(N),N)=G(t;g(N),m_q=0,N).$$

$$\check{G}(\tau ;g(N),N)={\left(\frac{N{}^{\prime }}{N}\right)}^{-2\gamma }\check{G}(\tau ;g(N{}^{\prime }),N{}^{\prime }).$$

$$\check{\mathfrak{m}}(\tau ,g(N),N)=-{\partial }_{\tau }\ln \check{G}(\tau ,g(N),N)$$

$${\check{\mathfrak{m}}}_t(\tau ,g(N),N)={\check{\mathfrak{m}}}_t(\tau ,g(N{}^{\prime }),N{}^{\prime })$$

$${\check{\mathfrak{m}}}_s(\tau ,g(N),N)={\check{\mathfrak{m}}}_s(\tau ,g(N{}^{\prime }),N{}^{\prime }).$$

$${\check{\mathfrak{m}}}_{\mathfrak{s}}(\tau ,g(N),N)={\mathfrak{F}}_{\mathfrak{s}}(\tau ),$$

$$⟨\bar{\psi }\psi {⟩}_{\text{sub}}=2m_{q}a(2K{)}^2\sum _t{G}_{\text{PS}}(0,t),$$

$$\underset{m\to 0}{\lim }⟨\bar{\psi }\psi (m)⟩\simeq m^{\eta }$$

$$\underset{\lambda \to 0}{\lim }\rho (\lambda )\simeq {\lambda }^{\eta }$$

$$6/g_1^2-6/g_2^2=\hat{\mathfrak{B}}(N_1,N_2)\ln ({\mu }_1/{\mu }_2)$$

$$\hat{\mathfrak{B}}(32,16)=0.72$$

$$\hat{\mathfrak{B}}(32,24)=0.69$$

$$\hat{\mathfrak{B}}(24,16)=0.74$$

$$\hat{\mathfrak{B}}(\ln ({\mu }_1/{\mu }_2)+\ln ({\mu }_2/{\mu }_3)+\ln ({\mu }_3/{\mu }_1))=0$$

Full text

††preprint: APS/123-QED

I introduction

The properties of quarks and gluons at high temperatures and the nature of the transition from the quark-gluon state to the hadronic state are key ingredients for understanding the evolution of the Universe. Yet, even the order of the chiral phase transition is still under debate. To make the discussions more transparent, let us consider the idealistic situation in which we have two massless quarks in the SU(3) gauge theory, and ask the question whether the finite temperature phase transition of the chiral symmetry is first order or second order. Pisarski:1983ms - Nakayama:2014sba

Lattice QCD is probably the most reliable constructive formulation for the investigation of non-perturbative properties of quarks and gluons at the chiral phase transition point (see e.g. review1 for a review). In this article we investigate the chiral phase transition in the massless two-flavor QCD, using the renormalization group (RG) improved gauge action and the Wilson quark action with two degenerate quarks on $32^{3}\times 16$, $24^{3}\times 12$, and $16^{3}\times 8$ lattices. The simulation itself may be a standard one, but we are going to employ a new method based on the renormalization group to investigate the long distance physics of the chiral phase transition.

The key idea is to test the RG scaling relation of the meson propagator. Based on Wilson’s idea of RG, we derive the RG equation for the effective masses of mesons. The RG equation we derive contains only two relevant (including marginal relevant) operators (i.e. the gauge coupling and the quark mass), and the RG equation has an interesting scaling solution at the chiral phase transition if it is second order. Alternatively, when the chiral phase transition is first order, then we do not expect the existence of such a solution, or even the RG equation as it is may not hold because of the emergence of other relevant operators.

We show that effective masses calculated on the three lattices for the pseudo scalar (PS) mesons and vector (V) mesons agree with the predicted scaling relations in an excellent manner. While we cannot exclude the possibility that the behavior changes at the larger lattice size or in the thermodynamic limit, we argue that the numerically verified scaling behavior supports the claim that the chiral phase transition in the massless two-flavor QCD is second order.

The organization of the paper is as follows. After describing our setup in section 2, we revisit RG equations in section 3 and derive RG scaling relations in section 4. After giving the job parameters in section 5, we identify the chiral phase transition points in section 6. We show our numerical results with the verification of the scaling relation in section 7. The characteristic of the universal scaling curves is discussed in section 8. In section 9, the chiral phase transition temperature is estimated. Finally summary and discussion are given in section.10.

II Action and Observables

We define continuous gauge theories as the continuum limit of lattice gauge theories, defined on the Euclidean lattice of the size $N_{x}=N_{y}=N_{z}=N_{s}$ and $N_{t}$. We impose an anti-periodic boundary condition in the temporal direction for fermion fields and periodic boundary conditions otherwise.

The continuum limit of QCD is taken at or toward a UV critical point (in this article we consider only the case where the UV fixed point is $g_{0}=0$ and $m_{0}=0$), with the lattice spacing $a\rightarrow 0$ together with $N_{s}\rightarrow\infty$ and $N_{t}\rightarrow\infty,$ changing the coupling $g(a)$ in a such way that physical quantities with physical dimension such as hadron masses $m_{H}/a$ fixed. We may use the aspect ratio $r=N_{s}/N_{t}$ and $N=N_{t}$ instead of $N_{s}$ and $N_{t}.$ In this article, we keep the aspect ratio $r$ fixed when we change the lattice size. The necessary change of the coupling constant is governed by the renormalization group equation.

When $N_{s}$ and $N_{t}$ are finite, the formulation of finite temperature QCD on a lattice is equivalent to a Euclidean path integral defined on a discrete three dimensional lattice cite with a transfer matrix for a discrete time. However, at finite temperature we have to take the thermodynamic limit $r\rightarrow\infty$ for interpretation of physical quantities at finite temperature.

Although our methodology can be applied to any gauge theories with fermions in arbitrary representation, we focus on $SU(3)$ gauge theories with degenerate $N_{f}=2$ fundamental fermions (“quarks”) in this article. We employ the Wilson quark action and the RG improved gauge action RG-improved . The theory is defined by two parameters; the bare coupling constant $g_{0}$ and the bare degenerate quark mass $m_{0}$ at ultraviolet (UV) cutoff. We also use, instead of $g_{0}$ and $m_{0}$, $\beta={6}/{g_{0}^{2}}$ and the hopping parameter $K=1/2(m_{0}a+4)$.

We measure the plaquette and the Polyakov loop in each space-time direction, the mass of hadrons such as the pseudo-scalar meson mass $m_{\text{PS}}$. The quark mass $m_{q}$ is defined Bochicchio:1985xa ,Itoh:1986gy .as the large $t$ value of $m_{q}(t)$ obtained through Ward-Takahashi identities by the ratio of thermal propagators=

[TABLE]

where $P(x,t)$ is the pseudo-scalar density and $A_{4}(x,t)$ the fourth component of the local axial vector current, renormalization constants being suppressed. This is also the same definition that we use in the zero-temperature lattice QCD. The quark mass $m_{q}$ thus defined does not depend on whether the system is confining or deconfining, and depends on only $\beta$ and $K$ up to order $1/N_{s}$ and $1/N_{t}$ corrections.

The propagator is defined by

[TABLE]

for the temporal one and

[TABLE]

for the spatial one. $H$ here stands for the channel, and in this article, we focus on the pseudo-scalar $(H=PS)$ and the vector channel $(H=V)$, and we sometimes omit the subscript $H$ in the following.

We define the effective mass $m(t)$ (or similarly $m(x)$) through

[TABLE]

which reduces to

[TABLE]

when boundary effects can be neglected. In the case of exponential-type decay the effective mass approaches a constant value in the large $t$ regime, which we call a plateau.

In our earlier studies Refs. Ishikawa:2013wf - Ishikawa:2013tua , we argued that there is a region where the propagator exhibits a power-law corrected Yukawa type decay when the quark mass is small on the $16^{3}\times 64$ lattice (with aspect ratio $r=N_{s}/N_{t}=1/4$). We observed the similar behavior for the aspect ratio $r=2$ which is larger than unity in Ishikawa:2015nox . In this article, we mainly discuss the $r=2$ case with the thermodynamic limit in mind.

III RG equations

Let us recall the RG equations on a finite lattice, which represent Wilson’s basic idea Wilson:1973jj . In the vicinity of a critical point where the correlation length diverges, short-range fluctuations can be integrated out while keeping the long distance behavior. When we change the size of the lattice, the action must be changed to keep the physics intact, but the key observation by Wilson is that near a critical point, it is sufficient to take into account relevant (including marginal relevant) operators.

Now let us consider the case when the chiral phase transition is second order. The RG equation in the vicinity of the critical point is given  DelDebbio:2010ze - DelDebbio:2010hu

[TABLE]

for the temporal propagator and

[TABLE]

for the spatial propagator. The suffix $s$ and $t$ represents spatial and temporal, respectively. In the following we often suppress the suffix except when it is better to specify which we refer to.

Here $n=(n_{x},n_{y},n_{z},n_{t})$, and $\mu{{}^{\prime}}=\mu/s$ and $N_{s}{{}^{\prime}}=N_{s}/s,$ $N_{t}{{}^{\prime}}=N_{t}/s,$ and $n{{}^{\prime}}=n/s$ with $s$ being the change of the scale under the renormalization. The UV renormalization scale $\mu$ in lattice theories is set by the inverse lattice spacing $a^{-1}.$ Note $N_{s}a=L_{s}$ and $N_{t}a=L_{t}$ are kept constant. The relation between $g^{\prime}$ and $g$ and $m^{\prime}_{q}$ and $m_{q}$ are determined by the RG beta function $\mathcal{B}$ and the mass anomalous dimension $\gamma$:

[TABLE]

and

[TABLE]

To make the expression simpler, it is assumed that $\gamma$ does not change very much as the ratio $N^{\prime}/N$ changes, which is always valid when we are at the RG fixed point, or when $N^{\prime}/N\sim 1$ even if we are not at the fixed point. With these caveats in mind, we are going to verify the RG equation numerically.

As a technical remark, we note that the RG equations (7) and (8) for the temporal propagator and for the spatial propagator define the same RG beta function and the mass anomalous dimension, since both are derived from the RG equation for the common two-point functions. However, the RG beta function $\mathcal{B}$ and the anomalous dimension $\gamma$ may depend on the aspect ratio $r$ and in addition the latter takes a different value for a different channel.

IV RG Scaling relations

Our goal is to test the RG equation to study a long distance behavior of the meson propagator. For this purpose, let us rewrite the RG equation in a more convenient way. We first define the scaled effective mass by

[TABLE]

suppressing $\mu$. It reduces in the continuum limit $N\rightarrow\infty$ to the form

[TABLE]

Here $\tau=n_{t}/N_{t}$. The variable $n_{t}$ takes $0,1,2,\cdots,N_{t}-1$ so that the scaled time always lie in the range $0\leq\tau\leq 1$.

For the scaled effective mass, the RG equation takes the form

[TABLE]

IV.1 Solution at an IR fixed point

The RG equations (7) and (8) derived at the chiral phase transition point are also applied to the IR fixed point of conformal QCD at zero temperature. Let us consider when the theory is in the conformal window ($N_{f}*<N_{f}<16.5$) at zero temperature. Suppose we are at an IR fixed point, then we have the fixed parameters $g^{\prime}=g=g^{*}$ and $m_{q}^{\prime}=m_{q}=0$ so that $\mathcal{B}=0$ and $\gamma=\gamma^{*}$ in the RG equation Eq.(7). Hence, the propagator can be expressed by a simplified notation suppressing $g=g^{*},m_{q}=0,$

[TABLE]

and the RG relation eq.(7) reduces to

[TABLE]

Combining Eqs.(12) and (15), we derive the RG scaling relations for effective masses

[TABLE]

at the IR fixed point. It means that the scaled effective mass does not depend on $N$ as a function of $\tau$. Therefore, the agreement of the scaled effective mass as a function of $\tau$ is a stringent test of the fixed point. We applied in Ref.Ishikawa:2015iwa the scaling relation (16) in order to identify the IR fixed point for $N_{f}=16,12,8,7$ in $SU(3)$ gauge theories.

IV.2 Solution at chiral phase transition

Let us assume that we are at a chiral phase transition point. The chiral phase transition point depends on the lattice size $g(N)$, where the precise criterion to call the chiral phase transition point on a finite lattice will be discussed in section  VI. Using the simplified expression

[TABLE]

Eq. (7) reduces to

[TABLE]

Correspondingly, the scaled effective mass becomes

[TABLE]

as in the previous subsection, and we obtain the RG equation

[TABLE]

and

[TABLE]

Eqs. (20) and (21) are key scaling relations of this article, which are valid when the chiral phase transition is second order. Since the number of the parameters is reduced by setting $g=g(N)$, one may solve the RG equation as

[TABLE]

where $\mathfrak{F}_{s}(\tau)$ is renormalization invariant.

The distinct point of Eqs. (20) and (21) from Eq.(16) is that the coupling $g$ is fixed in the former, while $g$ moves according to the RG beta function in the latter. Nevertheless, the both (16) and (22) show the universal scaling behavior irrespective of $N$. Furthermore, while the RG equation does not say anything about the shape of the scaling curve $\mathfrak{F_{s}(\tau)}$, we will see the numerical result for scaled effective masses $\mathfrak{F_{s}(\tau)}$ at the chiral phase transition shows the very similar behavior to the one at the zero temperature in the conformal QCD.

When the RG equation is evaluated in the vicinity of the UV fixed point $g_{0}=0$ and $m_{0}=0$, the quark mass term is relevant and the gauge coupling is marginal. Along the RG trajectory from the UV fixed point to the IR critical point, the beta function possesses a zero at a certain value $g=g^{*}$ in the conformal QCD, at which the quark mass term is relevant and the gauge coupling is irrelevant. In the case of the chiral phase transition, it does not possess a zero, that is, the beta function is negative along the RG trajectory. The gauge coupling constant (or temperature) is relevant, and we have to tune it to obtain the criticality. More generally, we may consider an infinite parameter space of possible interactions $H=\sum g_{i}O_{i}$ and we may define the beta function $\beta_{i}(\mu)$ for each coupling $g_{i}.$ By RG transformation the effective Hamiltonian possess many possible interactions like a four-fermion interaction. These operators are irrelevant operators at the UV fixed point and retain irrelevant at the chiral phase transition point. We have implicitly assumed that the number of relevant operators does not increase along the RG trajectory, which was at the core of the derivation of our RG equation.

In contrast, let us see what would happen if the chiral phase transition were the first order phase transition. Then, typically, we would expect an extra relevant operator that will induce the chiral symmetry breaking. The search for such an operator in the effective Landau-Ginzburg model was a key diagnosis for determining the order of phase transition in the study of Pisarski and Wilczek Pisarski:1983ms . Accordingly, if this scenario is the mechanism of the first order phase transition, the assumption of our RG scaling relation is invalidated, and there is no generic reason to believe that the scaling relation that we proposed is observed.

V Job parameters

We perform simulations with two degenerate quarks on $32^{3}\times 16,$ $24^{3}\times 12$ and $16^{3}\times 8$ lattices to investigate the scaling of the effective masses of mesons. The algorithm we employ is the blocked HMC algorithm Hayakawa:2010gm . We choose the run-parameters in such a way that the acceptance of the HMC Metropolis test is about $70\%\sim 90\%.$ The statistics are 1,000 MD trajectories for thermalization and $1000\sim 5000$ MD trajectories for the measurement. We estimate the errors by the jack-knife method with a bin size corresponding to 100 HMC trajectories (see Table LABEL:job_parameter).

VI Identification of chiral phase transition points

The first task in our numerical simulations is to determine the chiral phase transition point on the parameter space $K$ and $\beta$. For this purpose, let us discuss the nature of the chiral phase transition and its realization on a finite lattice. At high temperature the chiral symmetry is restored and the order parameter (i.e. chiral condensate) vanishes: $\,\langle\bar{\psi}\psi\rangle=0,$ while below the chiral phase transition temperature $\langle\bar{\psi}\psi\rangle\neq 0$. In the Wilson fermion formalism a properly subtracted $\langle\bar{\psi}\psi\rangle$ can be defined via an axial Ward-Takahashi identity Bochicchio:1985xa

[TABLE]

where $G_{\text{PS}}(0,t)$ is the propagator of the pseudo-scalar meson.

A direct approach to identify the (would-be) chiral phase transition points in lattice simulations is to check whether the oder parameter of the chiral symmetry vanishes or does not. We should note, however, that there is no real phase transition on a finite lattice in this strict sense. The long-time average of order parameter is always zero in finite systems, and if we try to practically measure it over a finite time average, the change of the order parameter around the (would-be) phase transition will be rounded.

Rather than directly measuring the chiral condensate, we take the following alternative approach. Let us recall that the Banks-Casher relation Banks:1979yr DelDebbio:2010ze states that the critical exponent $\eta$

[TABLE]

is identical with the critical exponent of eigenvalue density of the massless Wilson-Dirac operator.

[TABLE]

Based on this fact, we proposed the “on Kc method” Iwasaki:1996zt , which monitors the number of iteration of CG inversion” to determine the chiral phase transition points. In the chiral symmetry broken phase, $\eta=0$ and it implies that the zero eigenvalues accumulate. Hence it is impossible to invert the Dirac-Wilson operator in the chiral symmetry broken phase, while in the symmetric phase there are no theoretical obstruction for the inversion albeit many iterations may be needed in numerical methods. We, therefore, proposed to identify the chiral phase transition as the point where we cannot invert the Dirac-Wilson operator to continue the numerical simulation.

Let us summarize our strategy in our setup. We first determine the hopping parameter $K_{c}$ for the massless quark from $\beta=2.3$ to $3.0$ (Table LABEL:job_parameter) Then we estimate the chiral phase transition points on the lattices of the three sizes: $32^{3}\times 16$, $24^{3}\times 12$ and $16^{3}\times 8$ by using this “on Kc method”,

As mentioned above, on a finite lattice the phase transition is rounded and the change of the number of iterations is not as drastic as it would be in the thermodynamic limit. Therefore there is certain arbitrariness in the definition of the phase transition on a finite lattice. We adopt the following as our definition of the (would-be) chiral phase transition point. We begin with simulations on the massless lines of the $16^{3}\times 8,$ $24^{3}\times 24$ and $32^{3}\times 16$ lattices at some $\beta$ (we take $\beta=3.0$) which is apparently in the chiral symmetric phase and decrease $\beta$ by a step $\Delta\beta=-0.1$. As $\beta$ decreases, the number of the iterations for the inversion of the Wilson-Dirac operator in molecular dynamics steps increases and the acceptance ratio decreases. Finally at some points we cannot continue the simulations without decreasing the molecular steps.

We identify such $\beta$’s as the chiral phase transition points. Our simulation shows that these three points are:

• •

$\beta_{*}\simeq 2.8;K_{*}=0.1455$ on the $32^{3}\times 16$ lattice;

• •

$\beta_{*}\simeq 2.6;K_{*}=0.1480$ on the $24^{3}\times 12$ lattice;

• •

$\beta_{*}\simeq 2.3;K_{*}=0.1547$ on the $16^{3}\times 8$ lattice.

To double-check our criterion, we have also computed the chiral condensates on the $K_{c}$ in the chiral symmetric phase as given in Table 1: They are of order $O(10^{-2})$ in lattice unit. On the other hand, in the confining phase we are unable to compute the order parameter directly at the massless point. For reference we quote some of them at small (not very small) quark mass: At $\beta=2.7$, $K=0145$ on the $32^{3}\times 16$ lattice, where $m_{q}=0.049(1)$, we have $\,\langle\bar{\psi}\psi\rangle=0.161(1)$. At $\beta=2,5$, $K=0.149$ on the $24^{3}\times 12$ lattice, where $m_{q}=0.031(1)$, we have $\,\langle\bar{\psi}\psi\rangle=0.103(1)$. At $\beta=2.2$, $K=0.155$ on the $8^{3}\times 16$ lattice, where $m_{q}=0.068(1),$ we have $\,\langle\bar{\psi}\psi\rangle=0.296(5).$ They are of order $O(10^{-1})$ and are larger than those in the symmetric phase by a factor of ten. This is consistent with our identification of the (would-be) chiral phase transition points.

Now we introduce a discrete beta function $\mathfrak{\hat{B}}(N_{1},N_{2})$ defined by

[TABLE]

which relates the critical values $\beta^{*}$ to each other. Note $N_{i}=L/a_{i}$ and $\mu_{i}=1/a_{i}.$ Therefore $\ln{(\mu_{1}/\mu_{2})}=\ln{(N_{1}/N_{2})}.$ This definition may be more convenient because the one-loop result shows that $\mathfrak{\hat{B}}$ is just a constant independent of $g$. From our results on the chiral phase transition point, we obtain

[TABLE]

[TABLE]

[TABLE]

Thus the beta function $\mathfrak{\hat{B}}(N_{1},N_{2})$ is roughly a constant $0.7.$ It means the three equations are mutually consistent. The sum of the left side of the three equations (26) vanishes and therefore the sum of the right sides also should vanish. If we can regard $\mathfrak{\hat{B}}$ is a constant around here, the consistency is automatically satisfied

[TABLE]

Converting the result to the standard one,

[TABLE]

we obtain the beta function $\mathfrak{B}=-0.183$ at $\beta=2.8$ and $-0.246$ at $\beta=2.3$. We curiously note that this number is close to the one-loop beta function of $\mathcal{B}=-\frac{33-2N_{f}}{48\pi^{2}}g^{3}$.

VII Numerical results for RG scaling

VII.1 spatial effective masses

Now let us show numerical results of the spatial effective masses measured at the critical points in Fig.1. We plot the scaled spatial effective masses defined in eq.(11) in terms of $\tau=n_{s}/N_{s}$ in order to test the scaling relation (21). We overlay the data on the three lattices of $32^{3}\times 16$, $24^{3}\times 12$ and $16^{3}\times 8$.

On the left panel we show the effective masses for the pseudo-scalar channel, while on the right panel we show those for the vector channel. We see that all data are excellently on the scaling curve except for three points at short distance ($n_{s}=0,1,2$) on each of the lattices. In particular, the scaling curve for both pseudo-scalar channel and the vector channel on the $32^{3}\times 16$ and $24^{3}\times 12$ are within one-standard deviation around $\tau=0.4$. The scaling curve for the vector channel on the $16^{3}\times 8$ lattice is also almost within one standard deviation, while that for the pseudo-scalar channel on the $16^{3}\times 8$ lattice is several standard deviations off the other two values. We interpret that the $N=16$ is slightly small to regard it as $N\rightarrow\infty$, in which limit the RG scaling relation holds strictly.

The deviation in short distance behavior is expected to be due to the finite lattice spacing effect. Since it probes the UV behavior, the short distance physics is more sensitive to $1/a$ corrections. We expect that the data at short distance will envelop the limiting scaling curve with increasing $N_{t}$.

VII.2 temporal effective masses

Next, let us discuss the temporal effective mass $m_{t}(n_{t})$. The available data points are half of the spatial effective masses. Even in the case of the largest lattice we take, $32^{3}\times 16$, the number of the data point is only $8$. For the other cases they are $6$ and $4$ on $24^{3}t\times 12$ and $16^{3}\times 8$ lattices. Furthermore since the three points at short distance are affected by the lattice spacing effects, the effective data points we can use are $5$, $3$ and $1$, respectively. Therefore we do not expect to see a clear scaling.

The data are shown in Fig.(2). On the left panel we present the spatial effective masses of the pseudo-scalar channel, while on the right panel those of the vector channel. In the pseudo-scalar channel the difference between the $32^{3}\times 16$ lattice and the $24^{3}\times 12$ lattice is order of $1\%,$ which is order of one standard deviation. The difference between the $16^{3}\times 8$ lattice and the $32^{3}\times 16$ lattice is of order $5\%$, which is of order several standard deviations. In the vector channel the differences are of order $2.5\%$. Although the clear scaling is yet to be seen, we find the tendency that the larger $N$ curves will converge to the scaling curve.

VIII Screening mass and power-law corrected Yukawa type decay

As shown in the preceding section, we find that our numerical results on the three sizes of lattices are excellently on the scaling curves. In particular the data on the $32^{3}\times 16$ lattice well represents the universal curve $\mathfrak{F_{s}}(\tau)$. Note the spatial scaling curve corresponds to the screening mass. It is time to ask the physical meaning of the shape of the scaling curve $\mathfrak{F_{s}}(\tau)$.

The renormalization group equation itself does not determine the shape of the curve, so it should contain the more detailed physics of the chiral phase transition. Fig.1 shows that the effective masses on the $32^{3}\times 16$ lattice is monotonously decreasing and does not show a plateau up to the maximum point $n_{s}=15$ (with the periodic boundary condition). The asymptotic behavior of propagator $G(n_{s})$ in the limit $n_{s}\rightarrow 15$ is well represented by

[TABLE]

We fit the data (from $n_{s}=10$ to $15$) on the $32^{3}\times 16$ lattice to the power-law corrected Yukawa type form. The results of the fit are shown with the data in Fig. 3 for the pseudo-scalar channel on the left panel and the vector channel on the right panel, respectively. The best fitting parameters are $\hat{m}=0.314(6)$ and $\alpha=0.80(8)$ for the pseudo-scalar channel and $\hat{m}=0.351(11)$ and $\alpha=1.11(11)$ for the vector channel. We present them in terms of $\tau=n_{s}/N_{s}$ (x-axis) and $N_{s}\times m_{s}$ (y-axis).

If we fitted the data with a constant assuming a plateau at $n_{s}=14$ and $15$, we would get $m=0.366(1)$ and $m=0.429(14)$ for the pseudo-scalar channel and the vector channel respectively. These are about $15-20\%$ larger than the one estimated from the power-law corrected Yukawa type form.

In previous articles  Ishikawa:2013wf - Ishikawa:2013tua , we pointed out that the propagators show the power-law corrected Yukawa type decay at zero temperature in the $N_{f}=7,8,12$ and $16$ cases and at finite temperature $N_{f}=2$ case on a $16^{3}\times 64$ lattice. There, we argued that the power-corrections to the exponential decay may be originated from the unparticle nature of the (conformal) fixed point in agreement with the holographic analysis. Based on this theoretical picture, we claimed that the power-law corrected Yukawa-type decay observed not only for the special aspect ratio $r=1/4$, but for any aspect ratio. Our new results support this claim.

IX Critical phase transition temperature

Before we conclude, let us estimate the chiral phase transition temperature $T^{*}$ given by $1/(aN_{t})$. In order to fix the scale, a standard way we also apply is extrapolating the $\rho$ meson mass (in lattice units) to the chiral limit and comparing it with the experimental value. The spatial lattice size should be the same as that of the finite temperature for a given $\beta$ in a such way that the Hamiltonian in the three dimensional space is identical. While the lattice size $N_{t}$ should be large such that $r=N_{s}/N_{t}=1/4.$

Due to the CPU time constraint, we focus on the case with the $16^{3}\times 64$ lattice at $\beta=2.3$. The gauge coupling constant $\beta$ here is fixed by the chiral phase transition point for the $16^{3}\times 8$ lattice that we determined in section 6.

The simulations are done at nine hopping parameters from $K=0.140$ to $K=0.152$ (from $m_{q}=0.372$ to $m_{q}=0.067$) listed in Table LABEL:16x64_spectroscopy. We fit the rho meson mass in the linear way $m_{\rho}=am_{q}+c$ to obtain the vector meson mass in the chiral limit (more precisely; the physical point). The fit depends on the choice of the fit range of the quark mass. For example, when the fit range is from $K=0.145$ to $K=0.152$, the rho meson mass extrapolated takes $m_{\rho}=0.42(1)$ with $\chi^{2}/nf=4.4$. The fit with the data are shown in Fig.4. If we restrict the fit with the condition $\chi^{2}\leq 6.0$, $m_{\rho}=0.400\sim 0.445\,\,1/a.$ Using the experimental value $m_{\rho}=770$ MeV as an input, we obtain $1/a=1.73\sim 1.83$ GeV at $\beta=2.3.$ This implies that $a=0.104\sim 0.115$ fm and the lattice size in spatial direction is $L_{s}=1.66\sim 1.85$ fm. The transition temperature is estimated to be $T^{*}\simeq 215\sim 240$ MeV.

This value is larger than the value $T_{c}\sim 170-190$ MeV quoted in Karsch:1999vy . See also the more recent review Kanaya:2010qd and reference therein. However it is notoriously difficult to obtain the common value for the chiral phase transition temperature. It is important to perform a similar analysis on larger lattices and take the continuum limit carefully.

X summary and discussion

We have investigated the nature of the chiral phase transition in the degenerate two flavor QCD using the RG improved gauge action and the Wilson quark action on the $32^{3}\times 16$, $24^{3}\times 12$, and $16^{3}\times 8$ lattices.

We derived the RG scaling relations at the chiral phase transition point based on Wilson’s idea by assuming that the chiral phase transition is second order

[TABLE]

Our numerical results show that three effective spatial masses for the pseudo scalar mesons and vector mesons excellently scale. While the number of the data points are small in the temporal masses to see a clear scaling behavior, the data are consistent with the scaling. The RG scaling relations are derived with the assumption that the chiral phase transition is second order, and there is no reason to hold if the chiral phase transition is the first order. Our results therefore support the claim that the chiral phase transition is the second order. Nevertheless, we are not able to exclude a possibility that accidentally the correlation length is very large in a first order transition.

The RG approach is very powerful in investigating the properties of the critical points. It is our future mission to verify the RG scaling relations on larger lattices such as $48^{3}\times 24$ to see if the scaling relation persists in agreement with that the chiral phase transition is second order. In addition, it is important to derive critical exponents and compare them with the theoretically predicted ones to conclude the issue of the order of the chiral phase transition. In particular, the precise knowledge on the critical exponents will reveal whether the $U(1)_{A}$ symmetry is restored at the chiral phase transition point or not because the (non-)existence of the $U(1)_{A}$ symmetry affects the critical exponents.

We would also like to thank A. Ukawa for reading through the article and K. Kanaya for useful discussion. The calculations were performed on Hitachi SR16000 at KEK under its Large-Scale Simulation Program and HA-PACS computer at CCS, University of Tsukuba under HA-PACS Project for advanced interdisciplinary computational sciences by exa-scale computing technology.