Symmetry restoration at high-temperature in two-color and two-flavor lattice gauge theories

We consider the $SU(2)$ gauge theory with $N_f=2$ flavors of Dirac fundamental fermions. We study the high-temperature behavior of the spectra of mesons, discretizing the theory on anisotropic lattices, and measuring the two-point correlation functions in the temporal direction as well as screening masses in various channels. We identify the (pseudo-)critical temperature as the temperature at which the susceptibility associated with the Polyakov loop has a maximum. At high temperature both the spin-1 and spin-0 sectors of the light meson spectra exhibit enhanced symmetry properties, indicating the restoration of both the global $SU(4)$ and the axial $U(1)_A$ symmetries of the model.


Introduction
We consider the SU (2) gauge theory with N f = 2 flavors of Dirac fundamental fermions, and study the finite-temperature behavior by using numerical methods based on formulating the theory on anisotropic lattices. The main purpose of this work is to collect evidence that the global symmetries of the model are implemented à la Wigner at high-temperature, where the condensate breaking global symmetry is expected to melt and the global symmetries to be linearly realized.
This model has been considered before in three different contexts, as it represents the prototype of non-trivial gauge theory in which lattice numerical methods have concrete potential to provide useful information about the dynamics of the underlying theory. First of all, it is a useful toy model for the study of generalizations of Quantum Chromo-Dynamics (QCD) at finite temperature T and finite chemical potential µ. One trivial reason for this is that the number of fundamental degrees of freedom is smaller than for two-flavor QCD, making the numerical treatment easier. Most importantly though, the fundamental representation of SU (2) is pseudo-real, and hence there is no sign problem. It is then possible to study the phase diagram of the model in the (T, µ)-plane, and to apply numerical techniques to extract its detailed structure. For an incomplete list of useful references on the subject see [1].
A second context in which this model is important is that of traditional technicolor (TC) [2,3]. The choice of SU (2) with 2 fundamental Dirac fermions yields the minimal model such that one can embed the electro-weak SU (2) L × U (1) Y group of the Standard Model of particle physics (SM) within the global symmetries of the matter field content. One expects spontaneous symmetry breaking to arise dynamically at the scale Λ, hence providing a natural way to implement the Higgs mechanism for giving mass to the electroweak bosons within a fundamental theory. Aside from the fact that, once more, the small number of degrees of freedom makes practical applications amenable to numerical treatment, the fact that the field content is minimal also minimizes the potentially problematic contributions to precision parameters such as the oblique S and T as defined by Peskin and Takeuchi [4], that on the basis of perturbative arguments one expects to grow with N f and N c , and that are not dynamically suppressed when one identifies Λ with the electroweak scale v W ∼ 246 GeV. The dynamics preserves a custodial SU (2) that further suppresses the T parameter, as the underlying masses of the fermions vanish.
The model has received some attention in a third context [5,6], as a concrete realization of the idea of Higgs compositeness [7]. This is a quite distinct framework in respect to traditional TC. The underlying dynamics is the same, being based upon a gauge theory with a given global symmetry, for which one expects the formation of a non-trivial symmetrybreaking condensate. Yet, one chooses to embed the electroweak gauge group into the global symmetry group of the theory in such a way that the fermion condensate does not break it. 1 The long-distance behavior of the theory is hence captured by an Effective Field Theory (EFT) that includes the SM gauge theory, supplemented by a set of light, composite pseudo-Goldstone bosons arising at the scale Λ, a subset of which is interpreted as the Higgs doublet field.
The gauging of the SM group explicitly breaks the global symmetries, and hence provides a potential for the Higgs fields. Additional ingredients, not arising from the SU (2) fundamental gauge theory, are invoked in order to drive spontaneous symmetry breaking in the Higgs sector, which ultimately yields electro-weak symmetry breaking (EWSB) at the scale v W Λ. For example, one has to introduce a mechanism to give mass to the SM fermions, which requires coupling the Higgs field to the quarks and leptons. It is well known that, as a byproduct of doing so, the theory yields radiative corrections to the Higgs potential due to loops of the top quark, naive estimates of which show that they can destabilize the minimum of the Higgs potential. In the following we will not discuss any of these points, related to realistic model-building in the electro-weak sector.
The reason why composite scenarios are viable within this model originates from the pseudo-real nature of the fundamental representation of SU (2). In particular, in the presence of two Dirac fermions, the global symmetry of the Lagrangian is enhanced from the U (1) A × U (1) t B × SU (2) t L × SU (2) t R global symmetry of QCD and TC to a U (1) A × SU (4) global symmetry, and the condensate breaks it to the Sp(4) subgroup. Excluding for the time being the anomalous U (1) A from the discussion, this yields 5 (pseudo-)Goldstone bosons, that form a multiplet of the unbroken Sp(4) ∼ SO (5). The gauging of SU (2) L × U (1) Y ⊂ SO(4) ⊂ Sp(4) splits the 5 into a 4 of SO (4), which is identified with the Higgs doublet, and an additional singlet, that may have important phenomenological implications.
In this paper, we compute the masses of the composite (meson) states created and annihilated by operators of the formQΓQ, with Γ = 1, γ 5 · · · , and discuss their dependence on temperature T . 2 In particular we track how the mass-splittings between parity partners change by going to high-temperature. In order to do so, we formulate the theory on anisotropic lattices, and use Monte Carlo methods to extract the spectral masses as a function of T . We are looking for clear signals of the restoration in the thermal bath of the much larger global symmetry of the underlying theory. This is the first step of a more ambitious and long-term program, which we envision will include also the study of the effects due to the presence of explicit symmetry-breaking terms, in particular due to the chemical potential µ, and to the weakly-coupled gauging of the SM electroweak group.
The paper is organized as follows. In Section 2 we describe the model and summarize effective field theory and symmetry arguments that play a role in the rest of the paper. In Section 3 we describe the lattice set-up used in the numerical calculations, particularly by explaining in details how the bare parameters are tuned in the presence of anisotropic lattices. In Section 4 we report our results, which we critically discuss in Section 5. Appendix A contains some useful notation about spinors, and we show explicitly how the enhanced global symmetry emerges. In Appendix B we summarize the algebraic properties of SU (4) and Sp(4), by providing an explicit example of generators for SU (4). Examples of the renormalized versus bare parameters are given in Appendix C.

