Present and future prospects for lattice QCD calculations of matrix elements for nEDM

TL;DR

This paper reviews the current status and future prospects of lattice QCD calculations for neutron EDM contributions from various CP-violating operators, highlighting progress, challenges, and expected timelines for improved results.

Contribution

It provides a comprehensive status report on lattice QCD calculations of neutron EDM contributions from four CP-violating operators, emphasizing recent progress and future challenges.

Findings

Precise results currently only for quark EDM from PNDME

Initial signals observed for $ heta$-term and chromo EDM contributions

Progress with ~50% uncertainty expected in five years

Abstract

A status report on the calculations of the contribution of four CP violating operators, the $\Theta$-term, the quark EDM, the chromo EDM and the Weinberg operator to the neutron EDM are presented. At this time, there exit precise physical results only for the quark EDM operator by the PNDME collaboration. First results showing signal in the contributions of the $\Theta$-term and the connected part of the chromo EDM operator have been presented. The challenge of divergent mixing in the chromo EDM and Weinberg operators has motivated calculations in the gradient flow scheme. While there has been steady progress, the challenges remaining are large. Results with $O(50\%)$ uncertainty with control over all systematic errors can be expected for the $\Theta$-term over the next five years. Prediction of a timeline for progress on the chromo EDM and the Weinberg operators will depend on when the renormalization and divergent mixing of these operators is brought under control. The most optimistic scenario is that the gradient flow scheme provides a solution to the numerical signal and mixing problems for both the gluonic and quark operators.

Equations

$${\mathcal{L}}_{\mathrm{QCD}}⟶{\mathcal{L}}_{\mathrm{QCD}}^{\cancel{\text{CP}}}={\mathcal{L}}_{\mathrm{QCD}}$$

$${⟨n\mid J_{\mu }^{\mathrm{EM}}\mid n⟩|}_{/\mathrm{CP}}^{\Theta }$$

$${⟨n\mid J_{\mu }^{\mathrm{EM}}\mid n⟩|}_{/\mathrm{CP}}^{\mathrm{qEDM}}$$

$${⟨n\mid J_{\mu }^{\mathrm{EM}}\mid n⟩|}_{/\mathrm{CP}}^{\mathrm{cEDM}}$$

$${⟨n\mid J_{\mu }^{\mathrm{EM}}\mid n⟩|}_{/\mathrm{CP}}^G$$

$$⟨n|J_{\mu }^{\mathrm{EM}}|n{⟩}_{/\mathrm{CP}}=\frac{F_3(q^2)}{2M_n}{\bar{u}}_n{q}_{\nu }{\sigma }^{\mu \nu }{\gamma }_5{u}_n\mathrm{and}{d}_n=\underset{q^2\to 0}{\lim }\frac{F_3(q^2)}{2M_n}$$

$$e\sum _q\overline{q}{\gamma }^{\mu }q⟶e\sum _q\overline{q}{\gamma }^{\mu }q+{ϵ}^{\mu \nu \rho \sigma }{p}^{\nu }\sum _q{d}_q^{\gamma }\overline{q}{\Sigma }^{\rho \sigma }q,$$

$$d_n=d_u^{\gamma }{g}_T^u+d_d^{\gamma }{g}_T^d+d_s^{\gamma }{g}_T^s,\cdots ,$$

$$g_T^u=0.784(28)(10);g_T^d=-0.204(11)(10);g_T^s=-0.0027(16).$$

$$d_N^{\Theta }=a{M}_{\pi }^2+b{M}_{\pi }^2\ln M_{\pi }^2+\cdots .$$

Full text

1 Introduction

One of the deepest mysteries of the observed universe is the matter-antimatter asymmetry. The observed universe has $6.1^{+0.3}_{-0.2}\times 10^{-10}$ baryons for every black body photon [1], whereas in a baryon symmetric universe, we expect no more that about $10^{-20}$ baryons for every photon [2]. It is difficult to include such a large excess of baryons as an initial condition in an inflationary cosmological scenario [3]. The way out of the impasse lies in generating the baryon excess dynamically (baryogenesis, leptogenesis, …) during the evolution of the universe.

In the early history of the universe, if the matter-antimatter symmetry was broken post inflation and reheating, then one is faced with Sakharov’s three necessary conditions [4]: the process has to violate baryon number, evolution has to occur out of equilibrium, and CP (or equivalently time reversal invariance if CPT remains unbroken) has to be violated.

