The θ -dependence of the SU( N ) critical temperature at large N

,


INTRODUCTION
Non-abelian SU(N ) gauge theories exhibit a nontrivial dependence on the parameter θ, coupling the integer-valued topological charge to the standard pure-gauge Yang-Mills action S YM = 1 2g 2 d 4 x Tr {G µν (x)G µν (x)}, whose behavior in the large-N limit is of particular theoretical and phenomenological interest, due to its relation with the solution of the U(1) A puzzle via the Witten-Veneziano mechanism [1][2][3][4][5][6][7][8].Among the several intriguing implications of θ-dependence, it is particularly interesting to investigate how the phase structure of Yang-Mills theories is modified in the presence of a non-zero θ [9][10][11][12][13][14].
At θ = 0, it is well known that SU(N ) pure Yang-Mills theories with N ≥ 3 undergo a genuine first-order phase transition (for N = 2, instead, it is of second order [15]) when crossing the critical temperature T c [16][17][18][19][20][21][22][23][24][25][26], which is associated to the spontaneous breaking of the Z N center symmetry enjoyed by S YM , with the Polyakov loop as an order parameter.The critical temperature T c is required to stay finite in the large-N limit, since the Witten-Veneziano solution to the U(1) A problem assumes the confining phase to survive up to N = ∞.This behavior has been confirmed by lattice results [19][20][21][22][23].
In the presence of a non-vanishing θ, since Q preserves the center symmetry too, one expects this picture to remain true (at least as long as θ < π), with the critical temperature changing as a function of θ.On general grounds, T c (θ) is an even function of θ, being Q a CPodd quantity, whose Taylor expansion around θ = 0 is usually parameterized, at leading order, as: where T c = T c (0) is the θ = 0 critical temperature, and where the O(θ 2 ) coefficient R represents the curvature of the phase diagram in the T − θ plane.This picture has been well-verified for N = 3 from lattice simulations [27][28][29], and a numerical computation of R, based on the imaginary-θ dependence of the critical temperature, can be found in [27,28].On the other hand, no lattice investigation of the large-N behavior of R can be found in the literature.
In this respect, in Ref. [27] a prediction for the large-N behavior of R is given, based on the known properties of the θ-dependence of the Yang-Mills free-energy density (here V is the four-dimensional volume) in the broken and the symmetric phases.In particular, the argument of [27] assumes the presence of a first order transition (which is indeed the case for N ≥ 3) and relies on the fact that, by definition, at the critical point the free energies in the confined (c) and deconfined (d) phases are equivalent: f c (T c (θ), θ) = f d (T c (θ), θ).Introducing the customary parameterization of the Taylor expansion of f (T, θ) around θ = 0 up to O(θ) 2 , with the topological susceptibility, maintaining the equivalence condition as θ changes, f c (T c (θ), θ) = f d (T c (θ), θ), implies: arXiv:2312.12202v2 [hep-lat] 8 Feb 2024 where ∆χ = χ c − χ d is the jump of the topological susceptibility at the critical point, while ∆ϵ = ϵ d − ϵ c is the latent heat of the transition, with both quantities computed at θ = 0.The relation in Eq. ( 6) implies that R ∼ O(1/N 2 ), as shown by the following argument [27].In the confined phase the topological susceptibility approaches a finite large-N limit, as required by the Witten-Veneziano argument, and as it has been extensively checked by lattice simulations [30][31][32][33][34][35][36][37][38][39][40][41][42], while higher-order terms in the θexpansion are expected [43], and have been verified from the lattice [39][40][41]44], to be suppressed as powers of 1/N 2 , or faster.In the deconfined phase, instead, lattice simulations have shown the free energy to approach, as expected, the functional behavior predicted by the Dilute Instanton Gas Approximation (DIGA) [45,46]: where χ(T ) ∼ T −N is exponentially suppressed in N [21,[47][48][49][50][51][52][53][54][55].Putting these facts together, it is thus clear that ∆χ ∼ O(N 0 ), and since the latent heat is proportional to the number of degrees of freedom, ∆ϵ ∼ O(N 2 ), we thus expect, from Eq. ( 6), R = ∆χ/(2∆ϵ) ∼ O(1/N 2 ).However, although this relation between R and the features of the deconfinement transition has been well verified from lattice simulations for N = 3 [56], no direct probe of this prediction at large N , nor of the large-N scaling of R, can be found in the literature.Given the lack of lattice studies of the θ-dependence of T c apart from the N = 3 case, the goal of this work is to go beyond the state of the art discussed so far, and pursue a direct lattice investigation of the curvature R of the critical temperature at large N , with the aim of checking the prediction of [27] in Eq. ( 6) in the large-N limit.To do so, we will perform simulations for N = 4 and N = 6 and for several imaginary values of θ, and we will determine R from two different and complementary strategies.The first one relies on a θ = 0 calculation of ∆χ and ∆ϵ to determine R via Eq.( 6).Our second strategy, instead, relies on determining the critical temperature as a function of the imaginary θ, so that R can be directly obtained from its definition in Eq. ( 2).Finally, we will discuss the large-N behavior of R combining our new results with the pre-existing ones for N = 3.
This manuscript is organized as follows.In Sec. 2 we will present our numerical setup.In Sec. 3 we will present our numerical results for R for SU(4) and SU(6), and we will discuss the large-N limit of the curvature combining our determinations with the N = 3 one.Finally, in Sec. 4 we will draw our conclusions.

NUMERICAL SETUP
In this section we will describe our lattice setup, namely, the adopted discretization of the action, the adopted algorithms used for Monte Carlo (MC) simulations, and the strategy pursued to identify the critical point and calculate the critical temperature and other relevant observables.

A. Lattice action
We discretize the SU(N ) pure-gauge action on a N 3 s × N t lattice with periodic boundary conditions in all directions, where the temperature is given by T = (aN t ) −1 , with a the lattice spacing.Concerning the θ = 0 theory, we adopted the standard Wilson lattice action: where β = 2N/g 2 is the bare coupling and is the product of the SU(N ) gauge link matrices along the plaquette in the (µ, ν) plane based on the site x of the lattice.
Concerning the theory at non-zero θ, as it is well known, doing direct MC simulations at real values of θ is impossible due to the sign problem introduced by the topological term, which is purely imaginary: A common strategy to deal with this issue, which has been widely employed in lattice simulations targeting the study of θ-dependence [27, 28, 36, 39-41, 44, 57-66], is to perform imaginary-θ simulations.Assuming analyticity around θ = 0, it is possible to perform an analytic continuation towards imaginary values of θ with the change of variables θ I = iθ, so that the action in (9) becomes purely real and the sign problem due to the θ-term is avoided.
Concerning the topological term, we define it in terms of the so-called clover discretization, which is the simplest lattice discretization of the topological charge Q, defined in Eq. ( 1), with a definite parity: where the Levi-Civita tensor with a negative index is defined as ε (−µ)νρσ = −ε µνρσ .In the end, our total lattice action looks like: The discretization of the charge in Eq. ( 10) is not integer valued, and in particular its value is related to the physical integer one, on average, by a multiplicative renormalization, which stems from UV fluctuations at the scale of the lattice spacing [35,67].This, in particular, leads to a multiplicative renormalization of the lattice parameter θ L as follows: where θ I is the physical imaginary-θ.Concerning the topological susceptibility, instead, a naive definition based on Q clov would suffer also for an additional additive renormalization due to short-distance contact terms [31,68]: Such additive renormalization would eventually become dominant in the continuum limit compared to the physical signal.
A customary choice to get rid of renormalization effects is to compute the topological charge on smoothened configurations.Smoothing is a procedure to dampen UV fluctuations at the scale of the lattice cut-off ∼ 1/a while leaving the relevant IR topological fluctuations intact.Many different smoothing algorithms have been proposed, such as cooling [69][70][71][72][73][74][75], gradient flow [76,77] or stout smearing [78,79], all giving agreeing results when properly matched to each other [75,80,81].
In this work we adopted cooling for its simplicity and numerical cheapness.A single step of cooling consists of a lattice sweep where each link is aligned to its local staple (computed always at θ = 0), in order to locally minimize the Wilson action.After cooling, Z ≃ 1, and the distribution of Q (cool) clov shows sharp peaks very close to integer values.Thus, we assign to each configuration an integer topological charge via the following procedure [30]: where "round{x}" denotes rounding x to the closest integer, and α is defined as the minimum of ⟨(αQ The typical values we found for α were of the order of ∼ 1.05-1.1.Since, after smoothing, one also has M ≃ 0, the quantity χ = ⟨Q 2 ⟩ /V is a proper lattice definition of the topological susceptibility, with a proper scaling towards the continuum limit.Moreover, smoothing allows us to perform a non-perturbative computation of the constant Z needed to renormalize the lattice parameter θ L : cf. Eq. (12).In all explored cases, we found that the obtained results for Q, for χ and Z were stable after ∼ 10 cooling steps, thus we always computed the lattice rounded charge after 30 cooling steps.

B. Monte Carlo algorithms
In most cases, to generate gauge configurations we adopted the standard local over-relaxation (OR) [82] and heat-bath (HB) [83,84] algorithms.More precisely, our single MC updating step is given by 4 lattice sweeps of OR, followed by a lattice sweep of HB, where both local updates have been implemented à la Cabibbo-Marinari [85], i.e., updating all the N (N − 1)/2 diagonal SU(2) subgroups of SU(N ).
Local algorithms for SU(N ) gauge theories are well known to suffer, when approaching the continuum limit, from a severe increase of the autocorrelation time of topological quantities, usually known as "topological critical slowing down".Such problem further worsens as N is increased, eventually leading to a complete freezing of the MC evolution of the topological charge density even for moderate coarse lattice spacings when N is sufficiently large [18,30,42,[86][87][88][89][90][91][92][93][94][95].Since we are interested in measuring the θ-dependence of the critical temperature, topological freezing can potentially introduce sizeable systematic effects in R.Moreover, topological freezing makes it very difficult to obtain a reliable estimation of the topological susceptibility at large N , which is necessary to accurately check Eq. ( 6).Since, as we will show in the following, for the finest lattice spacings explored at our largest value of N , the MC histories of the topological charge are almost frozen, it is clear that in those cases a strategy to mitigate topological freezing is needed to keep possible systematics affecting R under control.
To this end, for those simulations points, we adopted the Parallel Tempering on Boundary Conditions (PTBC) algorithm.This algorithm was first proposed [96] and employed [66,97] in 2d CP N −1 models, and later implemented for 4d SU(N ) gauge theories too [41,98,99].
The PTBC algorithm involves the simultaneous simulation of N r copies of the system differing only for the boundary conditions imposed on a small 3-dimensional cubic L 3 d region, orthogonal to a spatial direction of the lattice, referred to as the "defect".One of the replicas, say the one labelled by r = 0, has standard periodic boundary conditions, while, for the other replicas (r > 0), each link crossing the defect gets multiplied in the Wilson action by a factor c(r), interpolating between c(0) = 1 and c(N r − 1) = 0.The topological term is instead not altered in any replica.
The last replica (r = N r − 1) has open boundary conditions along the defect, thus suffering for much smaller autocorrelation times of the topological charge [90,100].The fast MC evolution of the topological charge achieved in the open replica is then transferred towards the periodic one, since, after each replica has been independently updated via the standard MC step earlier introduced, the swap of each couple of adjacent replicas (r, r + 1) is proposed.Swaps are proposed via a Metropolis test whose acceptance probability is: where S (r) W denotes the Wilson action in the presence of the boundary conditions of the replica r, and U (r) denotes the gauge fields of the replica r.
To ensure the efficiency of the swapping procedure, the c(r) coefficients are tuned so that the probability of accepting the swaps during the Metropolis test is approximately constant for each pair of exchanged replicas (r, r+ 1), i.e., ⟨p(0, 1)⟩ ≈ ⟨p(1, 2)⟩ ≈ • • • ≈ ⟨p(N r − 2, N r − 1)⟩, which allows a given configuration to perform a random walk among the different replicas with different boundary conditions.
The advantage of using the PTBC algorithm is that it allows to exploit the much smaller autocorrelation times of the topological charge achieved with open boundary conditions, while at the same time bypassing several complications introduced by them.Indeed, this algorithm allows to perform measures directly in the periodic replica, where no issue related to the unphysical contributions of the boundary is present (which would require larger lattices to be kept under control), and where translation invariance is unbroken, so that a notion of global topological charge is preserved.Thus, all observables will always be computed on the periodic lattice.

C. Identification of the critical point
The order parameter of the phase transition is given by the expectation value of the Polyakov loop, which we discretize as follows: where ⃗ x is the space position in each 3-dimensional slice of the lattice at fixed Euclidean time t.As usual, we will study the expectation value of its modulus, ⟨|L|⟩.
To look for the critical point, we will follow the socalled "fixed-N t "approach.Namely, for a given lattice, we change the temperature T by changing the bare coupling β according to 1/T = a(β)N t , and we identify he critical point by looking at the behavior of the susceptibility of the Polyakov loop as a function of β.Indeed, on a finite volume, χ L presents a peak of finite width in correspondence of the critical temperature, associated to the critical coupling β c .
To determine more precisely such value of β, we interpolated the values of χ L (β) computed with direct simulations using the standard "multihistogram method" [101].After the multihistogram analysis, we obtain a smooth curve for χ L (β).We then fitted the peak of the interpolated curve with a Lorentzian function centered around β c .The uncertainties on the fit parameters were estimated from O(1000) binned bootstrap resamplings of the original data.A proper systematic uncertainty has been also associated to the variation observed upon varying the fit range around the peak.
To pass from the critical β to the critical temperature, we need to set the scale.To this end we used the determinations of the string tension in lattice units a(β) √ σ reported in [22].In particular, we determined the critical temperature in units of the string tension as follows: Once the critical point is identified, we will be interested in two observables computed at β c (which will be computed only for θ = 0): the latent heat and the discontinuity of the topological susceptibility.In both cases, we will need to compute differences between expectation values in the confined and deconfined phases.To this end, one needs a clear separation between the two phases, i.e., two clearly distinct peaks in the histogram of |L|, in order to assign a configuration to either a phase or the other.In practice, thus, this requires to work on large enough volumes, and to compute the quantities of interest only considering the gauge configurations with |L| lying in the range of values that are typical of each phase.More details about this procedure will be given in Sec. 3.
Concerning the jump of the topological susceptibility, we computed the discontinuity as follows: where Q is the rounded charge defined in Eq. ( 15), and where ⟨O⟩ c/d denote the expectation values in the confined (c) and deconfined (d) phases.
Regarding the latent heat, instead, we computed the jump at the critical temperature of the trace anomaly (ϵ − 3p)/T 4 c .Indeed, this coincides with the discontinuity of the internal energy ϵ, since the pressure p is continuous.The trace anomaly can be computed on the lattice using the thermodynamical identity: where ) is the partition function.Rewriting Eq. ( 22) in terms of lattice quantities yields: The latent heat ∆ε can, thus, be related to the difference in the expectation values of the mean plaquette between the two phases: where

RESULTS
This section is devoted to the discussion of our results for R obtained for N = 4 and 6 both from the imaginaryθ dependence of the critical temperature and from the determination of the latent heat of the transition and of the gap in the topological susceptibility measured at θ = 0.In both cases we will discuss both the extrapolation towards the continuum limit, and the extrapolation towards the large-N limit.
A. Results for SU (4) We start our discussion from the determination of R in the case N = 4.We considered 3 values of the temporal extent, namely, N t = 5, 6, 8.For each lattice, we performed simulations for a few values of θ L , and for each value of the imaginary-θ we performed simulations for several values of β, in order to identify the critical point.For this value of N , we observed a reasonable number of topological fluctuations in all explored points, thus, all simulations were performed with the standard local algorithm, without relying on the parallel tempering.Examples of MC histories of the topological charge are shown in Fig. 1 For each value of θ L , we observe that the Monte Carlo histories of the Polyakov loop, when the critical point is approached, exhibit the typical bistability which characterizes the coexistence of the broken and the symmetric phases, cf.Fig. 2. This, in turn, results in a peak of the Polyakov susceptibility χ L as a function of β in correspondence of the critical coupling β c , which is computed by means of the multihistogram analysis explained in Sec. 2 C. As expected, the peak of χ L (β), corresponding to the critical β c , moves towards higher couplings as θ L is increased, see Fig. 3.
All determinations of β c (θ L ) were done on lattices with aspect ratios N s /N t = 4, except for the simulations performed for N t = 8, where we used N s /N t = 3.Since we checked that, for our lattice with N t = 6, and for the smallest and largest values of θ L explored, the results obtained for β c (θ L ) with aspect ratios N s /N t = 3 and N s /N t = 4 were perfectly compatible, we conclude that finite size effects are well under control.All determinations of the critical values of β found for N = 4 are reported in Tab.I.
After determining the critical coupling β c as a function of θ L , we determined T c / √ σ according to Eq. ( 20), and the physical imaginary-θ parameter θ I = Z(β c )θ L according to the determinations of Z(β c ) obtained from Eq. (17).We remark that the renormalization constants Z(β c (θ L ) were obtained from dedicated θ = 0 simulations performed on hypercubic N 4 s lattices (i.e., deep in the confined phase) for β ≃ β c .All our results for Z(β c ), a(β c ) (σ) and T c / √ σ, θ I are reported in Tab.I.   I: Summary of results obtained for N = 4.The renormalization constants Z(βc(θL)) were obtained from dedicated θ = 0 simulations performed for β ≃ βc(θL) on a 16 4 lattice.The lattice spacings a(βc) √ σ were obtained interpolating determinations of [22].The critical temperatures are reported with two statistical errors: the first one is due to the uncertainty on βc, the second one combines in quadrature the uncertainties on βc and a(βc) √ σ.
We are now ready to compute the curvature R. Let us start from the computation of this quantity from the imaginary-θ dependence of the critical temperature.We fitted our data for T c / √ σ as a function of θ I according to the analytic continuation for imaginary values of θ of Eq. ( 2), which reads: An example of the imaginary-θ fit of the critical temperature is shown in Fig. 4.
Concerning statistical errors on R, it is important to stress here that such best fit was carried out taking into account the errors on θ I , inherited from the ones on Z(β c ), and the errors on the data of [22] for the string tension).In particular, a bootstrap analysis done to take into account correlations among a(β c ) √ σ reveals that the joint fluctuations of a(β c ) √ σ just have the effect of overall shifting the data for T c (θ I ), with a negligible impact on their curvature, eventually leading to a small uncertainty on R.
Concerning instead systematic errors, since we are dealing with a truncated Taylor expansion, one needs to check that systematic effects coming from the neglected terms are under control.Thus, to check that contributions from higher-order terms are indeed small, we performed both parabolic and quartic fits in θ I for various choices of the fit range, cf.Fig. 5.In all cases, systematics due to higher-order powers of θ I are well under control within our typical statistical errors, see Fig. 5 and Tab.II.Our final results for R from the imaginary-θ fits of the critical temperature are reported for each N t in Tab.III.
We now move to the determinations of R from θ L = 0 simulations alone, using Eq.(6).In order to reliably compute the jump of the topological susceptibility according to Eq. ( 21) and the latent heat according to Eq. ( 24), which require to compute differences of expectation values between the confined and the deconfined phases, we performed (for all values of N t but the smaller) dedicated θ = 0 simulations at β ≃ β c on lattices with a larger spatial volume compared to those employed for the identification of the critical point, the purpose being of obtaining a sharper separation between the two phases.
An example of the histories of the Polyakov loop obtained from such simulations on a 48 × 8 lattice is shown in the bottom panel of Fig. 6.As it can be observed, by starting from ordered/random configurations we obtain histories which explore only a single phase, thus allowing an unambiguous calculation of the expectation values ⟨O⟩ c,d in each of the two phases.TABLE III: Summary of our final results for the curvature R for N = 4 obtained from dedicated θ = 0 simulations aiming at determining the jump of the topological susceptibility ∆χ and the latent heat ∆ϵ.For completeness we also report the final results for R obtained from the imaginary-θ fit.
On the other hand, for the coarsest lattice spacing explored (corresponding to N t = 5) the latent heat turned out to be large enough to allow an unambiguous separation of the two phases already on the volume employed to compute β c , as it can be observed by inspecting the histogram of the Polyakov loop in that case, shown in the top panel of Fig. 6.Therefore, in that case we just fixed two cuts, L  cut , to, respectively, the confined and the deconfined phase.We also verified that our results for ∆ϵ(N t = 5) were largely insensitive to the specific choice of these cuts as long as their values were varied within the depleted region separating the two peaks, cf.top panel of Fig. 6.
Our results for the jump of the topological susceptibility ∆χ/T 4 c , the latent heat ∆ϵ/T 4 c and the curvature R = ∆χ/(2∆ϵ) according to Eq. ( 6) for all explored points are reported in Tab.III, along with the final results for R obtained from the imaginary-θ fit of T c .We now proceed to perform the extrapolation of both determinations of R towards the continuum limit.In both cases, we took the continuum limit assuming O(a 2 ) corrections, i.e., since 1/N t = aT c , according to the law: where R represents the continuum-extrapolated result.
The continuum limits of R from both data sets are shown and compared in the top panel of Fig. 7. Our final continuum results are: R(N = 4) = 0.0095 (11) (from imaginary-θ), ( 28) R(N = 4) = 0.0119(16) (from latent heat).(29) As it can be observed, obtained determinations are perfectly compatible within the errors once the continuum limit is taken, confirming that, also for N > 3, the predicted relation between the curvature of the T − θ phase diagram and the jump at the critical point of the internal energy and the topological susceptibility holds.
Given the very good agreement between these two results, we also tried to take a common continuum limit, where a common fit parameter was used for R in Eq. ( 27).The two data set can be nicely fitted jointly with such a  fit function, and the best fit yields a reduced chi-squared of 2.3/3, corresponding to a p-value of 52%, and a combined continuum limit: R(N = 4) = 0.01025(92) (combined). ( The combined fit is shown in the bottom panel of Fig. 7.

B. Results for SU(6)
The computation of R for N = 6 goes essentially along the same lines outlined for SU(4), meaning that we looked for the critical couplings β c (θ L ) for several values of the imaginary-θ parameter, and for 3 values of the temporal extent N t = 5, 6, 8.However, compared to the N = 4, there are two important differences that we want to stress before presenting our numerical results.
The first difference relies on the adopted volumes, which are smaller than those chosen for N = 4.As a matter of fact, to determine β c (θ L ) for N = 6 we adopted in all cases lattices with aspect ratios N s /N t = 2. Indeed, being ∆ϵ ∼ N 2 , the strength of the first order for N = 6 was such that, on lattices with larger volumes, we were not able to unambiguously identify the critical point.Nonetheless, on the basis of general theoretical arguments relying on large-N volume reduction [102], we expect finite size effects to become milder at large N .This fact has also been verified in previous lattice studies of the SU(N ) phase transition [22,23].Moreover, the curvature R is essentially the leading-order deviation from 1 of T c (θ)/T c (0), thus, finite size effects for this quantity are expected to largely drop in the ratio.In conclusion, we are confident that, even though smaller spatial volumes were used compared to the N = 4 case, our determinations of R in the N = 6 case are still reliable.
The second difference with respect to the N = 4 case, is that, as already discussed in Sec. 2 B, for N = 6 the simulations performed for larger values of β suffer from a severe topological freezing.In particular, we found this problem to affect all simulations performed on the 16 3 ×8 lattice, for all values of β and θ L employed.Thus, all simulations performed on this lattice for N = 6 were done using the PTBC algorithm, which provides a tremendous gain in terms of the autocorrelation time of the topological charge.To be more precise, for each value of β and θ L explored at N t = 8 we simulated N r = 30 replicas of the lattice, and boundary conditions were changed on a cubic defect of size L d = 3.The tempering parameters used to change boundary conditions were tuned through short test runs to ensure an approximately constant replica swap rate (18) of ∼ 30%.
Let us start exactly by the latter point.In Fig. 8, we show the history of the topological charge obtained from the standard combination of local algorithms, and from the parallel tempering.As it can be observed, the PTBC algorithm allows to obtain a substantial gain in terms of the observed number of fluctuations of Q with the same numerical effort1 .Remarkably enough, we also found that, with the PTBC algorithm, it is possible to achieve a reduction of about a factor of 3 of the autocorrelation time of L. This is probably due to the large correlation between L and Q. Examples of MC histories of the Polyakov loop are shown in Fig. 9, where for the sake of clarity we just report the MC evolution of L obtained with parallel tempering.
We now move to the determination of R from the imaginary-θ fit.In Fig. 10   TABLE IV: Summary of results obtained for N = 6.The renormalization constant Z(βc(θL)) were obtained from interpolating the results for Z(β) reported in Refs.[40,41].The lattice spacings a(βc) √ σ were obtained interpolating determinations of [22].The critical temperatures are reported with two statistical errors: the first one is due to the uncertainty on βc, the second one combines in quadrature the uncertainties on βc and a(βc) √ σ.  obtained from the imaginary-θ fit, in this case the situation is different.Indeed, we estimated from a bootstrap analysis that, for our N = 6 data, these two sources of uncertainty on T c / √ σ gave comparable contributions to the final errors on R. Therefore, the uncertainties on R reported in Tab.V have been computed from the sum of these two errors, which were added in quadrature.
Concerning instead the determination of R from the latent heat, again we relied on dedicated θ = 0 simulations performed at β ≃ β c and using the same volumes adopted for the identification of the critical point.As already discussed in the previous section, we assigned configurations satisfying suitable cuts on the Polyakov loop to either one of the two phases.Also in this case we verified that changing the values of these cuts within the depleted region between the two peaks has a negligible impact on the latent heat, as any corresponding variation of ∆ϵ is much smaller than the obtained statistical errors on this quantity.(37) 0.00800(50) 16 8 0.1463(46) 12.17(76) 0.00608 (39) 0.00598 (39) TABLE VI: Summary of our final results for the curvature R for N = 6 obtained from dedicated θ = 0 simulations aiming at determining the jump of the topological susceptibility ∆χ and the latent heat ∆ϵ.For completeness we also report the final results for R obtained from the imaginary-θ fit.
An example of the histogram of the Polyakov loop sand of the cuts used to separate the two phases is shown in Fig. 13, while our determinations of ∆χ/T 4 c , ∆ϵ/T 4 c and R = ∆χ/(2∆ϵ) for all explored points are collected in Tab.VI.
The continuum limits of R(N = 6) from both data sets, again taken assuming leading O(a 2 ) corrections, is shown in the top panel of Fig. 14.Our final continuum results are: R(N = 6) = 0.00385(75) (from imaginary-θ), (31) R(N = 6) = 0.00524(61) (from latent heat).(32) Again, as it can be observed, obtained determinations for R from the two different strategies are perfectly compatible among them in the continuum limit, confirming again the validity of the prediction in Eq. ( 6) for N = 6 too.
Given the very good agreement between these two continuum extrapolations, also in this case we consider a combined continuum limit of the two data sets, which yields a reduced chi-square of 3.18/3, corresponding to a p-value of 37%, and a combined limit: R(N = 6) = 0.00469(47) (combined).( The combined limit is shown in the bottom panel of Fig. 14.

C. The large-N limit
This section is devoted to the discussion of the large-N limit of R, which we extrapolate using our determinations for N = 4, 6 and the pre-existing N = 3 results.Our goal is in particular to compare the large-N scaling obtained from the imaginary-θ and the latent heat determinations.
The results for R as a function of N are summarized in Tab.VII.Concerning the N = 3 determinations, we considered the value of R(N = 3) obtained from the imaginary-θ fit of Ref. [27], and the values of ∆ϵ(N = 3) and ∆χ(N = 3) obtained, respectively, in Refs.[24] and [56] to compute R = ∆χ/(2∆ϵ) for N = 3. TABLE VII: Summary of our final results for the curvature R for N = 4 and 6, along with the N = 3 results obtained in Refs.[28] (imaginary-θ) and [24,56] (latent heat and jump of the susceptibility).
As explained earlier in the introduction, we expect the following asymptotic large-N scaling: In order to verify this prediction, we considered the following 3 ansätze, which we used to perform a best fit of the two data sets reported in Tab.VII: We start our investigation by comparing the best fits obtained from these ansätze using the imaginary-θ data set.Assuming a pure power-law behavior R ∼ 1/N c , with c left as a free parameter, the best fit of R(N ) as a function of 1/N yields c = 2.20 (24), which is in perfect agreement with the 1/N 2 predicted leading large-N scaling for R.This finding thus justifies a best fit fixing c = 2, which gives a perfectly compatible results for the pre-factor R, both when the determination for N = 3 is included and excluded from the fit.Finally, fitting our data assuming the presence of further 1/N 4 higher-order corrections also gives perfectly agreeing results with the pure 1/N 2 fit, as we find a compatible result for R, and a value for the O(1/N 4 ) coefficient R(1) = 0.22 (25) which is compatible with zero within the error.
The perfect agreement among the best fitting curves obtained from the 3 different ansätze employed can also be clearly seen from the top panel of Fig. 15, where one can clearly observe the 3 curves and their corresponding error band to be perfectly overlapping bands, and also from the comparison reported in Fig. 16 of the different values of R obtained from the various employed ansätze, which are all in very good agreement.
Perfectly agreeing results and perfectly analogous considerations can be drawn by repeating this systematic study of the N -dependence of R from the latent heat data set, as it can be seen from the comparison in the bottom panel of Fig. 15, and from the obtained results for R from these determinations reported in Fig. 16.
As a matter of fact, from the comparison of obtained results for R from the various employed ansätze, cf.Fig. 16, we quote as our final results: R = 0.159(4), (from imaginary-θ), R = 0.177 (14), (from latent heat).
Summarizing, we can conclude that our numerical results for N = 4, 6, along with the previously-obtained N = 3 determinations, show that R(N ) follows the predicted leading large-N scaling, thus confirming the conjectured relation between the curvature of the phase diagram in the T − θ plane and the topological features of the θ = 0 critical deconfinement transition.
One intriguing outcome of our results regarding the leading large-N scaling of R is the confirmation that ∆χ = χ c − χ d ∼ O(N 0 ), as expected on the basis of general theoretical arguments, since χ c ∼ O(N 0 ), according to the Witten-Veneziano mechanism, and ), according to the DIGA prediction.However, while the former is a very well established fact from lattice simulations [38,40,41,103], the former has not been directly verified, even though there is numerical evidence that χ d /χ c = χ(T + c )/χ(T − c ) is rapidly suppressed as a function of N [21,49].Since to obtain ∆χ we computed χ d = χ(T + c ) at θ = 0, a byproduct of our large-N investigation of R is the first direct lattice study of the N -dependence of χ d , which we are able to perform combining our continuum results for N = 4 and 6 with those obtained for N = 3 in Ref. [24].The extrapolation towards the continuum limit of our values of χ d /T 4 c for SU(4) and SU(6), as well as that of the SU(3) results reported in Tab. 2 of Ref. [56] are shown in the left panel of Fig. 17 The result for N = 3 was obtained extrapolating towards the continuum limit the determinations in Tab. 2 of Ref. [56].lent description of these data.The best fit parameters turn out to be A = 0.446 (38) and b = 0.451 (25).
The one-loop DIGA prediction, however, overestimates χ d /T 4 c by almost an order of magnitude, as it can be observed from the comparison in Fig. 17.On the other hand, the decay width of the exponential, although not exactly compatible, is in qualitative good agreement with lattice results.As a matter of fact, performing an exponential best of the numerical DIGA prediction according to A DIGA exp(−b DIGA N ) yields A DIGA = 4.7(3) and b DIGA = 0.65(2) (technical details about the numerical calculation of the DIGA prediction can be found in App.A).This is not surprising, as similar conclusions have been reached by comparing the lattice results for the temperature dependence of the topological susceptibility in the deconfined phase with the DIGA prediction both in quenched QCD and in the full theory with dynamical fermions [51-53, 55, 103, 104]

CONCLUSIONS
In this work we have presented a large-N lattice investigation of the O(θ 2 ) coefficient R, parametrizing the θ-dependence of the critical deconfinement transition of SU(N ) gauge theories at leading order: This quantity can be related to two features of the deconfinement transition at θ = 0, namely, the latent heat and the jump of the topological susceptibility, R = ∆χ/(2∆ϵ), implying also R ∼ O(1/N 2 ).
In order to verify this prediction, we performed a direct lattice calculation of R for N = 4 and 6 using two different and complementary strategies.The first one relied on the fit of the imaginary-θ dependence of the critical temperature, which was achieved by performing simulations in the presence of a non-vanishing imaginary-θ.The second strategy instead is based on the calculation for vanishing θ of the latent heat and of the jump of χ at the critical point, and thus relied just on θ = 0 simulations.
In both cases, for the finest lattice spacings explored at N = 6, due to the significant topological freezing suffered from standard simulations, we adopted the recentlyproposed state-of-the-art PTBC algorithm to obtain a sizable reduction of the auto-correlation time of the topological charge.
For example, in Ref. [51] one-loop DIGA result for χ(N = 3) at T = Tc were found to underestimate the lattice one by an order of magnitude.The origin of this discrepancy can be traced back to the fact that the authors employed the four-loop renormalized coupling, while here we considered just the one-loop expression.
Our numerical results show that, after the continuum limit is taken, the two different determinations of R perfectly agree among each other, thus confirming the predicted relation between R, ∆ϵ and ∆χ for several values of N .Moreover, combining our new N = 4, 6 results with the pre-existing N = 3 determination of R, ∆χ and ∆ϵ, we were able to study the large-N behavior of R(N ), finding perfect agreeing results with the predicted leading 1/N 2 scaling already for N ≥ 3.As a byproduct of this investigation, we were also able to show that the topological susceptibility in the deconfined phase computed at the critical temperature is exponentially suppressed as a function of N , as predicted by the DIGA.While the DIGA decay width qualitatively agrees with lattice data, the DIGA pre-factor of the exponential overestimates the lattice result by an order of magnitude.
Centro de Excelencia Severo Ochoa CEX2020-001007-S and, partially, by grant PID2021-127526NB-I00, both funded by MCIN/AEI/10.13039/501100011033.This work has been partially supported by the project "Non-perturbative aspects of fundamental interactions, in the Standard Model and beyond" funded by MUR, Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN), Bando 2022, grant 2022TJFCYB (CUP I53D23001440006).Numerical simulations have been performed on the MARCONI machine at Cineca, based on the agreement between INFN and Cineca, under projects INF22 npqcd and INF23 npqcd.

10000 FIG. 1 : 1000 FIG. 2 :N = 4 , 20 FIG. 3 :
FIG. 1: Evolution of the topological charge Q defined in Eq. (15) for the largest β explored at θL = 0 and θL = 20 on a 24 3 × 8 lattice for N = 4.The horizontal MC time is expressed in terms of the number of single MC steps defined at the beginning of Sec. 2 B.

FIG. 5 :
FIG.5: Study of the systematic effects on R(N = 4) due to the truncation of the Taylor expansion in θI .The solid lines and the shaded areas represent our final results.Empty points were ignored, as they are associated with best fits with a p-value smaller than 5%.

FIG. 6 :
FIG. 6: Top panel: history and histogram of the Polyakov loop L for N = 4 obtained on a 20 3 ×5 for β ≃ βc(θ = 0).The two solid lines correspond to the cuts L < L (c) cut and L > L used to assign a configuration to the confined or deconfined phase.Bottom panel: histories of the Polyakov loop L for N = 4 obtained on a 48 3 × 8 for β ≃ βc(θ = 0) and starting respectively from an ordered/random configuration.In both cases, the horizontal MC time is expressed in terms of the number of single MC steps defined at the beginning of Sec. 2 B.

FIG. 7 :
FIG.7: Extrapolation towards the continuum limit according to(27) of the determinations of the curvature R(N = 4) in Tab.III by, respectively, fitting the two data sets separately (top), and imposing a common continuum limit (bottom).
we show how the peaks of χ L (β) shift towards larger couplings as θ L is increased.Our determinations of β c (θ L ) for N = 6 are collected in Tab.IV, along with the corresponding values of T c / √ σ.A few examples of the imaginary-θ fit of T c / √ σ for N = 6 results are shown instead in Fig. 11.Again, in order to assess possible systematic effects on the curvature due to the truncation of the Taylor series, we performed both parabolic and quadratic fits in θ I for various choices of the fit range.All obtained results are collected in Tab.V.It is worth recalling that, as already explained in Sec. 3 A, the critical temperature T c (θ L )/ √ σ has two sources of uncertainties: the error on β c (θ L ), and the error on a(β c (θ L )) due to the scale setting, which were treated separately in our analysis due to the existing correlations among a(β c (θ L )).Unlike what was found for N = 4, where the first source of uncertainty gave by far the dominant contribution to the final error on

= 8 ,= 8 , 500 FIG. 8 := 8 ,= 8 , 100 FIG. 9 :
FIG. 8:Comparison of the evolutions of the topological charge Q defined in Eq. (15) obtained with PTBC and with the standard algorithm for the largest β explored at θL = 0 and θL = 15 on a 16 3 × 8 lattice for N = 6.The horizontal MC time is expressed in both cases in terms of the number of stadard MC steps defined at the beginning of Sec. 2 B. For the PTBC algorithm this means, in practice, to multiply the number of parallel tempering steps by the number of replicas.

N = 6 ,= 15 FIG. 10 : 15 FIG. 11 :
FIG. 10: Peaks of the Polyakov loop susceptibility χL(β) as a function of the bare coupling β observed on a 12 3 × 6 lattice for 3 different values of θL = 0, 15, 20 and N = 6.The lighter shaded areas represent the multihstogram interpolations of the MC data, while the darker solid lines are the results of the best fits of the interpolated data according to a Lorentzian ansatz f (β) = A1/[1 + (β − βc) 2 /A2].Filled points represent the position of the critical points.

FIG. 12 :
FIG.12: Study of the systematic effects on R(N = 6) due to the truncation of the Taylor expansion in θI .The solid lines and the shaded are represent our final results.Empty points were ignored, as they are associated with best fits with a p-value smaller than 5%.

FIG. 13 :
FIG. 13: History and histogram of the Polyakov loop L for N = 6 obtained on a 10 3 × 5 for β ≃ βc(θ = 0).The two dashed lines correspond to the cuts L (c) cut and L (d) cut used to assign a configuration to the confined or deconfined phase.The horizontal MC time is expressed in terms of the number of single MC steps defined at the beginning of Sec. 2 B.

FinalFIG. 16 :
FIG. 16:Comparison of the values obtained for the prefactor of the O(1/N 2 ) coefficient of R(N ) using fit functions(35), (36) and (37) to fit our results for R(N ) in Tab.VII.

FIG. 17 :
FIG. 17: Left panel: continuum extrapolations of our determinations of χ d /T 4 c = χ(T + c )/T 4 c for N = 4 and 6, along with the extrapolation of the N = 3 ones found in Tab. 2 of Ref. [56].Right panel: behavior of the continuum-extrapolated determinations of χ d /T 4 c as a function of N .The dashed line represents an exponential best fit of χ d /T 4 c as a function of N , as predicted by the DIGA.

TABLE II :
Results of the imaginary-θ fit of Tc/ √ σ as a function of θI for N = 4 and for various fit ranges, and considering both a parabolic and a quartic fit function in θI .Results marked with * were not considered when assessing the systematic error on R because of the unreliability of the fits. R

TABLE V :
Results of the imaginary-θ fit of Tc/ √ σ as a function of θI for N = 6 and for various fit ranges, and considering both a parabolic and a quartic fit function in θI .Results marked with * were not considered when assessing the systematic error on R due to the truncation of the Taylor series because of the unreliability of the obtained best fits.
(36)om panel: comparison of the large-N scaling of the determinations of R(N ) obtained from the imaginary-θ fit and from the latent heat.Dashed and dotdashed lines represent the best fit of both data sets according to fit function(36).
(37)15: Study of the large-N scaling of R. Top panel: Solid, dashed and dotted lines represent, respectively, the best fit of our data for R(N ) in Tab.VII according to fit functions(35),(36)and(37).