Lattice QCD without topology barriers

As the continuum limit is approached, lattice QCD simulations tend to get trapped in the topological charge sectors of field space and may consequently give biased results in practice. We propose to bypass this problem by imposing open (Neumann) boundary conditions on the gauge field in the time direction. The topological charge can then flow in and out of the lattice, while many properties of the theory (the hadron spectrum, for example) are not affected. Extensive simulations of the SU(3) gauge theory, using the HMC and the closely related SMD algorithm, confirm the absence of topology barriers if these boundary conditions are chosen. Moreover, the calculated autocorrelation times are found to scale approximately like the square of the inverse lattice spacing, thus supporting the conjecture that the HMC algorithm is in the universality class of the Langevin equation.


Introduction
One of the central goals in numerical lattice QCD is the computation of the properties of the light mesons and baryons with controlled errors. While the most important systematic errors in these calculations (finite volume and lattice spacing effects) are theoretically well understood, the relevant time scales in QCD simulations remain unpredictable. In practice, the correctness of the simulations within the quoted statistical errors can therefore only be established through empirical tests and thus only to a limited level of confidence.
In order to preserve the translation symmetry, the lattice theory is usually set up with periodic boundary conditions in all space-time directions. A side-effect of this choice of boundary conditions is the emergence of disconnected topological sectors in the continuum limit. On the lattice the sectors are not strictly separated from each other, but the relative weight in the functional integral of the gauge fields "between the sectors" decreases with a high power of the lattice spacing [1]. As a consequence, transitions from one sector to another tend to be suppressed in the simulations and may eventually become so rare that a proper sampling of the sectors would require far longer runs than are practically feasible [2][3][4].
In this paper we address both issues, the very long autocorrelation times caused by the emergence of the topological sectors and the lack of theoretical control over the simulations. The first of them we propose to avoid by imposing open boundary conditions on the gauge field in the time direction (see sect. 2). With this choice, the field space in the continuum theory becomes connected and the topological charge can smoothly flow in and out of space-time through its boundaries. All statistically relevant parts of field space are therefore expected be accessible to the simulation algorithms without having to cross higher and higher topology barriers as the lattice spacing is reduced.
When properly renormalized, some algorithms may even converge to a well-defined stochastic process in the continuum limit (sect. 3). In asymptotically free theories, such algorithms have a predictable scaling behaviour as a function of the lattice spacing and are thus theoretically controlled to some extent. The HMC algorithm [5] recently turned out to be non-renormalizable in perturbation theory and is therefore not of this kind [6], but the algorithm may conceivably fall in the universality class of the Langevin equation (whose renormalizablity was established long ago [7,8]). The empirical tests reported in sections 4 and 5 partly serve to verify that the topology barriers are indeed absent if open boundary conditions are imposed and partly to find out whether the HMC algorithm scales like an algorithm that integrates the Langevin equation.

QCD with open boundary conditions
Open boundary conditions are easily imposed in QCD and do not give rise to important theoretical complications. While the discussion in this section is more generally valid, the gauge group is taken to be SU(3) from the beginning and we assume there is a multiplet of quarks in the fundamental representation of SU (3). Our notational conventions are summarized in appendix A.

Boundary conditions in the continuum theory
In the continuum limit, the gauge and quark fields live on a four-dimensional spacetime M with Euclidean metric, time extent T and spatial size L × L × L. Time thus runs from 0 to T , while space is taken to be a three-dimensional torus, i.e. all fields are required to satisfy periodic boundary conditions in the space directions. At time 0 and T , the boundary conditions imposed on the gauge potential A µ (x) are F 0k (x)| x 0 =0 = F 0k (x)| x 0 =T = 0 for all k = 1, 2, 3, (2.1) where F µν (x) denotes the gauge-field tensor. Note that these conditions preserve the gauge symmetry and therefore do not constrain the gauge degrees of freedom of the field. In perturbation theory, the boundary conditions on the latter instead derive from the gauge-fixing procedure. If the usual Lorentz-covariant gauge is chosen, for example, the time and space components of the gauge potential are found to satisfy Dirichlet and Neumann boundary conditions, respectively †.
In the case of the quark and antiquark fields ψ(x) and ψ(x), we require that These boundary conditions are familiar from the discussion of the QCD Schrödinger functional [9,10]. Many of the theoretical results obtained in that context can actually be reused here. In particular, as explained in ref. [11], one is practically forced to choose the boundary conditions (2.2),(2.3) if parity and the time reflection symmetry are to be preserved. The action of the theory (without gauge fixing terms) is then given as usual by where g 0 and M 0 are the bare coupling and quark mass matrix.

