Ruling out the massless up-quark solution to the strong CP problem by computing the topological mass contribution with lattice QCD

The infamous strong CP problem in particle physics can in principle be solved by a massless up quark. In particular, it was hypothesised that topological effects could substantially contribute to the observed non-zero up-quark mass without reintroducing CP violation. Alternatively to previous work using fits to chiral perturbation theory, in this letter we bound the strength of the topological mass contribution with direct lattice QCD simulations by computing the dependence of the pion mass on the dynamical strange-quark mass. We find that the size of the topological mass contribution is inconsistent with the massless up-quark solution to the strong CP problem.

Introduction.-One of the unsolved puzzles in particle physics is the so-called strong CP problem, where CP stands for the combined charge conjugation and parity symmetry. In quantum chromodynamics (QCD), the theory of strong interactions, the non-trivial topological vacuum structure generates a CP-violating term ∝ θ G µνG µν in the Lagrangian, where θ is an a priori unknown parameter, G is the gluon field strength tensor andG its dual. However, experimentally there is no sign of CP violation in QCD. Instead, the strong upper bound θ < 10 −10 [1] from measurements of the neutron electric dipole moment leads to a severe fine tuning problem.
There are several proposals to overcome this problem, for instance by postulating the existence of an axion [2][3][4]. A simple alternative could be the vanishing of the up-quark mass m u , which at first sight seems inconsistent with results of current algebra. However, Refs. [5][6][7][8] pointed out that the up-quark mass in the chiral Lagrangian has two different contributions: a CP-violating perturbative contribution m u and a CP-conserving nonperturbative contribution m eff from topological effects, such as instantons. While m u = 0 could be easily ensured by an accidental symmetry [8][9][10][11][12], m eff does not contribute to the neutron electric dipole moment and is parametrically of order m eff ∼ m d m s /Λ QCD , plausibly as large as the total required up-quark mass. Testing this simple solution to the strong CP problem is particularly important because the other proposed solutions, including the axion [2][3][4] and Nelson-Barr [13,14] mechanisms, face several theoretical challenges [15].
As the only tool to reliably test the m u = 0 proposal [8], lattice gauge theory has determined the up-quark mass to m u (2 GeV) ∼ 2 MeV by fitting the light meson spectrum with errors around 5% (see Ref. [16] for a review). As proposed in Refs. [17,18], it would be beneficial to perform a complementary analysis by calculating the dependence of the pion mass on the dynamical strange-quark mass while keeping the light quark masses fixed. This direct calculation would have the advantage of avoiding any fitting procedures.
Lattice QCD simulations are now being performed taking into account the first two quark generations as dynamical degrees of freedom. In addition, simulations are performed at (or very close to) the physical values of the pion, kaon, and D-meson masses [19] and at various values of the lattice spacing and volumes, such that systematic effects can be studied and eventually controlled [20,21]. Finally, the theoretically sound definitions of the topological charge and susceptibility on the lattice (see Ref. [22] for a review) allow for directly accessing topological effects related to m eff .
In this paper, we perform a cross-check of the m u > 0 hypothesis based on the proposals of Refs. [17,18]. In particular, we compute the parameter β 2 /β 1 which measures the strength of m eff and probes the contribution of small instantons and other topological effects to the chiral Lagrangian. While β 2 /β 1 is usually obtained from a combination of low-energy constants [23], this indirect lattice method requires chiral perturbation theory (χPT). Using direct lattice computations instead, we obtain the result β 2 /β 1 = 0.63(39) GeV −1 by computing the dependence of the pion mass on the strange-quark mass. Since a bound significantly smaller than 5 GeV −1 provides an exclusion of the massless up-quark hypothesis [8,18], our result rules out this hypothesis, in accordance with previous fits of chiral perturbation theory to lattice data [16].
Method.-We test the m u = 0 proposal by investigating the variation of the pion mass with respect to the strange-quark mass. The general form of the quark-mass dependence of the pion mass reads [24] where the first term is the first-order contribution of the light quark masses in χPT. The second term receives contributions both from small instantons that could mimic a non-zero m u and from higher-order terms in χPT that are proportional to m s , including logarithmic corrections. In order to let topological effects explain the observed value for m u and to allow for a solution of the strong CP problem, β 2 /β 1 ≈ 5 GeV −1 at renormalization scalē µ = 2 GeV in the MS scheme is required [18].
The most precise and computationally challenging test of the ratio β 2 /β 1 is to vary either the strange quark mass or the light quark mass, m u = m d ≡ m . For example, by varying m s while keeping m fixed, we obtain [17,18] where M π,i = M π (m s,i ) is the average pion mass as a function of the varied strange-quark mass m s,i at fixed m . Note that the approximate result for β 2 /β 1 in Eq. (2) is independent of the up/down quark masses. Crucially, this allows us to reliably compute β 2 /β 1 even at larger than physical quark masses. The higher-order corrections in Eq. (1) reintroduce a small residual pion-mass dependence for β 2 /β 1 that finally needs to be cancelled by a chiral extrapolation. While this challenging direct method to compute the ratio β 2 /β 1 is independent of χPT, the more common indirect method is to use the chiral Lagrangian. For example, Ref. [18] used lattice data from the FLAG report of 2013 [25] to estimate β 2 /β 1 (1±1) GeV −1 neglecting chiral logarithms and higher order terms in the chiral Lagrangian. To check the consistency of our computations with the results of Ref. [18], we have also computed β 2 /β 1 indirectly by using chiral fits and measuring M 2 K (m s ), obtaining excellent agreement with Ref. [18].
Lattice Computation.-In this paper we use gauge configurations generated by the Extended Twisted Mass collaboration (ETMC) with the Iwasaki gauge action [26] and Wilson twisted mass fermions at maximal twist [27,28] with up, down, strange and charm dynamical quark flavours. Up and down quarks are mass degenerate. All the gauge configuration ensembles we used are listed together with the corresponding pion and strange quark mass values in Table I. For details on how these values are obtained, we refer to the supplemental material [29].
We first perform the analysis using three sets each with a pair of ensembles (AX and AXs with X = 60, 80, 100) without the so-called clover term in the action. Details on the production of these ensembles can be found in Ref. [30]. Each pair with X = 60, 80, 100 has identical parameters apart from strange and charm quark mass values, which are close to their physical values. The three pairs have equal strange and charm quark masses within errors but differ in the light quark mass value corresponding to unphysically large pion mass values of about 386 MeV, 444 MeV and 494 MeV, respectively. The lattice spacing value corresponds to a = 0.0885(36) fm [31] determined from the pion decay constant f π . In addition, we use one ensemble (cA211.30.32) that includes the clover term in the action [21]. While Wilson twisted mass fermions at maximal twist automatically remove discretization effects linear in the lattice spacing a [32] and thus leave only lattice artifacts at O(a 2 ), the clover term reduces these O(a 2 ) effects even further [33]. The cA211.30.32 ensemble has a smaller pion mass value of about 270 MeV and strange and charm quark mass values again close to their physical values. The lattice spacing value is a = 0.0896(10) fm determined using the nucleon mass dependence on the pion mass. This estimation is done by employing χPT at O(p 3 ) [34,35], where p is a typical meson momentum. Similarly to Ref. [21] the N f = 2 + 1 + 1 non-clover twisted mass ensembles [36] at different lattice spacings, which also include the AX ensembles, were used to control the chiral extrapolations.
Since the pion mass of the cA211.30.32 ensemble is significantly smaller than the ones of the AX ensembles and thus closer to the physical value, we consider this ensemble as the most appropriate one to compute our final value of β 2 /β 1 . Moreover, cA211.30.32 uses the same action as ensembles that are currently under production with a physical value of the pion mass. These ensembles could be used, in principle, in future work to repeat the calculation presented here at the physical point. For cA211.30.32 we have simulations for only one dynamical strange quark mass; thus it is necessary to apply the so-called reweighting technique to investigate the strange quark mass dependence of M π , while keeping the charm and light quark masses constant [29]. We denote with cA211.30.32l (cA211.30.32h) the reweighted ensem- ble with 5% lower (higher) strange quark mass value than the original ensemble cA211.30.32. In contrast, for the AX(s) ensembles we have pairs of ensembles with different dynamical strange quark masses; thus we can use a direct approach to investigate the strange quark mass dependence of M π . Note that in this case also the charm quark mass differs slightly, but its value is so close to the cut-off that this difference will not affect our results. While the AX(s) ensembles have rather heavy pion masses (see Table I), they are ideal to test the robustness of the reweighting procedure that we apply to cA211.30.32. In fact, we use these ensembles to demonstrate that reweighting works successfully [29]. In addition, the β 2 /β 1 values from these ensembles provide an insight into the pion mass dependence of β 2 /β 1 .
Results.-Using the values of M π and m s from Table I as input, we compute β 2 /β 1 from Eq. (2). The results from this direct approach for the three pairs A60(s), A80(s) and A100(s) are compiled in Table II, where we quote besides β 2 /β 1 also the numerator and denominator of Eq. (2) separately. Since the pion mass differences are all zero within errors (see Table I), we find with this approach also β 2 /β 1 to be compatible with zero. Note that the errors of the observables compiled in Table I are correlated per ensemble. This correlation is taken into account in our analysis for β 2 /β 1 .
Finally, we use reweighting on the cA211.30.32 ensemble to vary the strange quark mass by ±5% around its original value. The change in the pion mass with the strange quark mass is not significant, see Table I. The corresponding values for β 2 /β 1 are again compiled in Table II In Fig. 1 we show the values of the ratio β 2 /β 1 at µ = 2 GeV in the MS-scheme as a function of the squared pion mass M 2 π , both in physical units. The three blue points at heavier pion mass values correspond to the three pairs of the AX(s) ensembles without clover term. The three red points at lower pion mass value correspond to the cA211.30.32 ensemble including the clover term. The The ratio β2/β1 as a function of the squared pion mass M 2 π , both in physical units. The solid line with the 1σ error band represents a linear extrapolation in M 2 π . We extrapolate to the chiral limit to eliminate higher-order corrections to β2/β1, see Eqs. (1) and (2).
latter three points are slightly displaced horizontally for better legibility. While all of the points are compatible with zero at the 1.5σ-level, we observe a slight trend towards larger β 2 /β 1 -values with decreasing pion mass values.
In addition, we show in Fig. 1 a linear extrapolation of β 2 /β 1 in M 2 π to the chiral limit. This linear dependence can be justified with χPT, which predicts [16] modulo logarithmic corrections, where α 1,2,3 are combinations of low-energy constants with Since the data points for the ensemble cA211.30.32 are highly correlated, we include only the combination of cA211.30.32h and cA211.30.32l in the fit denoted as cA211.30.32(h,l) in Table II. The fit has χ 2 /dof = 3.28/2, i.e., a p-value of 0.2 and the chirally extrapolated value reads β 2 /β 1 = 0.63(25) GeV −1 . As mentioned, we extrapolate to the chiral limit to cancel the residual pionmass dependence in Eq. (3), which stems from higherorder corrections in Eq. (1) and does not appear in the expression for β 2 /β 1 in Eq. (2). Our data thus confirms in hindsight that the approximation in Eq. (2) is justified.
Discussion.-All the estimates for the ratio β 2 /β 1 presented in this letter are consistent with zero at the 1.5σ level. With the chiral extrapolation explained above and 1σ statistical uncertainty, we exclude a value of 5 GeV −1 by an amount significantly larger than 10σ. The remaining question is whether there are additional systematic uncertainties that could potentially spoil this conclusion.
Let us first consider the discretization errors for β 2 /β 1 , which are of order (aΛ QCD ) 2 multiplied by an unknown coefficient, with Λ QCD = 341(12) MeV [37]. We can reliably estimate the coefficient by using the known continuum extrapolation values for M 2 π and m s for the AX ensembles [31]. By comparing these continuum values to our lattice results for M 2 π and m s , we can infer the size of discretization errors at our given lattice spacing. Depending on the scaling variable, the discretization errors in M 2 π and the strange quark mass are both of order 5 − 10%. If propagated generously, this implies a 10% uncertainty on the numerator, a 15% uncertainty on the denominator and thus about 20% on the ratio β 2 /β 1 . Note that this estimate is highly conservative, because most of the discretization effects cancel in the differences in both the numerator and the denominator. Due to the reduced lattice artifacts with the action including the clover term (see supplementary material [29]), we do not expect larger uncertainties on the ratio for the ensemble cA211.30.32 stemming from discretization effects.
In addition, there is a residual pion mass dependence of β 2 /β 1 , which we account for by extrapolating to the chiral limit. In this extrapolation, the errors stemming from different lattice artifacts of the AX and cA211.30.32 ensembles are taken into account by the above-mentioned 20% uncertainty. Lastly, there are finite-size effects for M π proportional to exp(−M π L), with L the spatial extend of the lattice, but no finite-size corrections to m s . Since the strange quark mass dependence of M π is so weak, these finite-size effects are equal for M 2 π,1 and M 2 π,2 and thus cancel in the ratio β 2 /β 1 . In summary, taking the chirally extrapolated value for β 2 /β 1 plus the 1σ statistical error and the 20% uncertainty for discretization effects, we arrive at the following conservative estimate: atμ = 2 GeV in the MS scheme. For the final estimate we have added the errors linearly. Note that our data is equally well compatible with a constant extrapolation in M 2 π , which would lead to significantly smaller value at the physical point. Thus, we consider Eq. (4) as a conservative estimate. Moreover, the logarithmic corrections from chiral perturbation theory contributing to β 2 /β 1 are of the same order as our value in Eq. (4), therefore the topological contribution to β 2 /β 1 should be even smaller.
Conclusion.-In this letter, we have tested the massless up-quark solution to the strong CP problem by directly investigating the strange quark mass dependence of the pion mass on the lattice. This allows us to determine the ratio β 2 /β 1 , which would need to be larger than 5 GeV −1 to solve the strong CP problem.
Since all our estimates of β 2 /β 1 are compatible with zero, we obtain a strong upper bound for β 2 /β 1 including residual uncertainties stemming from discretization errors and chiral extrapolation. The result in Eq. (4) is clearly incompatible with the massless up-quark solution to the strong CP problem. This exclusion of the m u = 0 solution is consistent with previous results using χPT and direct fits of the light meson spectrum.
Given our conservative error estimates, we consider it highly unlikely that the factor five needed to rescue the solution to the strong CP problem is hidden in the quoted uncertainties. A confirmation of this result using ensembles with physical pion mass values could be undertaken in the future, once different values for the lattice spacing become available for a continuum extrapolation.
Our direct lattice results also quantitatively support the large-N picture as a good description of QCD at low scales, because the coefficient of the non-perturbatively induced mass operator is known to be suppressed in the large-N limit [8,17]. Thus, our computations reliably demonstrate that the topological vacuum contributions to the chiral Lagrangian are negligible.
Supplemental Material: Applying mass reweighting to the strange quark mass keeping maximal twist We discuss here our procedure to determine the strange quark mass and explain the reweighting approach used for the cA211.30.32 ensemble to obtain the probability distribution at different values of the strange quark mass. In Table III we list the bare parameters of the gauge ensembles used in this work. They consist of N f = 2 + 1 + 1 twisted mass ensembles without a clover term for the AX(s) ensembles and with a clover term [49] for the cA211.30.32 ensemble. The ensembles without a clover term come in pairs having the same light quark mass, determined by the parameter µ , and different strange quark masses, related to the parameters µ σ and µ δ . They are denoted as A60 and A60s for the lightest pion mass, A80 and A80s for the intermediate and A100 and A100s for the heaviest pion mass. They all have the same lattice spacing a = 0.0885 (36) fm. For the ensemble with a clover term, denoted as cA211.30.30, and with the lattice spacing a = 0.0896(10), we use re-weighting in the strange quark mass of ±5%. We denote the resulting lighter and heavier mass ensembles as cA211.30.30l and cA211.30.32h, respectively.

Strange Quark Mass Determination
In order to evaluate the ratio β 2 /β 1 , we need to determine the renormalised strange quark mass. Two different approaches are used that provide a cross-check. The first approach uses the following relation of the bare parameters µ σ and µ δ of Table III to the renormalized strange and charm quark masses: Here, Z P and Z S are the pseudoscalar and scalar renormalization functions. For the ensembles without a clover term, they have the following values [31]: Z P = 0.529(07) and Z S = 0.747 (12) .
The renormalization functions are given atμ = 2 GeV in the MS-scheme. Knowing the values of Z P and Z S , we can directly compute m a s,c from the bare parameters µ σ and µ δ .
An alternative approach is to use the kaon mass M K computed in two ways: i) we compute M K within the same fermion discretization using the same bare strange quark mass for the valence as for the one in the sea. We denote this kaon mass by M unitary K since the valence and  (14) TABLE III. Parameters of the ensembles used in this work. β = 6/g 2 0 is the inverse squared gauge coupling, µ the bare light quark mass, and µσ and µ δ are parameters related to the renormalized strange and charm quark masses. All dimensionful quantities are given in units of the lattice spacing a. For the cA211 ensembles, the value of the clover improvement coefficient is csw = 1.74.
the sea quark masses are the same; ii) we compute the kaon mass using Osterwalder-Seiler (OS) valence strange quarks [51] with bare strange quark massμ s . The value ofμ s is then adjusted such that the resulting kaon mass The desired value of the bare strange quark mass µ s is given by yielding the renormalised strange quark mass where Z P is the same pseudoscalar renormalization constant as the one used in the first approach.
In order to carry out the matching as described above, we compute M OS K for three values of aμ s and interpolate (aM OS K ) 2 linearly in aμ s . The unitary kaon masses are taken from Ref. [52] or computed when not available. For the ensembles without the clover term we use aμ s = 0.017, 0.019, 0.0225 and for the cA211.30.32 ensemble we use aμ s = 0.0176, 0.0220, 0.0264.
In Table IV, we give β 2 /β 1 computed using m a s and m b s . As can be seen, the mean values agree very well, albeit with large statistical errors. This indicates that discretization artifacts largely cancel in the ratio. We thus conclude that either approach can be utilized to fix the strange quark mass and proceed with m b s to compute the ratio β 2 /β 1 .