$\cancel{\text{CP}}$  exists in the electroweak sector of the standard model (SM) of particle interactions due to a phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix [5], and possibly by a similar phase in the leptonic sector, given that the neutrinos are not massless [6]. The strength of the $\cancel{\text{CP}}$  in the CKM matrix is much too small to explain baryogenesis. Leptogenesis from the neutrino sector with large lepton-baryon conversion is another possible mechanism, however, no CP violation in the lepton sector has been observed so far.

The SM has an additional source of CP violation arising from the effect of QCD instantons. The presence of these finite action non-perturbative configurations in a non-Abelian theory often leads to inequivalent quantum theories defined over various ‘$\Theta$’-vacua [7]. However, because of asymptotic freedom, all non-perturbative configurations including instantons are strongly suppressed at high temperatures where rates of baryon number violating processes are sizable. Because of this, $\cancel{\text{CP}}$  due to such vacuum effects do not lead to appreciable baryon number production.

In short, the overriding consensus is that additional much larger $\cancel{\text{CP}}$ is needed from physics beyond the SM (BSM). Even though the BSM theory that describes nature above the TeV scale is not known, using the tools of effective field theory one can organize, by symmetry and dimension, possible $\cancel{\text{CP}}$ interactions at the hadronic scale. One then needs to quantify their contribution to the neutron electric dipole moment (nEDM). Each contribution consists of a product of the coupling (BSM model dependent) and the matrix element of the low-energy effective interaction (both defined at the hadronic scale0 within the neutron state (BSM model independent). In this review, I will discuss the status of lattice QCD calculation of the matrix elements of four of the leading, within the effective field theory framework, $\cancel{\text{CP}}$ operators.

2 Nucleon Matrix elements of $\cancel{\text{CP}}$  operators

Over the past few decades, many extensions of the SM have been proposed in the literature. At the hadronic scale (${}\sim 2$ GeV), the effects of BSM scenarios that involve heavy degrees of freedom at the mass scale $\Lambda_{\rm BSM}>M_{W}$ can be described in terms of effective local operators composed of quarks and gluons. Using tools of effective field theory, one can organize all possible effective $\cancel{\text{CP}}$ interactions of quarks and gluons based on symmetry and dimension, and independent of specific BSM theory [8, 9]. In general, operators with higher dimension are suppressed by increasing inverse powers of $\Lambda_{\rm BSM}$ where $\Lambda_{\rm BSM}$ is the scale of new physics. The couplings associated with these operators encode information about the BSM model at the scale $\Lambda_{\rm BSM}\sim$ TeV and the renormalization group evolution from $\Lambda_{\rm BSM}$ to the hadronic scale.

The current goal of lattice QCD calculations is to examine $\cancel{\text{CP}}$  operators with dimension six and lower. Of these, $\cancel{\text{CP}}$  four quark operators of dimension six have not been considered because they are sub-leading in many BSM scenarios and because the lattice methodology to compute their contribution to the neutron EDM has not yet been developed. Focus of the lattice community has been on the following four, which encode the leading $\cancel{\text{CP}}$  effects in a large class of BSM models:

[TABLE]

where the first term is the $\Theta-$interaction and the last is the dimension six three-gluon Weinberg operator. The $\Theta$-term is a part of the SM, but is usually neglected because the coupling $\Theta$ is constrained to be smaller than $10^{-10}$ by the current bound on the nEDM and/or it is assumed that some form of a Peccei-Quinn mechanism tunes $\Theta$ to zero [10]. Note that the $\Theta$-term can be rotated into a pseudoscalar mass term $im_{\ast}({\Theta})\ \overline{q}\gamma_{5}q$ under a chiral transformation [11]. The two dimension five operators are called the quark EDM (qEDM) and the quark chromo-EDM (cEDM). The couplings $d_{u,d,s}^{\gamma}$ are the quark EDMs, the $d_{u,d,s}^{g}$ are the quark chromo-EDMs, and $d_{G}$ is the strength of the Weinberg operator. They are generated directly by threshold effects at the scale $\Lambda_{\rm BSM}$ or by mixing under renormalization group evolution. They parameterize the strength of new CP violating interactions that a given BSM theory generates at the hadronic scale.

In an ideal world, the best way to calculate these matrix elements would be to simulate a lattice theory with these $\cancel{\text{CP}}$ interactions (Eq. (1)) added to say the Wilson-clover fermion action. This ideal approach does not work because these interactions are complex and we do not yet know how to efficiently simulate theories with a complex action. The approach that works for lattice QCD is to treat the small $d_{q}^{\gamma,g}$, $\Theta$ and $d_{G}$ as perturbations and expand the theory about the normal CP conserving action, such as the Wilson-clover action. Then, to lowest order in $\alpha_{em}$, the lattice calculation involves the product of the electromagnetic current $J^{\rm EM}_{\mu}$ and each of these operators:

[TABLE]

The qEDM is an exception as the leading contribution to it arises from the modification to $J^{\rm EM}_{\mu}$ as discussed in Sec. 4. These matrix elements of $J^{\rm EM}_{\mu}$ between neutron states in the presence of CP violation are model independent and provide the “connection” between the couplings and the nEDM as exemplified in Eq. (8). Note that, since each $\cancel{\text{CP}}$ interaction contributes to the nEDM, the value of (bound on) the nEDM provides a single constraint on the sum of all the contributions. Nevertheless, lowering the bound on the nEDM will provide increasingly tight, and amongst the most stringent constraint, on possible BSM models.

3 Challenges to the calculation

The contributions of the $\Theta$, cEDM and Weinberg terms to the nEDM are, in most lattice calculations, obtained from the $P$ and $T$ violating form factor $F_{3}$ that arises in the decomposition of the matrix elements defined in Eqs. (2), (4), and (5). For example, for the cEDM operator and ignoring the standard form factors $F_{1}$ and $F_{2}$ in the first term for brevity, one finds

[TABLE]

In practice, in the extraction of $F_{3}$ from the matrix element, one has to take into account mixing of $F_{3}$ with the $CP$ conserving form factor $F_{2}$ due to a phase $e^{i\alpha_{N}\gamma_{5}}$ that $\cancel{\text{CP}}$ interaction generates in the neutron interpolating operator. Note that there is a different $\alpha_{N}$ generated for each of the $\cancel{\text{CP}}$ interactions. This subtlety has been discussed and resolved in Ref. [12]. To keep this review brief, the reader is referred to the original papers cited for details. Here, I assume that the reader is familiar with the correct phase convention, the calculation of this phase $\alpha_{N}$ from the nucleon 2-point function, and its impact on the extraction of $F_{3}(Q^{2})$.

There are two very important challenges to getting results at the physical point for the contribution of the $\Theta$, cEDM and Weinberg operators to the nEDM:

• •

The first important challenge is the signal in $F_{3}$ is very small. In fact, in most calculations, it is barely significant. The matrix elements have to be calculate for non-zero momentum transfer, excited-state contamination (ESC) removed, $F_{3}$ extracted, mixing with $F_{2}$ resolved, and then extrapolated to $Q^{2}=0$. Each of these steps is non-trivial. In Secs. 5 and 6, I describe newly developed variance reduction methods that reduce the errors by almost a factor of ten.

• •

Even after a signal has been demonstrated, defining the renormalized cEDM and Weinberg operators and obtaining finite results in the continuum limit is non-trivial. The reason is they mix with the same and lower dimension operators under renormalization. The 1-loop analysis of the cEDM operator [13] shows that it mixes with the pseudoscalar $i\overline{q}\gamma_{5}q$ and the $\Theta$ operators. Of these, the mixing with $i\overline{q}\gamma_{5}q$ is quadratically divergent in all lattice formulations, i.e., both chiral symmetry preserving and violating. As a result, the mixing coefficients have to be determined very precisely nonperturbatively. In fermion formulations such as Wilson-clover that explicitly break chiral symmetry, there is an additional divergent mixing with $im\overline{q}\gamma_{5}q$. The mixing pattern of the Weinberg operator is even more complicated.

4 Quark EDM operator and its matrix element

In the presence of $\cancel{\text{CP}}$ interactions, the electromagnetic current, defined as $\delta{\cal L}/\delta A^{\mu}$ with $\cal L$ given in Eq. (1), gets an additional term

[TABLE]

which, at $\vec{p}=0$ is the quark bilinear operator with tensor structure. There is also a term analogous to Eq. (4) that is generated. It is, however, suppressed by $\alpha_{\rm em}$ and therefore not considered. Thus, to leading order the contribution to neutron EDM is given by

[TABLE]

where $g_{T}^{u,d,s}$ are the flavor diagonal tensor charges given by the matrix elements $\langle N|\overline{q}\sigma^{\mu\nu}q|N\rangle$ and $d_{u,d,s}^{\gamma}$ are the corresponding BSM couplings.

The lattice methodology for the calculation of $g_{T}^{u,d,s}$ is robust and reliable results at the physical pion mass $M_{\pi}=135$ MeV and in the continuum limit have been presented in Ref. [14]. A brief description of the methodology of lattice QCD calculations of nucleon 3-point functions for the analogous axial operator is given in a companion paper in the proceedings of this conference, PoS(Spin2018)018. It involves calculating the three correlations functions illustrated in Fig. 1, with the qEDM (tensor) operator inserted at time $t$ in the connected and disconnected 3-point diagrams.

Lattice calculations of nucleon charges have been reviewed in the FLAG 2019 report [15]. Data for the tensor charges exhibit small discretization or finite volume corrections and have been stable over time. The current best estimates, in the $\overline{MS}$ scheme at 2 GeV, are

[TABLE]

Both the connected and disconnected contributions were obtained at the physical point by fitting data at multiple values of $a$ and $M_{\pi}$ and removing the leading continuum-chiral corrections. In short, flavor-diagonal tensor charges have been calculated with control over all systematics [14, 15].

Using these results, the authors of Ref. [14] analyze constraints on the split SUSY model [16, 17, 18]. This model is pertinent because in it qEDM is the dominant $\cancel{\text{CP}}$  operator. Using the experimental bounds $d_{n}\leq 2.9\texttimes 10^{−26}$ e cm [19] and $d_{e}\leq 1.1\texttimes 10^{−29}$ e cm [20], gave the upper bound $d_{n}\leq 4.1\texttimes 10^{−29}$ e cm for the split-SUSY model [14]. More BSM theories can be analyzed as results for other $\cancel{\text{CP}}$ operators become available.

5 Calculation of the $\Theta$-term

The $\Theta$-term breaks $P$ and $T$ invarinace and thus $CP$ by the CPT theorem. In the absence of a Peccei-Quinn mechanism, the $\Theta$-term arises naturally in the SM. Operator mixing under renormalization group flow of the cEDM and Weinberg operators between $\Lambda_{\rm BSM}$ and the hadronic scale also generates it, i.e., the cEDM and Weinberg operators mix with the $\Theta$-term under renormalization. Therefore, calculating its contribution to nEDM is essential. Attributing the cause of a non-zero nEDM to an intrinsic $\Theta$-term versus one generated from BSM interactions will be challenging and will require knowing the matrix elements of all [leading] $\cancel{\text{CP}}$ operators.

A number of calculations of the contribution of the $\Theta$-term to nEDM using the $F_{3}$ form factor method have been done [21, 22, 23, 24, 25, 12, 26, 27]. Of these, calculations done prior to Ref. [12] have poor [no] statistical signal and did not properly include the phase induced by $\cancel{\text{CP}}$ operators in the nucleon state as pointed out in Ref. [12].

The calculation involves the correlation of the purely gluonic operator $\int d^{4}xG^{\mu\nu}\tilde{G}^{\mu\nu}$ (topological charge) with the nucleon 3-point function $\langle N(0)J^{\rm EM}_{\mu}(t)N(\tau)\rangle$ as shown in Fig. 2. Depending on the lattice generation methodology, the value of the gluonic term, the topological charge, can exhibit very long time auto-correlations. The signal in the fermionic part is very good and typically $\langle N(0)J^{\rm EM}_{4}(t)N(\tau)\rangle$ is used. The correlation between the two is, therefore, the fermionic part weighted by the topological charge. Consequently, if the topological charge is frozen during lattice generation, then the correlation and the projection on to $F_{3}$ will have a poor/biased signal.

Preliminary analysis with evidence of a signal has been presented in Ref [26]. This study used only one 2+1 flavor domain wall ensemble with $a=0.1105$ fm and $M_{\pi}=340$ MeV; calculated only the connected 3-point nucleon correlation function; and did not remove the ESC. The topological charge was defined in the gradient flow scheme to reduce noise in it. Most of the effort has been devoted to developing/testing the following variance redution method to get a signal. In the correlation of $\int d^{4}xG^{\mu\nu}\tilde{G}^{\mu\nu}$ with the nucleon 3-point function, the authors propose to use a 4-d volume centered about the nucleon correlator over which to sum $G^{\mu\nu}\tilde{G}^{\mu\nu}$ rather than the whole lattice which gives the topological charge. This setup is illustrated in Fig. 3. The motivation is that $G^{\mu\nu}\tilde{G}^{\mu\nu}$ on points outside this volume contribute only noise. While there is some evidence for variance reduction as a result, the concomitant issue of introducing a possible bias has not been settled. One could apply a variant of the standard bias correction method [28, 29] by calculating the correlation with both the full and subvolume sum for a few of the nucleon 3-point functions on each configuration. Such a bias correction method has not yet been explored.

Furthermore, the authors contend that a reliable signal at the physical pion mass will need a new level of precision or alternate methods. Chiral perturbation theory predicts that the contribution of the $\Theta$-term to $d_{n}$ vanishes in the chiral limit as [30, 31, 32]

[TABLE]

Thus, as one tunes $M_{\pi}\to 135$ MeV in lattice simulations, the precision of the calculation will have to be increased significantly to keep the fractional error the same. For example, on going from $340\to 135$ MeV, the value of $d_{N}^{\Theta}\propto M_{\pi}^{2}$ is expected to decrease by a factor of about six. Thus, the statistics will have to be increased by $O(100)$ to keep the fractional error the same [26].

Results from six ensembles with $M_{\pi}\geq 410$ MeV have just been published in Ref. [27]. The topological charge density is again calculated in the gradient flow scheme. The central values are all negative, lie between $-0.0021$ and $-0.0070$, but five of them differ by less than 2$\sigma$ from zero. There is no evidence of a $d_{N}^{\Theta}\propto M_{\pi}^{2}$ behavior at these heavy pion masses, and a less than convincing chiral-continuum fit gave $d_{n}=-0.00186(59)\ \Theta$ e fm. Clearly, far more precise calculations near the physical pion mass are needed before a result can be presented with confidence.

My bottom line conclusion on the contribution of the $\Theta$-term is that, while, the methodology for the calculation of $F_{3}$ from the $\Theta$-term is now established and there are no show stopping issues of renormalization, more work needs to be done to demonstrate a $5\sigma$ signal on $M_{\pi}\gtrsim 350$ MeV ensembles, and $O(100)$ more to obtain results at $M_{\pi}=135$ MeV and in the $a\to 0$ limit.

6 Enhancing the signal in $F_{3}$ from the cEDM operator

Calculations of the cEDM operator need to address both the signal and the renormalization, including the divergent mixing with the pseudoscalar operator $i\overline{q}\gamma_{5}q$. The two calculations reported in Refs. [33, 26] are both based on extracting $d_{n}$ from $F_{3}$. So far, the focus in these works has been on getting a signal, and that too in only the connected part. No results have been reported for the mixing coefficients calculated non-perturbatively in either of the two approaches.

The work described in Ref. [33] uses the Schwinger source method to include the $\cancel{\text{CP}}$ interaction in the Dirac action. The authors show that the phase $\alpha_{N}$ associated with the neutron ground state and generated by the $\cancel{\text{CP}}$ interaction can be extracted reliably from the 2-point function. To obtain a signal in $F_{3}$, which $a\ priori$ was poor, they developed a general variance reduction method. It exploits correlations between an observable $O$ and a number of other quantities $R_{i}$ with $\langle R_{i}\rangle=0$. For the connected part of cEDM, the variance of the combination $O+\xi_{i}R_{i}$ was shown to be reduced by almost a factor of ten by a suitable choice of $R_{i}$ and optimizing the parameters $\xi_{i}$.

The extraction of the three form factors, $F_{1}$, $F_{2}$ and $F_{3}$, from the 8 matrix elements (real and imaginary parts of the four components of the vector current) forms an over complete set. The authors find that the results for the unrenormalized connected contributions to $F_{3}$ from different combinations give different estimates. These differences point to possible large excited state and/or discretization effects that need to be resolved. The real hurdle in these calculations, in schemes such as RI-sMOM, is controlling the divergent mixing with the lower dimension pseudoscalar operator. Work has, therefore, been initiated to do the calculation in the gradient flow scheme.

Results for the unrenormalized connected contributions to both the cEDM and pseudoscalar operator it mixes with using the direct 4-point method have been presented in Ref. [26]. Having demonstrated a signal in $F_{3}$ for both the cEDM and pseudoscalar operator, they are now investigating the position-space renormalization scheme to control the mixing problem, and on including the disconnected contributions.

7 The Weinberg operator

The lattice calculation of the Weinberg operator is similar to that of the $\Theta$-term except for the additional serious complication of a divergent mixing with the $\Theta$-term under renormalization. The method currently being explored that would control this mixing on the lattice is gradient flow [34, 33]. Fig. 4 shows a comparison of the susceptibilities of the topological charge and the Weinberg operator as a function of the flow time. The data for the topological susceptibility quickly flattens out as expected while that for the scale-dependent Weinberg susceptibility continues to evolve. The next step is to calculate the 4-point correlation function as a function of the flow time. Simultaneously, our extended collaboration at LANL is working to relate the gradient flow scheme to the continuum $\overline{MS}$ scheme. With this matching in hand, we will need to demonstrate that there exists a window in flow time in which results in the continuum are independent of the flow time. In short, a number of developments have to pan out before physical results can be obtained.

Fig. 5 shows data for the Weinberg operator and the topological charge versus Monte Carlo configuration number on $a12m310$ ensemble at flow time $t_{\rm WF}=5$. The two are highly correlated. So, any technique for improving the signal in the $\Theta$-term will also be applicable to the Weinberg operator. It also highlights a large divergent mixing with the $\Theta$-term that needs to be addressed.

8 Conclusions

The prospect of reducing the upper bound on the nEDM from current and future experiments, and possibly finding a value is exciting. It will signal $T$ (and $CP$ assuming CPT) violation larger than in the standard model, and will put stringent constraints on BSM theories provided the matrix elements of $\cancel{\text{CP}}$ operators can be calculated reliably. A number of groups are using large scale simulations of lattice QCD to calculate the matrix elements of the four operators reviewed, and their contributions to the nEDM. These calculations are extremely hard except for the quark EDM operator. Clearly new ideas and algorithms are needed! My summary of the status and prospects for obtaining results of dimension six and less operators, in order of the most likely within the next five years, is the following:

• •

The quark EDM operator. The leading contribution of this operator has been calculated and results with $O(5\%)$ errors, after extrapolation to the continuum limit and at the physical pion mass, have been obtained as discussed in Sec. 4.

• •

The $\Theta$-term: Calculations of its contribution to nEDM have the longest history. Evidence of a signal is growing, and I expect validated results with $O(50\%)$ uncertainty to be available over the next five years.

• •

The chromo EDM operator: Methods to get a $5\sigma$ signal are being developed and tested. Results including both the connected and disconnected contributions at multiple values of the lattice spacing require $O(100)$ or more in computing resources. Thereafter, tests and control of the divergent mixing with the pseudoscalar operator can begin.

• •

The Weinberg operator: The mixing problem under renormalization is severe. Working in the gradient flow scheme is the approach of choice being investigated [34, 33]. First tests of the methodology and the numerical signal are being performed. For both the Weinberg and chromo EDM operators, estimates with $O(1)$ uncertainty are possible in five years if the gradient flow method works.

• •

4-fermion operators: So far there is no published work on a renormalization framework for these operators that will be suitable for lattice calculations nor have exploratory lattice calculations begun. We are unlikely to see significant calculations over the next five years, if for no other reason than due to the limited access to computer time. Most groups will likely channel focus and resources to the previous three operators first.

Acknowledgments.

I thank Andreas Wirzba, Paolo Lenisa and the organizers of Spin 2018 for inviting me to give this review and their hospitality. I thank Tanmoy Bhattacharya, Vincenzo Cirigliano, Andrew Kobach, Emanuelle Mereghetti and Boram Yoon, for a fun and fruitful collaboration in which we have seen the statistical signal come and go, and the methods evolve in response. We thank the MILC collaboration for sharing the $2+1+1$-flavor HISQ ensembles generated by them. We gratefully acknowledge the computing facilities at, and resources provided by, NERSC, OLCF at Oak Ridge, USQCD and LANL Institutional Computing.