Possible scenario of dynamical chiral symmetry breaking in the interacting instanton liquid model

We briefly explain the anomaly-driven breaking of chiral symmetry which was proposed by the previous work Kono:2021 based on the three flavor NJL model with the axial anomaly term. In Ref. Kono:2021 , the following Lagrangian is considered:

	$\displaystyle{\mathcal{L}}_{\rm NJL}$	$\displaystyle=$	$\displaystyle\bar{\psi}(i\gamma^{\mu}\partial_{\mu}-{\mathcal{M}})\psi$		(1)
			$\displaystyle+g_{S}\sum_{a=0}^{8}\left[\left(\bar{\psi}\frac{\lambda_{a}}{2}\psi\right)^{2}+\left(\bar{\psi}i\gamma_{5}\frac{\lambda_{a}}{2}\psi\right)^{2}\right]$
			$\displaystyle+\frac{g_{D}}{2}\left\{\det_{i,j}\left[\bar{\psi}_{i}(1-\gamma_{5})\psi_{j}\right]+{\rm H.c.}\right\},$

for the quark fields $\psi=(u,d,s)^{T}$, where ${\mathcal{M}}$ is the quark mass matrix given by ${\mathcal{M}}={\rm diag}(m_{q},m_{q},m_{s})$ with assuming isospin symmetry $m_{q}=m_{u}=m_{d}$, $\lambda_{a}\ (a=0,1,\dots,8)$ represent the Gell-Mann matrices in the flavor space with $\lambda_{0}=\sqrt{2/3}\bm{1}$, $\det_{i,j}$ is understood as determinant operation over the flavor indices of the quark fields $\psi_{i},\bar{\psi}_{j}\ (i,j=u,d,s)$, and $g_{S},g_{D}$ are the coupling constants for the four-point vertex interaction and the determinant-type ${\rm U}(1)_{A}$ breaking interaction, respectively. In the mean field approximation, the effective potential is obtained from the Lagrangian (1) as a function of the dynamical quark mass $M$ as follows

			$\displaystyle V_{\rm eff}(M)$		(2)
		$\displaystyle=$	$\displaystyle iN_{c}\sum_{f=q,q,s}\int\frac{d^{4}p}{(2\pi)^{4}}\mathrm{Tr}\left[\ln(\gamma\cdot p-M_{f})+\frac{\gamma\cdot p-m_{f}}{\gamma\cdot p-M_{f}}\right]$
			$\displaystyle-\frac{g_{S}}{2}(2\braket{\bar{q}q}^{2}+\braket{\bar{s}s})-g_{D}\braket{\bar{q}q}^{2}\braket{\bar{s}s},$

where $\braket{\bar{q}q}=\braket{\bar{u}u}=\braket{\bar{d}d}$ is the quark condensate with isospin symmetry and is given as a function of $M$. Their specific forms are given by Eq. (A.1) and Eq. (A.2) in Ref. Kono:2021 .

Following Ref. Kono:2021 , we start with the case of the chiral limit and no anomaly term, i.e., $m_{q}=m_{s}=0$ and $g_{D}=0$ in Eq. (1). In this situation, chiral symmetry is dynamically broken at the vacuum with a finite dynamical quark mass when the dimensionless coupling constant $G_{S}\equiv g_{S}/g_{S}^{\rm crit}$ is larger than $1$. Here, $g_{S}^{\rm crit}$ is the critical coupling constant defined by $g_{S}^{\rm crit}=2\pi^{2}/(3\Lambda_{3}^{2})$ with a three-momentum cutoff $\Lambda_{3}$ for the quark loop function. In other words, if there exists sufficiently strong four-quark interaction, chiral symmetry is dynamically broken in the vacuum. This is the well-known scenario of the dynamical chiral symmetry breaking in the NJL model.

Next, let us take the anomaly term into account with $g_{D}\neq 0$ in Eq. (1) while keeping the chiral limit. In this situation, even though $G_{S}$ is less than $1$, chiral symmetry can be broken dynamically in the vacuum due to the existence of the axial anomaly term. The previous work Kono:2021 demonstrated that such a situation is realized with a sufficiently large contribution from the anomaly term. We refer to such chiral symmetry breaking as the anomaly-driven breaking of chiral symmetry or shortly the anomaly-driven breaking in this paper.

The above patterns of the chiral symmetry breaking are related to the hadronic properties, such as the mass of the $\sigma$ meson Kono:2021 . In the literature, in order to study the relationship more quantitatively, the explicit chiral symmetry breaking by finite current quark mass is introduced and the $\sigma$ meson is assumed to be the chiral partner of the pion. The authors calculated the $\sigma$ meson mass with varying values of the dimensionless couplings $G_{S}$ and $G_{D}\equiv\Lambda_{3}g_{D}/(g_{S}^{\rm crit})^{2}$ so as to reproduce the $\eta^{\prime}$ mass. As a result, they found that the $\sigma$ mass is smaller than about $800~{}{\rm MeV}$ when the anomaly-driven breaking is realized, i.e., $G_{S}<1$, and larger than about $800~{}{\rm MeV}$ if the normal breaking is done, i.e., $G_{S}>1$.

According to the previous work Kono:2021 , the definition of the anomaly-driven breaking in the NJL model is that the chiral symmetry is dynamically broken even though the dimensionless coupling constant $G_{S}$ is less than $1$. However, the definition of this determination procedure relies on the model-specific parameter $G_{S}$, which makes its application to other models nontrivial. To discuss the anomaly-driven breaking in other models and systems, in the next section, we generalize the definition of the determination procedure for the anomaly-driven breaking.

In this subsection, we generalize the definition of the determination procedure for anomaly-driven breaking based on arguments in the NJL model. The anomaly-driven breaking of chiral symmetry was initially introduced in the chiral effective models using the model-specific coupling constants. Here, we present a key quantity that links the model-specific and the model-independent determination procedures for the anomaly-driven breaking. That is the sign of the curvature of the effective potential at the point with zero quark condensate. We use the latter as the definition of determination procedure for the anomaly-driven breaking in this paper.

We first consider the NJL model with no anomaly term in the vanishing quark mass limit $(m\to 0)$. In this case, as we have mentioned in the last section, chiral symmetry is dynamically broken in the vacuum when a dimensionless four-quark coupling $G_{S}$ is greater than $1$. Here, the evaluation of the second derivative of the effective potential (2) with respect to the quark condensate at the point $\braket{\bar{q}q}=0$ yields the following result:

$\displaystyle\left.\frac{\partial^{2}V_{\rm eff}}{\partial\braket{\bar{q}q}^{2}}\right|_{\braket{\bar{q}q}=0}=g_{S}^{\rm crit}-g_{S}=g_{S}^{\rm crit}\left(1-G_{S}\right),$

(3)

where we use $G_{S}=g_{S}/g_{S}^{\rm crit}$. From Eq. (3), we find that the parameter region where chiral symmetry is dynamically broken is linked to the negativity of the curvature of the effective potential at the origin. We refer to such chiral symmetry breaking as the normal breaking in contrast to the anomaly-driven one.

Next, let us turn on the anomaly term by $g_{D}\neq 0$ while keeping the chiral limit. Again, chiral symmetry can be dynamically broken and that leads to the nonzero quark condensate in the vacuum. The main difference compared to the case without the anomaly term is that chiral symmetry can be broken dynamically even though the dimensionless coupling $G_{S}$ is less than $1$. Calculating the second derivative of the effective potential (2) with the anomaly term, we obtain the same result to Eq. (3). Thus, the pattern of dynamical chiral symmetry breaking which is distinguished whether $G_{S}$ is greater than $1$ or not corresponds directly to the sign of the curvature (3).

In Fig. 1, we show the effective potentials as a function of the absolute value of quark condensate for four parameter sets in the chiral limit. In order to set the origin to zero, irrelevant constants are subtracted from the effective potentials. We can see that when the dimensionless coupling $G_{S}$ is greater than $1$, the curvature of effective potential at the origin is negative (red dotted line). On the other hand, when chiral symmetry is dynamically broken even though $G_{S}$ is less than $1$ (green solid line), the curvature is positive at the origin. For the remaining two parameter sets (magenta dash-dotted and blue dashed line), the effective potential has a minimum value at the origin and thus chiral symmetry is not broken in the vacuum. In this way, when the chiral symmetry is dynamically broken, the dimensionless coupling constant $G_{S}$ and the curvature of the effective potential at the origin are linked in the NJL model including the anomaly term for the chiral limit.

We also confirm that such a relationship between the dimensionless coupling $G_{S}$ and the curvature of the effective potential remains unchanged even if small current quark masses are introduced. The current quark mass dependence appears in the effective potential as a product with the linear term of the quark condensate for small quark masses. That term vanishes by performing the second derivative with respect to the quark condensate. Equation (3) thus remains the same in the presence of small current quark masses.

Based on the arguments in the NJL model, we use the curvature of the effective potential to determine whether the anomaly-driven breaking or not, instead of the model specific coupling constant. Since the effective potential is proportional to the vacuum energy density $\epsilon$ except for a constant in field theory, we take the inequality

$\displaystyle\left.\frac{\partial^{2}\epsilon}{\partial\braket{\bar{q}q}^{2}}\right|_{\braket{\bar{q}q}=0}>0,$

(4)

to be the definition of the determination procedure for the anomaly-driven breaking also for finite quark masses. We apply this definition to the instanton liquid model and discuss the feasibility of the anomaly-driven chiral symmetry breaking.

The QCD Euclid partition function is approximated by the superposition of instantons Schafer:1996 . The canonical partition function is given as

	$\displaystyle Z$	$\displaystyle=$	$\displaystyle\frac{1}{N_{+}!N_{-}!}\int\prod_{i=1}^{N_{+}+N_{-}}\left[d\Omega_{i}f(\rho_{i})\right]\exp(-S_{\rm int})$		(5)
			$\displaystyle\times\prod_{f=1}^{N_{f}}\mathrm{Det}({D}/+m_{f}).$

Here, $N_{+}$ and $N_{-}$ are the numbers of instantons and anti-instantons, respectively, $d\Omega_{i}=dU_{i}d\rho_{i}d^{4}z_{i}$ is the measure of the path integral over the collective coordinate space for the color orientation in the color ${\rm SU}(N_{c})$ group, size and position associated with the $i$-th instanton, and $f(\rho)$ is semiclassical instanton amplitude. The interactions between instanton-instanton and instanton-quark are included in $S_{\rm int}$ and $\mathrm{Det}({D}/+m_{f})$, respectively. $D_{\mu}$ and $m_{f}$ represent the covariant derivative for the quark fields and the current quark mass of flavor $f$, respectively.

The explicit expression of the semiclassical instanton amplitude calculated by ’t Hooft tHooft:1976_b reads

	$\displaystyle f(\rho)$	$\displaystyle=$	$\displaystyle C_{N_{c}}\left[\frac{8\pi^{2}}{g^{2}(\rho)}\right]^{2N_{c}}\exp\left[-\frac{8\pi^{2}}{g^{2}(\rho)}\right]\frac{1}{\rho^{5}}$		(6)
		$\displaystyle=$	$\displaystyle C_{N_{c}}\frac{1}{\rho^{5}}\beta_{1}(\rho)^{2N_{c}}$
			$\displaystyle\times\exp\left[-\beta_{2}(\rho)+\left(2N_{c}-\frac{b^{\prime}}{2b}\right)\frac{b^{\prime}}{2b}\frac{1}{\beta_{1}(\rho)}\ln\beta_{1}(\rho)\right]$
			$\displaystyle+{\mathcal{O}}(3\mathchar 45\relax{\rm loop}),$
	$\displaystyle C_{N_{c}}$	$\displaystyle=$	$\displaystyle\frac{0.466e^{-1.679N_{c}}}{(N_{c}-1)!(N_{c}-2)!},$		(7)

where the gauge coupling $g^{2}$ is given as a function of the instanton size $\rho$, the beta functions $\beta_{1},\beta_{2}$ include up to 2-loop order and the Gell-Mann–Low coefficients are given as follows:

	$\displaystyle\beta_{1}(\rho)$	$\displaystyle=$	$\displaystyle-b\ln(\rho\Lambda),\hskip 10.00002pt\beta_{2}(\rho)=\beta_{1}(\rho)+\frac{b^{\prime}}{2b}\ln\left[\frac{2}{b}\beta_{1}(\rho)\right],$
	$\displaystyle b$	$\displaystyle=$	$\displaystyle\frac{11}{3}N_{c}-\frac{2}{3}N_{f},\hskip 10.00002ptb^{\prime}=\frac{34}{3}N_{c}^{2}-\frac{13}{3}N_{c}N_{f}+\frac{N_{f}}{N_{c}}.$

Here, $N_{c}$ and $N_{f}$ denote the number of colors and of the dynamical quark flavors, respectively, and $\Lambda$ is the scale parameter which is introduced for calculating $\beta$-function by using the Pauli-Villars regularization scheme Schafer:1996 .

The instanton-instanton interaction $S_{\rm int}$ in Eq. (5) is expressed by the sum of all possible pairs of instantons

$\displaystyle S_{\rm int}=\sum_{l<m,l\neq m}^{N_{+}+N_{-}}S^{(2)}_{\rm int}(l,m),$

(10)

where $l,m$ refer to each instanton or anti-instanton in the ensemble. The two-body interaction between instantons (or anti-instantons) $S^{(2)}_{\rm int}(l,m)$ is defined as the difference between the action of two-instanton configuration and twice of the single instanton action $S_{0}(\rho)=8\pi/g^{2}(\rho)$:

