Open-Boundary Conditions in the Deconfined Phase

TL;DR

This paper investigates open-boundary conditions at high temperatures to measure boundary effects, screening masses, and topological susceptibility in lattice gauge theories, revealing temperature-dependent boundary zone sizes and challenges due to topological freezing.

Contribution

It introduces the use of open-boundary conditions at high temperatures to study boundary effects and topological susceptibility, providing new measurements and insights into temperature-dependent phenomena.

Findings

Boundary effects are larger at T=1.5T_c than at T=0, smaller at T=2.7T_c.

Scalar and pseudo-scalar screening masses are measured, showing a mass gap at T=1.5T_c but not at T=2.7T_c.

Topological susceptibility estimates are consistent with literature at T=1.5T_c but limited at T=2.7T_c due to topological freezing.

Abstract

In this work, we consider open-boundary conditions at high temperatures, as they can potentially be of help to measure the topological susceptibility. In particular, we measure the extent of the boundary effects at $T=1.5T_c$ and $T=2.7T_c$. In the first case, it is larger than at $T=0$ while we find it to be smaller in the second case. The length of this "boundary zone" is controlled by the screening masses. We use this fact to measure the scalar and pseudo-scalar screening masses at these two temperatures. We observe a mass gap at $T=1.5T_c$ but not at $T=2.7T_c$. Finally, we use our pseudo-scalar channel analysis to estimate the topological susceptibility. The results at $T=1.5T_c$ are in good agreement with the literature. At $T=2.7T_c$, they appear to suffer from topological freezing, impeding us from providing a precise determination of the topological susceptibility. It still provides us with a lower bound, which is already in mild tension with some of the existing results.

Figures

Tables

Table 1. Table 1: Lattices used in this study. The spatial size of the lattice is denoted l 𝑙 l .

$N_s^3\times N_t$	$\beta$	$a\left[\mathrm{fm}\right]$	$\frac{l}{\sqrt{8t_0}}$	$\frac{T}{T_c}$	$n_{OBC}$	$n_{PBC}$
${64}^3\times 6$	$6.139$	$0.074$	$10.0$	$1.5$	$257$	$442$
${88}^3\times 8$	$6.335$	$0.056$	$10.4$	$1.5$	$430$	$533$
${112}^3\times 10$	$6.498$	$0.045$	$10.7$	$1.5$	$246$	$155$
Continuum extrapolations for OBC				$1.5$	-
Continuum extrapolations for PBC				$1.5$	-
${64}^3\times 4$	$6.253$	$0.063$	$8.52$	$2.7$	$533$	$512$
${88}^3\times 6$	$6.55$	$0.042$	$7.81$	$2.7$	$318$	$532$
${120}^3\times 8$	$6.778$	$0.031$	$7.86$	$2.7$	$533$	$532$
Continuum extrapolations for OBC				$2.7$	-
Continuum extrapolations for PBC				$2.7$	-
${102}^3\times 6$	$6.64$	$0.037$	$7.98$	$3.0$	$479$	$472$

Equations

$$Q=\int {\mathrm{d}}^{4}xq(x)=\frac{1}{32{\pi }^2}\int {\mathrm{d}}^{4}x{ϵ}^{\mu \nu \rho \sigma }\mathrm{Tr}\left(G_{\mu \nu }{G}_{\rho \sigma }\right),$$

$$S^{OBC}=-\frac{\beta }{3}\sum _{P}w(P)\mathrm{Tr}\left(1-P\right),$$

$$P_{\mu \nu }=U_{\mu }(n)U_{\nu }(n+\hat{\mu })U_{\mu }(n+\hat{\nu }{)}^{†}{U}_{\nu }(n{)}^{†}$$