The model: symmetry considerations
The matter field content consists of two (massive) Dirac fermions Q i a , where a = 1, 2 is the SU (2) color index and i = 1, 2 the flavor index. The covariant derivative is with V A µ the gauge fields, g the coupling, and T A the generators of SU (2) obeying Tr where the summations over flavor index i = 1, 2 and color index a = 1, 2 are understood, and where the field-strength tensors are defined in terms of the gauge bosons as We collect in Appendix A and B several useful relations between 2-component spinors q and 4-component spinors Q, as well as details about the algebra of SU (4) and Sp(4) (see also [5,6]). The global symmetry acting on the matter fields is U (1) A × SU (4), and we explicitly list the transformation properties of the fields in Table 1. It is convenient to define: and to write the mass explicitly as a matrix M ≡ m Ω, with Ω the symplectic matrix in Eq. (B.1). The index n, m = 1, · · · , 4 andC = −iτ 2 acts on spinor indexes. In the lower half of Table 1 we list the transformation properties of the composite field Σ 0 , as well as the (symmetry-breaking) spurion M . In the body of the paper, we will describe the finite-temperature properties of composite states that we identify with the pions π, ρ vector, a 1 axial-vector, and a 0 scalar mesons. In the rest of this Section, we summarize the basic properties of these objects, using the language of effective field theory (EFT). What results is a Lagrangian density that includes potentially heavy and strongly-coupled degrees of freedom, and hence does not yield a calculable weakly-coupled low-energy EFT in the usual sense. We use this language to guide our book-keeping exercise, focused on classifying the physical particles, their quantum numbers, and the degeneracies-in particular the difference of mass between the ρ and a 1 vectors and between the π and a 0 scalars-that are consequences only of the symmetry structure of the theory and its vacuum.

Composite states: scalars
In the low-energy EFT description, the real antisymmetric field Σ transforms as under the action of an element U of SU (4). The VEV Σ ∝ Ω breaks SU (4) to the Sp(4) subgroup. The generators T A with A = 1, · · · , 5 are broken, while T A with A = 6, · · · , 15 are unbroken. For instance, see Eq. (B.4) in Appendix B.
In terms of the matrix-valued π(x) = 5 A=1 π A (x)T A , the convenient parameterization automatically satisfies the non-linear constraint Σ † Σ = I 4 . The leading-order term of the low-energy EFT is (2.7) The pion fields are canonically normalized and hence f = f π is the pion decay constant. The quark mass is incorporated in the EFT by adding the symmetry-breaking term The expansion in pion fields confirms that the 5 pions are still degenerate, if not massless, in the presence of the explicit breaking given by the Dirac mass for the fermions, with The degeneracy of the five pions is a consequence of the unbroken Sp(4) ∼ SO(5) symmetry.
The spurion M formally transforms as M → U * M U † , so that if it were promoted to a field then L m would be manifestly invariant under the full SU (4) symmetry.
Here we pause to make two general observations. In the context of composite-Higgs models, the presence of a (small) mass term for the quarks is allowed, contrary to the TC case. While in the latter the quark mass explicitly breaks the gauge symmetries, in the composite-Higgs case the SM gauge group is a subgroup of Sp(4), and hence the term in Eq. (2.8) does not break it. The distinction between TC and composite-Higgs cases reduces (in the massless case) to a vacuum alignment issue driven by the weak gauging of the SU (2) L × U (1) Y symmetry. In the presence of a mass of the form in Eq. (2.8), this problem has a trivial solution: the mass m stabilizes the composite-Higgs vacuum. Yet, some caution is in order: if m π is large, it might become impossible to induce electro-weak symmetry breaking. We leave these and similar issue out of this study (see [12]), as in our numerical work all calculations are done with the SU (2) theory in isolation.
To describe the regime in which the symmetry is restored, which is expected to be realized at high temperature, we remove the non-linear constraint, and hence replace Σ by the field H ∼ 6, that transforms as a complete antisymmetric representation of SU (4). The kinetic term is where the normalizations are chosen so that all the fields have canonical kinetic terms. Unconstrained by symmetry considerations, the scalar σ (singlet of Sp (4)) is expected to have a large mass m σ , and in general decay fast to pions. Besides the SU (4) → Sp(4) breaking, the vacuum also induces the breaking of the (anomalous) U (1) A . To discuss it, we need to promote H a complex field, hence doubling the field content. We defineH with H a second real antisymmetric representation of SU (4). The action of U (1) A is where θ is the parameter of the U (1) A transformation. The field H introduces an additional Sp(4) singlet that is the analog of the η in QCD and 5 additional scalars that form a multiplet of the SO(5) unbroken symmetry, and are the analogue of the a 0 isovectors of QCD. The treatment presented here is indeed a generalization of what done in the context of the linear-sigma-model description of low-energy QCD [13]. The presence of the anomaly produces a large mass for η . At high temperatures both the fermion condensate and the effect of the anomaly are suppressed. Hence, the mass splitting between a 0 and π provides a measure of the level of breaking of U (1) A in addition to global SU (4), and can be used to look for SU (4) × U (1) A thermal restoration. Similar arguments hold in the case of QCD (see for example [14,15] and references therein).
Because the σ and η are flavor singlets, and the flavor-singlet sector of the spectrum is more difficult to study numerically than the flavored channels, we will study the a 0 -π mass splitting in order to discuss the restoration of the axial U (1) A at high temperatures. We will do so in the body of the paper, using numerical techniques based on the formulation of the theory on anisotropic lattices.

Composite states: vectors
The full set of spin-1 vector and axial-vector mesons spans the adjoint representation of the SU (4) global symmetry. A cartoon representing the EFT description of their longdistance dynamics is depicted in Fig. 1, and represents a generalization of hidden local symmetry [16][17][18]. One extends the symmetry from SU (4) to SU (4) A × SU (4) B , with SU (4) A weakly gauged, with coupling g ρ . Then one enlarges the field content to include two non-linear sigma-model fields S and Σ. The non-linear sigma-model S transforms as the bifundamental of SU (4) B × SU (4) A , while the field Σ transforms on the antisymmetric of SU (4) A : (2.14) In a composite-Higgs model, the SM gauge group SU The gauging of the SU (4) A symmetry means that (for global SU (4) B ) one has to introduce the covariant derivatives The quark masses also contribute to the masses of the spin-1 states in a more complicated way, that will be discussed elsewhere [19].
In the absence of the antisymmetric condensate (for Σ = 0), ρ and a 1 mesons would be exactly degenerate. Their mass splitting is hence a measure of the amount of breaking SU (4) → Sp(4). In the main body of the paper we use the mass splitting between ρ (vector) and a 1 (axial-vector) as a way to test whether the global symmetry is restored at high temperatures. The generalization to the case in which Σ is replaced byH does not require any new ingredients. In particular the restoration of the axial U (1) A and of the Table 2. Interpolating operators, and corresponding flavored particles (i.e. i = j in the interpolating operators), studied in the body of the paper. Color and spinor indexes (summed over) are understood.
global SU (4) can, at least in principle, be treated independently. We summarize in Table 2 the properties of the states discussed in the body of the paper. One of the purposes of this paper is to make the first steps towards a quantitative assessment of the relation between the two phenomena at high temperature, in the specific theory of interest here.

Lattice action
In this Section, we describe the discretized Euclidean lattice action used for our numerical study. For the gauge sector, we modify the standard plaquette action by treating the operators containing temporal gauge links separately from those solely containing spatial links, where β = 2N/g 2 and ξ 0 g are the lattice bare gauge coupling and the bare gauge anisotropy, respectively. The plaquette P is defined by where U µ (x) denotes the link variables. For the fermion sector, we use the Wilson action for fermions in the fundamental represention with the massive Wilson-Dirac operator given by where ∇ and ∇ * denote the forward and backward covariant derivatives, respectively: The ratio v µ of the bare fermion to gauge anisotropy is introduced as it can be different to unity. From the redefinition of the fermion field (Q → √ v t Q and m 0 → m 0 /v t ), along with the introduction of the fermion anisotropy ξ 0 For the rest of this paper we do not explicitly show the lattice spacings for convenience, i.e. a t = 1, except when we need to distinguish the spatial and temporal lattice spacings and to discuss the finite temperature. The bare anisotropy parameters, ξ 0 g and ξ 0 f , are renormalized such that physical probes at scales well below the cut-off ∼ 1/a exhibit Euclidean symmetry, i.e. ξ g = ξ f = ξ. For the input quark mass, M q , we parameterize the renormalized parameters (ξ g , ξ f , M q ) as functions of bare parameters (ξ 0 g , ξ 0 f , m 0 ). For a small region in the parameter space, we assume that the renormalized parameters are linear in the bare parameters. We further assume that we are in the region of light quark masses, i.e. M 2 P S ∼ M q , and arrive at the form [20] For each set of bare parameters, nonperturbative determinations of ξ g and ξ f are carried out through the interquark potential and the relativistic meson dispersion relation, respectively, which will be discussed in details in the following subsections.

Simulation details
We consider the lattice action in Eq. (3.1) and Eq. (3.3) with two mass-degenerate Wilson fermions. Configurations are generated using the Hybrid Monte Carlo(HMC) algorithms with the second order Omelyan integrator for Molecular Dynamics(MD) evolution, where different lengths of MD time steps δτ µ are used for gauge and fermion actions such that the acceptance rate is in the range of 75 − 85%. The simulation codes are developed from the HiRep code [21] modified by implementing the gauge and fermion anisotropies described  (7) 0.1382(11) 6.13(6) 6.42(13) · Table 3. Simulation parameters and results for the tuning of the lattice bare parameters of an anisotropic lattice. The masses of pseudoscalar (PS) and vector(V) mesons are measured in units of a t .
in Section 3.1. To optimise the acceptance rate, we also treat the variance of temporal and spatial conjugate momenta differently by introducing a new tunable parameter [22], which is essentially equivalent to the multiscale anisotropic molecular dynamics update [20]. Without changing the validity of the algorithm, such a setup is helpful for the anisotropic lattice calculations through balancing the temporal and spatial MD forces: typically the former is larger than the latter approximately by the anisotropy in the lattice spacings.
Except the lattice of N t × N 3 s = 128 × 10 3 for the investigation of finite volume effects, all of the numerical calculations for the tuning of bare parameters are performed on N t × N 3 s = 128 × 12 3 lattices. We use periodic boundary conditions in each direction of both link variables and fermion fields. 5 Twelve ensembles are created with different bare quark masses, gauge and fermion anisotropies at β = 2.0, where the details are found in Table 3. Thermalization and autocorrelation times are estimated by monitoring the average plaquette expectation values. For each ensemble N conf = 138 − 300 configurations are accumulated after 200 trajectories for thermalization, where every two adjacent configurations are separated by one auto-correlation length of which the typical size is 8 ∼ 12 trajectories. The statistical errors for all quantities extracted in this work are obtained using the standard bootstrapping technique.

Gauge anisotropy
The gauge anisotropy ξ g is determined from the static potential using Klassen's method [23]. We first define the ratios of spatial-spatial and spatial-temporal Wilson loops by respectively. In an asymptotic region, these ratios fall exponentially with the linear interquark potential and do not depend on r, R s (r, y) ∼ e −asVs(yas) and R t (r, t) ∼ e −asVs(tat) . Finite volume effects are expected to be suppressed since they are canceled out in the ratios [23,24]. As the interquark potential at the same physical distance should yield the same value, one can extract the anisotropy ξ g by imposing R s (r, y) ≡ R t (r, t = ξ g y).
In practice, we determine ξ g by minimizing [24] L(ξ g ) = r,y (ξ g ; r, y), where ∆R s and ∆R t are the statistical errors of R s and R t , respectively. In the original Klassen's approach, the planar Wilson loops are considered where r is either x or z. A typical difficulty in this approach is the limited number of data points as one quickly encounters a severe signal-to-noise problem in the calculations of the large Wilson loops. By noting that r can be any two-dimensional path in the x-z plane with r = √ x 2 + z 2 , we extend the Klassen's method by including nonplanar Wilson loops along the closed path C y (x, z, y) and C t (x, z, t) with x ≥ z. To maximize the overlap with the physical ground state the shortest paths in the x-z plane are considered using the Bresenham algorithm which has been applied for the lattice study of quark antiquark potential, i.e. see [25]. Analogous to the planar case, we define r = (x, z) and r + 1 = (x + 1, z) for a fixed value of z.
Using the generalized Klassen's method, we are able to secure enough data points having reasonable statistical errors. As a consequence, not only do we find the clean signal of an asymptotic region in which ξ g converges, but also reduce the statistical error of the gauge anisotropy ξ g . However, due to the breaking of rotational symmetry on the lattice, results obtained mixing on-axis and off-axis loops might be affected by a large systematics. To investigate this issue, we calculate ξ g by minimizing the function (ξ g ; r, y) with z = 0, · · · , 3, corresponding to the different shapes of the 2-dimensional paths. The results are shown in Fig. 2. We find no significant deviations between colored data, suggesting that any potential effect of the breaking of rotational symmetry cancels in the ratios of Wilson loops in Eq. (3.8). For r * y ≥ 5 all data points are statistically consistent with one another. The measured value of ξ g is denoted by the blue band in the figure, where its extraction is discussed in the following.
In Fig. 3 we plot ξ g , obtained by using Eq.   y = 1 (blue circle), planar and nonplanar Wilson loops except y = 1 (red square), and planar and nonplanar Wilson loops except y = 1 and r = 1 (green diamond). The largest value of r * y is the one before we encounter significant numerical noise.
For all ensembles we find that ξ g converges to the asymptotic value at around min(r * y) = 4 ∼ 6 and thus we choose min(r * y) = 6, as for this value we expect the size of systematic errors to be small compared to the statistical error. Since the inclusion of y = 1 Wilson loops causes significant systematic effects due to short-range lattice artefacts, as can be seen in the plots (see also the discussion in [20,24], in the case of QCD), we exclude these Wilson loops for the determination of ξ g . In summary, we calculate the asymptotic value of ξ g using planar and nonplanar Wilson loops, except the ones having y = 1, at min(r * y) = 6 and the results are reported in Table 3.

Fermion anisotropy
The fermion anisotropy ξ f is determined through the leading-order relativistic dispersion relation of mesons where N s is the spatial lattice size. The energy E and the mass m are in units of a t , while the momentum p is in units of a s . In the Euclidean formulation, meson two-point correlation functions exponentially fall off with the lowest energy at an asymptotically large time. In practice, it is useful to define an effective mass, where C(t) is the ensemble average of meson correlators. Then, ground state energies are obtained from a constant fit to the plateau of m eff in the asymptotic region of large t. In the case of zero momentum these energies are nothing but the meson masses. The measured masses of pseudoscalar and vector mesons are reported in Table 3.
As an example, in Fig. 4 we show the effective mass plots for pseudoscalar and vector mesons with m 0 = 0.2, ξ 0 g = ξ 0 f = 4.7 and β = 2.0. We construct the meson interpolating operators at source and sink using point sources. Various momentum projections with | n| = 0, 1, 2, 3 are denoted by red, green, yellow and brown colors, respectively, while the measured ground state energies are denoted by the blue bands.
In Fig. 5 we plot the resulting squared energy E 2 as a function of | n| 2 and find a good linearity, consistent with Eq. (3.11). In the determination of ξ f , to minimize the systematic effects due to excited state contamination at higher momenta, we only use the lowest four momentum vectors n = (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1) in the linear fit of E 2 (| n| 2 ) to Eq. (3.11). As seen in the figures, the fit results denoted by blue bands explain the data very well. The extracted value of ξ f = 6.41(11) from a pseudoscalar meson is in good agreement with the one from a vector meson, ξ f = 6.36 (14), and shows better precision. Therefore, for the tuning of lattice bare parameters we use ξ f from pseudoscalar mesons which are summarized in Table 3.

Tuning results
To determine the coefficients, a i , b i , and c i , we perform the simultaneous χ 2 fit of the numerical data in Table 3 to the functions in Eq. where the values of χ 2 per degrees of freedom are 1.72, 0.72, 0.23, respectively. In Appendix C we show some examples of the results of the fit in the two-dimensional spaces of the renormalized and bare parameters. Our interpretation of the above results requires that we comment on a few important features. First of all, renormalized anisotropies are somewhat larger than the bare anisotropies, which we interpret as a signal of the fact that the calculations are performed far from the weak coupling limit. Secondly, we find that the coefficients a 2 and b 1 are small, in particular, b 1 is zero within the statistical errors. In the quenched approximation, one would expect that the gauge and fermion anisotropies can be determined independently. The mild dependences of ξ f on ξ 0 g and ξ g on ξ 0 f are consistent with the fact that this part of the numerical study is performed in the regime of heavy quarks. Yet, we note that over the range of considered lattice parameters our results show a good linear dependence of the squared mass of a pseudoscalar meson M 2 P S on the bare quark mass m 0 , which is consistent with our use of Eq. (3.7) to extrapolate to the limit of vanishing physical mass for the quarks.

(3.15)
We will use these choices for the lattice parameters in measuring the physical properties of the field theory. Note that m * 0 falls slightly outside the range of masses used in this part of the study (see Table 3), and hence we expect some (small) residual quark mass and symmetry-breaking effects to be present in our physical simulations.

Numerical results: Finite temperature
From now on, the lattice bare parameters are fixed by the central values in Eq. (3.15) along with β = 2.0. We perform finite temperature calculations on the anisotropic lattices of N t × 16 3 and N t × 16 2 × 24. Simulation details and numerical results for these two lattices are summarized in Tables 4, 5 and 6. Two different values of N z are considered to estimate the systematic errors due to excited state contaminations in the calculations of screening masses. The algorithms for the generation of gauge ensembles have been discussed in Sec. 3.2.
Before we discuss the numerical results of finite temperature calculations in details, we perform a zero temperature calculation in order to check how well the tuned bare parameters are working. Using the ensemble of 128×16 3 in Table 5, we obtain ξ g = 6.29(4), ξ f = 6.1(2), and M 2 ps = 0.00517 (14). These results are compatible with the renormalized parameters of ξ = 6.3 and m 2 ps = 0.005, where the largest uncertainty occurs in the detemination of ξ f with ∼ 3%. Finite volume effects are expected to be negligible as the lattice volume is much larger than the size of the pseudo-scalar meson, m ps L ∼ 7.  Adopting anti-periodic boundary condition along the temporal direction, temperature is defined by T ≡ 1 Ntat . We will find it convenient to measure the temperature in units of the (pseudo-)critical temperature T c , discussed and measured in the next section.

Deconfinement crossover
As is the case for QCD with small number of quarks, our model is also expected to exhibit confinement at low temperature and form a quark-gluon plasma across the (pseudo-)critical temperature T c . Although the Polyakov loop is not an exact order parameter when the number of quarks is finite, it is widely used as an indicator of deconfinement. Following the method used in [26,27], we define the expectation value of the renormalized Polyakov where the bare Polyakov loop L 0 (T ) is related to the bare free energy F 0 (T ) as L 0 (T ) = exp(−F 0 (T )). The multiplicative renormaliztion constant is defined by Z L = exp(−∆F 0 ), which only captures the short distant physics and thus is independent on the temperature. As different choices of Z L denote different renormalization schemes, to incorporate the scheme dependence on the detemination of T c we impose a renormalization condition for a given temperature T R by L R (T R ) ≡ constant. We consider three renormalization schemes, defined by the conditions L R (N t = 24) = 0.9, L R (N t = 24) = 0.5, and L R (N t = 20) = 0.9 respectively. The results are shown in Fig. 6. The temperature T c is determined from the peak of the susceptibility of the Polyakov loop, χ(L R ) = ∂L R /∂T , denoted by dashed lines in the figure. Combining the statistical uncertainty and the systematic uncertainty of scheme dependences in quadrature, we find that T c a t = 0.0255 (25), or equivalently that N c t = 39 (4). As anticipated, we will measure temperatures in units of this T c in the following.

Temporal correlation functions
At zero temperature, the Euclidean two-point correlation functions of mesonic observables fall off with a single exponential at a large time so that the ground state energy of mesons can be extracted in a clear way in principle. In the finite temperature lattice calculations this process is affected by some limitations. Firstly, the maximum available physical temporal extent is limited by the inverse of the temperature. In addition, a single exponential analysis becomes subtle as the spectral function of mesons no longer exhibits a sharp peak at the  Figure 6. Renormalized Polyakov loop and their susceptibility. The renormalized Polyakov loops L R denoted by empty squares are obtained from the ensembles of N t × 16 2 × 24 with N t ranged over [16,56]. The solid curves are the interpolation of L R connected by cubic splines, while the dashed curves are the corresponding susceptibility χ(L R ), the derivatives of L R with N −1 t . Different colors are associated with different renormalization conditions, while the blue band denotes the (pseudo-)critical temperature T c with uncertainties as described in the text. mass of mesons. In this case, it is more desirable to investigate the correlation functions by themselves.
We introduce the normalized correlation function with the reference choice t = N t /2: We consider isovector pseudo-scalar, scalar, vector, and axial-vector mesons, where the corresponding interpolating fields are defined by respectively (flavour indices selecting non-singlet states are understood). In order to improve the statistics, we use stochastic wall sources [30] for the study of meson spectrum at finite temperature. Using these mesonic operators we compute the function C Nt/2 (t). In Fig. 7 we show the results of log C Nt/2 (t) for N t = 48 and 40, which exemplify the typical behaviors of C Nt/2 (t) below and near T c respectively. By comparing the two plots in Fig. 7 one can see that while at low temperature (N t = 48) the vector and axial-vector correlators are different, they become hard to distinguish from one another in proximity of T c (N t = 40). The overlap of C Nt/2 (t) between vector  and axial-vector mesons can be considered as an indication of the parity doubling in the vector channel and thus the restoration of the global SU (4) symmetry. By contrast, the situation for scalar and pseudo-scalar correlators is quite different, as we will discuss better by looking at spatial correlation functions in the next subsection, and indicates that at this temperature we do not yet see evidence of the restoration of the U (1) A symmetry. Notice that the correlation functions still satisfy the Weingarten's mass inequalities [31].