$\displaystyle S^{(2)}_{\rm int}(l,m)=S[A_{\mu}(l,m)]-2S_{0}.$

(11)

Two-instanton configuration is no longer an exact solution to the classical Yang-Mills equations. To calculate the two-body interaction, one uses an Ansatz for the gauge configurations. This idea was developed by Schäfer and Shuryak Schafer:1996 ; Schafer:1998 . We use the streamline Ansatz Yung:1988 ; Verbaarschot:1991 and its form of $S_{\rm int}$ is given by Eq. (A5) in Ref. Schafer:1996 .

The instanton-quark interaction represented by the determinant in Eq. (5) is evaluated by factorizing it into two parts, a high and a low momentum parts as

$\displaystyle\mathrm{Det}({D}/+m_{f})$

$\displaystyle=$

$\displaystyle\mathrm{Det}_{\rm high}({D}/+m_{f})\mathrm{Det}_{\rm low}({D}/+m_{f}).~{}~{}~{}~{}$

(12)

The high momentum part is evaluated as the product of contributions of each instanton calculated by using the Gaussian approximation. The low momentum part is evaluated by using the quark zero-mode wave functions in the instanton and anti-instanton backgrounds Schafer:1996 . As a result, the instanton-quark interaction takes the following form:

			$\displaystyle\mathrm{Det}({D}/+m_{f})=\left(\prod_{i=1}^{N_{+}+N_{-}}1.34\rho_{i}\right)\mathrm{Det}_{{\rm I},{\rm\bar{I}}}(-iT+m_{f}\bm{1}),$		(13)
			$\displaystyle T=\begin{pmatrix}\bm{0}_{N_{+}\times N_{+}}&({\mathcal{T}})_{N_{+}\times N_{-}}\vspace{1em}\\ ({\mathcal{T}}^{\dagger})_{N_{-}\times N_{+}}&\bm{0}_{N_{-}\times N_{-}}\\ \end{pmatrix},$
			$\displaystyle({\mathcal{T}})_{IJ}=\int d^{4}x\,\psi^{*}_{0,I}(x;U,\rho,z)i\gamma_{\mu}D_{\mu}\psi_{0,J}(x;U,\rho,z).~{}~{}~{}~{}$		(14)

Here $T$ is called the overlap matrix with a size of $(N_{+}+N_{-})\times(N_{+}+N_{-})$. This matrix is spanned by the quark zero-mode wave functions $\psi_{0,I}(x;U,\rho,z)$ in the instanton (${\rm I}$) and anti-instanton $({\rm\bar{I}})$ backgrounds which are labeled by the collective coordinate $\{U,\rho,z\}$. The specific forms of the quark zero-mode wavefunctions are summarized in Appendix 2 of Ref. Schafer:1998 . The determinant operation $\mathrm{Det}_{{\rm I},{\rm\bar{I}}}$ is taken over the space represented in Eq. (13) and $m_{f}\bm{1}$ is understood as a diagonal matrix of $(N_{+}+N_{-})\times(N_{+}+N_{-})$. The explicit expression of ${\mathcal{T}}$ is given by Eq. (B2) in Ref. Schafer:1996 .

In this paper, we work on the ${\rm SU(3)}$ flavor symmetric limit, i.e., the number of quark flavors is set to $N_{f}=3$ and the current quark masses are equally set to $m=m_{f}~{}(f=1,2,3)$. The quark condensate $\braket{\bar{q}q}$ represents the one-flavor amount for our calculations.

The vacuum energy density is identified the free energy density as $\epsilon=F$ at zero temperature and we simply call it the free energy denoted as $F$ in what follows. This relationship is derived from the standard thermodynamics relation. The free energy is expressed as

$\displaystyle F=-\frac{1}{V}\ln Z,$

(15)

with the four-dimensional space-time volume $V=L^{4}$ and the partition function $Z$ of the system considered.

We explain the method to compute the value of the partition function $Z$. The method that we use is well-known as the thermodynamics integration method, and it has been applied to the IILM calculation in the previous work Schafer:1996 . In this method, one writes the effective action as

$\displaystyle S_{\rm eff}(\alpha)=S_{\rm eff}(0)+\alpha S_{1},$

(16)

