$\eta'$ and $\eta$ mesons at high T when the U_A(1) and chiral symmetry breaking are tied

The approach to the eta'-eta complex employing chirally well-behaved quark-antiquark bound states and incorporating the non-Abelian axial anomaly of QCD through the generalization of the Witten-Veneziano relation, is extended to finite temperatures. Employing the chiral condensate has led to a sharp chiral and U_A(1) symmetry restoration, but with the condensates of quarks with realistic explicit chiral symmetry breaking, which exhibit a smooth, crossover chiral symmetry restoration in qualitative agreement with lattice QCD results, we get a crossover U_A(1) transition, with smooth and gradual melting of anomalous mass contributions. This way we obtain a substantial drop of the eta' mass around the chiral transition temperature, but no eta mass drop. This is consistent with the present empirical evidence.


I. INTRODUCTION
The experiments at heavy-ion collider facilities-such as RHIC, LHC, FAIR, and NICA-aim to produce a new form of hot and/or dense QCD matter [1,2]. Clear signatures of its production are thus very much needed. The most compelling such signal would be a change in the pertinent symmetries, i.e., the restoration (in hot and/or dense matter) of the symmetries of the QCD Lagrangian which are broken in the vacuum, notably the [SU A ðN f Þ flavor] chiral symmetry for N f ¼ 3 ¼ 2 þ 1 light quark flavors q, and the U A ð1Þ symmetry. This provides much motivation to establish that experiment indeed shows this, as well as to give theoretical explanations of such phenomena.
The first signs of a (partial) restoration of the U A ð1Þ symmetry were claimed to be seen in 200 GeV Au þ Au collisions [3,4] at RHIC by Csörgő et al. [5]. They analyzed the η 0 -meson data of the PHENIX [3] and STAR [4] collaborations through several models for hadron multiplicities, and found that the η 0 mass (M η 0 ¼ 957.8 MeV in vacuum) decreases by at least 200 MeV inside the fireball. The vacuum η 0 is, comparatively, so very massive since it is predominantly the SU V ðN f Þ-flavor singlet state η 0 . Its mass M η 0 receives a sizable anomalous contribution ΔM η 0 due to the U A ð1Þ symmetry violation by the non-Abelian axial Adler-Bell-Jackiw anomaly ["gluon anomaly," or "U A ð1Þ anomaly" for short], which makes the divergence of the singlet axial quark currentqγ μ γ 5 1 2 λ 0 q nonvanishing even in the chiral limit of vanishing current masses of quarks, m q → 0. This mass decrease is then a sign of a partial U A ð1Þ symmetry restoration in the sense of a diminishing contribution of the U A ð1Þ anomaly to the η 0 mass, which would decrease to a value readily understood in the same way [6] as the masses of the octet of the light pseudoscalar mesons P ¼ π 0;AE , K 0;AE ,K 0 , η, which are exceptionally light almost-Goldstone bosons of dynamical chiral symmetry breaking (DChSB).
A recent experimental paper studied 200 GeV Au þ Au collisions [7]. Although a new analysis of the limits on the η 0 and η masses was beyond the scope of Ref. [7], the data contained therein make it possible, and preliminary considerations [8] confirm the findings of Ref. [5].
The first explanation [9] of these original findings [5] was offered by conjecturing that the Yang-Mills (YM) topological susceptibility, which leads to the anomalously high η 0 mass, should be viewed through the Leutwyler-Smilga (LS) [10] relation (12). This ultimately implies that the anomalous part of the η 0 mass decreases together with the quark-antiquark (qq) chiral-limit condensate hqqi 0 ðTÞ as the temperature T grows towards the chiral restoration temperature T Ch and beyond. This connection between the U A ð1Þ symmetry restoration and the chiral symmetry restoration was just a conjecture until our more recent paper [11] strengthened the support for this scenario. Nevertheless, there was also a weakness: our approach predicted the decrease of not only the η 0 mass, but also an even more drastic decrease of the η mass M η , and signs for that have not been seen in any currently available data [7,12]. In the present paper, we show that the predicted Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP 3 . decrease of M η [9] was the consequence of employing the chiral-limit condensate hqqi 0 ðTÞ, since it decreases too fast with T after approaching T ∼ T Ch . We then perform T > 0 calculations in the framework of the more recent work by Benić et al. [11], where the LS relation (12) is replaced by the full-QCD topological charge parameter (18) [13][14][15]. There, one can employ qq condensates for realistically massive u, d, and s quarks, with a much smoother T dependence. As a result, the description of the η-η 0 complex of Ref. [9] is significantly improved, since our new T dependences of the pseudoscalar meson masses do not exhibit a decrease of the η mass, while a considerable decrease of the η 0 mass still exists, which is consistent with the empirical findings [5].
The light pseudoscalar mesons are simultaneously qq 0 bound states (q, q 0 ¼ u, d, s) and (almost-)Goldstone bosons of the DChSB of nonperturbative QCD. We can implement both simultaneously by using the Dyson-Schwinger (DS) equations as Green functions of QCD (see, e.g., Refs. [16][17][18][19] for reviews). Particularly pertinent are the gap equation for dressed quark propagators S q ðpÞ with DChSB-generated self-energies Σ q ðpÞ, (while S free q are free ones), and the Bethe-Salpeter equation (BSE) for the qq 0 meson bound-state vertices Γ qq 0 , where K is the interaction kernel, and e, f, g, h represent (schematically) the collective spinor, color, and flavor indices. This nonperturbative and covariant bound-state DS approach can be applied for various degrees of truncations, assumptions, and approximations, ranging from ab initio QCD calculations and sophisticated truncations (see, e.g., Refs. [16][17][18][19][20][21][22] and references therein) to very simplified modeling of hadron phenomenology, such as utilizing Nambu-Jona-Lasinio point interactions. For applications in involved contexts such as nonzero temperature or density, strong simplifications are especially needed for tractability. This is why the separable approximation [23] is adopted in this paper [see the discussion between Eqs. (4) and (5)]. However, when describing pseudoscalar mesons (including η and η 0 ) reproducing the correct chiral behavior of QCD is much more important than the dynamicsdependent details of their internal bound-state structure.
A rarity among bound-state approaches, the DS approach can also achieve the correct QCD chiral behavior regardless of the details of model dynamics, but under the condition of a consistent truncation of DS equations, respecting pertinent Ward-Takahashi identities [16][17][18][19]. A consistent DS truncation, where DChSB is very well understood, is the rainbow-ladder approximation (RLA). Since it also enables tractable calculations, it is still the most used approximation in phenomenological applications, and we also adopt it here. In the RLA, the BSE (2) employs the dressed quark propagator solution SðpÞ from the gap equation (1) and (4), which in turn employs the same effective interaction kernel as the BSE. It has a simple gluon-exchange form, where both quark-gluon vertices are bare, so that the quark self-energy in the gap equation is where D ab μν ðkÞ eff is an effective gluon propagator. These simplifications should be compensated by modeling the effective gluon propagator D ab μν ðkÞ eff in order to reproduce well the relevant phenomenology; here, pseudoscalar (P) meson masses M P , decay constants f P , and condensates hqqi, including T-dependence of all these. In the present paper, we use the same model as in Ref. [9] and attempt to improve their approach to the T dependence of the U A ð1Þ anomaly. All of the details on the functional form and parameters of this model interaction can be found in the Sec. II A of Ref. [24]. Such models-socalled rank-2 separable models-are phenomenologically successful (see, e.g., Refs. [23][24][25][26][27]). However, they have the well-known drawback of predicting a somewhat too low transition temperature: the model we use in this paper and that was used in Refs. [9,24,26,27] has T Ch ¼ 128 MeV, i.e., some 17% below the now widely accepted central value of 154 AE 9 MeV [28][29][30]. But, rather than quantitative predictions at specific absolute temperatures, we are interested in the relative connection between the chiral restoration temperature T Ch and the temperature scales characterizing signs of the effective disappearance of the U A ð1Þ anomaly, for which the present model is adequate. In addition, Ref. [31] showed that coupling to the Polyakov loop can increase T Ch , while the qualitative features of the T dependence of the model are preserved. Thus, separable model results at T > 0 are most meaningfully presented as functions of the relative temperature T=T Ch , as in Refs. [9,24].
Anyway, regardless of the details of the model dynamics [i.e., the choice of D ab μν ðkÞ eff ] and thanks to the consistent truncation of DS equations, the BSE (2) yields the masses M qq 0 of pseudoscalar P ∼ qq 0 mesons which satisfy the Gell-Mann-Oakes-Renner-type relation with the current masses m q , m q 0 of the corresponding quarks: While this guarantees that all M qq 0 → 0 in the chiral limit, it also shows that the RLA cannot lead to any U A ð1Þanomalous contribution responsible for ΔM η 0 . That is, the RLA gives us only the nonanomalous partM 2 NA of the squared-mass matrixM 2 ¼M 2 NA þM 2 A of the hidden-flavor (q ¼ q 0 ) light (q ¼ u, d, s) pseudoscalar mesons. In the basis fuū; dd; ssg,M 2 NA is simplyM 2

