An exploratory study of heavy domain wall fermions on the lattice

We report on an exploratory study of domain wall fermions (DWF) as a lattice regularisation for heavy quarks. Within the framework of quenched QCD with the tree-level improved Symanzik gauge action we identify the DWF parameters which minimise discretisation effects. We find the corresponding effective 4$d$ overlap operator to be exponentially local, independent of the quark mass. We determine a maximum bare heavy quark mass of $am_h\approx 0.4$, below which the approximate chiral symmetry and O(a)-improvement of DWF are sustained. This threshold appears to be largely independent of the lattice spacing. Based on these findings, we carried out a detailed scaling study for the heavy-strange meson dispersion relation and decay constant on four ensembles with lattice spacings in the range $2.0-5.7\,\mathrm{GeV}$. We observe very mild $a^2$ scaling towards the continuum limit. Our findings establish a sound basis for heavy DWF in dynamical simulations of lattice QCD with relevance to Standard Model phenomenology.


Introduction
With LHCb and BESIII generating data, and Belle-II soon to start production, increasingly accurate Standard Model (SM) predictions for heavy flavour physics are dearly needed to constrain or hopefully identify new physics. These predictions typically involve matrix elements of the operators of the Weak Effective Hamiltonian among hadronic states. As a result they require a non-perturbative approach, making lattice QCD simulations crucial. This is why in the last few years several approaches to implement heavy quarks in simulations of lattice QCD have been proposed. Some of these are based on an effective description of the heavy degrees of freedom, such as the Non Relativistic treatment of the heavy quark (NRQCD) [1,2] or Heavy Quark Effective Theory (HQET) [3,4], or on a non relativistic re-interpretation of relativistic discretisation [5][6][7][8]. More recently collaborations have started treating the charm and bottom quarks in the same relativistic framework used to discretise the light quarks, e.g. [9,10].
However, simulations of full lattice QCD where both the physical light quarks (u, d and s) and the charm or heavier quarks are represented by the same discretisation are still rather scarce. The main reason is certainly that the relevant energy scales, associated with the pion and the heavy quark masses respectively, are computationally costly to reconcile. This is particularly true in a fully relativistic and dynamical setup with controlled finite volume and discretisation errors. Simulations in which all quarks are discretised in the same way have a number of advantages, though. For instance, continuum flavour symmetries at finite lattice spacing simplify many calculations. Moreover, only such a setup seems suitable for the study of GIM-cancellation, which is an important ingredient in a number of phenomenological applications [11].
This paper is the second [12] in a series towards a lattice phenomenology program with domain wall fermions (DWF) [13,14], in particular Möbius DWF (MDWF) [15][16][17], as the discretisation for light as well as heavy quarks. Compared to Twisted Mass [18], DWF offer the attractive properties of conserving both chiral and parity symmetries at finite lattice spacing. Compared to HISQ fermions [9], a single quark can be simulated without the need of taking the root of the determinant to eliminate the different tastes, thus providing a theoretically clean regularisation.
Since we enter mostly uncharted territory with simulations of heavy DWF (see also [12,19,20]), we dedicate this paper to the investigation of its basic properties. We are particularly interested in studying the approach of heavy-light meson observables to the continuum limit. Our simulations have been carried out within quenched QCD. This is computationally much cheaper than dynamical QCD and therefore allows us to access a much wider range of lattice spacings (a −1 ≈ 2.0−5.7 GeV). While the quenched approximation is certainly not suited for making phenomenologically relevant predictions, we expect it to share a number of properties with the unquenched case. Most importantly, we expect that the continuum limit scaling observed in the quenched theory over a large range of lattice spacings will be qualitatively the same as in the dynamical theory. Such information is particularly valuable given that for phenomenologically relevant simulations only dynamical ensembles at coarser lattice spacings are currently available.
The rest of the paper is organised as follows: in section 2 we outline the overall computational strategy followed in this paper, report on the properties of the generated quenched gauge field ensembles, define the quantities that we compute and discuss several more technical aspects of our computation. In section 3 we describe the tuning of the MDWF parameters. This is followed in section 4 by a study of the continuum limit scaling of the dispersion relation and decay constants. In section 5 we draw our conclusions. In the appendix we provide supplementary material, in particular the numerical values for all data underlying the analysis.

