Simulation of QCD with N_f=2+1 flavors of non-perturbatively improved Wilson fermions

We describe a new set of gauge configurations generated within the CLS effort. These ensembles have N_f=2+1 flavors of non-perturbatively improved Wilson fermions in the sea with the Luescher-Weisz action used for the gluons. Open boundary conditions in time are used to address the problem of topological freezing at small lattice spacings and twisted-mass reweighting for improved stability of the simulations. We give the bare parameters at which the ensembles have been generated and how these parameters have been chosen. Details of the algorithmic setup and its performance are presented as well as measurements of the pion and kaon masses alongside the scale parameter t_0.


Introduction
What can be achieved in lattice QCD computations in terms of observables and their accuracy depends to a large extent on the availability of suitable ensembles of gauge field configurations. For reliable results, many sources of error have to be controlled: Fine lattices are needed for minimal discretization effects, the quark masses have to be close to their physical values and the volume of the lattices has to be large enough for finite size effects to be small. The final precision also depends on the quark flavor content of the sea. On top of these systematic effects come statistical uncertainties: Simulations have to be long enough such that the statistical errors can be estimated reliably, and with statistical uncertainties getting smaller, the need for control over systematic effects increases.
Since the generation of gauge field configurations is computationally the most demanding part of the whole computation, a careful evaluation of the physics parameters needs to be made in view of the target precision of the observables -as far as this is possible at this stage. The goal is to balance the various sources of systematic and statistical uncertainties in the final result: In light of the findings of Ref. [1], for example, we do not include a dynamical charm quark in the sea as we do not anticipate to be able to reach an accuracy comparable to its effect on typical low-energy observables after taking the continuum limit and the chiral extrapolation. On the contrary, including the charm might introduce large lattice artifacts and would make the tuning procedure more difficult.
Recent year's advances have led to a re-evaluation of the requirements for a reliable lattice computation regarding the control over statistical errors. Notably, it has been known for a while that the global topological charge freezes on fine lattices with periodic boundary conditions [2][3][4]. However, with the advent of the gradient flow in lattice computations [5,6] it has been discovered that at moderate lattice spacing other quantities constructed from smoothed fields evolve even slower in Monte Carlo time [7,8]. To exclude uncontrolled biases in any observable, Monte Carlo histories much longer than the exponential autocorrelation time are required, i.e. much longer than the times observed in these smoothed observables and therefore longer than previously thought.
In this paper, we give an overview of the first round of the CLS (Coordinated Lattice Simulations) effort to generate configurations with N f = 2 + 1 flavors of nonperturbatively improved Wilson fermions. In some of its aspects it is a continuation of the N f = 2 flavor project: We use a non-perturbatively improved Wilson fermion action, we do not employ link smearing, simulations are done using a public code and we focus on small lattice spacings for a controlled continuum limit [9].
By adding the additional flavor to the sea, one naturally aims at higher accuracy than with two flavor simulations. In order to achieve this, there are also improvements over the previous project. We use open boundary conditions in the time direction which prevent the topological charge from freezing [5,7] and twistedmass reweighting to avoid the sector formation due to zero eigenvalues of the Wilson fermions and the resulting instabilities in the simulation [10].
The simulations are performed using the openQCD code version 1.2 [11] whose general algorithmic setup is described in Ref. [12]. The code provides broad flexibility with respect to the algorithms used, starting from the determinant decomposition, the molecular dynamics integration and the methods employed for solving the Dirac equation. In this paper we give the physics and algorithmic parameters of these simulations and report on our experiences with this new setup. Furthermore we present first measurements of basic physics observables: the masses of the pion and the kaon as well as the scale parameter t 0 [6] on which we base the tuning of the runs.
Similar large-scale simulations of QCD have recently been performed by the PACS-CS simulating improved Wilson fermions [13], the QCDSF collaboration with N f = 2 + 1 flavors of NP improved, smeared Wilson fermions [14], the Hadron Spectrum collaboration using N f = 2 + 1 flavors of tree-level improved, smeared Wilson fermions on anisotropic lattices [15], as well as the ETM collaboration using twistedmass fermions [16] and the BMW collaboration with tree-level improved smeared fermions [17]. Also domain wall fermions are employed by RBC-UKQCD [18] and overlap fermions by JLQCD [19] as well as smeared rooted staggered fermions by MILC [20]. Our simulations are unique by their use of open boundary conditions and, among the simulations with standard Wilson fermions, twisted-mass reweighting as a safeguard against the effects of near-zero modes of the Dirac operator.
The paper is organized as follows: In Section 2 we give the details of the action, the tuning strategy and the parameters of the runs. The algorithmic setup is described in Section 3. Autocorrelations observed in the simulations are the subject of Section 4, while the two types of reweighting used in the light and the strange quark sector are discussed in Section 5. This is followed in Section 6 by the measurement of the pseudoscalar masses and the scale parameter t 0 and a discussion of discretization effects in Section 7.

