This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

thanks: huangll22@mails.jlu.edu.cnthanks: yuanyuanw23@jlu.edu.cnthanks: hxzhang18@163.comthanks: ishidah@pu-toyama.ac.jpthanks: mamiya@ucas.ac.cnthanks: synya@jlu.edu.cnthanks: akio@yukawa.kyoto-u.ac.jp

Impact of local CP-odd domain in hot QCD on axionic domain-wall interpretation
for NANOGrav 15-year Data

Linlin Huang Center for Theoretical Physics and College of Physics, Jilin University, Changchun, 130012, China    Yuanyuan Wang Center for Theoretical Physics and College of Physics, Jilin University, Changchun, 130012, China    He-Xu Zhang Center for Theoretical Physics and College of Physics, Jilin University, Changchun, 130012, China    Hiroyuki Ishida Center for Liberal Arts and Sciences, Toyama Prefectural University, Toyama 939-0398, Japan    Mamiya Kawaguchi School of Nuclear Science and Technology, University of Chinese Academy of Sciences, Beijing 100049, China    Shinya Matsuzaki Center for Theoretical Physics and College of Physics, Jilin University, Changchun, 130012, China    Akio Tomiya Department of Information Technology, International Professional University of Technology in Osaka, 3-3-1 Umeda, Kita-Ku, Osaka, 530-0001, Japan
Abstract

We argue that the axionic domain-wall with a QCD bias may be incompatible with the NANOGrav 15-year data on a stochastic gravitational wave (GW) background, when the domain wall network collapses in the hot-QCD induced local CP-odd domain. This is due to the drastic suppression of the QCD bias set by the QCD topological susceptibility in the presence of the CP-odd domain with nonzero θ\theta parameter of order one which the QCD sphaleron could generate. We quantify the effect on the GW signals by working on a low-energy effective model of Nambu-Jona-Lasinio type in the mean field approximation. We find that only at θ=π\theta=\pi, the QCD bias tends to get significantly large enough due to the criticality of the thermal CP restoration, which would, however, give too big signal strengths to be consistent with the NANOGrav 15-year data and would also be subject to the strength of the phase transition at the criticality.

I Introduction

The observation of a stochastic GW background has recently reported from the NANOGrav pulsar timing array collaboration in 15 years of data NANOGrav:2023gor ; NANOGrav:2023hfp . Possible origins of the detected nano-Hz peak frequency in a view of particle physics have been investigated also by the NANOGrav collaboration NANOGrav:2023hvm . Other recent pulsar timing array (PTA) data, such as those from the European PTA (EPTA) EPTA:2023sfo ; EPTA:2023akd ; EPTA:2023fyk , Parkes PTA (PPTA) Reardon:2023gzh ; Reardon:2023zen , and Chinese PTA (CPTA) Xu:2023wog have also supported the presence of consistent nano-Hz stochastic GWs. Thus this nano-Hz GW evidence might provide us with a hint on the new aspect of the thermal history of the universe in terms of beyond the standard model of particle physics.

Among the new-physics candidate interpretations, the axionlike particle (ALP)-domain wall annihilation triggered by a QCD-induced bias has been considered as an attractive model Li:2024psa ; Ellis:2023oxs ; Lozanov:2023rcd ; Gelmini:2023kvo ; Kitajima:2023cek ; Geller:2023shn ; Bai:2023cqj ; Blasi:2023sej , which can naturally be realized in the QCD-phase transition epoch at the temperature T𝒪T\sim{\cal O}(100) MeV consistent with the produced peak frequency of nano Hz NANOGrav:2023hvm . The QCD-induced bias is supplied by the topological susceptibility χtop\chi_{\rm top}, and its TT-dependence and the value at T=0T=0 have already been measured in a lattice simulation at the physical point for quark masses with the continuum limit is properly taken Borsanyi:2016ksw . This is how the ALP-domain wall prediction gets definite and unambiguous except the domain-wall network formulation and annihilation analysis, even though the system that the ALP acts in is nonperturbative QCD.

However, the thermal history of the QCD phase transition epoch may not be so simple: a local CP-odd domain may be created in hot QCD plasma due to the presence of the QCD sphaleron Manton:1983nd ; Klinkhamer:1984di , so that the QCD vacuum characterized by the strong CP phase θ\theta and its fluctuation (in the spatial-homogeneous direction) gets significantly sizable Kharzeev:2007tn ; Kharzeev:2007jp ; Fukushima:2008xe within the QCD time scale McLerran:1990de ; Moore:1997im ; Moore:1999fs ; Bodeker:1999gx  #1#1#1 The QCD sphaleron transition rate is not suppressed by the thermal effect in contrast to the QCD instanton’s one Moore:1997im ; Moore:1999fs ; Bodeker:1999gx . The topological charge fluctuation (within the QCD time scale =𝒪(110)={\cal O}(1\mathchar 45\relax 10) fm), Δqtop(t)=t𝑑t𝑑x3Qtop(t,x)\Delta q_{\rm top}(t)=\int^{t}dt^{\prime}d\vec{x}^{3}Q_{\rm top}(t^{\prime},\vec{x}), will be non-vanishing. This implies that the corresponding source θ(x)\theta(x) in the generating functional ZQCDZ_{\rm QCD} needs to be (at least) time-dependent and fluctuate: δZQCD/δΔθ(t)qtop(t)0\delta Z_{\rm QCD}/\delta\Delta\theta(t)\sim q_{\rm top}(t)\neq 0. The time and/or thermal average of Δθ(t)\Delta\theta(t) thus acts as what we call the theta parameter in the QCD thermal plasma, which namely means θΔθ(t)\theta\equiv\Delta\theta(t) in there. See also, e.g., the literature Andrianov:2012hq ; Andrianov:2012dj . Besides, the time fluctuation tθ(t)\partial_{t}\theta(t), to be referred to as the chiral chemical potential Kharzeev:2007tn ; Kharzeev:2007jp ; Fukushima:2008xe ; Andrianov:2012hq ; Andrianov:2012dj , will be significant as well when the non-conservation law of the U(1)U(1) axial symmetry is addressed. See also Summary and Discussions. . Though θ\theta should be tiny enough (<1010<10^{-10}) at present, nonzero θ\theta contribution to the QCD bias χtop\chi_{\rm top} may be non-negligible when the ALP domain wall starts to collapse in the QCD phase transition epoch. This is irrespective to whether or not the ALP acts as the QCD axion which relaxes the θ\theta to zero at present.

In this paper, we argue that the ALP-domain wall with the QCD bias becomes incompatible with the NANOGrav 15-year data on a stochastic GW background, when the domain wall network collapses in the hot-QCD induced local-CP odd domain. This happens due to the drastic suppression of the QCD bias set by the QCD topological susceptibility in the presence of the CP-odd domain with nonzero θ\theta parameter of order one which the QCD sphaleron could generate.

This paper is structured as follows. We first make a generic argument. It is based only on the anomalous Ward-Takahashi identities in QCD and the mixing structure between the scalar quark condensate (q¯q\langle\bar{q}q\rangle) and pseudoscalar quark condensate (q¯iγ5q\langle\bar{q}i\gamma_{5}q\rangle) bilinear operators via the U(1)U(1) axial transformation with nonzero θ\theta. χtop\chi_{\rm top} is shown to generically get small when θ\theta takes the value of order one, because of the dramatic suppression of the scalar quark condensate at any temperature.

We next implement this suppression effect into the produced GW signals by working on a low-energy effective model of Nambu-Jona-Lasinio (NJL) type in the mean field approximation (MFA) and the random phase approximation (RPA). In accord with the generic argument, χtop\chi_{\rm top} as well as the scalar quark condensate get highly suppressed when θ=𝒪(1)\theta={\cal O}(1), instead, the pseudoscalar condensate as the CP/axial partner develops.

We also find that only at θ=π\theta=\pi, the QCD bias tends to get significantly large enough due to the criticality of the thermal CP restoration of the second order, which is signaled when the pseudoscalar condensate reaches zero. This, however, gives too big signal strengths to be consistent with the NANOGrav 15-year data. Going beyond the MFA or other types of effective models could yield a different strength of the phase transition at the criticality. In summary of the present paper, we also briefly address the outlook along this possibility and the prospective impact on the QCD bias for the ALP-domain wall annihilation in light of nano Hz GW signals.

II General argument

We begin with presenting a general consequence on the suppression of χtop\chi_{\rm top} when θ=𝒪(1)\theta={\cal O}(1). First of all, we note that the anomalous U(1)U(1) axial transformation (qeiγ5θ4qqq\to e^{-i\gamma_{5}\frac{\theta}{4}}q\equiv q^{\prime}) can make the θ\theta dependence in QCD present only in the quark mass term. In two-flavor QCD, for simplicity, the quark mass term then takes the form

q=u,d[mcosθ2(q¯q)+msinθ2(q¯iγ5q)],\displaystyle\sum_{q=u,d}\left[m\cos\frac{\theta}{2}(\bar{q}^{\prime}q^{\prime})+m\sin\frac{\theta}{2}(\bar{q}^{\prime}i\gamma_{5}q^{\prime})\right]\,, (1)

where the primed quark bilinear fields are related to the original ones via the orthogonal rotation,

(q¯qq¯iγ5q)=(cosθ2sinθ2sinθ2cosθ2)(q¯qq¯iγ5q).\displaystyle\begin{pmatrix}\bar{q}q\\ \bar{q}i\gamma_{5}q\end{pmatrix}=\begin{pmatrix}\cos\frac{\theta}{2}&\sin\frac{\theta}{2}\vspace{2mm}\\ -\sin\frac{\theta}{2}&\cos\frac{\theta}{2}\end{pmatrix}\begin{pmatrix}\bar{q}^{\prime}q^{\prime}\\ \bar{q}^{\prime}i\gamma_{5}q^{\prime}\end{pmatrix}\,. (2)

The scalar condensate q¯q\langle\bar{q}q\rangle thus mixes with the pseudoscalar one q¯iγ5q\langle\bar{q}i\gamma_{5}q\rangle by nonzero θ\theta.

χtop\chi_{\rm top} is given as the functional derivative of the generating functional of QCD twice with respect to θ\theta evaluated at θ0\theta\neq 0. Taking into account the quark mass term in Eq.(1) we thus get

iχtop(T,θ)\displaystyle i\chi_{\rm top}(T,\theta) =Td4xQtop(x)Qtop(0)θ,\displaystyle=\int_{T}d^{4}x\langle Q_{\rm top}(x)Q_{\rm top}(0)\rangle_{\theta}\,,
Qtop(x)\displaystyle Q_{\rm top}(x) =gs232π2GμνG~μν,\displaystyle=\frac{g_{s}^{2}}{32\pi^{2}}G_{\mu\nu}\widetilde{G}^{\mu\nu}\,, (3)

where QtopQ_{\rm top} denotes the topological charge with the QCD gauge coupling gsg_{s} and the (dual) field strength GμνG_{\mu\nu} (G~μνϵμνρσ2Gρσ\tilde{G}_{\mu\nu}\equiv\frac{\epsilon_{\mu\nu\rho\sigma}}{2}G^{\rho\sigma}); Td4x01/T𝑑τd3x\int_{T}d^{4}x\equiv\int^{1/T}_{0}d\tau\int d^{3}x with the imaginary time τ=ix0\tau=ix_{0}; the subscript “θ\theta” attached on the vacuum (thermodynamic ground) states stands for the implicit θ\theta parameter dependence. Following the procedure in the literature Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr χtop\chi_{\rm top} in two-flavor QCD is expressed in terms of the original-basis fields as

χtop(T,θ)=14[mq=u,dq¯qθ,T+im2χη(T,θ)],\displaystyle\chi_{\rm top}(T,\theta)=-\frac{1}{4}\left[m\sum_{q=u,d}\langle\bar{q}q\rangle_{\theta,T}+im^{2}\chi_{\eta}(T,\theta)\right]\,, (4)

where we have taken the isospin symmetric mass for up (uu) and down (dd) quarks, and

χη(T,θ)=Td4xq=u,d(q¯(0)iγ5q(0))(q¯(x)iγ5q(x))θ,\displaystyle\chi_{\eta}(T,\theta)=\int_{T}d^{4}x\sum_{q=u,d}\langle(\bar{q}(0)i\gamma_{5}q(0))(\bar{q}(x)i\gamma_{5}q(x))\rangle_{\theta}\,, (5)

which is Tq=u,dq¯iγ5qθ,T\sim\frac{\partial}{\partial T}\sum_{q=u,d}\langle\bar{q}i\gamma_{5}q\rangle_{\theta,T}. Since the quark mass is perturbatively small enough and q¯q\langle\bar{q}q\rangle develops a nonzero value at the normal QCD vacuum with θ=0\theta=0, the chiral perturbation conventionally works well for χtop\chi_{\rm top}, so that the q¯q\langle\bar{q}q\rangle term will dominate over the χη\chi_{\eta} term, to saturate the measured χtop\chi_{\rm top} value (75.6MeV)4\simeq(75.6\,{\rm MeV})^{4}, i.e.,

|χtop(T,θ=0)|14mq=u,d|q¯qθ=0,T|(75.6MeV)4.\displaystyle|\chi_{\rm top}(T,\theta=0)|\approx\frac{1}{4}m\sum_{q=u,d}|\langle\bar{q}q\rangle_{\theta=0,T}|\simeq(75.6\,{\rm MeV})^{4}\,. (6)

This is sort of the well-known formula called the Leutwyler-Smiluga formula Leutwyler:1992yt , and this feature has also been observed in chiral effective model approaches at any TT Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr . When θ0\theta\neq 0, however, the value of q¯q\langle\bar{q}q\rangle will be shared with the pseudoscalar condensate q¯iγ5q\langle\bar{q}i\gamma_{5}q\rangle through Eq.(2), so that q¯q\langle\bar{q}q\rangle term in χtop\chi_{\rm top} of Eq.(4) gets small, to be as small as or smaller than the χη\chi_{\eta} term, as θ\theta gets sizable. Thus, we expect the dramatic suppression of χtop\chi_{\rm top} at a sizable θ\theta #2#2#2 Even when the strange quark contributions are incorporated in χtop\chi_{\rm top}, the form of Eq.(4) is intact, as has been discussed in the literature Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr and also reviewed in Appendix A, hence χtop\chi_{\rm top} can be evaluated only via the lightest quark condensate and χη\chi_{\eta}, in which the strange-quark loop contributions are implicitly incorporated. Thus the lightest quark condensate term (q=u,dq¯q\sum_{q=u,d}\langle\bar{q}q\rangle) in |χtop(T,θ=𝒪(1))||\chi_{\rm top}(T,\theta={\cal O}(1))| still persists being suppressed due to the sizable CP violation, while the χη\chi_{\eta} term keeps sizable with the large CP violation, where the CP-violating strange quark contribution almost decouples simply because of its heaviness (see also Summary and Discussions). Therefore, the inequality in Eq.(7) holds even in the case of three-flavor QCD. :

|χtop(T,θ=𝒪(1))|14m2|χη(T,θ)||χtop(T,θ=0)|.\displaystyle|\chi_{\rm top}(T,\theta={\cal O}(1))|\approx\frac{1}{4}m^{2}|\chi_{\eta}(T,\theta)|\ll|\chi_{\rm top}(T,\theta=0)|\,. (7)