Strategy
The main purpose of this work is to gain a qualitative understanding of discretisation effects of heavy MDWF. To this end, we study the MDWF parameter space and the heavy quark mass dependence of basic heavy-heavy and heavy-strange meson matrix elements and the energy as the cutoff is varied. Simulations of the quenched theory allow us to adopt algorithms (over-relaxed [21,22] heat-bath [23]) that are, compared to the algorithms used with dynamical quarks (Hybrid Monte Carlo [24]), computationally much cheaper. To some extent the problem of critical slowing down [25][26][27] can therefore be circumvented by brute force. This enables us to probe finer lattice spacings than those affordable in dynamical simulations and check the scaling of the theory towards the continuum limit in more detail. In order to reduce simulation costs further, a relatively small physical lattice volume of L ≈ 1.6 fm was considered. The volume was kept approximately constant while decreasing the lattice spacing. Since the finite size effects in physical quantities are then constant across all simulated lattice spacings, cut-off effects can be studied in detail.
An important point addressed in this study concerns the residual chiral symmetry breaking of MDWF. The restoration of chiral symmetry in the massless limit is crucial to the simulation of QCD on the lattice, and is also responsible for automatic O(a)-improvement, which is especially important when studying heavy quark physics. In our notation, the five dimensional MDWF action is and we define with the usual chiral projectors P ± = 1 2 (1 ± γ 5 ) and the Wilson matrix Besides the bare quark mass am, MDWF have two further input parameters that need to be specified in each simulation: the extent of the fifth dimension L s and the domain wall height parameter M = −M 5 , respectively. More specifically, M 5 is the negative mass parameter in the 4-dimensional Wilson Dirac operator that resides in the 5-dimensional MDWF Dirac operator. Since both L s and M 5 are parameters of the discretisation rather than of QCD we have some freedom in varying them. In the limit L s → ∞ and with the Wilson kernel this formalism coincides with the overlap formulation [28,29] and allows for the simulation of a four-dimensional chirally symmetric theory (in the limit of massless quarks) that is free of doublers. When L s is finite however, chiral symmetry remains broken by a small amount. 1 This can be quantified by measuring the amount of additive quark mass renormalisation, also known as residual mass m res (defined later in 2.3). For a given extent of the fifth dimension, the parameter M 5 sets the scale for the exponential localisation of the chiral modes of the fermionic fields at the boundaries of the 5th dimension. The decay rate of the physical mode away from the boundary is however also modified by the presence of an explicit quark mass term, and care must be taken in order to maintain the localisation of the physical modes on the boundary [13,20,[31][32][33][34]. As we will see, this becomes particularly crucial for heavy input quark masses. We will study how the choice of a heavy quark mass am = am h and M 5 changes the ultra-violet properties of the discretisation. In the following we chose an extent L s = 12 of the fifth dimension, which guarantees a small value of m res for light quarks [35]. The particular choice of MDWF is the same implementation as the one used in [35] with a Möbius scale of α = b + c = 2.

Ensemble generation
We generated ensembles based on the tree-level Symanzik improved [36,37] gauge action with lattice spacings in the range of 0.034-0.1 fm. The gauge configurations have been produced with the heatbath algorithm [21][22][23]. The coarser three ensembles were generated using CHROMA [38] Table 3. Autocorrelation time of topological charge (τint (Qtop)) and of charge squared (τint Q 2 top ) in units of sweep steps; number of sweeps separating each configuration included in the measured ensemble (Nsep), and total number of gauge configurations considered.
for the finest lattice spacing (which involved the highest computational cost) we recurred to a faster implementation, especially optimised for IBM BG/Q [39]. In tables 1, 2 and 3 we summarise the simulation parameters used and basic ensemble properties.
Lattice spacings have been determined at each simulated β by enforcing the Wilson-flow scale w 0 [40,41] to take its "physical" value, which we assumed to be w phys 0 = 0.17245(99) fm as recently determined in [35]. 3 We kept the physical volume fixed such that the spatial extent remained at about 1.6 fm (cf. table 2).
The evolution of the topological charge Q (measured with the GLU package [42]) is illustrated in figure 8 in appendix A. These quantities are expected to couple strongly to the slowest evolving mode in the evolution of the algorithm [26]. We obtain sets of decorrelated measurements by choosing only configurations for further processing that are separated by N sep intermittent update steps with N sep larger than twice its autocorrelation time τ int Q 2 [43], as detailed in table 3.

Observables
The pseudoscalar decay constant f X is defined as the matrix element of the conserved MDWF axial vector current [44] between a pseudoscalar meson state X and the vacuum, We determine the decay constant f X and the energy E X (p) of the pseudoscalar state X from fits to the time dependence of Euclidean QCD two-point correlation functions projected onto momentum p, (2.4) The operator O si M is an interpolating operator with the quantum numbers of the meson, i.e. O s M = q 2 ω s Γ M q 1 , where we consider the pseudoscalar case Γ P = γ 5 and the axial vector case Γ A = γ 0 γ 5 , respectively. The superscript s indicates the smearing type induced via the spacial smearing kernel ω, which in the simulations presented here is either local (s = L, ω(x, y) = δ x,y ) or Gaussian via Jacobi iteration [45][46][47] (see table 4 for our choice of smearing radii). The constants The fits leading to the extraction of masses and decay constants are multi-channel fits to combinations of the two-point correlation functions C AA , C AP , C P A and C P P . We note the relation between the conserved MDWF axial current [35,44]

and the renormalised local axial current
A further quantity that we wish to monitor during our simulations is the residual quark mass am res [44], which provides an estimate of residual chiral symmetry breaking in the MDWF formalism. It is defined in terms of the axial Ward identity (AWI) where ∆ − µ is the lattice backward derivative and am is the bare quark mass in lattice units in the Lagrangian. It motivates the definition Here, J 5q is the pseudoscalar density in the centre of the 5th dimension. We compute the correlation functions in eq. (2.4) with two types of quark sources. The analysis of the decay constant and the residual mass is based on Z 2 noise sources and the one-end-trick [48][49][50] (in this case we only consider p = 0) while the analysis of the dispersion relation is based on point source data. The computation of heavy quark propagators by means of conjugate gradient type algorithms can be affected by roundoff errors [51]. We monitor proper convergence during the computation of the quark propagators by checking that the desired solver residual is fulfilled on all time slices using the time slice residual [51] defined as where |x| t is the norm of the vector x restricted to time slice t. We determined statistical errors using the bootstrap method with 500 samples.  Table 4. Simulated strange and heavy input quark masses am h and the choices of smearing radii for heavy quark masses. The simulated bare quark masses are quoted in lattice units for the MDWF action. The heavy quark masses starting from "start" with a step of "step" and ending at "end" are simulated. r P sm and r Z 2 sm refer to the choice of the smearing parameter for the Gaussian smearing of the source/sink of the propagators for the point and Z2 noise sources, respectively. For the Gaussian smearing we employed 400 Jacobi iterations. All measurements are carried out with MDWF with parameters Ls = 12.

Tuning MDWF for charm
In this section we present results for the am h and M 5 dependence of the heavy-heavy meson decay constant f hh and the residual mass am res .

M 5 dependence
The left hand panel of figure 1 shows the dependence of the heavy-heavy decay constant on the heavyheavy inverse pseudoscalar mass m hh observed on the coarsest (β = 4.41) ensemble. We normalise the results for a given M 5 by the value of the decay constant at m hh = 1.5 GeV as obtained from a polynomial interpolation. For small values of m hh the decay constant shows little dependence on the value of M 5 , but as m hh is increased a strong dependence is observed.
The right hand panel of figure 1 shows the same data for

Residual mass
Next we quantify how the residual chiral symmetry breaking is affected by M 5 by observing the response of the size of the residual mass to variations in am h and M 5 . In the left panel of figure 2 we show the ratio of correlation functions eq. (2.6) from which we determine am res as a function of time for several values of the quark mass at M 5 = 1.6. Note that for large t the time dependence in ratio eq. (2.6) is expected to cancel between the numerator and denominator. While the expected (constant) behaviour in time is observed for small quark masses, this is strikingly not the case for values of am bare h 0.4. In these cases it is difficult to interpret the operator J 5q 's matrix element as a constant, residual additive mass correction in the chiral Ward identity. The effect is of course rather small compared to the explicit chiral symmetry breaking, but there is a risk that the physical modes no longer remain bound to the walls of the fifth dimension in this large mass limit. To be more quantitative, we define am res (t = T /2) as the value of this correlator ratio in the (temporal) middle of the lattice. Note, however, that above am h ≈ 0.4 the meaning of am res as a unique measure of residual chiral symmetry breaking is no longer clear, only indicative. The right hand panel in figure 2 shows am res (t = T /2) as  a function of the quark mass. We observe the same qualitative behaviour for all values of M 5 : as the input quark mass is increased beyond am h ≈ 0.4 the residual mass am res (t = T /2) starts to increase drastically. Although this quantity is L s dependent, it is likely unsafe to use domain wall fermions at masses where the physical modes become unbound from the walls and the matrix elements of J 5q have such non-trivial behaviour. The impact on 4d observables will be studied later in this paper.

Locality of the effective 4d Dirac operator
Given the above observation indicating the reduced binding of surface states of MDWF above am h ≈ 0.4, a further concern one might have is that we should check the locality property of the corresponding effective 4d MDWF Dirac operator. The connection of the 5d MDWF operator D 5 M DW F defined in eq. (2.1) to a four dimensional effective theory is well established in the literature, [15-17, 35, 52-54]. We identify D ov as an approximation to the overlap operator with approximate sign function where the Möbius kernel is The transfer matrix in the fifth dimension can be identified as The effective overlap operator may be simply found as where this is known to reduce to the standard overlap formalism in the limit L s → ∞ and when b = c, and the projection matrix P is Following eq. (3.5), we may place the mass dependence of D ov (am) at non zero mass in the following form: We see that the kinetic term in the four dimensional effective action should remain unaltered as the mass is changed up to a trivial rescaling factor (1−am) affecting the surface field renormalisation. The induced overlap bare mass is therefore better interpreted as the combination am 1−am , which of course varies non-linearly and diverges as we take the domain wall mass towards the Pauli-Villars mass of unity. The exponential locality [55] is fully encoded in the massless operator, and is independent of the quark mass. So, from this perspective there should be no locality issues as we take the mass large, since the kinetic term is trivially rescaled compared to the light mass case.
We demonstrate this with a second use of eq. (3.5). The effective operator may be constructed by the simple application of the inverse of the Pauli Villars operator. Following the methodology of ref. [55] we now study the locality properties of this operator.
We start by defining a point source ξ, ξ α,a (x) = 1 x = y, α = a = 0 (spin, colour) 0 otherwise, (3.10) where y is the source location, and ψ is the result of the multiplication of the effective 4d Dirac operator with ξ, ψ = D ov ξ . (3.11) We say D ov is strictly local (or "ultralocal") if the only non-zero contributions to (D ov ξ) (x) come from a finite set of terms D ov (x, y) ξ (y) with y in the vicinity of x [55]. We collect all lattice points {x} r separated by r hoppings from the origin, such that x ∈ {x} r if |x| 1 = r. Here |x| 1 is the "taxi driver" (or "Manhattan") norm of x, defined by 12) where N µ is the number of lattice sites along the µ axis. This definition accounts for the periodicity of the lattice. Finally, for each value of r we define the maximum of the norm of ψ at the set of points In figure 3 we show the function f (r) for two bare quark masses on all three ensembles. As expected, we observe that the slope of f (r) is independent of the bare quark mass as well as of the lattice spacing, indicating that locality is recovered in the continuum limit.
We can make a more quantitative statement for the mass independence of the locality of D ov (am): motivated by eq. (3.8) we define the functionf : where we have introduced a term to subtract the additive mass term in eq. (3.8). We can then define the ratio where the subscripts indicate the bare quark masses at which the function f was evaluated (am 1 = 0.1 and am 2 = 0.5). According to eq. (3.8), we expect R(r) = 1, which is confirmed by our data to the level of arithmetic precision used in the computation. This provides a strong consistency check of our setup and our understanding of the locality of the MDWF operator.

Continuum limit of the decay constant and the dispersion relation
The results in the previous section provide the first evidence for a region in parameter space where MDWF can be used as a suitable discretisation for heavy quarks. To further substantiate this picture we now fix M 5 = 1.6 and study the continuum scaling of a basic heavy-strange pseudoscalar meson matrix element, the decay constant, and the corresponding dispersion relation, as a function of the mass of the heavy quark.

Choice of strange and heavy quark masses
We study the continuum limit along lines of constant strange and heavy quark mass. We fix the s-quark by considering a fictitious meson η s composed of two different quarks, s and s , of degenerate mass m s . This meson differs from the physical η −η mesons by quark-disconnected Wick contractions. We tuned the strange quark mass to its "physical value", by imposing at each lattice spacing the mass of the simulated η s meson to reproduce m ηs = 0.6858(40) GeV [56]. This sets a common renormalised strange quark mass on all the ensembles. In table 4 we report on the values of the corresponding bare strange quark mass and on our choices for the simulated heavy quark masses.

Decay constants for heavy-strange mesons
We consider the renormalised ratio where we introduce f norm sh m norm sh , interpolated to m sh = 1 GeV, to cancel the axial current renormalisation constant. We also include in both the numerator and denominator a factor of √ m sh to make both of these quantities individually finite in the limit am h → ∞. We interpolate R sh to the reference pseudoscalar masses 1.3, 1.6, m Ds = 1.9685 [57] and 2.4 GeV on all ensembles. To fulfil the constraint am h ≤ 0.4 we are forced to drop the coarsest lattice spacing for the heaviest mass considered. A first visual inspection (see figure 4) suggests the absence of cutoff effects beyond O(a 2 ). Moreover, cutoff effects are observed to be very mild for the choice M 5 = 1.6, in agreement with the observation made in section 3.
To obtain a more quantitative understanding we perform continuum limit extrapolations by considering two different fit ansätze, namely  Table 5. Results of the continuum limit extrapolation for the heavy-strange decay constants. The first block summarises the results for the linear extrapolation in a 2 , the second block the quadratic extrapolation in a 2 . We also show corresponding results for the χ 2 /dof and p-values.
The results are illustrated in figure 4 as solid and dashed lines with error bands, respectively, and the resulting fit coefficients are listed in table 5. For the two lightest reference masses, 1.3 and 1.6 GeV, the slope of the continuum limit is compatible with zero. For higher masses the continuum limit starts exhibiting a significant slope. In fact, the dimensionless term D 2 a 2 /R (a = 0), which indicates the fractional amount of discretisation errors, is around 2% for the physical D s meson on the coarsest ensemble (a −1 ≈ 2 GeV), and of O (1%) on the next finest one (a −1 ≈ 2.9 GeV). At the level of statistical precision achieved here the fits reveal only a very mild sensitivity to higher order (O(a 4 )) coefficients: E 4 is compatible with zero within one standard deviation.