Physical parameters
The simulations are done on lattices of size N t × N 3 s , with open boundary conditions imposed on time slice 0 and N t −1. Lattices with N t points in the temporal direction therefore have a physical time extent of T = (N t − 1)a in conventional notation, with a the lattice spacing.

Action
The general setup of the lattice actions which can be simulated with the openQCD code has already been given in detail in Ref. [12]. In particular it is described there how the boundary conditions are imposed. Therefore, here we only give details of the bulk action. Throughout, the coefficients of the boundary improvement terms are set to their tree level values.
For the gauge fields, we use the Lüscher-Weisz action [21] with tree level coefficients -which is different from the earlier reference [12] where the Iwasaki action has been employed. In the bulk, the plaquette and rectangle terms are multiplied by their respective coefficients c 0 = 5/3 and c 1 = −1/12 where the sums run over the plaquettes p and the rectangles r contained in the lattice and β = 6/g 2 0 with the bare gauge coupling g 0 . For the fermions, the Wilson Dirac operator [22] including the Sheikholeslami-Wohlert term needed for O(a) improvement of the action [23] is used with ∇ µ and ∇ * µ the covariant forward and backward derivatives, respectively. The improvement term containing the standard discretization of the field strength tensor F µν [24] comes with the coefficient c SW whose value has been determined nonperturbatively in Ref. [25].
The three flavor fermion action then reads where we take the up and down quark masses to be degenerate m 0,ud ≡ m 0,u = m 0,d . The strange-quark mass m 0,s is tuned as a function of the light quark mass. In the following, we frequently quote the hopping parameters κ f instead of the bare quark masses

Choice of parameters
Since we restrict ourselves to N f = 2 + 1 flavor QCD, which is different from the full Standard Model, the point of "physical" quark masses is not unique even in the continuum and we therefore have to fix observables which define it. For the tuning of our runs, we set the scale through t 0 defined by the Wilson flow [6], see Section 6.3. The quark masses are set using the masses of the pion and the kaon. While this choice is convenient during the tuning of the runs, it can be changed in the future once more observables are available. The lattices at different cutoff are matched via the dimensionless parameters where all quantities are the ones measured at the parameter values of the ensemble in question. Note that in leading order of Chiral Perturbation Theory (ChPT) they are proportional to the sum of the quark masses, φ 2 ∝ (m u + m d ) and φ 4 ∝ (m u +m d +m s ) [26,27]. The advantage of this strategy is that we obtain all quantities involved with high statistical accuracy from the simulated ensembles, without further need of renormalization constants or chiral extrapolation. Particular drawbacks of this strategy are the significant cutoff effects which we observe in the various definitions of t 0 /a 2 on our largest lattice spacings, as discussed in Section 7.1. Furthermore, the value of t 0 is not an experimentally accessible observable and only known from other lattice simulations. In the literature one finds √ 8t 0 = 0.4341(33) fm by the ALPHA collaboration using Wilson fermions in two-flavor QCD [28] and √ 8t 0 = 0.4144(59)(37) fm by the BMW collaboration using N f =2+1 flavors [29]. In a 2 + 1 + 1 flavor setup with rooted staggered fermions, the HPQCD collaboration finds √ 8t 0 = 0.4016(23) fm [30]. As has been observed in Ref. [28], these numbers exhibit a significant flavor content effect, which however is monotonic in the number of flavors. Since our simulation setup is also with N f = 2 + 1 flavors, we choose the value of Ref. [29] and the QCD values of m π = 134.8(3) MeV and m K = 494.2(4) MeV in the isospin limit [31] which leads to a physical point estimate where errors have been added in quadrature. From this choice and our measurements of t 0 /a 2 presented below, we estimate for our three values of β = 3.4, 3.55 and 3.7 lattice spacings a of a ≈ 0.086 fm, 0.064 fm and 0.05 fm, respectively.

