Temperature Dependence of the Axion Mass in a Scenario Where the Restoration of Chiral Symmetry Drives the Restoration of the $U_A(1)$ Symmetry

The temperature ($T$) dependence of the axion mass is predicted for $T$'s up to $\sim 2.3 \times$ the chiral restoration temperature of QCD. The axion is related to the $U_A(1)$ anomaly. The squared axion mass $m_a(T)^2$ is, modulo the presently undetermined scale of spontaneous breaking of Peccei-Quinn symmetry $f_a$ (squared), equal to QCD topological susceptibility $\chi(T)$ for all $T$. We obtain $\chi(T)$ by using quark condensates calculated in two effective Dyson-Schwinger models of nonperturbative QCD. They exhibit the correct chiral behavior, including the dynamical breaking of chiral symmetry and its restoration at high $T$. This is reflected in the $U_A(1)$ symmetry breaking and restoration through $\chi(T)$. In our previous studies, such $\chi(T)$ yields the $T$-dependence of the $U_A(1)$-anomaly-influenced masses of $\eta'$ and $\eta$ mesons consistent with experiment. This in turn supports our prediction for the $T$-dependence of the axion mass. Another support is a rather good agreement with the pertinent lattice results. This agreement is not spoiled by our varying $u$ and $d$ quark mass parameters out of the isospin limit.


Introduction
The axion, one of the oldest hypothetical particles beyond the Standard Model, intensely sought for by many experimentalists already for 40 years now, still escapes detection [1]. It was introduced theoretically [2][3][4][5] to solve the so-called Strong CP problem of QCD. The problem is that no experimental evidence of CP-symmetry violation has been found in strong interactions, although the QCD Lagrangian L QCD (x) can include the so-called θ-term L θ (x) = θ Q(x) where gluon field strengths F b µν (x) form the CP-violating combination Q(x) named the topological charge density: Whereas Q(x) can be re-cast in the form of a total divergence ∂ µ K µ , discarding L θ is not justified even if F b µν (x) vanish sufficiently fast as x → ∞. Specifically, FF = ∂ µ K µ can anyway contribute to the action integral, since in QCD there are topologically nontrivial field configurations such as instantons. They are important for, e.g., obtaining the anomalously large mass of the η ′ meson. Also, precisely the form (1) from the θ-term appears in the axial anomaly, breaking the U A (1) symmetry of QCD -see Equation (2).
For these reasons, one needs L θ (x) = θ Q(x) in the QCD Lagrangian, as reviewed briefly in Section 1 of Ref. [6]. Moreover, the Strong CP problem cannot be removed by requiring that the coefficient θ = 0, since QCD is an integral part of the Standard Model, where weak interactions break the CP symmetry. This CP violation comes from the complex Yukawa couplings, yielding the complex CKM matrix [7,8] and the quark-mass matrix M which is complex in general. To go to the arXiv:1909.09879v2 [hep-ph] 28 Oct 2019 The breaking of U A (1) by the anomaly makes the flavor singlet (a = 0) axial current of quarks, A µ 0 (x) = ∑ q=u,d,sq (x) γ µ γ 5 q(x), not conserved even in the chiral limit: unlike the corresponding octet currents A µ a (x), a = 1, 2, ..., N 2 f − 1. In the chiral limit, the current masses of non-heavy quarks all vanish, m q → 0 (q = u, d, s), but the divergence of the singlet current A µ 0 is not vanishing due to the U A (1) anomaly contributing no other but the topological charge density operator Q(x) (1)-that is, precisely the quantity responsible for the strong CP problem.
The quantity related to the U A (1)-anomalous mass in the η ′ -η complex is the QCD topological susceptibility χ, where T denotes the time-ordered product. Figure 1 shows how anomaly contributes to the mass matrix (in the basis of quark-antiquark (qq) pseudoscalar bound states P) by depicting how hidden-flavor qq pseudoscalars mix, transiting through anomaly-dominated gluonic intermediate states. anomaly-induced, hidden-flavor-mixing transitions from pseudoscalar quark-antiquark states P = qq to P ′ = q ′q′ include both possibilities q = q ′ and q ≠ q ′ . Springs symbolize gluons. All lines and vertices are dressed in accord with the nonperturbative QCD. Nonperturbative configurations are essential for nonvanishing anomalous mass [20] contribution to η 0 ∼ η ′ , since Q(x) is a total divergence. The gray blob symbolizes the infinity of all intermediate gluon states enabling such transitions, so that the three bold dots represent any even [21] number of additional gluons. Just one of infinitely many, but certainly the simplest realization thereof, is when such a transition is mediated by just two gluons (and no additional intermediate states), whereby the above figure reduces to the so-called "diamond graph".