with the partition function $Z(\alpha)=\int d\Omega\exp\left[-S_{\rm eff}(\alpha)\right]$. The desired partition function $Z$ is reproduced as $Z(\alpha=1)$. This form of the effective action (16) interpolates between a known solvable action $S_{\rm eff}(0)\equiv S_{\rm eff}^{0}$ and the full action $S_{\rm eff}(1)=S_{\rm eff}^{0}+S_{1}$. One obtains the partition function $Z(\alpha=1)$ straightforwardly as

$\displaystyle\ln\left[Z(\alpha=1)\right]=\ln\left[Z(\alpha=0)\right]-\int_{0}^{1}d\alpha\braket{0}{S_{1}}{0}_{\alpha},~{}~{}$

(17)

where the expectation value $\braket{0}{O}{0}_{\alpha}$ is defined by

$\displaystyle\braket{0}{O}{0}_{\alpha}\equiv\frac{1}{Z(\alpha)}\int d\Omega~{}O(\Omega)e^{-S_{\rm eff}(\alpha)},$

(18)

with configurations according to $p(\Omega_{i})\propto\exp\left[-S_{\rm eff}(\alpha)\right]$. Therefore, if we know the decomposition (16) and the values of $S_{\rm eff}^{0}$ and $Z(\alpha=0)$, we can compute the partition function $Z(\alpha=1)$ from Eq. (17).

For our computation of the partition function in the instanton liquid model, we use the following decomposition of the effective action $S_{\rm eff}$ as Ref. Schafer:1996 :

	$\displaystyle S_{\rm eff}(\alpha)$	$\displaystyle=$	$\displaystyle\sum_{i=1}^{N_{+}+N_{-}}\left\{\ln\left[f(\rho_{i})\right]+(1-\alpha)\nu\frac{\rho_{i}^{2}}{\bar{\rho}^{2}}\right\}$		(19)
			$\displaystyle+\alpha\left[S_{\rm int}-\sum_{f=1}^{N_{f}}\ln\mathrm{Det}({D}/+m_{f})\right],$

where $\nu=(b-4)/2$ is a coefficient with the Gell-Mann–Low coefficient $b$ given by Eq. (III.1), and $\bar{\rho}^{2}$ is the average size squared of instantons in the full ensembles including the instanton-instanton and instanton-quark interactions. The variational single instanton distribution is used as $Z(\alpha=0)\equiv Z_{0}$. Its form is given by

$\displaystyle Z_{0}$

$\displaystyle=$

$\displaystyle\frac{1}{N_{+}!N_{-}!}\left(V\mu_{0}\right)^{N_{+}+N_{-}},$

(20)

with

$\displaystyle\mu_{0}$

$\displaystyle=$

$\displaystyle\int_{0}^{\infty}d\rho~{}f(\rho)\exp\left(-\nu\frac{\rho^{2}}{\bar{\rho}^{2}}\right).$

(21)

In the instanton liquid model, the quark condensate is evaluated as the expectation value of the traced quark propagator at the same space-time coordinate as follows:

	$\displaystyle\braket{\bar{q}q}$	$\displaystyle=$	$\displaystyle\sum_{A,\alpha}\braket{q^{\dagger}(x)^{A}_{\alpha}q(x)^{A}_{\alpha}}$
		$\displaystyle=$	$\displaystyle-\lim_{y\to x}\sum_{A,\alpha}\braket{q(x)^{A}_{\alpha}q^{\dagger}(y)^{A}_{\alpha}}$
		$\displaystyle=$	$\displaystyle-\lim_{y\to x}\frac{1}{Z}\int D\Omega~{}\mathrm{Tr}\left[S(x,y;m)\right]e^{-S_{\rm int}}\mathrm{Det}({D}/+m).$

Here we write the measure of the path integral as $D\Omega$ in short that is given in the partition function (5), $A=1,\dots,N_{c}$ and $\alpha=1,\dots,4$ represent the color and the Dirac indices, respectively and $\mathrm{Tr}$ taken for the both indices. The quark propagator is approximated as a sum of contributions from the free and the zero-mode propagators by inverting the Dirac operator in the basis spanned by the quark zero-mode wave functions in instantons background as following Schafer:1998

			$\displaystyle S(x,y;m)\approx S_{0}(x,y)+S^{\rm ZM}(x,y;m),$
			$\displaystyle S_{0}(x,y)=\frac{i}{2\pi^{2}}\frac{\gamma\cdot(x-y)}{(x-y)^{4}},$
			$\displaystyle S^{\rm ZM}(x,y;m)=\sum_{I,J}\left[\psi_{0,I}(x)\left[(-iT+m)^{-1}\right]_{IJ}\psi^{\dagger}_{0,J}(y)\right],$