Quark mass trajectory
In order to achieve O(a) improvement, the bare coupling -as all bare parameters -has to be improved with a mass-dependent term [24] with m cr the critical quark mass whose precise value is not known at this stage. To keep the lattice spacing constant as we change the sea quark masses, this modified coupling constantg 0 has to be kept constant. While the coefficient b g is small at one loop in perturbation theory [32], b g = 0.012 N f g 2 0 , a non-perturbative result is not known for any action. To keepg 0 fixed, we therefore keep the sum over the subtracted quark masses fixed, a strategy already proposed in Ref. [14]. Note that this is equivalent to keeping the sum over the bare quark masses m 0,f fixed Up to effects of order O(am ud ), this also implies a constant sum of improved PCAC quark masses [33]. We can therefore define chiral trajectories by a point in the φ 2 -φ 4 plane, at which different lattice spacings are matched and the requirement that the sum of the bare quark masses is constant. For each value of β, the lattice spacing is constant along these lines and a continuum limit can be performed for each value of φ 2 .
As explained in the next section, we match lattices with different lattice spacings at m π = m K ≈ 420 MeV, where we will also show first results concerning the size of the O(am) cutoff effects introduced by this choice.

Tuning strategy
By choosing the chiral trajectories of Eq. (2.8), the tuning process can be highly simplified: Keeping β fixed, for each chiral trajectory we match the lattices at the flavor symmetric point, i.e. where all quarks have equal masses.
Determining the slope of φ 4 as a function of φ 2 at β = 3.4 from a set of preliminary runs, not shown here, we estimate the target value on the symmetric line φ 4 m ud =ms = 1.15 . (2.9) With the final statistics, we are able to reach better than 1% accuracy in this quantity and a matching of the target value within one standard deviation. In the chiral limit this translates into an accuracy of about 1 MeV in the strange-quark mass. In the future, we plan to have more chiral trajectories which will allow us to study the consequences of the remaining mistuning. The result of the tuning effort and the resulting trajectory in the φ 2 -φ 4 plane is shown in Figure 1 with results from the ensembles given in Table 1. Within the statistical accuracy, we do not observe significant cutoff effects. The one point at β = 3.7 is still under production and its error therefore not yet trustworthy. We observe, the quark mass effect on φ 4 along this trajectory is moderate, around 5% between the chiral limit and the symmetric point, as expected from ChPT.

Algorithmic parameters
The basic algorithmic setup has already been described in detail in Ref. [12], but since we are presenting simulations with larger lattices, statistics and a different action, the various settings needed to be reconsidered. Here we give the parameters at which the runs were performed and the reasoning behind the various choices.

Twisted-mass reweighting
Since the Wilson Dirac operator is not protected against eigenvalues below the quark mass, field space is divided by surfaces of zero eigenvalues. These barriers of infinite  action cannot be crossed during the molecular dynamics evolution, at least if the equations of motion are integrated exactly. While at a sufficiently large volume and quark mass this might not be a problem in practice [34], it can lead to instabilities during the simulation and meta-stabilities in the thermalization phase. Lüscher and Palombi [10] therefore suggested to introduce a small twisted-mass term into the action during the simulation and compensate for this by reweighting.
In the present simulations, we use the second version of the reweighting suggested in Ref. [10], which is less affected by fluctuations in the reweighting factor from the ultraviolet part of the spectrum of the Dirac operator. Contrary to the original proposal, we do not apply it to the Hermitian Dirac operator Q = γ 5 D W but to the Schur complementQ = Q ee − Q eo Q −1 oo Q oe of the asymmetric even-odd preconditioning [35]. This amounts to replacing the determinant of the light quark pair by The reweighting factor which needs to be included in the measurement of primary observables then reads The choice of the parameter µ 0 will be discussed in Section 5.1.  Table 1: List of the ensembles. In the id, the letter gives the geometry, the first digit the coupling and the final two label the quark mass combination. We give rounded values of m π and m K using the t 0 /a 2 of the ensemble and √ 8t 0 = 0.4144 fm. Using t 0 /a 2 extrapolated to the physical light quark masses, we estimate lattice spacing of a ≈ 0.086 fm, 0.064 fm and 0.05 fm for β = 3.4, 3.55 and 3.7, respectively.

Determinant factorization
The fluctuations in the forces have to be reduced further than what can be achieved by introducing an infrared cutoff by the twisted mass µ 0 . To this end we use Hasenbusch's mass factorization [36] with a twisted mass [37] applied to the last term in Eq. (3.1) [38] with a tower of increasing values of µ 0 < µ 1 < · · · < µ N mf . The values of these masses can significantly influence the performance of the algorithm. Here we roughly set them at equal distances on a logarithmic scale as suggested in Ref. [12]. The precise values of the µ i are listed in Table 2, which implicitly gives also the number of factors in Eq. (3.3). The combination of twisted-mass reweighting and mass factorization leads to an effective action for the light quark pair with N mf + 2 terms  Table 2: Parameters of the algorithm: We give the twisted masses used in the reweighting and mass factorization, the N mf,2 lightest of which are integrated on the coarsest time scale, the number of poles N p and the range used in the RHMC, with N p put on single pseudofermions , N p,2 of which are integrated on the outer level. N s,2 is the number of steps of the outer level of the MD integrator used for one trajectory. The total length of the Markov chain and the acceptance rate are also given.
The single term with the largest twisted mass and the one from the diagonal determinant det Q oo are always integrated together and are therefore counted as one term.