On Some Possibilities of Modeling the U A (1) Anomaly Influence
Light pseudoscalar mesons can be studied by various methods. We have preferred using [17,[21][22][23][24][25][26][27][28][29][30][31][32] the relativistic bound-state approach to modeling nonperturbative QCD through Dyson-Schwinger equations (DSE), where, if approximations are consistently formulated, model DSE calculations also reproduce the correct chiral behavior of QCD. This is of paramount importance for descriptions of the light pseudoscalar mesons, which are quark-antiquark bound states but simultaneously also the (almost-)Goldstone bosons of DChSB of QCD. (For general reviews of the DSE approach, see, e.g., Refs. [33][34][35][36]. About our model choice at T = 0 and T > 0, in the further text see especially the Appendix.) Figure 1 illustrates how hard computing would be "in full glory" the U A (1)-anomalous mass and related quantities, such as the presently all-important topological susceptibility (3) in the DSE approach with realistically modeled QCD interactions -especially if the calculation should be performed in a consistent approximation with the calculation of the light pseudoscalar bound states, to preserve their correct chiral behavior. (For this reason, they have most often been studied in the rainbow-ladder approximation of DSE, which is inadequate for the anomalous contributions [33][34][35][36], as also Figure 1 shows.) However, our DSE studies of pseudoscalar mesons have been able to address not only pions and kaons, but also η ′ and η mesons, for which it is essential to include the anomalous U A (1) symmetry breaking at least at the level of the masses. This was done as described in Refs. [21,24,29,30,37], namely exploiting the fact that the U A (1) anomaly is suppressed in the limit of large number of QCD colors N c [38,39]. This allows treating the anomaly contribution formally as a perturbation with respect to the non-anomalous contributions to the η and η ′ masses [21,24,29]. This way we avoid the need to compute the anomalous mass contribution together, and consistently, with the non-anomalous, chiral-limit-vanishing parts of the masses. The latter must be evaluated by some appropriate, chirally correct method, and our preferred tool-the relativistic bound-state DSE approach [33][34][35][36]-is just one such possibility. The point is that they comprise the non-anomalous part of the η ′ -η mass matrix, to which one can add, as a first-order perturbation, the U A (1)-anomalous mass contribution M U A (1) -and it does not have to be modeled, but taken from lattice QCD [29]. Specifically, at T = 0, M U A (1) can be obtained from χ YM , the topological susceptibility of the (pure-gauge) Yang-Mills theory, for which reliable lattice results have already existed for a long time 3 [40][41][42].
This can be seen from the remarkable Witten-Veneziano relation (WVR) [38,39] which in a very good approximation relates the full-QCD quantities (η ′ , η and K-meson masses M η ′ , M η and M K respectively, and the pion decay constant f π ), to the pure-gauge quantity χ YM : The right-hand-side must be the total U A (1)-anomalous mass contribution in the η ′ -η complex, since in the combination on the left-hand-side everything else cancels at least to the second order, O(m 2 q ), in the current quark masses of the three light flavors q = u, d, s. This is because the non-anomalous, chiral-limit-vanishing parts M qq ′ of the masses of pseudoscalar mesons 4 P ∼ qq ′ composed of sufficiently light quarks, satisfy the Gell-Mann-Oakes-Renner (GMOR) relation with their decay constants f qq ′ and the quark-antiquark (qq) condensate signaling DChSB: Here f ch.lim qq ′ = f qq ′ (m q , m q ′ → 0), and ⟨q q⟩ 0 denotes the massless-quark condensate, i.e., the qq chiral-limit condensate, or "massless" condensate for short. (In the absence of electroweak interactions, the "massless" condensates have equal values for all flavors: ⟨q q⟩ 0 = ⟨q ′ q ′ ⟩ 0 .) It turns out that even s-flavor is sufficiently light for Equation (5) to provide reasonable approximations.
Using WVR and χ YM to get the anomalous part of the η ′ and η masses is successful [29] only for T ∼ 0, or at any rate, T's well below T c , the pseudocritical temperature of the chiral transition. In the absence of a systematic re-derivation of WVR (4) at T > 0, its straightforward extension (simply replacing all quantities by their T-dependent versions) is tempting, but was found [31] unreliable and with predictions in a drastic conflict with experiment when T starts approaching T c . This is because 3 This is in contrast with χ = χ QCD , the full-QCD topological susceptibility (3), which is much harder to find on the lattice because of the light quark flavors. χ = χ QCD approaches χ YM only if one takes the quenched limit of infinitely massive quarks, χ YM = χ quench , since quarks then disappear from the loops of Equation (3). 4 The combinations P ∼ qq ′ need not always pertain to physical mesons. The pseudoscalar hidden-flavor states uū, dd, ss are not physical as long as the U A (1) symmetry is not restored (i.e., the anomaly effectively turned off, see around Equation (2.6) in Ref. [43] for example), but build the SU(3) states η 0 , η 8 and π 0 . the full-QCD quantities M η ′ (T), M η (T), M K (T) and f π (T) have very different T-dependences from the remaining quantity χ YM (T), which is pure-gauge and thus much more resilient to increasing temperature: the critical temperature of the pure-gauge, Yang-Mills theory, T YM , is more than 100 MeV higher than QCD's T c = (154 ± 9) MeV [44,45]. The early lattice result T YM ≈ 260 MeV [40,41] is still accepted today [46], and lattice groups finding a different T YM claim only it is even higher, for example T YM = (300 ± 3) MeV of Gattringer et al. [42]. (There are even some claims about experimentally established T YM = 270 MeV [47]. ) We thus proposed in 2011. [48] that the above mismatch of the T-dependences in WVR (4) can be removed if one invokes another relation between χ YM and full-QCD quantities, to eliminate χ YM , i.e., substitute pertinent full-QCD quantities instead of χ YM at T > 0. This is the Leutwyler-Smilga (LS) relation, Equation (11.16) of Ref. [49], which we used [17,37,48] in the inverted form (and in our notation): to express (at T = 0) pure-gauge χ YM in terms of the full-QCD topological susceptibility χ ≡ χ QCD , the current quark masses m q , and ⟨qq⟩ 0 , the condensate of massless, chiral-limit quarks. The combination which these full-QCD quantities comprise, i.e., the right-hand-side of the LS relation (6), we denote (for all T) by the new symbolχ for later convenience -that is, for usage at high T, where the equality (6) with χ YM does not hold. The remarkable LS relation (6) holds for all values of the current quark masses. In the limit of very heavy quarks, it correctly yields χ → χ quenched = χ YM for m q → ∞, it but it also holds for the light m q . In the light-quark sector, the QCD topological susceptibility χ can be expressed as [49][50][51]: where C m represents corrections of higher orders in light-quark masses m q . Thus, it is small and often neglected, leaving just the leading term as the widely used [52] expression for χ in the light-quark, N f = 3 sector. However, setting C m = 0 in the light-quark χ (7) returns us χ YM = ∞ through Equation (6) [48]. Or conversely, setting χ YM = ∞ in the LS relation (6), gives the leading term of χ (7) (see also Ref. [43]). This can be a reasonable, useful limit considering that in reality χ YM χ ≳ 40. Nevertheless, in our previous works on the η ′ -η complex at T > 0 [17,48] we had to fit C m (and parameterize it with Ansätze at T > 0), since we needed the realistic value of χ YM (T = 0) =χ(T = 0) from lattice to reproduce the well-known masses of η and η ′ at T = 0. (However, just for χ(T) this is not necessary.) Replacing [48] χ YM (T) by the full-QCD quantityχ(T) obviously keeps WVR at T = 0, but avoids the 'YM vs. QCD' T-dependence mismatch with f π (T) and the LHS of Equation (4), so it is much more plausible to assume the straightforward extension of T-dependences. The T-dependences of χ(T) (7) andχ(T), and thus also of the anomalous parts of the η and η ′ masses, are then obviously dictated by ⟨qq⟩ 0 (T), the "massless" condensate.
General renormalization group arguments suggest [53] that QCD with three degenerate light-quark flavors has a first-order phase transition in the chiral limit, whereas in QCD with (2+1) flavors (where s-quark is kept significantly more massive) a second-order chiral-limit transition 5 is also possible and even more likely [62][63][64]. What is important here, is that in any case the chiral-limit condensate ⟨qq⟩ 0 (T) drops sharply to zero at T = T c . (The dotted curve in Figure 2 is just a special example thereof, namely ⟨qq⟩ 0 (T) calculated in Ref. [48] using the same model as in [17] and here.) This causes a similarly sharp drop ofχ(T) and χ(T). (We may be permitted to preview similar dotted 5 This is a feature exhibited by DSE models, or at least by most of them, through the characteristic drop of their chiral-limit, massless qq condensates [54][55][56][57][58][59][60][61]. curves in Figures 5 and 7 in the next section and anticipate that in this case one would get a massless axion at T = T c .) This was also the reason, besides the expected [65,66] drop of the η ′ mass, Ref. [48] also predicted so drastic drop of the η mass at T = T c , that it would become degenerate with the pion. However, no experimental indication whatsoever for a decreasing behavior of the η mass, and much less for such a conspicuous sharp mass drop, has been noticed to this day, which seems to favor theoretical descriptions with a smooth crossover. Also, recent lattice QCD results (see [44,67,68] and their references) show that the chiral symmetry restoration is a crossover around the pseudocritical transition temperature T c .

All values in [MeV]
A 1/4 The relative-temperature T Tc dependences (where T Ch ≡ Tc) of (the 3 rd root of the absolute value of) the qq condensates, and of (the 4 th root of) the topological susceptibility χ(T) and the full QCD topological charge parameter A(T). Everything was calculated in the isosymmetric limit using the separable rank-2 DSE model which we had already used in Refs. [17,31,48]. (See the Appendix for the model interaction form and parameters, and the first line of Table 1 for the numerical values of the condensates and χ at T = 0.) Only the chiral-limit condensate ⟨qq⟩ 0 (T) falls steeply to zero at T = T Ch , indicative of the second-order phase transition. This sharp transition would through (7) and (6) be transmitted to, respectively, the chiral-limit χ(T) andχ(T), and ultimately to the η and η ′ masses in Ref. [48]. The highest curve (dash-dotted) and the second one from above (dashed) are (3 rd roots of the absolute values of) the condensates ⟨ss⟩(T) and ⟨ūu⟩(T), respectively. Their smooth crossover behaviors carry over to χ(T) and A(T) through Equations (9) and (8), respectively, leading to the empirically acceptable predictions [17] for the T-dependence of the η and η ′ masses. It turns out that χ(T) (9) also gives the smooth crossover behavior also to the T-dependence of the axion mass.
To describe such a crossover behavior of the chiral transition, we incorporated [17,37] into our approach Shore's generalization [69] of the Witten-Veneziano relation. To be precise, we studied it at T = 0 already in 2008 [30] and adapted it to our DSE bound-state context by applying some very plausible simplifications [20]. The more recent reference [37] presented the analytic, closed-form solutions to Shore's equations for the pseudoscalar meson masses. These solutions showed that Shore's approach is then actually quite similar to the original WVR, leading to a similar η-η ′ mass matrix [37].
Presently, the most important advantage is that Shore's generalization leads to the crossover T-dependence. As obtained in Section 3 for two specific interactions modeling the nonperturbative QCD interaction, these condensates exhibit a smooth, crossover chiral symmetry transition around T c . Here Figure 2 illustrates this generic behavior by displaying the results obtained in Section 3 for a specific DSE model: the higher current quark mass, the smoother the crossover behavior, which then results in the crossover behavior also of other quantities, like the presently all-important quantity, the QCD topological susceptibility χ(T).
This comes about as follows: the quantity which in Shore's mass relations [69] has the role of χ YM in the Witten-Veneziano relation, is called the full-QCD topological charge parameter A. (Shore basically took over this quantity from Di Vecchia and Veneziano [50].) At T = 0, it is approximately equal to χ YM in the sense of 1 N c expansion. Shore uses A to express the QCD susceptibility χ through a relation similar to the Leutwyler-Smilga relation (see Equation (2.11) and (2.12) in Ref. [69]), but using the condensates ⟨ūu⟩, ⟨dd⟩, ⟨ss⟩ of realistically massive u, d, s quarks. The inverse relation, yielding A (with the opposite sign convention), is the most illustrative for us: Obviously, it is analogous to the inverted LS relation (6) definingχ , except that A is expressed through "massive" condensates. (If they are all replaced by ⟨qq⟩ 0 , then A →χ.) They are in principle different for each flavor, but in the limit of usually excellent isospin symmetry, ⟨ūu⟩ = ⟨dd⟩.
One can examine the limiting assumption A = ∞ in analogy with taking the limit χ YM = ∞ compared to χ = χ QCD . Then, for A = ∞ (be it in our Equation (8) or Shore's Equations (2.11) and (2.12) for χ), one recovers the leading term of the QCD topological susceptibility expressed by the "massive" condensates. However, if one needs a finite A, as in η-η ′ calculations [17] where one needs to reproduce A ≈ χ YM , one also needs the appropriate correction term C ′ m , just as C m in Equation (7), so that: Again, C ′ m is a very small correction term of higher orders in the small current quark masses m q (q = u, d, s), and we can neglect it in the present context, where we actually have a simpler task than finding T-dependence of the η and η ′ masses in Ref. [17]. Since it turns out that for determining the T-dependence of the mass of the QCD axion we do not need to find A, we set C ′ m = 0 in this paper throughout. One needs just the topological susceptibility χ(T) for that, and just the leading term of (9) will suffice to yield the crossover behavior found on lattice (e.g., in Refs. [70][71][72]).

The Axion Mass from the Non-Abelian Axial Anomaly of QCD
Peccei and Quinn (PQ) introduced [2,3] a new global symmetry U(1) PQ which is broken spontaneously at some very large, but otherwise still unknown scale f a > 10 8 GeV [1,73], which determines the absolute value of the axion mass m a . Nevertheless, this constant cancels from ratios such as m a (T) m a (0), where T is temperature. Thus, useful insights and applications, such as those involving the nontrivial part of axion T-dependence, are possible in spite of f a being presently unknown.
The factor in the axion mass which carries the nontrivial T-dependence, is the QCD topological susceptibility χ(T). This quantity is also essential for our description of the η ′ -η complex at T > 0, since it relates the T-dependence of the anomalous breaking of U A (1) symmetry.

The Axion as the Almost-Goldstone Boson of the Peccei-Quinn Symmetry
The pseudoscalar axion field a(x) arises as the (would-be massless) Goldstone boson of the spontaneous breaking of the PQ symmetry U(1) PQ [4,5]. The axion contributes to the total Lagrangian its kinetic term and its interaction with fermions of the Standard model. Nevertheless, what is important for the resolution of the strong CP problem, is that the axion also couples to the topological charge density operator Q(x) defined in Equation (1) and generating the U A (1)-anomalous term in Equation (2). Theθ-term in L QCD thus changes into Because of this axion-gluon coupling, the U(1) PQ symmetry is also broken explicitly by the U A (1) anomaly (gluon axial anomaly). This gives the axion a nonvanishing mass, m a ≠ 0 [4,5].
Gluons generate an effective axion potential, and its minimization leads to the axion expectation value ⟨a⟩ which makes the modified coefficient of Q(x) in Equation (10) Obviously, the experiments excluding the strong CP violation, such as [9], have in fact been finding that consistent with zero isθ ′ , the coefficient of Q(x) in the QCD Lagrangian when the axion exists. The strong CP problem is thereby solved, irrespective of the initial value ofθ. (Relaxation from anyθ-value in the early Universe towards the minimum atθ = − ⟨a⟩ f a is called misalignment production. The resulting axion oscillation energy is a good candidate for cold dark matter [10][11][12][13][14][15].)

Axion Mass from the Topological Susceptibility from Condensates of Massive Quarks
Modulo the (squared) Peccei-Quinn scale f 2 a , the axion mass squared is at all temperatures T given by the QCD topological susceptibility [11,43,[70][71][72]74] very accurately [75,76] (up to negligible corrections of the order (pion mass) 2 f 2 a ): as revealed by the quadratic term of the expansion of the effective axion potential.
We explained in Section 2 how the U A (1) symmetry-breaking quantity χ(T) can be obtained through Equation (9) as a prediction of any method which can provide the quark condensates ⟨qq⟩(T) (q = u, d, s). Thus, one can get the T-dependence of the axion mass (11) through the mechanism where DChSB drives the U A (1) symmetry breaking. (And conversely, of course: the chiral restoration then drives the restoration of U A (1) symmetry.) An excellent tool to study DChSB, and in fact "produce" it in the theoretical sense, is one of the basic equations of the DSE approach -the gap equation. The most interesting thing it does for nonperturbative QCD is explaining the notion of the constituent quark mass around 1 3 of the nucleon mass M N by generating them via DChSB, in the same process which produces the qq condensates. Thanks to this, qq condensates can be evaluated from dressed quark propagators. Specifically, hadronic-scale large (∼ M N 3 at small momenta p) dressed quark-mass functions M q (p 2 ) ≡ B q (p 2 ) A q (p 2 ) are generated despite two orders of magnitude lighter current quark masses m q , and in fact even in the chiral limit, when m q = 0! This happens in low-energy QCD thanks to nonperturbative dressing via strong dynamics, making strongly dressed quark propagators S q (p) out of the free quark propagators S free q : The solution for the dressed quark propagator S q (p) of the flavor q, i.e., the dressing functions A q (p 2 ) and B q (p 2 ), are found by solving the gap equation where Σ q (p) is the corresponding DChSB-generated self-energy, for example, Equation (14) if the rainbow-ladder truncation is adopted.
In the present work, all we want to model of nonperturbative QCD are the condensates ⟨qq⟩ at all temperatures T, and for that the solutions of the quark-propagator gap Equation (13) are sufficient, i.e., we do not need the Bethe-Salpeter equation (BSE) for the qq ′ pseudoscalar bound states. However, we want the same condensates, and basically the same (leading term of) the topological susceptibility as we had in our related η-η ′ paper [17], and in a number of earlier papers, such as Refs. [31,48]. Thus, we now use the same model interaction we have been using then in the consistent rainbow-ladder (RL) truncation of DSE's to produce chirally correctly behaving pseudoscalar mesons -that is, with the non-anomalous parts of their masses given by GMOR (5).
Thus, the quark self-energy in the gap Equation (13) in the RL truncation is where D ab µν (k) eff is an effective gluon propagator, which should be chosen to model the nonperturbative, low-energy domain of QCD. This can be done in varying degrees of DSE modeling, depending on the variety of problems one wants to treat [33][34][35][36]77]. For example, in the context of low-energy meson phenomenology, if one does not aim to address problems of perturbative QCD, it is better not to include the perturbative part of the QCD interaction. Otherwise, in the words of very authoritative DSE practitioners, "the logarithmic tail and its associated renormalization represent an unnecessary obfuscation." [78] In medium, the original O(4) symmetry is broken to O(3) symmetry. The most general form of the dressed quark propagator then has four independent tensor structures and four corresponding dressing functions. At nonvanishing temperature, T > 0, we use the Matsubara formalism, where four-momenta decompose into three-momenta and Matsubara frequencies: p = (p 0 , ⃗ p) → p n = (ω n , ⃗ p). Therefore, the (inverted) dressed quark propagator S q (p) (13) becomes (15) (The T-dependence of the propagator dressing functions is understood and, to save space, is not indicated explicitly, except in the Appendix.) Nevertheless, the last dressing function D q ( ⃗ p 2 , ω n ) is so very small that it is quite safe and customary to neglect it -e.g., see Refs. [34,79]. Thus, also we set D q ≡ 0, leaving only A q , C q and B q .
For applications in involved contexts, such as calculations at T > 0, appropriate simplifications are very welcome for tractability. This is why in Refs. [17,31,48] and presently, we adopted relatively simple, but phenomenologically successful [30][31][32]77,80] separable approximation [77]. The details on the functional form and parameters of the presently used model interaction can be found in the Appendix.
As already pointed out in the original Ref. [77], the model Ansätze for the nonperturbative low-energy interaction ("interaction form factors") are such that they provide sufficient ultraviolet suppression. Therefore, as noted already in Ref. [77], no renormalization is needed and the multiplicative renormalization constants, which would otherwise be needed in the gap Equation (13) with (14), are 1. The usual expression for the condensate of the flavor q then becomes where Tr is the trace in Dirac space, and the combined integral-sum symbol indicates that when the calculation is at T > 0, the four-momentum integration decomposes into the three-momentum integration and summation over fermionic Matsubara frequencies ω n = (2n + 1)πT, n ∈ Z.  (17) and introduced by Ref. [81]. The lattice data points are from Figure  6 of Ref. [82], but scaled for the critical temperatures Tχ from their Table 2, which is different for the "crosses" (data points [82] for mπ ≈ 370 MeV) and "bars" (data points [82] for mπ ≈ 210 MeV). The lower, green curve results from the R ⟨ψψ⟩ (17) subtraction of our u-quark condensate. The upper, red curve is the T-dependence of our u-quark condensate when regularized in the usual way (see text).
As is well known, the condensates (16) are finite only for massless quarks, m q = 0, i.e., only ⟨qq⟩ 0 is finite, while the "massive" condensates are badly divergent, and must be regularized, i.e., divergences must be subtracted. Since the subtraction procedure is not uniquely defined, the chiral condensate at nonvanishing quark mass is also not uniquely defined. However, the arbitrariness is in practice slight and should rather be classified as fuzziness. It should not be given too large importance in the light of small differences between the results of various sensible procedures.
Our regularization procedure is subtracting the divergence-causing m q (∼ several MeV) from the scalar quark-dressing function B q (p 2 ) (∼ several hundred MeV) whenever it is found in the numerator of the condensate integrand. To justify our particular regularization of massive condensates as physically meaningful and sensible, we have examined its consistency with two different subtractions used on lattice [81][82][83] and in a recent DSE-approach paper [84].
We shall now test our massive condensates obtained from the separable rank-2 DSE model (see the Appendix), whose regularized versions have already been shown in Figure 2.
Let us first consider the subtraction on lattice (normalized to 1 for T = 0) first proposed in Ref. [81] in their Equation (17), rewritten in our notation and applied to our condensate of u-quarks (and of course d-quarks in the isospin limit): In Figure 3, the upper, red curve shows (normalized) u-quark condensate ⟨ū u⟩(T) ⟨ū u⟩(0) when regularized in the usual way, by subtracting m q from B q (p 2 ) in the numerator of the condensate integrand. It agrees very well with the lattice regularization R ⟨ūu⟩ (17) of our condensate ⟨ū u⟩(T), represented by the green curve. The agreement with the lattice data points taken (if pertinent) from Table 6 of Ref. [82] is also rather good.
Next, we examine the consistency of our subtraction with the most usual condensate subtraction on the lattice, which combines the light and strange quark condensates and their masses like this: with the lattice data of Ref. [83]. The agreement is very good, which implies also the agreement with the subtracted condensates in the recent DSE paper [84], which made this successful comparison first (in its Figure 3).
To conclude: results shown in Figures 3 and 4 demonstrate that certain arbitrariness in the choice of regularization does not disqualify our massive condensates from useful applications, such as using them in Equation (9) to make predictions on the topological susceptibility.

