Updated evaluation of $\varepsilon_K$ in the Standard Model with lattice QCD inputs

We report a strong tension in $\varepsilon_K$ at the $4\sigma$ level between the experimental value and the theoretical value calculated directly from the standard model using lattice QCD inputs such as $\hat{B}_K$, $|V_{cb}|$, $|V_{us}|$, $\xi_0$, $\xi_2$, $\xi_\text{LD}$, $F_K$, and $m_c$. The standard model with lattice QCD inputs describes only 70% of the experimental value of $\varepsilon_K$, and does not explain its remaining 30%. We also find that this tension disappears when we use the inclusive value of $|V_{cb}|$ (results of the heavy quark expansion based on QCD sum rules) to determine $\varepsilon_K$. This tension is highly correlated with the present discrepancy between the exclusive and inclusive values of $|V_{cb}|$. In order to resolve, in part, the issue with $|V_{cb}|$, it would be highly desirable to have a comprehensive re-analysis over the entire set of experimental data on the $\bar{B} \to D^* \ell \bar{\nu}$ decays using an alternative parametrization of the form factors, such as the BGL parametrization, and a comparison with results of the CLN method.


I. INTRODUCTION
CP violation serves as a natural place to search for new physics [1,2]. The CP violation in the neutral kaon system is, in particular, attractive to us, because the experimental results are already extremely precise [3], and lattice QCD allows us to perform a high precision calculation in kaon physics [4]. In this paper, we focus on the indirect CP violation parameter ε K , which we want to determine using lattice QCD inputs.
Indirect CP violation in the neutral kaon system is parametrized by ε K , ε K ≡ A(K L → ππ(I = 0)) A(K S → ππ(I = 0)) , where K L and K S are the neutral kaon states in nature, and I = 0 represents the isospin of the final two-pion state. In experiment [3], Here, the ε K value represents an ≈ 0.2% impurity of the CP even eigenstate in the K L state, which contains ≈ 99.8% of the CP odd eigenstate. The standard model (SM) describes the CP violation using a single phase in the CKM matrix elements. Hence, if there exists another phase coming from new physics, ε K is a natural place to find it, since ε K is highly sensitive * E-mail: wlee@snu.ac.kr to it. Therefore, it has been one of the top priorities in lattice QCD to calculate ε K to the highest possible precision [4].
In order to evaluate ε K directly from the SM, we need to know 18 input parameters [5]. Out of them, we can, in principle, obtain 7 parameters from lattice QCD:B K , |V cb |, ξ 0 , ξ 2 , ξ LD , |V us |, m c (m c ), and F K . 1 During the last decade, lattice QCD has made such remarkable progress in calculatingB K that its error is only 1.3% at present [6]. At present, the largest error in theoretical calculation of ε K comes from |V cb | [5,[7][8][9].
Here, we would like to report the final results to draw your attention to the key issues. Evaluating ε K directly from the SM with lattice QCD inputs, we find that it has 4.2σ ∼ 3.9σ tension with the experimental result when we use exclusive |V cb |. 2 We also find that this tension disappears with inclusive |V cb |. Hence, it is clear that the key issue is the input value of |V cb |; the 4σ tension in ε K is highly correlated with the 3σ tension between exclusive and inclusive |V cb | [5].
At present, there are two independent methods to determine |V cb |: one is the exclusive method and the other is the inclusive method. In the exclusive method, the experimentalists use the exclusive decaysB → D ( * ) ν to determine |V cb |F(1), and then combine them with lattice QCD results for the form factor F(w) to deter- 1 In this number count, ξ 0 and ξ 2 are redundant. We need to know only ξ 0 , but it is possible to obtain ξ 0 from ξ 2 using ε /ε. For more details, refer to Ref. [5]. 2 Here, the 4.2σ tension is obtained with the estimate of RBC-UKQCD for ξ LD , while the 3.9σ tension is obtained with the BGI estimate. For more details, refer to Section III E.
mine |V cb | [10]. In the final analysis, they also include results for |V cb |/|V ub | obtained by combining the LHCb results for the ratio of the branching fractions between the Λ b → Λ c ν and Λ b → p ν decays with lattice QCD form factors [10]. In the inclusive method, one use the heavy quark expansion (HQE) as the theoretical framework to perform the data analysis onB → X c ν decay processes [10]. The current status of |V cb | is, in units of 1.0 × 10 −3 , exclusive |V cb | = 39.13 ± 0.59 from Ref. [10] (3) inclusive |V cb | = 41.98 ± 0.45 from Ref. [10] (4) where the result in Eq. (4) is obtained in the 1S scheme. The difference between (3) and (4) is 3.8σ. This gap between exclusive and inclusive |V cb | has not been resolved yet. However, a number of interesting ideas have been proposed in order to resolve this issue [11,12]. We review them in Section III C 1 and Appendix A when we discuss |V cb |.
The main goal of this paper is to present the most upto-date results for ε K obtained directly from the SM by using lattice QCD and experimental inputs. In Section II, we review the master formula for ε K and describe each term in detail, including the physical meaning. In Section III, we explain how to obtain the 18 input parameters one by one. In the case of |V cb |, caveats in various methods for the form factor parametrization are addressed in some detail. In Section IV, we present results for ε K obtained using various combinations of input parameters. In Section V, we conclude.

II. REVIEW OF εK
A. Master Formula: εK In the standard model (SM), the direct CP violation parameter ε K in the neutral kaon system can be reexpressed in terms of the well-known SM parameters as follows, This is the master formula, and its derivation is well explained in Ref. [5]. Here, we use the same notation and convention as in Ref. [5].

B. Short Distance Contribution to εK
In the master formula of Eq. (5), the dominant leadingorder effect (≈ +107%) comes from the short distance (SD) contribution proportional toB K . Here, C ε is a dimensionless parameter defined as: Here, X SD represents the short distance effect from the Inami-Lim functions [13]: where λ i = V * is V id is a product of the CKM matrix elements with i = u, c, t, and η ij with i, j = c, t represent the QCD corrections of higher order in α s [14]. There exists a potential issue with poor convergence of perturbation theory for η cc at the charm scale, which is discussed properly in Ref. [5]. Here, S 0 's are Inami-Lim functions [13] defined as where i = c, t, x i = m 2 i /M 2 W , and m i = m i (m i ) is the scale invariant MS quark mass. In X SD of Eq. (7), the S 0 (x t ) term from the top-top contribution in the box diagrams describes about +72.4% of X SD , the S 0 (x c , x t ) term from the top-charm contribution takes over about +45.4% of X SD , and the S 0 (x c ) term from the charmcharm contribution depicts about −17.8% of X SD .
Here, the kaon bag parameterB K is defined aŝ where b(µ) is the renormalization group (RG) running factor to makeB K invariant with respect to the renormalization scale and scheme: Here, details on K + (µ) are given in Ref. [5].

C. Long Distance Contribution to εK
There are two kinds of long distance (LD) contributions on ε K : one is the absorptive LD effect from ξ 0 and the other is the dispersive LD effect from ξ LD . The absorptive LD effects are defined as tan They are related with each other through ε : The overall contribution of the ξ 0 term to ε K is about −7%.
The dispersive LD effect is defined as where if the CPT invariance is well respected. The overall contribution of the ξ LD to ε K is about ±2%.

D. Erratum
There were two pure typos in Ref. [5]. One typo is found in Eq. (50) of Ref. [5]. The correct equations for S 0 (x i ) and S 0 (x i , x j ) are given in Eq. (9) of this paper. The other typo is found in Eq. (62) of Ref. [5]. The correct equation for δm LD is if the CPT invariance is well respected. The −i factor is missing in Ref. [5]. In our actual calculation of ε K , we used the correct equations with no mistake, even through we introduced the above two typos in writing up the paper of Ref. [5].

III. INPUT PARAMETERS
We need to know values of 18 parameters defined in the standard model (SM) in order to evaluate ε K directly from the SM. Out of the 18 parameters, we can obtain, in principle, 7 parameters such asB K , |V cb |, ξ 0 , ξ 2 , ξ LD , |V us |, F K , and m c (m c ) directly from lattice QCD. Here, we describe how to obtain the 18 input parameters from the experiments and from lattice QCD results in detail.

A. Wolfenstein Parameters
The CKMfitter [15] and UTfit [16] collaborations provide the Wolfenstein parameters [17] (λ,ρ,η) determined by the global unitarity triangle (UT) fit. The 2017 results are summarized in Table I. As pointed out in Ref. [5], the Wolfenstein parameters extracted by the global UT fit have unwanted correlation with ε K , because ε K is used as an input to obtain them. Hence, in order to avoid this correlation, we take another set of the Wolfenstein parameters determined from the angle-only-fit (AOF) suggested in Ref. [18]. In the AOF, ε K ,B K , and |V cb | are not used as inputs to determine the UT apex (ρ,η). Then, we determine λ from |V us | which is obtained from the K 2 and K 3 decays using the lattice QCD results. The Wolfenstein parameter A is determined directly from |V cb |, which will be discussed later in Section III C. The Wolfenstein parameters from the AOF are summarized in Table I.

B.BK
In the FLAG review [6], they present lattice QCD results forB K with N f = 2, N f = 2+1, and N f = 2+1+1. Recent calculations ofB K in lattice QCD have been done with N f = 2 + 1 + 1 dynamical quarks [20]. We do not prefer it for two physical reasons: 1. The master formula in Eq. (5) is derived by integrating out the heavy particles including the charm quark. Hence, it is complicated and inconvenient to use the results forB K with dynamical charm quarks [6,20].

2.
A proper and systematic procedure to incorporate the charm quark effect in the calculation of ε K is available in Refs. [21,22]. However, this new method is in the stage of exploratory study and has not reached the stage of precision measurement yet.
Similarly, we prefer using the results forB K with N f = 2 + 1 to those with N f = 2 because they are obtained by quenching the vacuum polarization contributions of the strange quark.
In Table II, we present the FLAG results forB K with N f = 2 + 1. Here, they take a global average over the four data points from BMW 11 [26], Laiho 11 [25], RBC/UKQCD 14 [24], and SWME 15 [23]. The FLAG 17 in the table represents the final result forB K . Here, we use this for our evaluation of ε K .
C. |V cb | In Table III, we summarize updated results for both exclusive |V cb | and inclusive |V cb |. Recently HFLAV reported them in Ref. [10]. The results for exclusive |V cb | depend on the lattice QCD calculations of form factors of Refs. [27][28][29]. Here, when we obtain the in-combined results in Table III (b), we neglect the hidden correlation of the inclusive |V cb | between the kinetic scheme and the 1S scheme, even though there must be some correlation because they share some experimental data with each other. Hence, we prefer using results of the 1S scheme to the results of in-combined here. We use the combined results (ex-combined) for the exclusive |V cb | and the results of the 1S scheme for inclusive |V cb | when we evaluate ε K . a In this analysis, they use the lattice QCD results for the semileptonic form factors in Ref. [27]. b In this analysis, they use the lattice QCD results for the semileptonic form factors in Ref. [28]. c In this analysis, they use the lattice QCD results for the semileptonic form factors in Ref. [29].
In Fig. 1, we present results for |V cb | and |V ub |. The big change is that, as of Lattice 2016, the result for exclusive |V cb | fromB → D ν was about one sigma away from that fromB → D * ν (refer to Ref. [8,9] for more details), but in 2017, they are on top of each other, as shown in Fig. 1. The 2017 results forB → D * ν are not visibly different, but those forB → D ν shift downward by about 1σ. The difference is due to several factors acting in concert: The 2017 results of HFLAV include all results from the B factories, BABAR and BELLE, as well as the older results from CLEO and the LEP experiments ALEPH, OPAL, and DELPHI. Before the results are averaged, they are rescaled by HFLAV to updated values of the inputs, and the averages include the effects of correlations.
The light-blue band represents |V cb | determined from theB → D * ν decay mode. The light-green band represents |V cb | determined from theB → D ν decay mode. The yellow band represents |V ub | determined from theB → π ν decay mode. The magenta band represents |V ub /V cb | determined from the LHCb data of the Λ b → Λc ν and Λ b → p ν decay modes. The orange circle represents the combined results for exclusive |V cb | and |V ub | from the B meson and Λ b decays within 1.0σ. The black cross (×) represents the inclusive |V cb | and |V ub | from the heavy quark expansion. The details are given in Table III.
To obtain the 2017 results for exclusive |V cb | (and |V ub |), HFLAV performed a combined fit to all results for the decaysB → D ( * ) ν, B → π ν, and the ratio of the branching fractions for Λ b → p ν and Λ b → Λ c ν. With lattice QCD results for the form factors, the B → π ν decay yields |V ub |, while the ratio of the branching fractions of the Λ b decays yields |V ub |/|V cb |. Due to the addition of more data to the HFLAV analysis, the results for |V ub |/|V cb | shift downward by about 3 4 σ in 2017, while those for |V ub | shift downward by about 0.1σ. For more details, refer to Ref. [10].

Caveats on CLN and BGL
In order to extract a value of |V cb |F(1) (F(w) is a form factor at a recoil point w) from the experiment of B → D ( * ) ν we need to know the functional form of the form factors as a function of w. There have been two kinds of parametrization methods developed to do this job: one is an HQET-dependent method, and the other is an HQET-independent method. 3 The former is a method of Caprini, Lellouch, and Neubert (CLN) in Ref. [30] and its sibling paper [31], and the latter is a method of Boyd, Grinstein, and Lebed (BGL) in Ref. [32] and its sibling papers [33][34][35]. 4 Both of them have been developed on top of the building blocks designed for K 3 decays in Refs. [37][38][39].
Recently, in Refs. [11,12,40], they claim that the gap between the inclusive |V cb | and exclusive |V cb | might be explained in part by the observation of both groups that CLN consistently underestimates the value of exclusive |V cb | compared with that of BGL. In this claim, they refer to the numbers of HFLAV in Ref. [10] which are obtained using the CLN method. Certainly, this claim is interesting enough to deserve our full and careful investigation on it.
Let us first describe the key points of the claim in Refs. [11,12]. In the CLN parametrization, they introduce the form factor h A1 (w), and the ratios of R 1 (w) and R 2 (w) to describe the form factors forB → D * ν decays. Their definition (Eq. (A4)) and detailed explanation are given in Appendix A. Let us directly address the problematic part in CLN. CLN is constructed based on HQET and its perturbative application to the slope and curvature of h A1 (w), R 1 (w), and R 2 (w). CLN was originally designed to have its error in the level of about 2% precision [30]. At present, the trouble is that the experimental precision goes below the 2% level. The lattice QCD results have precision better than that of the 2% level. The typical size of errors from the slope and curvature in R i (w) obtained using the perturbation theory in HQET is about 10% which has a potential to cause 1 ∼ 2% errors in |V cb |. They ( [11,12]) observed that the CLN method consistently underestimates exclusive |V cb | compared with that of the BGL method which is model-independent by construction. Details on BGL are summarized in Appendix A. To support their claim, they used a preliminary unfolded data of BELLE in Ref. [41]. In their conclusion, they recommended comprehensive reanalysis of old experimental data used in Ref. [10] using the BGL method as well as some suggestions to the lattice QCD community.
In Ref. [42], they incorporate all the O(Λ QCD /m c,b ) and O(α s ) contributions in the HQET framework into their analysis for |V cb | based on the CLN method. They find that their results for |V cb | with improved precision agree with those of HFLAV [10]. In Ref. [43], they use the same kind of CLN method as in Ref. [42] and its variations as well as the BGL method to determine |V cb | and semi-leptonic form factors. In this study they find that the slope of the form factor ratio R 1 (w) at zero recoil obtained using the BGL method has potentially large deviation from heavy quark symmetry, and, in addition, has significant tension with the preliminary lattice QCD results of FNAL/MILC [44] and JLQCD [45,46]. They point out that the tensions between the exclusive and inclusive determinations of |V cb | are far away from being considered resolved at present. In Ref. [40], however, they claim that the conclusions previously reached in Ref. [11] are not changed by taking into account heavy quark symmetry. The extraction of |V cb | using the CLN and BGL parametrizations and preliminary BELLE data has been further investigated in Ref. [47].
The claim in Refs. [11,12] is interesting, but far away from conclusive or decisive in that they used only a preliminary subset of the BELLE data, and the BGL results for the R 1 (w) slope has significant violation of heavy quark symmetry and is disfavored by preliminary lattice QCD results [43]. This issue might well be resolved one way or the other, once the next round of comprehensive reanalysis by HFLAV on the old experimental data used in Ref. [10] becomes available. Lattice QCD calculation of semi-leptonic form factors for theB → D ( * ) ν decays at non-zero recoil will be helpful [45]. Hence, please stay tuned for this coming update.

D. ξ0
The absorptive part of long distance effects in ε K is parametrized into ξ 0 . We can express ε /ε in terms of ξ 0 and ξ 2 as follows, There are two independent methods to determine ξ 0 in lattice QCD: one is the indirect method and the other is the direct method. In the indirect method, we determine ξ 0 using Eq. (22) with lattice QCD input ξ 2 and with experimental results for ε /ε, ε K , and ω. In the direct method, we can determine ξ 0 directly using lattice QCD results for Im A 0 combined with experimental results for Re A 0 .
Recently, RBC-UKQCD reported results for ξ 2 in Ref. [48]. Using the indirect method, we can obtain the result for ξ 0 as in Table IV. Recently, RBC-UKQCD also reported results for Im A 0 in Ref. [49]. Using the experimental value of Re A 0 , we can obtain ξ 0 directly from Im A 0 , which is summarized in Table IV. In Ref. [49] RBC-UKQCD also reported the S-wave π − π scattering phase shift for the I = 0 channel: δ 0 = 23.8 (49)(12). This value is 3.0σ lower than the conventional results for δ 0 in Refs. [50] (KPY-2011) and [51,52] (CGL-2001), and [53]. The values for δ 0 are summarized in Table V. In Fig. 2, we show the experimental results for δ 0 with the fitting results of KPY-2011. They (KPY-2011) used a singly subtracted Roy-like equation to do the interpolation around √ s = m K (the physical kaon mass). Their fitting to the experimental data works well from the threshold to √ s ∼ = 800MeV. In this range they use the singly subtracted Roy-like equation to do the fitting.
In Fig. 3, we show the fitting results of both KPY-2011 and CGL-2001 as well as the results of RBC-UKQCD. There is essentially no difference between KPY-2011 and CGL-2001 in the region near √ s = m K . As one can see in Fig. 3 (a), we observe the 3.0σ tension for δ 0 between RBC-UKQCD and KPY-2011. In the case of δ 2 , there is no difference between RBC-UKQCD and KPY-2011 within statistical uncertainty as one can see in Fig. 3 (b). Taking into account all the aspects, we conclude that the direct calculation of Im A 0 and ξ 0 by RBC-UKQCD in Ref. [49] might have unresolved issues. Indeed, preliminary results presented by RBC-UKQCD in Lattice 2018 suggests that this discrepancy might disappear with improved analysis [54].
Therefore, we prefer the indirect method to the direct method for the following two reasons. The first reason is that the lattice QCD calculation of Im A 0 is much noisier than that of Im A 2 thanks to many disconnected diagrams. The second reason is that the S-wave phase shift δ 0 of the π −π scattering in Ref. [49] is lower by 3.0σ than the conventional determination of δ 0 in Ref. [50,51,53], which indicates that there might be some issues unresolved at present. In Table IV, we present results of ξ 0 determined using both indirect and direct methods.
Here, we use the value of ξ 0 obtained using the indirect method.
One remaining caveat is that the ξ 0 and ξ 2 in Ref. [48] is calculated in the isospin symmetric limit. The isospin breaking effects on ε /ε are studied in Refs. [55,56]. These studies conclude that the isospin violation correction in the CP violation correction for ε is below 15% within the uncertainties of large N c estimates for the low energy constants. Since ξ 0 has an effect of about −7% on ε K , the isospin violation effect maximum, 15% of ξ 0 amounts to ±1% correction for ε K . Here, we neglect this effect completely without loss of generality in our conclusion.

E. ξ LD
The long distance (LD) effects on ε K are explained in Sec. II C. Hence, here we would like to summarize the recent progress in calculating the LD effects in lattice QCD.
Lattice QCD tools to calculate the dispersive LD effect, ξ LD are well established in Ref. [21,22,57]. In addition, recently, there have been a number of attempts to calculate ξ LD on the lattice [58,59]. In these attempts, RBC-UKQCD used pion mass of 329MeV and kaon mass of 591MeV. Hence, the energy of the two pion state and three pion states are heavier than the kaon mass. Therefore, the sign of the denominator in Eq. (18) is opposite to that of the physical contribution in which the two and three pion state energy is lighter than the kaon mass. Therefore, this attempt in Refs. [58,59] belongs to the category of exploratory study rather than to that of precision measurement.
The net contribution of ξ LD to ε K in Eqs. (17) and (18) turns out to be of the same order of magnitude as ξ 0 using chiral perturbation theory [60]. They claim that where we use the indirect results of ξ 0 given in Table IV, including its error. Here, we call this method the BGI estimate for ξ LD . This also indicates that ξ LD is at most a 4% correction to ε K . This claim is highly consistent with the estimate of about 2% in Ref. [21,61]: Here, we call this method the RBC-UKQCD estimate for ξ LD .
In this paper, we use both of the above estimates of ξ LD with the BGI and RBC-UKQCD methods to determine ε K .

F. Top quark mass
The pole mass of top quarks coming from Ref. [3] is The pole and MS masses are related as follows, where m t (µ) is the MS mass renormalized at scale µ.
Here, Z OS is the renormalization factor in the on-shell scheme, and Z MS is the renormalization factor in the MS scheme. The top scale-invariant quark mass µ t is the MS mass m t (µ) with the scale µ set equal to the scaleinvariant mass, where we use the four-loop perturbation formula for z(µ t ). Details on the four-loop conversion formula are described in Appendix B. In Eq. (27), the first error comes from the error of the top pole mass, and the second error represents the uncertainty due to truncation of higher loops in the conversion formula which is estimated as the difference in µ t between the 3-loop and 4-loop formulas.
We have neglected the renormalon ambiguity and corrections due to the three-loop fermion mass such as m b (bottom quark mass) and m c (charm quark mass).

G. Other Input Parameters
For the higher order QCD corrections η cc , η ct , and η tt , we use the same values as in Ref. [5]. They are summarized in Table VI. Other input parameters are summarized in Table VII. They are the same as Ref. [5] except for charm quark mass m c (m c ), the kaon mass m K 0 , and the kaon decay constant F K . For the charm quark mass, we use the HPQCD results of m c (m c ) reported in Ref. [64]. For the kaon mass, we use the updated results of Particle Data Group (PDG) in Ref. [3]. For the kaon decay constant, we use the updated results of PDG reported in Ref. [3], which are obtained from the FLAG data [6].

IV. RESULTS
A. RBC-UKQCD estimate for ξ LD In Fig. 4, we present results for ε K calculated directly from the standard model with the lattice QCD inputs described in Section III. In Fig. 4 (a), the blue curve which encircles the histogram represents the theoretical evaluation of ε K using the FLAG-2017B K , AOF for Wolfenstein parameters, and exclusive |V cb | which corresponds to ex-combined in Table III (a), and the RBC-UKQCD estimate for ξ LD . The red curve in Fig. 4 represents the experimental result for ε K . In Fig. 4 (b), the blue curve represents the same as in Fig. 4 (a) except for using the inclusive |V cb | which corresponds to 1S scheme in Table  III The updated results for |ε K | are, in units of 1.0 × 10 −3 , presented in Table VIII. From Table VIII, we observe that the theoretical evaluation of |ε K | with lattice QCD inputs (with exclusive |V cb |), which corresponds to |ε K | SM excl , has 4.2σ tension with the experimental result |ε K | Exp , while there is no tension in the inclusive |V cb | channel (heavy quark expansion based on the OPE and QCD sum rules). In Fig. 5 (a), we plot the ∆ε K ≡ |ε K | Exp − |ε K | SM excl in units of σ (which is the total error of ∆ε K ) as the time evolves starting from 2012. We began to monitor ∆ε K in 2012 when several lattice QCD results forB K obtained using different discretization methods for the light and strange quarks became consistent with one another within one sigma. In 2012, ∆ε K was 2.5σ, but now it is

4.2σ.
To understand the change of ∆ε K /σ with respect to time, we have performed an additional analysis on the average and error.
In Fig. 5 (b), we plot the time evolution of the average ∆ε K and the error σ ∆ε K . Here, we find that the average of ∆ε K has increased with some fluctuations by 27% during the period of 2012-2018, and its error σ ∆ε K has decreased monotonically by 25% in the same period. These two effects interfere constructively to produce the 4.2σ tension in ∆ε K in 2018. We can understand the monotonic decrease in σ ∆ε K in the following way. As time goes on, the lattice QCD calculations are becoming more precise and the experimental results also are becoming more accurate, which constructively leads to the monotonic decrease in σ ∆ε K .
In Table IX, we present the error budget for |ε K | SM excl . Here, we find that the largest error in |ε K | SM excl comes from |V cb |, while the errors coming fromη and η ct are subdominant. Hence, if we are to see a gap ∆ε K greater than 5.0σ, it is essential to reduce the error in |V cb | significantly.
In exclusive |V cb |, there are two major error sources: one is experimental and the other is theoretical. The experimental error is discussed in Section III C, and the resolution is beyond the scope of this paper. The largest error in the theoretical part of |V cb | comes from the heavy quark discretization error (HQDE) for the charm quark in lattice QCD. If one use the Fermilab action, the HQDE is about 1.0%, which is significantly larger than any other error in the theoretical side. In order to reduce the HQDE by a factor of ∼ 1/5, there are on-going efforts to use the OK action to calculate the B → D ( * ) semileptonic form factors in Refs. [65][66][67].  Here, σ = σ∆ε K represents the error of ∆εK . The ∆εK is obtained using the same input parameters as in Fig. 4 (a) for the exclusive |V cb | channel.

B. BGI estimate for ξ LD
Here, we present the results obtained using the BGI estimate in Eq. (23) for ξ LD .
In Fig. 6 (a), the blue curve represents the theoretical evaluation of |ε K | directly from the standard model (SM) using the same input parameters as in Fig. 4 (a) except for the BGI estimate in Eq. (23) for ξ LD . The red curve in Fig. 6 represents the experimental result for |ε K |. In Fig. 6 (b), the blue curve represents the same as in Fig. 6 (a) except for using inclusive |V cb | (1S scheme in Table III (b)).
Results for |ε K | in Fig. 6 are summarized in Table X. From Table X, we find that the value for |ε K | SM excl (the theoretical evaluation of |ε K | with lattice QCD inputs such as exclusive |V cb |) has 3.9σ tension with the experimental result |ε K | Exp , whereas there is no tension in the inclusive |V cb | channel (with heavy quark expansion and QCD sum rules).
In Fig. 7 (a), we plot ∆ε K ≡ |ε K | Exp − |ε K | SM excl in units of σ (the total error of ∆ε K ) as a function of time starting from 2012. In 2012, ∆ε K was 2.3σ, but now it is 3.9σ. To understand this transition, we have done an additional  analysis on the average and error. In Fig. 7 (b), we plot the time evolution of the average and error for ∆ε K . Here, we find that the average of ∆ε K has increased by 29% with some fluctuations during the period of 2012-2018, and its error has decreased by 25% monotonically in the same period. These two effect has produced, constructively, the 3.9σ tension in ∆ε K .
In Table XI, we present the error budget for |ε K | SM excl . Here, we find that the largest error in |ε K | SM excl still comes from |V cb |. Here, note that the error from ξ LD (the BGI estimate) is larger than that fromρ, which is different from Table IX. In summary, if we are to observe the gap ∆ε K greater than 5.0σ, it is essential to reduce the error in |V cb | significantly.

V. CONCLUSION
In this paper, we find that there exists a remarkable gap of 4.2σ ∼ 3.9σ in ε K between experiment and the SM theory with lattice QCD inputs. The upper bound of 4.2σ tension is obtained with the RBC-UKQCD estimate for ξ LD . The lower bound of the 3.9σ tension is obtained when we use the BGI estimate for ξ LD . In the BGI estimate [60], they added 50% more error to be on the safe side and more conservative. Even if we remove this 50% bubble in the error of the BGI estimate, we end up with the same tension of 3.9σ. To obtain this result, we choose the angle-only-fit (AOF), exclusive |V cb | from lattice QCD, and FLAGB K (N f = 2 + 1) from lattice QCD, to determine the theoretical value for ε K directly from the SM. In 2015, we reported a 3.4σ tension between |ε K | SM excl and |ε K | Exp [5], and the tension is 4.2σ  Here, σ = σ∆ε K represents the error of ∆εK . The ∆εK is obtained using the same input parameters as in Fig. 6 (a) for the exclusive |V cb | channel.
at present. 5 We find that the tension between |ε K | SM excl and |ε K | Exp continues to increase during the period of 2012-2018. Part of the reason is that the uncertainties of results for the SM input parameters continues to decrease monotonically. In Table XII, we present how the values of ∆ε K have changed from 2015 to 2018. Here, we find that the positive shift of ∆ε K is about the same for the inclusive and exclusive values of |V cb |. This reflects the changes of other input parameters since 2015. We also note that there is no significant tension observed yet for inclusive |V cb |, which is obtained using the heavy quark expansion based on the QCD sum rules.
There has been an interesting claim [11,12] which has potential to resolve the issue of the inconsistency between the exclusive and inclusive |V cb |. However, this claim is far away from conclusive yet since it is based on an analysis over a preliminary and specific subset of experimental data. We find that it would be highly desirable if an experimental group were to perform a comprehensive reanalysis over the entire set of experimental data for thē B → D * ν decays using an alternative parametrization method for the form factors, and compare results with those of CLN. 5 In 2015, we used the RBC-UKQCD estimate for ξ LD .

ACKNOWLEDGMENTS
We would like to express our sincere gratitude to Carleton Detar, Aida El-Khadra, and Andreas Kronfeld for helpful discussion. We also would like to express sin- Appendix A: Brief summary on CLN and BGL Let us considerB → D * ν decays. The recoil variable w is defined as w ≡ v B ·v D * , where v B is the four velocity of the mother particle (B meson) and v D * is that of the daughter particle (D * meson). The differential decay rate [68] is given by where G F is Fermi's constant, η EW is a small electroweak correction, and F(w) is the form factor. The kinematic factor χ(w) is where r ≡ m D * /m B . So far the formalism is quite general. We may express the form factor as follows , without loss of generality, In the CLN method [30], the form factor functions are parametrized as follows, where z is a typical conformal mapping variable defined as The basic idea of CLN is a zero-recoil expansion around w = 1. The slope and curvature of R 1 (w) and R 2 (w) are determined by perturbation theory and O(1/M ) corrections at the leading order using Heavy Quark Effective Theory (HQET) [68]. The original claim of CLN [30] is that the accuracy of h A1 (w) is better than 2%, which makes us become apprehensive as we get |V cb | in the precision level below 2%. In addition, the slope and curvature of R 1 (w) and R 2 (w) contain truncation errors coming from O(Λ 2 /m 2 c ) and O(α s Λ/m c ) corrections as well as those uncertainties due to the QCD sum rules on which it is based [11]. Typically, Λ/m c ≈ 1/3 and α s ≈ 0.3, which implies that the accuracy of the slope and curvature in the ratio R i (w) is only in the 10% level.
Using the CLN method, experimentalists perform a four parameter fit of η EW F(1)|V cb |, ρ 2 , R 1 (1) and R 2 (1) to some unfolded data in experiment [10]. Since η EW is very well known and lattice QCD can determine F(1) very precisely, we can determine |V cb | from the experimental fits.
Let us switch the gear to BGL. In the case of CLN, it is built on the basis of HQET and its perturbative expansion. Unlike CLN, BGL is an HQET-independent approach to the form factor parametrization. The basic idea of BGL is composed of three building blocks: dispersion relationship, analytic continuation, and crossing symmetry.
Let us begin with the first building block: dispersion relation. In QCD, consider the two point function of flavor changing current In general, Π T,L J (q 2 ) is not finite. Hence, in order to obtain finite dispersion relations, we need to make one or two subtractions as follows, Let us introduce the Källen-Lehmann spectral decomposition by inserting a complete set of states X into the two point function.
where the sum includes an integral over the phase space allowed to each state X which has the same quantum number as the current J. The positivity of Im Π T J (q 2 ) and Im Π L J (q 2 ) follows from Eq. (A12) [39]. In other words, for any complex 4-vector ξ µ . This implies that where the sum is over polarizations of H b and H c states, and the ellipsis denotes strictly positive contributions from the higher resonances and multi-particle states (three-body or higher multi-body states). Here, we assume that H b = B, B * meson states, and H c = D, D * meson states. Since the right-hand side (RHS) of Eq. (A16) is a sum of positive contributions, we can obtain the following simple inequality.
where t = q 2 , k(t) is a calculable kinematic function arising from two-body phase space, and F (t) is the form factor associated with a specific decay of our interest. For example, in the case ofB → D * ν decay channel, At this stage, we need to use the second building block: crossing symmetry [33]. Let us define t ± ≡ (M H b ± M Hc ) 2 . The crossing symmetry insures that the 0|J i |H b H c amplitude which shows up in pair production of the H b and H c mesons from a virtual W boson shares the same form factor F (t) as the H c |J i |H b amplitude, while we can connect the pair production region t + ≤ t < ∞ with the semi-leptonic region m 2 ≤ t ≤ t − through the analytic continuation (the third building block).
Let us define the hadronic moments χ (n) J as in Ref. [33,34], where χ (0) J = χ J . Hence, from the inequality in Eq. (A17), we can obtain the following inequality: where t + represents the pair production threshold. At this stage, we need to introduce a key idea of quarkhadron duality in QCD sum rules which claims that the hadronic moments χ (n) J can be calculated in perturbative QCD at q 2 = 0 [34]. Then, we can rewrite the inequality in Eq. (A21) as follows, The inequality in Eq. (A22) imposes an upper bound on the form factor F (t) in the pair production region (t + ≤ t < ∞).
In order to turn Eq. (A22) into a constraint in the semileptonic region (m 2 ≤ t ≤ t − ), we need to use the third building block: analyticity which allows us to extend the analytic region of the integrand to the region below the pair-production threshold (t < t + ). To do this, it is convenient to introduce a conformal mapping function z is real for t < t + , z = −1 at t = t + , zero at t = t s , and a U(1) phase z = e iθ for t > t + . t s ≥ t > −∞ maps into 0 ≤ z < 1 along the real axis. The upper contour of t + ≤ t < ∞ along the real axis maps into the upper half of a unit circle (π ≥ θ > 0 for z = e iθ ). Similarly, the lower contour of t + ≤ t < ∞ along the real axis maps into the lower half of a unit circle (π ≤ θ < 2π for z = e iθ ). All the poles in the integrand of Eq. (A22) can be removed by multiplying by various powers of z(t, t s ), if we know the positions t s of the sub-threshold poles in F (t) and h (n) (t). Each pole has a distinct value of t s , and the product z(t, t s1 )z(t, t s2 )z(t, t s3 ) · · · can remove all of them. For example, such poles include the contribution of B c resonances to the form factor F (t) as well as singularities in the kinematic function h (n) (t).
Once we determine the pole positions phenomenologically, we can rewrite the inequality in Eq. (A22) as follows, where the outer function φ is The factorP (t) is a product of z(t, t s )'s and z(t, t s )'s such that t s is chosen to remove the sub-threshold poles and branch cuts in the kinematic function h (n) (t). The Blaschke factor P (t) is where t Pi = M 2 Pi represents the pole positions of F (t) below the threshold (t Pi < t + ), and N is the number of the sub-threshold poles in F (t). Here, note that z Pi is real (z * Pi = z Pi ) for the sub-threshold poles. In addition, note that P (t) is unimodular (|P (t)| = 1), if z is unimodular (|z| = 1). We have a full freedom to choose t 0 . Here, we set t 0 = t − for convenience and simplicity, and without loss of generality. 6 Since φ(t, t 0 )P (t)F (t) is analytic even in the subthreshold region, it is possible to expand this in powers of z(t, t 0 ). Hence, This is called the BGL method of form factor parametrization [32].
In addition, the integral in Eq. (A25) can be rewritten as follows, 1 π where z = e iθ above the threshold (t + ≤ t < ∞). Hence, the final version of the inequality after the Fourier analysis is ∞ n=0 a 2 n ≤ 1 .
This is called the unitarity condition (the weak version). 7 For theB → D * ν decay process that are the main subject of this paper, z(t, t 0 ) is in the physical region of 0 ≤ z ≤ 0.056 for any physical momentum transfer m 2 ≤ t ≤ t − . Hence, in practice, it is possible to truncate the expansion after the first two or three terms. Since the BGL method does not use any model to constrain a n , it is model-independent by construction. We follow the terminology in the literature. The pole and MS masses are related by the ratio of (mass) renormalization factors z, 6 For other choices of t 0 , refer to Ref. [34]. 7 A stronger version can be obtained simply by adding more decay channels in the right-hand side of the inequality in Eq. (A17) [11].
where m(µ) is the MS mass renormalized at scale µ, M is the pole mass, Z OS is the renormalization factor in the on-shell scheme, and Z MS is the renormalization factor in the MS scheme. The ratio z depends on the scale µ via the strong coupling α s (µ) and ln(µ/M ). The top scaleinvariant (SI) mass µ t is the MS mass m t (µ) with the scale µ set equal to the scale-invariant mass, Given the perturbative expansion of z(µ) in powers of α s (µ), Eq. (B2) gives the SI mass in terms of α s (µ t ) and ln(µ t /M t ), i.e., in terms of the SI mass and the pole mass.
To obtain a formula for the SI mass in terms of the pole mass alone, one can iterate the perturbative expansion. The result is an expansion of µ t in powers of α s (M t ).

Three-loop result
For the three-loop conversion, we use Eq. (16) of Ref. [69]. For the top mass conversion, this equation is an expansion of the ratio z(µ t ) in powers of the coupling α (6) s (M t ), with six active quark flavors. We obtain the coupling α  [69] and using Eq. (15). We have verified that doing so yields results consistent with Refs. [69,70] and the RunDec3 code [71]. Our result for the top SI mass is µ t = 163.82 ± 1.05 GeV , where the uncertainty is propagated from the uncertainty in the top pole mass, and all other uncertainties are neglected. Below we detail the inputs and the steps of the calculation.