Spatial correlation functions
In contrast to the temporal correlation function, the spatial correlation function at finite temperature exhibits a single exponential decay at large time. The decay rate is called screening mass, as it defines the effective length scale associated with the excitation of mesonic operators in the medium [32]. At zero temperature the screening mass is equivalent to the meson mass, as the temporal and spatial correlation functions share the same spectral function. By using the meson interpolating fields in Eq. (4.3), we calculate the ensemble average of spatial correlators C(z) along the z-direction and extract the masses in units of a s using the analysis method described in Sec. 3.4. Notice that in our anisotropic lattice calculations the spatial and temporal lengths are measured differently. To have the consistent lattice unit of mass in a t , we therefore define the screening mass M S by multiplying ξ −1 to the measured spatial masses. In addition to the screening masses, we define the following normalized mass ratios for the vector channel, and for the scalar channel. These quantities are useful to quantify the level of parity doubling in the mass spectrum.
Our main results are presented in Table 5 and 6, as well as in Fig. 8 and 9. The error bar of each data point only represents the statistical uncertainty. We show explicitly the comparison between N t × 16 3 (black) and N t × 16 2 × 24 (red) lattices. The level of agreement of the two ensembles implies that there is no significant systematic uncertainty due to excited state contaminations.
By looking first at the vector and axial-vector masses, we see a plateau in R V above T c , which together with the change of behavior of the masses above T c strongly suggests that parity partners are degenerate and the global symmetry is effectively restored. There is small deviation from zero in the mass ratio at asymptotically large values of T , that may be the result of finite spacing, finite mass and possibly other small lattice artefacts.
In the case of the scalar channel, the plateau in R S appears at somewhat larger temperature, ∼ 1.5T c . This result may imply that the axial U (1) A and global SU (4) symmetries are restored at different temperatures. However, this is not conclusive, for several reasons. First of all, because we do not know what kind of transition is appearing in the underlying dynamical model, and it is likely that T c actually identifies a cross-over. But also because we do not know how much each of the lattice artefacts affects the results, and it might be that different observables are affected in different amounts by the finite quark mass, or the finite value of the coupling. A relevant discussion in the context of two-flavor QCD can be found in [14], for instance, where the numerical results strongly suggest that the symmetry restorations occur simultaneously in the massless limit.  In the finite temperature calculations, it is often suggested to plot the screening mass divided by the temperature as it shows linear dependency above T c . The results are shown in Fig. 10 and Fig. 11. The black dashed line corresponds to 2π which is associated with the Matsubara frequency for massless free quarks. For all mesonic channels, data points  approach the dashed line as the temperature increases and seem to form a plateau. However, they start to deviate from the plateau above 2T c , possibly as a consequence of the finite lattice spacing. 7 This suggests that in looking at R V and R S (and in general in discussing parity-doubling) one should not include in the physical very high temperatures, but rather restrict attention to T < ∼ 2 T c .
artefacts can significantly be reduced if highly improved lattice fermions being used.

Discussion
We collected numerical evidence of the fact that the high-temperature behavior of the SU (2) theory with N f = 2 Dirac fundamental fermions differs in three respects from the lowtemperature one. The numerical study of the Polyakov loop and its fit shows the existence of a pronounced peak in the susceptibility. Its position identifies a temperature T c , that we interpret in terms of the deconfinement (cross-over) temperature. While the study of the details of the transition would require a dedicated program, this result is accurate enough to allow us to clearly separate the high-T and low-T regimes, and concentrate on the symmetry properties of the physical spectrum above T c .
The study of temporal correlation functions shows that for T > T c the vector and axial-vector 2-point functions have compatible t-dependence, supporting the hypothesis that parity doubling is emerging at T c , and global symmetry is restored. This is confirmed by the study of spatial correlation functions, in which one clearly sees that the behavior of the screening masses of vectors and axial-vector mesons changes at T ∼ T c : while the two masses are different and depend on T in two different ways when T < T c , for T > T c the masses come close to one another, and, most importantly, show the same T -dependence.
This last observation suggests that the small splitting in the masses we observe is due to a combination of lattice artefacts (in particular finite spacing and finite quark mass). To confirm or disprove this statement, one would need to extend the study in this paper, and consider more than one value of the bare coupling and of the bare quark mass, in order to extrapolate them both to the physically relevant regime. By doing so, one might not only be able to show that the mass difference between vectors and axial-vectors vanishes, but also to study other properties of the transition itself, such as its order.
The numerical study of the scalar and pseudo-scalar masses, in which we focused on cleaner states that form a fundamental of SO(5), shows qualitative features that are in broad agreement with the restoration also of the axial U (1) A symmetry at high temperature. Our data on spatial correlation functions seems to suggest that this is taking place at a larger temperature T ∼ 1.5T c . This is also supported by the fact that in the temporal correlation functions we do not see the effect of parity doubling in the spin-0 correlators, for the same choice N t = 40 for which the vector and axial-vector correlators do agree with one another. This is the most striking element of novelty of this study, although it must be considered as preliminary.
This paper is to be understood as a first step in what is a potentially broad and extensive research program. The results obtained are in good agreement with what expected on field-theory grounds about the non-trivial behavior of this theory at high temperatures: it deconfines, and both the global SU (4) and axial U (1) A symmetries are restored. Two main sets of explorations are interesting to pursue in the future. On the one side, it is interesting to perform precision studies of this system, in which larger statistics, and a broader set of values of the lattice parameters, are used in order to establish whether the three transitions we identified are distinct (and in this case how to classify them, and precisely measure the critical temperatures), or whether they are just three manifestations of the broader phenomenology related to a cross-over.
On the other hand, it is also interesting to understand how the system reacts to the introduction of additional sources of symmetry breaking at the Lagrangian level. For example, the weak gauging of a subgroup of Sp(4) (as in phenomenological composite-Higgs models), is going to break the global symmetry of the model, and with it the large degeneracies of states. It would be useful to know how these phenomena depend on finite temperature. Closely related, although possibly simpler, is the question of what happens at finite µ: given that SU (2) is pseudo-real, this model is free of the traditional sign problem of similar models with larger gauge groups. It should hence be possible to attempt a more general study of the phase diagram as a function of both T and µ.
The richness of the field theory behavior of this model, the wide variety of its possible applications and the fact that this study shows that its thermal features are amenable to quantitative numerical studies, all contribute to making it an ideal environment in which to study highly non-trivial phenomena, which might shed light on many aspects of direct relevance to QCD, TC and composite-Higgs scenarios. In this paper we performed a first study along these lines, mainly aimed at collecting evidence of symmetry restoration at high temperature. We also discussed ways to improve our results, and suggested avenues for further investigation, which we will pursue in the future.
Given a 2-component spinor u we can build a 4-component Majorana spinor as In 4-component notation this ensures that λ L = Cλ R T and λ L = λ T R C = −λ T R C −1 . With these definitions in place, and making use of the fact that Grassmann variables anticommute, after some algebra one finds that which implies that the kinetic term can be written equivalently in terms of λ R as of λ L (the total derivative can be dropped), or equivalently one can write it in terms of the 4-component Majorana spinor λ (with an overall factor of 1 2 to avoid double counting). We specify now the model of interest in this paper, with SU (2) gauge symmetry. Starting from the 2-component spinors q i a , with i = 1 , · · · , 4 the flavor index and a = 1, 2 the color index, we can build four LH and four right-handed(RH) 4-component spinors as with j = 1 , · · · , 4. Notice that the charge-conjugation used for the RH spinors implies to lower the SU (2) indexes, as it turns a fundamental of SU (2) in its conjugate. The essential property of SU (2) is that this can be compensated by the ab antisymmetric tensor. One can define two Dirac spinors Q i a = q i a L + q i+2 a R , with i = 1, 2. We identify such Q i a with the Dirac spinors that form the fundamental matter fields of the SU (2) gauge theory. Because of the structure of the gamma matrices, the kinetic terms do not couple different chiralities, and hence we can write and hence we can write which is manifestly SU (4) × U (1) A invariant. Notice that besides the SU (2) t × SU (2) t , the SU (4) ∼ SO(6) group includes also the U (1) t B associated with baryon number. For completeness, we can explicitly verify that where the fact that Ω is antisymmetric comes from the antisymmetric ab . This is a Majorana mass, with M = m Ω, which breaks explicitly the symmetry to Sp(4). The U (1) t B is a subgroup of Sp(4) ∼ SO (5), hence the spectrum of composite states cannot be classified in terms of baryon number, as mesons and baryons are in common Sp(4) multiplets. In the case one gauges the baryon number, then the symmetry would be explicitly broken back to the familiar U (2) 2 .
B SU (4) and Sp(4) algebra. The 5 Goldstone bosons can be written as π(x) = 5 A=1 π A (x)T A , or explicitly as The maximal SO(4) ∼ SU (2) L × SU (2) R subgroup of the unbroken Sp(4) can be chosen to be generated by The T L generators satisfy the SU (2) L algebra T i L , T j L = i ijk T k L , and similarly T i R , T j R = i ijk T k R , while T A L , T B R = 0. These generators being all unbroken (in a vacuum aligned with Ω), this is the natural choice of embedding of the SO(4) symmetries of the Higgs field in the context of composite Higgs. The same model can be used also to describe traditional technicolor. In this case, the embedding of the Standard Model symmetries is based on the natural choice of generators of SO(4) t ∼ SU (2) t L × SU (2) t R as follows: In this case, one finds that (with the vacuum aligned with Ω) the breaking SU (2) t L × SU (2) t R → SU (2) t V emerges, and the unbroken generators are t A V = (t A L +t A T R ), or explicitly: The normalization is Tr t A V t B V = δ AB , as in this case we are writing the generators in the bifundamental representation.
The unbroken U (1) t B associated with baryon number is generated by T 15 = 1 √ 2 (T 3 L + T 3 R ), while the anomalous axial U (1) A is generated by

C Fit results of renormalized parameters
In this Appendix, we demonstrate how the fits in Section 3.5 work by showing the renormalized parameters in Eq. (3.7) along with the fit results of Eq. (3.13) in the two-dimentional slices of the measured and lattice parameters. See Fig. 12 and Fig. 13 for the fermion and gauge anisotropies, and see