The anomalous partM 2
A arises because the pseudoscalar hidden-flavor states qq are not protected from the flavormixing QCD transitions (through anomaly-dominated pseudoscalar gluonic intermediate states), as depicted in Fig. 1. They are obviously beyond the reach of the RLA and horrendously hard to calculate. Nevertheless, they cannot be neglected, as can be seen in the Witten-Veneziano relation (WVR) [32,33], which remarkably relates the full-QCD quantities (η 0 , η, the K-meson masses M η 0 ;η;K , and the pion decay constant f π ) to the topological susceptibility χ YM of the (pure-gauge) YM theory: Its chiral-limit-nonvanishing rhs is large (roughly 0.8 to 0.9 GeV 2 ), while Eq. (5) basically leads to the cancellation of all chiral-limit-vanishing contributions on the lhs [9]. The rhs is the WVR result for the total mass contribution of the U A ð1Þ anomaly to the η-η 0 complex, M U A ð1Þ . TheM 2 A matrix elements generated by the U A ð1Þanomaly-dominated transitions qq → q 0q0 (see Fig. 1) can be written [35] in the flavor basis fuū; dd; ssg as Here b q ¼ ffiffi ffi β p for both q ¼ u, d, since we assume m u ¼ m d ≡ m l [i.e., isospin SUð2Þ symmetry] which is an excellent approximation for most purposes in hadronic physics. For example, M uū ¼ M dd ≡ M ll ¼ M ud ≡ M π obtained from the BSE (2) is our RLA model pion mass for π þ ðπ − Þ ¼ udðdūÞ and It still contains M ss , the mass of the unphysical (but theoretically very useful) ss pseudoscalar obtained in the RLA. However, thanks to Eq. (5), it can also be expressed through the masses of physical mesons, to a very good approximation [24,27,[34][35][36][37]. Its decay constant f ss is calculated in the same way as f π and f K .
Since the s quark is much heavier than the u and d quarks, in Eq. (7) we have b q ¼ X ffiffiffi β p for q ¼ s, with X < 1. Transitions to and from more massive s quarks are suppressed, and the quantity X expresses this influence of the SUð3Þ flavor symmetry breaking. The most common choice for the flavor-breaking parameter has been the estimate X ¼ f π =f ss [9,24,27,[34][35][36][37], but we found [11] that it necessarily arises in the variant of our approach relying on Shore's generalization of the WVR (6) [13,14] (see Sec. III).
The anomalous mass matrixM 2 A [which is of the pairing form (7) in the hidden-flavor basis fuū; dd; ssg] in the octet-singlet basis fπ 0 ; η 8 ; η 0 g of hidden-flavor pseudoscalars becomeŝ which shows that the SUð3Þ flavor breaking [X ≠ 1] is necessary for the anomalous contribution to the η 8 mass squared, ΔM 2 In the flavor SUð3Þsymmetric case (X ¼ 1), only the η 0 mass receives a U A ð1Þ-anomaly contribution: (8)] to be off diagonal, but in this basis the fη 8 ; η 0 g submatrix ofM 2 NA also gets strong, negative off-diagonal elements, g., Ref. [35]). Equation (8) thus shows that the interplay of the flavor symmetry breaking (X < 1) with the anomaly is necessary for the partial cancellation of the off-diagonal (8,0) elements in the complete mass matrixM 2 ¼M 2 NA þM 2 A , i.e., to obtain the physical isoscalars in a rough approximation as η ≈ η 8 and η 0 ≈ η 0 . How this changes with diminishing U A ð1Þ-anomaly contributions is exhibited in Secs. IV and V.
Since the isospin-limit π 0 decouples from the anomaly and mixing, only the isoscalar-subspace 2 × 2 mass matrix FIG. 1. Axial-anomaly-induced, flavor-mixing transitions from hidden-flavor pseudoscalar states P ¼ qq to P 0 ¼ q 0q0 including both possibilities q ¼ q 0 and q ≠ q 0 . All lines and vertices are dressed. The gray blob symbolizes all possible intermediate states enabling this transition. The three bold dots symbolize an even [34] but otherwise unlimited number of additional gluons. As pointed out in Ref. [34], the diamond graph is just the simplest example of such a contribution. M 2 needs to be considered. Even thoughM 2 is strongly off diagonal in the isoscalar basis fη NS ; η S g (the NS-S basis), in this basis it has the simple form which also shows that when the U A ð1Þ-anomaly contributions vanish (i.e., β → 0) the NS-S scenario is realized. This means that not only do the physical isoscalars become η → η NS and η 0 → η S , but also that their respective masses become M π and M ss . Our experience with various dynamical models (at T ¼ 0) shows [27,[34][35][36][37] that after pions and kaons are correctly described, a good determination of the anomalous mass shift parameter is sufficient for Eq. (10) to give good η 0 and η masses, since M 2 Nevertheless, calculating the anomalous contributions (∝ β) in DS approaches is a very difficult task. Reference [38] explored this by taking the calculation beyond the RLA, but they had to adopt extremely schematic model interactions (proportional to δ functions in momenta) for both the ladder-truncation part (3) and the anomaly-producing part. Another approach [39] obtained qualitative agreement with the lattice on χ YM (and, consequently, acceptable masses for η 0 and η) by assuming that the contributions to Fig. 1 are dominated by the simplest one-the diamond graph-if it is appropriately dressed (in particular, by an appropriately singular quark-gluon vertex).
However, we take a different route, since our goal is not to figure out how the breaking of U A ð1Þ comes about on a microscopic level, but rather to phenomenologically model and study the high-T behavior of the masses of the realistic η 0 and η, along with other light pseudoscalar mesons. In the DS context, the most suitable approach is then the one developed in Refs. [27,[34][35][36][37] and extended to T > 0 in Refs. [9,24].
The key is that the U A ð1Þ anomaly is suppressed in the limit of large number of QCD colors N c [32,33]. So, in the sense of the 1=N c expansion, it is a controlled approximation to view the anomaly contribution as a perturbation with respect to the (nonsuppressed) results obtained through the RLA (3)-(4). While considering meson masses, it is thus not necessary to look for anomalyinduced corrections to the RLA Bethe-Salpeter wave functions, 1 which are consistent with DChSB and with the chiral QCD behavior (5) that is essential for describing pions and kaons. The breaking of nonet symmetry by the U A ð1Þ anomaly can be introduced just at the level of the masses in the η 0 -η complex, by adding the anomalous contributionM 2 A to the RLA-calculatedM 2 NA . Its anomaly mass parameter β can be obtained by fitting [34] the empirical masses of η and η 0 or, preferably, from lattice results on the YM topological susceptibility χ YM (because then no new fitting parameters are introduced). Employing the WVR (6) yields [9,35] β ¼ β WV , while Shore's generalization gives (see Sec. III) β ¼ β Sho [11], where A is the QCD topological charge parameter, given below by Eq. (18) in terms of qq condensates of massive quarks, which turns out to be crucial for a realistic T dependence of the masses in the η 0 -η complex.

III. EXTENSION TO T ≥ 0
Extending our treatment [27,[34][35][36][37] of the η 0 − η complex to T > 0 is clearly more complicated. Since to the best of our knowledge there is no systematic derivation of the T > 0 version of either the WVR (6) or its generalization by Shore [13,14], it is tempting to try to straightforwardly replace all quantities by their T-dependent versions. In the WVR, these are the full-QCD quantities M η 0 ðTÞ, M η ðTÞ, M K ðTÞ, and f π ðTÞ, but also χ YM ðTÞ, which is a puregauge, YM quantity and thus much more resistant to high temperatures than QCD quantities that also contain quark degrees of freedom. Indeed, lattice calculations indicate that the decrease of χ YM ðTÞ (from which one would expect the decrease of the anomalous η 0 mass) only starts at a T some 100 MeV (or even more) above the (pseudo)critical temperature T Ch for the chiral symmetry restoration of full QCD, near where decay constants already decrease appreciably. It was then shown [24] that the straightforward extension of the T dependence of the YM susceptibility would even predict an increase of the η 0 mass around and beyond T Ch , contrary to experiment [5].
It could be expected that at high T, the original WVR (6) will not work since it relates the full-QCD quantities with a much more temperature-resistant YM quantity, χ YM ðTÞ. 1 It is instructive to recall [36,40] that nonet symmetry (or a broken version thereof) is in fact assumed (explicitly or implicitly) by all approaches using the simple hidden-flavor basis qq, e.g., to construct the SUð3Þ pseudoscalar meson states η 0 and η 8 without distinguishing between the qq states belonging to the singlet and those belonging to the octet. An independent a posteriori support for our approach is also that η and η 0 → γγ ðÃÞ processes are described well [34][35][36][37].
However, this problem can be eliminated [9] by using, at T ¼ 0, the (inverted) Leutwyler-Smilga (LS) relation [10] to express χ YM in the WVR (6) through the full-QCD topological susceptibility χ and the chiral-limit condensate hqqi 0 . Thus the zero-temperature WVR is retained, while the full-QCD quantities inχ do not have the T dependence mismatch with the rest of Eq. (6). Thus, instead of χ YM ðTÞ, Ref. [9] used the combinationχðTÞ [Eq. (12)] at T > 0, where the QCD topological susceptibility χ in the lightquark sector can be expressed as [10,15,41] This implies that the (partial) restoration of U A ð1Þ symmetry is strongly tied to the chiral symmetry restoration, since it is not χ YM ðTÞ but rather hqqi 0 ðTÞ [throughχðTÞ] that determines the T dependence of the anomalous parts of the masses in the η-η 0 complex [9]. The dotted curve in Fig. 2 illustrates how hqqi 0 ðTÞ decreases steeply to zero as T → T Ch , indicative of the second-order phase transition.
This behavior is followed closely byχðTÞ, and therefore also by the anomaly parameter β WV ðTÞ [Eq. (11)]. This makes the mass matrix (10) diagonal immediately after T ¼ T Ch , which marks the abrupt onset of the NS-S scenario M η 0 ðTÞ → M ss ðTÞ, M η ðTÞ → M π ðTÞ [9]. In Eq. (13), C m denotes corrections of higher orders in small m q , but it should not be neglected as C m ≠ 0 is needed to have a finite χ YM with Eqs. (12) and (13). They in turn give us the value C m at T ¼ 0 in terms of the qq condensate and the YM topological susceptibility χ YM . However, to the best of our knowledge, the functional form of C m is not known. Reference [9] thus tried various parametrizations covering reasonably possible T dependences of C m ðTÞ, but this did not greatly affect the results for the T dependence of the masses in the η 0 -η complex.
An alternative to the WVR (6) is its generalization by Shore [13,14]. There, relations containing the masses of the pseudoscalar nonet mesons take into account that η and η 0 should have two decay constants each [42]. If one chooses to use the η 8 -η 0 basis, they are f 8 η , f 8 η 0 , f 0 η , f 0 η 0 , and can be equivalently expressed through purely octet and singlet decay constants (f 8 , f 0 ) and two mixing angles (θ 8 , θ 0 ). This may seem better suited for use with effective meson Lagrangians than with qq 0 substructure calculations starting from the (flavor-broken) nonet symmetry, such as ours. Nevertheless, Shore's approach was also adapted for the latter bound-state context, and successfully applied there (in particular, to our DS approach in the RLA [27]). This was thanks to the simplifying scheme of Feldmann, Kroll, and Stech (FKS) [43,44]. They showed that this "two mixing angles for four decay constants" formulation in the NS-S basis, although in principle equivalent to the η 8 -η 0 basis formulation, can in practice be simplified further to a one-mixing-angle scheme using plausible approximations based on the Okubo-Zweig-Iizuka (OZI) rule. The decayconstant mixing angles in this basis are mutually close, ϕ S ≈ ϕ NS , and both are approximately equal to the state mixing angle ϕ rotating the NS-S basis states into the physical η and η 0 mesons, which diagonalizes the mass (squared) matrix (10). So, Ref. [27] numerically solved Shore's equations (combined with the FKS approximation scheme) for meson masses for several dynamical DS bound-state models [24,34,35]. Then, Ref. [11] presented analytic solutions thereof, for the masses of η and η 0 and the state NS-S mixing angle ϕ. These are rather long but closed-form expressions in terms of nonanomalous meson masses M π , M K and their decay constants f π , f K , as well as f NS and f S (the decay constants of the unphysical η NS and η S ), and, most notably, the full-QCD topological charge parameter A. This quantity (taken [13,14] from Di Vecchia and Veneziano [15]) plays the role of χ YM in the WVR in the mass relations of Shore's FIG. 2. The relative-temperature T=T Ch dependences of the pertinent order parameters calculated in our usual [9,24] separable interaction model. The odd man out is the (third root of the absolute value of the) chiral condensate hqqi 0 ðTÞ, which decreases steeply at T ¼ T Ch and dictates similar behavior [9] tõ χðTÞ. All of the other displayed quantities exhibit smooth, crossover behaviors, which are smoother for heavier flavors: the dash-dotted and dashed curves are the (third roots of the absolute values of the) condensates hssiðTÞ and hūuiðTÞ, respectively, the thin solid curve is the resulting topological susceptibility χðTÞ 1=4 , and the thick solid curve is the topological charge parameter AðTÞ 1=4 . The decay constants f π ðTÞ and f ss ðTÞ are, respectively, the lower dashed and dash-dotted curves. generalization. A will be considered in detail for the T > 0 extension, but now let us note that although Shore's generalization is in principle valid to all orders in 1=N c [13,14], Shore himself took advantage of and approximated A (as we shall at T ¼ 0) by the lattice result χ YM ¼ ð0.191 GeVÞ 4 [45]. Further, one should note that since the FKS scheme neglects OZI-violating contributions (that is, gluonium admixtures in η NS and η S ) it is consistent to treat them as pure qq states, accessible by our BSE (2) in the RLA. Then f NS ¼ f π , and f S ¼ f ss (the decay constant of the aforementioned "auxiliary" RLA ss pseudoscalar). We calculate its mass M ss with the BSE, but at T ¼ 0 it can also be related to the measurable pion and kaon masses, M 2 ss ≈ 2M 2 K − M 2 π , due to Eq. (5). Similarly, f ss can also be approximately expressed with these measurable quantities as f ss ≈ 2f K − f π . Thus, after taking A ≈ χ YM from lattice data, Ref. [11] calculated the η-η 0 complex using both the model-calculated and the empirical M π , M K , f π , and f K in their analytic solutions. This serves as a check (independently of any model) of the soundness of our approach at T ¼ 0.
Since the adopted DS model also enables the calculation of nonanomalous qq masses and decay constants for T > 0, the only thing still missing is the T dependence of the full-QCD topological charge parameter A, as χ YM ðTÞ is inadequate. But, A is used to express the QCD susceptibility χ through the "massive" condensates hūui, hddi, and hssi, i.e., away from the chiral limit, in contrast to Eqs. (12) and (13) [see, e.g., Eq. (2.12) in Ref. [13]]. Its inverse (expressing A) thus also contains the qq condensates out of the chiral limit for all light flavors q ¼ u, d, s, and so should χ in Eq. (18). That is, the light-quark expression for the QCD topological susceptibility in the context of Shore's approach should be expressed in terms of the current masses m q multiplied by their respective condensates hqqi realistically out of the chiral limit: As before [9], the small and necessarily negative correction term C m is found by assuming A ¼ χ YM at T ¼ 0. This large-N c approximation also easily recovers the LS relation (12): by approximating the realistically massive condensates with hqqi 0 everywhere in Eq. (18), the QCD topological charge parameter A reduces toχ, justifying the conjecture of Ref. [9] that connects the U A ð1Þ symmetry restoration with the chiral symmetry restoration. This connection between the two symmetries is still present. However, with the massive condensates we also get a more realistic, crossover T dependence of the masses, depicted in Figs. 3 and 4, and presented in Sec. IV. Figures 3 and 4 correspond to two variations of the unknown T dependence C m ðTÞ of the correction term in Eq. (19). As in Ref. [9], the simplest ansatz is a constant, C m ðTÞ ¼ C m ð0Þ, which is most reasonable for T < T Ch , where the condensates [and thus also the leading term in FIG. 3. T dependence, relative to T Ch , of various η 0 -η complex masses described in the text, the π mass (thick, dash-dotted curve) for reference, the halved (to maintain clarity) total U A ð1Þanomaly-induced mass 1 2 M U A ð1Þ (short-dashed curve), and the topological charge parameter A 1=4 (solid curve). The straight line is 2 times the lowest fermion Matsubara frequency 2πT. χðTÞ] change little. But above some higher T, the negative C m ð0Þ-although initially much smaller in magnitude than the leading term-will make χðTÞ [and therefore also AðTÞ] change sign. Concretely, this limiting T above which there is no meaningful description is found a little above 1.6T Ch .
For another, nonconstant C m ðTÞ that would not have such a limiting temperature, we now have a lead from lattice data where the high-T asymptotic behavior of the QCD topological susceptibility has been found to be a power law, χðTÞ ∝ T −b [46,47]. The high-T dependence of our model-calculated condensates is also (without fitting) such that the leading term of our χðTÞ in Eq. (19) has a similar power-law behavior, with b ¼ 5.17. Also, the values of our leading terms are, qualitatively, for all T roughly in the same ballpark as the lattice results [46,47]. We thus fit the quickly decreasing power-law C m ðTÞ for high T by requiring that (i) this more or less rough consistency with lattice χðTÞ values is preserved, (ii) the whole χðTÞ has the high-T power-law dependence as the leading term (with b ¼ 5.17), and (iii) C m ðTÞ joins smoothly with the low-T value C m ð0Þ determined from χ YM at T ¼ 0.
Our nonconstant choice of C m ðTÞ yields the masses in Fig. 3 [and χðTÞ and AðTÞ in Fig. 2], but these results are very similar to the ones with C m ðTÞ ¼ C m ð0Þ (of course, only up to the limiting T a little above 1.6T Ch ) in Fig. 4. Thus, Fig. 4 uses a different scale than Fig. 3, i.e., only the mass interval between 0.55 and 1.05 GeV, so as to zoom in on the η-η 0 complex and better discern its various overlapping curves, including M U A ð1Þ ðTÞ.
The second choice of C m ðTÞ enables in principle the calculation of χðTÞ and AðTÞ without any limiting T. Nevertheless, Fig. 3 does not reach higher than T ¼1.8T Ch , because the model chosen for the RLA part of our calculations seems to become unreliable at higher T's: the mass eigenvalues seem increasingly too high, since they tend to cross the sum of the lowest q þq Matsubara frequencies. Fortunately, by T=T Ch ¼ 1.8 the asymptotic scenario for the anomaly has been reached, as we explain in the next section where we give a detailed description of all pertinent results at T ≥ 0. Figure 2 shows how various magnitudes of current-quark masses m q influence the T dependence and size of qq condensates hqqi and pseudoscalar decay constants f qq calculated in our adopted model. Defined, e.g., in Sec. II A of Ref. [24], it employs the parameter values m u ¼ m d ≡ m l ¼ 5.49 MeV and m s ¼ 115 MeV.

IV. RESULTS AT T ≥ 0 IN DETAIL
For both condensates and decay constants, larger current-quark masses lead to larger "initial" (i.e., T ¼ 0) magnitudes and, what is even more important for the present work, to smoother and slower falloffs with T. The magnitude of (the third root of) the strange-quark condensate is the top dash-dotted curve in Fig. 2. Its T ¼ 0 value jhssij 1=3 ¼ 238.81 MeV remains almost unchanged until T ¼ T Ch , and falls below 200 MeV (i.e., by some 20%) only for T ≈ 1.5T Ch . On the other hand, the T ¼ 0 value of the isosymmetric condensates of the lightest flavors, hūui ¼ hddi ≡ hlli ¼ ð−218.69 MeVÞ 3 , is quite close to the chiral one, hqqi 0 ¼ ð−216.25 MeVÞ 3 , showing how well the chiral limit works for u and d flavors in this respect. Still, the small current masses of u and d quarks are sufficient to lead to a very different T dependence of the lightest condensates, depicted by the dashed curve. It exhibits a typical smooth crossover behavior around T ¼ T Ch , and while the decrease is much more pronounced than in the case of hssi, it differs qualitatively from the sharp decrease to zero exhibited by the chiral condensate [and thus also by the anomaly-related quantityχðTÞ defined by the LS relation (12)].
The isosymmetric pion decay constant f π ðTÞ ≡ f ll ðTÞ is the lower dashed curve in Fig. 2, starting at T ¼ 0 from our model-calculated value f π ¼ 92 MeV. It decreases rather quickly, in contrast to f ss ðTÞ [starting at f ss ðT ¼ 0Þ ¼ 119 MeV], the decay constant of the unphysical RLAss pseudoscalar. It exhibits a much "slower" T dependence, in accordance with the s-quark condensate hssiðTÞ.
The smooth, monotonic decrease of AðTÞ after T ∼ 0.7T Ch reflects the degree of gradual, crossover restoration of the U A ð1Þ symmetry with T. How this is reflected in the masses in the η-η 0 complex also depends on the ratios of AðTÞ with f 2 π ðTÞ, f π f ss ðTÞ, and f 2 ss ðTÞ in Eqs. (16) and (17). M 2 NS S ∝ AðTÞ=½f π ðTÞf ss ðTÞ decreases comparably to AðTÞ 1=2 , and 2AðTÞ=f ss ðTÞ 2 decreases even faster. Thus, M S ðTÞ [Eq. (17)] monotonically becomes the anomaly-free M ss ðTÞ in basically the same way as in Ref. [9], except now this process is not completed at T ¼ T Ch but rather [due to the AðTÞ crossover] drawn out until T ≈ 1.15T Ch .
These two limited increases of AðTÞ=f 2 π ðTÞ may be model dependent and are not important, but what is systematic and thus important is that the "light" decay constant f π ðTÞ makes ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi AðTÞ=f 2 π ðTÞ p more resilient to T than not only AðTÞ 1=4 itself, but also other anomalous mass contributions in Eqs. (16) and (17). Indeed, β Sho ðTÞ ¼ 2AðTÞ=f 2 π ðTÞ decreases only after T ≈ 0.95T Ch (contributing over a half of the η 0 mass decrease), and then again increases somewhat after T ≈ 1.15T Ch , to start definitively decreasing only after T ≈ 1.25T Ch , but even then slower than other anomalous contributions. This makes M NS ðTÞ larger enough than M S ðTÞ to increase ϕðTÞ to around 80°, and keep it there up to T ∼ 1.5T Ch (see Fig. 5).
This explains how the masses of the physical mesons η 0 and η (thick and thin solid curves in Figs. 3 and 4), exhibit the decrease of the mass of the heavier partner η 0 which is almost as strong as in the case [9] of the abrupt disappearance of the anomaly contribution, while on the contrary the lighter partner η now shows no sign of a decrease in mass around T ¼ T Ch , let alone an abrupt degeneracy with the pion. The latter happens in the case with the sharp phase transition because the fast disappearance of the whole M U A ð1Þ around T Ch can be accommodated only by a sharp change of the state mixing (ϕ → 0) to fulfill the asymptotic NS-S scenario immediately after T Ch . (See in particular Fig. 2 in Ref. [9]. Note that in our approach M η 0 ðTÞ cannot decrease by much more than a third of M U A ð1Þ , since the RLA M ss ðTÞ is the lower limit of M η 0 ðTÞ both in Ref. [9] and here.) In the present crossover case, however, T ¼ T Ch does not mark a drastic change in the mixing of the isoscalar states, but η 0 stays mostly η 0 and η stays mostly η 8 . Then, ΔM 2 Nevertheless, in M η 0 [Eq. (20)], the anomalous contributions from Eqs. (16) and (17) are all added together. The partial restoration of U A ð1Þ symmetry around T Ch , where around a third of the total U A ð1Þ-anomalous mass M U A ð1Þ goes away, is consumed almost entirely by the decrease of the η 0 mass over the crossover.
After T ≈ 1.15T Ch , M η 0 ðTÞ starts rising again, but this is expected since after T ≈ T Ch light pseudoscalar mesons start their thermal increase towards 2πT, which is twice the lowest Matsubara frequency of the free quark and antiquark. This rather steep joint increase brings all of the mass curves M P ðTÞ quite close after T ∼ 1.5T Ch . The kaon mass M K ðTÞ is not shown in Figs. 3 and 4 to maintain clarity by avoiding crowded curves, but at this temperature of the characteristic η-η 0 anticrossing, M K ðTÞ is roughly in between M π ðTÞ and the η mass, and is soon crossed by M η ðTÞ which tends to become degenerate with M π ðTÞ (as detailed below).
The rest of M U A ð1Þ ðTÞ [melting as 2 ffiffiffiffiffiffiffiffiffiffi AðTÞ p =f π ðTÞ] under 1.5T Ch is sufficiently large to keep M NS ðTÞ > M S ðTÞ and ϕ ≈ 80°. So a large ϕ makes θ positive, but not very far from zero, so that there we still have η 0 ≈ η 0 and η ≈ η 8 . This is also a fairly good approximation for T > 1.25T Ch , but there an even better approximation is η 0 ≈ η NS , M η 0 ðTÞ ≈ M NS ðTÞ and η ≈ η S , M η ðTÞ ≈ M S ðTÞ.

V. SUMMARY, DISCUSSION, AND CONCLUSIONS
We have studied the temperature dependence of the masses in the η 0 − η complex in the regime of the crossover restoration of chiral and U A ð1Þ symmetry. We relied on the approach of Ref. [11], which demonstrated the soundness of the approximate way in which the U A ð1Þ-anomaly effects on pseudoscalar masses were introduced and combined [24,27,[34][35][36][37] with chirally well-behaved DS RLA calculations in order to study η 0 and η. For T ¼ 0, this was demonstrated [11] model independently, with the only inputs being the experimental values of pion and kaon masses and decay constants, and the lattice value of the YM topological susceptibility. However, at T > 0 dynamical models are still needed to generate the temperature dependence of nonanomalous quantities through DS RLA calculations, and in this paper we used the same chirally correct and phenomenologically well-tested model as in numerous earlier T ≥ 0 studies (see, e.g., Refs. [9,24,31] and references therein).
Following Ref. [11], we assumed that the anomalous contribution to the masses is related to the full-QCD topological charge parameter (18), which contains the massive quark condensates. They give us the chiral crossover behavior for high T. This is crucial, since lattice QCD calculations have established that for the physical quark masses, the restoration of the chiral symmetry occurs as a crossover (see, e.g., Refs. [29,48,49] and references therein) characterized by the pseudocritical transition temperature T Ch .
Nevertheless, what happens with the U A ð1Þ restoration is still not clear [48,[50][51][52]. Whereas, e.g., Ref. [29] found its breaking as high as T ∼ 1.5T Ch , Ref. [53] found that above the critical temperature U A ð1Þ is restored in the chiral limit, and the JLQCD Collaboration [52] discussed the possible disappearance of the U A ð1Þ anomaly and pointed out the tight connection with the chiral symmetry restoration. Hence, there is a need to clarify "if, how (much), and when" [48] U A ð1Þ symmetry is restored. In such a situation, we believe instructive insight can be found in our study of how an anomaly-generated mass influences the η-η 0 complex, although this study is not done at the microscopic level.
Since the JLQCD Collaboration [52] has recently stressed that the chiral symmetry breaking and U A ð1Þ anomaly are tied for quark bilinear operators (as, e.g., in our Eqs. (12), (13), (18) and (19), where the chiral symmetry breaking drives the U A ð1Þ one through qq condensates), we again recall how Ref. [11] provided support for the earlier proposal of Ref. [9] relating DChSB to the U A ð1Þ-anomalous mass contributions in the η 0 -η complex. This adds to the motivation to determine the full-QCD topological charge parameter (18) on the lattice from simulations in full QCD with massive, dynamical quarks [besides the original motivation [13,14] to remove the systematic Oð1=N c Þ uncertainty of Eq. (15)]. More importantly, this connects the U A ð1Þ symmetry breaking and restoration to those of chiral symmetry. It connects them in basically the same way in both Refs. [9,11] (and here), except that the full-QCD topological charge parameter (18) enables the crossover U A ð1Þ restoration by allowing the use of the massive quark condensates. But, if the chiral condensate (i.e., of massless quarks) is used to extend the approach of Ref. [11] to finite temperatures, the T > 0 results are, in essence, very similar to those of Ref. [9]: the quick chiral phase transition leads to quick U A ð1Þ symmetry restoration at T Ch (consistent with Ref. [53]), which causes not only the empirically supported [5] decrease of the η 0 mass but also an even larger η mass decrease; if M 2 U A ð1Þ ðTÞ ∝ βðTÞ → 0 abruptly when T → T Ch , Eq. (10) mandates that M η ðT → T Ch Þ → M π ðT Ch Þ equally abruptly (as in Ref. [9]). However, no experimental indication for this has ever been seen, although this is a more drastic decrease than for the η 0 meson.
The present paper predicts a more realistic behavior of M η ðTÞ thanks to the smooth chiral restoration, which in turn yields the smooth, partial U A ð1Þ symmetry restoration (as far as the masses are concerned) making various actors in the η-η 0 complex behave quite differently from the abrupt phase transition (such as that in Ref. [9]). In particular, the η mass is now not predicted to decrease, but to only increase after T ≈ T Ch , just like the masses of other (almost-) Goldstone pseudoscalars, which are free of the U A ð1Þ anomaly influence. Similarly to T ¼ 0, η agrees rather well with the SUð3Þ flavor state η 8 until the anticrossing temperature, which marks the beginning of the asymptotic NS-S regime, where the anomalous mass contributions become increasingly negligible and η → η NS .
In contrast to η, the η 0 mass M η 0 ðTÞ does decrease similarly to the case of the sharp phase transition, where its lower limit [namely, M ss ðTÞ] is reached at T Ch [9]. Now, M η 0 ðTÞ at its minimum (which is only around 1.13T Ch because of the rather extended crossover) is some 20 to 30 MeV above M ss ðTÞ, after which they both start to grow appreciably, and M η 0 ðTÞ is reasonably approximated by M η 0 ðTÞ up to the anticrossing. The effective restoration of U A ð1Þ regarding the η-η 0 masses only occurs beyond the anticrossing at T ≈ 1.5T Ch , in the sense of reaching the asymptotic regime M η 0 ðTÞ → M ss ðTÞ. Another, less qualitatively illustrative but more quantitative criterion for the degree of U A ð1Þ restoration is that there, at T ≈ 1.5T Ch , M U A ð1Þ is still slightly above 40%, and at T ≈ 1.8T Ch still around 14% of its T ¼ 0 value. Thus, the decrease to the minimum of M η 0 ðTÞ around 1.13T Ch in any case signals only a partial U A ð1Þ restoration.
This M η 0 ðTÞ decrease is around 250 MeV, which is consistent with the current empirical evidence claiming that it is at least 200 MeV [5]. For comparison with some other approaches that explore the interplay of the chiral phase transition and axial anomaly, note that the η 0 mass decrease around 150 MeV is found in the functional renormalization group approach [54]. A very η 0 AND η MESONS AT HIGH T WHEN THE … PHYS. REV. D 99, 014007 (2019) recent analysis within the framework of the Uð3Þ chiral perturbation theory found that the (small) increase of the masses of π, K, and η after around T ∼ 120 MeV, is accompanied by the decrease of the η 0 mass, but only by some 15 MeV [55]. Admittedly, the crossover transition leaves more space for model dependence, since some model changes that would make the crossover even smoother would reduce our η 0 mass decrease. Nevertheless, there are also changes that would make it steeper, and those may, for example, help M η 0 ðTÞ saturate the M ss ðTÞ limit. Exploring such model dependences, as well as attempts to further reduce them at T > 0 by including more lattice QCD results, must be relegated to future work. However, here we can already note a motivation for varying the presently isosymmetric model current u-and d-quark mass of 5.49 MeV. Since it is essentially a phenomenological model parameter, it cannot be quite unambiguously and precisely related to the somewhat lower Particle Data Group values m u ¼ 2.2 þ0.5 −0.4 MeV and m d ¼ 4.70 þ0.5 −0.3 MeV [56]. Still, their ratio m u =m d ¼ 0.48 þ0.07 −0.08 is quite instructive in the present context, since the QCD topological susceptibility χ [Eq. (19)] and charge parameter A [Eq. (18)] contain the current-quark masses in the form of harmonic averages of m q hqqi (q ¼ u, d, s). Since a harmonic average is dominated by its smallest argument, our χ and A are dominated by the lightest flavor, providing the motivation to venture beyond the precision of the isospin limit and in future work explore the maximal isospin violation scenario [57] within the present treatment of the η-η 0 complex.