Axion Mass and Topological Susceptibility-Results from the Rank-2 Separable Model in the Isosymmetric Limit
Our result for χ(T) 1 4 = m a (T) f a is presented in Figure 5 as a solid curve and compared, up to T ≈ 2.3 T c , with the corresponding results of two lattice groups [70,71], rescaled to the relative-temperature T T c . (Table 1 gives numerical values of our results at T = 0.) In our case, the results for χ(T) and condensates ⟨ūu⟩(T), ⟨dd⟩(T) and ⟨ss⟩(T) needed to obtain it, are predictions of the dynamical DSE model used in the T > 0 study of η ′ -η [17]. This is the same modeling of the low-energy, nonperturbative QCD interactions as we have already employed in our earlier studies of light pseudoscalar mesons at T ≥ 0 [30][31][32]48]: the separable model interactionsee, e.g., [77,85], and references therein. We have adopted the so-called rank-2 variant from Ref. [77]. The adopted model with our choice of parameters is defined in detail in the Appendix of the present work, after the subsection II.A of Ref. [31] and Ref. [86]. It employs the model current Contrary to, e.g., Ref. [48], where the condensate of massless quarks ⟨qq⟩ 0 (T) was used, in Ref. [17] and here we follow Shore [69]   . The relative-temperature T Tc dependence of the (normalized) subtracted quark condensate (19) from the lattice [83] (blue squares) and from our condensates. Slightly lower, green curve results from our unsubtracted condensates plugged in Equation (19), while the very slightly higher, red curve is from our already subtracted condensates.
masses. The smooth, crossover behavior around the pseudocritical temperature T c for the chiral transition (now confirmed at vanishing baryon density by lattice studies such as [44,67,68,70,71]), is obtained thanks to the DChSB condensates of realistically massive light quarks -i.e., the quarks with realistic explicit chiral symmetry breaking [17]. In contrast, using in Equation (9) the massless-quark condensate ⟨qq⟩ 0 (which drops sharply to zero at T c ) instead of the "massive" ones, would dictate a sharp transition of the second order at T c [17,48] also for χ(T), illustrated in Figure 5 by the dotted curve. This would of course imply that axions are massless [87] for T > T c . It is of academic interest to know what consequences would be thereof for cosmology, but now it is clear that only crossover is realistic [70,71].
The rather good agreement with lattice in Figure 5 resulted without any refitting of this model, either in Ref. [17] for η ′ and η, or in this subsection. The model is in the isosymmetric limit, m u = m d ≡ m l , which is perfectly adequate for most purposes in hadronic physics. Nevertheless, the QCD topological susceptibility χ in its version (9) contains the current quark masses in the form of harmonic averages of m q ⟨qq⟩ (q = u, d, s). A harmonic average is dominated by its smallest argument, and presently this is the lightest current-quark-mass parameter, motivating us to investigate the changes occurring beyond the isospin symmetric point.
This seems not to bode well for the attempts out of the isosymmetric limit, because lowering the values of the current masses seems to threaten yielding unacceptably low values of the topological susceptibility. Indeed, taking the central values from the current quark masses m u = 2. s . The rest of model parameters, namely those in the Ansatz functions F 0 (p 2 ) and F 1 (p 2 ) modeling the strength of the rank-2 nonperturbative interaction (see the Appendix), are also not varied.) Table 1. For the both variants of the DSE separable model (with the rank-2 and rank-1 interaction Ansatz) used in the present paper, various sets of values of the model quark-mass parameters mq (q = u, d, s) are related to the model results for the topological susceptibility χ and the "massive" condensates ⟨q q⟩ at T = 0. The topological susceptibility χ varies because of varying mq and (to a much lesser extent) because of the changes of "massive" condensates induced by changes of the quark-mass parameters mq. The massless-quark condensate, ⟨q q⟩ 0 , depends only on the dynamical DSE model: always ⟨q q⟩ 0 = −216.25 3 MeV 3 for the rank-2 model, and ⟨q q⟩ 0 = −248.47 3 MeV 3 for the rank-1 model. Thus, the topological susceptibility χ 0 calculated with the chiral-limit condensate, varies for a given model only because of varying values of mq. All values are in MeV (or the indicated 3rd or 4th powers of MeV). The T-dependence of the resulting χ(T) 1 4 is given by the short-dashed black curve in Figure 6. Except its better agreement with the lattice results [70,71] at low T, the new (dashed) χ(T) 1 4 curve is very close to the isosymmetric (solid) curve. Now we will check the model dependence by comparing our results presented so far (obtained in the rank-2 model) with those we get in the separable rank-1 model of Ref. [77]. It is similar to the previously considered rank-2 one by modeling the low-energy, nonperturbative QCD interaction with an Ansatz separating the momenta p a , p b of interacting constituents, but is of a simpler form, proportional to just F 0 (p 2 a ) F 0 (p 2 b ). Its presently interesting feature is that for similar quark-mass parameters, it yields significantly larger condensates than those in the separable rank-2 model, and thus also larger χ. (This also holds at low and vanishing T even for χ(T) calculated using only the "massless" condensate ⟨qq⟩ 0 (T). This case is depicted in Figure 7 as the dotted curve.) The original rank-1 model employs the light-quark current mass parameters in the isosymmetric limit: m u = m d ≡ m l = 6.6 MeV [77]. However, in Equation (9) 3 . This gives too large topological susceptibility at T = 0, namely χ(0) = (84.08 MeV) 4 . Nevertheless, for large T, it also falls with T somewhat faster than the rank-2 χ(T), since rank-1 condensates fall with T somewhat faster than the rank-2 ones.
The isosymmetric rank-1 χ(T) is depicted by the solid black curve in Figure 7, showing that it actually falls with T faster even than χ(T)'s from lattice [70,71] for practically all T's high enough to induce changes. Then, comparing Figures 6 and 7 shows that the lattice high-T results are in between high-T results of the two separable models.  Figure 5 to help resolve the short-dashed curve from the solid curve representing again the isosymmetric case of the same separable rank-2 model. Also, the lattice results [70,71] are again depicted as in Figure 5.
To go out of the isospin limit with the rank-1 model, we again require that the changed parameters m fit u and m fit d (together with the condensates resulting from them) fit the currently most precise T = 0 value of the topological susceptibility, χ = (75.44 MeV) 4 [76], while also obeying m fit u m fit d = 0.48, i.e., the central value of the PDG [1] m u m d ratio. (Again, other model parameters including m s are not varied: m s ≡ m fit s .) In our rank-1 model, these requirements yield m fit d = 6.56 MeV, i.e., it practically remained the same as in the originally fitted model [77]. Of course, its condensate ⟨dd⟩ also remains the same. The lightest flavor has the mass parameter lowered to m fit u = 3.15 MeV. Now it has only slightly lower condensate ⟨ūu⟩(T = 0) = (−248.87 MeV) 3 , which is even closer to the "massless" ⟨qq⟩ 0 (T = 0). Nevertheless, ⟨ūu⟩(T) retains the crossover behavior for T > 0, although it falls with T steeper than more "massive" condensates.
The resulting non-isosymmetric χ(T) 1 4 = m a (T) f a is in Figure 7 shown as the short-dashed black curve, which is everywhere consistently the lowest (among the "massive", crossover curves).

Summary and Discussion
In the DSE framework, we have obtained predictions for the nontrivial part of the T-dependence of the axion mass m a (T) = χ(T) f a , Equation (11) Figure 7. From the calculation in the separable rank-1 DSE model [77], the relative-temperature T Tc dependence of (the leading term of) χ(T) 1 4 is represented by: (i) the solid curve for the isosymmetric case with mu = m d = 6.6 MeV and ms = 142 MeV, (ii) the dotted curve is for the same mass parameters, but with all condensates approximated by the "massless" condensate ⟨qq⟩ 0 (T), and (iii) the short-dashed curve for the non-isosymmetric case m fit u = 3.15 MeV, m fit d = 6.56 MeV, while m fit s ≡ ms = 142 MeV. The pertinent lattice results are presented in the same way as in the previous two Figures: the dash-dotted and long-dashed curves extracted, respectively, from Refs. [70] and [71]. Colors online. empirically successful dynamical models of the separable type [77] to model nonperturbative QCD at T > 0. We also studied the effects of varying the mass parameters of the lightest flavors out of the isospin limit, and found that our χ(T), and consequently m a (T), are robust with respect to the non-isosymmetric refitting of m u and m d = m u 0.48 .
All these results of ours on χ(T), and consequently the related axion mass, are in satisfactory agreement with the pertinent lattice results [70,71], and in qualitative agreement with those obtained in the NJL model [88]. Everyone obtains qualitatively similar crossover of χ(T) around T c , but it would be interesting to speculate what consequences for cosmology could be if χ(T), and thus also m a (T), would abruptly fall to zero at T = T c due to a sharp phase transition of the "massless" condensate ⟨qq⟩ 0 (T). Of course, dynamical models of QCD can access only much smaller range of temperatures than lattice, where T ∼ 20 T c has already been reached [71]. (On the other hand, the thermal behavior of the U A (1) anomaly could not be accessed in chiral perturbation theory [89].) Since it is now established that (at vanishing and low density) the chiral transition is a crossover, it is important that one can use massive-quark condensates, which exhibit crossover behavior around T ∼ T c . In the present work, they give us directly, through Equation (9), the crossover behavior of χ(T). However, these are regularized condensates, because a nonvanishing current quark mass m q makes the condensate ⟨qq⟩ plagued by divergences, which must be subtracted. In Section 3, we have shown that our regularization procedure is reasonable and in good agreement with at least two widely used subtractions on the lattice.
To discuss our approach from a broader perspective, it is useful to recall that JLQCD collaboration [90] has recently pointed out how the chiral symmetry breaking and U A (1) anomaly are tied, and stressed the importance of the qq chiral condensate in that. The axion mass presently provides a simple example thereof: through Equation (11), m a (T) is at all temperatures directly expressed by the QCD topological susceptibility χ(T), which is a measure of U A (1) breaking by the axial anomaly. We calculate χ(T) through Equation (9) from the quark condensates, which in turn arise from DChSB. In addition, conversely: melting of condensates around T ∼ T c signals the restoration of the chiral symmetry. Therefore, the U A (1) symmetry breaking and restoration being driven by the chiral ones is straightforward.
The relation of χ(T) to the η ′ mass is, however, a little less straightforward [17] because it involves several other elements, but the topological susceptibility remains the main one. Since our present results on the axion are, in a way, a by-product of the framework which was initially formulated to understand better the T-dependence of η ′ and η masses, we have explained it in detail in Subsection 2.2. Therefore, here in the Summary, we should just stress that the topological susceptibility χ(T) is the strong link between the QCD axion and the η-η ′ complex. It so also in the case of the present paper regarding our η-η ′ reference [17]: specifically, we should note that the actual T-dependence of η ′ and η is rather sensitive to the behavior of χ(T), and rather accurate χ(T) is needed to get acceptable M η (T) and M η ′ (T). Thanks to its crossover behavior, our χ(T) gives in Ref. [17] empirically allowed T-dependence of the masses in the η-η ′ complex. However, even a crossover, if it were too steep, would lead to the unwanted (experimentally never seen) drop of the η mass, just as a too slow one would not yield the drop of the η ′ mass required according to some experimental analyses [65,66].
In that sense, our present predictions on m a (T) are thus supported by the fact that our calculated topological susceptibility χ(T) gives the T-dependence of the U A (1) anomaly-influenced masses of η ′ and η mesons [17] which is consistent with experimental evidence [65,66]. Acknowledgments: This work was supported in part by STSM grants from COST Actions CA15213 THOR and CA16214 PHAROS.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: At T = 0, the Dyson-Schwinger equation (DSE) approach in the rainbow-ladder approximation (RLA) tackles efficiently solving Dyson-Schwinger gap equation and Bethe-Salpeter equations, but extending this to T > 0 is technically quite difficult. We thus adopt a simple model for the strong dynamics from Ref. [77], namely the model we already used in Refs. [17,30,31,48]. For the effective gluon propagator in a Feynman-like gauge, we use the separable Ansatz: g 2 D ab µν (p − ) eff = δ ab g 2 D eff µν (p − ) → δ µν D(p 2 , 2 , p ⋅ ) δ ab , whereby the dressed quark-propagator gap Equation (13) with (14) yields B q (p 2 ) = m q + 16 3 d 4 (2π) 4 D(p 2 , 2 , p ⋅ ) More specifically, the so-called rank-2 separable interaction entails: Then, the solutions of Equations (A2)-(A3) for the dressing functions are of the form B q (p 2 ) = m q + b q F 0 (p 2 ) and A q (p 2 ) = 1 + a q F 1 (p 2 ), reducing Equations (A2)-(A3) to the nonlinear system of equations for the constants b q and a q : If one chooses that the second term in the interaction (A4) is vanishing, by simply setting to zero the second strength constant, D 1 = 0, one has a still simpler rank-1 separable Ansatz, where A q (p 2 ) = 1.
For separable interactions, the dressing functions A q , C q and B q depend only on the sum p 2 n = ω 2 n + ⃗ p 2 . In the separable models (A4), with their characteristic form (A5) of the propagator solutions at T = 0, the dressing functions obtained as solutions of the gap equation at T > 0 are: