Open-boundary conditions in the deconfined phase

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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=1.5T_c$$\end{document}T=1.5Tc and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=2.7T_c$$\end{document}T=2.7Tc. In the first case, it is larger than at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=0$$\end{document}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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=1.5T_c$$\end{document}T=1.5Tc but not at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=2.7T_c$$\end{document}T=2.7Tc. Finally, we use our pseudo-scalar channel analysis to estimate the topological susceptibility. The results at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=1.5T_c$$\end{document}T=1.5Tc are in good agreement with the literature. At \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=2.7T_c$$\end{document}T=2.7Tc, they appear to suffer from topological freezing, which prevents us from providing a precise determination of the topological susceptibility.


Introduction
In general, finite-size systems differ from their infinitevolume 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 a e-mail: adrien.florio@epfl.ch b e-mail: okacz@physik.uni-bielefeld.de c e-mail: lmazur@physik.uni-bielefeld.de typical example of this is the use of open-boundary conditions (OBC) in lattice QCD, which have been introduced in [1] 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 (1) 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 [2,3]. A poor sampling of the topological charge affects in principle all observables, leading to finite volume effects (see [4,5] 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 χ = Q 2 V , decreases with T . At asymptotically-high temperatures, it is even suppressed as T −7 [6]. Nonetheless, for moderate temperatures, Q = 0 configurations still contribute in a nonnegligible way to the path-integral. In this context, OBC may also be of interest at non-zero temperatures. 1 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 Sect. 2, we recall known facts about OBC and discuss our datasets and methodology. Then, in the spirit of the zero temperature analysis of [8], we investigate in Sect. 3 the typical length over which the boundary effects propagate, the "boundary 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 Sect. 4. As a byproduct, we report in Sect. 5 an extraction of the topological susceptibility from our rather large volumes simulations. We finally discuss our results in Sect. 6.

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 [1] where the sum runs over all the plaquettes 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 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 [1], 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 [9][10][11].
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 matrixT = e −Ĥ withĤ the Hamiltonian, the Euclidean equivalent of the evolution operator. It evolves states between temporal slices. In particular, going from the state |γ i at t = 0 to the state |γ f at t = T , and given an operator O inserted at t, we can write To label our basis of states, we use the lattice version of the translation operators and get a basis consisting of |E n (p) , with n labelling extra quantum numbers and p the momentum eigenstates. Inserting a complete basis of states, we can then write with γ i, f n,p = E n (p) γ i, f . 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 Z is with α 1 and β 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 Sect. 4. Note that this argument can be generalised to two-point functions [10] 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 andP x exchange roles, withP 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 [12]. The fundamental gauge fields A μ (x) are evolved to finite flow-time τ , B(x, τ ), using the flow equatioṅ The associated smearing radius is √ 8τ . 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 [13] 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 where AH is the projection on the traceless antihermitian part and Q μν is defined as with the plaquette discretisation U μ,ν .

Boundary zone
As explained in Sect. 2, the presence of a boundary affects observables in the bulk, at least close to the boundary. The length of this "boundary 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 tempera-ture case, we adopted the method of [8]. 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 ) × N y × N z × N t from the boundary, by with 0 ≤ r < N x /2 − 1. For r = 0, we expect the strongest dependence on the boundary excitations. By studying the rdependence, 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 [14,15] a "magnetic condensate" E ss and an "electric condensate" E st . In Fig. 1, we show both densities at the reference flowtime 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 Fig. 2, we show this continuum extrapolation for T = 1.5T c and three different radii. As reported in [12], 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 Fig. 3, together with the zero temperature result of [8]. 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   Fig. 2 Continuum extrapolation of the magnetic clover density for three different radii r = 0.6, 0.8, 1.8 √ 8t 0 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.

Screening masses
As explained in Sect. 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 [16].
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 Sect. 3.

Scalar screening mass
The strong boundary contamination seen in the E ss channel in Fig. 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 [17,18] 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 with r a radius in the boundary zone (see Sect. 3). The constant β has to reproduce the continuum bulk value and the γ factor encodes the a 2 finite lattice spacing corrections. We look for an intermediate range r ∈ [r min , r max ] of values where we can extract a candidate mass m 0+ . On the left-hand side of Fig. 4 we show the behaviour of E r ss t 2 0 for different flow-times (top panel) together with the extraction of the screening mass for different r min and different flowtimes. 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 α 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   . 3 Comparison of the normalised clover action between different temperatures. We report in this figure the zero temperature results of [8] in blue. We observe the length of the boundary zone to depend on temperature. At 1.5T c , we see that the boundary effects propagate over a larger distance than at zero-temperature. We take as a conservative estimate of the boundary zone at 1.5T c a length of l 1 Fig. 4 Top left: Continuum extrapolated (E r ss − E bulk ss ) for different flow-times together with its τ → 0 limit. Already qualitatively, one can see that there is an exponential decay, whose exponent does not seem to be sensitive to the flow-time, whilst its prefactor does. Bottom left: Extraction of the effective mass as a function of the minimal radius used in the exponential fit. We see that when the parameter saturates to a plateau, different flow-times lead to the same prediction, as expected. Note that our errors seem to be overestimated for large r min ; we do not correct for this. Top right: Normalised prefactor of the exponential. This quantifies the interactions with the boundary states and increases with flow-time. This is due to the smoothing effect of the flow evolution; generically it increases the overlap between states. Bottom right: Corresponding effect on the boundary zone, its length increases with the flow-time as the bulk states interact more and more strongly with the boundary states on the right-hand side of Fig. 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 Fig. 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 μν , 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 α. More precisely let us consider [19] We can write The flow dependence then reads, using the expansions of [19], with T R μν the renormalised field strength tensor and G μν G μν R the renormalised version of G μν G μν . The coefficients can be expanded perturbatively as with g 0 the bare coupling and g MS the running coupling in the MS scheme (see also [20]). The coefficient c 1 is a mixing with unity and is set to c 1 (τ ) = E(τ, x bulk ) 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 [8]. Equation (19) allows to obtain a renormalised quantity to study the screening mass,  The zero flow-time extrapolation is shown in the top left plot of Fig. 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 Fig. 5. As expected, the extracted screening mass is compatible with the one obtained at other flow-times, as is shown in Fig. 6.
We also extracted the screening mass at T = 2.7T c , but did not extrapolate to zero flow-time; this is shown in Fig. 7. Note that the mass is noticeably larger at 2.7T c than at 1.5T c ; see Sect. 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 .

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 Unfortunately, we cannot proceed with its integrated average, as we did in Sect. 4.1 with the energy density, as Q = 0 in our case, with Q the topological charge (Eq. (1)); in other words we are in the sector θ = 0, with θ the QCD θ -angle. To circumvent this issue, we consider the two-point function of q over different sub-volumes, where we defined an averaged q, In other words, χ 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 Sect. 5. 2 We show the r dependence of this quantity in Fig. 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 Fig. 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, χ r displays a characteristic "bump". This feature is inherited from the behaviour of the correlator q(x)q(0) around x = 0 (see [22] 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: χ 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 Sect. 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 pseudoscalar 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 Figs. 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 χ r .

Discussion
All masses determined in this study are shown in physical units in Fig. 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 [24].
On the same plot, we also show the values obtained in [23] 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 [23] are at fixed lattice spacings and other systematics, an intrinsic difference between the two methods cannot be excluded.
Setting this aside, the data of [23] 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 Fig. 3 about the length of the boundary zone.

Topological susceptibility
Actually, the topological charge density square presented in Sect. 4.2 also allows us to extract a value for the topological susceptibility. Indeed, as the continuum topological susceptibility is defined to be Fig. 8 Topological charge density square. The legend's labels correspond to N t . For readability, we show only a subset of our data. At T = 1.5T c (left-hand side), everything behaves as expected. The topological charge density square converges when integrated from the bulk and saturates to a constant value, which we can identify with the topological susceptibility. The OBC have the same bulk behaviour but suffer from exponentially suppressed contributions from the boundary states.
The T = 2.7T c case (right-hand side) is more interesting. We see that even the PBC charge density does not saturate. It can be interpreted as an indication of topological freezing, as it is known that the charge density over a sub-volume is less autocorrelated than the total charge [3]. The OBC presents a similar pattern, calling for a more careful analysis of their autocorrelation time Fig. 9 Extraction of the pseudo-scalar screening mass from the boundary pollution at T = 1.5T c . The x-axis is the radius at which we start our single-exponential fit. We extract the mass from the plateau value with the E(θ ) the vacuum energy at non-zero θ , 3 we expect the plateau value of χ r to give the topological susceptibility.
We show the obtained values in Fig. 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 · a 2 + χ to remove the discretisation effect and extract the constant χ . For PBC, we fit from the boundary up to r max = 2 √ 8t 0 . In the case of OBC, we excluded the data in the boundary zone. Correspondingly, we used values in the 3 We recall that Euclidean SU (3) gauge theory at non-zero θ can be described by the Lagrangian L = L Fig. 10 Extraction of the pseudo-scalar screening mass from the boundary pollution at T = 2.7T c . We observe a milder temperature dependence than in the scalar sector, see Fig. 9 range [1.3 √ 8t 0 , 2 √ 8t 0 ]. It gives us measurements for the topological susceptibility which are in good agreement with χ(1.5T c )t 2 0 = 2.25(12) · 10 −5 of [22] and χ t (1.5T c )t 2 0 ∈ [1.5 · 10 −5 , 4.4 · 10 −5 ], the global fit of reference [25]. They are also consistent with the fixed lattice spacing results of [26].
At 2.7T c the situation is more intricate, as it is not clear that χ r saturates to a plateau. The trouble comes from two reasons. First, as it was already discussed in [27] in the context of a toy model with OBC, χ 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.  Fig. 11 Summary plot for the screening masses in (MeV). Only statistical errors are shown. Firstly, as expected, the scalar screening mass is the lightest of the two states. Then, the behaviour of the masses is consistent with the behaviour of the boundary zone. At T = 1.5T c , the scalar screening mass is lighter than the T = 0 lightest glueball. At T = 2.7T c it is heavier. The behaviour of the pseudo-scalar mass is also consistent; from a large mass gap between the two channels at T = 1.5T c , we move to an almost degeneracy at T = 2.7T c , which is a signal of dimensional reduction. On this figure, we also show the fixed lattice spacing results of [23]. The 15% discrepancy can most likely be attributed to systematic uncertainties (fixed lattice spacings, finite volume effects and conversion to physical units), even though a systematic difference between our methods cannot be excluded First, in Figs. 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 Figs. 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 χ 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 [3], 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 [22,28], 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 "unfrozen" 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 χ(2.7T c ) can be extracted from Fig. 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 [21] 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 Fig. 16, we show χ r and χ 2 pt (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  16 Topological charge density square χ r versus topological charge integrated two-point function χ 2 pt . Both methods are consistent in their determination of the topological susceptibility. In this figure, χ 2 pt was averaged over 24 source points associated to χ 2 pt 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.

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 Sect. 3, we investigated the typical propagation length of these effects and compared it to the zero temperature results of [8]. 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 Sect. 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 ∈ [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 [29]. 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 pseudoscalar 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 [22]. 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 [22]. Even so, it shows again that even with rather large volumes, the determination of the topological susceptibility is delicate. This supports the recent efforts [22,30], which have been undertaken to reassess the robustness of hightemperature 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.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The data will be shared upon request to one of the authors.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .