Thermal behavior of effective U A (1) anomaly couplings in reflection of higher topological sectors

Thermal behavior of effective, chiral condensate-dependent U A (1) anomaly couplings is investigated using the functional renormalization group approach in the N f = 3 flavor meson model. We derive flow equations for anomaly couplings that arise from instantons of higher topological charge, dependent also on the chiral condensate. These flow equations are solved numerically for the | Q | = 1 , 2 topological sectors at finite temperature. Assuming that the anomaly couplings at the ultraviolet scale may also exhibit explicit temperature dependence, we calculate the thermal behavior of the effective potential. In accordance with our earlier study, [G. Fejos and A. Patkos, Phys. Rev. D 105 , 096007 (2022)], we find that for increasing temperatures, the anomalous breaking of chiral symmetry tends to strengthen toward the pseudocritical temperature ( T C ) of chiral symmetry breaking. It is revealed that below T C , around ∼ 10% of the U A (1) breaking arises from the | Q | = 2 topological sector. Correspondingly, a detailed analysis on the thermal behavior of the mass spectrum is also presented.


I. INTRODUCTION
The axial U A (1) subgroup of the approximate U V (N f ) × U A (N f ) chiral symmetry of quantum chromodynamics (QCD) with N f quark flavors is known to be broken anomalously in the ground state of the system [1,2].Microscopic origin of this anomaly, at least at sufficiently high temperatures, can be well described in the dilute instanton gas picture.However, at lower temperatures, toward the pseudocritical point (T c ), the eigenvalue spectrum of the Dirac operator is rather different to its high T counterpart, revealing that the instanton gas approximation definitely breaks down [3,4].The correct picture in such regimes is assumed to be more of an instanton liquid [5], where the instanton density and radius play a central role.Though the understanding of the underlying mechanism of the U A (1) breaking has significantly improved over time, the actual fate of the anomaly along the complete temperature axis still remains an open question.
Since the isotriplet pseudoscalar meson (π) and the isotriplet scalar meson (a 0 ) are related through the axial U A (1) transformation, it is natural to quantify the degree of the anomaly breaking by either the m a0 − m π mass or χ π − χ a0 susceptibility differences.In a chirally symmetric background, where the π and σ excitations degenerate, one notes that the former difference is equivalent to the disconnected part of the total chiral susceptibility, χ π − χ a0 = χ π − χ σ + χ disc → χ disc [6], leaving the quantification of the anomaly to determine χ disc .
Lattice QCD studies have not yet reached a firm con-clusion on the presence of the anomaly at the pseudocritical point.Refs.[7][8][9], with the use of domain-wall fermions, concluded that the anomaly is present even beyond the pseudocritical temperature.Ref. [10], using Wilson fermions, shows that U A (1) symmetry is effectively restored at T c in the chiral limit (for N f = 2).Eigenvalue spectrum analyses argue that the anomaly is present even beyond the critical temperature [11], and using highly improved staggered quark action also shows that U A (1) is broken at 1.6T c [12].The same finding is presented in Ref. [13] just above T c , using analyses of eigenvalue densities.Ensembles generated by twoflavor (Möbius) domain wall sea quarks with the eigenvalue analysis show, however, that the obtained results are consistent with the U A (1) symmetry being restored at T c , at least in the chiral limit [14,15].Predictions of the density of state method in the SU (3) pure gauge theory show that even the |Q| = 2 topological sector has a non-negligible contribution to the topological susceptibility even at 1.5T c [16].For a comprehensive review on recent developments along these directions, the reader is referred to Ref. [17].
The order of the chiral transition for zero quark masses may also contain information on the U A (1) restoration at the critical temperature.Recently, there have been indications that the chiral phase transition is of second order for both N f = 2, 3 in the chiral limit [32][33][34] in contradiction with earlier studies [35], especially those using the perturbative renormalization group in the ϵ expansion [36].The latter shows that for N f = 3, irrespectively of the fate of the U A (1) breaking at the critical point, the transition can only be of first order.Recently, using the functional renormalization group (FRG) technique, however, it was shown that the transition can appear to be second order after all, but only if the anomaly is very weak at T c [37].
In our earlier studies [38,39] we argued that in the N f = 3 low energy effective meson model, fluctuations tend to make the anomaly stronger with respect to the temperature.The newly found mechanism is related to the resummation of an infinite number of anomaly breaking operators in the functional renormalization group (FRG) formalism, which effectively made the standard Kobayashi-Maskawa-'t Hooft (KMT) coupling chiral condensate dependent.Through such dependence one observed a strengthening of the effective anomaly coupling as the condensate evaporated, showing that mesonic fluctuations are working against instanton contributions, which, at least at asymptotically large T definitely lead to the restoration of U A (1).
The standard KMT coupling and the corresponding determinant term arise from instanton contributions of the underlying theory that carry Q = ±1 topological charge.In [40] it is explicitly shown for N f = 2 that instantons with higher winding numbers lead to higher powers of the determinant term, which are typically considered nonrenormalizable, and thus entirely dropped in the effective description.Even though one might expect that the value of these couplings in the ultraviolet (even being zero) does not effect significantly the physics of the infrared, the corresponding terms do get generated at low energies, and in principle should not be neglected in the effective action.The main goal of this study is to derive scale evolution equations for those anomalous terms that arise from instantons with any winding number, and on top of that, realize a resummation in terms of chirally invariant operators to make the former condensate dependent.Our aim is also to quantify the importance of the higher topological sectors in the effective model setting.
We wish to emphasize that a consistent treatment of the U A (1) breaking, at least considering the physical point, can only be dealt with via the N f = 3 scenario.In case of N f = 2, the strange sector is completely decoupled as the s-quark mass is formally taken to be infinity.In such a case, the determinant in effect becomes a mass term, and the corresponding anomaly parameter also needs to be set to infinity in order for one of the O(4) multiplets can be dropped, as required by, e.g. the m K < m η ′ mass relation in the physical point.(Keep in mind that because of suppressing the strange content, by consistency m K → ∞, thus m η ′ → ∞.)Even if one wishes to keep both multiplets [40,41], one faces the fact that beyond T c , by the effective evaporation of the (non-strange) chiral condensate, keeping track of anomalous contributions in the effective potential is no longer possible.The N f = 3 model, however, does, through the very presence of the strange condensate, allow to investigate the effective KMT coupling even beyond T c , which proves to be crucial in understanding characteristic features of the high temperature meson spectrum.
The paper is organized as follows.In Sec.II, we discuss our model setup and introduce the corresponding ansatz for the effective potential, with particular emphasis on the resummation that leads to condensate dependent effective couplings.In Sec.III, we discuss how to obtain the scale evolution of the complete effective potential via the use of various background fields for the construction of the FRG flow equations (presented explicitly in the Appendix).Section IV contains the numerical results, where we used a model parameter set that corresponds to the physical point.The reader also finds a detailed analysis on the temperature dependence of the mesonic spectrum in light of the anomaly evolution.Section V is devoted for summary.