Dispersion relation
On the lattice, the continuum dispersion relation for pseudoscalar mesons is modified: all even powers of the lattice spacing with p-dependent coefficients, invariant under hypercubic group transformations (e.g. p 2 , µ p 2n µ ...), are allowed. Here we investigate whether the continuum expression is correctly reproduced after taking the continuum limit of the lattice data for the heavy-strange meson energy at various momenta p = 2π L n. In particular, we consider the cases n ∈ {(0, 0, 0) , (1, 0, 0) , (1, 1, 0) , (1, 1, 1)}. The measured meson energies are sufficiently precise to be sensitive to the slight mistunings in the physical volume of our ensembles (cf. table 2). In particular, for any given n the simulated lattice momenta p sim in physical units only agree approximately amongst the different ensembles.
We correct for this by defining a reference volume with spatial extent L ref = 1.648 fm and therefore reference momenta The meson energies E sim are interpolated to this by taking advantage of the lattice dispersion relation Considering eq. (4.4) for a meson of momentum p ref on two different volumes we obtain the interpolated energy: (4.5) In figure 5 we show an example of the interpolation to the chosen reference momenta for the ensemble requiring the largest corrections (cf. table 2) for the case of the physical-mass D s meson.
We now proceed to perform the continuum limit extrapolation of the meson energy. In figure 6 we illustrate the extrapolation of the physical D s meson energies for two different momenta. In both cases the extrapolated result is compatible with the energy predicted by the continuum dispersion relation (4.3).
This procedure was repeated for all momenta and reference masses. In figure 7 we show the results for the energies after the continuum limit extrapolation for the different momenta and choices of reference rest masses. The expected continuum dispersion relation eq. (4.3) is recovered, indicating a good control over the continuum limit.

