Chiral observables and topology in hot QCD with two families of quarks

We present results on QCD with four dynamical flavors in the temperature range $150$ MeV $\lesssim T \lesssim 500$ MeV. We have performed lattice simulations with Wilson fermions at maximal twist and measured Polyakov loop, chiral condensate and disconnected susceptibility, on lattices with spacings as fine as 0.065 fm. For most observables spacing effects are below statistical errors, which enables us to identify lattice results with continuum estimates. Our estimate of the pseudocritical temperature compares favorably with continuum results from staggered and domain wall fermions, confirming that a dynamical charm does not contribute in the transition region. From the high temperature behaviour of the disconnected chiral susceptibility we infer the topological susceptibility, which encodes relevant properties of the QCD axion, a plausible Dark Matter candidate. The topological susceptibility thus measured exhibits a power-law decay for $T/T_c \gtrsim 2$, with an exponent close to the one predicted by the Dilute Instanton Gas Approximation (DIGA). Close to $T_c$ the temperature dependent effective exponent seems to approach the DIGA result from above, a behaviour which would support recent analytic calculations based on an Instantons-dyons model. These results constrain the mass of a hypothetic QCD post-inflationary axion, once an assumption concerning the relative contribution of axions to Dark Matter is made.


I. INTRODUCTION
The properties of strong interactions at high temperatures are under active scrutiny both theoretically and experimentally. From an experimental viewpoint, collisions of ultrarelativistic ions at the LHC create matter whose temperature has been estimated to reach 500 MeV and more. Theoretically, the behaviour of the theory at high temperature offers important insights into the mechanisms of chiral symmetry, confinement and the degrees of freedom active in the different phases. For instance, both experiments and lattice results indicate that the system remains very strongly coupled above the critical temperature T c , and that the perturbative regime, and eventually a free field behavior, will only set in at much higher temperatures, see e.g. Ref. [1].
As a consequence of this, fundamental degrees of freedom manifest their influence in a sequential way when temperature increases: particles of the light spectrum dissolve first, heavy bound states later. Among the quarks themselves, the light ones, which receive most of their masses at the quark-gluon transition, have a significant influence on the transition itself: the dependence of the value of the pseudocritical temperature(s) as well as of the nature of the transition on the masses of light current quark masses, as well as the role of the strange quark mass around T c have been studied at a great length, and the current status of our understanding of the sensitivity of the phase diagram to the light degrees of freedom is for instance summarized in the so-called Columbia plot [2][3][4]. Closely related with the issue of nature of the phase transition is the (effective) restoration of the U A (1) symmetry [5,6].
By contrast, we know that the charm quark mass does not carry any influence on the transition: if we were to draw a multidimensional Columbia plot including the charm mass m c the plot will just look the same for values of the charm mass ranging from infinity down to the experimental value. This can be anticipated just by considering that the charm does not acquire its mass across the chiral phase transition, and our numerical results will confirm this. However, as the temperature increases the role of the charm becomes more important, and HTL calculations predict a temperature threshold for sizeable charm effects at about 300 MeV [7].
This suggests that precise calculations of QCD thermodynamics in the LHC working region, where temperature as high as 500 MeV may be reached, require a dynamical charm. However, in comparison with studies with 2 + 1 flavors, lattice results with two dynamical families are much less developed -lattice discretization effects are more severe due to the charm large mass. Only the MILC collaboration [8] and the Wuppertal-Budapest collaboration [9] have published results for the Equation of State with N f = 2 + 1 + 1 flavors based on the staggered action. Our own study of the Equation of State is underway and only preliminary and partial results have appeared [10], which will not be presented here.
In this paper we study the chiral and topological properties of QCD for temperatures ranging from about 150 MeV till about 500 MeV, well above the deconfinement transition, and close to the charm threshold. We use four flavors of maximally twisted mass Wilson fermions in the isospin limit, i. e. with degenerate up and down quark masses, and physical strange and charm masses, relying on the zero temperature results of the ETM collaboration for scale setting and tuning of the algorithm to simulate at maximal twist [11].
The purpose of this study is twofold: firstly, we present the first results around the pseudocritical temperature obtained with twisted mass Wilson fermions. We successfully cross check our results with earlier ones obtained with N f = 2 + 1 staggered and domain wall fermions [12][13][14], thus confirming that a dynamical charm has no influence around the transition. Secondly, we will use the results for the disconnected chiral susceptibility to calculate the topological susceptibility for temperatures as high as 500 MeV, apparently reaching the onset of the Dilute Instanton Gas Approximation. Since this second aspect is at the moment under active investigation by several groups, we conclude this introduction with a brief review of the current status of topology in hot QCD.
QCD topology is an eminently non perturbative problem which has been addressed on the lattice since long. It is well known that topological studies are hampered by technical difficulties on a discrete lattice [15,16]. However, recent methodological progress, together with more adequate computer resources, have to some extent reopened the field, leading to the first results on topology at high temperature with dynamical fermions [9,[17][18][19][20]. Although these studies exhibit some common features, which we will review in the following together with our own results, a quantitative agreement has not been reached yet. Particularly significant -and still under debate -is the onset of the Dilute Instanton Gas behaviour: once this is reached, the results could be safely extrapolated to temperatures T = O(1) GeV of cosmological relevance.
In the calculation of the topological susceptibility presented here we follow an early proposal of Kogut, Lagae, and Sinclair [21] which has also been investigated by other groups [19,22,23]. In a nutshell, we will use well known identities in the fermionic sector based on the Atiyah-Singer theorem -as we will review in the followingto infer the topological susceptibility at high temperatures from the results for the chiral susceptibility. Let us mention from the start that this approach requires the approximate restoration of the U A (1) symmetry. We will not discuss in detail this important problem itself [6], rather we rely on the fact that the residual breaking above T c , if any, is negligible. Let us finally mention that a seminal paper [24] (still limited to the quenched approximation) has proposed to use lattice results on topology and the Equation of State to constrain the mass of the post-inflationary QCD axion, paving the way to more quantitative and realistic results. We will return to this point in Section V. This paper is organised as follows. In the next Section we review the lattice action and the setup of the simulations. Observables are introduced in Section III. The following two Sections contain our results: Section IV is devoted to the analysis of the pseudocritical region, while in Section V we discuss topology and its implication on the bounds on the post-inflationary QCD axion. We close with a brief discussion on the present status of N f = 2 + 1 + 1 thermodynamics and topology, and future steps.
Some of the results presented here have been presented in a preliminary form at Conferences [10,25].

II. THE ACTION AND THE SIMULATION SETUP
We performed simulations with four flavors of maximally twisted mass Wilson fermions in the isospin limit, i. e. with degenerate, larger than physical, up and down quark masses, and physical strange and charm quark masses. In terms of the twisted fields χ l,h = exp (−iπγ 5 τ 3 /4)ψ l,h the light and heavy quark twisted mass actions have the following form: and similarly: where D W [U ] is the usual Wilson operator, a is the lattice spacing, and χ l,h are quark spinors in the twisted basis. The hopping parameter κ is set to its coupling dependent critical value κ c (β) leading to the so-called "maximal twist" of the action (1)-(2) with the property of automatic O(a) improvement for expectation values  of any operator [26,27]. The parameter µ l describes the mass of the degenerate light quark doublet, which is still unphysically large in our study: the charged pion mass values m π ± considered at present are 210, 260, 370 and 470 MeV. The heavy twisted mass parameters µ σ and µ δ have been tuned in the unitary approach to reproduce approximately the physical K and D meson mass values within the accuracy of 10%, thus allowing for a realistic treatment of s and c quarks. For the gauge sector the Iwasaki action is used (c 0 = 3.648 and c 1 = −0.331): The two sums extend over all possible plaquettes (P ) and planar rectangles (R), respectively. Our finite temperature simulations have been performed for three values of β = (1.90, 1.95, 2.10). Using the nucleon mass to fix the scale, this gives a = 0.0646 fm, a = 0.0823 fm and a = 0.0936 fm [28]. For each lattice spacing we explored temperatures ranging from 150 MeV to 500 MeV by varying the temporal size of the lattice N τ . So far we have generated finite temperature configurations for eight sets of parameters that correspond to four values of the charged pion mass of about 470, 370, 260 and 210 MeV, for which two, three, two and one value(s) of the lattice spacing have been considered, respectively.
For each lattice spacing the advantage of this setup ("fixed-scale approach" [29]) is that we are allowed to rely on the setup of T = 0 simulations of ETMC, thus exploring a wide set of temperatures. Different lattice spacings allow us to study the approach to the continuum limit. In principle, we have to deal with the mismatch of different temperatures obtained for different choices of β. In practice, the temperature scans are fine enough to overcome this potential disadvantage, as we will see when presenting the results. The full set of parameters, as well as the indicative statistics, is reported in Table I.

III. OBSERVABLES
We concentrate first on the study of the critical region: our primary observables here are the quasi order parameter for deconfinement and chiral symmetry breaking -the Polyakov loop and the chiral condensate in the light sector. We also consider the disconnected chiral susceptibility: as it is merely the fluctuations of the order parameter, it has the same meaning of ordinary susceptibilities in spin models, hence it carries the relevant information on the pseudocritical behaviour. We will also discuss later how this observable acts as a proxy for the topological susceptibility in the symmetric phase.

A. Polyakov loop
The Polyakov loop is defined as a Wilson loop of gauge fields winding once around the temporal (thermal) direction. The important quasi order parameter is the real part of it, We may only consider the real part since when quarks are present the center symmetry is explicitly broken towards the real sector. In pure gauge theory this quantity is strictly the order parameter for the deconfinement transition, and since the pure gauge limit is approached for large quark masses the Polyakov loop acts as a quasi order parameter also with dynamical quarks. It may be renormalized as follows [30] Re where V (r 0 ) denotes the static quark-antiquark potential at the distance of the Sommer scale r = r 0 [31]. The latter is to be determined at zero temperature. The static quark-antiquark potential V (r) has been evaluated using 20 steps of APE smearing with a smearing parameter α APE = 0.5. For r 0 we use the values determined by ETMC. Subsequently V (r) was interpolated using cubic splines and its value V (r 0 ) has been extracted.

B. Chiral condensate
The chiral condensate ψ ψ in the light sector serves as an (approximate) order parameter for the SU (2) × SU (2) symmetry, which is explicitly broken by the quark mass: The trace of the inverse is evaluated by means of the technique of noisy estimators using multiple quark matrix inversions acting on 24 Gaussian noise vectors. At maximal twist the leading ultraviolet divergent part is proportional to µ/a 2 which in our early N f = 2 study was removed by subtracting the corresponding condensate at the same mass at T = 0 [32]. The multiplicative renormalization factors are then canceled by dividing by the condensate at zero temperature in the chiral limit: Having now the strange quark at our disposal, we are using another prescription based on a suitable subtraction (involving the light and strange masses) of the strange quark condensate. It eliminates the additive µ-proportional divergence contained in the light condensate being estimated by means of the strange quark condensate. This procedure has been proposed for N f = 2 + 1 flavors in the literature [33]: where ψ ψ s is the strange quark condensate calculated in the Osterwalder-Seiler setup [26,34] to avoid mixing in the heavy quark sector and where the strange mass µ s has been determined as to reproduce the physical mass ofsγ µ s. The corresponding expression at T = 0 serves as normalization factor.

C. Chiral susceptibility
The chiral susceptibility, defined as the mass derivative of the chiral condensate reads: The disconnected chiral susceptibility, which carries the relevant information on the critical behaviour, is the quadratic fluctuation of the chiral condensate: for which the traces arising from the path integral have been evaluated using the stochastic noise method and we made sure that no net connected piece enters the final result. Ultraviolet divergences cancel in the difference, and we are left only with a multiplicative renormalization.

IV. THERMAL TRANSITION TEMPERATURE(S)
We have considered the renormalized Polyakov loop and two chiral observables in order to extract three different pseudo-critical temperatures. We will discuss the used fit strategies in what follows.

A. From the Polyakov loop
The deconfinement transition temperature T deconf is read off the inflection point of a hyperbolic tangent function fit to the renormalized Polyakov loop data The data becomes more and more noisy for larger N τ which mostly reduces the data quality for the small mass ensembles.
We have restricted the fits to temperatures below T = 310 MeV for our central fits in order to avoid the region of visible discretisation artifacts at high temperatures. We have always included into the fits all available data at the low temperature end since data is very limited there.
A second fit was performed including higher temperature data up to T = 380 MeV (or up to T = 400 MeV in the case of B260, respectively) in order to estimate a systematic error for T deconf taking half of the deviation from the central fit results above.
In the middle panel of Figure 1 we have added the Polyakov loop data of our finite size test 24 3 ensemble B370.24. There is no difference in the data as compared to the data obtained in the larger volume 32 3 .

B. From the renormalized chiral condensate
Here the first estimate T ∆ for the crossover temperature T c is determined from the renormalized chiral condensate. The data of ∆ ls follows an S-shape curve for all considered ensembles. We therefore used as a fit ansatz the following sigmoid function: We always used all available data at low temperatures and used two upper limits for the fit ranges in temperature.
The main fits were obtained with data with T < 350 MeV and another extended fit has been done with T < 450 MeV. Half the deviation of the latter from the main fit results was used to estimate the systematic error. Here the second estimate for the crossover temperature T χ is determined from the chiral susceptibility. We noted that by considering the T χ disc ψψ rather than χ disc ψψ the signal is sharper, and we used this in our estimates. The modest shift in the pseudocritical temperature associated with this different normalization is about 10 MeV, negligible with respect to other errors.
For estimating the position of the peak in the (modified) disconnected chiral susceptibility Eq. (10) we fit an ansatz quadratic in the temperature to the data: a definition that does not assume a model of the T -dependence of the susceptibility data but is a generic fit function that should be valid in the close vicinity of a peak. We have fitted it to the peak regions of our χ disc ψψ data. Each ensemble was fitted separately and we have considered several fit regions in temperature for each ensemble to address the systematic uncertainty of its choice. We take the central values and the statistic errors from the best fit result (in terms of χ 2 /dof) and use fits with χ 2 /dof < 2.5 or with χ 2 /dof = ∞ (i. e. the parameters are fully constrained by the data). Among the allowed fit results imposing the just mentioned conditions on the fit quality we picked the one with the maximum deviation from the central fit and took half of the deviation as the systematic error.
Since we adopted a fixed scale approach the number of data points we are able to measure in the peak region is restricted by construction. It is even more restricted for the smaller pion masses since there intermediate points would be at large odd N τ which are prohibitively expensive in terms of computing time and therefore out of reach at our present possibilities. Consequently, the peak region for the small pion masses is not too well covered by data and the fits turned out to be more difficult than those to the renormalized quark condensate discussed previously.
In Figures 3-5 we show for each ensemble the best fit obtained upon varying the fit range in T . Figure 3 shows the chiral susceptibility for the three ensembles with the lightest pion. While the peak is nicely visible for the ensemble at ∼ 210 MeV pion mass the maximum of the data is much more weakly pronounced for the two ensembles at about 260 MeV pion mass which in both cases is mainly due to the smallest temperature points. While for the A260 ensemble the best fit is accordingly found discarding the smallest temperature point, the fits for the B260 ensemble favour to include it.   Figure 6 shows the obtained results for T χ and T ∆ . The estimates of T χ at physical quark masses were obtained by the Wuppertal-Budapest [12] and the HotQCD [13] collaborations with staggered quark discretisations. There is also the HotQCD result from the independent study with chiral domain wall fermions [14], which is almost identical to the value from their staggered analysis. We fitted our values of T χ (taken from fits to the chiral susceptibility) obtained at larger than physical light quark masses with the two ansaetze that have been under discussion in the context of the order of the phase transition in the chiral limit of the two-flavour theory [2,4].   A few comments are in order. First, we confirm that for these largish masses we cannot discriminate between different critical scenario. On the positive side, either plausible parameterizations interpolates the data, and fares well through the results for physical pion masses obtained with N f = 2 + 1 staggered and domain wall fermions. We thus cross check the results for the pseudocritical temperature provided by other groups, and, at the same time, confirm that a dynamical charm does not contribute to the chiral dynamics in the pseudocritical region.  6: Pseudocritical temperatures as a function of the pion mass, superimposed by the curves corresponding to first and second order chiral scenarios. The staggered value from Wuppertal-Budapest collaboration [12] and the coinciding result from both staggered [13] and domain wall [14] studies of HotQCD collaboration for physical pion mass are shown as well. Data from the disconnected chiral susceptibility T χ disc ψψ have been used in the fits. See text for discussions.

V. TOPOLOGY AND AXION'S PROPERTIES
Let us first briefly remind how axion physics and topology are intimately connected. We will then describe the results for the topological susceptibility, comment on the scaling with the pion mass, and discuss the implications on the post-inflationary axion's mass bound.

The QCD Lagrangian admits a CP violating term
where we recognize that g 2 32π 2 F a µνF µν a is the topological charge density Q(x). The θ term gives an electric dipole moment to the neutron, which is strongly constrained experimentally [35] leading to the bound θ < 10 −10 . The strong CP problem consists in explaining this unnaturally small value.
An elegant solution to the strong CP problem postulates the existence of an extra particle [36][37][38], a pseudo-Goldstone boson of the spontaneously broken Peccei-Quinn symmetry, which couples to the QCD topological charge, with a coupling suppressed by a scale f a . The thermal Grand Canonical partition function of QCD is now a function of θ and the temperature T : and the free energy F (θ, T ) is the axion potential. At leading order in 1/f a -well justified as f a 4 × 10 8 GeV -the axion can be treated as an external source, and its mass is given by The cumulants C n of the topological charge distribution are related to the Taylor coefficients of the expansion of the free energy around θ = 0 hence higher order cumulants and their ratios carry information on the axion's interactions. We will consider higher order cumulants in a companion paper based on the gluonic measurements [17,39], while in this paper we will use our results on topological susceptibility exclusively to constrain the (post-inflationary) axion mass.

A. Topological susceptibility
We measure the topological susceptibility following Refs. [19,21,22]. These authors start from the continuum indentity and note that this remains true on a smooth gauge field configuration (i.e. if lattice artifacts are small) hence, by squaring Eq. (18), It is well known that χ disc 5 suffers from huge fluctuations. Rather than attempting its measurement, one considers that when both SU (2) × SU (2) and U A (1) are restored, χ disc 5 ≈ χ disc ψψ . We are well aware that the issue of the (effective) restoration of the axial symmetry is a subtle and unsettled one. On the other hand, by making a virtue of this, we can with some confidence assume that a residual U A (1) breaking above T c , if any, is small [5,6]. More precisely, even those studies advocating a residual breaking at T c [5] suggest that the breaking disappears at T 1.5T c . We will take this as the threshold beyond which χψ ψ /T 2 500 400 300 200 100  The main quantity of interest for our topological analysis thus turns out to be the disconnected susceptibility of the chiral condensate, χ disc ψψ = which we have already discussed within the pseudocritical region. For our topological study we will consider results in the extended temperature range, and based on these results we evaluate the topological susceptibility which we show in Figure 7. We have grouped the results according the mass value, and in each plot we show results for the different lattice spacing. Within our statistical errors we cannot see any lattice spacing dependence, and we will take our results as continuum estimates.
In Figure 8 we show all the results for the topological susceptibility on a log-log scale in the high temperature region. Superimposed we show the fits to a simple power-law fall off χ top (T ) = AT −d , which describe the data for T > 350 MeV, with a power roughly independent on the mass. As already stated, we could not detect any lattice spacing dependence within our accuracy, hence we globally fit all data belonging to the same mass for all values of lattice spacing.
The fit parameters A and d are strongly correlated and fit errors have little meaning -we merely quote the central values for d = (6.26, 6.88, 7.52, 7.48) for decreasing values of the mass m = (470, 370, 260, 220) MeV. To get a feeling on the error of the exponent d we may explore the functional dependence of the topological susceptibility on the temperature by defining a local effective power [40] d eff (T ) = T d log χ top (T )/dT which we show in Figure 9. In doing so, we do not distinguish different masses and spacings: as the log derivative is rather noisy, it is not possible to detect any trend with mass and spacing and we simply plot all the results together. The average result is broadly consistent with the fits, and apparently approaches a constant value above T ∼ 350 MeV, while the spread of the exponents could be taken as an estimate of the error on the power d of the fall-off.
Concerning the apparent constant asymptotic behaviour, we expect that at very large temperatures the partition function is described by a dilute gas of instantons and anti-instantons (DIGA), which leads to a well defined prediction for the fall off of the topological susceptibility, also shown in Figure 9. Apparently the DIGA result is approached already for these temperatures: however, it is not clear from the plot whether such decreasing trend with temperature has reached its nearly asymptotic value (modulo small corrections), or, rather, whether the apparent coincidence with the DIGA value is accidental, and restricted to a limited range of temperatures. Only simulations for larger temperatures can settle this issue.
Further, we note that the DIGA result is approached from above: the exponent has a larger value for smaller temperature -see again Figure 8. Very interestingly, it has recently been proposed [41] that close to T c the dilute instanton gas changes to an ensemble of instantons and dyons, and the signature for this should be a faster decrease of the topological susceptibility close to T c . It remains to be seen whether a hypothetical residual breaking of U A (1) is affecting the results close to T c -in which case the decreasing trend would reflect the progressive U A (1) restoration, rather than the dyons dynamics suggested in [41]. Clearly clarifying the fate of the U A (1) symmetry is important also in this context -an early restoration of U A (1) would validate the results for the topological susceptibility also around T c .

B. Pion mass dependence and rescaling
Our results have been obtained with pion masses ranging from 220 till 470 MeV, above the physical pion mass. The same dilute instanton gas model has a prediction for the mass scaling which may be used to extrapolate our results to the physical pion mass: χ top ∝ m 4 π . We note that this leading mass dependence is more general than DIGA and simply follows from the analyticity of the chiral condensate in the chiral limit above T c . In fact, taking ψ ψ = n=0 a n m 2n+1 l in the symmetric phase, the total susceptibility is an even series in the quark mass Barring unexpected cancellation we may assume that the same holds for the connected and disconnected susceptibilities separately, hence An exact DIGA form would imply that the leading order is exact, i.e. that the disconnected chiral susceptibility does not depend on the pion mass in the mass range considered. Our results do show a mass dependence, not surprisingly, given the largish masses which we are using -smaller masses are probably needed to get rid of the subleading mass corrections. In this first study, rather than attempting an extrapolation in mass we content ourselves with the simple rescaling dictated by the leading order, as done first in [18]. In Figure 10 we present the results for the topological susceptibility obtained by rescaling the data for different pion masses to the physical pion mass according to the leading scaling prescription. In the same diagram we reproduce the results from Table S7 of Ref. [9] obtained with physical quark masses for a comparison. As we have already discussed, the trend close to T c is different, while there is a broad agreement at high temperatures -we will discuss below the implications of the residual differences on the axion mass. A similar agreement holds with the results of Ref. [19].

C. Axions
A power-law decay for the topological susceptibility as noted by many authors opens a very interesting possibility: a safe extrapolation to very large temperatures. This feature has been exploited in applications to axion physics, briefly reminded in the Introduction. In a nutshell (see e.g Refs. [42]), when the axion mass is of the order of the inverse of the Hubble parameter, the axion starts to oscillate: 3H(T ) = m a (T ) = χ top (T )/f a . f a can be traded for the zero temperature topological susceptibility and the axion mass via χ top = f 2 a m 2 a , hence the knowledge of the temperature dependence of the topological susceptibility suffices to determine the time of  π . We also superimpose the tabulated results (Table S7) from [9].
the beginning of oscillation. At this time, the energy density ρ a (T ) of the oscillating axion field is the same as a collection of axions at rest ρ a (T ) ≈ 1/2m 2 a (T )f 2 a θ 2 , and the number density n a (T ) = ρ a (T )/m a (T ) can be estimated as n a (T ) ≈ 1/2m a (T )f 2 a θ 2 . The axion-to-entropy ratio remains constant after the beginning of the oscillations, so the present mass density of axions is ρ a,0 = n a s m a s 0 , where s, s 0 are the entropies at time T , and of today, and Ω a = ρ a,0 ρ c , with ρ c the critical density. To obtain simple expressions in closed form we employ the power-law parameterization g * (T ) = 50.8 (T /(MeV)) 0.053 for the number of relativistic degrees of freedom g * (T ) entering the Hubble parameter and entropy density, which reproduces the results up to a few percent in the temperature interval 800 < T < 1500 MeV. Following e.g. [42] we finally arrive at ρ a (m a ) ∝ m − 3.053+d/2 2.027+d/2 a , where d > 0 defines the power law decay of the topological susceptibility at high temperatures discussed above, χ top A T −d , and we use the latest PDG results for the required astrophysical constants [44]. We present the results in a graphic form in Figure 11: we plot the axion's fractional contribution to Dark Matter versus the axion mass for various situations. Similar analyses have been presented in previous works, and we confirm to a large extent their conclusions [9,19,43]. The first three lines are obtained from our results, rescaled to the physical pion mass. We have included errors for the lowest mass, which are not visible in the graph. The discrepancy between the results for the two lowest pion masses, and for a pion of 370 MeV may be ascribed to violations to the mass scaling discussed above, and are anyway small at practical level. Clearly, as we know, further uncertainties might be hidden in the results -to estimate their impact we plot a few mock curves on the basis of 210 MeV data: first, we keep the same amplitudes as the one we have measured, but choose the DIGA exponent. Second, we choose a small exponent similar to the one reported in [18]. To study the effect of the amplitude, we fix the exponent to the one we have measured at 210 MeV, and consider the fourth root of the amplitude ten times larger or smaller. In all cases the intercept with the abscissa axis (overclosure bound) defines the absolute lower bound for the axion mass: for our results, this gives m a 20µeV. This limit becomes more stringent if we assume that the axions only contribute to a fraction of DM -the results can be read off the plot from the intercept of the desired fraction and the corresponding curves. Clearly the bounds are robust against 'small' changes of parameters. However a significant variability remains. In particular, as also noted in Ref. [45], slower decays (as those observed by us using the gluonic definition of the topological susceptibility) would considerably lower the axion bound.

VI. DISCUSSION
We have studied chiral and topological properties of N f = 2 + 1 + 1 QCD with twisted mass Wilson fermions. The strange and charm mass have their physical values, while for the light pion we have considered four different masses.
We have identified the pseudocritical temperature as a function of the pion mass. We have found that an extrapolation to the pseudocritical temperature for physical pion is robust with respect to different assumption for the universality class of the two flavor massless theory, and agrees well with previous estimates. This serves as a sanity check of for the twisted mass chiral dynamics around T c and in addition it confirms the irrelevance of a dynamical charm in the transition region.
We have measured the topological susceptibility in the range 150−500 MeV and found a rather fast decrease with temperature in the range T c < T 350 MeV. This feature is predicted in recent instanton-dyon models [41], however it has not been observed in other recent lattice studies. Such fast decrease may also been understood within the framework of the QCD magnetic Equation of State [46], as it is known that around T c the disconnected susceptibility almost saturates the total susceptibility ∂ ψ ψ /∂m. Since the total susceptibility in the symmetric phase behaves as 1/(T − T c ) γ with γ 1 (with a weak mass dependence due to Griffith analyticity) it is rather natural to expect a fast decrease of the topological susceptibility leading to an apparently large exponent in a simple power-law parameterization T α . This scenario should be checked by directly computing the regular contribution to the susceptibility, and most importantly the γ 5 susceptibility, which would allow the computation of the topological susceptibility without relying on the restoration of the U (1) axial symmetry. At higher temperatures the contribution from the regular part becomes significant and the results approach the DIGA behaviour from above -although, as we have already mentioned it is hard to exclude a continual decreasing trend which would bring the results well below that.
This work might be extended along several directions: firstly -and obviously -the results should be further pushed towards lower masses. The first zero temperature results by the ETM collaboration with N f = 2 + 1 + 1 and a physical pion mass have appeared recently [47], and we hope to be able to extend them to high temperature. This will be particularly relevant in the context of the restoration of the U A (1) symmetry. Second, we would