II. MODEL SETUP
Our investigations involve a mesonic dynamical variable M = (s i + iπ i )T i , where T i (i = 0, ..., 8) are the U (3) generators, and s i (π i ) refer to the scalar (pseudoscalar) fields.A chiral transformation is represented as M → LM R † , where L (R) are left (right) handed transformations.The effective potential, V , of our model can only depend on chirally invariant combinations of M .For N f = 3, a possible set of independent operators are the usual combinations.Higher order invariants can be expressed in terms of (1).The U A (1) breaking can be included via the Kobayashi-Maskawa-'t Hooft (KMT) term, where it is important to mention that the plus sign in the rhs of ( 2) is due to parity reasons (a minus sign would lead to a parity odd combination), and also that ∆ ≡ (det M † − det M )2 is not independent from (2) and (1).That is, the potential, both classical or quantum, can only depend on ρ, τ, ρ 3 , ∆.Note that in the broken phase chiral symmetry shows the U (3) × U (3) → U (3) breaking pattern, which is realized by the condensate ⟨M ⟩ ∼ 1 , i.e. it is proportional to the unit matrix.In such backgrounds τ = 0 = ρ 3 , i.e., even if one includes explicitly symmetry breaking terms (arising from nonzero quark masses) that somewhat shift the vacuum expectation value of M , it is expected that neither τ nor ρ 3 carry significant contributions in V .Therefore, it is expected to make sense to perform the following chiral invariant expansion [42]: where α 1 , α 2 ∈ Z.Note that in the classical version of the model, based on perturbative renormalizability, ρ 3 is dropped, and an expansion in terms of ρ and ∆ is also performed: which is the usual potential of the three flavor linear sigma model with real parameters m 2 , λ 1 , λ 2 , a.In the original variable M , (4) yields with λ1 = λ 1 −λ 2 /3.Note that at high momentum scales the form of ( 5) is certainly correct1 , but for low scales it is very restrictive, and thus if one includes quantum and thermal fluctuations, the more general form, (3) should be considered, at least via some approximation.From here onward, we entirely drop the ρ 3 dependence due to the reason described above, but we do keep τ , as it corresponds to a renormalizable operator that should be present at the UV scale.Its coefficient is, in principle, ρ and ∆ dependent, but we consider only the former: which will be the main starting point of our investigations.Even though (6) looks very simple, it is only due to our compact notations.Note that, the infinite resummation (6) realizes is threefold: Here U (α,β) and C (α) can be thought of as usual coupling constants.This demonstrates that the simple form of ( 6) contains an infinite number of couplings.Our goal is to determine these couplings through obtaining the functions U (ρ, ∆) and C(ρ).
A comment on the anomalous terms is now in order.It was shown by Pisarski and Rennecke that the ∼ ∆ α terms in the potential arise from multi instanton contributions of the underlying theory, with |Q| = α topological charges.Emergence of these interactions for N f = 2 can be found in [40].In what follows, we investigate the importance of the interactions generated by multi instanton contributions, in terms of an effective mesonic description at finite temperature.In the next section, using the functional renormalization group (FRG) formalism, we derive and then solve scale evolution equations for U (ρ, ∆) and C(ρ).

