The QCD topological charge and its thermal dependence: the role of the $η'$
Angel Gómez Nicola Jacobo Ruiz De Elvira Andrea Vioque-Rodríguez
hep-ph






gomez@ucm.es
Departamento de Física
Teórica and IPARCOS. Univ. Complutense. 28040 Madrid. Spain

elvira@itp.unibe.ch
Albert Einstein Center for Fundamental Physics, Institute for Theoretical Physics,
University of Bern, Sidlerstrasse 5, CH–3012 Bern, Switzerland

avioque@ucm.es
Departamento de Física
Teórica and IPARCOS. Univ. Complutense. 28040 Madrid. Spain




We analyze the contribution of the η'(958) meson in the first two non-trivial moments of the QCD topological charge distribution, namely, the topological susceptibility and the fourth-order cumulant of the vacuum energy density.
We perform our study within U(3) Chiral Perturbation Theory up to next-to-next-to-leading order in the combined chiral and large-N_c expansion. 
We also describe the temperature dependence of these two quantities and compare them with previous analyses in the literature. 
In particular, we discuss the validity of the thermal scaling of the topological susceptibility with the quark condensate, which is intimately connected with a Ward Identity relating both quantities.  
We also consider isospin breaking corrections from the vacuum misalignment at leading order in the U(3) framework. 


 11.30.Rd, 
 11.10.Wx, 
  12.39.Fe, 
  25.75.Nq. 
 12.38.Gc. 
 14.80.Va	
 11.30.Er	
 
 
The QCD topological charge and its thermal dependence: the role of the η'
    A. Vioque-Rodríguez
    
=========================================================================





§ INTRODUCTION AND MOTIVATION



The rich topological structure of the QCD  vacuum is encoded in the θ-angle dependence of the vacuum energy density 

    e_vac=-1/β Vln Z(θ)  with   Z_QCD(θ)=∫[ G][ψ̅][ψ] e^S_QCD(θ)

the QCD generating functional in a θ-vacuum, β=1/T is the inverse of the temperature, V the volume of the system and the QCD action is written as

    S_QCD(θ)=∫^4 x [ L_QCD-θ(x)ω(x)]

with L_QCD the QCD Lagrangian at θ=0 and 

    ω(x)   =g^2/32π^2_c G_μνG̃^μν

the winding number topological charge density, responsible for the U_A(1) anomaly.

The expansion of the vacuum energy density around θ=0 can be expressed as 

    ϵ_vac (θ)=∑_n=1^∞c_2n/(2n)!θ^2n

 with c_2=χ_top the topological susceptibility and c_4 the fourth-order cumulant, which in Euclidean space-time read

    χ_top   =∫_T  x  ⟨𝒯ω(x) ω(0) ⟩
    
    c_4   =-∫_T  x  y  z[ ⟨𝒯ω(x) ω(y)ω(z) ω(0)⟩ -3⟨𝒯ω(x) ω(0)⟩^2]
 
with ∫_T  x=∫_0^β dτ∫ d^3x⃗. 

The topological susceptibility is meant to be  connected with the η' mass through the U_A(1) anomaly. More specifically, in the quenched approximations where quarks loops are absent, 
χ_top^quenched is related to the mass of the η' (or rather its singlet part η_0) for N_f massless quarks as <cit.>

    χ_top^quenched=1/2N_fF^2 M_0^2=F^2[M^m_q=0_η_0]^2/2N_f,

where M_0 denotes the anomalous contribution to the η'  mass and F is the pion decay constant in the chiral m_u=m_d=m_s=0 limit. 
The quenched approximation is formally valid in the N_c→∞ limit, where meson-loop contributions are suppressed and the mass of the η_0 becomes of the same order than the other members of the Nambu-Goldstone boson (NGB) meson octet, since the U(1)_A anomaly scales with 1/N_c. 

However, even when the fluctuation of the winding number is directly linked to the U(1)_A anomaly, meson-loop corrections are indeed very relevant and one of the reasons why the QCD topological charge can be analyzed using low-energy effective field theories. 
This  can be seen by looking at the leading-order (LO) low-energy chiral prediction for the topological susceptibility, 
a well-known result in Chiral Perturbation Theory (ChPT) that for three light flavors N_f=2+1 reads <cit.>

    χ_top^SU(3),LO=Σ[1/m_u+1/m_d+1/m_s]^-1≡Σm̅,

where Σ=B_0F^2=-  is the single-flavor quark condensate in the chiral limit and B_0=M_0π^±^2/(m_u+m_d), with M_0π^± the mass of the charged pions. 
The result in (<ref>) can be easily extended to a larger number of flavors <cit.>, while the two-flavor result can be recovered by taking the m_s>>m_u,d limit.  

One of the main consequences of (<ref>) is that the topological susceptibility is linearly proportional to the quark mass and hence, unlike the pure gluonic quenched result (<ref>), it vanishes in the chiral limit. Thus, meson loop corrections generate terms that cancel the contribution in (<ref>). This fact is actually crucial to understand why the QCD topological charge can be reliable described just by including the lightest degrees of freedom.

In addition, (<ref>) also shows that at LO in the chiral expansion the topological susceptibility is proportional to the quark condensate. 
This property has been used in several lattice analyses to extract the value of  in the chiral limit from χ_top, both for N_f=2 <cit.> and N_f=2+1 <cit.>. 
In these analyses, χ_top is determined at different quark-mass values and then fitted using ChPT predictions to extract the quark condensate. 
This relation has been emphasized in <cit.> to study the strange mass paramagnetic suppression of the three-flavor condensate. 
Recent values of the quark condensate obtained with this method can be found in <cit.>, while direct lattice measurements of χ_top are provided in <cit.> for N_f=2+1, in <cit.> for N_f=2 and in <cit.> within the framework of two fermion families.  
Therefore, obtaining explicit analytic expressions for χ_top and studying its quark mass dependence within the effective Lagrangian framework under different approximations is of the utmost importance from the point of view of lattice analyses. 

Nevertheless, higher order corrections change the simple linear dependence of χ_top with the quark condensate given in (<ref>). In this sense, an important result that hints towards more complicated dependencies is the existence of the following family of Ward Identities (WI) connecting quark condensates with the topological and pseudoscalar susceptibilities in the isospin limit <cit.>
[We acknowledge a misprint in eq (31) of <cit.> where  should read .]:

    χ_top   =-1/4[m_q+m_q^2 χ_P^ll] 
       =-[m_s+m_s^2χ_P^ss]

 where =⟨u̅ u + d̅ d⟩, m_q=m_u=m_d and χ_P^ll, χ_P^ss are the pseudoscalar susceptibilities corresponding to the quark bilinear operators η_l=i(u̅γ_5 u+d̅γ_5 d) and η_s=i s̅ s, respectively. 

 The identities (<ref>) and (<ref>) have been verified in U(3) ChPT up to next-to-next-leading order (NNLO) in <cit.>. Furthermore, the identity (<ref>) has been used in lattice works to determine χ_top indirectly at finite temperature <cit.> and in <cit.> to justify that U(1)_A restoration approaches the O(4) phase transition when the chiral symmetry is exactly restored. 

Now, note that for N_f=2 in the isospin limit, the first term in the r.h.s. of (<ref>) corresponds to the LO ChPT expression (<ref>) when m_s>>m_q. Therefore, the term proportional to χ_P^ll necessarily includes  higher order corrections in the chiral series and/or terms suppressed as m_q/m_s for N_f=2+1 flavors.  An immediate conclusion is that χ_top is not necessarily proportional to the light quark condensate to all orders. In fact, as we will see here in detail, the dependence with temperature of the two terms in the r.h.s. of (<ref>) is completely different. 
Namely, the condensate term drops with T signaling chiral restoration while the χ_P^ll term shows a much smoother behavior.  

Higher order corrections to the topological susceptibility within the chiral Lagrangian framework have been obtained in <cit.> to next-to-leading order (NLO) (one-loop) in  SU(N_ f) ChPT. 
This result has been incorporated in the lattice analysis of <cit.> to extract the quark condensate and in <cit.> to provide a numerical estimate of χ_top in terms of the low-energy constants (LECs) of the (p^4) effective Lagrangian. Higher order isospin-breaking corrections as well as an analysis of χ_top within the  so called resummed ChPT can be found in <cit.>, whereas in the recent work <cit.> NNLO and electromagnetic corrections in SU(2) ChPT have been computed. 

Another important motivation for the study of the topological susceptibility is its relation with the Peccei-Quinn axion <cit.> and hence with various cosmological and astrophysical implications. The axion mass is directly proportional to χ_top, which allows for  numerical estimates based on chiral Lagrangians <cit.>. In addition, the fourth-order self-coupling of the axion field can be obtained from the fourth-order cumulant c_4 of the ϵ_vac(θ) expansion in (<ref>), which has been thoroughly analyzed  within SU(N_ f) ChPT up to NLO <cit.> and computed in the lattice <cit.>.  

In addition, the large-N_c behavior of the topological susceptibility and the fourth-order cumulant have been recently analyzed in <cit.> and compared with the lattice results provided at different large-N_c values in <cit.>.

The thermal dependence of the vacuum energy density  ϵ(θ) is important for several reasons:  as stated above χ_top(T) plays an important role in the relation between chiral and U(1)_A restoration. In addition, the connection of χ_top with the light- and strange-quark condensate as a function of temperature is relevant for lattice analyses and for the understanding of the temperature dependence of the different contributions in the WI (<ref>)-(<ref>).  Thermal corrections to the axion potential and its mass are also of importance in the cosmological and astrophysical context <cit.>.  
A ChPT analysis for χ_top(T) for N_f=2 has been performed in <cit.>. In that paper, the fact that for N_f=2 the one-loop corrections to χ_top can be  encoded in the physical pion mass and  decay constant is used to establish the scaling χ_top(T)/χ_top(0)= (T)/ (0), valid  at that order. Actually, that scaling law is nothing but the first term in the r.h.s. of (<ref>), which opens the question of how relevant are the additional corrections provided by the second term in that identity. As we have already mentioned, we will discuss that particular aspect in detail in the present  work.  



As for lattice results at finite temperature, direct measurements of the topological susceptibility have typically large errors. Results can be found  e.g. in <cit.>. As commented above, indirect measurements can be obtained precisely through (<ref>). Higher order cumulants at finite temperature in the lattice have been analyzed in <cit.>.