Topology of the classical field space
Since M is contractible to a three-dimensional torus, all SU(3) principal bundles over M are trivializable. Smooth classical gauge potentials may therefore be assumed to † As is the case with periodic boundary conditions, perturbation theory in finite volume is complicated by the presence of non-trivial gauge-field configurations with vanishing action (the constant Abelian fields). The remarks made here and below on perturbation theory refer to the situation at L = ∞ and finite T , where the minimum of the action is unique up to gauge transformations.
be globally defined differentiable fields. In view of the non-linearity of the boundary conditions (2.1), the classical field space is however not a linear space. We now show that the field space is connected. Evidently, any given gauge potential A µ (x) satisfying A 0 (x) = ∂ 0 A k (x) = 0 at x 0 = 0 and x 0 = T can be smoothly contracted to zero, without violating the boundary conditions, by multiplication with a scale factor. These fields are therefore continuously connected to the classical vacuum configuration. On the other hand, if one starts from an arbitrary field A µ (x) in the classical field space, a smooth curve of gauge transformations Λ s (x), 0 ≤ s ≤ 1, may be defined through the differential equation When applied to the potential A µ (x), the transformation generates a curve in field space (parametrized by s) along which the field is continuously deformed to another field at s = 1 with vanishing time component. The transformed field can then be contracted to zero, as explained above, which proves that the field space is connected. The absence of disconnected topological sectors goes along with the fact that the topological charge is not quantized. When an instanton on M is contracted to the vacuum configuration, for example, the charge flows away through the boundaries and Q smoothly varies from 1 to 0. It may be worth noting here that the massless Dirac operator has no zero modes in the space of complex quark fields satisfying the boundary conditions (2.2). Moreover, the eigenvalues λ of the Hermitian Dirac operator γ 5 (γ µ D µ + m) are all in the range |λ| > m (see ref. [11], sect. 2.2, for a proof of these statements).
As long as the quark masses are non-negative, the quark determinant has therefore a definite sign and never passes through zero even if some masses vanish.

Renormalization and stability of the boundary conditions
The renormalization of quantum field theories on space-time manifolds with boundaries in general requires the usual (bulk) counterterms to be added to the action as well as further counterterms that are localized at the boundaries [12]. In the present case, the symmetries of the theory, power counting and the fact that ψψ vanishes at the boundaries however exclude such boundary counterterms. The renormalization of the theory thus proceeds as in infinite volume by renormalizing the coupling, the quark masses and the fields in the correlation functions considered.
Boundary conditions are subject to renormalization too and sometimes require a fine-tuning of boundary counterterms. Neumann boundary conditions in scalar field theories, for example, are known to be unstable under quantum fluctuations [12]. The situation in QCD is safe from this point of view, because there are no relevant or marginal boundary counterterms with the required symmetries. In particular, the boundary conditions (2.1)-(2.3) are stable under quantum fluctuations (see ref. [11], sect. 3, for a broader discussion of the subject).

Lattice formulation
The lattice theory is set up on a hypercubic lattice with spacing a, time-like extent T + a and spatial size L × L × L, where T and L are integer multiples of a. Periodic boundary conditions are imposed on all fields in the space directions, while time runs from 0 to T inclusive, the terminal time-slices being the boundaries of the lattice.
As usual the gauge and quark fields reside on the links and points of the lattice. In particular, the link variables U (x, µ) ∈ SU(3) live on all links (x, µ) with both endpoints in the specified range of time. The Wilson gauge action is then given by the sum over all oriented plaquettes p on the lattice, U (p) being the ordered product of the link variables around p. Only those plaquettes are included in the sum whose corners are in the time interval [0, T ]. The weight w(p) is equal to 1 except for the spatial plaquettes at time 0 and T , which have weight 1 2 . In the case of the fermion fields, a possible choice of the action is where denotes the Wilson-Dirac operator and the fields are assumed to satisfy the boundary conditions (2.2),(2.3). Since the action (2.8) depends on the quark fields at times 0 < x 0 < T only, one may then just as well set all components of the fields at time 0 and T to zero. The dynamical components of the quark fields are thus those residing in the interior of the lattice.
The functional integral and the basic correlation functions are now defined in the standard manner. Evidently, only the dynamical components of the fermion fields are integrated over in the functional integral. Note that the quark determinant is a product of real factors, one for each quark flavour, since the Wilson-Dirac operator is γ 5 -hermitian with the chosen boundary conditions. The established QCD simulation algorithms can therefore be applied straightforwardly.
It may not be completely obvious at this point that the fields satisfy the boundary conditions (2.1)-(2.3) in the continuum limit. As already mentioned in subsect. 2.3, these boundary conditions are stable under quantum fluctuations, i.e. it suffices to check that they emerge at tree-level of perturbation theory when the lattice spacing is sent to zero. The explicit expression for the quark propagator obtained in ref. [13] and a similar computation of the gluon propagator in the standard covariant gauge actually show this to be so.

Quantum-mechanical representation
The formulation of the lattice theory described above admits a quantum mechanical description in terms of a Hilbert space H of physical states and a bounded, positivedefinite transfer matrix T [14]. In particular, the partition function of the theory is given by the expectation value of a power of the transfer matrix in a state Ω ∈ H that encodes the chosen boundary conditions at time 0 and T . A relatively direct way of introducing the transfer matrix formalism starts from a representation of the physical states through wave functions depending on a gauge field V ( x, k) and the components of a quark field on the spatial lattice (see ref. [10], for example). The boundary state Ω is represented by the wave function in this language, the covariant derivatives being evaluated in presence of the gauge field V . Note that this expression is manifestly invariant under the gauge symmetry, the lattice symmetries and the (vector-like) flavour transformations. In other words, Ω has the quantum numbers of the vacuum state.
Correlation functions of gauge-invariant fields have a quantum mechanical interpretation as well. The two-point function of a local scalar field φ(x), for example, is given by whereφ( x) denotes the operator field associated to φ(x) [14]. Since the transfer matrix and the space of physical states are independent of the boundary conditions at time 0 and T , the hadron masses and many other physical quantities can, in principle, be extracted from such correlation functions in much the same way as on lattices with periodic boundary conditions in time.

On-shell O(a) improvement
The O(a) improvement of the lattice theory follows the lines of refs. [15,16]. There is actually very little difference with respect to the case of the Schrödinger functional discussed in the second of these papers. In particular, all bulk O(a) counterterms and their coefficients are exactly the same as those required for the on-shell improvement of the theory on the infinite lattice. We wish to emphasize at this point that a further improvement is not needed if one is exclusively interested in the correlation functions of fields localized far away from the boundaries of the lattice, where the effects of the latter are exponentially suppressed. The improvement of correlation functions involving fields close to or at the boundaries however requires the addition of the O(a) boundary counterterms In these equations, p s runs over all space-like oriented plaquettes at the boundaries of the lattice and the coefficients must be adjusted so as to cancel the boundary effects of order a (since the boundary conditions are not the same, there is no reason to expect these coefficients to coincide with those needed to improve the Schrödinger functional).

Other lattice formulations of the theory
Lattice QCD with open boundary conditions can be set up in many different ways. Universality actually suggests that the details of the lattice theory become irrelevant in the continuum limit if the gauge, space-time and flavour symmetries are respected. Lattice formulations that preserve chiral symmetry away from the boundaries exist as well, but some care is required in this case in order to guarantee the locality of the lattice Dirac operator near the boundaries [17,11].

Dynamical properties of QCD simulations
The interpretation of simulation data requires good control over the simulation dynamics. In this section, the relevant notions are briefly discussed and some specific issues are addressed, which arise when studying the scaling behaviour of QCD simulations.

Autocorrelations
QCD simulation algorithms produce random sequences of gauge-field configurations recursively, where the next configuration is obtained from the current one according to some transition probability. The simulation algorithms considered in this paper are the HMC algorithm [5] and the closely related SMD (stochastic molecular dynamics, or generalized HMC) algorithm [18]. In both cases the simulation time is proportional to the molecular-dynamics time in lattice units and will, for simplicity, be identified with the latter in the following.
Let O i be a set of real-valued unbiased observables labeled by an index i. Their values O i (t) measured at simulation time t are statistically correlated to some extent, i.e. the connected parts of the n-point autocorrelation functions in general do not vanish. In this equation, the bracket . . . stands for the average over infinitely many statistically independent parallel simulations, which is the same as the average over time translations if the simulation is ergodic. The connected parts of the autocorrelation functions tend to fall off exponentially at large separations in simulation time. In particular, the two-point autocorrelation can be shown to have a spectral decomposition of the form † where τ 0 ≥ τ 1 ≥ . . . are the so-called exponential autocorrelation times of the algorithm. While these are independent of the observables considered, the coefficients c in measure how strongly the observables O i couple to the eigenmode number n of the transition probability. Note that neither the spectral values λ n nor the coefficients c in are guaranteed to be real, except in the case of the HMC algorithm and the Langevin limit of the SMD algorithm.
In practice, the integrated autocorrelation times of the observables of interest play an important rôle, where ∆t is the separation in simulation time of the observable measurements. The sum in eq. (3.4) amounts to a numerical integration of the normalized autocorrelation function ρ i (t) using the trapezoidal rule. In particular, in the Langevin limit of the SMD algorithm or if the HMC algorithm is used, the formula and thus the bound τ int (O i ) ≤ τ 0 hold up to integration errors.

Topology-changing transitions
On lattices with periodic boundary conditions, the probability per unit simulation time for a HMC or an SMD trajectory to pass from one topological sector to another is a rapidly decreasing function of the lattice spacing [2][3][4]. Such topology-tunneling † The HMC and the SMD algorithm both evolve the gauge field U (x, µ) together with its canonical momentum π(x, µ) (cf. subsect. 4.1). Equation transitions are non-perturbative lattice artifacts that may informally be described as "an instanton falling through the lattice". The integrated autocorrelation time of the topological charge consequently tends to become very large, sometimes to the extent that the correctness of the simulation is compromised.
With open boundary conditions, the situation is different, because the topological charge can change smoothly along a molecular-dynamics trajectory by flowing in and out of the lattice via its boundaries. A catastrophic slowdown of the algorithms as in the case of periodic boundary conditions is therefore not expected.

Renormalizable algorithms
The n-point autocorrelation functions of gauge-invariant local fields formally look like the correlation functions in a field theory in five dimensions, where the simulation time is the fifth space-time coordinate. When the lattice spacing is taken to zero, the autocorrelation functions may then conceivably have a continuum limit, provided the fields and the parameters (of both the theory and the simulation algorithm) are properly renormalized.
Algorithms that integrate the Langevin equation are known to be renormalizable in this sense to all orders of perturbation theory [7,8]. An example of an algorithm of this kind is provided by the SMD γ algorithm (cf. subsect. 4.2). The simulation time has physical dimension [length] 2 in this case and must be renormalized according to where t is the simulation time in lattice units, Z t (g 0 ) a renormalization constant and t R the renormalized simulation time in some physical units. Further renormalization is not required apart from the usual field and parameter renormalization. Beyond perturbation theory, the renormalizability of an algorithm (and thus the associated scaling laws) can break down as a result of non-perturbative lattice artifacts. On lattices with periodic boundary conditions, topology-changing transitions have this effect in the case of the SMD γ algorithm. However, if open boundary conditions are chosen, there is currently no reason to expect that the renormalizability of the algorithm does not extend to the non-perturbative level.

Scaling behaviour of the HMC algorithm
Free-field studies of the HMC algorithm suggest that the exponential autocorrelation times scale linearly (like a −1 ) if the length of the molecular-dynamics trajectories is scaled accordingly [19]. The algorithm however turns out to be non-renormalizable in perturbation theory [6] and its scaling behaviour in the presence of interactions may therefore be completely different.
The empirical studies reported later actually show that the HMC algorithm (on lattices with open boundary conditions) appears to fall into the universality class of the Langevin equation. In particular, the autocorrelation times scale approximately like a −2 rather than linearly. From this point of view, the non-renormalizability of the HMC algorithm in perturbation theory merely reflects the fact that the leadingorder theory is in the wrong dynamical universality class and therefore not a suitable starting point for the perturbation expansion.

Making QCD simulations safer
In practice, numerical simulations should be much longer (by, say, a factor 100 at least) than the longest exponential autocorrelation time τ 0 , as otherwise a proper sampling of the functional integral is not guaranteed and the simulation may consequently be biased in an unpredictable way. Usually the integrated autocorrelation times of the quantities of interest are monitored, but it should be noted that the correctness of the simulation results (within the estimated statistical errors) cannot be taken for granted if only these autocorrelation times are much smaller than the total simulation time.
Integrated autocorrelation times of both physical and other observables can in fact be very much smaller than τ 0 . In particular, the autocorrelation times of noisy quantities (large Wilson loops, for example) tend to be practically unrelated to the exponential autocorrelation times. To illustrate this point, consider an observable where c is a constant and η a statistically independent Gaussian noise with mean zero and unit variance. O 0 has the same expectation value as O 1 and its autocorrelation function is given by At large c, i.e. when the added noise term is large, the integrated autocorrelation time τ int (O 0 ) decreases like 1/c 2 and can therefore be made arbitrarily small. Nothing is gained in this way, but the example shows that integrated autocorrelation times may not be representative of the true autocorrelations in the simulation. Exponential autocorrelation times are difficult to determine reliably if very long simulations are impractical. In this case, a pragmatic way to proceed is to look for observables with large integrated autocorrelation times and to take the maximum of the latter as an estimate of τ 0 . The observables that provide the best probes for autocorrelations should be sensitive to the smooth modes of the gauge field, since these tend to be updated least efficiently. Moreover, for the reasons given above, good probes are likely to have small statistical fluctuations. Quantities obtained by integrating the Wilson flow [1], such as the average action density E at positive flow time, satisfy both criteria and are therefore recommended probes (see fig. 1).

Numerical studies
In order to verify and complement the theoretical discussion in the previous sections, we performed extensive simulations of the SU (3) gauge theory with open boundary conditions. The algorithms, observables and simulation parameters used in these studies are described in this section.

Simulation algorithms
Both the HMC and the SMD algorithm operate in phase space, i.e. on the gauge field U (x, µ) and its canonical su(3)-valued momentum field π(x, µ). The O(a) boundary counterterm (2.14) is not included in the Hamilton function H(π, U ) = 1 2 (π, π) + S G (U ) (4.1) of the system, partly for simplicity and partly because the term is unimportant in the present context. The HMC algorithm proceeds in cycles, where in each cycle one first chooses the momentum field randomly, with normal distribution, and then evolves the fields according to the molecular-dynamics equations that derive from the Hamilton function (4.1). In our simulations, the equations were integrated from molecular-dynamics time 0 to τ using n 0 iterations of the 4th-order Omelyan-Mryglod-Folk (OMF) integrator defined through eqs. (63) and (71) in ref. [20]. At the end of the evolution, the fields are submitted to an acceptance-rejection step that corrects for the integration errors. This algorithm has two parameters, τ and n 0 , and requires the derivative of the gauge action to be calculated 5n 0 times per cycle.
In the case of the SMD algorithm, one proceeds in essentially the same way, but the momentum field is only partially refreshed according to π(x, µ) → c 1 π(x, µ) + c 2 υ(x, µ), (4.2) where υ(x, µ) is a randomly chosen momentum field with normal distribution, while γ and δτ are parameters of the algorithm. The molecular-dynamics equations are then integrated from 0 to δτ by applying a single iteration of the 4th-order OMF integrator and the fields are finally submitted to an acceptance-rejection step. When rejected, the fields are reset to their values before the integration, except for a change of sign of the momentum field (see ref. [21] for a straightforward proof of the correctness of the algorithm). Note that the simulation time t elapsed after n SMD cycles is, by definition, equal to nδτ , irrespectively of the rejection rate.
Since the OMF integrator is applied only once, the molecular-dynamics evolution time δτ is usually set to a value much smaller than 1 in order to guarantee a high acceptance rate P acc . Otherwise the SMD algorithm is frequently backtracking, on average after every period of time equal to t acc = δτ P acc 1 − P acc , (4.5) and thus tends to become inefficient. With respect to the leapfrog and the 2nd-order OMF integrator, the 4th-order OMF integrator has the advantage that very high acceptance rates can be achieved with a moderate computational effort.

Stochastic equation, parameter scaling and the SMD γ algorithm
In the limit δτ → 0, the SMD algorithm amounts to solving the stochastic moleculardynamics equations where η s a Gaussian random noise with mean zero and variance In these equations, the evolution time s and the mass µ 0 are related to the simulation time t and the parameter γ through s = ta, µ 0 = γ/2a, (4.9) respectively. Evidently, eqs. (4.6),(4.7) reduce to the standard molecular-dynamics equations if µ 0 is set to zero (see appendix A for the definition of the derivative of the gauge action). When the continuum limit is approached, the scaling behaviour of the simulation algorithms depends on how their parameters are scaled. The fact that the evolution time in eqs. (4.6),(4.7) has dimension [length] suggests to scale the HMC trajectory length τ proportionally to 1/a [19]. For the same reason, one can argue that µ 0 should be scaled like a physical mass up to a logarithmically varying renormalization factor perhaps. This choice of the parameter scaling (which, however, leads to nonremovable ultra-violet singularities in perturbation theory [6]) will be referred to as free-field scaling.
Alternatively, if γ is held fixed, and if δτ is such that the continuous-evolution time t acc is on the order of the exponential autocorrelation times (or larger), the SMD algorithm effectively performs a numerical integration of the Langevin equation [6]. For clarity, we use the acronym SMD γ for the SMD algorithm with this parameter scaling.

Observables
As already noted in subsect. 3.5, observables based on the Wilson flow probe the slow modes of the gauge field and are therefore well suited for studying autocorrelations in QCD simulations. A review of the Wilson flow is beyond the scope of this paper and we merely write down the differential equation that generates the flow V t (x, µ), t ≥ 0, in the space of gauge fields (see refs. [1,4,22] for a comprehensive discussion of the flow and some of its surprising properties).
In the course of the simulations, the observables are evaluated at fixed separations in simulation time. Starting from the current gauge-field configuration U (x, µ), we first integrate the flow equation (4.10) numerically up to some flow time t. The field tensor G µν (x) of the gauge field V t (x, µ) generated in this way is defined through the clover formula, i.e. through the four plaquette Wilson loops in the (µ, ν)-plane that start and end at x (at the boundaries x 0 = 0 and x 0 = T we set G 0k (x) = 0). The primary observables considered are then the time-slice averages of the action density and the time-slice sums of the topological charge density. Evidently, the autocorrelations of the total charge are studied as well. In all these equations, the dependence on the flow time has been suppressed for simplicity. At positive flow time t, the expectation values of arbitrary (finite) products of the observables E(x 0 ), Q(x 0 ) and Q do not require renormalization and are expected to  (3) 17.49(4) * Calculated at physical time L/2 on the (L/a) 4 lattices have a well-defined limit when the lattice spacing is taken to zero [1,22]. While these expectation values do not have any obvious interpretation in terms of glueballs or colour flux tubes, for example, they are properties of the continuum theory which reflect the dynamics of the smooth modes of the gauge field (the smoothing radius being roughly equal to √ 8t). In particular, as explained in ref. [1], on lattices with periodic boundary conditions, the topological charge Q (as defined here) converges to an integer-valued observable in the continuum limit, which labels the topological sectors of field space.

Lattice and simulation parameters
In table 1 we list the spatial sizes and the inverse gauge couplings β = 6/g 2 0 of the lattices that we have simulated. The number of lattice points in the time direction (which is equal to T /a + 1) coincides with L/a in most cases, but lattices with larger time extent have been considered too. For the conversion to physical units, we use the Sommer radius r 0 = 0.5 fm [23] and the results obtained for r 0 /a by Necco and Sommer [24]. The values of the lattice spacing determined in this way (3rd column of table 1) are such that all lattices have the same physical size L, as is desirable for a scaling study.
As a reference for the Wilson flow time t, we prefer to use the scale t 0 determined through the implicit equation [1] t 2 E(L/2) t=t 0 = 0.3. (4.14) At flow time t 0 , the Wilson flow has a smoothing range approximately equal to r 0 , i.e. this point in flow time is about where the non-perturbative regime sets in. Since L is quite small in physical units, the values of t 0 /a 2 quoted in table 1 are probably affected by finite-volume effects and they are, in fact, a few percent lower than those  previously obtained in ref. [1] at L ≃ 2.4 fm. In the present context, the effect can however be safely ignored since L is the same on the lattices simulated. The trajectory length τ in the HMC simulations was scaled according to the freefield parameter scaling, and the number n 0 of integration steps (each consisting of one iteration of the 4th-order OMF integrator) was then tuned to achieve fairly high acceptance rates P acc (see table 2). In a second set of simulations, we used the SMD γ algorithm. Some experimenting suggests that the autocorrelation times of the observables considered have a flat minimum near γ = 0.3 and we therefore decided to stick to this value of γ. The other parameter of the algorithm, δτ , was adjusted to ensure acceptance over an average simulation time t acc significantly larger than the exponential autocorrelation times (see table 3).
The observables were measured at the separations ∆t in simulation time quoted in tables 2 and 3. On each lattice, a fairly large number N cnfg of configurations were analyzed, the total length of the simulations thus being equal to N cnfg ∆t.

Scaling properties of the autocorrelation functions
While the chosen observables do not require renormalization, the flow time at which they are evaluated should be scaled like a physical quantity of dimension [length] 2 in the continuum limit. In the following, the flow time is set to the reference time t 0 , the results at other values of the flow time being similar as long as the short-time regime (where lattice effects are large) is avoided.
The SMD 0.3 algorithm is renormalizable to all orders of perturbation theory since it effectively integrates the Langevin equation (cf. subsects. 3.3 and 4.2). Moreover, with open boundary conditions, the topology barriers that otherwise slow down the algorithm are absent. It is therefore not unreasonable to expect that the normalized autocorrelation functions of the selected observables converge to universal functions in the continuum limit, provided the simulation time is scaled according to eq. (3.6).
The autocorrelation functions plotted in fig. 2 in fact behave as expected if one assumes that the renormalization constant Z t varies only slightly on the lattices considered. Note that all points obtained on a given lattice are statistically correlated. In particular, the seemingly systematic deviation of the measured autocorrelation functions on the 40 4 lattice from those on the 32 4 and 24 4 lattices may very well be a statistical fluctuation. Large deviations are however seen in the case of the time-slice and the total topological charge on the coarser lattices, where topology-tunneling transitions are not totally suppressed and thus reduce the autocorrelations. Langevin scaling then sets in as expected once these lattice artifacts become unimportant.
On physically large lattices, the four-point autocorrelation function of the topological charge Q is dominated by its disconnected parts. The normalized two-point autocorrelation function of Q 2 is then related to the one of Q through Although the simulated lattices are not very large in physical units, we found that eq. (5.1) is accurately satisfied. In particular, the autocorrelation functions of Q and Q 2 scale in practically the same way.

Autocorrelation times
Similarly to the energy spectrum in finite volume, the exponential autocorrelation times depend on the symmetry sector considered. In particular, eq. (5.1) suggests that the longest autocorrelation time in the odd-parity sector is larger, by a factor 2 perhaps, than the one in the even-parity sector. On the basis of the data shown in fig. 2, we estimate that the latter is about 1.2×(r 0 /a) 2 (thus ranging from 30 to 187) in the case of the SMD 0.  usual, such estimates should be taken with a grain of salt, because the slowest modes in the system may not couple sufficiently strongly to the measured observables for their effects to be seen in the available data.
The integrated autocorrelation times plotted in fig. 3 and their errors were calculated following the lines of appendix B. As is evident from the figure, the autocorrelation times all scale linearly in 1/a 2 and thus as expected for algorithms that integrate the Langevin equation. From the point of view of the continuum limit, the intercepts at L/a = 0 of the straight lines in fig. 3 are O(a 2 ) lattice corrections to the Langevin scaling, while the ratios of their slopes are universal properties of the simulation dynamics. Figure 3 also shows that the HMC algorithm (with free-field parameter scaling) scales like the SMD 0.3 algorithm. The matching of the autocorrelation times requires a renormalization of the simulation time by the factor Z ≃ 1.32, but in terms of computer time, HMC simulations tend to be faster than SMD 0.3 simulations, because a very accurate integration of the molecular-dynamics equations is not needed.