RHMC
The strange quark is simulated using the RHMC algorithm [39,40], where the matrix square root is approximated by a rational function Zolotarev's optimal approximation to the inverse square root in the interval [r a , r b ] with a given number of poles N p determines the parameters A and {μ i ,ν i }. The strange-quark mass as argument of Q andQ has been suppressed for readability. W 1 is the reweighting factor, implicitly defined by Eq. (3.5), which has to be included in the measurement. The values used in the various runs are specified in Table 2.
The openQCD code gives the option to split the determinant of the rational function in Eq. (3.5) into several factors. In our simulations, we represent the N p terms with the smallestμ i of the product Eq. (3.5) by single pseudofermions, whereas the determinant of the remaining factors is expressed as a single pseudofermion integral Here again, the contribution from the two final terms is always considered together. This decomposition has several advantages. First of all, the small residues frequently can be integrated on a larger time scale, due to a small coefficient decreasing the forces. Furthermore, while the multi-shift conjugate gradient algorithm [41] is efficient for the combined solution of the systems in the last factor with the large shifts, it turns out to be advantageous to employ the deflated solver for the terms involving the smallerμ i . In this case it is no longer necessary to use a single pseudofermion field for all shifts.
The range of the rational approximation is given by the smallest and the largest eigenvalue ofQ 2 over typical gauge field configurations. On thermalized configurations, estimates of these numbers can be obtained in openQCD by the power method applied toQ −2 andQ 2 , respectively. Typically, O(20) iterations proved sufficient for the lower bound, whereas the largest eigenvalue required O(100) iterations. In particular the smallest eigenvalue turned out to be sensitive to thermalization effects and exhibit larger fluctuations than expected. This made it necessary to monitor it carefully at the beginning of each production run.

HMC and the integration of the molecular dynamics
In the algorithm the action is split into different components: the gauge action, the determinants from the Hasenbusch splitting for the light quarks and the various contributions to the strange-quark determinant from the rational approximation described above, N mf + N p + 4 components in total. The complete action is simulated with the Hybrid Monte Carlo (HMC) algorithm [42]; the classical equations of motion are solved numerically for trajectories of length τ = 2 in all simulations. This leads to Metropolis proposals which are accepted with an acceptance rate P acc , given for our runs in Table 2.
The goal of the splitting of the action, and the forces deriving from it, is the reduction of the computational cost needed to obtain a high acceptance rate at the end of the trajectory. The gauge forces are much cheaper to compute than the fermion forces, whose components differ by orders of magnitude in size and fluctuations. It is therefore natural to use a hierarchical integration scheme for the molecular dynamics of the HMC to reflect this [43].
We use the setup described in Ref. [12], i.e., a three-level scheme with the gauge fields integrated on the innermost level with the fourth order integrator suggested by Omelyan, Mryglod, and Folk (OMF) [44] and implemented in the openQCD code. Most of the fermion forces are on the intermediate level, again integrated with the fourth order integration scheme. Only particularly small components of the fermion forces, that contribute little to energy violation, are integrated on a larger scale with the second order OMF integrator [44], whose parameter λ is set to 1/6.
Since one step of an inner level integration scheme is done for each outer step, there are three parameters which define the scheme: the number of outermost steps per trajectory N s,2 and the number of poles N p,2 as well as the number N mf,2 of terms of Eq. (3.3) integrated on the outermost level. In the latter two cases the numbers refer to the terms with the smallest twisted-mass shifts. The values chosen in our runs can be found in Table 2.
The choice of the trajectory length affects the autocorrelation times and is therefore not easily studied. In general, longer trajectories have proven to be beneficial [4], but in particular with dynamical fermions one might prefer shorter trajectories because of instabilities of the integrator. As a compromise, we use τ = 2. Asymptotically, this leads to autocorrelations growing with τ int ∝ a −2 . Note, however, that this scaling behavior is also expected if the length of the trajectory is scaled [45].

Solver
The extensive use of the locally deflated solver [46][47][48] is an important part of the progress that made the presented simulations possible. It removes the largest part of the cost increase as the quark mass is lowered, thereby circumventing the significant slowing down observed in the past. The increase in performance of the solver comes at the price of a more complex setup and many additional parameters which have to be chosen.
Fortunately, relatively little tuning of the local deflation subspace was necessary here and we therefore do not list the parameters in detail. For most runs, we used deflation blocks of size 4 4 . The parallelization of the N s = 48 lattices required one or two dimensions to be set to 6; also blocks of 8 × 4 3 have been used.
The number of deflation modes per block has been chosen between 20 and 32, in order to balance the higher efficiency provided by the larger subspace and the cost associated with the application of the preconditioner.
For the smaller lattices with L/a = 32, we set the solver accuracy (the ratio between the norm of the residue to the norm of the right hand side of the equation) to 10 −11 in the action and 10 −10 in the force computation. To ensure the value of the action and the reversibility of the integration of the equations of motion are sufficiently precise, more stringent residues have been used for the lattices of larger volume.

Autocorrelations
Markov Chain Monte Carlo algorithms, like the Hybrid Monte Carlo used here, produce field configurations which exhibit autocorrelations characterized for an observable A by the autocorrelation function where t is the Monte Carlo time. The integral over the normalized autocorrelation function ρ(t) enters the error analysis. This is the integrated autocorrelation time To estimate τ int (A) with a finite variance, it is necessary to cut the summation at a window W [49,50]. In order to choose the window for our final error estimates and to account for the thereby neglected tail, we employ the method described in Ref. [4]. Its essential input is an estimate for the exponential autocorrelation time, which we discuss in the following.

Scaling of the autocorrelations
As we approach the continuum limit, the autocorrelation times are expected to grow due to critical slowing down. The open boundary conditions used in our setup should prevent catastrophic scaling due to the freezing of the topological charge. Since we have chosen the trajectory length constant in all our runs, we expect Langevin scaling τ int ∝ a −2 .
In Figure 2 we show autocorrelation times of notoriously slow observables: the global topological charge and the action density averaged over the plateau region, They are both defined in Eq. (6.2). We find a situation similar to that encountered in pure gauge theory [7]. While the energy density shows very good scaling, the topological charge decorrelates faster on coarser lattices. The fast growth of the integrated autocorrelation time of the charge does not mean that the 1/a 2 scaling is not valid. In pure gauge theory, the behavior could very well be fitted with τ int ∝ a −2 (c + da 2 ). In this picture, there are significant cutoff effects to the scaling, but no catastrophic behavior in the a → 0 limit. This is expected when simulating with open boundary conditions.

Cost of the simulation
Note, in two-flavor QCD with periodic boundary conditions at a lattice spacing of roughly a = 0.05 fm [8] the topological charge does not decorrelate slower than the smoothed energy. Rather, it shows similar autocorrelations for quark masses around 400 MeV. This means that we are not yet in the position to fully profit from the effect of the open boundary conditions, however, going to finer lattice spacings the freezing observed in two-flavor QCD at a ≈ 0.03 fm [51] will be avoided.
In the sense of fast decorrelations and minimal requirements on the number of units of molecular dynamics time (MDU), the presented simulations are not cheap, nevertheless. The exponential autocorrelation time of τ exp ≈ 14(3) t 0 /a 2 is consistent with what is found in Fig. 2. For biases to be small and a simulation to be reliable we need a total Monte-Carlo history of at least O(50)×τ exp . For β = 3.4 this translates to 2000 MDU, whereas for β = 3.55 and β = 3.7 a statistics of 3600 MDU and 6000 MDU, respectively, is necessary. For most of our ensembles listed in Table 1, we exceed these numbers, but for some, which are still in production, they are not yet reached. Those quoted results therefore have to be taken with care in these cases.

Reweighting factors
The simulations are not done with the exact QCD action as given by Eqs. (2.1) and (2.3), but differ due to the twisted-mass reweighting Eq. (3.4) and the inaccurate rational function in the RHMC Eq. (3.6). The observables are reweighted to the target theory, for which the reweighting factor W = W 0 W 1 needs to be computed. The factors W 0 and W 1 , as defined in Eqs. (3.2) and (3.5) respectively, contain ratios of determinants which are estimated stochastically as described below.
Expectation values A of primary observables A can then be computed from expectation values in the theory with the modified action · · · W , according to