III. RG FLOWS
In the core of the FRG formalism lies the scale dependent effective action (Γ), which, throughout this paper will be considered in the local potential approximation (LPA).That is, Γ only contains a standard kinetic term besides the effective potential, V .Denoting the scale variable by k, the latter obeys the Wetterich equation [44,45], which, in d spatial dimensions, using Litim's regularization [46] takes the form of where T is the temperature, Ω d = dΩ d /(2π) d comes from the angular integral, ω n = 2πnT are bosonic Matsubara frequencies, V ′′ k is the mass matrix, meaning the second derivative of V k with respect to the dynamical variables (in our case, i.e., for N f = 3 flavors its size is 18 × 18), and ∂k is a differential operator acting by definition only on the explicit k-dependence.Note that, the physical meaning of V k is an effective potential, where only fluctuations with momenta q > k are taken into account.Practically ( 8) is solved by starting from a typical hadronization scale Λ [47], where V k→Λ is considered to be the fluctuationless classical potential, and we integrate down toward k → 0. We also note that ( 8) is derived with an infrared regulator that only cuts off fluctuations in the spatial directions, thus the Matsubara sum runs over all possible frequencies.
Since we are working with the ansatz of ( 6), we need to transform the right-hand side (rhs) of ( 8) such that it becomes compatible with (6).We follow the same procedure developed in [38,39].The most important observation is that in order to derive flow equations for U k (ρ, ∆) and C k (ρ), one is allowed to work in any background field of ⟨M ⟩ that allows for a unique reconstruction of the dependence on the invariants, not necessarily the one, which will be realized physically in the ground state.
We start by calculating the flow of U k (ρ, ∆).Since for any ⟨M ⟩ ∼ 1 , τ = 0, in such a background the lefthand side (lhs) of ( 8) equals the flow of U k (ρ, ∆).(Note that, the derivatives of τ are, in principle, nonzero in this background.)We, therefore, choose ⟨M ⟩ = (s 0 + iπ 0 )T 0 , T 0 = 1/ √ 61 .In this case, the invariants become ρ = 6).Note that, this background respects the SU (3) × SU (3) subgroup of chiral symmetry, therefore, the mass matrix V ′′ k has 8+8 degenerate eigenmodes in the "planes" (s i , π i )(i = 1, ..., 8), and a doublet in the (s 0 , π 0 ) sector.It is important to note that throughout the calculation each matrix element depends on the background components of s 0 , π 0 , but in the 8+8 degenerate eigenmodes they naturally combine into the ρ and ∆ invariants, see their respective expressions above.This is not true for the doublet, but after calculating the corresponding 2 × 2 determinant, its contribution to the flow equation happens to depend again only on the invariant combinations ρ and ∆, as it should.Summing up all terms yields the flow for U k (ρ, ∆), see (A1) in the Appendix.
For determining the scale evolution of the coefficient function of the τ invariant, that is C k (ρ), the purely imaginary background ⟨M ⟩ = i(π 0 T 0 + π 8 T 8 ) is particularly convenient.In this background the cubic invariant ∆ vanishes, but ρ = 1 2 (π 2 0 +π 2 8 ) and τ = That is, the lhs of (8) now becomes and since we already obtained ∂ k U k , by calculating the rhs of (8) and subtracting the ∂ k U k part, we obtain the flow of C k (ρ).In accordance with [39], the mass matrix breaks up into three degenerate {σ i , π i } doublets (i = 1, 2, 3), plus another four degenerate doublets (i = 4, 5, 6, 7), with a fully coupled quartet in the subspace {s 0 , s 8 , π 0 , π 8 }.One way to proceed with the calculations is that one inverts the ρ and τ equations and expresses π 0 , π 8 in terms of the latter.Then all terms in the rhs of the flow equation can be expressed with the help of the invariants ρ and τ , but at first sight they turn out to be non analytic in the latter, i.e. they contain terms ∼ √ τ .When expanding the contribution of each sector in terms of τ , the piece that comes from the {s 0 , s 8 , π 0 , π 8 }-sector completes the two sets of degenerate doublets into an analytic function of both ρ and τ , meaning that all the peculiar ∼ √ τ terms vanish.After a long and tedious calculation, the first nontrivial order in τ provides the flow for C k (ρ), see (A2) in the Appendix.