With the above motivation in mind, we will carry out here a ChPT-based analysis of the topological susceptibility and the fourth-order cumulant, concentrating in particular in the following aspects: 


  * We will provide a NNLO calculation of the topological susceptibility and the fourth-order cumulant within the formalism of U(3) ChPT, in which the singlet η_0 is incorporated  as a ninth pseudo-Goldstone boson within the large-N_c framework  <cit.>. This approach allows us to study η' meson effects and to assess its contribution. For instance, the inclusion of η' will help to understand the relevance of meson-loop corrections to the quenched result (<ref>) and their role on the vanishing of χ_top in the chiral limit, as already noticed in <cit.>.
 Our calculation will also allow us to estimate numerically the effect of the additional U(3) corrections in terms of the LECs involved. As mentioned above, many lattice analyses of the quark condensate rely on the chiral expansion of χ_top. In that respect, studying the influence of an additional heavier degree of freedom is important. The η' case is especially significant due to its direct connection with the axial anomaly and its relevant role in the Ward identities described above. In addition, our U(3) analysis will provide a natural way to establish contact with recent large-N_c analysis of the topological susceptibility and the fourth-order cumulant. Our study will also have the advantage of including explicitly the dependence of those quantities on the η-η' mixing angle. 


  * We will also calculate the leading isospin-breaking m_u≠ m_d corrections within the U(3) formalism, thus extending previous SU(3) works. The importance of isospin breaking corrections to the topological susceptibility and the fourth-order cumulant stems from the vacuum misalignment induced by the combined m_u≠ m_d and θ≠ 0 effects <cit.>. It implies corrections proportional to (m_u-m_d)/(m_u+m_d), hence much larger than the typical isospin-breaking correction in other quantities, such as quark condensates or the π^0η mixing <cit.>, which are proportional to (m_u-m_d)/m_s.
Recent estimates within the SU(2) formalism show that these isospin contributions give rise to around a 4% correction to χ_top^1/4. 
As a natural extension to those analysis, we will include isospin breaking in the LO U(3) correction to ϵ_vac(θ) to estimate its numerical effect.


  * We will extend our analysis to finite temperature. In this way, we will account for corrections both within the SU(3) and U(3) formalisms for the different contributions in the Ward identities (<ref>) and (<ref>). This will also allow us to test the robustness of the N_f=2 scaling  performed in <cit.> when corrections from m_q /m_s and η,η' loops are properly incorporated. As commented before, the evolution of  quark condensates and susceptibilities towards chiral restoration makes it interesting to clarify their relation with the topological susceptibility as the temperature grows, within the context of chiral and (1)_A restoration.  In addition, we recall that the finite temperature dependence of the η' mass has been analyzed within the U(3) ChPT formalism in <cit.>,  in fermion models <cit.> and in the lattice <cit.>, confirming the U(1)_A restoring behavior, as in the recent analysis of the η-η' mixing angle <cit.>.



For that purpose, the paper is organized as follows. In section <ref> we will discuss the U(3) ChPT calculation of the topological susceptibility, providing explicit analytic expressions up to NNLO. 
We will also analyze its various limits of interest. In section <ref> we will extend the calculation to the fourth order cumulant. In Section <ref> we will provide numerical estimates and compare to previous approaches. The isospin breaking corrections to the U(3) susceptibility and the fourth-order cumulant will be calculated in section <ref>. In section <ref> we will study in detail their finite temperature dependence and  their connection with chiral and U(1)_A restoration. Some of the explicit analytic U(3) expressions will be collected in Appendix <ref>.

 


§ THE  TOPOLOGICAL SUSCEPTIBILITY TO NNLO IN U(3) CHIRAL PERTURBATION THEORY





Within U(3) ChPT one follows a similar approach as in standard ChPT <cit.>. The chiral power counting in terms of momenta and quark masses is used to construct the most general effective Lagrangian for SU(N_ f) pseudo-Goldstone bosons up to a given order, which ensures renormalizability order by order in the expansion. 
Nevertheless, in the U(3) formalism the singlet η_0 is also included as the ninth pseudo-Goldstone boson. Given the large η' mass value, this can be done consistently only in the large-N_c framework, 
since the winding number charge density ω(θ) defined in (<ref>), responsible for the anomalous contribution of the η' mass, is suppressed within the large-N_c counting <cit.>. 

Thus, the expansion is performed in terms of a small parameter δ such that {M^2, p^2, T^2, m_q,m_s,1/N_c}=(δ), where M and p denote typical meson masses and momenta. 
In this counting, the tree-level pion decay constant F^2=(N_c)=(1/δ), which hence suppresses loop diagrams. 
The counting of the different LECs according to their (N_c) trace structure is given in detail in <cit.>. 

Following the same steps as in <cit.>, the topological susceptibility can be calculated by taking functional derivatives with respect to the vacuum angle θ. 
Thus, taking into account the θ-vacuum coupling in the QCD action defined in (<ref>) and (<ref>), 
one can derive expectation values or thermal correlators involving the winding number density. In this way, the topological susceptibility reads
 
    χ_top   =∫_T  x .δ/δθ (x)δ/δθ (0)log Z(θ)|_θ=0
       =∫_T  x {⟨δ L_ eff (x)/δθ (x)δ L_ eff (0)/δθ (0)⟩_θ=0 +⟨δ/δθ (x)δ/δθ (0) L_ eff (x)  ⟩_θ=0δ^4(x)}

 where ⟨·⟩ denotes Euclidean vacuum expectation values for T=0 or thermal correlators for T≠ 0 and where we have used that ⟨ω(0)⟩=0. 
 
 In the effective Lagrangian L_ eff, θ(x) appears through the operator 

    X(x)=log[ U(x)] + iθ(x),

with U=exp(iΦ/F)=exp(iλ^aϕ_a/F) the NGB matrix field including the singlet contribution (i.e., U≠ 1), ϕ_i the Goldstone fields,  λ^a=1,… 8 the Gell-Mann matrix and λ^0=√(2/3) <cit.>,

    Φ = ( [  π^0+1/√(3)η_8+√(2/3)η_0                  √(2)π^+                  √(2)K^+;                  √(2)π^- -π^0+1/√(3)η_8+√(2/3)η_0                  √(2)K^0;                  √(2)K^-                 √(2)K̅^0     -2/√(3)η_8+√(2/3)η_0 ]) .

Nevertheless, due to the η-η' mixing the flavor eigenstates η_8 and η_0 are not mass eigenstates even at LO in the chiral Lagrangian. 
Thus, we use the angle θ̂ to describe their mixing at LO

    η_8=  η + η'  η_0= - η + η',

where =cosθ̂ and =sinθ̂. 

The combination (<ref>) ensures that L_ eff will be invariant under local U(N_f)_L× U(N_f)_R transformations. 
The effective Lagrangians containing the X field entering in χ_top and the fourth-order cumulant up to the order needed for our purposes here read
 
    (δ^0):        L_ eff^(0)=F^2/12 M_0^2  X^2 
    (δ):        L_ eff^(1)=-F^2/12Λ_2 X (U^†χ-χ^† U)
    (δ^2):        L_ eff^(2)= F^2/4[v_4^(0)X^4+v_2^(2) X^2 (U^†χ+χ^† U)] +L_25 X ( U^†χ U^†χ-χ^† U χ^† U)

where  M_0^2 is  the contribution to the tree-level mass of the singlet η_0 in the chiral limit, i.e., its anomalous contribution as given in (<ref>),  and χ=2B_0 M with M=(m_u,m_d,m_s) the quark mass matrix. The constants Λ_2=(δ), v_2^(2)=(δ^2), v_4^(0)=(δ^3)  and L_25=(δ^0) are LECs associated to the η-η' mixing. As in the standard SU(3) formalism, the LECs  are renormalized to absorb divergences coming from the loops. 

The first non-vanishing contribution to the topological susceptibility is (δ^0). On the one hand, the second term in the r.h.s. of (<ref>) gives rise at this order to a constant term, 
which is nothing but the contribution (<ref>) with [M^m_q=0_η_0]^2=M_0^2 and N_f=3. 
On the other hand, the first term in the r.h.s of (<ref>)  gives rise to terms of the type  M_0^4∫ x⟨η^(')(x)η^(')(0)⟩∼ M_0^4/M_η^(')^2, 
i.e., tree level two-point functions at p=0,  where η^(') denotes generically η or η' fields.  Gathering the two types of contribution yields in the isospin limit

    χ_top^U(3),  LO,  IL=F^2M_0^2/6[1-M_0^2/M_0η^2^2-M_0^2/M_0η'^2^2  ]

where M_0η and M_0η' are respectively the η and η' masses at tree level. 
They depend on the quark (or meson) masses and on M_0. Their explicit expressions in the isospin limit m_u=m_d=m_q can be found  e.g. in <cit.>. 
Thus, in terms of quark masses, (<ref>) can be recast as:

    χ_top^U(3),  LO,  IL=F^2M_0^2 B_0 m_q m_s/(2m_s+ m_q)M_0^2+6B_0 m_q m_s =Σ M_0^2 m̅_IL/M_0^2+6B_0 m̅_IL=Σm̂_IL

where m̅_IL denotes the isospin limit value of  m̅  in (<ref>), i.e.,  
 
    m̅_IL=m_q m_s/2m_s+m_q

   and
 
    m̂_IL=  M_0^2 m̅_IL/M_0^2+6B_0 m̅_IL

 
 
 Note that in this work we use the symbol m̂, as in (<ref>), with a different meaning than the light-quark mass, which we denote here by  m_q.  The result in (<ref>) and (<ref>) is the extension to N_f=2+1 flavors of the  result given in <cit.> for N_f degenerated flavors in the large-N_c limit. 
Note that our U(3) LO result can be obtained from the SU(3) one by replacing m̅ with m̂, which also holds when including isospin breaking corrections at LO, as we will see in section <ref>. 
In the M_0→∞ limit one recovers the SU(3) expression in (<ref>) for m_u=m_d=m_q, since m̂_IL→m̅_IL in this limit.
Furthermore, χ_top^U(3),LO vanishes in the chiral limit, as it also does in the SU(3) case.  
Actually, in this limit the second and third terms in the r.h.s. of (<ref>) (coming from η,η' propagators) cancel the first term, as one can check from the chiral limit behavior of s_θ,M_0η,M_0η'.
 
Our present U(3) calculation has also the advantage that one can formally recover the quenched result in (<ref>) by taking the limit M_0≪m̅. 
This limit can also be achieved by re-expanding (<ref>) in the 1/N_c expansion. Taking into account the 1/N_c scaling of the different constants involved, namely F^2=(N_c), M_0^2=(1/N_c) and so on, the 1/N_c LO contribution to (<ref>) gives
 
    χ_top^U(3),IL=F^2 M_0^2/6   + (1/N_c)=  F^2/6(M_0η'^2+M_0η^2-2M_0K^2)  + (1/N_c)

 which coincides with the result in (<ref>) for N_f=3.
 The last equality in (<ref>) reproduces the result given in <cit.>, where M_0i denote tree-level meson masses. 
 
The NLO ((δ)) and NNLO ((δ^2)) results require including the higher-order effective Lagrangians (<ref>)-(<ref>) in (<ref>). They involve one-loop corrections to the η- and η'-meson propagators at zero momentum, which include the LECs L_6, L_7, L_25, C_19, C_31, Λ_2 and v_2^(2) <cit.>. The renormalization of these LECs and the constant B_0 <cit.> allows one to absorb all one-loop divergences, rendering the result finite and independent of the renormalization scale μ. In addition, the calculation of χ_top, χ_P^ll, χ_P^ss and  up to NNLO in the U(3) expansion allow one to verify the Ward Identities (<ref>)-(<ref>).  The explicit results for the topological susceptibility at NLO and NNLO in U(3) ChPT are given explicitly in Appendix <ref> in the isospin limit. Note that the NLO and NNLO U(3) results do not correspond to simply perform the replacement m̅→m̂ as it happened at LO.  

In order to compare with the SU(3) calculation in ChPT, as given for instance in <cit.>, we recall that the NLO order ChPT result is included distributed among the NLO and NNLO U(3) outcome. On the one hand, we have checked that the M_0→∞ limit of χ_top^U(3),NLO in (<ref>) yields the contribution proportional to the renormalized LEC L_8^r in <cit.> while  χ_top^U(3),NNLO in (<ref>) for M_0→∞ provides the rest of the SU(3) contributions, proportional to L_6^r, L_7 and log(M^2_π,K,η/μ^2).  On the other hand, the surviving term proportional to log(M^2_η'/μ^2) is absorbed in  L_6^r. 

Finally, we remark that our present U(3)  formalism allows us to study  systematically the large N_c corrections in (<ref>), corresponding to the Witten-Veneziano result. Performing the 1/N_c expansion on the different orders in the U(3) ChPT expansion we obtain the (1/N_c) correction to (<ref>), namely

    χ_top^U(3),IL   =   F^2 M_0^2/6{1+1/3 M_0^2 (1/M_0π^2-2 M_0K^2-2/M_0π^2)+16  M_0^2 L_8^r/F^2-2 Λ _2 + 2 (2 M_0K^2+M_0π^2)/3 F^4 M_0^2[24  C_19 F^2 M_0^4 ..
       +   .. 16  C_31 F^2 M_0^4+9  F^4 v_2^(2)+16  F^2 Λ _2 M_0^2 L_8^r+24  F^2 M_0^2 L_25^r-128  M_0^4 (L_8^r)^2]} + (1/N_c^2)+ (δ^3)









The above result is consistent with the large-N_c  scaling analysis of  the topological susceptibility provided in <cit.>. 



§ FOURTH-ORDER CUMULANT TO NNLO IN U(3) CHPT



The fourth-order cumulant is defined in (<ref>) and involves the difference between the four- and two-point function square of the winding number density. 
Similarly to χ_top, it can be computed by taking functional derivatives with respect to the θ-vacuum angle. 
Taking into account again the θ-angle coupling in the QCD action as defined in (<ref>) and (<ref>), the fourth-order cumulant can be written as
 
    c_4=   -∫_T  x  y  z [.δ/δθ (x)δ/δθ (y)δ/δθ (z)δ/δθ (0)log Z(θ)|_θ=0-3(.δ/δθ (x)δ/δθ (0)log Z(θ)|_θ=0)^2]=
    
    =   -[⟨δ^4  L_ eff (0)/δθ (0)^4⟩_θ=0+∫_T  x{4⟨δ^3  L_ eff (x)/δθ (x)^3δ L_ eff (0)/δθ (0)⟩_θ=0+3⟨δ^2  L_ eff (x)/δθ (x)^2δ^2 L_ eff (0)/δθ (0)^2⟩_θ=0-3⟨δ^2  L_ eff (x)/θ (x)^2⟩_θ=0^2}
         +6∫_T  x  y {⟨δ^2 L_ eff (x)/δθ (x)^2δ L_ eff (y)/δθ (y)δ L_ eff (0)/δθ (0)⟩_θ=0-⟨δ^2  L_ eff (x)/δθ (x)^2⟩_θ=0⟨δ L_ eff (y)/δθ (y)δ L_ eff (0)/δθ (0)⟩_θ=0}
         +∫_T  x  y  z{⟨δ L_ eff (x)/δθ (x)δ L_ eff (y)/δθ (y)δ L_ eff (z)/δθ (z)δ L_ eff (0)/δθ (0)⟩_θ=0-3⟨δ L_ eff (x)/δθ (x)δ L_ eff (0)/δθ (0)⟩_θ=0^2}]

where we have used once more that  ⟨ω(0)⟩=0. The last terms in all rows are associated to the square of χ_top and hence they provide disconnected contributions, i.e., terms proportional to the Euclidean four-dimensional volume that should cancel out exactly with the disconnected contributions coming from the four-point function.
Thus, the calculation of the fourth-order cumulant involves five connected contributions associated to all possible combinations of a total even number of derivatives with respect to θ  applied to the effective lagrangian.


The LO topologies for each of these five connected contributions are depicted in Fig. <ref>[Note that we only refer about the LO topology for each connected contribution in (<ref>) independently of their corresponding counting in the U(3) expansion.].

Following the notation given in <cit.>, the θ-induced vertices are classified in terms of their number of derivatives with respect to the θ-vacuum angle.
Meson lines are always coupled to one θ-induced source coming from a single, second or third derivative with respect to the vacuum angle [The fourth derivative is just a contact term without any meson line up to NNLO in the δ expansion.], which in the isospin limit involve η and η' meson propagators at zero momentum. The four different θ-induced vertices contributing to the topologies in Fig. <ref> read

    .δ L_ eff (x)/δθ (x)|_θ=0=   -F/√(6)[M_0^2(  η'(x)- η(x))+2/3Λ_2B_0((2 m_q+m_s)(  η'(x)- η(x))
         +√(2)(m_q-m_s)(  η'(x)+ η(x)))+12B_0v_2^(2)(2m_q+m_s)(  η'(x)- η(x))]
         +16√(2/3)B_0^2/FL_25[(2m_q^2+m_s^2)(  η'(x)- η(x))+√(2)(m_q^2-m_s^2)(  η'(x)- η(x))]
         +1/9√(3)FB_0Λ_2[1/3((m_q - 4 m_s)^3 -3 √(2) (m_q + 2 m_s) ^2  + 6 (m_q - m_s) ^2 - √(2) (2 m_q + m_s) ^3))η(x)^3
         +(√(2) (m_q + 2 m_s)^3  - 3  m_q^2 -3 √(2)  m_s ^2 + 2 (m_q - m_s) ^3))η(x)^2η'(x)
         +(2 (m_q - m_s)^3  + 3  m_s^2 -3 √(2)  m_q ^2 - √(2) (m_q + 2 m_s) ^3))η(x)η'(x)^2 
         +1/3(√(2) (2 m_q + m_s) ^3 +6 (m_q - m_s) ^2  + 3 √(2) (m_q + 2 m_s) ^2 + (m_q - 4 m_s)^3))η'(x)^3]
         +2/F√(2/3)B_0v_2^(2)(  η'(x)- η(x))[((m_q + 2 m_s)^2- 2√(2)(m_q - m_s)+(2 m_q + m_s)^2))η(x)^2
          +2(m_q - m_s)(√(2)^2--√(2)^2))η(x)η'(x)+((2m_q + m_s)^2+2√(2)(m_q - m_s)
         +(m_q + 2m_s)^2))η'(x)^2]+6√(6)/F v_4^(4)(  η'(x)- η(x))^3+⋯+(δ^9/2)
    .δ^2  L_ eff (x)/δθ (x)^2|_θ=0=   -F^2[M_0^2/6+ (2m_q +m_s) B_0 v_2^(2)]+2/3v_2^(2)B_0[((m_q + 2 m_s)^2- 2√(2)(m_q - m_s)+(2 m_q + m_s)^2))η(x)^2
          +2(m_q - m_s)(√(2)^2--√(2)^2))η(x)η'(x)+((2m_q + m_s)^2+2√(2)(m_q - m_s)
         +(m_q + 2m_s)^2))η'(x)^2]+18v_4^(0)(η'(x)-η(x))^2+⋯+(δ^4),
    .δ^3  L_ eff (x)/δθ (x)^3|_θ=0=   6√(6) F v_4^(4)(  η'(x)- η(x))+(δ^7/2),
    .δ^4  L_ eff (x)/δθ (x)^4|_θ=0=   6 F^2 v_4^(0)+(δ^3),
 
where the ellipses denote further terms involving π's and K's fields, which do not contribute to connected contributions. 
The diagram displayed in Fig. <ref> (a) is generated from the fourth-order derivative induced vertex in (<ref>). 
It implies a contact term proportional to F^2v_4^(0) and hence it only contributes at (δ^2). 
The diagram (b) in Fig. <ref> comes from the third-order derivative induced vertex (<ref>) and the LO contribution in the δ expansion to the single-derivative vertex, i.e., first term in (<ref>).
It involves a single η^(') propagator at vanishing momentum and hence terms proportional to F^2M_0^2 v_4^(0)/M^2_0η^('), where η^(') stands for a η or η' field, that contribute at (δ^2).
Diagram (c) is coming from the product of two second-order derivative vertices at NLO in the δ expansion, hence involving two η^(') propagators. Nevertheless, it contributes only at (δ^4) and it will not enter in the NNLO U(3) calculation for the fourth-order cumulant. The topology shown in Fig. <ref> (d) is produced from the NLO U(3) contribution of the second derivative vertex in (<ref>) and two single-derivative vertices at LO in the δ expansion. 
It involves two η^(') propagators and terms proportional to M_0^4 v_2^(2)/M^4_0η^(') and M_0^4 v_4^(0)/M^4_0η^('). Thus, it contributes at (δ^2). 
Finally, diagrams (e) and (f) involve four single-derivative vertices. Diagram (e) requires one of the vertices to emit three η^(') lines, i.e., the last three term in brackets in (<ref>), 
and the three remaining producing only one meson line, which in total entails three  η^(') propagators. It involves terms proportional to M_0^6 Λ_2/M^6_0η^('), which are (δ), and terms multiplying M_0^6 v_2^(2)/M^6_0η^('),  M_0^6 v_4^(0)/M^6_0η^(') and M_0^4 Λ^2_2/M^6_0η^('), contributing at (δ^2). Finally, at LO in U(3) diagram (f) implies an interaction vertex with four internal legs, which in turn involves four η^(') propagators evaluated at zero momentum. It implies that only mass terms contribute to the interaction vertex, leading to a total contribution proportional to F^2 M_0^8 m_q,s/M^8_0η^(')  at (1), F^2 M_0^6 Λ_2 m_q,s^2/M^8_0η^(') at  (δ) or terms proportional F^2 M_0^6 v_2^(2) m_q,s^2/M^8_0η^(') and F^2 M_0^4 Λ_2^2 m_q,s^3/M^8_0η^(') at  (δ^2). 


Furthermore, since diagram (f) is the only one contributing at (1) and (δ), it can be also dressed at higher orders in two different ways. On the one hand, the four-point interaction vertex can be dressed by including the higher-order diagrams depicted in Fig. <ref>. 
Namely, the LO four-point vertex can be replaced by its NLO counterpart in the U(3) expansion, i.e., by including a LEC, Fig. <ref> (a).  Since the propagators are evaluated at p=0 this diagram only involves the LECs L_6, L_7, L_8 from the L_4 Lagrangian and C_19, C_31 from L_6.  In addition, one might add a loop to the vertex either by including a π, K or η^(') tadpole, Fig. <ref> (b), or a rescattering diagram Fig. <ref> (c). Note that in these two cases off-shell momentum-dependent terms enter. In addition, the one-loop contribution in Fig. <ref> (c) is evaluated at p=0, which leads to the function M_0i^2nlogM_0i^2/μ^2, where M_0i denote a π, K, η or η' and n=0,1,2 depending on the fields and derivatives running in the loop. 
Terms proportional to the LEC L_8 are (δ). The remaining LECs and loop topologies enter only at (δ^2).

On the other hand, one or more of the η^(') propagators can be dressed at NLO in the U(3) expansion. The η^(') self energies were computed in <cit.> in U(3) ChPT. At zero momentum they involve the LECs Λ_2, L_8 at  (δ) and L_6, L_7, L_25, v_2^2, L_8^2, L_8 Λ_2, Λ^2_2 at (δ^2), Fig. <ref> (a), as well as further π, K or η^(') tadpole diagrams entering at (δ^2), Fig. <ref> (b). Finally, the product of the topology in Fig. <ref> (a) proportional to L_8 and in Fig. <ref> (a) multiplying Λ_2 and L_8 also contribute at  (δ^2). 

All together, the LO contribution in the δ expansion to the fourth-order cumulant in the isospin limit reads

    c_4^U(3),  LO,  IL    =F^2M_0^8B_0[^4/162 M_0η^8(m_q(^4-4√(2)^3+12^2^2-8√(2)^3+4^4).
         .+2m_s(4^4+8√(2)^3+12^2^2+4√(2)^3+^4))
         -2^3/81 M_0η^6 M_0η'^2(m_q(√(2)^4-5^3+3√(2)^2^2+2^3-2√(2)^4).
         .-2m_s(2√(2)^4+2^3-3√(2)^2^2-5^3-√(2)^4))
         +^2^2(m_q+2m_s)/27 M_0η^4 M_0η'^4(2^4-2√(2)^3-3^2^2+2√(2)^3+2^4)
         +2^3/81 M_0η^2 M_0η'^6(m_q(-2√(2)^4-2^3+3√(2)^2^2+5^3+√(2)^4).
         .+2m_s(√(2)^4-5^3+3√(2)^2^2+2^3-2√(2)^4))
         +^4/162 M_0η'^8(m_q(4^4+8√(2)^3+12^2^2+4√(2)^3+^4).
         .+2m_s(^4-4√(2)^3+12^2^2-8√(2)^3+4^4))].


Expressing once more the tree-level η and η' masses and mixing angle in terms of M_0, m_q and m_s, (<ref>) simplifies remarkably leading to

    c_4^U(3),  LO,  IL =-B_0 F^2 M_0^8 m_q m_s(m_q^3 + 2 m_s^3))/(M_0^2(m_q+ 2m_s) + 6B_0 m_qm_s)^4=-Σm̅_IL^4/m̅^[3]_IL(M_0^2+6B_0m̅_IL)^4=-Σm̂_IL^4/m̅^[3]_IL,

 where m̅_IL and m̂_IL were defined in (<ref>) and (<ref>), respectively, and 

    m̅^[n]_IL= [2/m_q^n+1/m_s^n]^-1.


Like the topological susceptibility, c_4 vanishes in the chiral limit m_q→ 0 pointing out once more that it is a chiral quantity. The SU(3) result for the fourth-order cumulant at LO in the chiral expansion can be recovered by taking the limit M_0→∞. In that case m̂_IL→m̅_IL and one retrieves the results in <cit.>.
In addition, one can also study the opposite limit, i.e., the quenched approximation for which M_0<<m̂. In that case one obtains

    c_4^U(3), LO, IL= -F^2M_0^8/1296 B_0^3 m̅_IL^[3]+(1/N_c^4).
  

The NLO (δ) results are always proportional to the (N_c) LECs Λ_2 and L_8, while NNLO (δ^2) results involve the remaining pieces. Namely, terms proportional to the LECs L_6, L_7, L_25, C_29, C_31, v_4^(0), v_2^(2), L_8^2,  L_8 Λ_2, Λ^2_2 and meson logarithms. 
The renormalization of the LECs <cit.> render a finite and scale-independent result. The explicit formulas are too long to be displayed here and they are provided as supplementary material.
  In addition, in the SU(3) M_0→∞ limit we recover the NLO results in the chiral expansion given in <cit.>. 

Finally, we can also study the large-N_c expansion of the fourth-order cumulant. While the LO N_c behavior of χ_top was well established time ago <cit.>, the large-N_c behavior of c_4 is still under debate. On the one hand, in <cit.> it was argued from the large-N_c structure of the vacuum energy density that the fourth-order cumulant should scale as (1/N_c^2). On the other hand, an explicit calculation  based on U(N_f) ChPT at NLO in the chiral expansion for degenerate quark masses <cit.> suggests that it goes as (1/N_c^3).
Our LO  (1) and NLO  (δ) results do indeed reproduce the predictions in <cit.>, since we obtain that they are (1/N_c^3) and (1/N_c^4), respectively. 
Nevertheless, at NNLO (δ^2) the fourth-derivative contribution in (<ref>) involves the contact term -6 F^2v_4^0∼(1/N_c^2). In fact, one can show that it is the only term contributing to c_4 at this large-N_c order. Namely, the leading-N_c dependence of any U(N_f) chiral operator Ô(X) involving the field X in (<ref>) is given by <cit.>

    Ô(X)= N_c^2-#(Tr)-#(X),

where #(Tr) and #(X) denote the number of chiral flavor traces and powers of the operator X in (<ref>), respectively.
Thus, any higher order operators in X will be suppressed both by their large-N_c counting  and the pion decay constants F. 
It also implies that any nth-order cumulant of the topological charge distribution should scale as N_c^2-n. All together, we obtain

    c_4^U(3), IL=   -6F^2 v_4^0-F^2M_0^8/1296B_0^3m̅^[3]_IL
    
    -4F^2M_0^2/m̅_ILB_0[v_4^0+M_0^2/216(v_2^(2)-Λ_2^2)+M_0^6/27F^4(3C_19^r+4C_31^r+160/3L_8^r 2)]
       +(1/N_c^4) + (δ^3).

where, as explained, the first term in the r.h.s. is (1/N_c^2) and the rest of the displayed terms are (1/N_c^3). 



§ NUMERICAL RESULTS


In Table <ref> we provide the numerical results for χ_top^1/4, (-c_4)^1/4 and b_2=c_4/(12χ_top) (the latter is defined following standard lattice analyses, see below) calculated in ChPT for SU(2), SU(3) and U(3) at different orders in the chiral or δ expansion. The numerical values for  the parameters involved, i.e., F, M_0, meson masses and the LECs L_6^r, L_7^r, L_8^r, Λ_1, Λ_2, C_19 and C_31  are taken from <cit.>. Note that the U(3) LECs do not correspond to the usual SU(3) ChPT quantities, but larger differences might be expected between them.  More precisely, we consider the values of the NNLOFit-B, i.e., their best fit to lattice results for the η and η' masses. However, the constants L_25^r, v_4^(0) and v_2^(2) are not included in <cit.>. L_25^r and v_2^(2) were estimated in <cit.> in an additional fit to lattice data for the η and η' parameters and they enter both in χ_top^1/4 and c_4. 
On the contrary, v_4^(0) remains unknown and it only contributes to the fourth-order cumulant. 
In addition, while these LECs play a very small role on the topological susceptibility, v_2^(2) and v_4^(0) have a much more sizable effect on c_4.
In fact, taking  L_25^r and v_2^(2) from <cit.> and assuming the NNLO U(3) correction for c_4 to agree within uncertainties with the NLO estimate, one obtains for v_4^(0)=218(10).  Nevertheless, in order to avoid any bias in the behavior of the U(3) expansion, we will set instead te value of v_4^(0) value to zero.
Thus, since there are no current estimates for v_4^(0) and it would be inconsistent to include only v_2^(2), we simple set the values of L_25^r, v_4^(0) and v_2^(2) to zero. 
Furthermore, we will neglect the NNLO U(3) corrections coming from Λ_2^2 since they are of the same other than the  v_2^(2) and v_4^(0) effects that we are ignoring. 


We also include the LECs uncertainties quoted in <cit.>. One can see in Table <ref> that the uncertainties of the U(3) LECs are much larger than the standard SU(2) and SU(3) errors, the main source of error coming from F^2.  In that sense, 
let us remark that the SU(3) and SU(2) values quoted in Table <ref> are obtained from the U(3) expressions by taking the M_0→∞ limit, plus the m_s→∞ one in SU(2), but keeping the numerical values of the U(3) LECs and their uncertainties.
The reason for this is that our main purpose here is  to calibrate the numerical effect of the η' as compared with the rest of the light degrees of freedom. Recall that, according to our previous discussion, the NNLO U(3) contribution includes the NLO SU(3) one in the limit where the η' is decoupled and the NLO SU(2) results when also kaons and eta decouple.  This has however the drawback of losing numerical precision with respect to the corresponding  purely SU(2) and SU(3) LECs.
For instance, recent estimates of the pure SU(2) calculation are χ_top^1/4=75.5(5) MeV at NLO  <cit.> and χ_top^1/4=75.44(34) MeV at NNLO <cit.> for typical LECs values of SU(2) and SU(3) <cit.>. Note that the NNLO calculated for SU(2) in <cit.>   includes O(p^6)= O(δ^3, δ^4) corrections, which are beyond our present analysis. In particular, we do not recover those corrections when taking the M_0→∞ and m_s→∞ limits. For that reason, we do not include those results in Table <ref> where, as explained, we aim to discuss the different contributions that are obtained as limiting cases of our present approach. In any case, it is pointed out in <cit.>, despite the uncertainties in the O (p^6) LECs, the NNLO SU(2) numerical corrections are one order of magnitude smaller than the NLO ones.  In the case of the b_2 coefficient, the LO SU(2) and SU(3) expressions are LEC independent and therefore, they are given without theoretical uncertainty in Table <ref>. The cancellation of the F^2  dependence  in the U(3) LO for b_2, e.g. from (<ref>) and (<ref>),  explains also its smaller error compared to higher orders. 

The results in Table <ref> are obtained in the isospin limit. In the case of the topological susceptibility they have to be compared with the lattice result [χ_top^latt]^1/4=73(9) MeV in that case <cit.>. As for lattice results on c_4,  customarily given in terms of the b_2 coefficient,   they are provided only for pure gauge SU(N) theories in <cit.> and for domain-wall N_f=2+1 fermions with large light quark masses  m_l/m_s≥ 0.25 <cit.>. Remarkably, the value of b_2 seems to be quite stable under those different approximations and close to the simplest SU(2) ChPT value, as pointed out in <cit.>. We quote for reference the isospin-limit value for N_c=3 gluodynamics, b_2=-0.0216(15) <cit.>.  Nevertheless, more accurate lattice determinations for  c_4 and b_2 for the physical N_f=2+1 case would be needed to make further claims. 





From the results in Table <ref> we also observe the following features.  First, in the three theoretical frameworks, the perturbative corrections remain reasonably under control. Second, although the SU(2) approach already reproduces the main contribution, the η' meson and mixing angle corrections that we are including in the present work are actually comparable to the kaon and η ones introduced in the SU(3) approach. For χ_top, those corrections lower the central value and get closer to the lattice prediction and so on for | b_2|. Actually,  we  see that the full U(3) calculation for both observables remains compatible with the lattice results within the range provided by the LECs uncertainties, which in the case of χ_top holds also for all the different approximations collected in Table <ref>. The latter confirms that these are are good chiral quantities in the sense that they can be accurately described within ChPT.  







§ ISOSPIN BREAKING CORRECTIONS TO Χ_TOP AND C_4



As mentioned in the introduction, isospin breaking corrections can become important for the topological charge distribution.
The main reason is that for θ≠ 0 and m_u≠ m_d, the constant field configuration that minimizes the vacuum energy density is not U_0=, but generally U_0=(e^ϕ_1,…,e^ϕ_N_f)  <cit.>.
In the SU(3) framework, the constraint U_0=1 leads to ∑_j ϕ_j=0.
However, within our present U(3) formalism, such constraint does not hold, since the determinant of U_0 is an additional degree of freedom <cit.>. The U_0 configuration should be such that the vacuum energy density ϵ_vac (θ) defined in (<ref>) is minimized,
which is indeed achieved for a constant value of U_0.
Therefore, to LO in the chiral expansion, we have  to consider just the usual NGB mass term plus the M_0^2 term in (<ref>) for the Euclidean action. Namely,

    ϵ_vac^LO(θ)   =   -F^2 B_0/2[U_0 M^†+ M U_0^†]-F^2/12M_0^2[iθ+log U_0]^2
    
       =   -F^2 B_0 ∑_j=1^N_f m_jcos[ϕ_j (θ)]+F^2 M_0^2/12(θ+∑_j=1^N_fϕ_j (θ))^2

with M= (m_u,m_d,m_s) the quark mass matrix and ϕ_j(θ) are such that they minimize ϵ_vac^LO(θ), i.e., 

    B_0 m_ksin[ϕ_k(θ)] +M_0^2/6(θ+∑_j=1^N_fϕ_j (θ))=0      (k=1,…,N_f)


Note that the solution to (<ref>) is equivalent to encode the θ dependence in a complex quark mass matrix Mexp(-iθ/N_f) as in <cit.> with the change of variable ϕ_j→ϕ_j +θ/N_f. 

Now, following the same procedure as in <cit.>, in order to solve the minimization problem we expand cosϕ_j in powers of ϕ_j. The reason for this is that we are only interested in the power expansion of ϵ_vac(θ) around θ=0 and hence around the solution ϕ_j=0 of  (<ref>) for θ=0. For N_f=2 and M_0=0 the solution of (<ref>) can be found in <cit.>. The solution of the system (<ref>) to (ϕ_j) for N_f=3 is  ϕ_j(θ)=ϕ_j^0 θ+(θ^3) with

    ϕ_j^0=-m_j (6B_0/M_0^2+1/m̅)       (j=u,d,s)

with m̅ defined in (<ref>). Replacing the above linear order in the vacuum energy density (<ref>) yields

    ϵ_vac^LO(θ)=ϵ_vac(0)+1/2χ_top^U(3),LOθ^2 + (θ^4)

with

    χ_top^U(3),LO= Σm̂

and

    m̂=M_0^2 m̅/M_0^2+6B_0m̅


 The above result is the extension to m_u≠ m_d of the LO U(3) result in (<ref>), which amounts to the replacement m̅_IL→m̅ and reproduces the LO result with isospin breaking in SU(3) in (<ref>) by replacing m̅→m̂. 

In order to provide a numerical estimate of the isospin breaking effect, we write (<ref>) as

    χ_top^U(3),LO=F^2 M_0^2/6+(1+z)^2/zM_0^2/M_0π^0^2+(1+z)M_0^2/(1+z)M_0 K^0^2-M_0π^0^2

where z=m_u/m_d. Using the central value z=0.485 of the recent lattice analysis <cit.>, we get [χ_top^U(3),LO]^1/4=72 MeV, to be compared with the U(3) LO value in Table <ref>, which corresponds to z=1 in (<ref>). Thus, the isospin correction to LO U(3) is within the 5% range and lies within the theoretical LO uncertainty. It is therefore numerically safe to consider isospin breaking only for the LO U(3) result in our present analysis.


Following the same approach to the next order in the θ expansion allows us to calculate the isospin breaking corrections to the fourth-order cumulant. Thus, we expand (<ref>) up to (ϕ_j^3) and write its solution as ϕ_j(θ)=ϕ_j^0θ + ϕ_j^1θ^3 + (θ^5) and keep only up to (θ^3) in the equation, thus solving linearly for the  ϕ_j^1. Replacing then the solution in (<ref>) yields:


    ϵ_vac^LO(θ)=ϵ_vac(0)+1/2χ_top^U(3),LOθ^2 + 1/24 c_4^U(3),LOθ^4+(θ^6)

with the fourth order cumulant


    c_4^U(3),LO=-Σm̂^4/m̅^[3]

where, following the notation of <cit.>, we have defined, consistently with (<ref>),


    m̅^[3]= [1/m_u^3+1/m_d^3+1/m_s^3]^-1.


The result (<ref>) for the fourth-order cumulant corresponds again to the SU(3) one in <cit.> with the replacement m̅→m̂. As in the previous case, we also can write in terms of measurable meson parameters:



    c_4^U(3),LO=-F^2 M_0^8(1+z)^3 [1/[(1+z) M_0K^0^2-M_0π^0^2]^3+1+z^3/z^3 M_0π^0^6]/[6+(1+z)M_0^2/(1+z) M_0K^0^2-M_0π^0^2+(1+z)^2/zM_0^2/M_0π^0^2]^4

which yields (-c_4^U(3),LO)^1/4=53 MeV, so that the correction lies also within the 5% level when compared with the isospin-limit values in Table <ref>. 


Finally, we recall that at the order we are considering the isospin-breaking corrections, temperature dependence  is absent since it only enters through loop contributions. 



§ FINITE TEMPERATURE DEPENDENCE



From our present U(3) ChPT analysis, we can straightforwardly include the temperature dependence coming from  meson loops. In the case of χ_top, loop effects only arise at NNLO in the δ expansion from tadpole contributions coming from the Euclidean tree-level propagator G_i(x=0). Its finite part reads

    μ_i(T)=M_0i^2/32π^2 F^2logM_0i^2/μ^2+g_1(M_0i,T)/2F^2
    
    g_1(M,T)=T^2/2π^2∫_M/T^∞ x  √(x^2-(M/T)^2)/e^x-1
     
where  i=π,K,η,η', M_0i are the tree level masses and μ is the  renormalization scale.

In the case of c_4, in addition to tadpoles, which enter again from the Euclidean tree-level propagator G_i(x=0) in Fig <ref> (b) but also from the six-point interaction vertex in Fig. <ref> (b), one has to take into account the  one-loop function depicted in Fig. <ref> (c), which finite part  can be written in terms of

    ν_i (T)=   F^2d/d M_0i^2μ_i(T) = 1/32π^2[1+logM_0i^2/μ^2]-g_2(M_0i,T)/2,
    
    g_2(M,T)=   1/4π^2∫_M/T^∞ dx   1/x1/e^x-1.


That said, we want to remark that the U(3) thermal expansion is based on a perturbative calculation and hence, it is only expected to converge at very low temperatures.
The fact that thermal corrections arise at NNLO in the δ expansion implies that the series breaks down as soon as thermal effects are sizable. 
In Figure <ref>, we  show the temperature dependence of χ_top(T) and b_2(T) for the same parameter values and errors used for the T=0 results in Table <ref>.  In addition, we also plot the lattice data results for the topological susceptibility obtained in <cit.> and <cit.> and for b_2 in <cit.>. In the latter case, lattice errors are  larger than for the susceptibility and we have actually not considered the data set for a=0.0824 fm which has even  larger errors than those showed in the figure. 

We see that the present U(3) ChPT analysis is consistent with the lattice within uncertainties, even beyond its applicability range, which as mentioned before lies at low temperatures, well below the transition. Consequently, its extrapolation close and above T_c has to be taken with care.  
However, the good agreement with the lattice observed in Figure <ref> reveals once more that accounting properly for the lightest meson degrees of freedom is crucial for the description of the topological susceptibility and the cumulant.  



This implies an important difference with other thermodynamic observables like the quark condensate, which accurate description based on effective theories requires the contribution of many hadronic states,
like for instance in the Hadron Resonance Gas approach <cit.>.   
Actually, in figure <ref>  we compare different orders of the finite-T ChPT approach for the topological susceptibility and we see that the SU(3) and U(3) calculations represent rather small deviations from the SU(2) one. 



In the case of b_2, our theoretical U(3) result in Figure <ref> is quite flat with temperature, and lattice data are less accurate. This indicates  that both quantities χ(T) and c_4(T) decrease with T in a roughly similar way. Nevertheless, the agreement between theory and lattice results is also quite remarkable. 

As discussed in section <ref>, another important issue regarding the temperature dependence is to what extent it can be approximated by just the scaling of the quark condensate, i.e.,  whether or not the second term in the  Ward Identity in (<ref>) can be ignored.  This is actually the case if one sticks to SU(2) at NLO in the chiral expansion, what has been used in <cit.>. 
In our present work, we are calculating  χ_top(T) in U(3) at NNLO, including SU(3) and SU(2) NLO as special cases. 
Thus, we can provide a much more accurate analysis in that respect. 
That  issue can be also of relevance for lattice analyses. If the quark condensate terms dominate, the combined use of (<ref>) and (<ref>) may help to relate χ_top with quantities much better determined in the lattice. Namely,
 
    χ_top=m_q2(2-m_q^2/m_s^2)^-1[2Δ_l,s-m_q(χ_p^ll-2χ_p^ss)]

   where 

    Δ_l,s= (T)-2m_q/m_s (T)

is the reduced quark condensate used in lattice calculations to eliminate finite-size divergences appearing in individual condensates <cit.>. 
The relation (<ref>) offers a way to measure indirectly χ_top, which is alternative to the usual method based on the WI 

    χ_top=m_q^2/4(χ_P^π-χ_P^ll)
 
which stems directly from (<ref>), with χ_P^π=-/m_q the pion susceptibility <cit.>.
The reduced quark condensate is a particular example of a well determined quantity in the lattice. 

The possible dominance of the quark-condensate term in the WI (<ref>) is also relevant for a current topic of discussion, which has been actively studied both theoretically and in the lattice <cit.>.
Namely, whether the chiral and U(1)_A restoration temperatures are close enough.  Since χ_top and  are meant to vanish at exact U(1)_A (asymptotically) and chiral O(4) restoration, respectively, their difference, encoded in χ_P^ll in (<ref>), provides a direct measure of the separation between the two transitions. 
In this sense, it is useful to recall the behavior of these quantities near the light chiral limit m_q→ 0^+ (M_π→ 0^+), where the effects of chiral symmetry restoration are meant to be enhanced. 
In NNLO U(3) ChPT <cit.>, χ_top=(m_q), while (T) and χ_P^ll (T) are both (1) quantities in the m_q→ 0^+ limit. Thus, in the equation
 
    χ_top/m_q=-1/4[+m_q χ_P^ll]

only the quark condensate contribution survives in the right-hand side of (<ref>) in the chiral limit, which supports its dominance at low temperatures.  
However, near the transition (T)→ 0^+  in the chiral limit, while χ_P^ll(T) changes much more slowly, since it is controlled by  a term proportional to T^2/M_0K^2. 
Thus, it brings up the question as to whether near the transition the χ_P^ll  term  can become important enough for physical masses. 
 


In Figure <ref>,   we compare the temperature dependence of the full U(3) topological susceptibility with the scaling of the light quark condensate calculated in the same framework,
which would correspond to neglect the χ_ll term in the WI (<ref>) and (<ref>). 
The values of the LECs and other parameters involved are the same as those used in Figure <ref>. 
In the same figure, we also show the simple NLO SU(2) scaling used in <cit.>, corresponding just to the NLO ChPT quark condensate, namely

    χ_top^SU(2), NLO(T)/χ_top^SU(2), NLO(0)=^SU(2),NLO(T)/^SU(2),NLO(0)=1-3/2F_π^2 g_1(M_π,T)



The results in Figure <ref> show that the contribution from the additional χ_ll term in the WI, although not large, may be significant as T approaches the transition point. 
The simple SU(2) description remains also close to the full U(3) one, which is a test of its robustness despite its simplicity. 
However, it is very important to point out that, although the topological susceptibility may be well described within a ChPT analysis  including only the lightest degrees of freedom, that is definitively not the case for the quark condensate. To show this explicitly, we have plotted in Figure <ref>  the quark condensate resulting from the HRG approximation provided in <cit.>, which includes hadron states with masses up to 2 GeV and provides a very good fit to lattice condensate data. One can see clearly a much larger departure from the scaling of the topological susceptibility than the one observed with the  ChPT expressions for the quark condensate, which happen to remain close to the topological susceptibility. In this sense, we remark  that the addition of degrees of freedom is expected to reduce drastically the chiral condensate, as expected from approximate chiral restoration at T_c.
Nevertheless, that may not be the case for the topological susceptibility, whose behavior is not directly related to chiral symmetry restoration but includes U(1)_A restoration features, describing lattice data just with the light degrees of freedom. 
Thus, in a full description of the hadron gas we do expect large deviations from the quark condensate scaling and hence significant contributions from the second term in (<ref>), which becomes large close to the transition point.
This also indicates that the U(1)_A symmetry is still sizeably broken at the chiral transition for physical quark masses. 
This analysis should prevent from the use of the topological susceptibility to extract the quark condensate at finite temperature. 

 

§ CONCLUSIONS

 
The main conclusions achieved in this work are the following:


  * We have provided a full calculation of the topological susceptibility and the fourth-order cumulant up to NNLO in U(3) Chiral Perturbation Theory. Our result allows one to consider the effect of the η' meson consistently, as well as the η-η' mixing angle dependence. As limits of interest,  we recover the SU(2) and SU(3) results when M_0→∞. In addition, we  have discussed the large-N_c corrections to both quantities. In the case of the topological susceptibility, we have provided the (1/N_c) correction to the Witten-Veneziano formula up to (δ^2) in the U(3) ChPT expansion, and so on for the (1/N_c^2) and (1/N_c^3) corrections to the fourth-order cumulant. 


  * We have estimated the η' corrections to χ_top and c_4 at zero temperature, which turn out to be of the same order as the K and η contributions. 
Furthermore, it provides results compatible with lattice analyses, consistently with the idea that the QCD topological charge is an observable well described by the expansion around the chiral limit provided by ChPT. 


  * Including the dominant isospin breaking effect in the vacuum misalignment, we have provided both χ_top and the fourth cumulant c_4 of the vacuum energy density expansion in the θ parameter to LO in the U(3) ChPT expansion. Numerical corrections due to isospin breaking remain below the 5% level for χ_top^1/4 and (-c_4)^1/4. 


  * We have calculated χ_top and c_4 at finite temperature up to NNLO in U(3) ChPT. The temperature dependence obtained for the topological susceptibility is consistent with lattice data, supporting again that this quantity is well described by a gas made only of light mesons, unlike for instance the quark condensate. We have also discussed the relation between these two quantities, which are connected through a Ward Identity valid at all temperatures. 
Although the quark condensate calculated within the same ChPT formalism seems to scale quite similarly to χ_top, we argue that this cannot be the case for the full hadron gas.
It reveals a sizable gap between the chiral transition and the U(1)_A one for physical quark masses, even though recent theoretical analysis show that those transition tend to coincide near the chiral limit for exact chiral restoration. 

  
 

§ ACKNOWLEDGMENTS

 We are very grateful to   F. K. Guo, U. G. Meißner,  J. J. Sanz-Cillero, Z. H. Guo and T. Vonk for very useful comments. Work partially supported by  research contract FPA2016-75654-C2-2-P  (spanish “Ministerio de Economía y Competitividad") and the Swiss National Science Foundation, project No. PZ00P2_174228. This work has also received funding from the European Union Horizon 2020 research and innovation programme under grant agreement No 824093. A. V-R acknowledges support from a fellowship of the UCM predoctoral program. 
 
 

§ RESULTS FOR Χ_TOP AT NLO AND NNLO IN U(3) CHPT IN THE ISOSPIN LIMIT

 
 
 We provide here the full results for the topological susceptibility in  the U(3) ChPT formalism at NLO and NNLO in the δ expansion.  
 
    χ_top^U(3),NLO   =   1/9 M_0η '^4{
     F^2 Λ _2 M_0^2 [2 M_0K^2 (-√(2))+M_0π^2 (2 √(2)+)](M_0^2 ^2-M_0η
       '^2).
       +   .8 M_0^4 ^2 L_8^r [4 M_0K^4 (2 ^2+^2-2 √(2))..
       -   ..4 M_0π^2 M_0K^2 (2 ^2+^2-2
       √(2))
       +3 M_0π^4]}
       +   1/9 M_0η^4{F^2 Λ _2 M_0^4 ^3 [2 M_0K^2 (+√(2))+M_0π^2 (-2 √(2))].
       
       +   .
       8 M_0^4 ^2
       L_8^r [4 M_0K^4 (^2+2 ^2+2 √(2))..
       -   ..4 M_0π^2 M_0K^2 (^2+2 ^2+2 √(2))+3 M_0π^4]}
       +   1/9 M_0η^2 M_0η '^2{
       F^2 Λ _2 M_0^2[
       (2 M_0K^2 (+√(2))+M_0π^2 (-2 √(2))) (M_0^2 ^2-M_0η '^2)..
       +   ..M_0^2[2 M_0K^2 (-√(2))+M_0π^2 (2 √(2)+)]]
       .
       +   . 64 M_0^4 M_0K^2 (M_0K^2-M_0π^2) L_8^r (-√(2)^2+√(2)^2-)
       }



 
    χ_top^U(3),NNLO   =   
     -16 M_0^2 Λ _2 L_8^r /27 M_0η^6 M_0η '^6{8M_0K^6 [(^3 √(2)^2+6 ^2 ...
         
       -   ...2 √(2)^3) (2 ^2 M_0^2-M_0η '^2) M_0η^6M_0η '^2 (2 (^5+2
       √(2)^4....
       -   ....
     7 ^2 ^3-2 √(2)^3 ^2+7 ^4 -√(2)^5) M_0^2
       ...
       +   ...
       (-2 ^4-√(2)^3+6 ^2 ^2+√(2)^3 -2 ^4) M_0η '^2) M_0η^4
       ..
       +   ..
         M_0η '^4 (2 (√(2)^5+7 ^4+2 √(2)^2 ^3-7 ^3 ^2-2 √(2)^4 +^5) M_0^2
       ...
       -   ...
    (2 √(2)^3+6^2+3 √(2)^2 +^3) M_0η '^2) M_0η^2
      ..
       +   ..
    2 ^3(2 √(2)^3+6^2+3 √(2)^2+^3) M_0^2 M_0η '^6]
       .
       -   .
       4 M_0π^2 M_0K^4 [(^3-6 √(2)^2+18 ^2 -8 √(2)^3) (2 ^2 M_0^2-M_0η '^2) M_0η^6
       ..
       +   ..
       2
       M_0η '^2 ((4 ^5+2 √(2)^4-25 ^2 ^3+4 √(2)^3 ^2+22
       ^4 ....
       -   ....
    4 √(2)^5) M_0^2+(-4 ^4+√(2)^3+9 ^2 ^2-√(2)^3 -4 ^4) M_0η '^2) M_0η^4
       ..
       +   ..
    M_0η '^4 (2 (4 √(2)^5+22^4-4 √(2)^2 ^3-25 ^3 ^2-2 √(2)^4 +4 ^5) M_0^2
      ...
       -   ...
    (8 √(2)^3+18
      ^2+6 √(2)^2 +^3) M_0η '^2) M_0η^2
       ..
       +   ..
    2 ^3 (8 √(2)^3+18 ^2+6 √(2)^2 +^3) M_0^2 M_0η '^6] 
       .
       +   .
    2 M_0π^4 M_0K^2 [(^3-3 √(2)^2+15 ^2 -11 √(2)^3) (2 ^2 M_0^2-M_0η '^2) M_0η^6
      ..
       +   ..
       M_0η '^2 ((8 ^5-5 √(2)^4-38 ^2
       ^3+32 √(2)^3 ^2+38 ^4 ....
       -   ....
    11 √(2)^5) M_0^2-4 (2 ^4-2 √(2)^3-3 ^2 ^2+2 √(2)^3 +2 ^4) M_0η '^2) M_0η^4
      ..
       +   ..
     M_0η '^4 ((11
       √(2)^5+38^4-32 √(2)^2 ^3-38 ^3 ^2
       ....
       +   ....
    5 √(2)^4 +8 ^5) M_0^2-(11 √(2)^3+15^2+3 √(2)^2 +^3) M_0η
       '^2) M_0η^2
       ..
       +   ..
    2 ^3 (11 √(2)^3+15^2+3 √(2)^2 +^3) M_0^2 M_0η
       '^6] 
      .
       +   .
    3 M_0π^6 [(+2 √(2)) (2
       ^2 M_0^2-M_0η '^2) M_0η^6
       ..
       +   ..
          2 (-√(2)^2++√(2)^2) M_0^2 M_0η '^2
       M_0η^4
       ..
       +   ..
     M_0η '^4 (2 (-√(2)^2++√(2)^2) M_0^2+(2 √(2)-) M_0η '^2) M_0η^2+2 ^3 (-2 √(2)) M_0^2 M_0η '^6]
       }
       -   128 M_0^4 (L_8^r)^2/27 F^2 M_0η^6 M_0η '^6{
        16 M_0K^8[^2 (^4-4 √(2)^3+12 ^2 ^2-8 √(2)^3 +4 ^4) M_0η^6
        ..
       +   ..
       (2 ^2-^2)
       (^3-6 ^2 +4 √(2)^3)
       M_0η '^2 M_0η^4
        ..
       +   ..
      (^2-2 ^2) (4 √(2)^3+6^2-^3) M_0η '^4 M_0η^2
        ..
       +   ..
       ^2 (4 ^4+8 √(2)^3+12 ^2 ^2+4 √(2)^3 +^4) M_0η '^6]
        .
       -   .
       32 M_0π^2 M_0K^6 [^2 (^4-4 √(2)^3+12 ^2 ^2-8 √(2)^3 +4 ^4) M_0η^6
    ..
       +   ..
       (2 ^2-^2) (^3-6 ^2 +4 √(2)^3) M_0η '^2 M_0η^4
       ..
       +   ..
    (^2-2 ^2) (4 √(2)^3+6^2-^3) M_0η '^4 M_0η^2
       ..
       +   ..
       ^2 (4 ^4+8 √(2)^3+12 ^2 ^2+4
       √(2)^3 +^4) M_0η '^6] 
       .
       +   .
    8 M_0π^4 M_0K^4[^2 (5 ^4-14 √(2)^3+33 ^2 ^2-22 √(2)^3 +14 ^4) M_0η^6
        ..
       +   ..
       (4 ^5+3 √(2)^4-29 ^2 ^3+16 √(2)^3 ^2+9 ^4 -11 √(2)^5) M_0η '^2
       M_0η^4
        ..
       +   ..
      (11 √(2)^5+9^4-16 √(2)^2 ^3-29 ^3 ^2-3
       √(2)^4 +4 ^5) M_0η '^4 M_0η^2
        ..
       +   ..
          ^2 (14 ^4+22 √(2)^3+33 ^2^2+14 √(2)^3 +5 ^4) M_0η '^6]
       .
       -   .
        24 M_0π^6 M_0K^2[^2 (^2-2
       √(2)+2 ^2) M_0η^6
       ..
       +   ..
       (√(2)^2--√(2)^2) M_0η '^2 M_0η^4
        ..
       +   ..
       (√(2)^2--√(2)^2) M_0η '^4 M_0η^2+^2 (2 ^2+2 √(2)+^2) M_0η '^6] 
        .
       +   .
    9 M_0π^8 (^2 M_0η^6+^2
       M_0η '^6)  }
       +   1/M_0η^4 M_0η '^4{(2 M_0K^2+M_0π^2) F^2 v_2^(2)[ (^2
       M_0^2-M_0η '^2) M_0η^2+ ^2 M_0^2 M_0η '^2]^2
        .
       +   .
    8 L_6^r M_0^4 (2 M_0K^2+M_0π^2)/9[
      ^2 (2 (^2-2 √(2)+2 ^2) M_0K^2
       ...
       +   ...
       (^2+4 √(2)-^2) M_0π^2) M_0η^4
       ..
       +   ..
       4 (√(2)^2--√(2)^2) (M_0K^2-M_0π^2) M_0η '^2 M_0η^2
       ..
       +   ..   
       ^2 (2 (2 ^2+2 √(2)+^2) M_0K^2+(-^2-4 √(2)+^2) M_0π^2) M_0η '^4
    ]
    .
       +   .
    8 C_19 M_0^4/3[
       ^2 (8 (^2-2 √(2)+2 ^2) M_0K^6
    ...
       -   ...   
       12 (^2-2 √(2)+2 ^2) M_0π^2 M_0K^4
       ...
       +   ...
      6 (^2-2 √(2)+2 ^2) M_0π^4 M_0K^2+(^2+4
       √(2)-^2) M_0π^6) M_0η^4
       ..
       +   ..
    4 (√(2)^2--√(2)^2) (4 M_0K^6-6 M_0π^2 M_0K^4+3 M_0π^4 M_0K^2-M_0π^6) M_0η '^2 M_0η^2
         ..
       +   ..
       ^2 (8 (2 ^2+2 √(2)
      + ^2) M_0K^6-12 (2 ^2+2 √(2)+^2) M_0π^2 M_0K^4
        ...
       +   ...
       6 (2 ^2+2 √(2)+^2) M_0π^4 M_0K^2+(-^2-4 √(2)+^2) M_0π^6) M_0η '^4
    ]
     .
       +   .
    16 C_31 M_0^4/9[
      ^2 (8 (^2-2 √(2)+2 ^2) M_0K^6
     ...
       -   ...  
       12 (^2-2 √(2)+2 ^2) M_0π^2 M_0K^4
        ...
       +   ...
       6 (^2-2 √(2)+2 ^2) M_0π^4 M_0K^2+(^2+4 √(2)-^2) M_0π^6) M_0η^4
        ..
       +   ..
       4 (√(2)^2--√(2)^2) (4 M_0K^6-6 M_0π^2 M_0K^4+3 M_0π^4 M_0K^2-M_0π^6) M_0η '^2 M_0η^2
       ..
       +   ..
       ^2
       (8 (2 ^2+2 √(2)+^2)
       M_0K^6-12 (2 ^2+2 √(2)+^2) M_0π^2 M_0K^4
       ...
       +   ...
       6 (2 ^2+2 √(2)+^2) M_0π^4 M_0K^2+(-^2-4 √(2)+^2) M_0π^6) M_0η '^4]
      .
       +   .
    8 L_7^r M_0^4/9[
       ^2 (8 ^2 (M_0K^2-M_0π^2)^2+^2 (2 M_0K^2+M_0π^2)^2
       ...
       +   ...
       4 √(2)(-2 M_0K^4+M_0π^2 M_0K^2+M_0π^4)) M_0η^4
       ..
       +   ..
       2 (4 (√(2)^2--√(2)^2) M_0K^4
        ...
       +   ...
       2 (-√(2)^2+10 +√(2)^2) M_0π^2 M_0K^2
        ...
       +   ...
       (-2 √(2)^2-7
      +2 √(2)^2) M_0π^4) M_0η '^2 M_0η^2
        ..
       +   ..
       ^2 (8 ^2 (M_0K^2-M_0π^2)^2+^2 (2 M_0K^2+M_0π^2)^2+4 √(2)(2 M_0K^4-M_0π^2 M_0K^2-M_0π^4)) M_0η
       '^4]
       .
       -   .
    8 L_25^r M_0^2/3[
       (4 (-√(2)) M_0K^4-4 (-√(2)) M_0π^2 M_0K^2+3  M_0π^4) M_0η^2
       ..
       +   ..
      (4 (√(2)+) M_0K^4
       ...
       +   ...
       -4 (√(2)+) M_0π^2 M_0K^2+3  M_0π^4) M_0η '^2) ((^2 M_0^2-M_0η '^2) M_0η^2
    
    
      + ^2 M_0^2 M_0η '^2]
       }
       -   M_0^4/54 M_0η^4{
          4 ^2 (2 (2 μ
       _η (T)+μ _η '(T)) ^4
        +  2 √(2)(4 μ _η(T)-μ _η '(T)) ^3
       ..
       +   ..
       3 ^2 (4 μ _η(T)-μ _η '(T)) ^2
       +2 √(2)((2 μ _η(T)+μ _η '(T)) ^2+3 μ _K(T)) ..
       +   ..
       6 μ _K(T)
       ^2 +^4 (μ _η(T)+2 μ _η '(T))) M_0K^2
    .
       +   .
       ^2M_0π^2 (-(7 μ _η(T)+2 μ _η '(T)) ^4+2 √(2)(μ
       _η '(T)-10 μ _η(T)) ^3
    ..
       +   ..   
       3 ((μ _η '(T)-4 μ _η(T)) ^2-2 μ
       _K(T)+3 μ _π(T)) ^2
    ..
       -   ..      
       2 √(2)((8 μ _η(T)+μ _η '(T)) ^2+6 μ _K(T)+9 μ _π(T)) +18  ^2 μ _π(T)
     ..
       +   ..        
       2 ^4 (μ _η(T)-μ _η '(T)))
       }
       -   M_0^4/54 M_0η '^4{
        4 M_0K^2^2
       [(2 μ _η(T)+μ _η '(T)) ^4-2 √(2)(μ _η(T)+2 μ _η
       '(T)) ^3
       ..
       +   ..    
       (6 μ _K(T)-3 ^2 (μ _η(T)-4 μ _η '(T))) ^2+2 √(2)(^2 (μ _η(T)-4 μ _η '(T))-3 μ
       _K(T)) ..
       +   ..    
       2 ^4 (μ _η(T)+2 μ _η '(T))] +M_0π^2^2[2 (μ _η '(T)-μ _η(T)) ^4
       ..
       +   ..    
     2 √(2)(μ
       _η(T)+8 μ _η '(T)) ^3+3 ((μ _η(T)-4 μ _η '(T)) ^2+6 μ
       _π(T)) ^2
       ..
       +   ..    
       2 √(2)(-(μ _η(T)-10 μ _η '(T)) ^2+6 μ _K(T)+9 μ _π(T))
       ..
       -   ..    
       ^2 ((2 μ _η(T)+7 μ
       _η '(T)) ^2+6 μ _K(T)-9 μ _π(T))]
       }
       -   M_0^4/54 M_0η^2 M_0η '^2{
         8  M_0K^2[√(2)(2 μ _η(T)+μ _η '(T)) ^4+(2 μ _η(T)-5 μ _η '(T)) ^3
       ..
       +   ..
       3 √(2)((μ _η '(T)-μ _η(T)) ^2+μ _K(T)) ^2+((2 μ _η '(T)-5 μ _η(T))
       ^2+6 μ _K) 
       ..
       -   ..
       √(2)^2 ((μ _η(T)+2 μ
       _η '(T)) ^2+3 μ _K(T))] 
        .
       +   .
       2  M_0π^2
       [√(2)(4 μ _η(T)+5 μ _η '(T)) ^4+(8 μ _η(T)+μ _η
       '(T)) ^3
       ..
       +   ..
       3 √(2)((μ _η(T)-μ _η '(T)) ^2+2 μ _K(T)+3 μ _π(T)) ^2
       ..
       +   ..
     ((μ _η(T)+8 μ _η '(T)) ^2+6 
       μ _K(T)+9 μ _π(T)) 
        ..
       -   ..
       √(2)^2 ((5 μ _η(T)+4 μ _η '(T)) ^2+6 μ _K(T)+9 μ _π(T))]
      }



  
99

Witten:1979vv 
  E. Witten,
 
  Nucl. Phys. B 156, 269 (1979).
  
  
Veneziano:1979ec 
  G. Veneziano,
  
  Nucl. Phys. B 159, 213 (1979).
 
  
  
  
DiVecchia:1980yfw 
  P. Di Vecchia and G. Veneziano,
  
  Nucl. Phys. B 171, 253 (1980).
 
  
  
  
Gasser:1984gg 
  J. Gasser and H. Leutwyler,
  
  Nucl. Phys. B 250, 465 (1985).

  
  
  
Leutwyler:1992yt 
  H. Leutwyler and A. V. Smilga,
  
  Phys. Rev. D 46, 5607 (1992).
  
  
  
  
  
  Aoki:2007pw 
  S. Aoki et al. [JLQCD and TWQCD Collaborations],
  
  Phys. Lett. B 665, 294 (2008).

 
  
  
 
  
  
Bernardoni:2010nf 
  F. Bernardoni, P. Hernandez, N. Garron, S. Necco and C. Pena,
  
  Phys. Rev. D 83, 054503 (2011).


  

 
 Allton:2007hx 
  C. Allton et al. [RBC and UKQCD Collaborations],
  
  Phys. Rev. D 76, 014504 (2007).

 
  
 
 
  Chiu:2008jq 
  T. W. Chiu et al. [TWQCD Collaboration],
  
  Phys. Lett. B 671, 135 (2009).

 
 
  
Bernard:2012fw 
  V. Bernard, S. Descotes-Genon and G. Toucas,
  
  JHEP 1206, 051 (2012).

 
  
  
  Aoki:2016frl 
  S. Aoki et al.,
  
  Eur. Phys. J. C 77, no. 2, 112 (2017).

 
  
  
  Bonati:2015vqz 
  C. Bonati, M. D'Elia, M. Mariti, G. Martinelli, M. Mesiti, F. Negro, F. Sanfilippo and G. Villadoro,
  
  JHEP 1603, 155 (2016).

 
  
  
   
Petreczky:2016vrs 
  P. Petreczky, H. P. Schadler and S. Sharma,
  
  Phys. Lett. B 762, 498 (2016).

 
  
  
  
Borsanyi:2016ksw 
  S. Borsanyi et al.,
  
  Nature 539, no. 7627, 69 (2016); 

  [arXiv:1606.07494 [hep-lat]] for additional information.
  
  
  Dimopoulos:2018xkm 
  P. Dimopoulos et al.,
  
  Phys. Rev. D 99, 034511 (2019).

 
  
  
  
Burger:2018fvb 
  F. Burger, E. M. Ilgenfritz, M. P. Lombardo and A. Trunin,
  
  Phys. Rev. D 98, no. 9, 094501 (2018).

 
  
  
  
  
Hansen:1990yg 
  F. C. Hansen and H. Leutwyler,
  
  Nucl. Phys. B 350, 201 (1991).
 
  
  
  Nicola:2016jlj 
  A. Gómez Nicola and J. Ruiz de Elvira,
  
  JHEP 1603, 186 (2016).

 
  
  
  
  
Azcoiti:2016zbi 
  V. Azcoiti,
  
  Phys. Rev. D 94, no. 9, 094505 (2016).

 
  
     
  Nicola:2018vug 
  A. Gómez Nicola and J. Ruiz De Elvira,
  
  Phys. Rev. D 98, no. 1, 014020 (2018).

 
  
  
  

  Buchoff:2013nra 
  M. I. Buchoff, M. Cheng, N. H. Christ, H.-T. Ding, C. Jung, F. Karsch, Z. Lin and R. D. Mawhinney et al.,
  
  Phys. Rev. D 89, 054514 (2014).

  


Bhattacharya:2014ara 
  T. Bhattacharya, M. I. Buchoff, N. H. Christ, H.-T. Ding, R. Gupta, C. Jung, F. Karsch and Z. Lin et al.,
  
  Phys. Rev. Lett.  113, 082001 (2014).
 
  
  
  
  
GomezNicola:2017bhm 
  A. Gómez Nicola and J. Ruiz de Elvira,
  
  Phys. Rev. D 97, no. 7, 074016 (2018).


  
  
  
  
  
Mao:2009sy 
  Y. Y. Mao et al. [TWQCD Collaboration],
  
  Phys. Rev. D 80, 034502 (2009).

 
  




diCortona:2015ldu 
  G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro,
  
  JHEP 1601, 034 (2016).

 
  
  
  
  Gorghetto:2018ocs 
  M. Gorghetto and G. Villadoro,
  
  JHEP 1903, 033 (2019).


  
  
  
Peccei:1977hh 
  R. D. Peccei and H. R. Quinn,
  
  Phys. Rev. Lett.  38, 1440 (1977).

  
  
  
Weinberg:1977ma 
  S. Weinberg,
  
  Phys. Rev. Lett.  40, 223 (1978).

  
  
  
Spalinski:1988az 
  M. Spalinski,
  
  Z. Phys. C 41, 87 (1988).
 
  
  
  
Bernard:2012ci 
  V. Bernard, S. Descotes-Genon and G. Toucas,
  
  JHEP 1212, 080 (2012).
 
 
  
  
  
Guo:2015oxa 
  F. K. Guo and U. G. Meißner,
  
  Phys. Lett. B 749, 278 (2015).


  
   
  
Vicari:2008jw 
  E. Vicari and H. Panagopoulos,
  
  Phys. Rept.  470, 93 (2009).


  
  
  
Panagopoulos:2011rb 
  H. Panagopoulos and E. Vicari,
  
  JHEP 1111, 119 (2011).


  

  Bonati:2016tvi 
  C. Bonati, M. D'Elia, P. Rossi and E. Vicari,
  
  Phys. Rev. D 94, no. 8, 085017 (2016).

 
  
  
  
  
  
  
   
Vonk:2019kwv 
  T. Vonk, F. K. Guo and U. G. Meißner,
  
  JHEP 1906, 106 (2019).


  

  
  
Bazavov:2012qja 
  A. Bazavov et al. [HotQCD Collaboration],
  
  Phys. Rev. D 86, 094503 (2012).

 
  
  
 
  
  
  
Bonati:2013tt 
  C. Bonati, M. D'Elia, H. Panagopoulos and E. Vicari,
  
  Phys. Rev. Lett.  110, no. 25, 252003 (2013).
 
 
  
    
  
   
  HerreraSiklody:1996pm 
  P. Herrera-Siklody, J. I. Latorre, P. Pascual and J. Taron,
  
  Nucl. Phys. B 497, 345 (1997).
 
  
  

Kaiser:2000gs 
  R. Kaiser and H. Leutwyler,
  
  Eur. Phys. J. C 17, 623 (2000).
 
  
  
  
  
Nicola:2010xt 
  A. Gómez Nicola and R. Torres Andres,
  
  J. Phys. G 39, 015004 (2012).

 
  



Escribano:2000ju 
  R. Escribano, F. S. Ling and M. H. G. Tytgat,
  
  Phys. Rev. D 62, 056004 (2000).

 
  
  


  
  
  
Gu:2018swy 
  X. W. Gu, C. G. Duan and Z. H. Guo,
  
  Phys. Rev. D 98, no. 3, 034007 (2018).

 
  
  
  
  
Ishii:2016dln 
  M. Ishii, H. Kouno and M. Yahiro,
  
  Phys. Rev. D 95, no. 11, 114022 (2017).


  


Kotov:2019dby 
  A. Y. Kotov, M. P. Lombardo and A. M. Trunin,
  
  Phys. Lett. B 794, 83 (2019).

 
  
  


    Rosenzweig:1979ay 
 C. Rosenzweig, J. Schechter and C. G. Trahern,

  Phys. Rev. D 21, 3388 (1980).

        
  
Guo:2012ym 
  Z. H. Guo, J. A. Oller and J. Ruiz de Elvira,
  
  Phys. Lett. B 712, 407 (2012);

  
  
  
  


  
  Phys. Rev. D 86, 054006 (2012).


  
  

  Guo:2015xva 
  X. K. Guo, Z. H. Guo, J. A. Oller and J. J. Sanz-Cillero,
  
  JHEP 1506, 175 (2015).
 
 
  
  

GomezNicola:2010tb 
  A. Gomez Nicola, J. R. Pelaez and J. Ruiz de Elvira,
  
  Phys. Rev. D 82, 074012 (2010);


  
  



  
  Phys. Rev. D 87, 016001 (2013).


  
  

 

Leutwyler:1996np 
  H. Leutwyler,
  
  Phys. Lett. B 374, 181 (1996).

 
  
  

  
Aoki:2019cca 
  S. Aoki et al. [Flavour Lattice Averaging Group],
  
  arXiv:1902.08191 [hep-lat].
  
  
  
 Karsch:2003vd 
  F. Karsch, K. Redlich and A. Tawfik,
  
  Eur. Phys. J. C 29, 549 (2003).
 
 
  
  
Tawfik:2005qh 
  A. Tawfik and D. Toublan,
  
  Phys. Lett. B 623, 48 (2005).
  
  
  

Huovinen:2009yb 
  P. Huovinen and P. Petreczky,
  
  Nucl. Phys. A 837, 26 (2010).
 
 
  
 
Jankowski:2012ms 
  J. Jankowski, D. Blaschke and M. Spalinski,
  
  Phys. Rev. D 87, 105018 (2013).


  
 

Nicola:2013vma 
  A. Gomez Nicola, J. Ruiz de Elvira and R. Torres Andres,
  
  Phys. Rev. D 88, 076007 (2013).


  
  

 
Shuryak:1993ee 
  E. V. Shuryak,
  
  Comments Nucl. Part. Phys.  21, no. 4, 235 (1994).

  


Cohen:1996ng 
  T. D. Cohen,
  
  Phys. Rev. D 54, R1867 (1996).

  
  


Lee:1996zy 
  S. H. Lee and T. Hatsuda,
  
  Phys. Rev. D 54, R1871 (1996).
 
 
  

Meggiolaro:2013swa 
E. Meggiolaro and A. Morda,
  
Phys. Rev. D 88, no. 9, 096010 (2013).
 
 
  
  
 Aoki:2012yj
  S. Aoki, H. Fukaya and Y. Taniguchi,
  
  Phys. Rev. D 86 (2012) 114512.

 

Cossu:2013uua 
  G. Cossu, S. Aoki, H. Fukaya, S. Hashimoto, T. Kaneko, H. Matsufuru and J. I. Noaki,
  
  Phys. Rev. D 87, no. 11, 114514 (2013)
  Erratum: [Phys. Rev. D 88, no. 1, 019901 (2013)].
  
 







    
    
Dick:2015twa 
  V. Dick, F. Karsch, E. Laermann, S. Mukherjee and S. Sharma,
  
  Phys. Rev. D 91, no. 9, 094504 (2015).

 
  
    
  
Tomiya:2016jwr 
  A. Tomiya, G. Cossu, S. Aoki, H. Fukaya, S. Hashimoto, T. Kaneko and J. Noaki,
  
  Phys. Rev. D 96, no. 3, 034509 (2017).
  
 
  
    
  
Brandt:2016daq 
  B. B. Brandt, A. Francis, H. B. Meyer, O. Philipsen, D. Robaina and H. Wittig,
  
  JHEP 1612, 158 (2016).

  
  
  
  
  
  