Conclusion
This study is motivated by the need to explore new and alternative ways for discretising heavy flavours in simulations of lattice QCD: more independent predictions for heavy flavour hadronic quantities with a solid control of systematic uncertainties are urgently needed [58] to make reliable predictions for SM phenomenology.
To this end we explored the feasibility of Möbius domain wall fermions (MDWF) as a lattice regularisation for heavy quarks. DWF have so far been widely used as a discretisation for the light u, d and s quarks. Its desirable features are chiral symmetry to a good approximation and automatic O(a)-improvement.
From our simulations within quenched QCD with the tree-level improved Symanzik gauge action we have identified a point in MDWF parameter space, the domain wall height M 5 = 1.6, for which discretisation effects turn out to be particularly small. We demonstrated that the salient features of MDWF persist for heavy quarks as long as the bare input quark mass obeys the bound am h 0.4. Based on these findings we carried out a detailed scaling study of the heavy-strange dispersion relation and decay constant. Over the range of lattice cutoffs 2.0-5.7 GeV the observables were compatible with a linear scaling in a 2 . At the level of precision achieved in this work, coefficients of a 4 terms were found to be almost always compatible with zero, remaining remarkably small even for the heaviest quark mass (heavier than charm) simulated.
The results accumulated in this paper constitute a proof of concept for MDWF as a powerful discretisation to study charm and heavier quarks on current dynamical gauge field ensembles. This work constitutes a solid basis for RBC/UKQCD's heavy MDWF phenomenology program [59]. Nevertheless, we are considering ideas for how to improve the current setup, for example by link-smearing the MDWF kernel [12,59].

A Topological charge evolution
In figure 8 we show the Monte Carlo histories and histograms of the topological charge restricted to the configurations on which we also determined the decay constant and the meson energies.