A. Numerics
We do not intend to solve the flow equation for U k (ρ, ∆) in its full generality, mainly due to the high numerical cost.In what follows we are interested in the scale evolution of the |Q| = 1, 2 instanton contributions.By introducing the notations U (ρ) := U (0) (ρ), A(ρ) := U (1) (ρ), B(ρ) := U (2) (ρ), using the decomposition of (7), we consider the effective potential as The flow equation for C k (ρ) is given by (A2), while those for U k (ρ), A k (ρ) and B k (ρ) can be obtained by the zeroth, first, and second order terms in a ∆ expansion of (A1).The choice of ( 9) allows us to investigate the importance of the |Q| = 2 interaction [i.e.B k (ρ)] in light of the one with The coupled differential equations are solved for three spatial dimensions (d = 3) using the grid method.All grids are set up in ρ space with the spacing δρ = 50 MeV 2 .Field derivatives are calculated using the six-point formula, except those close to the boundaries, where the five-and four-point formulas are used.Initial conditions are chosen to correspond to k = Λ ≡ 1 GeV , where, based on perturbative renormalizability, we assume that with some real UV parameters m 2 , λ 1 , λ 2 , a.These initial values at k = Λ are determined by the requirement of optimally reproducing the pseudoscalar meson spectra once all fluctuations are taken into account, see Table I.Note that, the initial condition for the B k (ρ) coefficient function is set to be zero, which is due to the (perturbative) irrelevance of the ∆ 2 operator at high scales.We have checked the sensitivity of the results with respect to the variation of the UV value of the latter interaction in the natural range of O([0.1 − 10]/Λ 2 ), and we have found no significant difference between the final results, as expected.Note that, in the physical point one needs to subtract an explicit symmetry breaking term from the potential, i.e. h 0 T 0 + h 8 T 8 , where h 0 and h 8 are determined using the partially conserved axialvector current relations (PCAC), see details in [38,39].As discussed earlier, without explicit breaking, in the vacuum ⟨M ⟩ ∼ T 0 ∼ 1 , but by introducing h 0 and h 8 into the potential, the minimum point of V shifts, and the actual physical background becomes ⟨M ⟩ = v 0 T 0 + v 8 T 8 .
The differential equations are integrated via the fourthorder Runge-Kutta method.We typically expect a critical slowing down in the numerics when approaching k → 0, therefore, we stop the flows at k/Λ = 0.1.At this point we avoid numerical instabilities, while all functions are practically converged and the scale dependence is rather moderate.We also note that all Matsubara sums are performed numerically, with a cutoff of n max = 10 4 .