where the matrix $T$ is the overlap matrix given in Eq. (13). Here we omit to write the collective coordinates of instantons $\{U,\rho,z\}$ from the argument of the quark zero-mode wave functions. We obtain the quark condensate by averaging it over the configurations.

The simplest way to simulate the instanton liquid model described by the Euclid partition function (5) is to use the Monte Carlo method with a weight function $S_{\rm eff}$ given by

	$\displaystyle S_{\rm eff}=-\sum_{i=1}^{N_{+}+N_{-}}\ln[f(\rho_{i})]+S_{\rm int}-\sum_{f=1}^{N_{f}}\ln\mathrm{Det}({D}/+m_{f}).$
			(24)

To perform Monte Carlo simulations using the Markov Chain Monte Carlo (MCMC) method, it is crucial to understand the weight function $S_{\rm eff}$. This function represents the probability density $p(\Omega_{i})$ as $p(\Omega_{i})\propto\exp[-S_{\rm eff}(\Omega_{i})],$ where $\Omega_{i}$ denotes a configuration within the considered ensemble. Once the partition function is established, we derive the weight function, as illustrated in Eq. (24). Employing algorithms such as the Metropolis algorithm or Hybrid Monte Carlo (HMC) algorithm, we generate a series of configurations $\{\Omega\}=(\Omega_{1},\Omega_{2},\dots,\Omega_{N_{\rm conf}})$ that converge to the given probability density function because of the detailed balance condition of the algorithms. Here, $N_{\rm conf}$ represents the number of configurations.

The expectation value $\braket{O}$ of an operator $O$ that is expressed as a function of a configuration $\Omega_{i}$ is computed from these configurations using the formula:

$\displaystyle\braket{O}=\lim_{N_{\rm conf}\to\infty}\frac{1}{N_{\rm conf}}\sum_{i=1}^{N_{\rm conf}}O(\Omega_{i}).$

(25)

For more details on the Monte Carlo simulations using the MCMC, we referred to some textbooks and an introductory article Binder:2010 ; Randau:2014 ; Hanada:2018 .

In our simulations, each configuration $\Omega_{i}$ consists of $N_{+}$ instantons and $N_{-}$ anti-instantons labeled by their collective coordinate $\{U,\rho,z\}$ and the partition function for each instanton density is determined by generating 5000 configurations after 1000 initial sweeps with $N=N_{+}+N_{-}=16+16$ instantons and anti-instantons. The simulations with different instanton density $n=N/V$ are achieved by changing the simulation box size $V=L^{4}$ with the fixed number of instantons $N$. Whole simulations are performed under the periodic boundary condition for the coordinates of instantons and anti-instantons. All quantities in this calculation are nondimensionalized by $\Lambda$ which has been introduced through the $\beta$-functions (III.1). The value of the scale parameter $\Lambda$ is determined so that the free energy has a minimum value at the instanton density $n=1~{}{\rm fm^{-4}}$ Schafer:1996 . This density with the minimum free energy is the vacuum instanton density. We are interested in the quark condensate dependence of the free energy in the small quark mass regime, so we set the current quark masses $m$ to be as small as possible within a stable run of the simulations. The calculations below are performed with the current quark masses $m=0.1\Lambda,0.15\Lambda,0.2\Lambda$.

In Table 1, we show the bulk parameters, such as the current quark mass $m$, the instanton density in the vacuum $n$, the scale parameter $\Lambda$ and the instanton average size $\bar{\rho}$ in our simulations. These values are consistent with the previous work Schafer:1996 . By fixing the unit such that $n=1~{}{\rm fm}^{-4}$ at the vacuum, the values of the scale parameter are determined as $\Lambda=373,360,351~{}{\rm MeV}$ for each current quark mass. Using these values, the average instanton sizes in our calculations are evaluated as $\bar{\rho}=0.35$ to $0.37~{}{\rm fm}$ compared to $\bar{\rho}=0.42~{}{\rm fm}$ in Ref. Schafer:1996 .

To calculate the free energy by Eq. (15), we need the value of the partition function. The partition function is obtained through Eq. (17), so we first calculate $Z_{0}$ with the variational single instanton distribution (21). The single instanton distribution can be evaluated from initial 1000 sweeps with full interaction $\alpha=1$. We calculate the average size squared of instantons $\bar{\rho}^{2}$ appearing in Eq. (21) using the initial sweeps. The $\rho$-integral over an infinite interval $[0,\infty]$ appearing in Eq. (21) is practically performed over a finite interval $[0,\Lambda^{-1}]$. The remaining task for the calculation of the partition function is to perform the integral by summing the integrands. This integral is done at 10 different coupling values $\alpha$.

In the computation of the quark condensate, only the quark zero-mode propagator $S^{\rm ZM}(x,y;m)$ in Eq. (LABEL:eq:S_full) is evaluated to subtract the infinite contribution initiated by the free propagator $S_{0}(x,x)$ at the same space-time coordinate. In the actual calculations, for each instanton density, the trace of the quark zero-mode propagator is averaged over the 5000 configurations, and also averaged over 10,000 different space-time coordinates to reduce the effect of incomplete equilibration of the configurations.

In Fig. 2, we show our numerical results of the free energies (in the unit of ${\rm MeV/fm^{3}}$) versus the instanton density (in the unit of ${\rm fm^{-4}}$) for three values of the current quark mass. Our results are consistent with the previous work Schafer:1996 . The free energy is monotonically decreasing towards a minimum value and then rapidly increasing at higher instanton density. This shows that attractive interaction is dominant in lower instanton density regions, while repulsive interaction becomes important at higher instanton density.

In Table 2, we summarize our numerical results of the free energy and the quark condensate at the vacuum for three different current quark masses $m=37,54,70~{}{\rm MeV}$. By using the values of the scale parameter in Table 1, the free energies are evaluated as $F=-149,-187,-228~{}{\rm MeV/fm^{3}}$. For reference, we perform further calculations for the current quark masses with values close to those in the previous work, such as $m=96~{}{\rm MeV}$ Schafer:1996 . That shows a good agreement with the previous work, and we conclude that our simulations work well.

In Fig. 3, we show the instanton density dependence of cubic root of the quark condensate $(-\braket{\bar{q}q})^{1/3}$ in the unit of ${\rm MeV}$. We find that the absolute values of the quark condensate increase monotonically as the instanton density increases. We also find that the value of the quark condensates at the vacuum is insensitive to the value of the current quark masses. This means that the contribution from the explicit breaking of chiral symmetry due to the current quark mass is not so large for the value of the quark condensate.

In Table 2, we show the values of the quark condensate at the vacuum instanton density. Our results almost reproduce the empirical values obtained by various previous works GellMann:1968 ; Reinders:1985 ; Dosch:1998 ; Harnett:2021 ; Giusti:2001 ; McNeile:2013 ; Borsanyi:2013 ; Durr:2014 ; Cossu:2016 ; Davies:2019 ; FLAG:2021 ; Gasser:1985 ; Jamin:2001 ; Boyle:2016 ; Kneur:2020 . We obtain the values of the quark condensate as $|\braket{\bar{q}q}|^{1/3}=188,196,200~{}{\rm MeV}$ at the vacuum with the current quark masses $m=37,54,70~{}{\rm MeV}$, respectively. These values give a good estimate of the quark condensates although they are slightly smaller in the absolute value.

Combining the results of the free energy and the quark condensate, we obtain the quark condensate dependence of the free energy as shown in Fig. 4. The free energies monotonically decrease towards the minimum values as the quark condensate increases in magnitude. Once the free energies have the minimum value, they start to increase as expected from the instanton density dependence of them (also see Fig. 2). From this, we observe that the free energy has its minimum value at the point with the finite value of the quark condensate. This shows that chiral symmetry is dynamically broken in the vacuum of the IILM.

The behavior of the free energy near the origin appears to be decreasing in a downward convex trend. This trend is crucial for the sign of the curvature of the free energy at the origin. In other words, it provides one of the hints for discriminating patterns of chiral symmetry breaking through our definition of the determination procedure for the anomaly-driven breaking discussed in Sect. II.2. So we study the quark condensate dependence of the free energy more quantitatively.

We aim to evaluate the curvature of the free energy at the origin from our simulation results of the IILM. For that, we perform polynomial fits for our three data sets with different current quark masses. Each data set consists of $N_{\rm data}$ pairs of $\braket{\bar{q}q}$ and $F$, i.e., $(|\braket{\bar{q}q}|_{i},F_{i}),i=1,\dots,N_{\rm data}$ arranged in ascending order, where $N_{\rm data}$ is the number of data of the free energy and the quark condensate and it is equal to the number of the grid points of the instanton density, $N_{\rm data}=71$ in our calculations.

We use a polynomial function given by

$\displaystyle F(\braket{\bar{q}q})=\sum_{j=0}^{K}C_{j}\braket{\bar{q}q}^{j},$

(26)

for the fitting model in this analysis. We perform the fits for four orders $K=1,\dots,4$. For each order $K$, we optimize the parameters $\{C_{j}\}~{}(j=0,\dots,K)$ so that they minimize the reduced chi-square including errors of both axes as given in Ref. Orear:1984 by

$\displaystyle\chi^{2}_{\rm d.o.f.}=\frac{1}{N_{\rm d.o.f.}}\sum_{i=1}^{M}\frac{(y_{i}-f(x_{i}))^{2}}{\sigma_{y_{i}}^{2}+\sigma_{x_{i}}^{2}\left[f^{\prime}(x_{i})\right]^{2}},$

(27)

where $(x_{i},y_{i})$ with their errors $(\sigma_{x_{i}},\sigma_{y_{i}})$ correspond to our data set $(|\braket{\bar{q}q}|_{i},F_{i})$ with their errors, $f(x)$ and $f^{\prime}(x)$ represent the fit model (26) and its first derivative, $N_{\rm d.o.f.}$ is the number of degrees of freedom which is defined by $N_{\rm d.o.f.}\equiv M-(K+1)$, and $M$ is the number of data used in the fit. We use the data from $i=1$ to $i=M$ rather all data because we want to know the behavior of the free energy around the origin.

We determine an upper limit of fit range, $M$, as follows. We calculate the reduced chi-square (27) as we increase the number of data that we use and obtain the reduced chi-square as a function of the number of data. We then find the number of data for which the reduced chi-square has a local minimum value. We finally use this number of data as the value of $M$ for fitting. Here we skip a trivial local minimum obtained with a small number of data. The specific values of $M$ used in each fit, the corresponding quark condensate values and the reduced chi-square are summarized in Table 3.

In Table 4, we show the fit results of $C_{0}$. We find that the value of the coefficient $C_{0}$ is close to zero. For different fit orders and quark masses, the values of $C_{0}$ range from $-1.52$ to $-0.47$ in the unit of ${\rm MeV/fm^{3}}$. These values are sufficiently small compared to the typical value of the free energy $F$ and we conclude that the coefficient $C_{0}$ is consistent with zero.

In Table 5, we show the fit results of $C_{1}$. We find that the value of $C_{1}$ is insensitive to the fitting order $K$. This implies that the value of $C_{1}$ is well determined around the origin. We find also that $C_{1}$ has no clear quark mass dependence. The averaged values of $C_{1}$ over the different orders are $389~{}{\rm MeV},320~{}{\rm MeV}$ and $360~{}{\rm MeV}$ for $m=37,54$ and $70~{}{\rm MeV}$.

In Table 6, we show the fit results of $C_{2}$. We find that the coefficients $C_{2}$ appear to be positive. For each quark mass, the averaged values of $C_{2}$ over the different fit orders are $C_{2}=(2.94,0.98,1.17)\times 10^{-5}~{}{\rm MeV^{-2}}$ for $m=37,54$ and $70~{}{\rm MeV}$. These values of $C_{2}$ are located at the positive region for all orders and current quark masses. This suggests that the curvature of the free energy at the origin is positive for the IILM.

In Fig. 5, we show the current quark mass dependence of the coefficient $C_{2}$. The value of $C_{2}$ is not very sensitive to the value of the current quark mass. For the fit with the smallest quark mass, the values of $C_{2}$ are evaluated as slightly larger than the results of other two quark masses. These results suggest that the value of $C_{2}$ and the curvature of the free energy at the origin is positive for wide current quark mass region and thus we conclude that the anomaly-driven breaking of chiral symmetry is realized in the IILM by our definition (4)

Here we do not show the values of the coefficients $C_{3}$ and $C_{4}$ because these coefficients are introduced to confirm the stability of fits of $C_{0},C_{1}$ and $C_{2}$ with and without these higher-order terms. Furthermore, since we use the part of full data near the origin for fits, the values of $C_{3}$ and $C_{4}$ are not determined well and irrelevant to our discussion.

Interestingly we find the opposite sign of $C_{2}$ in the quenched calculations where the quark determinant part in the partition function (5) is set as $\prod_{f}\mathrm{Det}({D}/+m_{f})=1$ in generating the configurations. The details of the quenched calculation are shown in Appendix A. As we show in Table 7, the quenched calculations suggest the negative values of coefficient $C_{2}$. By our definition of the anomaly-driven breaking (4), this concludes that the normal breaking is realized in the quenched IILM. The opposite signs of $C_{2}$ obtained by the full and quenched calculations mean that the different scenarios of chiral symmetry breaking are possible depending on the presence of the quark determinant part in the IILM partition function. Since the quark determinant contributes to the instanton-quark interaction in the ensemble, this implies that the quarks play a crucial role in the anomaly-driven breaking for the IILM.