Reweighting in the strange quark mass
For the cA211.30.32 ensemble, we apply mass reweighting [53,54] to obtain the probability distribution at three   different values of the strange quark mass. This works if the mass shift is small compared to the original strange quark mass m s . Since the two above-mentioned approaches yield consistent values for β 2 /β 1 and the definition of m a s in Eq. (5) is technically easier to use for the re-weighting, we use the m a s definition to calculate the resulting shifts in the heavy quark parameters. To keep the notation tidy, we drop the index a from m a s and set the lattice spacing to unity in what follows.
In the Boltzmann weight W = e −Sg D f , a shift in the strange quark mass only affects the fermionic part where D and D ND are the light and non-degenerate twisted mass Wilson Dirac operators, respectively. The untwisted bare quark massm is an input parameter which receives an additive renormalization factor m cr and is, when subtracted by m cr , proportional to the current quark mass: m PCAC ∝m − m cr . The current quark mass can be determined via a suitable ratio of matrix elements via the partially-conserved axial current (PCAC) relation. It should be noted that m cr implicitly depends on all other bare parameters in the theory. In the twisted mass formulation, the condition m PCAC = 0 is referred to as maximal twist and results in automatic O(a)-improvement of all physical observables. Since we are using the N f = 2 + 1 + 1 twisted mass action, the mass-degenerate light quark part of D f is det D 2 (m, µ ) = D † (m)D (m) + µ 2 and depends onm, which acts as a constant shift to the diagonal of D , and on the light twisted mass µ , which acts as a twist in spin space via γ 5 to the diagonal of D . The non-degenerate heavy quark part is det D ND (m, m s , m c ) and depends oñ m and, of course, the strange and charm quark masses, m s and m c . The latter are functions of the bare quark mass parameters µ σ and µ δ as given in Eq. (5).
Our results presented in the main manuscript are obtained by performing a reweighting in m s keeping m c constant. From Eq. (5) one can rewrite µ σ = Z P m c − (Z P /Z S )µ δ . This means that a change of the strange quark mass m s keeping m c constant requires knowledge of the ratio of the pseudoscalar to scalar renormalization functions, Z P /Z S . However, a change in m s can result in a significant change in m cr , requiring a corresponding adjustment ofm to maintain maximal twist and an absence of discretization artifacts linear in the lattice spacing (see section 3.1 of Ref. [21] for more details). The AX, AXs and cA211.30.32 ensembles have been simulated at maximal twist while the reweighting procedure to obtain cA211.30.32(l,h) also takes into account the necessary shifts inm as the strange quark mass is varied.
We first confirm that our reweighting for cA211.30.32 works correctly by performing a corresponding reweighting between the parameters of the A60 and A60s ensembles with M π ∼ 400 MeV. We vary µ δ → µ δ + δµ δ = µ δ keeping µ σ constant. Note that both m s and m c change in this test case and we account for the change in m cr by suitably reweighting alsom →m + δm =m . Reweighting proceeds by correcting the Boltzmann weight with the factor where and .
(11) We show in Fig. 2 the results of the reweighting, starting from the ensemble A60. Since we have performed the reweighting only on a subset of the gauge configurations of A60, we include for comparison the value of the pion mass extracted when using the full ensemble. As expected, the two values are in agreement and the comparison is aimed at having an idea of the error using the subset. The determinants for a ratio of matrices A ∈ C V ×V with a positive definite partner matrix A † +A can be computed via a stochastic estimation of the integral representation of the determinant [54] given by with complex Gaussian distributed fields χ i , the total number of stochastic estimates N st and with a normalized integral measure Dη = i∈V dRe(η i )dIm(η i )/π, where η i is the ith entry of the complex vector η ∈ C V ×1 . Now, the reweighting factors in Eqs. (10)-(11) can be estimated by using the integral identity in Eq. (12), e.g. for R ND (m,m , µ δ , µ δ ) in Eq. (10) follows In order to reduce stochastic noise, we further split up the determinant ratios in Eq. (9) by introducing three intermediate steps with wherem j = ((4 − j)m + jm )/4 and µ δ,j = ((4 − j)µ δ + jµ δ )/4, and using N st = 16 sufficiently suppresses the stochastic fluctuations. Now, observables O(m , µ δ ), like correlations functions, can be evaluated at the new parameter set {m , µ δ } via using the ensemble generated at the old parameter set {m, µ δ }, where · · · is the the ensemble average. Splitting up the reweighting factor R w in several steps enables us to calculate the pion mass for the intermediate three steps, as shown in Fig. 2. After reweighting in steps, we reach the value of m s of the A60s ensemble and a value of M π which agrees with the one extracted by a direct analysis of A60s.
Having demonstrated that reweighting the strange quark mass while maintaining maximal twist works, we can apply it to the cA211.30.32 ensemble. For this ensemble, tuning to maximal twist is done at fixed strange and charm quark mass. Since we do not have a second ensemble at a different strange quark mass tuned to maximal twist to know the target parameters, we now proceed in two steps: we first change the strange quark mass and then re-adjust the bare light mass parameterm to maximal twist. This is done while keeping the charm quark mass constant. The change in the strange quark mass m s → m s + δm s = m s can be corrected by the reweighting factor R 1 = R ND (m s , m s ). Tuning to maximal twist at m s requires the change in the bare mass parameter m →m + δm =m , which can be taken into account via the reweighting factor R 2 = R (m,m )R ND (m,m ). For the cA211.30.32 ensemble, the stochastic fluctuations are suppressed by using N st = 64 for each reweighting factor and performing a single step. With the stochastic evaluation of Eq. (12), the fluctuations are where and the variance is given by The corresponding ratio matrix A is, for the specific example of R 1 , given by A = D ND (m, µ δ )D −1 ND (m , µ δ ). We find that the stochastic fluctuations are small compared to the statistical gauge fluctuations σ 2 ens = var(R)/ R 2 = R 2 / R 2 − 1. For the two different steps, we obtain σ 2 st (R 1 )/σ 2 ens (R 1 ) ∼ 0.02 and σ 2 st (R 2 )/σ 2 ens (R 2 ) ∼ 0.16, similar to what was obtained for clover Wilson fermions, see e.g. [55,56].
The parameters used for the cA211.30.32 ensemble for the different reweighting factors are listed in Table V. We vary the strange quark mass by ∼ 5 MeV. This leads to a significant shift of the PCAC mass (defined in Eq. (6) of Ref. [21]), as shown in Fig. 3. We then readjust by using reweighting in the bare quark mass parameters, following the tuning process outlined in Ref. [21], which establishes a dependence of the PCAC mass on the bare mass parameter given by am PCAC = 1.19(87) + 2.77(202)am .
This leads to δm given by which indeed re-adjusts the PCAC mass to zero, as shown in Fig. 3 and listed in Table V. After reweighting the ensemble cA211.30.32 by varying the strange quark mass by ±5% around its original value and keeping a constant and maximal twist, the resulting values of β 1 and β 2 can be derived as listed in Table II, where correlations between the data points reduce the total statistical error.