B. Ground state and anomaly
All results are split into two parts.On the one hand, we treat the initial a parameter (which corresponds to the bare parameter in the perturbative renormalization process) as a constant, but on the other hand, we are also interested in a scenario, where it explicitly depends on the temperature.The latter is motivated by the fact that since we are dealing with an effective theory cut off at the rather low Λ = 1 GeV scale, the initial parameters should inherit non-negligible environment dependence when deriving them from the underlying theory of QCD.Denoting by T c the pseudocritical temperature of the chiral transition, it is known that in the dilute instanton gas approximation, valid for T ≫ T c , the instanton density is exponentially suppressed, leading to recovery of the anomalously broken axial U A (1) symmetry.When lowering the temperature, however, an instanton liquid model is more appropriate, and it is sometimes argued that the instanton density is weakly dependent on the temperature below T c [48].Therefore, as a working hypothesis, and motivated by obtaining the correct large T behavior of the anomaly, we make the initial anomaly parameter of the meson model temperature dependent as This choice activates the high-temperature behavior of the anomaly beyond the critical temperature and maintains it constant below the critical temperature.For T > T c , the a ∼ T −5 behavior originates from integrating over the instanton size after applying the semi classical approximation for the instanton density 2 [40].
Note that, the remaining m 2 , λ 1 , λ 2 parameters can, in principle, also carry explicit T -dependence, but for the sake of our main motivation of investigating the anomaly evolution, these parameters will be treated as temperature independent quantities in the current study.In Fig. 1 and Fig. 2 we show the thermal evolution of the ground state (i.e. the minimum point of the effective potential V ) in terms of the nonstrange (v ns ) and strange (v s ) condensates for both scenarios described above. 3 The prediction for the pseudocritical temperature, defined as the inflection point of the v ns (T ) function, is T c ≈ 160 MeV , being very close to lattice simulations 2 Note that had one retained a nonzero initial value for B k=Λ ≡ b, for the same reasons as described in the main text, one would have had b(T ) = b(T = 0) 1 + (Tc/T ) 14 − 1 Θ(T − Tc) . 3The nonstrange and strange components are defined through the [9, 49].The break point in case of the second scenario has no physical meaning, it is due to the fact that a(T ) is a nondifferentiable function at T = T c .The choice we made for the activation point of the T -dependence is somewhat arbitrary, and furthermore, had we chosen a smoother interpolation between the two regimes, we would have obtained a curve without breaking points.Note that while the nonstrange condensate shows a significant evaporation, the strange component loses less than 20% of its T = 0 value at T ≈ 2T c .Now let us focus on the fate of the axial anomaly as the function of the temperature.In Fig. 3 we show, for a temperature independent initial a anomaly parameter, the absolute value of the dimensionless, Ā(ρ) = A(ρ)/k, B(ρ) = B(ρ)k 2 anomaly coefficient functions (note that both A and B are negative).First, we note that, the shape of each function does not depend very strongly on the temperature, but the actual anomaly strength does, as the evaporation of v ns and v s yields the chirally invariant combination ρ min to show a decreasing tendency with T .This is in agreement with our earlier findings [38,39] showing that mesonic fluctuation effects strengthens the anomaly.Now we see that on top of the |Q| = 1 sector, the same applies for interactions coming from instantons with |Q| = 2 topological charge.Second, Fig. 3 also tells us that in order to reproduce the correct high-T behavior of the axial anomaly, we absolutely must include the damping effects through the explicit temperature dependence of the initial anomaly coefficient, as already announced in the previous paragraph.This is illustrated FIG. 6. Ratio between the chirally symmetric and anomalous parts of the effective potential as a function of the temperature.
in Fig. 4, where the effective anomaly strength, is plotted against the temperature.The "min" subscripts refer to the actual minimum point of the complete effective potential, that is where the effective interaction is defined.We see that by neglecting the explicit temperature dependence of the initial a anomaly parameter, A eff would monotonically increase with T even beyond T c .To avoid such a scenario did we employ the assumption (11).
As a result, we still see an intermediate strengthening of the U A (1) anomaly toward T → T c , however, including the |Q| = 2 interactions somewhat moderates the effect compared to our earlier findings, which included only the |Q| = 1 sector [39].
One may now wonder about the importance of the |Q| = 2 interactions, and this is what is shown in Fig. 5.We plot the ratio between the |Q| = 2 contribution to the effective anomaly coupling and the latter itself.Interestingly, for this quantity both discussed scenarios show quite a rapid decrease (obviously the one that corresponds to an initially temperature dependent anomaly coupling is faster), but the most important observation is that for T ≲ T c the |Q| = 2 topological sector provides around 10% of the anomaly strength.This confirms that for T ≲ T c , instanton interactions with |Q| = 2 topological winding are definitely important, and we also see, however, that, for T ≳ T c it is safe to keep only the |Q| = 1 sector, as expected from a dilute instanton gas approximation in the underlying theory.Finally, we define as the absolute value of the ratio between the chirally symmetric and anomalous contributions in the effective potential 4 .Interestingly, as it can be seen in Fig. 6, for lower temperatures they are almost of the same order, but for T ≳ T c , the anomalous contributions rapidly disappear, much faster than the effective anomaly coupling itself.This is easily understood, since the numerator of r E contains an extra factor of v 2 ns evaporating fast at high T , compared to A eff .This comes from the cubic invariant being ∆ ∼ v s v 2 ns .Before moving on to the mass spectra we note that in [43] it is not entirely ruled out that in the chiral limit, at the critical point all anomaly couplings may simultaneously vanish.In the light of the present results this would mean a rather radical change in the behavior of the topological fluctuations compared to those at the physical point.

C. Mass spectrum
The mass spectrum is obtained by diagonalizing the second derivative tensor of V k=0 [see (9)] at the minimum point of the complete effective potential, computed in the background of the physical condensates v 0 (T ), v 8 (T ) [or, equivalently, v ns (T ) and v s (T )].For this we have to obtain numerically the first and second derivatives of the U (ρ), C(ρ), A(ρ), and B(ρ) functions, and make use of the derivatives of the ρ, τ and ∆ invariants with respect to all field variables.Analytical formulas for the latter can be found in [39].
The T -dependence of the spectrum for our two scenarios can be seen in Fig. 7 and Fig. 8, respectively.One notices the different high temperature behavior for each case.In Fig. 7, at high T the anomaly coupling 4 We note that U (ρ) is normalized such that U (0) = 0.A eff stays nonzero, though among the condensates, v ns practically disappears.The strange condensate v s is still significant.As for Fig. 8, A eff is damped according to (11) and it practically vanishes together with v ns , but v s is still alive.In what follows, we provide a semi quantitative analysis of the observed high temperature spectra based on these remarks.The basic multiplet structure is the same for both cases in the two component (v ns , v s ) background.We have both in the pseudoscalar and scalar sectors one triplet (π and a 0 ), one quartet (K and K * 0 ), and the mixing (η, η ′ ) pair in the pseudoscalar, and (σ, f 0 ) in the scalar sector, respectively.We do not list here exact but lengthy mass formulas separately, instead restrict ourselves on the presentation of useful relations, which can clarify the behaviors observed in Fig. 7 and Fig. 8, reflecting especially on the survival/suppression of the anomaly coupling.For the sake of transparency and simplicity, in the forthcoming analysis we neglect the ρ dependence of the A(ρ), B(ρ) and C(ρ) coefficient functions, and restrict ourselves to U (ρ) = m 2 ρ + λ 1 ρ 2 .Still, one can nicely in-terpret some of the characteristic features of the observed high-T spectra.
First, we have two compact expressions for the scalarpseudoscalar difference in the triplet and the quartet sectors: As we saw already, near and above the pseudocritical temperature, the nonstrange condensate quickly diminishes, while v s stays nearly constant.This means that in the triplet sector a nonzero mass difference in the high temperature region signals the persistence of the anomaly.
Below, in the high temperature region we set v ns ≈ 0, and then even the expressions of the actual masses become rather simple (see also [50]): The D discriminant constructed from the eigenvalue equation of the mixing matrix in the pseudoscalar (0 − 8) sector determines the η − η ′ mass difference, which, for a vanishing nonstrange condensate is rather transparent: This difference and also the above mass formulas can be discussed assuming two limiting situations, depending on the relative magnitude of the mass contribution from the nonanomalous and anomalous parts of the potential.When Cv 2 s ≫ −Av s , which is certainly true when A is suppressed by its explicit T -dependence, one finds asymptotically This means that in this case the levels of η ′ group separately.This can be clearly observed in Fig. 8.
On the other hand, when the anomaly survives and dominates, i.e. −Av s ≫ Cv 2 s , one finds In this case the lighter excitation of the pseudoscalar (0− 8) sector, i.e., η, becomes nearly degenerate with K and K * 0 , while the heavier η ′ will be located close to a 0 .Note that, the anomaly sustains a mass difference between π and a 0 .This characteristics is observed in Fig. 7, where one assumed no explicit temperature dependence for the initial anomaly coefficient a.
Based on the above analysis, by observing the variation of the meson masses as the temperature is gradually increased, such rearrangement of η and η ′ would clearly indicate the onset of anomaly suppression.We also note that the current setup shows that the drop in the η ′ mass at T c becomes larger compared to our earlier results, hinting that if this is maintained at the partial recovery of chiral symmetry at the nuclear liquid-gas transition, the η ′ -nucleon bound state might indeed be observable in a nuclear medium [51,52].
Also, we note that a similar analysis in the scalar (0−8) sector shows that, the σ and f 0 masses in the v ns ≈ 0 limit become Both scenarios Cv 2 s ≫ −Av s and −Av s ≫ Cv 2 s show that σ degenerates with π at high T , and it is also revealed that the f 0 mass is insensitive to the anomaly even when the latter dominates.
Finally we comment on the mixing angles in the scalar and pseudoscalar sectors.Note that, at first, one usually defines the scalar and pseudoscalar mixing angles, Θ s/π , in the 0 − 8 basis, which are then transformed into φ s/π , corresponding to the ns − s coordinate system, through the ideal mixing matrix (see, e.g., the Appendix of [19]).Depending on the orientation (i.e. the sign) of Θ s/π , one ends up with curves for the temperature dependence of φ s/π that are similar to either [18] or [19].Since our model and method are closer to that of [19], we adopt their convention, which, in the pseudoscalar sector leads to The second formula is especially convenient for calculating the angle, as it is not sensitive to either its orientation, or a possible level crossing between M 2 π,00 and M 2 π,88 .As the numerics show that M 2 π,08 < 0, in the above convention both Θ π > 0 and φ π > 0. Our results are very similar to that of Fig. 8 in [19].We see an almost completely identical behavior for φ s , but more interestingly, temperature dependence of φ π via our a(T ) ansatz is very close to that of the LPA'+Y approximation of [19], showing that φ π goes to zero at high temperature.It seems that the effect of the wave function renormalization is similar to a temperature dependent bare anomaly parameter, introduced in (11).

