On the strength of the $U_A(1)$ anomaly at the chiral phase transition in $N_f=2$ QCD

We study the thermal transition of QCD with two degenerate light flavours by lattice simulations using $O(a)$-improved Wilson quarks. Temperature scans are performed at a fixed value of $N_t = (aT)^{-1}=16$, where $a$ is the lattice spacing and $T$ the temperature, at three fixed zero-temperature pion masses between 200 MeV and 540 MeV. In this range we find that the transition is consistent with a broad crossover. As a probe of the restoration of chiral symmetry, we study the static screening spectrum. We observe a degeneracy between the transverse isovector vector and axial-vector channels starting from the transition temperature. Particularly striking is the strong reduction of the splitting between isovector scalar and pseudoscalar screening masses around the chiral phase transition by at least a factor of three compared to its value at zero temperature. In fact, the splitting is consistent with zero within our uncertainties. This disfavours a chiral phase transition in the $O(4)$ universality class.


Introduction
Nuclear matter under extreme conditions of high temperatures T and/or baryon chemical potential µ B is the subject of intense experimental and theoretical studies in nuclear, particle and astro-physics. One of the salient features of strongly interacting matter is the high-temperature transition from the hadronic phase to the deconfined quark-gluon plasma (QGP). The transition takes place in a temperature regime between 100 and 300 MeV, where the QCD running coupling is strong. Thus, a non-perturbative investigation of the transition and the properties of the QGP is necessary and a lot of effort has been invested by the lattice community in the past decades (for recent reviews see [1][2][3][4][5]). Because lattice studies of the QCD transition at finite baryon chemical potentials are severely hampered by the sign problem, the QCD phase diagram remains largely unknown. Even at zero baryon density, the nature of the thermal transition with light quark masses approaching the chiral limit is not yet determined in the continuum. Knowledge of this important limit would also help to constrain the phase diagram at non-zero µ B . Fig. 1 summarises the current knowledge about the order of the thermal transition for vanishing baryon density in the (m ud , m s )-plane, where m ud is the mass of the degenerate up and down quarks and m s the strange quark mass. In the opposite limits of pure gauge theory and QCD with three massless quarks, there are true first-order phase transitions associated with the breaking of centre symmetry [6], and the restoration of the SU(3) chiral symmetry [7], respectively. These get weakened by the explicit breaking of those symmetries by finite fermion masses, until they disappear along second order critical lines. For intermediate quark masses, the finite temperature transition is then merely an analytic crossover.
There is plenty of evidence that the physical quark mass configuration realised in nature is in the crossover region. Early results based on the staggered fermion discretisation [8,9] have been confirmed by domain wall fermions [10,11] and simulations with Wilson fermions are approaching the physical point as well [12][13][14][15]. The critical line separating the first order chiral transitions from the crossover region, the chiral critical line, has been mapped out on coarse lattices and is in the Z(2) universality class of the 3d Ising model [16,17]. The critical line in the heavy quark region, the deconfinement critical line, is in the same universality class and was mapped out on coarse lattices simulating a hopping expanded determinant [18] and a 3d effective lattice theory [19]. However, the location of the critical lines in the quark mass plane is subject to severe cut-off effects. With standard staggered fermions the N f = 3 critical pion mass is at around two to three times the physical quark mass [16,17], yet one finds that it shrinks to nearly half that value when going from N t = 4 to N t = 6 [20]. Improved staggered fermions can only give an upper bound for the critical mass which is around 0.1 times the physical quark mass [21,22]. First results with Wilson fermions on the other hand see the N f = 3 critical endpoint at around five times the physical quark mass [23]. While the latter result will still change when going towards finer lattices, it highlights the importance of taking the continuum limit before discussing critical behaviour.
While all current results indicate that the critical line passes the physical point to the left, its detailed continuation is still largely unknown [2]. There are two possible scenarios [7,24,25]: in scenario (1), depicted in the left panel of Fig. 1, the chiral critical line reaches the m ud = 0 axis at some tri-critical point m s = m tric s , implying a second order transition for N f = 2. In the alternative scenario (2) the chiral critical line never reaches the m ud = 0 axis, so that the chiral transition at m ud = 0 is first order for all values of the strange quark mass.  Many past studies have investigated the nature of the N f = 2 transition using staggered [26][27][28][29][30][31][32], O(a) improved [33][34][35] or twisted mass [36] Wilson fermions and domain wall and overlap fermions [37][38][39][40], without being able to provide a conclusive answer in the chiral limit. The main problem in investigating scaling properties is the similarity of the critical exponents of the universality classes in question (cf. section 2.4.1). A new method to exploit the known tricritical scaling when coming from the plane of imaginary chemical potential has been proposed in [41]. First results on coarse lattices with staggered [42] and Wilson fermions [43] agree on scenario (2), but show enormous quantitative discrepancies. A similar approach proposed recently is to look at the extension of the chiral critical line in the plane with N F additional degenerate heavy quarks [44,45].
Which of these scenarios is realised depends crucially on the strength of the anomalous breaking of the U A (1) symmetry at the critical temperature in the chiral limit when the mass of the strange quark is sent to infinity, i.e. in the N f = 2 case. If the breaking of U A (1) remains strong, the transition will be of second order in the O(4) universality class [7], so that scenario (1) will be realised. If, on the other hand, the symmetry is 'sufficiently' restored (see the discussion below about different restoration criteria), either scenario is possible [24,25]. For scenario (1) the breaking would then likely be in the U(2) × U(2) → U(2) universality class [24,25] (another symmetry breaking pattern/universality class of the form SU(2) × SU(2) × Z 4 → SU (2) has also been proposed [46,47]) and the O(4) universality class would be disfavoured. To be able to distinguish between the different scenarios by looking at the U A (1) symmetry, it is thus crucial to have a measure for the strength of the breaking.
First, we recall that the U A (1) classical symmetry does not imply the existence of a conserved current: the divergence of the singlet axial current A 0 µ (x) is proportional to the gluonic operator GG = µνρσ G µν G ρσ in the chiral limit, an equality valid in any on-shell correlation function. For instance, the static correlator (∂ 2 /∂x 2 3 ) dx 0 dx 1 dx 2 A 0 3 (x)A 0 3 (0) is proportional to the corresponding static two-point function of GG, which is certainly non-vanishing at any finite temperature. By contrast, the corresponding two-point function of the non-singlet axial current A a 3 (x) vanishes in the chiral limit. In this particular sense, the U A (1) symmetry is never restored.
In this paper we use static correlators and the associated screening spectrum to probe the U A (1) effects. A thermal state with an exactly restored U A (1) symmetry implies that the correlators ψ (x)τ a ψ(x)ψ(0)τ a ψ(0) and ψ (x)γ 5 τ a ψ(x)ψ(0)γ 5 τ a ψ(0) are equal in the massless theory. However, we expect the restoration of U A (1) symmetry in this sense to only be partial, improving as the temperature is increased (see e.g. [48]).
In the literature, the restoration of U A (1) symmetry has also been discussed in a different, more restrictive sense. The observables considered are e.g. the correlation functions mentioned in the previous paragraph projected onto zero-momentum ( d 4 x). Similar to the Banks-Casher relation for the chiral condensate, the difference of the zeromomentum correlators can be expressed in terms of the spectral density of the Dirac operator alone [49,50]. For a density ρ(λ) ∼ |λ| α with α > 1, the difference vanishes exactly in the chiral limit. In such a scenario, it is said that U A (1) symmetry is restored. More generally, if the SU(2) symmetry is restored at T C (which we will find to be the case later), the effective restoration of the U A (1) symmetry is indicated by the degeneracy of correlation functions belonging to a U(2) × U(2) multiplet (e.g. [49]). Using the restoration of the SU(2) symmetry, it was shown by Cohen [49,50] that the degeneracy of zero-momentum correlation functions in the multiplets is directly linked to the eigenvalue density ρ(λ) in the vicinity of zero. In particular, if the spectrum of the Dirac operator develops a gap around λ = 0 at T C , the U A (1) symmetry becomes effectively restored. Using QCD inequalities, it has been argued that the U A (1) symmetry is expected to be effectively restored as soon as the SU(2) symmetry is intact [49]. Later it was noted that the result using the inequalities was incorrect since the contributions from sectors with non-zero topology had not been taken into account properly [51,52]. However, the nonzero topology sectors only contribute away from the thermodynamic limit. Only a bit later it was shown that the eigenvalue density in the chiral limit behaves like ρ(λ) ∼ |λ| α with α > 1 [50]. Using mild assumptions and Ward-Takahashi identities for higher order susceptibilities in the framework of overlap fermions on the lattice, it has recently been shown that in fact α > 2 [53], meaning that not only the eigenvalue density, but also its first and second derivative vanish at the origin, speaking strongly in favour of a restoration of the U A (1) symmetry.
The relation between degeneracy of correlators and the behaviour of the low modes of the Dirac operator triggered a number of numerical studies of the eigenvalue spectrum in the vicinity of T C [10,[37][38][39][40]54]. Some groups [37,38] see a restoration of U A (1) at T C in the chiral limit (in particular, in [37] the behaviour of ρ(λ) ∼ |λ| 3 was observed), while others claim that U A (1) will still be broken in the chiral limit [10,54]. The difference between these studies is (apart from lattice spacings and volumes) the fermion action in use, in particular, how well the actions preserve chiral symmetry. In fact [53], to be able to show that the breaking of U A (1) still affects zero-momentum correlation functions via the low eigenvalues of the Dirac operator it is mandatory to fulfill the requirements of: (a) restored chiral symmetry on the lattice; (b) extrapolation to the infinite volume limit; (c) extrapolation to the chiral limit. In particular, condition (a) appears to play a crucial role. The reason is that any explicit breaking of the symmetry in the chiral limit by the lattice interferes with the effective restoration. Furthermore, it is mandatory to use the same fermion action also for the sea quarks, as shown in [39,40]. Here the authors looked at the small eigenvalues with domain-wall fermions and overlap fermions on the domain-wall ensembles and observed a broken U A (1) symmetry, actually made worse by the use of "quenched" overlap quarks. Only after reweighting of the configurations to the overlap ensemble an actual restoration of U A (1) in the chiral limit was observed. A possible explanation for this effect is that the "quenching" of the overlap operator leads to the appearance of non-physical near zero modes in the overlap spectrum, much like the appearance of exceptional configurations in quenched QCD. Indeed, the cases where a residual breaking was observed were those where either or both, valence and sea quarks, might not have a fully restored chiral symmetry on the lattice. Similar conclusions have been found using chiral susceptibilities, which can also be related to the eigenvalue spectrum, e.g. [10,40,55]. However, when computed on the lattice, the susceptibilities suffer from contact terms, which need to be carefully subtracted to obtain conclusive results.
In this article we present a study of the phase transition in two-flavour QCD using nonperturbatively O(a) improved Wilson fermions [56] and the Wilson plaquette action [57]. We work with a large temporal lattice extent of N t = 16 throughout, which at the chiral transition corresponds to a lattice spacing a ≈ 0.07 to 0.08 fm. Our pion masses range from about 200 to 500 MeV. In particular, we study the pseudo-critical temperatures defined by the change in the Polyakov loop and the chiral condensate, pertaining to deconfinement and chiral symmetry restoration, respectively, and check for the associated scaling in the approach to the chiral limit. As already discussed in [58], such a scaling analysis is not sufficient to distinguish between the universality classes in question. We thus direct our attention to the strength of the U A (1) breaking by investigating the degeneracy pattern of screening masses. This is complementary to other studies of the U A (1) symmetry in the literature described above, e.g. [10,11,[38][39][40], which are based on the eigenvalue structure of the Dirac operator. The screening masses probe the long-distance properties the correlators and are free of contaminations from contact terms, unlike chiral susceptibilities. We propose a measure for the strength of the U A (1)-breaking in the vicinity of T C and extrapolate it to the chiral limit. There we find it to be consistent with zero and 3 standard deviations away from its non-zero value at zero temperature. This suggests an effective restoration of the U A (1) symmetry around the critical temperature and thus a strengthening of the chiral transition for the lattice spacing considered.
As discussed above, at finite lattice spacing an exact chiral symmetry is mandatory in order to study the eigenvalue spectrum of the Dirac operator reliably. In this study we use O(a) improved Wilson fermions. While the action breaks chiral symmetry at finite lattice spacing, the static screening masses that we study approach their continuum limit with O(a 2 ) corrections. Therefore, as long as we work in a regime where these corrections are small compared to the physical mass splittings induced by the U A (1) anomaly, we should obtain qualitatively correct conclusions. If, at a given lattice spacing, a U A (1)-breaking mass splitting turns out to be small, a continuum extrapolation is required to determine how small exactly the splitting is.
Parts of our results have already been presented at conferences [58][59][60][61] and were used to investigate the properties of the pion quasiparticle in the vicinity of the transition [62][63][64][65].
The article is organised as follows: In the next section we introduce our observables, the details of our simulations and discuss the renormalisation and scale-setting procedures. In section 3, we present the numerical results. We first discuss the extraction of the pseudocritical temperatures in section 3.2 and try to compare the results to the scaling predictions in the approach to the chiral limit. We also compare our results with those from different fermion discretisations in the literature. In section 3.3 we discuss the screening masses in the different channels, before we come to the investigation of the strength of the breaking in the chiral limit in section 3.4. Finally we present our conclusions in section 4. Detailed tables collecting simulations parameters and results can be found in the appendices.

Simulation and scan setup
Our simulations are performed using two flavours of non-perturbatively O(a) improved Wilson fermions [56] and the unimproved Wilson plaquette action [57]. We use the clover coefficient determined non-perturbatively in Ref. [66]. The simulations are done employing deflation accelerated versions of the Schwarz [67,68] (DD) and mass [69] (MP) preconditioned algorithms, the latter in the implementation of Ref. [70]. Both algorithms make use of the Schwarz preconditioned and deflation accelerated solver introduced in [71,72]. As discussed in detail appendix A.1, the algorithms offer a significant speedup for large volumes and low quark masses, but also pose constraints on the available lattice sizes.
In general there are two known procedures to vary the temperature T = 1/(N t a). The first option is to vary the temporal extent whilst keeping the lattice spacing fixed. The advantage of this procedure is that all physical parameters and renormalisation constants remain fixed, making it the optimal tool for spectroscopy (see [73], for instance). The disadvantage of this procedure is that the resolution around T C is limited, made even worse by the use of improved algorithms (cf. appendix A.1). The second option, used in this study, is to vary the lattice spacing a by varying the coupling β = 6/g 2 0 , known as βscans. This offers the possibility to obtain a fine resolution around T C , but requires a good tuning of the bare quark mass to scan along lines of constant physics (LCPs), as well as the interpolation of quantities needed for scale setting and renormalisation. This is particularly demanding for Wilson fermions due to the additive quark mass renormalisation.
We will use throughout a comparatively large temporal extent of N t = 16 for two reasons. First, Wilson fermions break chiral symmetry explicitly at finite lattice spacing. Being as close as possible to the continuum helps to reduce the resulting effects as much as possible. Second, for our choice of lattice action the non-perturbative determination of the improvement coefficient c SW in the two-flavour theory extends only down to β = 5.2 [74]. Since a ≈ 0.08 fm at β = 5.2 this means that N t = 16 is necessary to allow for scans in the desired temperature range.

Scale setting and lines of constant physics
To convert our results to physical units we use the Sommer scale [75] r 0 with the interpolation of the CLS results from [74] discussed in appendix A.3. To convert to physical units we use the continuum result r 0 = 0.503(10) fm [74]. The temperature scans are done along LCPs, for which we estimate the values for the bare parameter κ corresponding to a particular quark mass by inverting the analytic relation m ud (β, κ) discussed in appendix A.4. We test the validity of this relation a posteriori by computing m ud along the β-scan. Conventionally, quark masses will be quoted in the MS scheme at a renormalisation scale of 2 GeV.
Whenever we quote pion masses for our temperature scans, we imply that these are zero temperature pion masses which correspond to the quark masses of the respective ensemble. We estimate the pion masses from our results for m ud using chiral perturbation theory to next-to-next-to leading order as given in [76]. For this we use the low-energy constants from [77] obtained by the fit denoted as 'NNLO F π , m 2 π ' with a mass cut of 560 MeV. The associated low-energy constants are given in table 6 of [77]. This procedure serves the purpose of enabling comparisons with the literature and should not be taken as a precision computation of the zero temperature pion mass.

Deconfinement and the Polyakov loop
To investigate the deconfinement properties of the transition we look at the associated order parameter, the Polyakov loop a nonzero value of which signals the spontaneous breaking of centre symmetry. Dynamical fermions explicitly break the centre symmetry of the gauge action and favour the centre sector of the Polyakov loop on the real axis, so that it is sufficient to look at Re(L) . In [78] it was found that the use of smeared links can enhance the signal in investigations of phase transitions. We have thus also computed the real part of the APE-smeared [79] Polyakov loop, Re(L) S , using 5 steps of APE smearing with a parameter of 0.5 multiplying the staples. The Polyakov loop susceptibility is given by and similarly for the smeared Polyakov loop. In order to have a quantity with a well defined continuum limit the Polyakov loop requires multiplicative renormalisation [80]. Here we will ignore this issue and work with the unrenormalised Polyakov loop. Since the renormalisation factor is expected to behave monotonically at the β values corresponding to the critical region, we do not expect the typical S-shape of the Polyakov loop vs. temperature graph to be affected.

The chiral condensate
A second aspect of the transition to the quark gluon plasma is the restoration of chiral symmetry. In the chiral limit, the associated order parameter is the chiral condensate It governs the response of the system with respect to the external 'field' which breaks the symmetry explicitly, i.e., the quark mass m ud . The bare chiral condensate is given by where D is the Dirac operator and the expectation value on the right-hand side is taken with respect to the gauge field. The associated susceptibility consists of a disconnected (the terms in the curly brackets) and a connected part, In the region around T c the disconnected part has been found to dominate the transition signal in the susceptibility [9,30,81]. Close to the chiral limit, however, this statement does not necessarily hold. The connected part only receives contributions from isovector states, while the disconnected part receives contributions both from isovector and isoscalar states.
Since an unbroken U A (1) symmetry would imply light isovector scalar states, the relevant magnitude of the two contributions is an important indicator of the nature of the transition.
Here we will focus on the disconnected part of the susceptibility and leave the comparison between the connected and the disconnected susceptibility for future publications. Due to mixings with operators of lower dimension, the condensate contains cubic, quadratic and linear divergences, and therefore requires additive renormalisation [82,83]. In addition, it renormalises multiplicative with the renormalisation factor associated with the scalar density, Z S , which is equivalent to the inverse of the mass renormalisation factor for Wilson fermions, Z S = Z −1 m [84], where Z m is defined in appendix A.4. Neither the additive nor the multiplicative renormalisation depends on the temperature, so that we can cancel the divergent terms by subtracting the chiral condensate at T = 0. The associated difference renormalises multiplicatively, For the determination of Z m we can use eqs. (A.5), (A.6) and (A.7) to obtain . Z PCAC (β) can be taken from the fit in appendix A.4. Using axial Ward identities (AWIs) one can also define an observable which is free of cubic divergences and reproduces the chiral condensate in the chiral limit [83,84]. The axial Ward identity in the form integrated over spacetime (and up to a contact term at y = x) is given by where m PCAC is the bare PCAC quark mass and b 0 represents the cubic divergence in the bare condensate. At finite quark mass the second term on the r.h.s. vanishes due to the absence of a true Goldstone boson [84] and we can use the first term as the definition of a bare subtracted condensate ψ ψ bare sub still suffers from additive quadratic and linear divergences and renormalises multiplicatively with Z P [84]. Again we can subtract the residual additive divergences using the T = 0 counterpart, so that we obtain an alternative renormalised chiral condensate We define the susceptibility of the quantity (2.10) as is not equivalent to the disconnected chiral susceptibility in eq. (2.6), it is expected to show a peak at the position of the chiral transition. The subtraction of the condensates at T = 0 requires the measurement of ψ ψ bare and P P on zero temperature ensembles. Here we use the set of N f = 2 ensembles generated within the CLS effort and the interpolation discussed in appendix A.5. Table 1. Bilinear operatorsψΓψ used for the screening correlators. Here i = 1, 2 for screening masses in x 3 -direction.

Mesonic correlation functions and screening masses
Mesonic correlation functions are a valuable probe of the properties of the QGP [85,86].
be the correlation function of two operators X and Y . The equality of two correlation functions in channels of different quantum numbers signals the restoration of the associated chiral symmetry. Here x µ is the coordinate of the direction in which the correlation function is evaluated and x ⊥ is the coordinate vector in the orthogonal subspace. The isovector correlation functions of interest for the chiral transition are the vector (V ) vs. axial-vector (A), and the pseudoscalar (P ) vs. scalar (S) channels, related by We choose the isovector channels as observables because they are free of disconnected diagrams, and the correlation functions can therefore be obtained with greater accuracy. The bilinear operators for the different channels are listed in table 1. While temporal correlation functions C XY (x 0 ) can be related to the real-time spectral densities (see [87]), here we are interested in spatial correlation functions C XY (x 3 ), which are related to the screening states of the plasma. In particular, the leading exponential decay of the correlator C XX (x 3 ) defines the lowest-lying 'screening mass' M X associated with the quantum numbers of the operator X. Screening masses can be interpreted as the inverse length scale over which a perturbation is screened by the plasma. If a symmetry imposes the equality of two correlation functions, it must also imply the degeneracy of the corresponding screening masses. The latter are thus quantities sensitive to the restoration of the symmetry. Consequently, the screening masses in the V and A channels provide an alternative way of defining the chiral symmetry restoration temperature. In contrast to susceptibilities, defined by the integrated correlation function, screening masses probe the long-distance properties of the correlation functions and thus do not suffer from contact terms.
Apart from their relation to chiral symmetry, mesonic screening masses are valuable quantities to probe the medium effects of the plasma and, at high temperatures, to test the applicability of perturbation theory. They have been studied in lattice QCD for a long time (for a review of early results see [88] and for more recent studies [89,90]). In the hightemperature limit, all screening masses approach the asymptotic value M ∞ = 2πT [91,92]. The leading order correction from the interactions has been computed in perturbation theory and is known to be positive [93,94]. Static and non-static screening masses can also be computed within an effective theory approach and provide an indirect probe for real-time physics in the Euclidean lattice setup [95,96].

Critical scaling
The main question driving the present study is the nature of the phase transition in the chiral limit. Simulations with vanishing quark masses are currently impossible; in order to extract information on the order of the transition, it is customary to investigate the scaling of various observables in the approach to the critical point (0, 0) in the parameter space of reduced temperature τ = (T − T c )/T c and external field h. The scaling laws can be derived from the scaling of the free energy F (see [97]). The variable playing the role of the external field depends on the particular scenario (see section 1). If the second order scenario (scenario (1)) is realised, no matter whether the universality class is O(4) or the one from the U(2)×U(2)→U(2) scenario 1 , the critical point is located in the chiral limit and the chiral condensate constitutes a true order parameter. In this case the external field h is proportional to the (renormalised) quark mass m ud . In the first-order scenario, depicted in the right panel of figure 1, the critical point is located at a finite quark mass m cr ud , so that chiral symmetry is broken explicitly at the critical point. One must in general expect that the external field is given by a linear combination of δm = m ud − m cr ud and τ . Furthermore, ψ ψ no longer constitutes a true order parameter. The situation is analogous to the approach of the chiral critical line along the N f = 3 axis [16].
Here we will perform an analysis based on the scaling of the order parameter or the transition temperature T C with the external field. In the vicinity of the critical point a true order parameter Θ satisfies the scaling relation (see e.g. [9,33,97,98]) Here f is a function depending on the universality class of the transition and s.v. stands for scaling violations which constitute terms that are regular in τ [9,97,98]. A number of studies have looked at the scaling of the chiral condensate in the approach to the chiral limit [9,30,33,36,98] and found consistency with O(4) scaling. From the scaling in eq. (2.16) one can also derive a scaling law for the critical temperature as a function of the external field [26]. The resulting scaling relation is where C is an unknown constant. It must be stressed that these scaling laws are only valid after the thermodynamic and continuum limits have been taken. Another problem for any study of the scaling laws is the similarity of the critical exponents in the three potentially relevant universality classes. They are summarised in  Table 2. Critical exponents for the universality classes (UC) relevant for the chiral transition. The critical exponents of the U(2)×U(2)→U(2) universality class, we have taken the results from [25] (from the MS scheme). The first error is statistical while the second quoted error denotes a systematic uncertainty arising from the scheme dependence. The critical exponents of the U(2)×U(2)→U(2) universality class have also been obtained recently using the bootstrap method [100].
is so small that very accurate results are needed to be able to distinguish between the two. Thus one cannot draw conclusions from the agreement of lattice data with the scaling of one universality class alone; instead one needs to demonstrate the ability to distinguish between the scenarios.

U A (1) symmetry restoration
The strength of the anomalous breaking of the U A (1) symmetry at the transition temperature in the chiral limit is thought to be crucial for the order of the chiral transition [7,101]. However, this raises the question of how to quantify the strength of the U A (1)-breaking. As a possible reference value we suggest the screening mass gap between the isovector pseudoscalar and scalar channels at T = 0, since the pion mass vanishes in the chiral limit. Ultimately, one would like to obtain the chirally extrapolated value of m a0 from lattice QCD, since this would give a result valid for the N f = 2 case in the range of relevant lattice spacings. Unfortunately, the scalar correlation function in the iso-vector channel is rather noisy, so that a reliable extraction is currently not possible. We will discuss a phenomenological estimate for the chiral limit in section 3.4. We note that in the two-flavour theory, the a 0 meson is expected to be stable or almost stable, since theKK and ηπ decay channels known from experiment are missing. Indeed, in N f = 2 QCD only a flavour-singlet pseudoscalar exists, sometimes called η 2 , whose nature is closer to the physical η meson, and whose mass has been estimated at about 800MeV at the physical pion mass [102]. The lightest isovector scalar state was found to lie between the physical a 0 (980) and a 0 (1450) states [103]. The splitting between the pion and the lightest isovector scalar state thus provides a convenient measure for the breaking of U A (1).

Ensembles and measurement setup
In this paper we present results for the chiral transition obtained from the scans on 16×32 3 (and first results from 16 × 48 3 ) lattices summarised in table 3. More details can be found in table 10 in appendix B. We consider three different values of the quark mass, corresponding to pion masses between 200 to 540 MeV. The scan corresponding to the largest pion mass at the critical point, denoted as scan B1 (in our naming convention the letter labels quark/pion masses while the number labels volumes), has been done at fixed hopping parameter, indicated by the subscript κ. Due to renormalisation and the change in the scale, the quark mass changes with the temperature in this scan. The scans at lighter quark masses, scans C1 and D1 with m π ≈ 300 and 220 MeV, are done along LCPs. To check the tuning of the quark masses we have measured the renormalised PCAC quark mass using the PCAC relation evaluated in the x 3 ≡ z-direction (see [104]). The simulation points in the (T, m ud ) parameter space are shown in figure 2. The plot shows that the tuning of the quark mass works well in the region below T C , while we see that we get smaller quark masses than expected above T C . It is unclear to us whether this is a cutoff effect or if our interpolation just becomes worse in this region (cf. appendix A.4). The quantities relevant for the transition temperature, i.e. the Polyakov loop and the chiral condensate, have been measured during the generation of the configurations with a separation of 4 MDUs. For the measurement of the condensate we have used 4 hits with a Z 2 × Z 2 volume source. The exception is the B1 κ scan, where the chiral condensate has only been measured on the stored configurations, using 100 hits. The screening masses have been measured on the stored configurations. For scan B1 κ , configurations have been saved every 200 MDUs, for scan C1 every 40 MDUs and for scan D1 every 20 MDUs. The results for the expectation values of Polyakov loops, the chiral condensates and screening masses are tabulated in tables 11 and 12 in appendix B.

The pseudocritical temperature
The first step of our investigation of the thermodynamics of QCD is the extraction of the pseudocritical temperatures. Since we are dealing with a crossover rather than a true phase transition there is no unique definition of the critical temperature and estimates for T C will depend on the defining observable.
To determine T C we will primarily look at the Polyakov loop and the chiral condensate. In particular, the inflection point of the Polyakov loop will define the deconfinement transition temperature T dc C , while the peak in the susceptibility of the chiral condensate defines the temperature of chiral symmetry restoration, which will be our main estimate for T C .

Polyakov loops
We start with the extraction of T dc C using the (unrenormalised) smeared Polyakov loop. We use Re(L) S since it typically shows stronger signals for the transition. We have, however, checked the agreement with the results for the unsmeared Polyakov loop explicitly (see also [59]). The results are shown in figure 3. The observable Re(L) S in scans C1 and D1 develops the S-shape characteristic for a phase transition, with some fluctuations in the vicinity of the inflection point. For scan B1 κ the S-shape is not as prominent, possibly due to the limited temperature range explored. To extract the inflection point we have fitted Re(L) S to the form of an arctangent. To check the model dependence of the results we have performed alternative fits using the Gaussian error function. Both fits tend to describe the data reasonably well and give similar χ 2 /dof values. The resulting curves from the arctangent fits are shown as black lines in figure 3. The results for the associated transition temperatures T dc C are given in table 4. Evidently the estimate for the uncertainty of the inflection point from the fit cannot be reliable due to the strong fluctuations in its vicinity. To account for this additional uncertainty we have assigned another systematic error of 10 MeV to the result from the fit, reflecting the size of the interval where we observe deviations from the smooth behaviour of the Polyakov loop. The shaded areas in figure 3 represent the estimates for the transition region. In the vicinity of T dc C the Polyakov loop susceptibility increases and shows fluctuations that can be interpreted as the onset of peak-like behaviour. Owing to this typical behaviour, we suspect that T dc  could be somewhat overestimated for scan D1.
taken into account. However, the position of the peak in the susceptibility should not be affected, since the additive renormalisation gives regular contributions around T C . Figure 4 displays the results for the disconnected susceptibilities for the unsubtracted and subtracted bare condensates. We define the pseudocritical temperature for chiral symmetry restoration through the position of the peak in the susceptibility of the subtracted condensate. To determine T C we fit the susceptibility to a Gaussian. Since the error estimate for T C from the fit will likely underestimate the true uncertainty we take the full spread of points included in the fit as a conservative error estimate. The resulting values for T C are given in table 4 and are shown as shaded areas in figure 4. The black curves correspond to the fit. Scan B1 κ is a problematic case, since, due to the change of the quark mass, the scan remains longer in the vicinity of the critical region, T C (m ud ). Consequently, the peak is broad and we obtain large uncertainties for both, T C and m C ud . Comparing the results for T dc C with T C , we see that they mostly agree within errors. The exception is D1, where we find that T dc C lies somewhat above T C . Like other studies in the literature we see that T C decreases with the quark mass [26][27][28][29][30][31][32][33][34][35][36], which is a general feature, persisting even when dynamical heavy quarks are included (e.g. [14,15]).
In figure 5 we show the results for the fully renormalised condensate and the fully renormalised subtracted version for scans C1 and D1. As expected, the condensates start close to zero and show a rapid decrease in the approach to T C . Around T C both condensates show fluctuations, especially the standard condensate fluctuates quite strongly. We note that T C , as defined by the inflection point of the renormalised condensate, does not necessarily have to agree with the peak of the susceptibility for a broad crossover.
We also have preliminary data for the condensates from simulations on 16×48 3 lattices. These were mainly used to check finite size effects on the screening masses below and are insufficient for a precise finite size scaling analysis. However, they are fully consistent with saturating susceptibilities and a crossover, as expected.

Scaling in the approach to the chiral limit
Following section 2.4.1 we will now try to extract information about the order of the transition in the chiral limit by looking at the scaling of the temperatures with the quark mass. With three transition temperatures at our disposal and their relatively large uncertainties we are not in the position to extract the critical exponents from a fit of the data for T C .  Instead, we fix the critical exponents to the ones from the different universality classes and check whether any particular scenario is favoured by our data. We start by fitting the results for T C from table 4 to eq. (2.17) with the critical exponents from the O(4) and U(2) scenarios. The results are shown in figure 6. The plot highlights the similarity of the two curves, indicating that the scaling of the transition temperatures alone will not be sufficient to distinguish between the two scenarios. We note that this is likely to remain true even when the error bars are reduced by an order of magnitude. Potentially, a scaling  analysis of the order parameter and its susceptibility (see [9]) might help in this respect. However, this would demand the knowledge of the scaling function from eq. (2.16) for the U(2) case. A distinction will only be possible if the scaling functions differ significantly.
In the first order scenario we also have to fix the value of the critical quark mass m cr ud . Since m cr ud is poorly constrained by the fit, we can try fits for different fixed values of m cr ud and look for minima in χ 2 /dof. We find that, not unexpectedly, χ 2 /dof is very flat and does not exhibit a minimum. Furthermore, χ 2 /dof is always of the same order as for the second order fits. As a typical case we show the curve obtained for m cr ud = 1.7 MeV in figure 7. As for the scaling with the O(4) and U(2) critical exponents, the Z(2) curve agrees very well with the data and is hardly distinct from the curves of figure 6.
Obviously, none of the scenarios is ruled out by the above analysis. These findings are in agreement with the results from the tmfT collaboration [36]. When we assume that the data shows consistency with one of the second order scenarios and extract the associated Both results are consistent with the findings from the tmfT collaboration for O(4) scaling, T C (0) = 152(26) MeV [36], and are on the lower side of the results for Wilson fermions at N t = 4, T C (0) = 171(4) MeV [33], and of the study of the QCDSF-DIK collaboration with different N t values T C (0) = 172 (7) MeV [35]. Calculations using staggered fermions only quote values for the critical coupling in the chiral limit, without providing results for the lattice spacing. The two results in eq. (3.1) indicate that the result for T C (0) is not sensitive to the universality class used for the extrapolation. This property is just another manifestation of the difficulty to distinguish between the two scenarios and shows that even a reduction of the error bars by an order of magnitude, in combination with results at much smaller quark masses, might not be sufficient using the scaling of the transition temperatures alone.

Screening masses and chiral symmetry restoration pattern
We now turn to the investigation of screening masses and the chiral symmetry restoration pattern. Since the behaviour of screening masses below and close to T C in general depends on the quark mass we will focus on the scans along LCPs in this section, i.e. on scans C1 and D1. The screening correlators have been measured in the x 3 direction on the stored configurations using unsmeared point sources. To make efficient use of the generated configurations we have computed the correlation functions for 48 randomly chosen source positions (see also [64]). Compared to Ref. [58], we have thus enlarged the statistics by another factor of three. Screening masses are extracted from the effective mass 2 , aM eff X , defined by the formula for the inverse hyperbolic cosine, where C = C XX . Since in scans C1 and D1 the spatial extent is rather small we could not find reasonable plateaus in most of the cases due to contaminations from excited states (see also [62]). To take the contaminations into account we have fitted the results for the effective mass to the form where A and ∆ are additional fit parameters. In figure 8 we show examples for the effective masses in the different channels above and below T C , together with the results from the associated fits. The formula for the effective mass follows to leading order when the contamination from the first excited state is included in the correlation function. In this case ∆ represents the energy gap between groundstate and first excited state in this channel. However, when contaminations from higher excited states become important ∆ (and A) will also contain contributions from those. In both cases a fit to the form (3.3) is known to improve the extraction of the groundstate energy significantly and to remove most of the contaminations of the excited states. Note that neither A nor ∆ are of direct importance for the following analysis. The fits to the form (3.3) typically work well for a variety of fitranges with starting points z min in the range from 4a to 9a for P and S channels and between 6a and 11a for V and A channels. 3 In these regions the results for M X do not depend significantly on the particular choice of fitrange and also the change of A and ∆ is only significant in a few cases. We show the results for M X , A and ∆ from two representative cases in P and S channels close to the critical temperature for D1 (where contaminations from excited states are typically expected to be most pronounced due to the small quark mass in the scan) in figure 9. Note, that the S channel is the more problematic one, both due to large statistical uncertainties and large contaminations from excited states (see also figures 8 and 12). In fact, the results shown in the right panel of figure 9 for the S channels are an example for the case where the fits do not work starting from z min /a = 8. Given these results, we conclude that for these values of z min the main contribution to ∆ comes from the first excited state. Within these regions where the fitparameters are insensitive to the particular choice for z min we can choose to work with any value of z min . We decided to use z min = 6a for P and S channels and z min = 8a for V and A channels. For these values all  Figure 10. Results for the screening masses in the different channels for scans C1 (left) and D1 (right). The screening masses are normalised to the asymptotic limit M ∞ .
of the fits give a reasonable χ 2 /dof around 1. 4 The results for the screening masses are shown in figure 10. At T /T C ≈ 0.7 the screening masses show the expected splitting from the zero temperature meson masses [88]. While the masses in the P and V channels initially remain constant, indicated by a slight decrease of M/T , the screening masses in S and A channels decrease drastically in the approach to T C . Around T /T C ≈ 0.9 all screening masses start to increase. In particular, the screening masses in P and S channels are drastically enhanced. Around T C the screening masses in the V and A channels are mostly degenerate and around 85 to 90% of the asymptotic 2πT limit, independent of the quark mass in the scan. This is consistent with the findings in simulations with staggered fermions [89,90].
The screening masses in P and S channels move closer together and fluctuate strongly around T C . For the P channel the screening masses are only around 40 to 50% of the 2πT limit for scan C1 with a pion mass of around 300 MeV and 35 to 40% for scan D1 with a pion mass around 200 MeV, indicating a quark mass dependence of the properties of pseudoscalar (and scalar) states around T C . The screening mass in the S channel is typically around 10% larger than M P . In the temperature interval covered by our calculations above T C , all screening masses are below the asymptotic high temperature limit. Note that weak-coupling calculations [93,94] predict the asymptotic approach to occur from above, implying that the screening masses must cross the value 2πT at a certain temperature.
Our findings are in good agreement with the results for screening masses obtained with staggered fermions [89,90]. In simulations within the quenched approximation, finite size effects have been found to be significant up to aspect ratios of N s /N t = 4 [105,106].
To get an idea about their magnitude, we compare the effective masses in the different  Figure 11. Results for the screening mass differences in the P and S channels, ∆M P S , and in V and A channels, ∆M V A , for scans C1 (left) and D1 (right). The differences are normalised to 2πT . channels at T = 150 MeV from scan C1 (with a volume of 32 3 ) with those obtained from a simulation with the same parameters but an increased spatial volume of 48 3 . The comparison is shown in figure 12, no significant finite size effects are visible in the data. Note that the comparison is done for the lattice with the smallest value of M P L (which is the relevant quantity governing the size of finite volume effects in the confined phase) in the scan C1, being smaller than M P L for all simulation points in the transition region. We thus conclude that our final results are not strongly affected by finite size effects. The chiral symmetry restoration pattern can be investigated by the degeneracies of the screening masses. To extract the differences, we have used the plateau in the effective masses of the ratios of the two correlation functions. In these ratios some of the fluctuations between different ensembles, evident in figure 10, and statistical fluctuations cancel. For these differences, some of the contaminations of the excited states cancel as well, so that we could use a fit to a constant, where we reduced the fitrange to the last few points. As for the fits to extract the screening masses, we have explicitly checked that our results do not depend significantly on the particular choice for the fitrange. This is in particular true when we consider that the main uncertainties in the following analysis are coming from the fluctuation between different ensembles in the region close to T C . On top of these checks for the fits, we also checked the extraction of the mass splittings using a fit to the plateau in the ratio of effective masses M X /M Y to obtain another estimate for the difference via In addition, we have also compared the results from the direct differences of screening masses. Note that in a very few cases the fit for ∆M P S failed, and in one case this also happened within the transition region (β = 5.37 for scan D1). We have excluded this datapoint from the following analysis, but we have checked with the mass difference of the independently determined screening masses that ∆M P S for this β value indeed lies within the uncertainty of the final estimate. In figure 11 we show the results for the mass differences. The plot indicates the degeneracy of M V and M A and the associated restoration of the SU A (2) symmetry around T C for both scans. As already noted above, this is in agreement with the results from staggered fermion simulations [89,90]. This also adds to the confidence concerning the extracted transition temperatures. In particular, both estimates for the mass difference agree very well over the whole temperature range for these channels (the ratio fit does not work for the lowest temperatures so that we have no reliable results in that region). For M P and M S there is still a significant mass splitting in the region around T C , which, however, appears to become weaker with decreasing quark mass. For the following analysis we will conventionally use the direct results for the mass differences, i.e. the results labeled by ∆M XY . The persistent mass splitting in P and S channels implies a residual breaking of U A (1) around T C in agreement with what has been found with staggered [89], domain wall [10,11,107] and overlap [38,39] fermion formulations. To make any statement about the fate and strength of the breaking in the chiral limit, however, we still need to perform  Table 5. Results for ∆M P S at T C for the different scans. The first error is statistical, the second reflects the systematic error owing to the fluctuations in the transition region. The value labeled with "m ud = 0 linear" is the result from the linear chiral extrapolation to all three points, as explained in the text, and the value labeled with "m ud = 0 sqrt" assumes a quark mass dependence proportional to √ m ud . The uncertainty on the result of the chiral extrapolations contains the statistical as well as the systematical uncertainty.
a chiral extrapolation, which is the topic of the next section.

On the relative size of the U A (1) breaking effects around T C
An important question is which amount of symmetry breaking is "strong" or "weak". When looking at the domain wall fermion results [10,11] for chiral susceptibilities, one might be led to the conclusion that the breaking is significant in the chiral limit. The same is true if one looks at the eigenvalues of the associated Dirac operator [107]. However, the chiral susceptibilities include contributions from contact terms which might give an additional contribution that overwhelms the effect of chiral symmetry breaking. The eigenvalues are a more sensitive probe of the U A (1) breaking. In this case, studies with the overlap operator [38][39][40] have indicated that the residual breaking in the chiral limit might be weak in contrast to the breaking implied by the eigenvalues of the domain wall operator at finite extent of the fifth dimension.
In contrast to the studies discussed above, we use non-perturbatively O(a)-improved Wilson fermions, which break chiral symmetry explicitly. The symmetry only becomes restored in the continuum limit. Since we use the mesonic screening spectrum as a probe, the results are expected to approach the continuum with O(a 2 ) corrections. On our relatively fine lattices (N t = 16), we expect these effects to be numerically small. We saw in the previous section that we do observe -to a good accuracy -the expected degeneracy between the vector and axial-vector screening above T C , signalling the restoration of the non-anomalous chiral symmetry. Similarly, if the U A (1) symmetry becomes effectively restored, we expect to obtain U A (1)-breaking mass splittings of O(a 2 ). In particular, since both the anomalous and the non-anomalous chiral symmetry are broken explicitly by the Wilson term, we expect the lattice artefacts to be of the same order of magnitude. Owing to the results for the difference between vector and axial-vector screening masses we expect this effect to contribute corrections of order 10 MeV. Furthermore, any accidental cancellation between lattice artefacts and a mass splitting in the continuum can only happen on this scale. Determining the strength of the breaking with a relative precision better than a few MeV requires taking the continuum limit.
We will now look at the chiral extrapolation of the mass difference ∆M P S in physical units. First of all, using the value m a0 = 980(20) MeV from the particle data group [108] we obtain as an estimate for the zero temperature reference value at the physical point This value might serve as an estimate for the effect in the mass difference when the breaking of U A (1) is substantial. However, this estimate is valid for physical, non-zero light quark mass, while we are interested in its value in the chiral limit. The pion mass vanishes in the chiral limit, so that ∆M T =0,m ud =0 . Next we need an estimate for m m ud =0 a0 , for which we use the following ansatz: we assume that the difference between chiral limit and physical point is solely due to the change of the constituent quark masses m const q in the meson and compare with another "iso-vector" scalar particle in the review of the Particle Data Group (PDG) [109], namely the K * 0 , where one of the u/d quarks is replaced by a strange quark. This results in the ansatz where C is a proportionality constant, which can be determined from eq. (3.7). Using C we can estimate the mass of the scalar in the chiral limit following where the factor 2 comes from the fact that we need to send the masses of two quarks to zero. Using the numbers from the PDG [109] we obtain the final estimate m m ud =0 a0 = 945(41) MeV. The error estimate follows from the uncertainties associated with the masses of the a0 and the K * 0 mesons. This is a rather crude estimate, but it is unlikely that it underestimates the effect by an order of magnitude (even then m m ud =0 a0 ≈ 600 MeV, which does not change the picture dramatically). Our final estimate for the chiral limit is (3.9) The width of the transition region must be taken into account when we extract an estimate for ∆M P S from our simulations. We thus compute the difference from a fit to a constant to the data points in the grey bands in figure 11. The spread of the results in the region is taken as a systematic uncertainty on top of the statistical uncertainty of the average. The results from this procedure are listed in table 5. Here we have also included a result for scan B1 κ to be able to perform a sensible chiral extrapolation. Unfortunately, B1 κ is not at fixed quark mass and thus remains longer in the vicinity of T C , since the latter increases with the quark mass. This accounts for the rather large error bars for the associated ∆M P S .
To perform the chiral extrapolation for ∆M P S we need to deduce its quark mass dependence. Since the pion is a Goldstone boson, its mass is expected to be proportional to √ m ud , at least at small temperatures, T < T C . On the other hand, the mass of the scalar should depend linearly on the quark mass which might also be the case for the pion at T C , where chiral perturbation theory breaks down. Given that we have only three data points at our disposal with relatively large uncertainties, our data clearly does not allow for a detailed investigation of the quark mass dependence of ∆M P S . We thus perform two types of fits; (i) linear in m ud , (ii) proportional to √ m ud . The results for the two different types of fits including all three data points are listed in table 5. We see that both results are consistent with zero within the relatively large error bars. As our final estimate we will thus use the linear fit. The associated result is shown in figure 13. We have also checked the robustness of the result with a linear fit using only the data from scans C1 and D1. For this chiral extrapolation the central value remains within the error bar of the result quoted in table 5, but the uncertainty of the result increases. Figure 13 indicates that the U A (1)-breaking screening-mass difference is a fairly small effect and has a mild quark-mass dependence only. In fact, at m ud = 0 the breaking effects are strongly suppressed compared to the effect in the same quantity at zero temperature and consistent with zero. Our result is thus in qualitative agreement with the results from the spectrum of overlap fermions [39,40].

Conclusions
In this paper we studied the finite temperature transition of two-flavour QCD on 16 × 32 3 (and 48 3 ) lattices. In particular, we have presented our results for the deconfinement and chiral symmetry restoration temperatures, extracted from the inflection point of the Polyakov loop and the peak in the susceptibility of the subtracted chiral condensate, respectively. In agreement with previous studies in the literature, we find both temperatures to decrease with the quark mass [26][27][28][29][30][31][32][33][34][35][36]. Our results for the chiral symmetry restoration temperatures, reported in table 4, are within uncertainties consistent with the values quoted by the tmfT collaboration [36,110].
The ultimate goal of our ongoing efforts is to determine the order of the transition in the chiral limit. In an attempt to extract information about the chiral limit, we tested for scaling of the transition temperatures in the approach to the chiral limit. As already reported in [58], the scaling behaviour alone is not sufficiently constrained to distinguish between a second order O(4) chiral transition or a first order transition with a Z(2) endpoint. This is consistent with earlier findings of the tmfT collaboration [36], even though we were able to reduce the pion masses in our study down to 200 MeV. Extrapolations of the critical temperatures to the chiral limit are not very sensitive to the universality class. Our results, when interpreted in terms of the O(4) and U(2) scenarios are consistent with the findings from the tmfT collaboration for O(4) scaling [36], and somewhat smaller than those for Wilson fermions at N t = 4 [33] and the QCDSF-DIK collaboration with different N t values [35]. The fact that the chiral transition temperature does not appear to be very sensitive to the universality class used for the chiral extrapolation is one particular manifestation of the difficulty to distinguish between different scenarios. Even a reduction of the error bars by an order of magnitude in combination with results from smaller quark masses might not be sufficient to allow to distinguish between the universality classes when using the scaling of the transition temperatures alone.
As an alternative, we have investigated the strength of the anomalous breaking of the U A (1) symmetry in the chiral limit by computing the symmetry restoration pattern of screening masses in various isovector channels. At T /T C ≈ 0.7 the screening masses assume values close to the zero temperature meson masses. Initially the masses in the scalar and axial-vector channels decrease, before at T /T C ≈ 0.9 all masses start to increase. Around T C the screening masses in the vector and axial-vector channels are degenerate and about 85 to 90% of the asymptotic 2πT limit, while the screening masses in the pseudoscalar channel are about 35 to 50% of this limit and exhibit a significant dependence on the quark mass. The screening mass in the scalar channel is typically around 10% larger. When going to higher temperatures all screening masses approach the 2πT limit from below, while leadingorder weak-coupling calculations [93,94] predict an asymptotic approach from above. All these findings are in qualitative agreement with the results found in N f = 2 + 1 simulations with staggered fermions [89,90].
To quantify the strength of the U A (1)-anomaly in the chiral limit, we use the difference between scalar and pseudoscalar screening masses, ∆M P S . Unlike susceptibilities, screening masses probe exclusively the long distance properties of the correlation functions and thus do not suffer from contact terms. The comparison of the chirally extrapolated value, ∆M m ud =0 P S = −81(282) MeV to its zero temperature analogue ∆M T =0,m ud =0 P S = −945(41) MeV (cf. eq. (3.9)) suggests that the U A (1)-breaking is strongly reduced at the transition temperature (see also figure 13). If this effect persists in the continuum limit, it disfavours a chiral transition in the O(4) universality class. through the John von Neumann Institute for Computing (NIC) within project HMZ21. The correlation functions and part of the configurations were computed on the dedicated QCD platforms "Wilson" at the Institute for Nuclear Physics, University of Mainz, and "Clover" at the Helmholtz-Institut Mainz. We also acknowledge computer time on the FUCHS cluster at the Centre for Scientific Computing of the University of Frankfurt. This work was supported by the Center for Computational Sciences in Mainz as part of the Rhineland-Palatine Research Initiative and by DFG grant ME 3622/2-1 Static and dynamic properties of QCD at finite temperature. B.B. has also received funding by the DFG via SFB/TRR 55 and the Emmy Noether Programme EN 1064/2-1.

A.1 Simulation algorithms and associated constraints
The simulations have been done using the deflation accelerated versions of the Schwarz [67,68] (DD) and mass [69,70] (MP) preconditioned HMC algorithms. Both algorithms employ the Schwarz preconditioned and deflation accelerated generalised conjugate residual (DFL-SAP-GCR) solver introduced by Lüscher in [71,72]. Due to the efficient solver, both algorithms exhibit an improved scaling behaviour when lowering the quark mass and the lattice spacing.
The block structures used in the solver and the preconditioning of the HMC algorithm impose constraints on the size of the local lattices allocated to the single processors. The main limitation follows from the constraint that Schwarz preconditioning blocks in the solver (SAP-blocks) have a minimal size of 4 4 [71]. At least two of these blocks need to fit into a sublattice, so that the minimal sublattice size when using a SAP solver is 8 × 4 3 . This restriction translates directly to the version of the MP-HMC algorithm from [70]. The restrictions for the sizes of the DD-blocks are essentially the same as for the SAP-blocks, but it also affects the efficiency of the HMC preconditioning, since the separation of modes works most efficiently if the DD-blocks have a minimal size of about half a fermi in each direction [67]. The inverse critical temperature corresponds to a temporal extent of about 1 fm, meaning that the size of the DD-blocks should be half the size of the temporal extent. For our scans the optimal sublattice size for the DD-HMC is thus given by 16 × 8 3 . Further constraints follow from the use of even/odd preconditioning in the different levels of the solver (see [71,72] for the details). Summarizing, the algorithm restricts one to use lattice sizes of N t/s = 8, 12, 16, 20, . . ..
In our initial runs (see [59][60][61]) we used the DD-HMC algorithm. A drawback of the algorithm is that a sizeable fraction R of the links remain fixed during the molecular dynamics evolution [67]. While the ergodicity can be restored by shifting the gauge field between two HMC trajectories, the autocorrelation between two trajectories increases significantly. For typical blocksizes of 4 4 and 8 4 one obtains R = 0.09 and 0.37. For N t = 16 the optimal block size is 8 4 , for which the autocorrelations are expected to be enhanced by about a factor of 3 compared to the MP-HMC algorithm. In particular on large scale machines, one is often forced to use a large number of processors and decrease the local β 5. 20 5.30 5.50 r 0 /a 6.15( 6) 7.26 ( 7) 10.00(11) κ c 0.136055(4) 0.136457(4) 0.1367749(8) Table 6. Input from zero temperature for scale setting and the estimation of LCPs taken from [74].
system size, meaning that autocorrelations increase drastically and the algorithm becomes inefficient. An example is the B3 κ lattice discussed in [58]. Due to the suppression of low modes above T C accompanying the effective restoration of the chiral symmetry, the main motivation for using deflation gradually disappears. The reduced benefit of deflation can be observed in practice. However, for the lattices and quark masses in this study, there is still a significant acceleration even above T C . In the region around and below T C we had rare appearances of problems with deflation, most likely stemming from the fact that the "little" Dirac operator (see [72] for the details) has large condition numbers. Mostly these issues could be solved by changing the parameters of the deflation subspace.

A.2 Error analysis
For the error analysis we employ the bootstrap procedure [111] with 1000 bins. The data for different time slices of correlation functions is correlated, which means that in fits the full covariance matrix has to be taken into account. In practice, however, the uncertainties for the entries of the correlation matrix are often not precise enough for a stable leastsquare minimisation (see [112] for instance). This is in particular true for screening masses that typically have a bad signal-to-noise ratio. All our fits to correlation functions have thus neglected the non-diagonal terms of the correlation matrices.
For scale setting and renormalisation we have used interpolations of quantities that have their own uncertainties and for which we have no access to the raw data. To take these uncertainties into account, we have generated uncorrelated pseudo-bootstrap distributions with 1000 bins, whose width reproduces the quoted uncertainties. The interpolations have then been done for each bin, giving bins for the desired quantities at each value of the coupling.

A.3 Interpolation of zero temperature quantities
In this appendix we discuss the interpolations of several quantities, such as lattice spacing and quark masses. For scale setting we perform an interpolation of r 0 /a obtained from the CLS lattices [74,113], summarised in table 6. The final result is obtained using the ansatz [114] ln The results for the coefficients are tabulated in table 7 and the interpolation is visualised in figure 14 (left). The plot displays the good agreement between the two interpolations, indicating that in the region 5.20 ≤ β ≤ 5.50, systematic errors due to the ansatz for r 0 /a(β) are negligible. Another important quantity is the critical hopping parameter κ c . We interpolate the result from [74], listed in table 6, using the two different ansätze, The results are also listed in table 7 and the interpolations are shown in figure 14 (right).
There is a slight model dependence in the region 5.35 β 5.45. However, the interpolation mostly concerns the estimation of LCPs and does not enter the final analysis.
For the other renormalisation factors and improvement coefficients we have used known interpolation formulas from the literature: for Z V we have used the interpolation formula from [115]; for c A the one from [116]; for Z A and Z P we have used the results from [74]; for b m and [b A − b P ] we have used the non-perturbative result from [117].

A.4 Estimating lines of constant physics
LCPs can be realised once we have an analytic relation between the bare parameter κ and the renormalised quark mass m ud . This relation can be obtained by using the two different definitions for m ud in the case of Wilson fermions. The first definition uses the bare quark massm [118], where Z m is the scheme dependent mass renormalisation factor and b m is an O(a) improvement coefficient. Another possibility is to use the PCAC quark mass [119], (17) 3.9(62) -0.21 (58) 184(120) -79 (45) (37) -20 (14) 1.9(13) 0.2 Table 7. Results of the fits for the interpolations of zero temperature results and the analytic relation am PCAC (β, κ). Note, that the for the interpolations there are as many parameters as data points, so that it is a parameterisation rather than a fit.
Here Z A and Z P are the renormalisation factors for the axial current and the pseudoscalar density (Z P is scheme dependent) and b A and b P are due to improvement. m PCAC is the bare quark mass defined by the PCAC relation (see [74,77]). Both quantities in eqs. (A.5) and (A.6) are estimates for the renormalized quark mass m ud and differ only in lattice artifacts. From eqs. (A.5) and (A.6) we can thus infer (see also [117]) Eq. (A.7) is valid up to O(a 2 ) corrections and Z PCAC (β) only depends on the regularisation scheme (i.e. the lattice discretisation), but not on the renormalisation scheme. Inserting the relation between am and κ we obtain For κ c (β) we take the known relation from appendix A.3. When m ud is small am is small as well. Since the factor Z 2 contains the combination of b-factors from eq. (A.7), which is dominated by b m , and is a number smaller than one, we can expect that the second term in eq. (A.8) is negligible for small quark masses. We thus have to obtain only the β-dependence of the factor Z PCAC . To this end we use the data for m PCAC at T = 0 from [74,77], listed in table 8. Since data is available only for three couplings, we will use a simple second order polynomial, The fit with this ansatz works reasonably well, giving a small χ 2 /dof of 0.2. In fact, χ 2 /dof does not make any statements about the goodness of the ansatz for Z PCAC (β), but it shows that for the given range of quark masses am PCAC and am are indeed linearly related. The results for the coefficients z i are listed in  [62,77]. Note, that the final error bars on m ud obtained from the relation m ud (β, κ) are much smaller than the uncertainties of the coefficients suggest, due to compensations of changes in one parameter by the others. Figure 15 shows a comparison between the predictions from m ud (β, κ) with the results obtained from the actual simulations. The plot displays the good agreement for temperatures below and up to T C . Note, that the tuning for scan C1 has been done using an early version of the matching with less information from the T = 0 side. This explains the fact that the black points do not lie on a constant m ud line in that case. The good agreement between the measured quark masses and the predictions from the matching reported here indicates that this updated form of the matching works well. Above T C , m ud is typically lower than expected. It is unclear to us whether this result is an effect due to markedly different cutoff effects above T C , or if the interpolation becomes worse at the higher β values. Note that our zero temperature ensembles are located at the β-values corresponding to 150, 180 and 245 MeV, so that the systematic error associated with the interpolation is expected to be most severe around 210 MeV.

A.5 Interpolation of the zero temperature chiral condensate
For the computation of the renormalised condensate the computation of the T = 0 subtraction terms is mandatory. Here we have computed ψ ψ bare and P P on the CLS lattices, summarised with the measurement details in table 8. For the computation we have used 12 Z 2 × Z 2 volume (for ψ ψ bare ) and wall (for P P ) sources on each configuration.
We have tested several different functional forms for the interpolation and found simple polynomials in am to work best. Note, that P P is expected to diverge linearly in am, since the condensate, which is finite at am = 0, is proportional tom P P , which can be described by a polynomial. We have thus used  Table 8. Set of CLS T = 0 lattices used for the determination of the relation m ud (β, κ) and the computation of the condensate at zero temperature and the associated results. For the measurements of am PCAC and more details on the lattices we refer to [74,77]. (similarly forp i with q →q) and as usual 2am = 1/κ − 1/κ c (β) with κ c (β) obtained from the standard interpolation (A.3). Altogether, each fit has 9 fit parameters. The results are given in table 9 and the interpolations for the different β-values are shown in figure 16. The χ 2 /dof values are 3.3 and 5.5 for ψ ψ bare and P P , respectively. The tiny error bars on ψ ψ bare and P P explain the bad values for χ 2 /dof, despite the fit providing a satisfactory description of the data. Once more, the uncertainties of the fit parameters do not reflect the uncertainty on the results for the condensate. Note that we have neglected finite volume effects for the condensates, which appears to be justified for m π L 4.    Table 10. Run parameters of the simulations. We list the bare lattice coupling β, the hopping parameter κ, the number of molecular dynamics units (MDU), temperature, renormalised quark mass in the MS-scheme at a renormalisation scale of 2 GeV and the autocorrelation time of the plaquette (τ Up ).  Table 11. Simulation results for the (smeared) real part of the Polyakov loop, the bare chiral condensate and its disconnected susceptibility ( ψ ψ bare andχ bare ψ ψ ), the associated renormalised condensate ( ψ ψ ren ) in units of r 0 , and the subtracted versions ( ψ ψ