This trend will be seen for any TT and χtop\chi_{\rm top} would simply get smaller and smaller as TT increases unless a sharp phase transition in χη\chi_{\eta} shows up at higher TT. Thus, the magnitude of the QCD-induced bias for the ALP domain wall annihilation, controlled by χtop\chi_{\rm top}, is expected to become dramatically small in the local CP-odd domain created in hot QCD.

Below we will explicitize this claim based on an NJL with nonzero θ\theta.

III NJL evaluation of χtop\chi_{\rm top} and ALP domain wall energy density with nonzero θ\theta

Since lattice data on χtop(T,θ)\chi_{\rm top}(T,\theta) with θ0\theta\neq 0 have not yet been available #3#3#3For more on the current status and future prospects, see also Summary and Discussions., we employ a low-energy chiral effective model, an NJL model, which matches with the underlying QCD via the consistent anomalous Ward-Takahashi identities associated with the chiral SU(2)SU(2) and U(1)U(1) axial symmetry breaking. The investigation of QCD with nonzero θ\theta has so far been carried out based on several chiral effective models  Dashen:1970et ; Witten:1980sp ; Pisarski:1996ne ; Creutz:2003xu ; Mizher:2008hf ; Boer:2008ct ; Creutz:2009kx ; Boomsma:2009eh ; Sakai:2011gs ; Chatterjee:2011yz ; Sasaki:2011cj ; Sasaki:2013ewa ; Aoki:2014moa ; Mameda:2014cxa ; Verbaarschot:2014upa ; Bai:2023cqj and also the recently developed ‘t Hooft-anomaly matching method Gaiotto:2017tne ; Gaiotto:2017yup extended from the original idea tHooft:1979rat ; Frishman:1980dq , so as to clarify the nature of the thermal chiral and strong CP phase transitions when θ=π\theta=\pi. Nevertheless, the TT and θ\theta dependence on χtop\chi_{\rm top} has never been addressed except Ref. Sasaki:2013ewa . In the reference χtop\chi_{\rm top} was computed based on the same NJL model as what we will work on below, however, consistency with the anomalous chiral Ward-Takahashi identities was not manifest, which is to be refined in the present study. We leave all the technical details in Appendices A and B and shall only give the results relevant to the discussion of the nano Hz GW signals.

The dependence of TT and θ\theta on χtop\chi_{\rm top} is plotted in Fig. 1. The present model yields |χtop1/4(T=0,θ=0)|77.5|\chi_{\rm top}^{1/4}(T=0,\theta=0)|\simeq 77.5 MeV, which is in good agreement with the lattice estimate, χtop1/4(T=0,θ=0)75.6\chi^{1/4}_{\rm top}(T=0,\theta=0)\simeq 75.6 MeV Borsanyi:2016ksw . As seen from the left panel of the figure, χtop\chi_{\rm top} prominently gets smaller when θ0.5π\theta\gtrsim 0.5\pi, which is due to the suppression of q¯q\langle\bar{q}q\rangle (as shown in Appendix B), instead q¯iγ5q\langle\bar{q}i\gamma_{5}q\rangle gets greater than q¯q\langle\bar{q}q\rangle, as was expected from the generic argument in Eq.(7). We also note from Appendix B that the chiral phase transition goes like crossover for any θ\theta, whereas the CP phase transition is of the second order type at θ=π\theta=\pi (see Fig. 1), which is also manifested as a spike structure of χtop\chi_{\rm top}. In particular, the CP symmetry is restored at the criticality (T221T\simeq 221 MeV) in accordance with the literature Boer:2008ct ; Boomsma:2009eh ; Sakai:2011gs .

Refer to caption
Figure 1: Plot of χtop\chi_{\rm top} (in magnitude) computed from the present two-flavor NJL model with nonzero θ\theta.

We assume that an ALP has already been present before the QCD phase transition epoch as in the literature Li:2024psa ; Ellis:2023oxs ; Lozanov:2023rcd ; Gelmini:2023kvo ; Kitajima:2023cek ; Geller:2023shn ; Bai:2023cqj ; Blasi:2023sej and developed the potential having the shift symmetry, aa+2πfaa\to a+2\pi f_{a} with nDWn_{\rm DW} being integer:

V0(a)=ma2fa2nDW2(1cos(nDWafa)),\displaystyle V_{0}(a)=\frac{m_{a}^{2}f_{a}^{2}}{n_{\rm DW}^{2}}\left(1-\cos\left(n_{\rm DW}\frac{a}{f_{a}}\right)\right)\,, (8)

where mam_{a} and faf_{a} are the ALP mass and decay constant, respectively. Then a domain wall profile exists as a soliton solution, which sweeps between the adjacent vacua, where nDWn_{\rm DW} in V0V_{0} above corresponds to the domain wall number. As the universe cools down to the QCD scale, assuming existence of the ALP coupling to gluon fields, the ALP develops another potential via the U(1)U(1) axial anomaly,

VQCD(a)=χtop(T,θ)cos(nafa+θ),\displaystyle V_{\rm QCD}(a)=-\chi_{\rm top}(T,\theta)\cos\left(n\frac{a}{f_{a}}+\theta\right)\,, (9)

which explicitly breaks the original shift symmetry, where nn is some factor related to and highly depending on the Peccei-Quinn charges of quarks. The original ALP vacuum-potential energy ma2fa2\sim m_{a}^{2}f_{a}^{2} is assumed to be |χtop|\gtrsim|\chi_{\rm top}|. Thus the domain wall configuration, which is supported by the original shift symmetry, becomes unstable in the total ALP potential Vtotal=V0+VQCDV_{\rm total}=V_{0}+V_{\rm QCD}:

Vtotal(a)=ma2fa2nDW2(1cos(nDWafa))+χtop(cosθcos(nafa+θ)),\displaystyle V_{\rm total}(a)=\frac{m_{a}^{2}f_{a}^{2}}{n_{\rm DW}^{2}}\left(1-\cos\left(n_{\rm DW}\frac{a}{f_{a}}\right)\right)+\chi_{\rm top}\left(\cos\theta-\cos\left(n\frac{a}{f_{a}}+\theta\right)\right)\,, (10)

with the normalization Vtotal(a=0)=0V_{\rm total}(a=0)=0. Meanwhile, the domain wall network (with nDW2n_{\rm DW}\geq 2) starts to collapse and release the latent heat into the universe. Here, VQCD(a)V_{\rm QCD}(a) plays the role of what is called the bias, ΔV\Delta V. The released energy density can be evaluated by the vacuum energy at the original vacua a=nDWπfaa=n_{\rm DW}\pi f_{a}. For instance, when nDW=2n_{\rm DW}=2, the original vacuum at a/fa=πa/f_{a}=\pi gets the energy shift by |Vtotal(a/fa=π)|=|V0VQCD|a/fa=π|V_{\rm total}(a/f_{a}=\pi)|=|V_{0}-V_{\rm QCD}|_{a/f_{a}=\pi}:

|V0VQCD|a/fa=π=|χtop(T,θ)(cos(nπ+θ)cosθ)||χtop(T,θ)|ΔV(T,θ).\displaystyle|V_{0}-V_{\rm QCD}|_{a/f_{a}=\pi}=|\chi_{\rm top}(T,\theta)(\cos(n\pi+\theta)-\cos\theta)|\sim|\chi_{\rm top}(T,\theta)|\equiv\Delta V(T,\theta)\,. (11)

Thus one may maximally have the released latent heat ρDW(T)\rho_{\rm DW}(T)

ρDWΔV(T,θ)=|χtop(T,θ)|.\displaystyle\rho_{\rm DW}\sim\Delta V(T,\theta)=|\chi_{\rm top}(T,\theta)|\,. (12)

Similar discussions have been made in the literature Li:2024psa ; Ellis:2023oxs ; Lozanov:2023rcd ; Gelmini:2023kvo ; Kitajima:2023cek ; Geller:2023shn ; Bai:2023cqj ; Blasi:2023sej . A more precise estimate based on the numerical simulations of the domain wall network suggests ρDW0.5ΔV\rho_{\rm DW}\simeq 0.5\Delta V Saikawa:2017hiv ; Hiramatsu:2010yz ; Hiramatsu:2012sc ; Hiramatsu:2013qaa .

The QCD bias |χtop||\chi_{\rm top}| generically does not significantly affect the ALP potential when the universe is hotter than the QCD scale of 𝒪(100){\cal O}(100) MeV, as seen from Fig. 1 and also from the generic formula in Eq.(4). This is essentially due to the correlation of the effective restoration of the chiral SU(2)SU(2) (via q¯q\langle\bar{q}q\rangle) and/or U(1)U(1) axial symmetry (χη\chi_{\eta}) at higher temperatures (See Eq.(4)). We assume that the ALP-domain wall annihilation takes place at T=TT=T_{*} by the QCD bias and fully provides the source of the GW #4#4#4 The ALP domain wall network will completely be decayed never to be left before the expected domain wall-dominated epoch arises. The critical temperature for the domain wall domination is estimated by using the present NJL model and evaluating the condition ρDW/ρrad>1\rho_{\rm DW}/\rho_{\rm rad}>1, to give T=Tdom19T=T_{\rm dom}\sim 19 MeV, which is indeed much lower than the QCD scale. . The produced GW power spectrum at the peak frequency is then assessed via the signal strength, free from the domain wall string tension σDWmafa2\sigma_{\rm DW}\sim m_{a}f_{a}^{2}, by the following signal strength (see, e.g., Blasi:2023sej ):

α(T,θ)\displaystyle\alpha_{*}(T_{*},\theta) =ρDW(T)ρrad(T)0.5ΔV(T,θ)ρrad(T)\displaystyle=\frac{\rho_{\rm DW}(T_{*})}{\rho_{\rm rad}(T_{*})}\simeq\frac{0.5\Delta V(T_{*},\theta)}{\rho_{\rm rad}(T_{*})}
0.15×(|χtop(T,θ)|1/4100MeV)4(T100MeV)4(g(T)10)1,\displaystyle\simeq 0.15\times\left(\frac{|\chi_{\rm top}(T_{*},\theta)|^{1/4}}{100\,{\rm MeV}}\right)^{4}\left(\frac{T_{*}}{100\,{\rm MeV}}\right)^{-4}\left(\frac{g_{*}(T_{*})}{10}\right)^{-1}\,, (13)

where we have used ρDW0.5ΔV\rho_{\rm DW}\simeq 0.5\Delta V Saikawa:2017hiv ; Hiramatsu:2010yz ; Hiramatsu:2012sc ; Hiramatsu:2013qaa and ρrad(T)=(π2/30)g(T)T4\rho_{\rm rad}(T)=(\pi^{2}/30)g_{*}(T)T^{4}. This α\alpha_{*} has been constrained by the NANOGrav 15yr data set as a function of TT_{*}. See Fig. 2, showing that the cases with θ/π0.5\theta/\pi\gtrsim 0.5 are incompatible with the interpretation of the NANOGrav 15-year data at around T=100T_{*}=100 MeV. This is essentially due to the suppression of the scalar quark condensate as observed in χtop\chi_{\rm top} as shown in Fig. 1. At θ=π\theta=\pi onset the criticality (at T221T\simeq 221 MeV), the signal gets dramatically enhanced by the singular spike structure reflecting the second order phase transition of the CP symmetry, which is, however, to be too big to interpret the data. This feature is insensitive to the deconfinement phase transition of QCD, which has not yet been incorporated in the figure.

Refer to caption
Figure 2: The contour plot for the signal strength α(θ)\alpha_{*}(\theta) in Eq.(13) versus the ALP domain wall annihilation temperature TT_{*} with χtop\chi_{\rm top} in Fig. 1 encoded. The 2σ\sigma contour from Ref. NANOGrav:2023hvm has also been displayed in blue, and the regimes inside the 2σ\sigma contour are thus allowed. α\alpha_{*} in Eq.(13) depends on the effective degrees of freedom g(T)g_{*}(T) around the QCD phase transition epoch, which we have adjusted g(T)g_{*}(T) available in Ref. ParticleDataGroup:2022pth , in such a way that the deconfinement transition happens at the same time as the chiral phase transition in the present NJL model takes place (at T=221T=221 MeV). The wiggles seen around a lower TT regime in the model prediction curves have thus been arisen from the thresholds encoded in g(T)g_{*}(T).

IV Summary and discussions

In summary, the ALP domain wall with a QCD bias may be incompatible with the NANOGrav 15-year data on a stochastic GW background, when the domain wall network collapses in the hot-QCD induced local-CP odd domain which could allow to have a sizable QCD θ\theta parameter. This is due to the drastic suppression of the QCD bias set by the QCD topological susceptibility, χtop\chi_{\rm top}, in the presence of the CP-odd domain with θ=𝒪(1)\theta={\cal O}(1) (see Eq.(7) and Fig. 1). An explicit model analysis of χtop\chi_{\rm top} with a large θ\theta based on the two-flavor NJL model implies that only at θ=π\theta=\pi, the QCD bias tends to get significantly large enough due to the criticality of the thermal CP restoration. However, this turned out to give too big signal strengths to be consistent with the NANOGrav 15-year data.

In closing, we give several comments on the issues to be pursed in the future.

  • The presently employed MFA of NJL model does not precisely reproduce the chiral crossover at θ=0\theta=0 as has been observed in the lattice simulations. The latter predicts the pseudocritical temperature Tpc155T_{\rm pc}\simeq 155 MeV Aoki:2009sc ; Borsanyi:2011bn ; Ding:2015ona ; Bazavov:2018mes ; Ding:2020rtq , while the present model yields Tpc220T_{\rm pc}\simeq 220 MeV. This may imply a simple shift of the TT-distribution of χtop\chi_{\rm top} in Fig. 1 and also α(θ)\alpha_{*}(\theta) toward lower TT with lower criticality, say, at T160T\sim 160 MeV. Even in that case, however, the trend of the dramatic suppression of χtop\chi_{\rm top} with θ=𝒪(1)\theta={\cal O}(1) should keep manifest for below the shifted criticality (160\sim 160 MeV) and the spike signal associated with the strong CP restoration at θ=π\theta=\pi will be left at the criticality.

  • In Appendix B we have also taken into account a sort of the QCD deconfinement phase transition by extending the NJL model into the Polyakov-loop NJL (PNJL) model Fukushima:2003fw (for a recent review, see, e.g., Fukushima:2017csk ). It turns out that the high suppression of χtop\chi_{\rm top} with θ=𝒪(1)\theta={\cal O}(1) is still manifest and the presence of a sharp spike in χtop\chi_{\rm top} at θ=π\theta=\pi still persists, in accord with the earlier work Sakai:2011gs . This implies the insensitivity of the deconfinement-confinement transition for the main conclusion addressed in the main text.

  • The local CP-odd domain would generate not only a large θ\theta, but also the fluctuation of θ\theta in the temporal direction, tθ(t)\partial_{t}\theta(t), which is identified as the so-called chiral chemical potential (often denoted as μ5\mu_{5}Kharzeev:2007tn ; Kharzeev:2007jp ; Fukushima:2008xe . The μ5\mu_{5} contribution to the thermal chiral phase transition as well as χtop\chi_{\rm top} has been discussed in Ref. Ruggieri:2020qtq based on a two-flavor NJL model with θ=0\theta=0. From the reference we can see that a part of χtop\chi_{\rm top}, which only includes the scalar condensate term q¯q\langle\bar{q}q\rangle as in Eq.(6) (without the χη\chi_{\eta} contribution in Eq.(4)), roughly gets enhanced by about a factor of (3/2)(3/2) at around T=𝒪(100)T={\cal O}(100) MeV. With θ=𝒪(1)\theta={\cal O}(1), the scalar condensate q¯q\langle\bar{q}q\rangle generically gets smaller due to the CP violation as in the generic argument (Sec. II), which still holds even in the presence of μ5\mu_{5} because μ5tθ(t)\mu_{5}\sim\partial_{t}\theta(t) does not modify the mixture structure between (q¯q)(\bar{q}q) and (q¯iγ5q)(\bar{q}i\gamma_{5}q) as in Eq.(2). Furthermore, the presence of μ5\mu_{5} as well as nonzero θ\theta does not change the form of all the anomalous chiral Ward-identities as clarified in Appendix A. Hence the high suppression of χtop\chi_{\rm top} would still be seen even with μ5\mu_{5}, so the conclusion as presently claimed would be intact. More precise discussions including the size of the spike at θ=π\theta=\pi would be worth pursing elsewhere.

  • The criticality associated with the thermal CP restoration at θ=π\theta=\pi would be crucial to precisely check if the ALP domain wall annihilation can still be viable to account for the NANOGrav 15-year data. The related spike structure, following the order of the phase transition, the first order or the second order, would be subject to effective model approaches for QCD Dashen:1970et ; Witten:1980sp ; Pisarski:1996ne ; Creutz:2003xu ; Mizher:2008hf ; Boer:2008ct ; Creutz:2009kx ; Boomsma:2009eh ; Sakai:2011gs ; Chatterjee:2011yz ; Sasaki:2011cj ; Sasaki:2013ewa ; Aoki:2014moa ; Gaiotto:2017tne ; Verbaarschot:2014upa ; Gaiotto:2017yup ; Bai:2023cqj . The presently employed NJL-type model with the MFA tends to predict the second order Boer:2008ct ; Boomsma:2009eh ; Sakai:2011gs , while the linear-sigma model type with the MFA Pisarski:1996ne ; Creutz:2003xu ; Mizher:2008hf ; Boer:2008ct ; Creutz:2009kx ; Boomsma:2009eh ; Bai:2023cqj and the ‘t Hooft anomaly matching argument Gaiotto:2017tne ; Gaiotto:2017yup supports the first order phase transition. Even going beyond the MFA and/or the RPA, including subleading order contributions in 1/Nc1/N_{c} expansion, the presently claimed incompatibility would be crucial to deduce a more definite conclusion of the compatibility of the ALP domain wall interpretation with NANOGrav 15-year data. The functional renormalization group analysis makes it possible to clarify the case, which deserves to another publication. At any rate, the presently claimed incompatibility would still pin down one benchmark point in the full parameter space of the QCD-biased ALP domain wall collapse.

  • As seen from Figs. 3 and 4 in Appendix B, at θ=π\theta=\pi the thermal CP restoration point coincides with the pseudocritical point (inflection point) for the chiral crossover. This is due to the present MFA, which involves the two intrinsic features at θ=π\theta=\pi: i) the chiral symmetry is not spontaneously broken via the scalar quark condensate q¯q\langle\bar{q}q\rangle, as long as the CP symmetry is broken; ii) no U(1)U(1) axial violating loop corrections are generated at this approximation level. Those lead to no thermal development of q¯q\langle\bar{q}q\rangle until the CP symmetry restores, hence q¯q\langle\bar{q}q\rangle gets kicked down instantaneously, i.e., undergoes the pseudocritical point at the same timing as the CP restoration. We have clarified those points in Appendix B (around Eq.(61)). Going beyond the MFA analysis, such as the functional renormalization group, would potentially generate a gap between two critical points for the (effective) chiral and CP symmetry restorations, which should come from the U(1)U(1) axial anomaly effects in loops. The further investigation along this line could also make a complementary benchmark compared to the prediction from the ‘t Hooft-anomaly matching in pure Yang-Mills theories Chen:2020syd , which suggests the CP restoration temperature is higher than or equal to the deconfinement phase transition one.

  • When strange quark contributions are incorporated in the analysis, the thermal CP restoration at θ=π\theta=\pi would be more involved. First of all, one notices that the CP phase contribution carried by the strange quark is highly suppressed by a factor of mu,d/msm_{u,d}/m_{s}, as reviewed in Appendix A (around Eqs.(23) - (25)). This observation comes from the robust flavor singlet nature of the θ\theta dependence in QCD that requires the strange quark field to carry the θ\theta phase with the form eimu,d2msθ1e^{i\frac{m_{u,d}}{2m_{s}}\theta}\simeq 1, i.e., almost free from θ\theta. Thus, as far as the CP violating effects are concerned, the three-flavor QCD is essentially decomposed into the 2+12+1 (the lightest two quarks and the strange quark) structure and the CP order parameter is almost controlled by the lightest two-flavor sector, i.e., q=u,dq¯iγ5q\sum_{q=u,d}\langle\bar{q}i\gamma_{5}q\rangle.

    This is the characteristic flavor violation and irrespective to the presence of the QCD topological charge fluctuation, which is flavor universal, or equivalently in the NJL framework, the ‘t Hooft-Kobayashi-Maskawa determinant term Kobayashi:1970ji ; Kobayashi:1971qz ; tHooft:1976rip ; tHooft:1976snw (as in Appendix B). In the three-flavor NJL with the MFA, thus the chiral crossover for s¯s\langle\bar{s}s\rangle would still persist even at θ=π\theta=\pi, simply because of presence of the finite strange quark mass, while q¯q\langle\bar{q}q\rangle would be subject to the approximation for QCD, or the (P)NJL. In the MFA of the three-flavor NJL, q¯q\langle\bar{q}q\rangle would be trapped to a constant value until the CP phase transition takes place at T=TcT=T_{c}, and then starts to drops when T>TcT>T_{c} (with discontinuity in the TT-derivative at T=TcT=T_{c}), as in the two-flavor case (Fig. 3 in Appendix B). As in the case with θ=0\theta=0, s¯s\langle\bar{s}s\rangle would drop down with TT more slowly than q¯q\langle\bar{q}q\rangle simply due to the larger strange quark mass. Thus, in the framework of the MFA of the NJL, the CP phase transition in the three-flavor case would be essentially identical to the one in the lightest two-flavor case, hence it would still be of the second order keeping the critical spike structure of χtop\chi_{\rm top}, hence α\alpha_{*} does as well, as seen from Fig. 2.

    However, as has been discussed in the recent literature Fejos:2016hbp ; Fejos:2023lvw ; Fejos:2021yod based on the functional renormalization group analysis on the three-flavor linear sigma model with θ=0\theta=0, the ‘t Hooft-Kobayashi-Maskawa determinant term beyond the MFA could be nonperturbatively enhanced at around T=TcT=T_{c} even in the case of θ=π\theta=\pi. This implies that the criticality of the CP phase transition as well as the chiral phase transition might significantly be altered and the entry of the strange quark contributions beyond the MFA might be nontrivial, though the carried CP phase is intrinsically highly suppressed as aforementioned. Thus, the detailed analysis in the three-flavor NJL, both within and beyond the MFA, is noteworthy to pursue in another publication.

Finally, we comment on impacts of lattice QCD calculations to this observation. Most of current lattice QCD calculations at the physical point have been done using the Monte-Carlo method. For nonzero theta, simulations are suffered from infamous sign problem. Namely, we cannot calculate expectation values with |θ|0|\theta|\gg 0. To avoid the problem in the realistic setup, we employ imaginary theta and analytic continue to the real theta to get expectation values for physical observables DelDebbio:2002xa ; DElia:2003zne ; DelDebbio:2006sbu ; Giusti:2007tu ; Izubuchi:2007rmy ; Vicari:2008jw ; DElia:2013uaf ; Hirasawa:2024vns . These calculations for χtop\chi_{\rm top} have been done away from the physical point. An interesting methods for nonzero theta is suggested Kitano:2021jho , but it is still hard to get physical results. Digital and analog quantum simulations and calculations with nonzero theta using tensor networks have been performed for toy models but not for the realistic setup like four dimensional QCD at the physical point Byrnes:2002gj ; Funcke:2019zna ; Kuramashi:2019cgs ; Chakraborty:2020uhf ; Honda:2021aum ; Honda:2022hyu .

Acknowledgments

We are grateful to Taishi Katsuragawa and Wen Yin for useful comments. This work was supported in part by the National Science Foundation of China (NSFC) under Grant No.11747308, 11975108, 12047569, and the Seeds Funding of Jilin University (S.M.), and Toyama First Bank, Ltd (H.I.). The work by M.K. was supported by the Fundamental Research Funds for the Central Universities and partially by the National Natural Science Foundation of China (NSFC) Grant No. 12235016, and the Strategic Priority Research Program of Chinese Academy of Sciences under Grant No. XDB34030000. The work of A.T. was partially supported by JSPS KAKENHI Grant Numbers 20K14479, 22K03539, 22H05112, and 22H05111, and MEXT as “Program for Promoting Researches on the Supercomputer Fugaku” (Simulation for basic science: approaching the new quantum era; Grant Number JPMXP1020230411, and Search for physics beyond the standard model using large-scale lattice QCD simulation and development of AI technology toward next-generation lattice QCD; Grant Number JPMXP1020230409).

Appendix A Generic properties of anomalous chiral Ward-Takahashi identities (ACWTIs) in QCD

A.1 ACWTIs with external gauge fields and θ\theta

We start with the QCD Lagrangian with NfN_{f} quarks including the external gauge fields,

QCD=q¯LiγμDμqL+q¯RiγμDμqRq¯L𝐦fqRq¯R𝐦fqL+12a=18tr[(GμνaTca)2]+θg264π2ϵμνρσGμνaGρσa,\displaystyle{\cal L}_{\rm QCD}=\bar{q}_{L}i\gamma^{\mu}D_{\mu}q_{L}+\bar{q}_{R}i\gamma^{\mu}D_{\mu}q_{R}-\bar{q}_{L}{\bf m}_{f}q_{R}-\bar{q}_{R}{\bf m}_{f}q_{L}+\frac{1}{2}\sum_{a=1}^{8}{\rm tr}\left[(G_{\mu\nu}^{a}T^{a}_{c})^{2}\right]+\theta\frac{g^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{a}G_{\rho\sigma}^{a}\,, (14)

where qL(R)q_{L(R)} denotes the left-(right-) handed quark fields belong to the fundamental representation of SU(Nf)SU(N_{f}); 𝐦f{\bf m}_{f} is the current quark mass for ff-quark taking the diagonal form like 𝐦f=diag(mu,md,){\bf m}_{f}={\rm diag}(m_{u},m_{d},\cdots); GμνaG_{\mu\nu}^{a} (a=0,1,,8)(a=0,1,\cdots,8) are the field strengths of the gluon fields GμaG_{\mu}^{a}; TcaT^{a}_{c} stand for the generators of the QCD color group SU(3)cSU(3)_{c}; gg denotes the QCD gauge coupling constant; θ=θ(x)\theta=\theta(x) plays the role of the source for the topological operator ig2ϵμνρσGμνaGρσaig^{2}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{a}G_{\rho\sigma}^{a}; The external gauge fields μA{\cal L}_{\mu}^{A} and μA{\cal R}_{\mu}^{A} (A=0,,Nf21A=0,\cdots,N^{2}_{f}-1) are introduced by gauging the global chiral U(Nf)LU(N_{f})_{L} and U(Nf)RU(N_{f})_{R} symmetry, which are embedded into the covariant derivatives as

DμqL\displaystyle D_{\mu}q_{L} =(μigGμaTcaiμATfA)qL,\displaystyle=\left(\partial_{\mu}-igG_{\mu}^{a}T^{a}_{c}-i{\cal L}_{\mu}^{A}T_{f}^{A}\right)q_{L}\,,
DμqR\displaystyle D_{\mu}q_{R} =(μigGμaTcaiμATfA)qR,\displaystyle=\left(\partial_{\mu}-igG_{\mu}^{a}T^{a}_{c}-i{\cal R}_{\mu}^{A}T_{f}^{A}\right)q_{R}\,, (15)

with TfAT_{f}^{A} being the generators of U(Nf)U(N_{f}) in the flavor space.

The vacuum energy of QCD is given by

VQCD(θ)=ilnZQCD,\displaystyle V_{{\rm QCD}}(\theta)=-i\ln Z_{{\rm QCD}}\,, (16)

where ZQCDZ_{{\rm QCD}} represents the generating functional of QCD in Minkowski spacetime,

ZQCD=[dqdq¯][dG]exp[id4xQCD].\displaystyle Z_{\rm QCD}=\int[dqd\bar{q}][dG]\exp\Biggl{[}i\int d^{4}x{\cal L}_{\rm QCD}\Biggl{]}\,. (17)

We define the topological susceptibility χtop\chi_{\rm top} including the θ\theta dependence:

χtop(θ)\displaystyle\chi_{\rm top}(\theta) =\displaystyle= d4xδ2VQCD(θ)δθ(x)δθ(0)id4xδVQCDδθ(x)δVQCDδθ(0)\displaystyle-\int d^{4}x\frac{\delta^{2}V_{{\rm QCD}}(\theta)}{\delta\theta(x)\delta\theta(0)}-i\int d^{4}x\frac{\delta V_{\rm QCD}}{\delta\theta(x)}\frac{\delta V_{\rm QCD}}{\delta\theta(0)} (18)
=\displaystyle= id4x0|T(g264π2ϵμνρσGμνaGρσa)(x)(g264π2ϵμνρσGμνbGρσb)(0)|0θ,\displaystyle-i\int d^{4}x\left\langle 0\Big{|}T\left(\frac{g^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{a}G_{\rho\sigma}^{a}\right)(x)\left(\frac{g^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{b}G_{\rho\sigma}^{b}\right)(0)\Big{|}0\right\rangle_{\theta}\,,

where we have introduced the subscript θ\theta to explicitize the θ\theta dependence on the vacuum.

The external fields, such as the electromagnetic field AμA_{\mu}, the baryon chemical potential μB\mu_{\rm B}, and the chiral chemical potential μ5\mu_{5}, are embedded in the external gauge fields as #5#5#5 The isospin chemical potential can also be incorporated, when Nf=2N_{f}=2, into the covariant derivative as an additional external gauge field term proportional to (σ3/2)(\sigma^{3}/2).

μATfA\displaystyle{\cal L}_{\mu}^{A}T_{f}^{A} =\displaystyle= [eQemAμ+μB3𝟏Nf×Nf][μ5𝟏Nf×Nf],\displaystyle\left[eQ_{\rm em}A_{\mu}+\frac{\mu_{\rm B}}{3}{\bm{1}}_{N_{f}\times N_{f}}\right]-\left[\mu_{5}{\bm{1}}_{N_{f}\times N_{f}}\right],
μaTfa\displaystyle{\cal R}_{\mu}^{a}T_{f}^{a} =\displaystyle= [eQemAμ+μB3𝟏Nf×Nf]+[μ5𝟏Nf×Nf],\displaystyle\left[eQ_{\rm em}A_{\mu}+\frac{\mu_{\rm B}}{3}{\bm{1}}_{N_{f}\times N_{f}}\right]+\left[\mu_{5}{\bm{1}}_{N_{f}\times N_{f}}\right]\,, (19)

where ee represents the electromagnetic coupling constant, and QemQ_{\rm em} denotes the electric charge matrix for quarks Qem=diag(Qemu,Qemd,)Q_{\rm em}={\rm diag}(Q_{\rm em}^{u},Q_{\rm em}^{d},\cdots). The chemical potentials have been introduced as constant fields. In QCD with the external gauge fields, the U(1)AU(1)_{A} symmetry is explicitly broken by the current quark mass term and the gluonic quantum anomaly, and also external electromagnetic field. This is reflected in the anomalous conservation law of the U(1)U(1) axial current for each quark in the NfN_{f}-plet quark field qq, labeled as qfq_{f},

μjA(f)μ\displaystyle\partial_{\mu}j^{(f)\mu}_{A} =\displaystyle= 2iq¯fmfγ5qf+g232π2ϵμνρσGμνaGρσa+Nce2[Qemf]232π2ϵμνρσFμνFρσ,\displaystyle 2i\bar{q}^{f}m_{f}\gamma_{5}q^{f}+\frac{g^{2}}{32\pi^{2}}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{a}G^{a}_{\rho\sigma}+N_{c}\frac{e^{2}[Q_{\rm em}^{f}]^{2}}{32\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\,, (20)

with the U(1)U(1) axial current jA(f)μ=q¯fγ5γμqfj_{A}^{(f)\mu}=\bar{q}^{f}\gamma_{5}\gamma^{\mu}q^{f}. Note that the constant chemical potentials do not contribute to the anomalous conservation law.

Under the U(1)AU(1)_{A} rotation with the rotation angle αA\alpha_{A}, the quark fields transform as

qLf\displaystyle q_{L}^{f} \displaystyle\to exp(iαAf/2)qLf,\displaystyle\exp(-i\alpha_{A}^{f}/2)q_{L}^{\prime f}\,,
qRf\displaystyle q_{R}^{f} \displaystyle\to exp(+iαAf/2)qRf.\displaystyle\exp(+i\alpha_{A}^{f}/2)q_{R}^{\prime f}\,. (21)

Then the QCD generating functional gets shifted as

[dqdq¯][dG]exp[id4x(q¯LiγμDμqL+q¯RiγμDμqRf(q¯LfmfeiαAfqRf+q¯RfmfeiαAfqLf)+12a=18tr[(GμνaTca)2]\displaystyle\int[dq^{\prime}d\bar{q}^{\prime}][dG]\exp\Biggl{[}i\int d^{4}x\Biggl{(}\bar{q}_{L}^{\prime}i\gamma^{\mu}D_{\mu}q_{L}^{\prime}+\bar{q}_{R}^{\prime}i\gamma^{\mu}D_{\mu}q_{R}^{\prime}-\sum_{f}\left(\bar{q}_{L}^{\prime f}m_{f}e^{i\alpha_{A}^{f}}q_{R}^{\prime f}+\bar{q}^{\prime f}_{R}m_{f}e^{-i\alpha_{A}^{f}}q^{\prime f}_{L}\right)+\frac{1}{2}\sum_{a=1}^{8}{\rm tr}\left[(G_{\mu\nu}^{a}T^{a}_{c})^{2}\right]
+(θfαAf)g264π2ϵμνρσGμνaGρσa(fαAfQf2)Nce264π2ϵμνρσFμνFρσ)].\displaystyle+\left(\theta-\sum_{f}\alpha_{A}^{f}\right)\frac{g^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}G_{\mu\nu}^{a}G_{\rho\sigma}^{a}-\left(\sum_{f}\alpha_{A}^{f}Q_{f}^{2}\right)N_{c}\frac{e^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\Biggl{)}\Biggl{]}\,. (22)

To rotate the QCD θ\theta term away, we take the following phase choice reflecting the flavor-singlet nature of the QCD vacuum Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr ,

αAf=m¯mfθ,\displaystyle\alpha_{A}^{f}=\frac{\bar{m}}{m_{f}}\theta\,, (23)

with

m¯=(f1mf)1.\displaystyle\bar{m}=\left(\sum_{f}\frac{1}{m_{f}}\right)^{-1}\,. (24)

Then the QCD generating functional goes like

ZQCD\displaystyle Z_{\rm QCD} =\displaystyle= [dqdq¯][dG]exp[id4x(q¯LiγμDμqL+q¯RiγμDμqR+12a=18tr[(GμνaTca)2]\displaystyle\int[dq^{\prime}d\bar{q}^{\prime}][dG]\exp\Biggl{[}i\int d^{4}x\Biggl{(}\bar{q}_{L}^{\prime}i\gamma^{\mu}D_{\mu}q_{L}^{\prime}+\bar{q}_{R}^{\prime}i\gamma^{\mu}D_{\mu}q_{R}^{\prime}+\frac{1}{2}\sum_{a=1}^{8}{\rm tr}\left[(G_{\mu\nu}^{a}T^{a}_{c})^{2}\right] (25)
f(q¯Lfmfexp(im¯mfθ)qRf+q¯Rfmfexp(im¯mfθ)qLf)\displaystyle-\sum_{f}\left(\bar{q}^{\prime f}_{L}m_{f}\exp\left({i\frac{\bar{m}}{m_{f}}\theta}\right)q^{\prime f}_{R}+\bar{q}^{\prime f}_{R}m_{f}\exp\left({-i\frac{\bar{m}}{m_{f}}\theta}\right)q^{\prime f}_{L}\right)
(fm¯mfθQf2)Nce264π2ϵμνρσFμνFρσ)].\displaystyle-\left(\sum_{f}\frac{\bar{m}}{m_{f}}\theta Q_{f}^{2}\right)N_{c}\frac{e^{2}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\Biggl{)}\Biggl{]}.

From this generating functional, the topological susceptibility is given as

χtop(θ)\displaystyle\chi_{\rm top}(\theta) =\displaystyle= χtop(0)(θ)+χtop(EM)(θ),\displaystyle\chi_{\rm top}^{(0)}(\theta)+\chi_{\rm top}^{\rm(EM)}(\theta)\,, (26)

where

χtop(0)(θ)\displaystyle\chi_{\rm top}^{(0)}(\theta) =\displaystyle= m¯2[f1mf0|(q¯fqfcosθ2+q¯fiγ5qfsinθ2)|0θ\displaystyle-\bar{m}^{2}\Biggl{[}\sum_{f}\frac{1}{m_{f}}\left\langle 0\left|\left(\bar{q}^{\prime f}q^{\prime f}\cos\frac{\theta}{2}+\bar{q}^{\prime f}i\gamma_{5}q^{\prime f}\sin\frac{\theta}{2}\right)\right|0\right\rangle_{\theta}
+id4x0|Tf(q¯fiγ5qfcosθ2q¯fqfsinθ2)(x)f(q¯fiγ5qfcosθ2q¯fqfsinθ2)(0)|0θ],\displaystyle+i\int d^{4}x\left\langle 0\left|T\sum_{f}\left(\bar{q}^{\prime f}i\gamma_{5}q^{\prime f}\cos\frac{\theta}{2}-\bar{q}^{\prime f}q^{\prime f}\sin\frac{\theta}{2}\right)(x)\sum_{f^{\prime}}\left(\bar{q}^{\prime f^{\prime}}i\gamma_{5}q^{\prime f^{\prime}}\cos\frac{\theta}{2}-\bar{q}^{\prime f^{\prime}}q^{\prime f^{\prime}}\sin\frac{\theta}{2}\right)(0)\right|0\right\rangle_{\theta}\Biggl{]},
χtop(EM)(θ)\displaystyle\chi_{\rm top}^{\rm(EM)}(\theta) =\displaystyle= i(fm¯Qf2mf)2d4x0|T(e2Nc64π2ϵμνρσFμνFρσ)(x)(e2Nc64π2ϵμνρσFμνFρσ)(0)|0θ\displaystyle-i\left(\sum_{f}\frac{\bar{m}Q_{f}^{2}}{m_{f}}\right)^{2}\int d^{4}x\left\langle 0\Big{|}T\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(x)\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(0)\Big{|}0\right\rangle_{\theta} (27)
i(fm¯2Qf2mf)d4x0|T(e2Nc64π2ϵμνρσFμνFρσ)(x)f(q¯fiγ5qfcosθ2q¯fqfsinθ2)(0)|0θ\displaystyle-i\left(\sum_{f}\frac{\bar{m}^{2}Q_{f}^{2}}{m_{f}}\right)\int d^{4}x\left\langle 0\Big{|}T\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(x)\sum_{f}\left(\bar{q}^{\prime f}i\gamma_{5}q^{\prime f}\cos\frac{\theta}{2}-\bar{q}^{\prime f}q^{\prime f}\sin\frac{\theta}{2}\right)(0)\Big{|}0\right\rangle_{\theta}
i(fm¯2Qf2mf)d4x0|T(q¯fiγ5qfcosθ2q¯fqfsinθ2)(x)(e2Nc64π2ϵμνρσFμνFρσ)(0)|0θ.\displaystyle-i\left(\sum_{f}\frac{\bar{m}^{2}Q_{f}^{2}}{m_{f}}\right)\int d^{4}x\left\langle 0\Big{|}T\left(\bar{q}^{\prime f}i\gamma_{5}q^{\prime f}\cos\frac{\theta}{2}-\bar{q}^{\prime f}q^{\prime f}\sin\frac{\theta}{2}\right)(x)\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(0)\Big{|}0\right\rangle_{\theta}\,.

These susceptibilities are written in terms of the primed quark field qq^{\prime}. They can be rewritten in terms of the original quark field qq using the following connection,

q=exp(iθγ5/4)q.\displaystyle q=\exp(i\theta\gamma_{5}/4)q^{\prime}. (28)

Then χtop(0)(θ)\chi_{\rm top}^{(0)}(\theta) and χtop(EM)(θ)\chi_{\rm top}^{\rm(EM)}(\theta) go like

χtop(0)(θ)\displaystyle\chi_{\rm top}^{(0)}(\theta) =m¯2[f1mf0|q¯fqf|0+id4x0|T(fq¯fiγ5qf)(x)(fq¯fiγ5qf)(0)|0θ],\displaystyle=-\bar{m}^{2}\left[\sum_{f}\frac{1}{m_{f}}\langle 0|\bar{q}^{f}q^{f}|0\rangle+i\int d^{4}x\left\langle 0\Big{|}T\left(\sum_{f}\bar{q}^{f}i\gamma_{5}q^{f}\right)(x)\left(\sum_{f^{\prime}}\bar{q}^{f^{\prime}}i\gamma_{5}q^{f^{\prime}}\right)(0)\Big{|}0\right\rangle_{\theta}\right],
χtop(EM)(θ)\displaystyle\chi_{\rm top}^{\rm(EM)}(\theta) =i(fm¯Qf2mf)2d4x0|T(e2Nc64π2ϵμνρσFμνFρσ)(x)(e2Nc64π2ϵμνρσFμνFρσ)(0)|0θ\displaystyle=-i\left(\sum_{f}\frac{\bar{m}Q_{f}^{2}}{m_{f}}\right)^{2}\int d^{4}x\left\langle 0\Big{|}T\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(x)\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(0)\Big{|}0\right\rangle_{\theta}
i(fm¯2Qf2mf)d4x0|T(e2Nc64π2ϵμνρσFμνFρσ)(x)(fq¯fiγ5qf)(0)|0θ\displaystyle-i\left(\sum_{f}\frac{\bar{m}^{2}Q_{f}^{2}}{m_{f}}\right)\int d^{4}x\left\langle 0\Big{|}T\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(x)\left(\sum_{f}\bar{q}^{f}i\gamma_{5}q^{f}\right)(0)\Big{|}0\right\rangle_{\theta}
i(fm¯2Qf2mf)d4x0|T(fq¯fiγ5qf)(x)(e2Nc64π2ϵμνρσFμνFρσ)(0)|0θ.\displaystyle-i\left(\sum_{f}\frac{\bar{m}^{2}Q_{f}^{2}}{m_{f}}\right)\int d^{4}x\left\langle 0\Big{|}T\left(\sum_{f}\bar{q}^{f}i\gamma_{5}q^{f}\right)(x)\left(\frac{e^{2}N_{c}}{64\pi^{2}}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}F_{\rho\sigma}\right)(0)\Big{|}0\right\rangle_{\theta}\,. (29)

The χtop\chi_{\rm top} term in Eq.(26) is precisely the generalization of the one in Eq.(4) with Eq.(5), which has also been addressed in the literature Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr without external gauge fields. Thus the presence of the electromagnetic axial anomaly yields additional topological susceptibility, χtop(EM)\chi_{\rm top}^{({\rm EM})}.

Next, we evaluate the ACWTIs for SU(Nf)L×SU(Nf)RSU(N_{f})_{L}\times SU(N_{f})_{R} transformation in QCD. Under the chiral SU(Nf)SU(N_{f}) transformation, the quark field qq transforms as

q(x)\displaystyle q(x) \displaystyle\to q(x)+iαA(x)TAγ5q(x),A=1,Nf21.\displaystyle q(x)+i\alpha^{A}(x)T^{A}\gamma_{5}q(x)\,,\qquad\,A=1,\cdots N_{f}^{2}-1\,. (30)

Then, the chiral SU(Nf)SU(N_{f}) transformation of the expectation value for an arbitrary local operator 𝒪(x1){\cal O}(x_{1}) in the path integral formalism yields the ACWTIs:

d4xlimαA0δ𝒪(x1)δαA(x)θ+d4x𝒪(x1)iDμj5Aμ(x)θ+d4x𝒪(x1)q¯{𝐦f,TA}γ5q(x)θ=0,\displaystyle\int d^{4}x\lim_{\alpha^{A}\to 0}\left\langle\frac{\delta{\cal O}^{\prime}(x_{1})}{\delta\alpha^{A}(x)}\right\rangle_{\theta}+\int d^{4}x\left\langle{\cal O}(x_{1})iD_{\mu}j_{5}^{A\mu}(x)\right\rangle_{\theta}+\int d^{4}x\left\langle{\cal O}(x_{1})\bar{q}\left\{{\bf m}_{f},T^{A}\right\}\gamma_{5}q(x)\right\rangle_{\theta}=0\,, (31)

where j5Aμj_{5}^{A\mu} denotes the chiral SU(Nf)SU(N_{f}) current, j5Aμ=q¯γμγ5TAqj_{5}^{A\mu}=\bar{q}\gamma^{\mu}\gamma_{5}T^{A}q, and

Dμj5Aμ\displaystyle D_{\mu}j_{5}^{A\mu} =\displaystyle= μj5Aμi[𝒱μ,j5Aμ],\displaystyle\partial_{\mu}j_{5}^{A\mu}-i[{\cal V}_{\mu},j_{5}^{A\mu}]\,, (32)

with 𝒱μ=(μ+μ)/2{\cal V}_{\mu}=({\cal R}_{\mu}+{\cal L}_{\mu})/2. Equation (31) is a generalization of the ACWTIs addressed in  Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr without external gauge fields. Unless external gauge fields possess a topologically nontrivial configuration, we can rewrite the covariant derivative term as

d4x𝒪(x1)iDμj5Aμ(x)θ=id4xDμ(x)0|T𝒪(x1)j5Aμ(x)|0θ=surface term.\displaystyle\int d^{4}x\left\langle{\cal O}(x_{1})iD_{\mu}j_{5}^{A\mu}(x)\right\rangle_{\theta}=i\int d^{4}xD_{\mu}^{(x)}\left\langle 0|T{\cal O}(x_{1})j_{5}^{A\mu}(x)|0\right\rangle_{\theta}=\textrm{surface term}\,. (33)

Thus we eventually have

d4xlimαA0δ𝒪(x1)δαA(x)θ+d4x𝒪(x1)q¯{𝐦f,TA}γ5q(x)θ=0.\displaystyle\int d^{4}x\lim_{\alpha^{A}\to 0}\left\langle\frac{\delta{\cal O}^{\prime}(x_{1})}{\delta\alpha^{A}(x)}\right\rangle_{\theta}+\int d^{4}x\left\langle{\cal O}(x_{1})\bar{q}\left\{{\bf m}_{f},T^{A}\right\}\gamma_{5}q(x)\right\rangle_{\theta}=0\,. (34)

This implies that the ACWTIs keep the same form as those in the case without external gauge fields Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr .

For instance, in the case of Nf=3N_{f}=3, we find the same form of the ACWTIs as in the literature GomezNicola:2016ssy ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr :

u¯uθ+d¯dθ\displaystyle\langle\bar{u}u\rangle_{\theta}+\langle\bar{d}d\rangle_{\theta} =mlχπ(θ),\displaystyle=-m_{l}\chi_{\pi}(\theta)\,,
u¯uθ+d¯dθ+4s¯sθ\displaystyle\langle\bar{u}u\rangle_{\theta}+\langle\bar{d}d\rangle_{\theta}+4\langle\bar{s}s\rangle_{\theta} =[ml(χPuu+χPdd+2χPud)2(ms+ml)(χPus+χPds)+4msχPss]θ,\displaystyle=-\left[m_{l}\left(\chi_{P}^{uu}+\chi_{P}^{dd}+2\chi_{P}^{ud}\right)-2(m_{s}+m_{l})\left(\chi_{P}^{us}+\chi_{P}^{ds}\right)+4m_{s}\chi_{P}^{ss}\right]_{\theta}\,,
u¯uθ+d¯dθ2s¯sθ\displaystyle\langle\bar{u}u\rangle_{\theta}+\langle\bar{d}d\rangle_{\theta}-2\langle\bar{s}s\rangle_{\theta} =[ml(χPuu+χPdd+2χPud)+(ml2ms)(χPus+χPds)2msχPss]θ,\displaystyle=-\left[m_{l}\left(\chi_{P}^{uu}+\chi_{P}^{dd}+2\chi_{P}^{ud}\right)+(m_{l}-2m_{s})\left(\chi_{P}^{us}+\chi_{P}^{ds}\right)-2m_{s}\chi_{P}^{ss}\right]_{\theta}\,, (35)

with ml=mu=mdm_{l}=m_{u}=m_{d}. Here χπ(θ)\chi_{\pi}(\theta) denotes the pion susceptibility defined as

χπ(θ)\displaystyle\chi_{\pi}(\theta) =Td4x[(u¯(0)iγ5u(0))(u¯(x)iγ5u(x))conn+(d¯(0)iγ5d(0))(d¯(x)iγ5d(x))conn]θ,\displaystyle=\int_{T}d^{4}x\left[\langle(\bar{u}(0)i\gamma_{5}u(0))(\bar{u}(x)i\gamma_{5}u(x))\rangle_{\rm conn}+\langle(\bar{d}(0)i\gamma_{5}d(0))(\bar{d}(x)i\gamma_{5}d(x))\rangle_{\rm conn}\right]_{\theta}\,, (36)

with conn\langle\cdot\cdot\cdot\rangle_{\rm conn} being the connected part of the correlation function, and the pseudoscalar susceptibilities χPuu,dd,ud\chi_{P}^{uu,dd,ud}, χPss\chi_{P}^{ss} and χPus,ds\chi_{P}^{us,ds} are defined as

χPf1f2(θ)\displaystyle\chi_{P}^{f_{1}f_{2}}(\theta) =Td4x(q¯f1(0)iγ5qf1(0))(q¯f2(x)iγ5qf2(x))θ,forqf1,2=u,d,s.\displaystyle=\int_{T}d^{4}x\langle(\bar{q}_{f_{1}}(0)i\gamma_{5}q_{f_{1}}(0))(\bar{q}_{f_{2}}(x)i\gamma_{5}q_{f_{2}}(x))\rangle_{\theta}\,,\qquad{\rm for}\quad{q_{{}_{f_{1,2}}}=u,d,s}\,. (37)

χtop(0)(θ)\chi_{\rm top}^{(0)}(\theta) in Eq.(29) then takes a couple of equivalent forms, related each other by the ACWTIs in Eq.(35Kawaguchi:2020kdl ; Kawaguchi:2020qvg ; Cui:2021bqf ; Cui:2022vsr

χtop(0)(θ)\displaystyle\chi_{\rm top}^{(0)}(\theta) =m¯2[u¯uml+d¯dml+s¯sms+χPuu+χPdd+χPss+2χPud+2χPus+2χPds]θ\displaystyle=\bar{m}^{2}\left[\frac{\langle\bar{u}u\rangle}{m_{l}}+\frac{\langle\bar{d}d\rangle}{m_{l}}+\frac{\langle\bar{s}s\rangle}{m_{s}}+\chi_{P}^{uu}+\chi_{P}^{dd}+\chi_{P}^{ss}+2\chi_{P}^{ud}+2\chi_{P}^{us}+2\chi_{P}^{ds}\right]_{\theta}
=14[ml(u¯u+d¯d)+ml2(χPuu+χPdd+2χPud)]θ\displaystyle=\frac{1}{4}\left[m_{l}\left(\langle\bar{u}u\rangle+\langle\bar{d}d\rangle\right)+m_{l}^{2}\left(\chi_{P}^{uu}+\chi_{P}^{dd}+2\chi_{P}^{ud}\right)\right]_{\theta}
=ml24(χπ(θ)χη(θ))\displaystyle=\frac{m_{l}^{2}}{4}\left(\chi_{\pi}(\theta)-\chi_{\eta}(\theta)\right)
=mss¯sθ+ms2χPss(θ),\displaystyle=m_{s}\langle\bar{s}s\rangle_{\theta}+m_{s}^{2}\chi_{P}^{ss}(\theta)\,, (38)

with m¯=(2ml+1ms)1\bar{m}=\left(\frac{2}{m_{l}}+\frac{1}{m_{s}}\right)^{-1}. It is interesting to note that χπ\chi_{\pi} and χη\chi_{\eta} are related not to the full topological susceptibility χtop\chi_{\rm top}, but χtop(0)\chi_{\rm top}^{(0)}, in Eq.(26).

A.2 Renormalization group invariance of χtop\chi_{\rm top}

The current quark mass parameter is multiplicatively renormalized as

mf0(Λ)=ZS1(Λ,μ)mfR(μ),\displaystyle m_{f}^{0}(\Lambda)=Z_{S}^{-1}(\Lambda,\mu)m_{f}^{R}(\mu)\,, (39)

where Λ\Lambda denotes the bare cutoff and the upper script RR stands for the quantity renormalized at the scale μ\mu. When the mass independent renormalization is thus applied. the renormalization factor ZS1Z_{S}^{-1} is flavor universal. Then we are allowed to drop the flavor label ff as

m¯0(Λ)=ZS1(Λ,μ)m¯R(μ).\displaystyle\bar{m}^{0}(\Lambda)=Z_{S}^{-1}(\Lambda,\mu)\bar{m}^{R}(\mu). (40)

This also implies that ZSZ_{S} is independent of the current quark mass as well.

The q¯fqf\bar{q}_{f}q_{f} bilinear operator is also multiplicatively renormalized by ZS(λ,μ)Z_{S}(\lambda,\mu),

(q¯fqf)R=ZS1(Λ,μ)(q¯fqf)Λ,\displaystyle\left(\bar{q}_{f}q_{f}\right)_{R}=Z_{S}^{-1}(\Lambda,\mu)\left(\bar{q}_{f}q_{f}\right)_{\Lambda}\,, (41)

so that we have the renormalization group invariant mass term like

mf0(Λ)(q¯fqf)Λ=mfR(μ)(q¯fqf)R.\displaystyle m_{f}^{0}(\Lambda)\left(\bar{q}_{f}q_{f}\right)_{\Lambda}=m_{f}^{R}(\mu)\left(\bar{q}_{f}q_{f}\right)_{R}\,. (42)

Consider also a (q¯fiγ5qf\bar{q}_{f}i\gamma_{5}q_{f}) operator to be renormalized in a similar way with the renormalization constant ZP1(λ,μ)Z_{P}^{-1}(\lambda,\mu)

(q¯fiγ5qf)R=ZP1(Λ,μ)(q¯fiγ5qf)Λ.\displaystyle(\bar{q}_{f}i\gamma_{5}q_{f})_{R}=Z_{P}^{-1}(\Lambda,\mu)(\bar{q}_{f}i\gamma_{5}q_{f})_{\Lambda}. (43)

Since none of quark mass dependence is generated in the mass independent renormalization, the U(1)U(1) axial invariance keeps manifest between renormalization of the (q¯fqf\bar{q}_{f}q_{f}) and (q¯fiγ5qf\bar{q}_{f}i\gamma_{5}q_{f}) operators. Hence we have

ZP1=ZS1.\displaystyle Z_{P}^{-1}=Z_{S}^{-1}\,. (44)

In that case we also find

(χPf1f2)Λ=ZS2(χPf1f2)R,\displaystyle\left(\chi_{P}^{f_{1}f_{2}}\right)_{\Lambda}=Z_{S}^{2}\left(\chi_{P}^{f_{1}f_{2}}\right)_{R}\,, (45)

Now, we apply the renormalization procedure as above to χtop(0)\chi_{\rm top}^{(0)} in the three-flavor case (Eq.(38)):

(χtop(0))Λ\displaystyle\left(\chi_{\rm top}^{(0)}\right)_{\Lambda} =\displaystyle= m¯02[(u¯u)Λmu0+(d¯d)Λmd0+(s¯s)Λms0+i(χPuu)Λ+i(χPdd)Λ+i(χPss)Λ+2i(χPud)Λ+2i(χPus)Λ+2i(χPds)Λ]\displaystyle-\bar{m}_{0}^{2}\left[\frac{\langle(\bar{u}u)_{\Lambda}\rangle}{m_{u}^{0}}+\frac{\langle(\bar{d}d)_{\Lambda}\rangle}{m_{d}^{0}}+\frac{\langle(\bar{s}s)_{\Lambda}\rangle}{m_{s}^{0}}+i(\chi_{P}^{uu})_{\Lambda}+i(\chi_{P}^{dd})_{\Lambda}+i(\chi_{P}^{ss})_{\Lambda}+2i(\chi_{P}^{ud})_{\Lambda}+2i(\chi_{P}^{us})_{\Lambda}+2i(\chi_{P}^{ds})_{\Lambda}\right] (46)
=\displaystyle= m¯R2[(u¯u)RmuR+(d¯d)RmdR+(s¯s)RmsR+i(χPuu)R+i(χPdd)R+i(χPss)R+2i(χPud)R+2i(χPus)R+2i(χPds)R]\displaystyle-\bar{m}_{R}^{2}\left[\frac{\langle(\bar{u}u)_{R}\rangle}{m_{u}^{R}}+\frac{\langle(\bar{d}d)_{R}\rangle}{m_{d}^{R}}+\frac{\langle(\bar{s}s)_{R}\rangle}{m_{s}^{R}}+i(\chi_{P}^{uu})_{R}+i(\chi_{P}^{dd})_{R}+i(\chi_{P}^{ss})_{R}+2i(\chi_{P}^{ud})_{R}+2i(\chi_{P}^{us})_{R}+2i(\chi_{P}^{ds})_{R}\right]
=\displaystyle= (χtop(0))R.\displaystyle\left(\chi_{\rm top}^{(0)}\right)_{R}\,.

Thus it has been proven that χtop(0)\chi_{\rm top}^{(0)} is renormalization group invariant. Note that this argument is also applicable to the case with nonzero θ\theta, because the QCD θ\theta is itself conventionally renormalization group invariant and no new divergent terms induced due to nonzero θ\theta will be generated, hence the renormalization factors will not be corrected. Furthermore, one can readily see that the external gauge-induced χtopEM\chi_{\rm top}^{\rm EM}, defined as in Eq.(27), is also manifestly renormalization group invariant.

Appendix B The details on the (P) NJL model analysis

B.1 NJL case

Our reference NJL Lagrangian with two flavors (up and down quarks) and nonzero θ\theta follows the literature Sakai:2011gs , which takes the form

=q¯(iγμμm)q+gs2a=03[(q¯τaq)2+(q¯iγ5τaq)2]+gd(eiθdet[q¯(1+γ5)q]+eiθdet[q¯(1γ5)q]),\displaystyle{\cal L}=\bar{q}(i\gamma_{\mu}\partial^{\mu}-m)q+\frac{g_{s}}{2}\sum_{a=0}^{3}\left[(\bar{q}\tau_{a}q)^{2}+(\bar{q}i\gamma_{5}\tau_{a}q)^{2}\right]+g_{d}\left(e^{i\theta}\mathop{\rm det}[\bar{q}(1+\gamma_{5})q]+e^{-i\theta}\mathop{\rm det}[\bar{q}(1-\gamma_{5})q]\right)\,, (47)

where q=(u,d)Tq=(u,d)^{T} and the determinant acts on the quark flavors. The gdg_{d} term, called the ‘t Hooft-Kobayashi-Maskawa determinant term Kobayashi:1970ji ; Kobayashi:1971qz ; tHooft:1976rip ; tHooft:1976snw , breaks U(1)U(1) axial symmetry, but keeps SU(2)L×SU(2)R×U(1)VSU(2)_{L}\times SU(2)_{R}\times U(1)_{V} symmetry, while others keep the full chiral U(2)L×U(2)RU(2)_{L}\times U(2)_{R} symmetry. The gdg_{d} term thus serves as the U(1)U(1) axial anomaly:

μjAμ=2iq¯mγ5q8gdIm[detq¯(1γ5)qeiθ].\displaystyle\partial_{\mu}j^{\mu}_{A}=2i\bar{q}m\gamma_{5}q-8g_{d}{\rm Im}\left[\mathop{\rm det}\bar{q}(1-\gamma_{5})q\cdot e^{-i\theta}\right]. (48)

Since the U(1)U(1) axial symmetry is broken only by gdg_{d} and the quark mass terms, one can move θ\theta in the gdg_{d} term to the mass term by a U(1)U(1) axial rotation as was done in the case of the general argument above:

qeiγ5θ2qq,\displaystyle q\to e^{-i\gamma_{5}\frac{\theta}{2}}q\equiv q^{\prime}\,, (49)

so that

=q¯(iγμμm(θ))q+gs2a=03[(q¯τaq)2+(q¯iγ5τaq)2]+gd(det[q¯(1+γ5)q]+h.c.),\displaystyle{\cal L}\to{\cal L}^{\prime}=\bar{q}^{\prime}(i\gamma_{\mu}\partial^{\mu}-m(\theta))q^{\prime}+\frac{g_{s}}{2}\sum_{a=0}^{3}\left[(\bar{q}^{\prime}\tau_{a}q^{\prime})^{2}+(\bar{q}^{\prime}i\gamma_{5}\tau_{a}q^{\prime})^{2}\right]+g_{d}\left(\mathop{\rm det}[\bar{q}^{\prime}(1+\gamma_{5})q^{\prime}]+{\rm h.c.}\right)\,, (50)

where

m(θ)=m[cosθ2+iγ5sinθ2],\displaystyle m(\theta)=m\left[\cos\frac{\theta}{2}+i\gamma_{5}\sin\frac{\theta}{2}\right]\,, (51)

in which nonzero θ\theta manifestly signals the CP violation. Use of this ”prime” basis is convenient to analyze the model because all the θ\theta dependence is transformed and collected into the complex mass m(θ)m(\theta) in the quark propagator. Similarly to Eq.(1), Eq.(49) relates the scalar and pseudoscalar bilinears between the original- and prime-base scalar and pseudoscalar bilinears (for each quark flavor ii) as

(q¯iqi)\displaystyle(\bar{q}_{i}q_{i}) =(q¯iqi)cosθ2+(q¯iiγ5qi)sinθ2,\displaystyle=(\bar{q}_{i}^{\prime}q_{i}^{\prime})\cos\frac{\theta}{2}+(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime})\sin\frac{\theta}{2}\,,
(q¯iiγ5qi)\displaystyle(\bar{q}_{i}i\gamma_{5}q_{i}) =(q¯iqi)sinθ2+(q¯iiγ5qi)cosθ2.\displaystyle=-(\bar{q}_{i}^{\prime}q_{i}^{\prime})\sin\frac{\theta}{2}+(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime})\cos\frac{\theta}{2}\,. (52)

We work in the MFA, so that the scalar and pseudoscalar bilinears (q¯iqi)(\bar{q}_{i}^{\prime}q_{i}^{\prime}) and (q¯iiγ5qi)(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime}) are expanded around the means fields S=q¯iqiS^{\prime}=\langle\bar{q}_{i}^{\prime}q_{i}^{\prime}\rangle and P=q¯iiγ5qiP^{\prime}=\langle\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime}\rangle, as q¯iqi=S+(:q¯iqi:)\bar{q}_{i}^{\prime}q_{i}^{\prime}=S^{\prime}+(:\bar{q}_{i}^{\prime}q_{i}^{\prime}:) and q¯iiγ5qi=P+(:q¯iiγ5qi:)\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime}=P^{\prime}+(:\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime}:), where the terms sandwiched by “::” stand for the normal ordered product, meaning that :𝒪:=0\langle:{\cal O}:\rangle=0 for 𝒪=S,P{\cal O}=S^{\prime},P^{\prime}. Then the interaction terms in Eq.(47) are replaced, up to the normal ordered terms, as

(q¯iqi)2\displaystyle(\bar{q}_{i}^{\prime}q_{i}^{\prime})^{2} 4Si(q¯iqi)4S2,\displaystyle\to 4S^{\prime}\sum_{i}(\bar{q}_{i}^{\prime}q_{i}^{\prime})-4S^{\prime 2}\,,
(q¯iiγ5qi)2\displaystyle(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime})^{2} 4Pi(q¯iiγ5qi)4P2,\displaystyle\to 4P^{\prime}\sum_{i}(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime})-4P^{\prime 2}\,,
det[q¯i(1+γ5)qi]+h.c.\displaystyle{\rm det}[\bar{q}_{i}^{\prime}(1+\gamma_{5})q_{i}^{\prime}]+{\rm h.c.} 2Si(q¯iqi)2Pi(q¯iiγ5qi)2(S2P2).\displaystyle\to 2S^{\prime}\sum_{i}(\bar{q}_{i}^{\prime}q_{i}^{\prime})-2P^{\prime}\sum_{i}(\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime})-2(S^{\prime 2}-P^{\prime 2})\,. (53)

In the MFA the NJL Lagrangian in Eq.(47) thus takes the form

MFA=iq¯i(iγμμ(S,P;θ))qi2gs(S2+P2)2gd(S2P2),\displaystyle{\cal L}_{\rm MFA}=\sum_{i}\bar{q}^{\prime}_{i}(i\gamma_{\mu}\partial^{\mu}-{\cal M}(S^{\prime},P^{\prime};\theta))q_{i}^{\prime}-2g_{s}(S^{\prime 2}+P^{\prime 2})-2g_{d}(S^{\prime 2}-P^{\prime 2})\,, (54)

with

(S,P;θ)\displaystyle{\cal M}(S^{\prime},P^{\prime};\theta) =α(S;θ)+iγ5β(P;θ),\displaystyle=\alpha(S^{\prime};\theta)+i\gamma_{5}\beta(P^{\prime};\theta)\,,
α(S;θ)\displaystyle\alpha(S^{\prime};\theta) =mcosθ22(gs+gd)S,\displaystyle=m\cos\frac{\theta}{2}-2(g_{s}+g_{d})S^{\prime}\,,
β(P;θ)\displaystyle\beta(P^{\prime};\theta) =msinθ22(gsgd)P.\displaystyle=m\sin\frac{\theta}{2}-2(g_{s}-g_{d})P^{\prime}\,. (55)

Integrating out quarks leads to the thermodynamic potential in the MFA:

Ω[S,P;θ]\displaystyle\Omega[S^{\prime},P^{\prime};\theta] =2gs(S2+P2)+2gd(S2P2)2NcNfd3k(2π)3[E+2Tln(1+eE/T)],\displaystyle=2g_{s}(S^{\prime 2}+P^{\prime 2})+2g_{d}(S^{\prime 2}-P^{\prime 2})-2N_{c}N_{f}\int\frac{d^{3}k}{(2\pi)^{3}}\left[E+2T\ln(1+e^{-E/T})\right]\,, (56)

where Nf=2N_{f}=2 and Nc=3N_{c}=3, and

E=M2+k2,M2=α2+β2,\displaystyle E=\sqrt{M^{2}+k^{2}}\,,\qquad M^{2}=\alpha^{2}+\beta^{2}\,, (57)

with k2=|k|2k^{2}=|\vec{k}|^{2}. Then q¯iqi=S\langle\bar{q}_{i}^{\prime}q_{i}^{\prime}\rangle=S^{\prime} and q¯iiγ5qi=P\langle\bar{q}_{i}^{\prime}i\gamma_{5}q_{i}^{\prime}\rangle=P^{\prime} are determined through the stationary condition,

ΩS=ΩP=0.\displaystyle\frac{\partial\Omega}{\partial S^{\prime}}=\frac{\partial\Omega}{\partial P^{\prime}}=0\,. (58)

Of particular interest is to see the thermodynamic potential at θ=π\theta=\pi,

Ω[α,β;θ=π]\displaystyle\Omega[\alpha,\beta;\theta=\pi] =\displaystyle= 12gs+2gDα2+12gs2gD(βm)22NcNfd3p(2π)3[E+2Tln(1+eE/T)],\displaystyle\frac{1}{2g_{s}+2g_{D}}\alpha^{2}+\frac{1}{2g_{s}-2g_{D}}(\beta-m)^{2}-2N_{c}N_{f}\int\frac{d^{3}p}{(2\pi)^{3}}\left[E+2T\ln\left(1+e^{-E/T}\right)\right]\,, (59)

where things have been written in terms of α\alpha and β\beta by means of Eq.(55). In the presently applied MFA, the loop corrections (corresponding to the last term in Eq.(59)) are U(1)U(1) axial invariant, i.e., do not separate α\alpha and β\beta due to no loop contributions involving vertices with the determinant coupling gdg_{d}. At θ=π\theta=\pi, the explicit-chiral SU(2)SU(2) breaking-effect (by mm) has completely been transported into the β\beta direction (the second term). The stationary condition then takes the form

α\displaystyle\alpha =\displaystyle= 2NcNf(gs+gD)αd3p(2π)31E[121+eE/T],\displaystyle 2N_{c}N_{f}(g_{s}+g_{D})\alpha\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{E}\left[1-\frac{2}{1+e^{E/T}}\right],
βm\displaystyle\beta-m =\displaystyle= 2NcNf(gsgD)βd3p(2π)31E[121+eE/T].\displaystyle 2N_{c}N_{f}(g_{s}-g_{D})\beta\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{E}\left[1-\frac{2}{1+e^{E/T}}\right]\,. (60)

When α0\alpha\neq 0, i.e., the CP symmetry is spontaneously broken, the gap equation for β\beta can be rewritten by eliminating the loop correction part as

βm=gsgDgs+gDβ\displaystyle\beta-m=\frac{g_{s}-g_{D}}{g_{s}+g_{D}}\beta (61)
\displaystyle\leftrightarrow β=gs+gD2gDm\displaystyle\beta=\frac{g_{s}+g_{D}}{2g_{D}}m
\displaystyle\leftrightarrow P=m4gD.\displaystyle P^{\prime}=-\frac{m}{4g_{D}}\,.

This implies that in the CP broken phase at θ=π\theta=\pi, the chiral SU(2)SU(2) symmetry is not spontaneously broken and the chiral order parameter does not evolve in TT.

We come back to the case with arbitrary θ\theta and derive relevant formulae for the susceptibilities χη\chi_{\eta} and χtop\chi_{\rm top} within the MFA. First, by noting the change of basis in Eq.(52), the η\eta susceptibility χη\chi_{\eta} in χtop\chi_{\rm top} of Eq.(4) is written in terms of the prime base as

χη=cos2θ2χη2sinθ2cosθ2χησ+sin2θ2χσ,\displaystyle\chi_{\eta}=\cos^{2}\frac{\theta}{2}\chi^{\prime}_{\eta}-2\sin\frac{\theta}{2}\cos\frac{\theta}{2}\chi_{\eta\sigma}^{\prime}+\sin^{2}\frac{\theta}{2}\chi^{\prime}_{\sigma}\,, (62)

The “primed” meson susceptibilities are evaluated in the so-called random phase approximation (the bubble-ring resummation ansatz) as

𝐗\displaystyle{\bf X}^{\prime} =𝚷1𝟏2×2+𝐆𝚷,\displaystyle={\bf\Pi}^{\prime}\cdot\frac{1}{{\bf 1}_{2\times 2}+{\bf G}^{\prime}\cdot{\bf\Pi}^{\prime}}\,, (63)

where

𝐗\displaystyle{\bf X}^{\prime} =(χσχησχησχη),\displaystyle=\left(\begin{array}[]{cc}\chi_{\sigma}^{\prime}&\chi_{\eta\sigma}^{\prime}\\ \chi_{\eta\sigma}^{\prime}&\chi_{\eta}^{\prime}\end{array}\right)\,, (66)
𝚷\displaystyle{\bf\Pi}^{\prime} =(ΠσΠησΠησΠη),\displaystyle=\left(\begin{array}[]{cc}\Pi_{\sigma}^{\prime}&\Pi_{\eta\sigma}^{\prime}\\ \Pi_{\eta\sigma}^{\prime}&\Pi_{\eta}^{\prime}\end{array}\right)\,, (69)
𝐆\displaystyle{\bf G}^{\prime} =(GσGησGησGη)=(gs+gd00gsgd),\displaystyle=\left(\begin{array}[]{cc}G_{\sigma}^{\prime}&G_{\eta\sigma}^{\prime}\\ G_{\eta\sigma}^{\prime}&G_{\eta}^{\prime}\end{array}\right)=\left(\begin{array}[]{cc}g_{s}+g_{d}&0\\ 0&g_{s}-g_{d}\end{array}\right)\,, (74)

with the vacuum polarization functions for each channel,

Πσ\displaystyle\Pi_{\sigma}^{\prime} =2IS,\displaystyle=2I_{S^{\prime}}\,,
Πησ\displaystyle\Pi_{\eta\sigma} =ISP,\displaystyle=I_{S^{\prime}P^{\prime}}\,,
Πη\displaystyle\Pi_{\eta}^{\prime} =2IP,\displaystyle=2I_{P^{\prime}}\,,
IS\displaystyle I_{S^{\prime}} =Ncπ20Λ𝑑kk2E2α2E3[12eE/T+1],\displaystyle=-\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{E^{2}-\alpha^{2}}{E^{3}}\left[1-\frac{2}{e^{E/T}+1}\right]\,,
ISP\displaystyle I_{S^{\prime}P^{\prime}} =Ncπ20Λ𝑑kk22αβE3[12eE/T+1].\displaystyle=\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{2\alpha\beta}{E^{3}}\left[1-\frac{2}{e^{E/T}+1}\right]\,.
IP\displaystyle I_{P^{\prime}} =Ncπ20Λ𝑑kk2E2β2E3[12eE/T+1].\displaystyle=-\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{E^{2}-\beta^{2}}{E^{3}}\left[1-\frac{2}{e^{E/T}+1}\right]\,. (75)

In evaluating the vacuum polarization functions, we have regularized the 3-momentum integrals by the cutoff Λ\Lambda.

Putting all the relevant things above into χtop\chi_{\rm top} in Eq.(4), we compute χtop\chi_{\rm top} as a function of TT and θ\theta. The model parameter setting follows the literature Boer:2008ct ; Boomsma:2009eh ; Sakai:2011gs :

m\displaystyle m =6MeV,\displaystyle=6\,{\rm MeV}\,,
Λ\displaystyle\Lambda =590MeV,\displaystyle=590\,{\rm MeV}\,,
gs\displaystyle g_{s} =2(1c)G0,gd=2cG0withG0Λ2=2.435andc=0.2.\displaystyle=2(1-c)G_{0}\,\,,\qquad g_{d}=2cG_{0}\,\qquad{\rm with}\qquad G_{0}\Lambda^{2}=2.435\qquad{\rm and}\qquad c=0.2\,. (76)

The gdg_{d} coupling has been related by cc to the gsg_{s} coupling just in a numerical manner, though the associated asymmetry features are different.

Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 3: Scalar and pseudoscalar condensates normalized to the scalar condensate at T=θ=0T=\theta=0 versus temperature at nonzero θ\theta.

See Fig. 3, where we observe several characteristic features: i) the scalar condensate becomes dramatically smaller with increasing θ\theta, so that χtop\chi_{\rm top} in Eq.(4) gigantically decreases; ii) at θ=π\theta=\pi, the pseudoscalar condensate undergoes the CP symmetry restoration of the second order type, which generates a significant spike at the criticality in χtop\chi_{\rm top}; iii) at θ=π\theta=\pi, the scalar condensate does not evolve in TT at all until the the CP symmetry is restored, in agreement with the analytic discussion around Eq.(61).

B.2 PNJL case

By extending the NJL detailed in the previous subsection into the PNJL model in the MFA following the procedure in the literature Fukushima:2003fw (for a recent review, see, e.g., Fukushima:2017csk ), we have the thermodynamic potential:

ΩPNJL[S,P,L,L]\displaystyle\Omega_{\rm PNJL}[S^{\prime},P^{\prime},L,L^{\dagger}] =2gs(S2+P2)+2gd(S2P2)\displaystyle=2g_{s}(S^{\prime 2}+P^{\prime 2})+2g_{d}(S^{\prime 2}-P^{\prime 2})
2Nfd3k(2π)3trcolor[E+T(ln(1+LeE/T)+h.c.)]+𝒰(L,L)\displaystyle-2N_{f}\int\frac{d^{3}k}{(2\pi)^{3}}{\rm tr}_{\rm color}\left[E+T\left(\ln(1+L\cdot e^{-E/T})+{\rm h.c.}\right)\right]+{\cal U}(L,L^{\dagger})
=2gs(S2+P2)+2gd(S2P2)\displaystyle=2g_{s}(S^{\prime 2}+P^{\prime 2})+2g_{d}(S^{\prime 2}-P^{\prime 2})
2NcNfd3k(2π)3[E+T3(ln(1+3ΦeE/T+3Φe2E/T+e3E/T)+h.c.)]+𝒰(Φ,Φ),\displaystyle-2N_{c}N_{f}\int\frac{d^{3}k}{(2\pi)^{3}}\left[E+\frac{T}{3}\left(\ln(1+3\Phi e^{-E/T}+3\Phi^{*}e^{-2E/T}+e^{-3E/T})+{\rm h.c.}\right)\right]+{\cal U}(\Phi,\Phi^{*})\,, (77)

where Φ=tr[L]/3\Phi={\rm tr}[L]/3, in which we have introduced the Polyakov loop field L=diag{eiϕ1,eiϕ2,eiϕ3}L={\rm diag}\{e^{i\phi_{1}},e^{i\phi_{2}},e^{i\phi_{3}}\} with the Z3Z_{3} charges ϕ1,2,3\phi_{1,2,3}; the SU(3)SU(3) constraint det[L]=1{\rm det}[L]=1 taken into account; EE takes the same form as in Eq.(57), but would depend on the Polyakov-loop fields Φ\Phi and Φ\Phi^{*} through the stationary condition as in Eq.(58). The Polyakov loop potential 𝒰(Φ,Φ){\cal U}(\Phi,\Phi^{*}) is assumed to take the form of polynomial type or logarithmic type:

𝒰poly(Φ,Φ)\displaystyle{\cal U}_{\rm poly}(\Phi,\Phi^{*}) =T4[b2(T)2ΦΦb36(Φ3+Φ3)+b44(ΦΦ)2],\displaystyle=T^{4}\left[-\frac{b_{2}(T)}{2}\Phi^{*}\Phi-\frac{b_{3}}{6}\left(\Phi^{*3}+\Phi^{3}\right)+\frac{b_{4}}{4}(\Phi^{*}\Phi)^{2}\right]\,,
𝒰log(Φ,Φ)\displaystyle{\cal U}_{\rm log}(\Phi,\Phi^{*}) =T4[c(T)2ΦΦ+d(T)ln(16ΦΦ+4(Φ3+Φ3)3(ΦΦ)2)],\displaystyle=T^{4}\left[-\frac{c(T)}{2}\Phi^{*}\Phi+d(T)\ln\left(1-6\Phi\Phi^{*}+4\left(\Phi^{3}+\Phi^{*3}\right)-3\left(\Phi\Phi^{*}\right)^{2}\right)\right]\,, (78)

with

b2(T)\displaystyle b_{2}(T) =a0+a1(T0T)+a2(T0T)2+a3(T0T)3,\displaystyle=a_{0}+a_{1}\left(\frac{T_{0}}{T}\right)+a_{2}\left(\frac{T_{0}}{T}\right)^{2}+a_{3}\left(\frac{T_{0}}{T}\right)^{3}\,,
c(T)\displaystyle c(T) =c0+c1(T0T)+c2(T0T)2,\displaystyle=c_{0}+c_{1}\left(\frac{T_{0}}{T}\right)+c_{2}\left(\frac{T_{0}}{T}\right)^{2}\,,
d(T)\displaystyle d(T) =d3(T0T)3.\displaystyle=d_{3}\left(\frac{T_{0}}{T}\right)^{3}\,. (79)

The Polyakov-loop potential parameters are fixed by fitting the lattice data in the pure Yang-Mills theory as Roessner:2006xn ; Ratti:2005jh

a0a1a2a3b3b4c0c1c2d3T0[MeV]6.751.952.6257.440.757.53.512.4715.21.75270.\displaystyle\begin{array}[]{|c|c|c|c|c|c|c|c|c|c|| c|}\hline\cr a_{0}&a_{1}&a_{2}&a_{3}&b_{3}&b_{4}&c_{0}&c_{1}&c_{2}&d_{3}&T_{0}[{\rm MeV}]\\ \hline\cr 6.75&-1.95&2.625&-7.44&0.75&7.5&3.51&-2.47&15.2&-1.75&270\\ \hline\cr\end{array}\,. (82)

The thermodynamic potential in Eq.(77) is thus minimized also with respect to Φ\Phi and Φ\Phi^{*}, in addition to the stationary condition along SS^{\prime} and PP^{\prime} directions as in the NJL case, Eq.(76).

In the RPA the following replacement rule for IPI_{P^{\prime}} in Eq.(75) due to the Polyakov-loop contribution is applied:

IPIP(Φ)=Ncπ20Λdkk2E2β2E3[1(eE/T(Φ+2ΦeE/T+e2E/T)1+3ΦeE/T+3Φe2E/T+e3E/T+h.c.)].\displaystyle I_{P^{\prime}}\to I_{P^{\prime}}(\Phi)=-\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{E^{2}-\beta^{2}}{E^{3}}\Bigg{[}1-\left(\frac{e^{-E/T}(\Phi+2\Phi^{*}e^{-E/T}+e^{-2E/T})}{1+3\Phi e^{-E/T}+3\Phi^{*}e^{-2E/T}+e^{-3E/T}}+{\rm h.c.}\right)\Bigg{]}\,. (83)

Similarly for ISI_{S^{\prime}} an ISPI_{S^{\prime}P^{\prime}} in Eq.(75), we have

ISIS(Φ)\displaystyle I_{S^{\prime}}\to I_{S^{\prime}}(\Phi) =Ncπ20Λdkk2E2α2E3[1(eE/T(Φ+2ΦeE/T+e2E/T)1+3ΦeE/T+3Φe2E/T+e3E/T+h.c.)],\displaystyle=-\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{E^{2}-\alpha^{2}}{E^{3}}\Bigg{[}1-\left(\frac{e^{-E/T}(\Phi+2\Phi^{*}e^{-E/T}+e^{-2E/T})}{1+3\Phi e^{-E/T}+3\Phi^{*}e^{-2E/T}+e^{-3E/T}}+{\rm h.c.}\right)\Bigg{]}\,,
ISPISP(Φ)\displaystyle I_{S^{\prime}P^{\prime}}\to I_{S^{\prime}P^{\prime}}(\Phi) =Ncπ20Λdkk22αβE3[1(eE/T(Φ+2ΦeE/T+e2E/T)1+3ΦeE/T+3Φe2E/T+e3E/T+h.c.)].\displaystyle=\frac{N_{c}}{\pi^{2}}\int_{0}^{\Lambda}dkk^{2}\frac{2\alpha\beta}{E^{3}}\Bigg{[}1-\left(\frac{e^{-E/T}(\Phi+2\Phi^{*}e^{-E/T}+e^{-2E/T})}{1+3\Phi e^{-E/T}+3\Phi^{*}e^{-2E/T}+e^{-3E/T}}+{\rm h.c.}\right)\Bigg{]}\,. (84)

Those modified integral functions are precisely reduced back to the ones in Eq.(75) when Φ=1\Phi=1 (without confinement).

In Fig. 4 we plot the scalar and pseudoscalar condensates, q¯q\langle\bar{q}q\rangle and q¯iγ5q\langle\bar{q}i\gamma_{5}q\rangle, as a function of temperature TT with θ/π=0,0.5,0.8\theta/\pi=0,0.5,0.8, and 1. The scalar condensate decreases sharply with increasing θ\theta. The pseudoscalar condensate undergoes a second-order phase transition at the critical temperature when θ=π\theta=\pi, and the CP symmetry is restored, which coincides with the NJL case (Fig. 3). At θ=π\theta=\pi, the scalar condensate does not evolve in TT at all until the CP symmetry is restored at the critical temperature, in the same way as in the NJL case (Fig. 3) due to the same form of the gap equations as in Eq.(61) because incorporation of the Polyakov loop field does not break the U(1)U(1) axial symmetry. This is the characteristic MFA feature, which persists both in the NJL and PNJL cases.

Those universal trends have also been seen in χtop\chi_{\rm top}, in Fig. 5, which shows no substantial difference from the case of the NJL as in Fig. 1, hence the dramatic suppression of α\alpha_{*} around T=𝒪(100)T={\cal O}(100) MeV still persists for θ=𝒪(1)\theta={\cal O}(1) even with the PL contribution taken into account (See Fig. 6). At θ=π\theta=\pi, the peak signal strengths in the PNJL model (whichever type L or P) deviates from the 2σ\sigma contour, as has also been observed in the NJL case (Fig. 2 in the main text).

Refer to caption
Refer to caption
Refer to caption
Refer to caption
Figure 4: Scalar and pseudoscalar condensates normalized to the scalar condensate at T=0T=0 and θ=0\theta=0 versus temperature with varying θ\theta from 0. The abbreviations L and P denote the PL potential of logarithmic and polynomial type, respectively.
Refer to caption
Refer to caption
Figure 5: Plot of χtop\chi_{\rm top} (in magnitude) computed from the two-flavor PNJL model (with two types of PL potentials) with ranging θ\theta from 0 to π\pi as in Fig. 4.
Refer to caption
Refer to caption
Figure 6: The signal strength (the normalized latent heat) α\alpha_{*} versus TT^{*} based on the formula in Eq.(13) with the two-flavor PNJL model estimate of χtop\chi_{\rm top} and θ\theta varied from 0 to π\pi as in Figs. 4 and 5.

References

  • (1) G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L8 (2023) doi:10.3847/2041-8213/acdac6 [arXiv:2306.16213 [astro-ph.HE]].
  • (2) G. Agazie et al. [NANOGrav], Astrophys. J. Lett. 952, no.2, L37 (2023) doi:10.3847/2041-8213/ace18b [arXiv:2306.16220 [astro-ph.HE]].
  • (3) A. Afzal et al. [NANOGrav], Astrophys. J. Lett. 951, no.1, L11 (2023) doi:10.3847/2041-8213/acdc91 [arXiv:2306.16219 [astro-ph.HE]].
  • (4) J. Antoniadis et al. [EPTA], Astron. Astrophys. 678, A48 (2023) doi:10.1051/0004-6361/202346841 [arXiv:2306.16224 [astro-ph.HE]].
  • (5) J. Antoniadis et al. [EPTA and InPTA], Astron. Astrophys. 678, A49 (2023) doi:10.1051/0004-6361/202346842 [arXiv:2306.16225 [astro-ph.HE]].
  • (6) J. Antoniadis et al. [EPTA and InPTA:], Astron. Astrophys. 678, A50 (2023) doi:10.1051/0004-6361/202346844 [arXiv:2306.16214 [astro-ph.HE]].
  • (7) D. J. Reardon, A. Zic, R. M. Shannon, G. B. Hobbs, M. Bailes, V. Di Marco, A. Kapur, A. F. Rogers, E. Thrane and J. Askew, et al. Astrophys. J. Lett. 951, no.1, L6 (2023) doi:10.3847/2041-8213/acdd02 [arXiv:2306.16215 [astro-ph.HE]].
  • (8) D. J. Reardon, A. Zic, R. M. Shannon, V. Di Marco, G. B. Hobbs, A. Kapur, M. E. Lower, R. Mandow, H. Middleton and M. T. Miles, et al. Astrophys. J. Lett. 951, no.1, L7 (2023) doi:10.3847/2041-8213/acdd03 [arXiv:2306.16229 [astro-ph.HE]].
  • (9) H. Xu, S. Chen, Y. Guo, J. Jiang, B. Wang, J. Xu, Z. Xue, R. N. Caballero, J. Yuan and Y. Xu, et al. Res. Astron. Astrophys. 23, no.7, 075024 (2023) doi:10.1088/1674-4527/acdfa5 [arXiv:2306.16216 [astro-ph.HE]].
  • (10) H. J. Li and Y. F. Zhou, [arXiv:2401.09138 [hep-ph]].
  • (11) J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki, M. Raidal, J. Urrutia, V. Vaskonen and H. Veermäe, Phys. Rev. D 109, no.2, 023522 (2024) doi:10.1103/PhysRevD.109.023522 [arXiv:2308.08546 [astro-ph.CO]].
  • (12) K. D. Lozanov, S. Pi, M. Sasaki, V. Takhistov and A. Wang, [arXiv:2310.03594 [astro-ph.CO]].
  • (13) G. B. Gelmini and J. Hyman, Phys. Lett. B 848, 138356 (2024) doi:10.1016/j.physletb.2023.138356 [arXiv:2307.07665 [hep-ph]].
  • (14) N. Kitajima, J. Lee, K. Murai, F. Takahashi and W. Yin, Phys. Lett. B 851, 138586 (2024) doi:10.1016/j.physletb.2024.138586 [arXiv:2306.17146 [hep-ph]].
  • (15) M. Geller, S. Ghosh, S. Lu and Y. Tsai, Phys. Rev. D 109, no.6, 063537 (2024) doi:10.1103/PhysRevD.109.063537 [arXiv:2307.03724 [hep-ph]].
  • (16) Y. Bai, T. K. Chen and M. Korwar, JHEP 12, 194 (2023) doi:10.1007/JHEP12(2023)194 [arXiv:2306.17160 [hep-ph]].
  • (17) S. Blasi, A. Mariotti, A. Rase and A. Sevrin, JHEP 11, 169 (2023) doi:10.1007/JHEP11(2023)169 [arXiv:2306.17830 [hep-ph]].
  • (18) S. Borsanyi, Z. Fodor, J. Guenther, K. H. Kampert, S. D. Katz, T. Kawanai, T. G. Kovacs, S. W. Mages, A. Pasztor and F. Pittler, et al. Nature 539, no.7627, 69-71 (2016) doi:10.1038/nature20115 [arXiv:1606.07494 [hep-lat]].
  • (19) N. S. Manton, Phys. Rev. D 28, 2019 (1983) doi:10.1103/PhysRevD.28.2019
  • (20) F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984) doi:10.1103/PhysRevD.30.2212
  • (21) D. Kharzeev and A. Zhitnitsky, Nucl. Phys. A 797, 67-79 (2007) doi:10.1016/j.nuclphysa.2007.10.001 [arXiv:0706.1026 [hep-ph]].
  • (22) D. E. Kharzeev, L. D. McLerran and H. J. Warringa, Nucl. Phys. A 803, 227-253 (2008) doi:10.1016/j.nuclphysa.2008.02.298 [arXiv:0711.0950 [hep-ph]].
  • (23) K. Fukushima, D. E. Kharzeev and H. J. Warringa, Phys. Rev. D 78, 074033 (2008) doi:10.1103/PhysRevD.78.074033 [arXiv:0808.3382 [hep-ph]].
  • (24) L. D. McLerran, E. Mottola and M. E. Shaposhnikov, Phys. Rev. D 43, 2027-2035 (1991) doi:10.1103/PhysRevD.43.2027
  • (25) G. D. Moore, Phys. Lett. B 412, 359-370 (1997) doi:10.1016/S0370-2693(97)01046-0 [arXiv:hep-ph/9705248 [hep-ph]].
  • (26) G. D. Moore and K. Rummukainen, Phys. Rev. D 61, 105008 (2000) doi:10.1103/PhysRevD.61.105008 [arXiv:hep-ph/9906259 [hep-ph]].
  • (27) D. Bodeker, G. D. Moore and K. Rummukainen, Phys. Rev. D 61, 056003 (2000) doi:10.1103/PhysRevD.61.056003 [arXiv:hep-ph/9907545 [hep-ph]].
  • (28) A. A. Andrianov, V. A. Andrianov, D. Espriu and X. Planells, Phys. Lett. B 710, 230-235 (2012) doi:10.1016/j.physletb.2012.02.072 [arXiv:1201.3485 [hep-ph]].
  • (29) A. A. Andrianov, D. Espriu and X. Planells, Eur. Phys. J. C 73, no.1, 2294 (2013) doi:10.1140/epjc/s10052-013-2294-0 [arXiv:1210.7712 [hep-ph]].
  • (30) M. Kawaguchi, S. Matsuzaki and A. Tomiya, Phys. Rev. D 103, no.5, 054034 (2021) doi:10.1103/PhysRevD.103.054034 [arXiv:2005.07003 [hep-ph]].
  • (31) C. X. Cui, J. Y. Li, S. Matsuzaki, M. Kawaguchi and A. Tomiya, Phys. Rev. D 105, no.11, 114031 (2022) doi:10.1103/PhysRevD.105.114031 [arXiv:2106.05674 [hep-ph]].
  • (32) C. X. Cui, J. Y. Li, S. Matsuzaki, M. Kawaguchi and A. Tomiya, Particles 7, no.1, 237-263 (2024) doi:10.3390/particles7010014 [arXiv:2205.12479 [hep-ph]].
  • (33) H. Leutwyler and A. V. Smilga, Phys. Rev. D 46, 5607-5632 (1992) doi:10.1103/PhysRevD.46.5607
  • (34) M. Kawaguchi, S. Matsuzaki and A. Tomiya, Phys. Lett. B 813, 136044 (2021) doi:10.1016/j.physletb.2020.136044 [arXiv:2003.11375 [hep-ph]].
  • (35) R. F. Dashen, Phys. Rev. D 3, 1879-1889 (1971) doi:10.1103/PhysRevD.3.1879
  • (36) E. Witten, Annals Phys. 128, 363 (1980) doi:10.1016/0003-4916(80)90325-5
  • (37) R. D. Pisarski, Phys. Rev. Lett. 76, 3084-3087 (1996) doi:10.1103/PhysRevLett.76.3084 [arXiv:hep-ph/9601316 [hep-ph]].
  • (38) M. Creutz, Phys. Rev. Lett. 92, 201601 (2004) doi:10.1103/PhysRevLett.92.201601 [arXiv:hep-lat/0312018 [hep-lat]].
  • (39) A. J. Mizher and E. S. Fraga, Nucl. Phys. A 831, 91-105 (2009) doi:10.1016/j.nuclphysa.2009.09.004 [arXiv:0810.5162 [hep-ph]].
  • (40) D. Boer and J. K. Boomsma, Phys. Rev. D 78, 054027 (2008) doi:10.1103/PhysRevD.78.054027 [arXiv:0806.1669 [hep-ph]].
  • (41) M. Creutz, Annals Phys. 324, 1573-1584 (2009) doi:10.1016/j.aop.2009.01.005 [arXiv:0901.0150 [hep-ph]].
  • (42) J. K. Boomsma and D. Boer, Phys. Rev. D 80, 034019 (2009) doi:10.1103/PhysRevD.80.034019 [arXiv:0905.4660 [hep-ph]].
  • (43) Y. Sakai, H. Kouno, T. Sasaki and M. Yahiro, Phys. Lett. B 705, 349-355 (2011) doi:10.1016/j.physletb.2011.10.032 [arXiv:1105.0413 [hep-ph]].
  • (44) B. Chatterjee, H. Mishra and A. Mishra, Phys. Rev. D 85, 114008 (2012) doi:10.1103/PhysRevD.85.114008 [arXiv:1111.4061 [hep-ph]].
  • (45) T. Sasaki, J. Takahashi, Y. Sakai, H. Kouno and M. Yahiro, Phys. Rev. D 85, 056009 (2012) doi:10.1103/PhysRevD.85.056009 [arXiv:1112.6086 [hep-ph]].
  • (46) T. Sasaki, H. Kouno and M. Yahiro, Phys. Rev. D 87, no.5, 056003 (2013) doi:10.1103/PhysRevD.87.056003 [arXiv:1208.0375 [hep-ph]].
  • (47) S. Aoki and M. Creutz, Phys. Rev. Lett. 112, no.14, 141603 (2014) doi:10.1103/PhysRevLett.112.141603 [arXiv:1402.1837 [hep-lat]].
  • (48) K. Mameda, Nucl. Phys. B 889, 712-726 (2014) doi:10.1016/j.nuclphysb.2014.11.002 [arXiv:1408.1189 [hep-ph]].
  • (49) J. J. M. Verbaarschot and T. Wettig, Phys. Rev. D 90, no.11, 116004 (2014) doi:10.1103/PhysRevD.90.116004 [arXiv:1407.8393 [hep-th]].
  • (50) D. Gaiotto, Z. Komargodski and N. Seiberg, JHEP 01, 110 (2018) doi:10.1007/JHEP01(2018)110 [arXiv:1708.06806 [hep-th]].
  • (51) D. Gaiotto, A. Kapustin, Z. Komargodski and N. Seiberg, JHEP 05, 091 (2017) doi:10.1007/JHEP05(2017)091 [arXiv:1703.00501 [hep-th]].
  • (52) G. ’t Hooft, NATO Sci. Ser. B 59, 135-157 (1980) doi:10.1007/978-1-4684-7571-5_9
  • (53) Y. Frishman, A. Schwimmer, T. Banks and S. Yankielowicz, Nucl. Phys. B 177, 157-171 (1981) doi:10.1016/0550-3213(81)90268-6
  • (54) K. Saikawa, Universe 3, no.2, 40 (2017) doi:10.3390/universe3020040 [arXiv:1703.02576 [hep-ph]].
  • (55) T. Hiramatsu, M. Kawasaki and K. Saikawa, JCAP 05, 032 (2010) doi:10.1088/1475-7516/2010/05/032 [arXiv:1002.1555 [astro-ph.CO]].
  • (56) T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, JCAP 01, 001 (2013) doi:10.1088/1475-7516/2013/01/001 [arXiv:1207.3166 [hep-ph]].
  • (57) T. Hiramatsu, M. Kawasaki and K. Saikawa, JCAP 02, 031 (2014) doi:10.1088/1475-7516/2014/02/031 [arXiv:1309.5001 [astro-ph.CO]].
  • (58) R. L. Workman et al. [Particle Data Group], PTEP 2022, 083C01 (2022) doi:10.1093/ptep/ptac097
  • (59) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 06, 088 (2009) doi:10.1088/1126-6708/2009/06/088 [arXiv:0903.4155 [hep-lat]].
  • (60) S. Borsanyi et al. [Wuppertal-Budapest], J. Phys. Conf. Ser. 316, 012020 (2011) doi:10.1088/1742-6596/316/1/012020 [arXiv:1109.5032 [hep-lat]].
  • (61) H. T. Ding, F. Karsch and S. Mukherjee, Int. J. Mod. Phys. E 24, no.10, 1530007 (2015) doi:10.1142/S0218301315300076 [arXiv:1504.05274 [hep-lat]].
  • (62) A. Bazavov et al. [HotQCD], Phys. Lett. B 795, 15-21 (2019) doi:10.1016/j.physletb.2019.05.013 [arXiv:1812.08235 [hep-lat]].
  • (63) H. T. Ding, Nucl. Phys. A 1005, 121940 (2021) doi:10.1016/j.nuclphysa.2020.121940 [arXiv:2002.11957 [hep-lat]].
  • (64) K. Fukushima, Phys. Lett. B 591, 277-284 (2004) doi:10.1016/j.physletb.2004.04.027 [arXiv:hep-ph/0310121 [hep-ph]].
  • (65) K. Fukushima and V. Skokov, Prog. Part. Nucl. Phys. 96, 154-199 (2017) doi:10.1016/j.ppnp.2017.05.002 [arXiv:1705.00718 [hep-ph]].
  • (66) M. Ruggieri, M. N. Chernodub and Z. Y. Lu, Phys. Rev. D 102, no.1, 014031 (2020) doi:10.1103/PhysRevD.102.014031 [arXiv:2004.09393 [hep-ph]].
  • (67) S. Chen, K. Fukushima, H. Nishimura and Y. Tanizaki, Phys. Rev. D 102, no.3, 034020 (2020) doi:10.1103/PhysRevD.102.034020 [arXiv:2006.01487 [hep-th]].
  • (68) M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 44, 1422-1424 (1970) doi:10.1143/PTP.44.1422
  • (69) M. Kobayashi, H. Kondo and T. Maskawa, Prog. Theor. Phys. 45, 1955-1959 (1971) doi:10.1143/PTP.45.1955
  • (70) G. ’t Hooft, Phys. Rev. Lett. 37, 8-11 (1976) doi:10.1103/PhysRevLett.37.8
  • (71) G. ’t Hooft, Phys. Rev. D 14, 3432-3450 (1976) [erratum: Phys. Rev. D 18, 2199 (1978)] doi:10.1103/PhysRevD.14.3432
  • (72) G. Fejos and A. Hosaka, Phys. Rev. D 94, no.3, 036005 (2016) doi:10.1103/PhysRevD.94.036005 [arXiv:1604.05982 [hep-ph]].
  • (73) G. Fejos and A. Patkos, Phys. Rev. D 109, no.3, 036035 (2024) doi:10.1103/PhysRevD.109.036035 [arXiv:2311.02186 [hep-ph]].
  • (74) G. Fejös and A. Patkos, Phys. Rev. D 105, no.9, 096007 (2022) doi:10.1103/PhysRevD.105.096007 [arXiv:2112.14903 [hep-ph]].
  • (75) L. Del Debbio, H. Panagopoulos and E. Vicari, JHEP 08, 044 (2002) doi:10.1088/1126-6708/2002/08/044 [arXiv:hep-th/0204125 [hep-th]].
  • (76) M. D’Elia, Nucl. Phys. B 661, 139-152 (2003) doi:10.1016/S0550-3213(03)00311-0 [arXiv:hep-lat/0302007 [hep-lat]].
  • (77) L. Del Debbio, G. M. Manca, H. Panagopoulos, A. Skouroupathis and E. Vicari, PoS LAT2006, 045 (2006) doi:10.22323/1.032.0045 [arXiv:hep-th/0610100 [hep-th]].
  • (78) L. Giusti, S. Petrarca and B. Taglienti, Phys. Rev. D 76, 094510 (2007) doi:10.1103/PhysRevD.76.094510 [arXiv:0705.2352 [hep-th]].
  • (79) T. Izubuchi, S. Aoki, K. Hashimoto, Y. Nakamura, T. Sekido and G. Schierholz, PoS LATTICE2007, 106 (2007) doi:10.22323/1.042.0106 [arXiv:0802.1470 [hep-lat]].
  • (80) E. Vicari and H. Panagopoulos, Phys. Rept. 470, 93-150 (2009) doi:10.1016/j.physrep.2008.10.001 [arXiv:0803.1593 [hep-th]].
  • (81) M. D’Elia and F. Negro, Phys. Rev. D 88, no.3, 034503 (2013) doi:10.1103/PhysRevD.88.034503 [arXiv:1306.2919 [hep-lat]].
  • (82) M. Hirasawa, K. Hatakeyama, M. Honda, A. Matsumoto, J. Nishimura and A. Yosprakob, PoS LATTICE2023, 193 (2024) doi:10.22323/1.453.0193 [arXiv:2401.05726 [hep-lat]].
  • (83) R. Kitano, R. Matsudo, N. Yamada and M. Yamazaki, Phys. Lett. B 822, 136657 (2021) doi:10.1016/j.physletb.2021.136657 [arXiv:2102.08784 [hep-lat]].
  • (84) T. Byrnes, P. Sriganesh, R. J. Bursill and C. J. Hamer, Nucl. Phys. B Proc. Suppl. 109, 202-206 (2002) doi:10.1016/S0920-5632(02)01416-0 [arXiv:hep-lat/0201007 [hep-lat]].
  • (85) L. Funcke, K. Jansen and S. Kühn, Phys. Rev. D 101, no.5, 054507 (2020) doi:10.1103/PhysRevD.101.054507 [arXiv:1908.00551 [hep-lat]].
  • (86) Y. Kuramashi and Y. Yoshimura, JHEP 04, 089 (2020) doi:10.1007/JHEP04(2020)089 [arXiv:1911.06480 [hep-lat]].
  • (87) B. Chakraborty, M. Honda, T. Izubuchi, Y. Kikuchi and A. Tomiya, Phys. Rev. D 105, no.9, 094503 (2022) doi:10.1103/PhysRevD.105.094503 [arXiv:2001.00485 [hep-lat]].
  • (88) M. Honda, E. Itou, Y. Kikuchi, L. Nagano and T. Okuda, Phys. Rev. D 105, no.1, 014504 (2022) doi:10.1103/PhysRevD.105.014504 [arXiv:2105.03276 [hep-lat]].
  • (89) M. Honda, doi:10.1142/9789811261633_0003
  • (90) A. Gómez Nicola and J. Ruiz de Elvira, JHEP 03, 186 (2016) doi:10.1007/JHEP03(2016)186 [arXiv:1602.01476 [hep-ph]].
  • (91) S. Roessner, C. Ratti and W. Weise, Phys. Rev. D 75, 034007 (2007) doi:10.1103/PhysRevD.75.034007 [arXiv:hep-ph/0609281 [hep-ph]].
  • (92) C. Ratti, M. A. Thaler and W. Weise, Phys. Rev. D 73, 014019 (2006) doi:10.1103/PhysRevD.73.014019 [arXiv:hep-ph/0506234 [hep-ph]].