V. SUMMARY
In this paper we analyzed the thermal behavior of the axial anomaly using the functional renormalization group method.Compared to our earlier study [39], we applied an enlarged set of terms in the effective potential of the N f = 3 meson model.We derived scale evolution equations for the effective potential that include interactions originated from instantons of the underlying theory of QCD with arbitrary topological charge [see Eq. (A1)].We discussed two distinct scenarios: (1) the initial |Q| = 1 anomaly coupling is an environment independent constant, and (2) the initial |Q| = 1 anomaly coupling inherits explicit temperature dependence from QCD interactions occurring above the initial momentum scale.
We numerically solved the flow equations for the chiral condensate dependent anomaly couplings coming from instantons with |Q| = 1, 2 topological charges and found that for low temperatures the |Q| = 2 interaction contributes around 10% to the effective U A (1) coupling A eff .The ratio between the |Q| = 2 and |Q| = 1 contributions to A eff rapidly decreases once the temperature hits T c , showing that at high temperatures it is sufficient to include instanton interactions with |Q| = 1.For T ≲ T c , however, |Q| = 2 interactions cannot be dropped.
In accordance with our earlier findings we saw that mesonic fluctuations tend to strengthen the anomaly couplings toward the critical point, which are then suppressed for T ≳ T c by instanton effects.We also calculated numerically the mass spectrum and provided semiquantitative analytical formulas to understand its high temperature behavior.In particular, we showed that the η and η ′ mass difference, on top of the anomalous terms, also contains one of the quartic couplings through the strange condensate.As a result, we argued that the η and η ′ clustering with other excitations occur differently for weak (temperature suppressed) and strong (temperature sustained) anomaly.Observing the finite tempera-ture spectrum in lattice simulations may, therefore, hint on the actual thermal behavior of the U A (1) anomaly coupling.
Our analysis could be extended in a straightforward fashion for couplings that correspond to |Q| > 2 instanton interactions by expanding (A1) to higher order in ∆.Our method also offers the possibility of treating all anomalous interactions on an equal footing by solving the flow equation for U (ρ, ∆) on a two dimensional grid.Though numerically this might be very demanding, it would certainly provide the most general answer to the question of the importance of all anomalous interactions.Furthermore, it could be interesting to apply our method to two-color QCD with N f = 2, where the global symmetry of the system turns into SU (4) (Pauli-Gursey symmetry) [31].The technique presented here, equipped with the corresponding lattice data could provide new insight into the behavior of the anomaly in such a system.

ACKNOWLEDGMENTS
This research was supported by the Hungarian National Research, Development, and Innovation Fund under Project No. FK142594.G.F. also received support from the Hungarian State Eötvös Scholarship.The authors would like to express their gratitude to Naoki Yamamoto for his insights regarding the explicit temperature dependence of the initial anomaly parameter in the ultraviolet.

APPENDIX: FLOW EQUATIONS
Here we list the explicit flow equations that were solved numerically.The flow equation for U k (ρ, ∆), determined in Sec.III, using the ⟨M ⟩ = (s 0 + iπ 0 )T 0 background is the following: where the subscripts ρ and/or ∆ refer to partial derivatives with respect to the corresponding variable.The flow equation for C k (ρ), which was calculated in the purely imaginary ⟨M ⟩ = i(π 0 T 0 + π 8 T 8 ) background, takes the form of

FIG. 1 .FIG. 2 .
FIG.1.Thermal behavior of the nonstrange and strange condensates with a temperature independent initial anomaly parameter.Prediction for the pseudocritical temperature is Tc ≈ 160 MeV .

FIG. 3 .FIG. 4 .
FIG.3.Dimensionless anomaly coefficient functions at T = 0 and T = Tc, as a function of the chirally invariant combination ρ, for a T -independent initial anomaly parameter.Both functions monotonically decrease showing that as the condensates evaporate, fluctuations tend to strengthen the axial anomaly.The scale is set to k/Λ = 0.1.

FIG. 5 .
FIG. 5. Importance of the |Q| = 2 contribution to the effective anomaly coupling for both scenarios.Ratio of the |Q| = 2 contribution to the full effective anomaly coupling is plotted.