Twisted-mass reweighting factor
The twisted-mass reweighting plays an important role in our setup. From a conceptual point of view, it removes barriers of infinite action created by zero eigenvalues of the Wilson Dirac operator. Together with Hasenbusch factorization, it also reduces the fluctuations of the forces which makes the simulations cheaper and more reliable in practice [38]. This situation is especially favoured for a large value of µ 0 in Eq. (3.2). It might, however, also lead to significant fluctuations in the reweighting factor and as a consequence to a larger statistical error of certain observables.
As a consequence, what constitutes the optimal choice of the parameter µ 0 will in general depend on the observable. As can be seen in Figure 3, W 0 is close to a constant for most configurations; only on some configurations the value will be much smaller. For observables with little or no correlation to the reweighting factor, like the gluonic ones we consider below, this effectively amounts to a reduction in statistics [52]. This reduction is negligible for our ensembles since var(W ) W 2 in all cases.
For observables with a strong correlation with W 0 , the situation is more delicate. Even after reweighting, this can lead to large fluctuations in the measurements and significantly increased statistical error. In particular in the case of anticorrelation, the situation is more problematic due to the stochastic estimation of W and, possibly, the observable. The cancellation between, e.g. a large value of the observable and a small value of W might require a rather precise determination of the two.

Reweighting and the pseudoscalar correlation function
The pseudoscalar correlation function is an observable showing a strong anticorrelation between its value and the reweighting factor. This can be easily understood by noting that at small quark masses both receive significant contributions from the smallest (in magnitude) eigenmodes of the Hermitian Dirac operator. It is precisely this region where the reweighting term has the largest effect. and two values of µ 0 is shown. As we can see, the larger µ 0 leads to larger fluctuations in W and f PP (x 0 ), as expected. In the product, however, they cancel and the average value W f PP (x 0 ) / W is then consistent within the statistical errors between the two ensembles.

Computation of W 0
Since the determinant ratios needed for the computation of W 0 (µ 0 ) cannot be computed directly, a stochastic estimator is taken instead. This can either be done by directly estimating the determinant ratio in Eq. (3.2) or by first splitting it up and then using stochastic estimates for the individual factors [53], a strategy already successful in Hasenbusch's mass factorization.

(5.2)
Now each of the factors is evaluated stochastically with N r complex-valued Gaussian random fields η of unit variancẽ such that up to an irrelevant constant factor the determinant is retrieved by averaging over the noise fields Following the initial proposal of Ref. [12], it is sufficient to use a single step N sp = 1 with a suitably chosen value of N r . Its value along with the other parameters of the reweighting can be found in Table 3. This is the method implemented in openQCD-1.2.
Once the fluctuations in the reweighting factor increase, it is advisable to use intermediateμ, a possibility given in openQCD-1.4. This is because the distribution of the results for the reweighting factors become long-tailed once exceptionally small eigenvalues of theQ 2 are encountered. In this situation it is very difficult to argue about the uncertainty of W 0 [54]. By splitting the estimate into smaller intervals iñ µ, the distribution of each of the factors becomes significantly more regular.
For ensemble H105r002 we find precisely such a situation. While with a single step inμ the smallest reweighting factors show a distribution which is far from Gaussian, using ten intermediateμ the individual factors can be computed reliably to O(15%) accuracy using 15 sources each.

RHMC reweighting factor
Since the rational approximation has been chosen to a good accuracy, the fluctuations in the reweighting factor are small and it turns out to be sufficient to estimate it with one stochastic source per configuration. The associated variances are given in Table 3. They are seen to receive a considerable contribution from the stochastic estimation of W 1 .
In order to study the effect of more sources, we observe using five instead of one stochastic estimate reduces the variance of W 1 by more than a factor 4, on ensemble  H105r005. The same is true for the H200 ensembles. Still, even with one source per gauge configuration the noise introduced by W 1 is negligible for all observables we investigated. Note that in some early runs we underestimated the upper bound of the interval in which the rational function is accurate. Since the accuracy does not deteriorate quickly outside the interval, the fluctuations of the reweighting factors nevertheless are sufficiently small.

Wilson flow
The Wilson flow can be a very useful tool in lattice QCD from which quantities with a finite continuum limit can be constructed [6,55,56]. The gauge fields U (x, µ) are subjected to the smoothing flow equation with S W the Wilson action. With clover-type discretization of the field strength tensorĜ µν (x, t) constructed from the smooth fields V t , the time slice energy E(x 0 , t) and the global topological charge Q top (t) can be constructed With the vacuum expectation value of the energy E(t) , the scale parameter t 0 is then defined by Throughout this paper we quote the observables of Eq. (6.2) at flow time t 0 .

Effects of the boundary
Due to the open boundary conditions in the temporal direction, time translational invariance is lost. Sufficiently far away from the boundaries, local observables are expected to assume their vacuum expectation values up to exponentially small corrections with a decay rate equal to the lightest excitation with vacuum quantum numbers.
On top of this continuum boundary effect, large discretization errors are observed close to the boundary. As an example we show in Figure 4 the behavior of the smoothed energy E(t, x 0 ), defined in Eq. (6.2). A further example for the pseudoscalar correlation function with the sink approaching the boundary can be found in Ref. [52].
In the case of the energy, it should be noted that it is at this point difficult to disentangle the discretization effects in the underlying gauge field from the ones introduced by the Wilson flow and the observable used to define the field strength tensor, but recent work by Ramos and Sint might clarify this issue [57]. Also the Dirac operator used in the measurement of the effective mass is only tree-level improved at the boundary. Still, the effects of the finite lattice spacing are very prominent at our coarser lattices, but become much less notable as the continuum limit is approached.
As can be seen in Figure 4, no sizable dependence on the quark mass is observed in t 2 0 E(x 0 , t 0 ). This is trivial for the bulk, since its value is equal to 0.3 by definition. But also the (cutoff) effects close to the boundaries show no quark mass dependence. Whether this is a generic feature of the sea quark contribution being small or it is due to our particular choice of chiral trajectory Eq. (2.8) cannot be judged from the data presented here.
In the context of the present paper we will not discuss these effects in detail, but perform the measurements in a region where they can be neglected. The determination of the plateau region is not always clear due to an effect already observed in Ref. [12]: is expected if the statistical analysis is done properly, however, they do make the plateau determination more difficult.
For meson correlation functions, these waves have been discussed previously [58], see Figure 6 for an example from our simulations. In other simulations they are typically not visible, because time translational invariance is used on the level of the correlation function by using different source positions in time and averaging them before computing effective masses. Again, this is not a principal problem. However, we need to ensure sufficient statistics and that errors are under control and require a procedure to deal with these waves.

Measurement of t 0
For the determination of t 0 , we need to determine the plateau in E(x 0 , t) for t = t 0 . Since we are looking at a smoothing radius of √ 8t 0 , the effects of the boundary visible in Figure 4 are at the expected length scale. Discretization effects are large though, and it is therefore difficult to argue about the expected functional form. In this situation, we use a two-stage procedure: First we fit in the range where this ansatz describes the data. This is only used to determine the fit range by the condition that in the whole range the statistical uncertainty δE(x 0 )  Table 4. In Figure 5 the quark mass dependence of t 0 is given for the three available lattice spacings. Recall, the values are given in terms of the symmetric point which defines the chiral trajectory at m ud = m s . From Taylor expansion around this point [14] as well as ChPT [27] one expects a constant behavior to the respective leading order. This is confirmed for the finer lattices to our level of accuracy. Only at the coarsest lattice spacing, cutoff effects seem to cause some deviation, albeit on a rather small scale.

Pseudoscalar masses
The masses of the pseudoscalar particles are computed from the pseudoscalar correlation function projected to zero momentum. With quark fields of flavor r and s,   and the pseudoscalar density P rs =ψ r γ 5 ψ s , it is given by x, y P rs (x)P sr (y) . (6.7) Due to the open boundary conditions in time, the translational invariance in the temporal direction is broken. However, we find that there is little to gain from using source fields at different time slices [52]. The U (1) stochastic source fields are therefore put only at y 0 = a and y 0 = T − a [59]. In the following we analyze In the continuum limit and for large volume and sink positions far away from the source and boundary, x 0 0 and x 0 T , the two-point function is expected to fall off as [12] f PP ( In line with Ref. [12],T is a free parameter. We follow a similar strategy as in Section 6.3 to make sure that in our final fit the excited state contribution is negligible. We show an example of an effective mass plot in Figure 6, where we can see that this fit works very well in a wide range of x 0 . The results for the masses are listed in Table 4.

Comparison to simulations with periodic boundary conditions
One possible concern regarding the open boundary conditions is that the region close to the boundaries is large and as a consequence one loses a sizeable fraction of the statistics. Which fraction of the lattice needs to be discarded depends on the observable and the statistical accuracy of the data, with boundary effects expected to decay close to the chiral limit as exp(−2m π x 0 ). However, it should be noted that also the systematic finite volume effects are parametrically similar with contributions proportional to exp(−m π L): A high accuracy requires large lattices in both, the temporal and spatial directions, which is also true for simulations with periodic boundary conditions in time.
In any case, since we observe large cutoff effects close to the boundary, arguments based on continuum physics are problematic at current lattice spacings. This is already visible in the x 0,min /a chosen in the measurement of E(t 0 ), which has only a minor dependence on the lattice spacing. On our smallest lattices at β = 3.4 with N t = 96, the x 0,min /a = 20 leads to a plateau average of almost 60% of the total time extent. On all other ensembles we have an even larger fraction over which we can take the plateau average.
In our measurements of the pseudoscalar masses we typically start the plateau at x 0,min ≈ T /4, from where on the effects from the excited states can be neglected. As noted before, moving the source away from the boundary has little effect, since the plateau is seen to start at the same position. The minimal distance x 0,max of the sink from the boundary is typically around T /6, such that we have in total a plateau stretching between 50% and 65% of the lattice. Even if the other half of the lattice was completely decorrelated, this would at most correspond to a factor of two in statistics.

Scaling violations
In the bulk, our action is fully O(a) improved, only for the boundary terms we use the tree-level values. This guarantees leading scaling violations close to the continuum limit to be of order a 2 , but at finite lattice spacings higher order terms will always be present as well. How large their contribution is and whether one can safely neglect them given the statistical uncertainties of the simulation is not a priori clear. In any case, once the higher order terms become important, they limit the value of coarser lattices in the continuum extrapolation.

Cutoff effects in t 0
A particularly precise way to study discretization effects is to look at observables which agree in the continuum limit but differ at finite lattice spacing. To this end, we take two slightly different definitions of t 0 : both are given by the implict relation Eq. (6.3), where in one case we use the conventional "clover" discretization of the field strength tensor for the the energy density E, Eq. (6.2), in the other case the plaquette definition is used as given in Ref. [6].
Since the continuum value of t 0 has to be the same, the ratio of t clov 0 and t plaq 0 has to be one up to cutoff effects. As they are evaluated on the same gauge field configurations, the two values of t 0 are highly correlated such that their ratio can be evaluated to exceedingly high accuracy.
As we can see in Figure 7, the ratios at β = 3.7 and β = 3.55 agree with the a 2 scaling hypothesis up to very high accuracy, with a total deviation of the ratio from its continuum value of 4% and 6%, respectively. With the assumption that higher order effects are negligible at β = 3.7, one concludes that at β = 3.55 higher orders contribute 0.4% to this observable, while at β = 3.4 an additional O(2%) effect can be attributed to higher orders on top of the 11% which come from the leading order scaling violation.
While this additional O(2%) effect originating from the higher order terms is not of concern for most observables in current lattice computations, it might impact studies of certain high accuracy observables.
In any case, the large discrepancies between the two definitions of t 0 are a problematic finding in view of the fact that our tuning strategy is entirely based on this quantity to set the scale. We therefore have to expect that for other scale setting strategies, the matching of the chiral trajectories will differ on the order of the ratio observed at the level of 10% at the coarsest lattice spacing.

Coarser lattices
In order to investigate the value of ensembles at coarser lattice spacings, we have also generated some β = 3.3 lattices at the m ud = m s matching point. After some tuning, κ ud = κ s = 0.136423 is found to match the φ 4 = 1.15 point to reasonable accuracy. However, the observation of large cutoff effects on 96 × 24 3 lattices, indicating this point no longer being in the assumed a 2 scaling region, has led us to abandon this coupling for now.
This decision is based on Figure 7, where for this parameter point we find a 2 /t 0 ≈ 0.5 and a ratio of t clov 0 /t plaq 0 ≈ 1.3, which is 16% above the leading order violations of 1.15. We are therefore clearly no longer in the scaling regime and since we aim with most observables at accuracies much below the 10% level, the points at this lattice spacing of roughly 0.1 fm, do not meet our precision goals. Furthermore, the autocorrelations observed in particular in the thin link plaquette were very significant and also large fluctuations in the lower spectrum of the Dirac operator have been observed. This makes these lattices difficult to simulate and poses another reason to refrain from considering this value of β at this time.

Conclusions
The generation of the gauge field configurations described here lays the ground for many future lattice QCD calculations. It is the first time that open boundary conditions in time and twisted-mass reweighting have been extensively used in such large scale calculations.
The two methods have been shown to work well. Keeping in mind simulations which are on our roadmap for the future, the experience gained with the use of open boundaries will prove very valuable as the lattice spacing is decreased, while not being strictly necessary at the lattice spacings under investigation.
We could show, that no particular obstacle is posed by the boundaries themselves, however, significant discretization effects are observed in their vicinity. Depending on the observable and its correlation with the reweighting factor, we observe the twisted-mass reweighting is under control. To study this in the future, we have generated ensembles with different values of the reweighting parameter µ 0 .
It is noteworthy that similar data sets previously needed to be accumulated over many years. However, due to advances in hardware and in algorithms, we could demonstrate the progress that has been made, by generating the current, new data set within a year and a half after the parameters of the action had been determined.
As of now, the covered parameter space is limited: We only have data on one chiral trajectory, a limited range of quark masses and lattice spacings and typically only one volume. In order to better control the associated systematic uncertainties, we therefore plan to extend the current set of ensembles. We are certain the configurations presented here will prove useful, and we are excited about the interesting physics results that will be obtained.