$$w(P)=\{\begin{array}{l}1\text{ if }P\in \text{bulk}\\ \frac{1}{2}\text{ if }P\in x\text{-face}\end{array}.$$

$$⟨O{⟩}_{OBC}=\frac{1}{\mathcal{Z}}⟨{\gamma }^f|{\hat{T}}^{-(T-t)}O(t){\hat{T}}^{-t}|{\gamma }^i⟩.$$

$$⟨O{⟩}_{OBC}$$

$$\cdot ⟨E_m(\vec{q})|{\hat{T}}^{-(T-t)}O(t){\hat{T}}^{-t}|E_n(\vec{p})⟩$$

$$=\frac{1}{\mathcal{Z}}\sum _{n,\vec{p},m,\vec{q}}{\gamma }_{n,\vec{p}}^i{\gamma }_{m,\vec{q}}^{f\ast }{e}^{-(T-t)E_m(\vec{q})}{e}^{-t{E}_n(\vec{p})}$$

$$\cdot ⟨E_m(\vec{q})|O(t)|E_n(\vec{p})⟩,$$

$$⟨O{⟩}_{OBC}=⟨0|O|0⟩$$

$$+{\beta }_1{e}^{-(E_1-E_0)(T-t)}+\dots .$$

$${\dot{B}}_{\mu }(x,\tau )=D_{\nu }{G}_{\nu \mu }(x,\tau ),B_{\mu }(x,\tau ){|}_{\tau =0}=A_{\mu }(x),$$

$$D_{\mu }={\partial }_{\mu }+\left[B_{\mu }(x,\tau ),\cdot \right].$$

$$G_{\mu \nu }(x)={\left(\frac{-i}{8a^2}\left(Q_{\mu \nu }(x)-Q_{\nu \mu }(x)\right)\right)}_{AH},$$

$$Q_{\mu \nu }(x)=U_{\mu ,\nu }(x)+$$

$$+U_{-\mu ,-\nu }(x)+U_{-\nu ,\mu }(x),$$

$$O^r=\frac{1}{N_y{N}_z{N}_t}$$

$$\sum _{y=0}^{N_y-1}\sum _{z=0}^{N_z-1}\sum _{t=0}^{N_t-1}\sum _{x=r}^{N_x-r-1}O(x,y,z,t),$$

$$E_{st}=\frac{1}{4}{G}_{0i}^a{G}_{0i}^a,E_{ss}=\frac{1}{4}{G}_{ij}^a{G}_{ij}^a;$$

$$E_{ss}(r)=\alpha \exp \left(-m_{0^{+}}r\right)+\beta +\gamma a^2$$

$$U_{\mu \nu }(x,\tau )$$

$$-\frac{1}{4}{\delta }_{\mu \nu }{G}_{\sigma \rho }^a(x,\tau )G_{\sigma \rho }^a(x,\tau )$$

$$E(x,\tau )$$

$$E_{ss}(x,\tau )=\frac{1}{4}\left(U_{ii}(x,\tau )-U_{00}(x,\tau )\right)-2E(x,\tau ).$$

$$E_{ss}(x,\tau )=\frac{c_T(\tau )}{4}$$

$$-\frac{c_E(\tau )}{2}{\left\{G_{\mu \nu }{G}_{\mu \nu }\right\}}^R-c_1(\tau ),$$

$$c_T(\tau )$$

$$c_E(\tau )$$

$$\Delta E_{ss}^r=\underset{\tau \to 0}{\lim }\left(E_{ss}^r(\tau )-E_{ss}^{r_{bulk}}(\tau )\right).$$

$$q(x)=\frac{1}{32{\pi }^2}{ϵ}^{\mu \nu \rho \sigma }\mathrm{Tr}\left(G_{\mu \nu }{G}_{\rho \sigma }\right)$$

$${\chi }^r$$

$$\overline{q}(r)=\sum _{x=r}^{N_x-r-1}\sum _{y=0}^{N_y-1}\sum _{z=0}^{N_z-1}\sum _{t=0}^{N_t-1}q(x,y,z,t).$$

$$\chi$$

$$=\frac{1}{V}\frac{1}{Z}\int \mathcal{D}A\left(\int {\mathrm{d}}^{4}x\int {\mathrm{d}}^{4}yq(x)q(y)\right)e^{-S[A]},$$

$$E(\theta )=-\frac{1}{V}\ln \mathcal{Z},$$

$${\chi }_{OBC}(1.5T_c)t_0^2$$

Full text

††journal: Eur. Phys. J. C

1 Introduction

In general, finite-size systems differ from their infinite-volume counterpart. One of the most simple examples is the ”particle-in-a-box” whose momenta are quantised. Not only the compactness, but also the boundary conditions affect the system. There, different choices lead to different quantisation conditions. The only restriction on such choices is that the infinite volume physics needs to be recovered in the thermodynamic limit. This requirement satisfied, the only remaining differences are related to the convergence to the infinite volume limit. When the system is discretised, discretisation effects may also vary between different types of boundary conditions.

In some circumstances, such differences may be used as algorithmic tools to improve numerical simulations. A typical example of this is the use of open-boundary conditions (OBC) in lattice QCD, which have been introduced in Luscher2011 as means to reduce autocorrelations of the topological charge. These autocorrelations become critical as the continuum is approached. and are signaled by the freezing of gauge field ensembles in given topological sectors. In this example, instead of considering QCD with periodic boundary conditions (PBC), which leads to a discrete topological charge

[TABLE]

the idea is to impose OBC in at least one of the directions. In this system, $Q$ spans a continuum range of value. This then lifts the topological barrier responsible for the topological freezing and improves the sampling of the configuration space.

Having small autocorrelations is crucial to keep control of the statistical errors in Monte Carlo simulations Wolff:2003sm ; Schaefer2011 . A poor sampling of the topological charge affects in principle all observables, leading to finite volume effects (see Aoki:2007ka ; Bietenholz:2016ymo for practical examples). The situation is partially improved when considering QCD in the deconfined phase. For $T>T_{c}$, the order parameter which quantifies the variance of the topological charge, i.e. the topological susceptibility $\chi=\frac{\langle Q^{2}\rangle}{V}$, decreases with $T$. At asymptotically-high temperatures, it is even suppressed as $T^{-7}$ Gross:1980br . Nonetheless, for moderate temperatures, $Q\neq 0$ configurations still contribute in a non-negligible way to the path-integral. In this context, OBC may also be of interest at non-zero temperatures111For a very exploratory study, see Burnier:2017osu .. However, before being able to use them systematically, an analysis of the influence of temperature on the boundary effects remains to be done. This is the content of this study, which focuses on pure $SU(3)$ gauge theory, as dynamical matter is not expected to drastically change the results.

In section 2, we recall known facts about OBC and discuss our datasets and methodology. Then, in the spirit of the zero temperature analysis of Husung:2017qjz , we investigate in section 3 the typical length over which the boundary effects propagate, the \sayboundary zone. We observe a noticeable temperature dependence. These differences can be understood in terms of the temperature dependence of the lightest propagating states’ screening masses, which we study in section 4. As a by-product, we report in section 5 an extraction of the topological susceptibility from our rather large volumes simulations. We finally discuss our results in section 6.

2 Open-Boundary Conditions and Setup

Conventional lattice QCD simulations use (anti-)periodic boundary conditions in all directions, for the obvious reason that they minimise boundary effects. In this study, we consider the use of OBC in one of the spatial directions (taken for definiteness to be the $x$ direction). This amounts setting the field-strength tensor to zero outside the lattice. In this case the Wilson action reads Luscher2011

[TABLE]

where the sum runs over all the plaquettes

[TABLE]

whose corners are in the interval $[x=0,x=N_{x}-1]$. The $U$’s are the usual link variables and the quantity $w(P)$ is an integration weight

[TABLE]

A bulk point is a point in the interval $[1,N_{x}-2]$. A plaquette is on a $x$-face if it is not oriented along $x$ and all of its corners are at $x=0$ or $x=N_{x}-1$. As shown in Luscher2011 , the continuum limit of this theory has a trivial topology in field space; all the admissible fields are connected by local gauge transformations.

Such boundary conditions break translational invariance and introduce boundary effects. These effects may be understood as the propagation of excited states from the boundary. Here we summarise the core of the argument, following Guagnelli:1999zf ; Bruno2016 ; GarciaVera2017 .

For the sake of clarity, let us first recall the argument for OBC in the time direction; it straightforwardly transposes to OBC in the $x$ direction. To quantise our Euclidean theory, we write down a transfer matrix $\hat{T}=e^{-\hat{H}}$ with $\hat{H}$ the Hamiltonian, the Euclidean equivalent of the evolution operator. It evolves states between temporal slices. In particular, going from the state $\left|\gamma_{i}\right>$ at $t=0$ to the state $\left|\gamma_{f}\right>$ at $t=T$, and given an operator $O$ inserted at $t$, we can write

[TABLE]

To label our basis of states, we use the lattice version of the translation operators and get a basis consisting of $\left|E_{n}(\vec{p})\right>$, with $n$ labelling extra quantum numbers and $\vec{p}$ the momentum eigenstates. Inserting a complete basis of states, we can then write

[TABLE]

with $\gamma^{i,f}_{n,\vec{p}}=\left<E_{n}(\vec{p})\right|\gamma^{i,f}\left.\right>$. Now we see that the main contribution comes from the state with smallest energy. We then have a tower of exponentially suppressed corrections. More explicitly, using the fact that the main contribution to $\mathcal{Z}$ is $\gamma^{i}_{0}\gamma^{f*}_{0}e^{-E_{0}T}$ (obtained by setting $O(t)=\mathbf{1}$ in our expansion), we find

[TABLE]

with $\alpha_{1}$ and $\beta_{1}$ some matrix elements.

In other words, OBC do not project out directly on the vacuum state but are affected by states which propagate from the boundary. We also see that, at least in some limits, the corrections should be dominated by an exponential decay in the lightest state. We will take advantage of this in section 4. Note that this argument can be generalised to two-point functions Bruno2016 and higher-point functions.

In the case of OBC in the $x$ direction, the previous analysis can be repeated by replacing the slicing in the $t$ direction by a slicing in the $x$ direction when quantising the system. Then $H$ and $\hat{P}_{x}$ exchange roles, with $\hat{P}_{x}$ the translation operator in $x$. Modulo this, the derivation goes through.

To measure the topological correlators, we used the gluonic definition of the topological charge density. It requires some smoothing of the gauge fields, which was performed by using the gradient flow Luscher:2010iy . The fundamental gauge fields $A_{\mu}(x)$ are evolved to finite flow-time $\tau$, $B(x,\tau)$, using the flow equation

[TABLE]

[TABLE]

The associated smearing radius is $\sqrt{8\tau}$. It is implemented on the lattice by using the standard Wilson gauge action (Wilson flow). The integration is done using a third order Runge-Kutta algorithm with a step-size of $0.01$, which was tested to be small enough for the lattice parameters of this study.

The configurations we used are listed in table 1. The quenched configurations were generated using a heat bath and an overrelaxation algorithm. One update consists of one heat bath and four overrelaxation steps. To make sure that the configurations are sufficiently thermalised we discard configurations from the first $4000$ iterations. Configurations are measured every $500$ Monte Carlo steps to minimise the autocorrelations. Working with flowed configurations, we use the scale $t_{0}$ with the interpolation given in Francis:2015lha to convert to physical units. The statistical errors were estimated by using Jackknife resampling.

To compute the topological charge and energy density we used the clover-shaped field strength tensor

[TABLE]

where $AH$ is the projection on the traceless antihermitian part and $Q_{\mu\nu}$ is defined as

[TABLE]

with the plaquette discretisation $U_{\mu,\nu}$.

3 Boundary Zone

As explained in section 2, the presence of a boundary affects observables in the bulk, at least close to the boundary. The length of this \sayboundary zone depends on how the observables couple to the lightest propagating states. To quantify this effect and in order to compare it to the zero temperature case, we adopted the method of Husung:2017qjz . We compute the value of the clover action density as a function of the distance to the boundary and extract the length of the boundary zone, i.e. the length over which this observable is significantly different from its bulk value. In more detail, for lattices with OBC in the $x$ direction and some operator $O$, we define its sub-average at a distance $r$, inside a sub-volume of size $(N_{x}-2r)\times N_{y}\times N_{z}\times N_{t}$ from the boundary, by

[TABLE]

with $0\leq r<N_{x}/2-1$. For $r=0$, we expect the strongest dependence on the boundary excitations. By studying the $r$-dependence, we can then characterise the typical size of the boundary contamination.

At non-zero temperature, the clover action density leads to two independent gluon condensates Datta:2015bzm ; Wandelt:2016oym

[TABLE]

a \saymagnetic condensate $E_{ss}$ and an \sayelectric condensate $E_{st}$.

In figure 1, we show both densities at the reference flow-time $t_{0}$ for our different configurations at $T=1.5T_{c}$. All temperatures used in our study behave in a qualitatively similar way. First, we see as expected the existence of a boundary zone and an agreement between OBC and PBC in the bulk of the lattices, i.e. when $r$ is sufficiently large to suppress the effects of the boundary on the sub-volume. Then we see that the component which displays the largest boundary zone is $E_{ss}$. The reason is that it couples to a lighter state than $E_{st}$.

To compare different results in all fairness, we proceed to a continuum extrapolation of both condensates. In figure 2, we show this continuum extrapolation for $T=1.5T_{c}$ and three different radii. As reported in Luscher:2010iy , the region close to the boundary is affected by linear lattice spacing artifacts when Wilson’s action is used without further improvements. We evade this complication by computing our continuum extrapolation only in the region where the $O(a)$ corrections are negligible. This region turned out to be large enough for all purposes of this study.

Different temperatures are compared in figure 3, together with the zero temperature result of Husung:2017qjz . In this plot, we show the energy density normalised to its bulk value. We see that the length of the boundary zone depends on temperature. At $1.5T_{c}$ we find it to be about $50\%$ larger than at zero temperature while we find it reduced by $20\%$ at $2.7T_{c}$, consistently with our fixed lattice spacing results at $3T_{c}$. This is also consistent with the temperature dependence of the screening masses. Actually, the behaviour of the observables in the boundary zone gives a handle on these screening masses, which will be discussed in the next section.

4 Screening Masses

As explained in section 3, the boundary effects are controlled by the masses of the propagating states in the theory. In pure $SU(3)$ gauge theory at finite temperature, these are the screening masses Nadkarni:1986cz .

In this section, we will take advantage of the boundary effects to extract the lightest scalar and pseudo-scalar screening masses. In particular, as the lightest scalar mass is expected to be the lightest state in our system, its value controls the length of the boundary zone of section 3.

4.1 Scalar Screening Mass

The strong boundary contamination seen in the $E_{ss}$ channel in figure 1 suggests that it might be an appropriate probe to extract the scalar screening mass $m_{0^{+}}$, which will correspond to the lowest screening mass of the state which couples to $E_{ss}$. At zero temperature, it would be the lowest glueball state. Such a strategy was used in Chowdhury:2014kfa ; Chowdhury:2014mra to extract glueball masses.

To make $E_{ss}$ ultraviolet (UV) finite and be able to take the continuum limit, we study it at some finite flow-time. To have good control of our errors, we perform a simultaneous fit of the type

[TABLE]

with $r$ a radius in the boundary zone (see section 3). The constant $\beta$ has to reproduce the continuum bulk value and the $\gamma$ factor encodes the $a^{2}$ finite lattice spacing corrections.

We look for an intermediate range $r\in[r_{min},r_{max}]$ of values where we can extract a candidate mass $m_{0+}$. On the left-hand side of figure 4 we show the behaviour of $\Delta E_{ss}^{r}t_{0}^{2}$ for different flow-times (top panel) together with the extraction of the screening mass for different $r_{min}$ and different flow-times. We also checked that the results were not sensitive to the choice of $r_{max}$.

The extracted screening mass should be flow-time independent, being the mass of some states (the flow evolution will mix different operators but not change the operator basis), and we see that within our precision it is. Outside the plateau region, the masses differ but they do match once a plateau is reached. Typically, small flow-times lead to a worse signal. The reason lies in the smoothing effect of the flow. For larger flow-times, the errors are reduced and generally speaking overlap between states increases, as do their matrix elements. We can verify this by looking at the prefactor $\alpha$ of our exponential fit, normalised by the bulk value. This quantifies the strength of the interaction with the $0^{+}$ boundary state. We extract it using the same procedure as for the screening mass and report its flow-time dependence on the top panel on the right-hand side of figure 4. As expected, we see it growing with the flow-time. This also explains the behaviour of the boundary region as a function of the flow-time, which is shown in the lower panel on the right-hand side of figure 4. The more we flow, the stronger the interaction with the boundary gets and the larger the boundary zone becomes. This suggests that upon a good knowledge of the flow dependence of the observable under consideration, smaller flow-times are advantageous with respect to the boundary contaminations.

In this spirit, it is also instructive to perform the same mass extraction in the limit of zero Wilson flow. It serves two purposes. First, it allows checking the robustness of our results. Then, since $E_{ss}$ is directly related to the energy-momentum tensor $T_{\mu\nu}$, taking the zero flow-time limit provides a properly renormalised observable. This would, for example, be required to extract any running quantities, such as the matrix elements encoded in $\alpha$. More precisely let us consider Suzuki:2013gza

[TABLE]

We can write

[TABLE]

The flow dependence then reads, using the expansions of Suzuki:2013gza ,

[TABLE]

with $T_{\mu\nu}^{R}$ the renormalised field strength tensor and $\left\{G_{\mu\nu}G_{\mu\nu}\right\}^{R}$ the renormalised version of $G_{\mu\nu}G_{\mu\nu}$. The coefficients can be expanded perturbatively as

[TABLE]

with $g_{0}$ the bare coupling and $g^{\overline{\rm{MS}}}$ the running coupling in the $\overline{\rm{MS}}$ scheme (see also DelDebbio:2013zaa ). The coefficient $c_{\mathbf{1}}$ is a mixing with unity and is set to $c_{\mathbf{1}}(\tau)=\langle E(\tau,x^{bulk})\rangle$ where by $x^{bulk}$ we mean the value in the centre of the lattice in the case of OBC. This sets the vacuum expectation of the trace of the energy-momentum tensor to zero Husung:2017qjz . Equation (19) allows to obtain a renormalised quantity to study the screening mass,

[TABLE]

The zero flow-time extrapolation is shown in the top left plot of figure 4. To perform the extrapolation, we used a quadratic fit and checked that the result was insensitive to higher order corrections. An example at fixed radii is shown on figure 5. As expected, the extracted screening mass is compatible with the one obtained at other flow-times, as is shown in figure 6.

We also extracted the screening mass at $T=2.7T_{c}$, but did not extrapolate to zero flow-time; this is shown in figure 7. Note that the mass is noticeably larger at $2.7T_{c}$ than at $1.5T_{c}$; see section 4.3 for a discussion. Consistently, the errors are also larger at $2.7T_{c}$. It also explains why we did not proceed to a zero flow-time extrapolation. As we may see, the signal quickly worsens at small flow-time and the noise reduction associated to the flow is crucial to extract the mass. It is thus extracted at $t_{0}$.

4.2 Pseudo-scalar Screening Mass

Upon considering different operators, this method allows us to access the mass of the screening states of different quantum numbers. In this section, we will proceed with the mass determination of the pseudo-scalar screening state. One of the first continuum operators which come to mind and couples to the pseudo-scalar sector is the topological charge density

[TABLE]

Unfortunately, we cannot proceed with its integrated average, as we did in section 4.1 with the energy density, as $\langle Q\rangle=0$ in our case, with $Q$ the topological charge (eq. (1)); in other words we are in the sector $\theta=0$, with $\theta$ the QCD $\theta$-angle. To circumvent this issue, we consider the two-point function of $q$ over different sub-volumes,

[TABLE]

where we defined an averaged $\overline{q}$,

[TABLE]

In other words, $\chi^{r}$ is the average of topological charge square over a sub-volume. As the notation suggests, this quantity is related to the topological susceptibility, see section 5. 222Note however that it is not the same quantity as is often used to study the topological susceptibility with OBC. Often, one considers the point-to-all two-point function of the topological charge, with the source far from the boundary (one can even average over different sources), see Bruno:2014ova for more details. We show the $r$ dependence of this quantity in figure 8 for various ensembles (some ensembles were omitted for the sake of clarity). Let us start discussing the ones at $1.5T_{c}$ (left-hand side of figure 8). As expected, we see again a boundary zone in the case of OBC and a saturation away from it. In the very centre of the lattice, $\chi^{r}$ displays a characteristic \saybump. This feature is inherited from the behaviour of the correlator $\langle q(x)q(0)\rangle$ around $x=0$ (see Giusti:2018cmp for a detailed discussion).

The results at $2.7T_{c}$ display the same global features as the ones at $T=1.5T_{c}$, with a notable exception: $\chi^{r}$ does not completely saturate; we observe a drift in its plateau value. We understand this effect as a manifestation of topological freezing (the lattices at $2.7T_{c}$ are finer than the one at $1.5T_{c}$), see section 5 for a discussion.

In all cases, to extract the screening masses, we are only interested in the exponential decay from the boundary. We use the same strategy as in the previous section. As the pseudo-scalar is heavier, we perform the extraction at flow-time $t_{0}$ to have a good signal to noise ratio; as in the case of the scalar mass at $2.7T_{c}$, the signal quickly deteriorates for smaller flow-times. We show the results in figure 9 and 10. The errors are comparable to the ones obtained for the scalar at $2.7T_{c}$, as the masses are of similar magnitude. We also checked that the masses are (within the statistical uncertainties) independent of the maximal radius used for the fit, as long as this radius is taken within the plateau region of $\chi^{r}$.

4.3 Discussion

All masses determined in this study are shown in physical units in figure 11. As expected, being less symmetric, the pseudo-scalar state is heavier than the scalar state. Whilst certainly present at $T=1.5T_{c}$, the difference is not statistically significant at $2.7T_{c}$. This is an indication of dimensional reduction; at high temperature, the scalar and pseudo-scalar are expected to become degenerate Datta:1998eb .

On the same plot, we also show the values obtained in Meyer:2008dt by measuring the asymptotic behaviour of the energy density. The qualitative behaviour is the same but we observe a shift of about $15\%$. Even if part of this discrepancy can presumably be explained by the fact that the results of Meyer:2008dt are at fixed lattice spacings and other systematics, an intrinsic difference between the two methods cannot be excluded.

Setting this aside, the data of Meyer:2008dt indicates that the main contribution to the pseudo-scalar mass is linear in $T$, as would be expected from perturbation theory at high temperatures. Taking this for granted, we can estimate that the scalar screening mass becomes heavier than the lightest glueball at around $2T_{c}$. This should correspond to the temperature at which the boundary zone becomes strictly smaller than the zero temperature one. And indeed, the fact that the scalar screening mass at $1.5T_{c}$ is lighter than the lightest $T=0$ glueballs and that the $T=2.7T_{c}$ scalar screening mass is heavier is consistent with what was reported in figure 3 about the length of the boundary zone.

5 Topological Susceptibility

Actually, the topological charge density square presented in section 4.2 also allows us to extract a value for the topological susceptibility. Indeed, as the continuum topological susceptibility is defined to be

[TABLE]

with the $E(\theta)$ the vacuum energy at non-zero $\theta$,333We recall that Euclidean $SU(3)$ gauge theory at non-zero $\theta$ can be described by the Lagrangian $\mathcal{L}=L^{SU(3)}_{\theta=0}+i\theta q$.

[TABLE]

we expect the plateau value of $\chi^{r}$ to give the topological susceptibility.

We show the obtained values in figure 8. The $1.5T_{c}$ case is the most straightforward and leads to a clean signal. We perform a global fit on our three ensembles of the type $f(a)=c_{1}\cdot a^{2}+\chi$ to remove the discretisation effect and extract the constant $\chi$. For PBC, we fit from the boundary up to $r_{max}=2\sqrt{8t_{0}}$. In the case of OBC, we excluded the data in the boundary zone. Correspondingly, we used values in the range $[1.3\sqrt{8t_{0}},2\sqrt{8t_{0}}]$. It gives us measurements for the topological susceptibility

[TABLE]

which are in good agreement with $\chi(1.5T_{c})t_{0}^{2}=2.25(12)\cdot 10^{-5}$ of Giusti:2018cmp and $\chi_{t}(1.5T_{c})t_{0}^{2}\in[1.5\cdot 10^{-5},4.4\cdot 10^{-5}]$, the global fit of reference Borsanyi:2015cka . They are also consistent with the fixed lattice spacing results of Berkowitz:2015aua .

At $2.7T_{c}$ the situation is more intricate, as it is not clear that $\chi^{r}$ saturates to a plateau. The trouble comes from two reasons. First, as it was already discussed in Bonati:2017woi in the context of a toy model with OBC, $\chi^{r}$ can present a slow convergence as a function of $r$. This is what we see in the centre of the lattice. Then, the behaviour of the PBC configurations seems to indicate that our ensembles are partially frozen.

First, in figures 12 and 13, we show the history of the topological charge for PBC. We observe a clear difference between the two temperatures. At $2.7T_{c}$, the topological transitions are highly correlated. Then, in figures 14 and 15, we focus on our finest configurations at those two temperatures. We also display the results obtained when restricting ourselves to the $Q=0$ sector, i.e. by artificially freezing our lattices. Of course, at $1.5T_{c}$, the effect of freezing is drastic, as a lot of topological transitions are still to be expected. What is more interesting is the qualitative behaviour of $\chi_{r}$. For small sub-volumes, the value for the topological susceptibility is not so far from the unfrozen value but decreases for larger sub-volumes. It is consistent with the observation reported in Schaefer2011 , where it was observed that the topological charge measured on sub-volumes is less autocorrelated than the total charge. It is also intimately tied to the fact that the freezing of the topological charge is only a finite volume effect, one of the key ideas behind master field simulations Giusti:2018cmp ; Luscher:2017cjh , where very large volumes are generated and the ensemble average is obtained by summing over decorrelated sub-volumes. The fact that we do not get the correct value for the topological susceptibility is only due to our volumes being too small to perform sub-volume averages in fixed sectors.

This discussion can also be applied to the $2.7T_{c}$ case. However, there, the same kind of behaviour is present in the \sayunfrozen case, which seems to indicate some partial freezing of our ensembles. This seems to be confirmed by the behaviour of the topological charge history of the PBC, which displays long correlations between jumps in topological sectors of the same sign; it shows correlations of at least $300$ configurations. Unfortunately, the OBC shows a similar kind of behaviour. This stresses the point that OBC are not a remedy to the topological freezing but only a potential improvement.

This discussion shows that no reliable estimate of $\chi(2.7T_{c})$ can be extracted from figure 8 without further investigations of the autocorrelations. In particular, it confirms that even with these relatively large volumes, an extraction of the topological susceptibility from $Q^{2}$ cannot be done reliably without a much larger statistics.

Note also that the way we extracted the topological susceptibility in this work is not the only way to do it. One can also consider the point-to-all integrated two-point function of the topological charge, with the source far-away from the boundary Bruno:2014ova

[TABLE]

with

[TABLE]

The first sum in (31) is an average over sources that are far enough from the boundary while the second sum is a genuine integration. The quantity $2l$ is the number of source points which are averaged over. While we did not systematically study this quantity, we did check that our method is consistent with this definition. In figure 16, we show $\chi^{r}$ and $\chi^{2pt}(r,12)$ for our largest configurations at $T=1.5T_{c}$. Note that the choice $l=12$ is presumably not optimal and the error bars associated to $\chi^{2pt}$ can presumably be reduced by tuning this parameter. We see that both methods are consistent; a careful study of their different systematics and how they relate is left as a potential interesting outlook.

6 Conclusion

In this study, we started a first systematic investigation of OBC at high temperature. The main difficulty in dealing with OBC is the presence of boundary effects. In section 3, we investigated the typical propagation length of these effects and compared it to the zero temperature results of Husung:2017qjz . At $T=1.5T_{c}$, the boundary zone is larger than at $T=0$, while it is smaller at $T=2.7T_{c}$ and $T=3.0T_{c}$. These differences can be understood in terms of the temperature dependence of the mass of the lightest state in our system, namely the scalar screening mass. Actually, the boundary contamination gives us means to measure this screening mass, giving results which are consistent with the already existing literature (see section 4.3). In particular, we predict that the scalar starts to be heavier than the $T=0$ lightest glueball at around $T=2T_{c}$. It tells us that the use of OBC in the region $T\in[T_{c},2T_{c}]$ is more delicate than at $T=0$ but becomes gradually easier at temperatures above $2T_{c}$. This is potentially useful as it is the interesting range of temperatures to measure the topological susceptibility, for example Moore:2017ond . Moreover, we do not expect the situation to change drastically in full QCD, in the deconfined phase.

We also used the boundary effects in the pseudo-scalar channel to estimate the corresponding screening mass. We measured a sizable mass gap between the scalar and pseudo-scalar at $T=1.5T_{c}$. Moreover, we could confirm that this gap reduces at higher temperature, which is an expected signal of the dimensional reduction taking place at high enough temperatures.

As a by-product of the pseudo-scalar analysis, we could extract a precise measurement of the topological susceptibility at $T=1.5T_{c}$, which is in good agreement with the recent results of Giusti:2018cmp . Finally, the same analysis at $T=2.7T_{c}$ exhibits some signs of topological freezing. A potential interesting outlook consists in studying quantitatively how the autocorrelation time depends on the lattice spacing and temperature and how it compares to the master field approach Giusti:2018cmp . Even so, it shows again that even with rather large volumes, the determination of the topological susceptibility is delicate. This supports the recent efforts Giusti:2018cmp ; Jahn:2018dke , which have been undertaken to reassess the robustness of high-temperature studies of the topological susceptibility. In particular, a careful reconsideration of the finite size effects on its determination, even in the quenched case, is called for.

Acknowledgements

The authors wish to thank M. Lombardo for an interesting discussion, A. Francis, M. Bruno and S. Sharma for comments on a previous version of this manuscript and particularly R. Sommer for his help and suggestions. A.F. thanks A. Rago, M, Hansen, F. Knechtli and N. Husung for insightful discussions and A. Jaquier for his technical support. The work of A.F. is supported by the Swiss National Science Foundation. The authors acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 \sayStrong-interaction matter under extreme conditions– project number 315477589 – TRR 211. Numerical calculations have been made possible through PRACE grants at CSCS, Switzerland, as well as grants at the GCS and at NIC-Jülich, Germany. These grants provided access to resources on Piz Daint at CSCS and JUWELS at NIC. Some of the numerical simulations have been performed on Fidis and Deneb, which are part of the EPFL HPC centre (SCITAS).