Dependence on the time-like extent of the lattice
In practice, the time extent of the simulated lattices will often have to be larger than the one of the (L/a) 4 lattices in our scaling studies. Autocorrelations in general depend on the physical situation and thus also on the lattice geometry. For illustration, the autocorrelation times of E and Q 2 calculated on three lattices with the same spacing and spatial size, but different time extent T , are plotted in fig. 4.
Close to boundaries of the lattice, the autocorrelation times shown in these plots are thus practically independent of T , while well inside the lattices they rapidly converge to a constant value when T is increased. The behaviour of the autocorrelation times of these observables discussed in subsection 5.2 is therefore expected to be representative of the situation on larger lattices as well.
The total topological charge Q is a special case, because it can only change (at small lattice spacings) by flowing in and out of the lattice. In the course of a simulation, the measured values of Q fluctuate around the origin with a standard deviation that increases proportionally to √ T L 3 on large lattices. The charge however flows through the boundaries with a rate proportional to √ L 3 only. The simulation time required for a significant change in Q must therefore be expected to grow with T (proportionally to T if Q performs a random walk). On the 24 4 , 48×24 3 and 80×24 3 lattices, we actually find that the autocorrelation times of Q 2 (42.1(2.5), 113(6) and 148(10), respectively) grow roughly linearly with T .
We wish to conclude this discussion by emphasizing that the autocorrelation times on lattices of a given physical size are expected to scale linearly in 1/a 2 . Independently of the chosen geometry, the computational effort for HMC simulations with the standard leapfrog integrator, for example, thus scales approximately like 1/a 7 .

Conclusions
The theoretical and empirical results presented in this paper show that the topology barriers in the SU(3) gauge theory can be avoided by choosing open boundary conditions in the time direction. Moreover, on lattices with these boundary conditions, the HMC and the SMD γ simulation algorithm both appear to fall in the dynamical universality class of the Langevin equation, i.e. simulations based on these algorithms slow down proportionally the square of the lattice spacing when the continuum limit is approached.
In our numerical studies, the autocorrelation times of the topological charge (as well as those of observables unrelated to the latter) went up to values greater than 100 in units of molecular-dynamics time. While such autocorrelations may be affordable in a given case, the experience suggests that there is ample room for algorithmic improvements. A separate treatment of the high-frequency and the smooth modes of the gauge field, for example, might be worth considering at this point.
Open boundary conditions can easily be imposed in QCD with a non-zero number of sea quarks. We do not foresee any technical issues when these boundary conditions are chosen, but an interesting theoretical question is whether the Langevin equation remains renormalizable in the presence of the pseudo-fermion fields that need to be introduced to be able to simulate the theory [27,28].
All simulations reported in this paper were performed on a dedicated PC cluster at CERN. We are grateful to the CERN management for funding this machine and to the CERN IT Department for technical support.

Appendix A. Notational conventions
The Lie algebra su(3) of SU(3) may be identified with the linear space of all traceless anti-hermitian 3 × 3 matrices. We choose the generators T a , a = 1, . . . , 8, of the Lie algebra to be such that The general element X of su(3) is then given by X = X a T a with real components X a (repeated indices are automatically summed over). The Euclidean Dirac matrices γ µ , µ = 0, . . . , 3, are assumed to be hermitian. Gauge potentials A µ (x) take values in su (3) and are normalized such that the field tensor and the covariant derivatives that appear in the Dirac operator are given by On the lattice, the gauge-covariant forward and backward difference operators in presence of a lattice gauge field U (x, µ) act on the quark field ψ(x) according to where a denotes the lattice spacing andμ the unit vector in direction µ.
The scalar product of any two vector fields ω(x, µ) and υ(x, µ) with values in su(3) is normalized such that (ω, υ) = −2a 4 x,µ tr{ω(x, µ)υ(x, µ)}. (A.6) If F(U ) is a differentiable function of the gauge field, its derivative with respect to the link variable U (x, µ) in the direction of the generator T a is defined by , U t (y, ν) = e tT a U (x, µ) if (y, ν) = (x, µ), U (y, ν) otherwise. (A.7) In particular, in the case of a scalar function F(U ), the combination T a ∂ a x,µ F(U ) is a vector field with values in su(3) that transforms under the adjoint representation of the gauge group.
On the (L/a) 4 lattices considered, the summation window for even-parity observables was set to Given the measured exponential autocorrelation times (subsect. 5.2), the systematic error that derives from the truncation of the sum (B.2) is estimated to be at most 3% with this choice. The statistical errors of the autocorrelation functions and the integrated autocorrelation times were determined using the Madras-Sokal approximation [25] (see ref. [26], appendix E, for a detailed description of the procedure).