Neutral B-meson mixing from full lattice QCD at the physical point

TL;DR

This paper presents the most precise lattice QCD calculations of bag parameters for neutral B-meson mixing, enabling improved determinations of CKM matrix elements and rare decay branching fractions within and beyond the Standard Model.

Contribution

First full four-flavour lattice QCD calculation of bag parameters for neutral B-meson mixing at the physical point, including beyond Standard Model operators.

Findings

Most accurate bag parameters to date for all five operators.

Determined CKM elements V_ts and V_td with high precision.

Predicted branching fractions for B_s and B_d to mu+ mu- decays.

Abstract

We calculate the bag parameters for neutral $B$-meson mixing in and beyond the Standard Model, in full four-flavour lattice QCD for the first time. We work on gluon field configurations that include the effect of $u$, $d$, $s$ and $c$ sea quarks with the Highly Improved Staggered Quark (HISQ) action at three values of the lattice spacing and with three $u/d$ quark masses going down to the physical value. The valence $b$ quarks use the improved NRQCD action and the valence light quarks, the HISQ action. Our analysis was blinded. Our results for the bag parameters for all five operators are the most accurate to date. For the Standard Model operator between $B_s$ and $B_d$ mesons we find: $\hat{B}_{B_s}=1.232(53)$, $\hat{B}_{B_d}=1.222(61)$. Combining our results with lattice QCD calculations of the decay constants using HISQ quarks from the Fermilab/MILC collaboration and with experimental values for $B_s$ and $B_d$ oscillation frequencies allows determination of the CKM elements $V_{ts}$ and $V_{td}$. We find $V_{ts} = 0.04189(93)$, $V_{td} = 0.00867(23)$ and $V_{ts}/V_{td} = 0.2071(27)$. Our results agree well (within $2\sigma$) with values determined from CKM unitarity constraints based on tree-level processes (only). Using a ratio to $\Delta M$ in which CKM elements cancel in the Standard Model, we determine the branching fractions ${\text{Br}}(B_s\rightarrow \mu^+\mu^-) = 3.81(18) \times 10^{-9}$ and ${\text{Br}}(B_d\rightarrow \mu^+\mu^-) = 1.031(54) \times 10^{-10}$. We also give results for matrix elements of the operators $R_0$, $R_1$ and $\tilde{R}_1$ that contribute to neutral $B$-meson width differences.

Tables

Table 1. Table 1: Normalizations η i q ​ ( m b ) superscript subscript 𝜂 𝑖 𝑞 subscript 𝑚 𝑏 \eta_{i}^{q}(m_{b}) for bag parameters (Eq. ( 4 )). These are calculated using M B s subscript 𝑀 subscript 𝐵 𝑠 M_{B_{s}} =5.3669(2) GeV and M B d subscript 𝑀 subscript 𝐵 𝑑 M_{B_{d}} =5.2796(2) GeV  Tanabashi  et al. ( 2018 ) , m ¯ b ​ ( m ¯ b ) subscript ¯ 𝑚 𝑏 subscript ¯ 𝑚 𝑏 \overline{m}_{b}(\overline{m}_{b}) = 4.162(48) GeV and m b / m s subscript 𝑚 𝑏 subscript 𝑚 𝑠 m_{b}/m_{s} =52.55(55)  Chakraborty  et al. ( 2015 ) , and m s / m l subscript 𝑚 𝑠 subscript 𝑚 𝑙 m_{s}/m_{l} =27.18(10)  Bazavov  et al. ( 2018 ) .

$B_q$	${\eta }_1^q$	${\eta }_2^q$	${\eta }_3^q$	${\eta }_4^q$	${\eta }_5^q$
$B_s$	2.667	$-2.669\left(62\right)$	0.534 (12)	3.536 (74)	2.068 (25)
$B_d$	2.667	$-2.678\left(62\right)$	0.536 (12)	3.547 (74)	2.071 (25)

Table 2. Table 2: Perturbative coefficients used in Eq. ( II.2 ) to convert matrix elements of lattice NRQCD-HISQ 4-quark operators into MS ¯ ¯ MS {\overline{\mathrm{MS}}}  matrix elements. Results are given for the NRQCD valence b 𝑏 b -quark masses (in lattice units) used with our different ensembles. The continuum scheme used is the MS ¯ NDR subscript ¯ MS NDR {\overline{\mathrm{MS}}}_{\mathrm{NDR}} scheme of  Beneke  et al. ( 1999 ) (BBGLN) with μ = m b 𝜇 subscript 𝑚 𝑏 \mu=m_{b} . The coefficients come from  Monahan  et al. ( 2014 ) , with z i ​ j ≡ ρ i ​ j − ζ i ​ j subscript 𝑧 𝑖 𝑗 subscript 𝜌 𝑖 𝑗 subscript 𝜁 𝑖 𝑗 z_{ij}\equiv\rho_{ij}-\zeta_{ij} where ρ i ​ j subscript 𝜌 𝑖 𝑗 \rho_{ij} and ζ i ​ j subscript 𝜁 𝑖 𝑗 \zeta_{ij} are listed in Tables III and IV of that paper 1 1 1 Note that we have corrected two typographical errors, for ρ 21 subscript 𝜌 21 \rho_{21} for a ​ m b = 2.66 𝑎 subscript 𝑚 𝑏 2.66 am_{b}=2.66 and ζ 22 subscript 𝜁 22 \zeta_{22} for a ​ m b = 2.62 𝑎 subscript 𝑚 𝑏 2.62 am_{b}=2.62 . . The perturbative coefficients  z A 0 subscript 𝑧 subscript 𝐴 0 z_{A_{0}} for the temporal axial current (Eq. ( 34 )) are also listed; these are from  Dowdall  et al. ( 2013a ) , which used results from  Monahan  et al. ( 2013 ) .

$a{m}_b$	$z_{11}$	$z_{12}$	$z_{22}$	$z_{21}$	$z_{33}$	$z_{31}$	$z_{44}$	$z_{45}$	$z_{55}$	$z_{54}$	$z_{A_0}$
3.297	$-0.472\left(2\right)$	$-0.299\left(2\right)$	0.440 (2)	0.041 (2)	0.036 (2)	0.092 (2)	0.646 (2)	$-0.252\left(2\right)$	$-0.141\left(2\right)$	$0.111\left(2\right)$	$0.024\left(2\right)$
3.263	$-0.469\left(2\right)$	$-0.296\left(2\right)$	0.438 (2)	0.041 (2)	0.038 (2)	0.091 (2)	0.640 (2)	$-0.251\left(2\right)$	$-0.140\left(2\right)$	$0.108\left(2\right)$	$0.022\left(2\right)$
3.25	$-0.469\left(2\right)$	$-0.294\left(2\right)$	0.438 (2)	0.041 (2)	0.040 (2)	0.091 (2)	0.639 (2)	$-0.252\left(2\right)$	$-0.139\left(2\right)$	$0.106\left(2\right)$	$0.022\left(2\right)$
2.66	$-0.429\left(2\right)$	$-0.235\left(2\right)$	0.394 (2)	0.044 (2)	0.101 (2)	0.080 (2)	0.514 (2)	$-0.254\left(2\right)$	$-0.127\left(2\right)$	$0.037\left(2\right)$	$0.006\left(2\right)$
2.62	$-0.427\left(2\right)$	$-0.229\left(2\right)$	0.388 (2)	0.044 (2)	0.105 (2)	0.080 (2)	0.501 (2)	$-0.254\left(2\right)$	$-0.128\left(2\right)$	$0.032\left(2\right)$	$0.001\left(2\right)$
1.91	$-0.296\left(2\right)$	$-0.108\left(2\right)$	0.340 (2)	0.045 (2)	0.259 (2)	0.053 (2)	0.299 (2)	$-0.243\left(2\right)$	$-0.063\left(2\right)$	$-0.084\left(2\right)$	$-0.007\left(2\right)$

Table 3. Table 3: Parameters of the gauge ensembles used in this calculation. β 𝛽 \beta is the gauge coupling, a Υ subscript 𝑎 Υ a_{\Upsilon} is the lattice spacing as determined by the Υ ​ ( 2 ​ S − 1 ​ S ) Υ 2 𝑆 1 𝑆 \Upsilon(2S-1S) splitting in Dowdall  et al. ( 2012a ) , where the three errors are statistics, NRQCD systematics and experiment. a ​ m l , a ​ m s 𝑎 subscript 𝑚 𝑙 𝑎 subscript 𝑚 𝑠 am_{l},am_{s} and a ​ m c 𝑎 subscript 𝑚 𝑐 am_{c} are the sea quark masses, L × T 𝐿 𝑇 L\times T gives the spatial and temporal extent of the lattices and n cfg subscript 𝑛 cfg n_{{\rm cfg}} is the number of configurations in each ensemble. We use 16 time sources on each configuration to improve statistics.

Set	$\beta$	$a_{{\rm Y}}$ (fm)	$a{m}_l$	$a{m}_s$	$a{m}_c$	$L\times T$	$n_{\mathrm{cfg}}$
1	5.8	0.1474(5)(14)(2)	0.013	0.065	0.838	16$\times$48	1000
2	5.8	0.1463(3)(14)(2)	0.0064	0.064	0.828	24$\times$48	1000
3	5.8	0.1450(3)(14)(2)	0.00235	0.0647	0.831	32$\times$48	1000
4	6.0	0.1219(2)(9)(2)	0.0102	0.0509	0.635	24$\times$64	1000
5	6.0	0.1195(3)(9)(2)	0.00507	0.0507	0.628	32$\times$64	1000
6	6.0	0.1189(2)(9)(2)	0.00184	0.0507	0.628	48$\times$64	1000
7	6.3	0.0884(3)(5)(1)	0.0074	0.037	0.440	32$\times$96	1007

Table 4. Table 4: Parameters used for the valence quarks. a ​ m b 𝑎 subscript 𝑚 𝑏 am_{b} is the bare b 𝑏 b quark mass in lattice units, u 0 ​ L subscript 𝑢 0 𝐿 u_{0L} is the Landau link value used for tadpole-improvement, and a ​ m l val 𝑎 superscript subscript 𝑚 𝑙 val am_{l}^{\rm val} , a ​ m s val 𝑎 superscript subscript 𝑚 𝑠 val am_{s}^{\rm val} are the HISQ light and strange quark valence masses. We also tabulate the values from  Colquhoun  et al. ( 2015 ) for the running coupling constant  α s subscript 𝛼 𝑠 \alpha_{s} to be used in Eq. ( II.2 ) for matching lattice 4-quark operators to the continuum. This is in the V-scheme at scale ( 2 / a ) 2 𝑎 (2/a) .

Set	$a{m}_b$	$u_{0L}$	$a{m}_l^{\mathrm{val}}$	$a{m}_s^{\mathrm{val}}$	${\alpha }_s$
1	3.297	0.8195	0.013	0.0641	0.346
2	3.263	0.82015	0.0064	0.0636	0.345
3	3.25	0.819467	0.00235	0.0628	0.343
4	2.66	0.834	0.01044	0.0522	0.311
5	2.62	0.8349	0.00507	0.0505	0.308
6	2.62	0.834083	0.00184	0.0507	0.307
7	1.91	0.8525	0.0074	0.0364	0.267

Table 5. Table 5: MS ¯ ¯ MS {\overline{\mathrm{MS}}} matrix elements obtained from Eq. ( II.2 ) together with simulation data from Table  8 . Results are given for each configuration set and both B s subscript 𝐵 𝑠 B_{s} (top) and B d subscript 𝐵 𝑑 B_{d}  (bottom) mesons. Values are also given (“phys.”) for our final results at physical quark masses.

	${⟨B_s\left|O_n\right|{\overline{B}}_s⟩}_{\overline{\mathrm{MS}}}^{\left(m_b\right)}/{\left(f_{B_s}{M}_{B_s}\right)}^2$
set	$O_1$	$O_2$	$O_3$	$O_4$	$O_5$
1	2.10 (21)	$-2.14\left(22\right)$	0.442 (59)	3.56 (41)	1.90 (14)
2	2.16 (21)	$-2.20\left(22\right)$	0.441 (59)	3.73 (43)	2.00 (14)
3	2.14 (21)	$-2.18\left(22\right)$	0.441 (58)	3.69 (41)	1.97 (14)
4	2.20 (15)	$-2.19\left(16\right)$	0.443 (48)	3.76 (28)	2.017 (97)
5	2.15 (14)	$-2.16\left(16\right)$	0.432 (47)	3.64 (27)	1.948 (91)
6	2.20 (14)	$-2.20\left(16\right)$	0.445 (48)	3.73 (27)	1.990 (93)
7	2.19 (10)	$-2.19\left(12\right)$	0.443 (37)	3.70 (16)	1.976 (93)
phys.	2.168 (93)	$-2.18\left(10\right)$	0.436 (29)	3.65 (15)	1.945 (76)
	${⟨B_d\left|O_n\right|{\overline{B}}_d⟩}_{\overline{\mathrm{MS}}}^{\left(m_b\right)}/{\left(f_{B_d}{M}_{B_d}\right)}^2$
set	$O_1$	$O_2$	$O_3$	$O_4$	$O_5$
1	2.07 (21)	$-2.12\left(22\right)$	0.438 (59)	3.57 (42)	1.89 (14)
2	2.10 (22)	$-2.13\left(22\right)$	0.421 (60)	3.77 (44)	2.04 (15)
3	2.06 (21)	$-2.10\left(21\right)$	0.398 (58)	3.77 (43)	1.98 (14)
4	2.20 (16)	$-2.15\left(16\right)$	0.403 (48)	3.93 (30)	2.11 (11)
5	2.16 (15)	$-2.06\left(15\right)$	0.396 (48)	3.70 (28)	1.965 (98)
6	2.11 (16)	$-2.17\left(17\right)$	0.447 (52)	3.87 (30)	2.04 (11)
7	2.20 (12)	$-2.14\left(13\right)$	0.413 (39)	3.83 (19)	2.03 (11)
phys.	2.15 (11)	$-2.06\left(11\right)$	0.400 (30)	3.82 (18)	2.015 (92)

Table 6. Table 6: MS ¯ ¯ MS {\overline{\mathrm{MS}}} bag parameters (Eq. ( 2 ) with μ = m b 𝜇 subscript 𝑚 𝑏 \mu=m_{b} ) for the five 4-quark operators. Results are given for both B s subscript 𝐵 𝑠 B_{s}  and B d subscript 𝐵 𝑑 B_{d}  mesons, and for the ratios of bag parameters.

	$B_{B_q}^{\left(1\right)}\left(m_b\right)$	$B_{B_q}^{\left(2\right)}\left(m_b\right)$	$B_{B_q}^{\left(3\right)}\left(m_b\right)$	$B_{B_q}^{\left(4\right)}\left(m_b\right)$	$B_{B_q}^{\left(5\right)}\left(m_b\right)$
$B_s$	0.813 (35)	0.817 (43)	0.816 (57)	1.033 (47)	0.941 (38)
$B_d$	0.806 (40)	0.769 (44)	0.747 (59)	1.077 (55)	0.973 (46)
$B_s/B_d$	1.008 (25)	1.063 (24)	1.092 (34)	0.959 (21)	0.967 (23)

Table 7. Table 7: Percent errors coming from different sources for the B s subscript 𝐵 𝑠 B_{s}  meson’s bag parameters  B B s ( n ) superscript subscript 𝐵 subscript 𝐵 𝑠 𝑛 B_{B_{s}}^{(n)} and B B s ( 1 ) / B B d ( 1 ) superscript subscript 𝐵 subscript 𝐵 𝑠 1 superscript subscript 𝐵 subscript 𝐵 𝑑 1 B_{B_{s}}^{(1)}/B_{B_{d}}^{(1)} (Table  6 ). The total error for each quantity is also shown. The error budgets for the B d subscript 𝐵 𝑑 B_{d}  meson’s bag parameters are very similar. Systematic errors from finite-volume, QED and strong-isospin breaking effects are estimated to be below 0.1% and hence negligible in Appendix  B.5 .

	$B_{B_s}^{\left(1\right)}$	$B_{B_s}^{\left(2\right)}$	$B_{B_s}^{\left(3\right)}$	$B_{B_s}^{\left(4\right)}$	$B_{B_s}^{\left(5\right)}$	$B_{B_s}^{\left(1\right)}/B_{B_d}^{\left(1\right)}$
lattice data	1.4	1.4	1.5	1.6	1.5	1.5
${\eta }_i^q$	0.0	2.3	2.3	2.1	1.2	0.0
${\alpha }_s^2$ terms	2.1	2.9	5.2	1.9	1.5	0.1
${\alpha }_s{\Lambda }_{\mathrm{QCD}}/m_b$ terms	2.9	2.8	2.9	2.8	2.7	0.0
${\left(a{\Lambda }_{\mathrm{QCD}}\right)}^{2n}$ terms	1.8	1.9	2.3	1.5	1.8	0.1
$m_l$ extrapolation	0.4	0.4	0.7	0.5	0.4	1.9
Total	4.3	5.3	7.0	4.6	4.1	2.5

Table 8. Table 8: Simulation results for ⟨ B q | O n | B ¯ q ⟩ latt / ( f B q ​ M B q ) 2 subscript quantum-operator-product subscript 𝐵 𝑞 subscript 𝑂 𝑛 subscript ¯ 𝐵 𝑞 latt superscript subscript 𝑓 subscript 𝐵 𝑞 subscript 𝑀 subscript 𝐵 𝑞 2 \langle{B}_{q}|O_{n}|\overline{B}_{q}\rangle_{\mathrm{latt}}/(f_{B_{q}}M_{B_{q}})^{2} for B q = B s , B d subscript 𝐵 𝑞 subscript 𝐵 𝑠 subscript 𝐵 𝑑 B_{q}=B_{s},B_{d} mesons. Results are presented for each of the configuration data sets described in Table  3 . Fit ranges for 2-point ( t min ≤ t ≤ t max subscript 𝑡 min 𝑡 subscript 𝑡 max t_{\mathrm{min}}\leq t\leq t_{\mathrm{max}} ) and 3-point ( t min ≤ t ≤ T − t min subscript 𝑡 min 𝑡 𝑇 subscript 𝑡 min t_{\mathrm{min}}\leq t\leq T-t_{\mathrm{min}} ) correlators are tabulated. Sample χ 2 superscript 𝜒 2 \chi^{2} s per degree of freedom from fits to both sets of correlators, with SVD and prior noise (see Appendix  D.4 ), are also listed.

meson	set	$t_{\min }$	$t_{\max }$	$T$	$⟨O_1⟩/{\left(fM\right)}^2$	$⟨O_2⟩/{\left(fM\right)}^2$	$⟨O_3⟩/{\left(fM\right)}^2$	$⟨O_4⟩/{\left(fM\right)}^2$	$⟨O_5⟩/{\left(fM\right)}^2$	${\chi }^2/\mathrm{dof}$ [dof]
$B_s$	1	4	17	8–12	$2.274\left(29\right)$	$-1.887\left(24\right)$	$0.3646\left(68\right)$	$3.041\left(35\right)$	$1.870\left(23\right)$	0.97 [258]
$B_s$	2	4	17	8–12	$2.336\left(29\right)$	$-1.939\left(21\right)$	$0.3625\left(75\right)$	$3.199\left(33\right)$	$1.978\left(22\right)$	1.11 [258]
$B_s$	3	4	17	8–12	$2.315\left(31\right)$	$-1.923\left(26\right)$	$0.3638\left(73\right)$	$3.162\left(39\right)$	$1.951\left(24\right)$	1.01 [258]
$B_s$	4	4	22	10–14, 17	$2.367\left(27\right)$	$-1.979\left(21\right)$	$0.3720\left(70\right)$	$3.382\left(34\right)$	$2.060\left(22\right)$	0.92 [324]
$B_s$	5	4	22	10–14, 17	$2.319\left(26\right)$	$-1.955\left(20\right)$	$0.3636\left(70\right)$	$3.288\left(34\right)$	$1.994\left(22\right)$	1.13 [324]
$B_s$	6	4	22	10–14, 17	$2.365\left(32\right)$	$-1.999\left(25\right)$	$0.3747\left(80\right)$	$3.366\left(40\right)$	$2.037\left(25\right)$	0.99 [324]
$B_s$	7	5	33	12–16, 19	$2.312\left(32\right)$	$-2.029\left(29\right)$	$0.3835\left(83\right)$	$3.548\left(47\right)$	$2.090\left(28\right)$	1.03 [399]
$B_d$	1	4	17	8–12	$2.238\left(56\right)$	$-1.869\left(39\right)$	$0.362\left(14\right)$	$3.054\left(59\right)$	$1.861\left(41\right)$	1.00 [258]
$B_d$	2	4	17	8–12	$2.283\left(70\right)$	$-1.875\left(43\right)$	$0.345\left(19\right)$	$3.232\left(69\right)$	$2.016\left(49\right)$	1.05 [258]
$B_d$	3	4	17	8–12	$2.228\left(74\right)$	$-1.854\left(49\right)$	$0.324\left(19\right)$	$3.232\left(77\right)$	$1.955\left(53\right)$	1.07 [258]
$B_d$	4	4	22	10–14, 17	$2.377\left(56\right)$	$-1.941\left(38\right)$	$0.333\left(14\right)$	$3.539\left(67\right)$	$2.158\left(44\right)$	1.06 [324]
$B_d$	5	4	22	10–14, 17	$2.332\left(59\right)$	$-1.870\left(38\right)$	$0.328\left(15\right)$	$3.341\left(61\right)$	$2.011\left(43\right)$	0.99 [324]
$B_d$	6	4	22	10–14, 17	$2.265\left(97\right)$	$-1.965\left(67\right)$	$0.379\left(25\right)$	$3.49\left(10\right)$	$2.086\left(71\right)$	1.08 [324]
$B_d$	7	5	33	12–16, 19	$2.332\left(70\right)$	$-1.984\left(68\right)$	$0.356\left(18\right)$	$3.68\left(10\right)$	$2.153\left(63\right)$	1.08 [399]

Table 9. Table 9: Sample error budgets from simulations on the 0.09 fm lattice (set 7) for ⟨ B q | O n | B q ⟩ / ( f B q ​ M B q ) 2 quantum-operator-product subscript 𝐵 𝑞 subscript 𝑂 𝑛 subscript 𝐵 𝑞 superscript subscript 𝑓 subscript 𝐵 𝑞 subscript 𝑀 subscript 𝐵 𝑞 2 \langle B_{q}|O_{n}|B_{q}\rangle/(f_{B_{q}}M_{B_{q}})^{2} for B s subscript 𝐵 𝑠 B_{s} and B d subscript 𝐵 𝑑 B_{d} mesons. Percentage errors coming from Monte Carlo statistics, the fit priors, and the SVD cuts are shown; these are added in quadrature to give the total error.

$B_s$	$⟨O_1⟩/{\left(fM\right)}^2$	$⟨O_2⟩/{\left(fM\right)}^2$	$⟨O_3⟩/{\left(fM\right)}^2$
statistics	1.25	1.15	2.00
prior	0.28	0.31	0.47
SVD	0.42	0.74	0.56
total	1.35	1.41	2.13

Table 10. Table 10: Correlations in the uncertainties of the R s n superscript subscript 𝑅 𝑠 𝑛 R_{s}^{n} (Eq. ( 51 )) for different values of  n 𝑛 n . Correlations are also shown for the bag parameters B B s ( n ) superscript subscript 𝐵 subscript 𝐵 𝑠 𝑛 B_{B_{s}}^{(n)} (Eq. ( 2 )), and for the ratio R s n / R d n subscript superscript 𝑅 𝑛 𝑠 subscript superscript 𝑅 𝑛 𝑑 R^{n}_{s}/R^{n}_{d} . Correlations for B B s ( n ) / B B d ( n ) superscript subscript 𝐵 subscript 𝐵 𝑠 𝑛 superscript subscript 𝐵 subscript 𝐵 𝑑 𝑛 B_{B_{s}}^{(n)}/B_{B_{d}}^{(n)} are the same as for R s n / R d n subscript superscript 𝑅 𝑛 𝑠 subscript superscript 𝑅 𝑛 𝑑 R^{n}_{s}/R^{n}_{d} .

	$R_s^1$	$R_s^2$	$R_s^3$	$R_s^4$	$R_s^5$
$R_s^1$	$1.000$	$-0.069$	$0.013$	$0.041$	$0.040$
$R_s^2$		$1.000$	$-0.039$	$-0.040$	$-0.026$
$R_s^3$			$1.000$	$0.023$	$0.012$
$R_s^4$				$1.000$	$0.144$
$R_s^5$					$1.000$
	$B_{B_s}^{\left(1\right)}$	$B_{B_s}^{\left(2\right)}$	$B_{B_s}^{\left(3\right)}$	$B_{B_s}^{\left(4\right)}$	$B_{B_s}^{\left(5\right)}$
$B_{B_s}^{\left(1\right)}$	$1.000$	$0.062$	$0.012$	$0.037$	$0.038$
$B_{B_s}^{\left(2\right)}$		$1.000$	$0.177$	$0.233$	$0.150$
$B_{B_s}^{\left(3\right)}$			$1.000$	$0.170$	$0.107$
$B_{B_s}^{\left(4\right)}$				$1.000$	$0.256$
$B_{B_s}^{\left(5\right)}$					$1.000$
	$R_s^1/R_d^1$	$R_s^2/R_d^2$	$R_s^3/R_d^3$	$R_s^4/R_d^4$	$R_s^5/R_d^5$
$R_s^1/R_d^1$	$1.000$	$0.296$	$-0.034$	$0.064$	$0.047$
$R_s^2/R_d^2$		$1.000$	$0.144$	$0.068$	$0.045$
$R_s^3/R_d^3$			$1.000$	$0.035$	$0.018$
$R_s^4/R_d^4$				$1.000$	$0.350$
$R_s^5/R_d^5$					$1.000$

Table 11. Table 11: Means and standard deviations for the reduced matrix elements R s n subscript superscript 𝑅 𝑛 𝑠 R^{n}_{s} and bag parameters B s ( n ) subscript superscript 𝐵 𝑛 𝑠 B^{(n)}_{s} for each 4-quark operator O n subscript 𝑂 𝑛 O_{n} , together with values for their B s / B d subscript 𝐵 𝑠 subscript 𝐵 𝑑 B_{s}/B_{d}  ratios.

	$R_s^n$	$R_s^n/R_d^n$	$B_s^{\left(n\right)}$	$B_s^{\left(n\right)}/B_d^{\left(n\right)}$
$O_1$	$2.1678\left(928\right)$	$1.0081\left(250\right)$	$0.8129\left(348\right)$	$1.0081\left(250\right)$
$O_2$	$-2.1801\left(1035\right)$	$1.0589\left(242\right)$	$0.8169\left(431\right)$	$1.0626\left(243\right)$
$O_3$	$0.4357\left(288\right)$	$1.0886\left(339\right)$	$0.8163\left(572\right)$	$1.0924\left(340\right)$
$O_4$	$3.6532\left(1480\right)$	$0.9558\left(213\right)$	$1.0332\left(471\right)$	$0.9589\left(214\right)$
$O_5$	$1.9448\left(759\right)$	$0.9650\left(232\right)$	$0.9406\left(384\right)$	$0.9668\left(233\right)$

Equations

$$O_1$$

$$O_2$$

$$O_3$$

$$O_4$$

$$O_5$$

$$B_{B_q}^{(i)}(\mu )\equiv \frac{⟨B_q|O_i^q|{\overline{B}}_q{⟩}_{\overline{\mathrm{MS}}}^{(\mu )}}{{\eta }_i^q(\mu )f_{B_q}^2{M}_{B_q}^2},$$

$$⟨0|{\overline{\Psi }}_q^i{\gamma }_0{\gamma }_5{\Psi }_b^i|B_q(\vec{p}=0)⟩=f_{B_q}{M}_{B_q}.$$

$${\eta }_1^q$$

$${\eta }_2^q$$

$${\eta }_3^q$$

$${\eta }_4^q$$

$${\eta }_5^q$$

$$⟨O_i{⟩}_{\overline{\mathrm{MS}}}^{(m_b)}$$

$$+\sum _{j\ne i}{\alpha }_s{z}_{ij}⟨O_j{⟩}_{\mathrm{latt}}$$

$$\displaystyle\quad+\mathcal{O}\Big{(}\alpha_{s}^{2},\frac{\alpha_{s}\Lambda_{\mathrm{QCD}}}{m_{b}},\alpha_{s}(a\Lambda_{\mathrm{QCD}})^{2}\Big{)}\,.$$

$$[{\overline{\Psi }}_b{\Gamma }_1{\Psi }_q][{\overline{\Psi }}_b{\Gamma }_2$$

$$\displaystyle+\frac{1}{2m_{b}}\big{[}\bm{\nabla}\overline{\Psi}_{b}^{\mathrm{NR}}\cdot\bm{\gamma}\Gamma_{1}\Psi_{q}\big{]}\big{[}\overline{\Psi}_{\overline{b}}^{\mathrm{NR}}\Gamma_{2}\Psi_{q}\big{]}$$

$$\displaystyle+\frac{1}{2m_{b}}\big{[}\overline{\Psi}_{b}^{\mathrm{NR}}\Gamma_{1}\Psi_{q}\big{]}\big{[}\bm{\nabla}\overline{\Psi}_{\overline{b}}^{\mathrm{NR}}\cdot\bm{\gamma}\Gamma_{2}\Psi_{q}\big{]}$$

$$\displaystyle+\big{(}\Gamma_{1}\leftrightarrow\Gamma_{2}\big{)},$$

$$\displaystyle{{O}}_{2}\Big{|}_{\mathrm{BJU}}$$

$$\displaystyle{{O}}_{3}\Big{|}_{\mathrm{BJU}}$$

$$\alpha_{s}\Big{(}c_{\alpha_{s}}^{n,0}+c_{\alpha_{s}}^{n,1}\delta_{a}+c_{\alpha_{s}}^{n,2}\delta^{2}_{a}\Big{)}\frac{\langle O_{n}\rangle_{\overline{\mathrm{MS}}}}{(fM)^{2}}$$

$${\delta }_a\equiv (a{m}_b-2.6)/1.4$$

$$\displaystyle\alpha_{s}\big{(}a\Lambda_{\mathrm{QCD}}\big{)}^{2}\Big{(}c_{a^{2}}^{n,0}+c_{a^{2}}^{n,1}\delta_{a}+c_{a^{2}}^{n,2}\delta^{2}_{a}\Big{)}\frac{\langle O_{n}\rangle_{\overline{\mathrm{MS}}}}{(fM)^{2}}$$

$$\displaystyle+\big{(}a\Lambda_{\mathrm{QCD}}\big{)}^{4}\Big{(}c_{a^{4}}^{n,0}+c_{a^{4}}^{n,1}\delta_{a}+c_{a^{4}}^{n,2}\delta^{2}_{a}\Big{)}\frac{\langle O_{n}\rangle_{\overline{\mathrm{MS}}}}{(fM)^{2}}$$

$${\hat{B}}_{B_q}^{(1)}=c_{\text{RGI}}{B}_{B_q}^{(1)}(m_b).$$

$${\hat{B}}_{B_s}^{(1)}$$

$${\hat{B}}_{B_d}^{(1)}$$

$$\frac{{\hat{B}}_{B_s}^{(1)}}{{\hat{B}}_{B_d}^{(1)}}$$

$$f_{B_s}\sqrt{{\hat{B}}_{B_s}^{(1)}}$$

$$f_{B_d}\sqrt{{\hat{B}}_{B_d}^{(1)}}$$

$$\xi$$

$$\Delta M_{s,\text{expt}}=17.757(21){\text{ps}}^{-1}$$

$$\Delta M_{d,\text{expt}}=0.5065(19){\text{ps}}^{-1}$$

$$\Delta M_q=\frac{G_F^2{M}_W^2{M}_{B_q}}{6{\pi }^2}{S}_0(x_t){\eta }_{2B}{\left|V_{tq}^{\ast }{V}_{tb}\right|}^2{f}_{B_q}^2{\hat{B}}_{B_q}^{(1)}.$$

$${\left|V_{ts}\right|}_{\text{CKMfitter, tree}}$$

$${\left|V_{td}\right|}_{\text{CKMfitter, tree}}$$

$${\left|V_{td}/V_{ts}\right|}_{\text{CKMfitter, tree}}$$

$${\left|V_{tb}\right|}_{\text{CKMfitter, tree}}$$

$$\Delta M_{s,\text{SM}}$$

$$\Delta M_{d,\text{SM}}$$

$${\left(\frac{\Delta M_d}{\Delta M_s}\right)}_{\text{SM}}$$

$$\left|V_{td}\right|$$

Full text

I Introduction

The Standard Model description of neutral $B_{d}$ and $B_{s}$ oscillations requires knowledge of hadronic parameters derived from the matrix elements of 4-quark operators between $B_{q}$ and $\overline{B}_{q}$ states. These 4-quark operators come from the effective electroweak Lagrangian at energy scales appropriate to $B$ physics and the matrix elements can only be determined by lattice QCD calculations, which are now able to include the full impact of QCD on such hadronic quantities Davies et al. (2004). The accuracy with which this can be done is the limiting factor in the constraint that can be obtained from the now very precise experimental results on the neutral meson mass difference (seen as an oscillation frequency). In the Standard Model this constraint leads to a determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements that accompany the 4-quark operators of the Standard Model. New physics models with extra heavy particles extend the effective Hamiltonian to include additional 4-quark operators. Constraints on the new physics from experiment then need accurate determination of the matrix elements of the new operators. Again this can come only from lattice QCD calculations.

Here we provide the first “second-generation” lattice QCD calculation of the matrix elements of all five $\Delta B$ = 2 operators of dimension six for the $B_{s}$ and $B_{d}$. We improve on earlier calculations by working on gluon field configurations generated by the MILC collaboration that include $u$, $d$, $s$ and $c$ quarks in the sea with $u/d$ quark masses going down to their physical values. Although this obviates the need for a chiral extrapolation, we also include heavier $u/d$ quark masses in our set of results so that we can map out the dependence on the light quark mass. The discretisation of QCD that we use is fully improved through $\mathcal{O}(\alpha_{s}a^{2})$ for both the gluons and the light quarks (including all of those in the sea) for which the Highly Improved Staggered Quark (HISQ) action is used. For the $b$ quarks we use improved NonRelativistic QCD, which includes $\mathcal{O}(\alpha_{s})$ corrections to terms at order $v^{4}$ (where $v$ is the heavy quark velocity). By linking this calculation directly to our earlier one for the $B$-meson decay constants Dowdall et al. (2013a) (that parameterize the amplitude to create a meson from the vacuum) we are able to give results directly for the “bag parameters” associated with each operator and take advantage of the cancellation of a number of systematic effects. The bag parameters (to be defined in Section II.1) encode the multiplicative factor by which the operator matrix element differs from that expected in the vacuum saturation approximation, which is related to the decay constant. To the extent to which this approximation works (and we will show here that it does work well) we expect the bag parameters to have very little dependence on light quark sea or valence masses and even on the lattice spacing. This enables improvements in accuracy over earlier work along with a much simpler picture of the extrapolation to the physical point.

The first unquenched lattice QCD calculations of the matrix elements for neutral $B$ meson mixing focussed purely on results for the Standard Model operators Dalgic et al. (2007); Gamiz et al. (2009) and ratios for $B_{s}$ to $B_{d}$ Bazavov et al. (2012a). Calculations have also been done in the infinite heavy quark mass limit Aoki et al. (2015). More recently calculations of matrix elements for the full set of SM and BSM operators have been done Carrasco et al. (2014); Bazavov et al. (2016). The calculations in Carrasco et al. (2014) use the twisted mass formalism for all quarks on gluon field configurations including $u$ and $d$ quarks in the sea. An extrapolation of results (renormalised using the RI-MOM scheme) is made from a heavy quark mass in the charm region up to the $b$ quark mass using ratios with a known infinite mass limit. The calculations in Bazavov et al. (2016) use the Fermilab formalism for the $b$ quark and the asqtad formalism for the light quarks on gluon field configurations that include $u$, $d$ and $s$ quarks in the sea with the asqtad formalism. Perturbatively renormalised 4-quark operator matrix elements (only) are calculated and so bag parameters must be derived using decay constant results from elsewhere. A very recent result in Boyle et al. (2018) uses domain-wall quarks on gluon field configurations including $u$, $d$ and $s$ in the sea and extrapolates in heavy quark mass to the $b$ quark mass for $B_{s}$ to $B_{d}$ ratios for SM mixing matrix elements and decay constants.

In Section II we discuss the 4-quark operators relevant to $B$ mixing and how they are implemented on the lattice. Section III describes our lattice calculation and results. We compare our results to previous work in Section IV, determine CKM elements $V_{ts}$ and $V_{td}$ using experimental results on $B$-meson mass differences, determine branching fractions for the rare decays of $B_{d}$ and $B_{s}$ to $\mu^{+}\mu^{-}$, and give matrix elements for derived operators that contribute to width differences. Section V gives our conclusions and discusses the prospects for future improvements. Details about our analysis are contained in four Appendices: Appendix A on the lattice QCD correlators that we calculate and how we fit them to obtain matrix elements and bag parameters; Appendix B on the chiral perturbation theory fits we use to combine results at physical and unphysical light quark masses; Appendix C on correlations in the uncertainties for our final results and lastly, with more general applications beyond this analysis, Appendix D on SVD cuts and fitting correlators.

II Background

II.1 Continuum 4-quark operators

Neutral $B$-meson mixing occurs at lowest order in the Standard Model through box diagrams involving the exchange of $W$ bosons and top quarks, see Figure 1. These box diagrams can be well approximated by an effective Lagrangian expressed in terms of 4-quark operators. Here we will examine all five of the independent local dimension-6 operators that could contribute to $\Delta B=2$ processes Gabbiani et al. (1996); Bazavov et al. (2016):

[TABLE]

where $V=\gamma_{\mu}$, $A=\gamma_{\mu}\gamma_{5}$, $S=1$ and $P=\gamma_{5}$, and sums over $\mu$ and color indices $i$ and $j$ are implicit. In the Standard Model, the most important of these for $B$-$\overline{B}$ mixing is $O_{1}$. This operator mixes with $O_{2}$ under renormalization. Operators $O_{4}$ and $O_{5}$ do not appear in the Standard Model, but do arise in various BSM scenarios.

It is conventional to parameterize matrix elements of these operators in terms of “bag parameters,”

[TABLE]

where here $\mu$ is the renormalization scale, and $M_{B_{q}}$ and $f_{B_{q}}$ are the mass and weak decay constant of the $B_{q}$ meson:

[TABLE]

The normalization parameter $\eta_{i}^{q}(\mu)$ is chosen so that the bag parameters equal 1 in the “vacuum saturation approximation,” where gluon (and other QCD) exchanges between the initial and final $\overline{B}_{q}$ and ${B}_{q}$ are ignored (see Gabbiani et al. (1996); Bazavov et al. (2016) for more details):

[TABLE]

We use renormalization scale $\mu=m_{b}(m_{b})$; the corresponding values for the normalization factors are given in Table 1.

The bag parameters provide both computational advantages and physical insights. The leading-order logarithms in chiral perturbation theory, coming from the matrix element of the 4-quark operator and $f_{B_{q}}^{2}$, partly cancel in the ratio; see Appendix B. In particular the coefficient of the chiral logarithm from the tadpole diagrams is reduced by a factor of 4. Therefore bag parameters should be less dependent upon the light-quark mass; we find very little mass dependence. Finite-volume effects will be correspondingly reduced. We also find that most of the dependence on lattice spacing cancels. Finally, as we will show, the bag parameters all turn out to be of order one, suggesting that vacuum saturation is a useful approximation. For these reasons, we focus here on bag parameters; values for the matrix elements are easily obtained from the bag parameters given values for the decay constants Dowdall et al. (2013a); Bazavov et al. (2018).

II.2 Lattice QCD 4-quark operators and matching

Matrix elements of the 4-quark operators are regulator dependent, and so we need to convert matrix elements calculated in our simulation (with the lattice regulator) into the corresponding matrix elements for the more conventional ${\overline{\mathrm{MS}}}$ scheme. The differences between the two schemes are ultraviolet and so can be calculated using QCD perturbation theory. To lowest and first order in $\alpha_{s}$ the relationship has the form (for $\mu=m_{b}$):

[TABLE]

The coefficients $z_{ij}$ relevant to our simulation were calculated in Monahan et al. (2014) and are summarized in Table 2. The scale for $\alpha_{s}$ depends on the lattice spacing; we use the same values for $\alpha_{s}$ used in Colquhoun et al. (2015) to calculate renormalizations for the axial-vector current that couples to $B_{q}$ mesons (see Table 4 for the values).

Our lattice analysis uses non-relativistic QCD (NRQCD) for the $b$ dynamics. Quarks and anti-quarks decouple in NRQCD, and so correspond to separate fields. As a result the lattice version of a 4-quark operator has the form Monahan et al. (2014)

[TABLE]

where $\overline{\Psi}_{b}^{\mathrm{NR}}$ creates a $b$ quark and $\overline{\Psi}_{\overline{b}}^{\mathrm{NR}}$ annihilates a $b$ anti-quark. These are the lattice operators we use on the right-hand side of Eq. (II.2). The $1/m_{b}$ terms are the $\mathcal{O}(\Lambda_{\mathrm{QCD}}/m_{b})$ corrections to the operator in NRQCD.

A complication for operators $O_{2}$ and $O_{3}$ is the treatment of “evanescent operators.” Our matching results use the ${\overline{\mathrm{MS}}}_{\mathrm{NDR}}$ scheme of Beneke et al. (1999) (BBGLN). These matrix elements are readily converted to the alternative scheme of Buras et al. (2001) (BJU) using the following equations (through $\mathcal{O}(\alpha_{s})$, with $\mu=m_{b}$) Becirevic et al. (2002); Ciuchini et al. (2003):

[TABLE]

The matching coefficients $z_{ij}$ from Eq. (II.2) are plotted against $am_{b}$ in Figure 2. These coefficients are not large and have a relatively benign dependence on the $b$-quark mass across the range that we use, although different coefficients behave differently. The diagonal coefficients $z_{ii}$ are much larger than the corresponding coefficients for the NRQCD-HISQ axial current ($z_{A_{0}}$ in Table 2), which are unusually small. Note that the only nonzero off-diagonal coefficients are for $ij$ equal to 12, 21, 31, 45, and 54, and that these tend to be smaller than the diagonal parameters.

It is worth remarking here on the similarities and differences between the perturbative matching we apply here and that used by the Fermilab/MILC collaborations Bazavov et al. (2016) in their determination of $B$ mixing matrix elements. They also make use of a perturbative calculation of the matching to $\mathcal{O}(\alpha_{s})$. They do this after a non-perturbative determination of factors that are needed to remove normalisation artifacts from the clover and asqtad actions that they use for heavy and light quarks respectively and without which they would have large $\mathcal{O}(\alpha_{s})$ coefficients. We do not need to apply this procedure because the NRQCD and HISQ actions are well-behaved in this respect Chakraborty et al. (2017). After applying their nonperturbative procedure, the Fermilab/MILC collaboration give results for their $\mathcal{O}(\alpha_{s})$ coefficients in Table III of Bazavov et al. (2016). Their coefficients differ from ours because they are using a different discretization of QCD for both heavy and light quarks. However qualitatively the coefficients show similar behavior in terms of magnitude and dependence on the lattice $b$ quark mass (given in their case by the parameter $\kappa_{b}^{\prime}$).

III Lattice Calculation

III.1 Simulations

We use seven ensembles of gluon field configurations recently generated by the MILC collaboration Bazavov et al. (2010, 2013). Details are given in Table 3. We use ensembles at three values of the lattice spacing, $a$, to control discretisation effects and at three values of the light quark mass down to the physical point to map out sea quark mass effects. Discretisation effects depend on $a^{2}$ and sea quark mass effects are approximately linear, so a range in $a^{2}$ of a factor of 3 and in sea light quark mass of a factor of 5 allows us substantial leverage to pin down these effects. The lattice spacing values were determined using the mass splitting between the $\Upsilon$ and $\Upsilon^{\prime}$, as described in Dowdall et al. (2012a) where a discussion of systematic errors can be found. The sea quarks use HPQCD’s HISQ action Follana et al. (2007) which we have shown to have small discretisation errors even for charm quarks Follana et al. (2008); Davies et al. (2010a); Donald et al. (2012). This enables four flavors of quarks to be included in the sea, with masses given in Table 3. The $u$ and $d$ quark masses are taken to be the same.

The valence $b$ quarks are implemented using lattice NonRelativistic QCD (NRQCD) Lepage et al. (1992). The action is described in detail in Dowdall et al. (2012a). It includes a number of improvements over earlier calculations, in particular one-loop radiative corrections (beyond tadpole-improvement Lepage and Mackenzie (1993)) to most of the coefficients of the $\mathcal{O}(v_{b}^{4})$ relativistic correction terms. The tadpole-improvement of the action is done using the Landau gauge-link, with $u_{0L}$ values given in Table 4. This action has been shown to give excellent agreement with experiment in recent calculations of the bottomonium Dowdall et al. (2012a); Daldrop et al. (2012); Dowdall et al. (2014) and $B$-meson spectra Dowdall et al. (2012b). The $b$ quark mass is tuned, giving the values in Table 4, by fixing the spin-averaged kinetic mass of the $\Upsilon$ and $\eta_{b}$ states to experiment Dowdall et al. (2012a). NRQCD breaks down as $am_{b}\rightarrow 0$ but all our values of $am_{b}$ are substantially larger than 1, where there is no problem.

The HISQ valence light quark masses are taken to be equal to the sea mass except on set 4 where there is a slight discrepancy. The $s$ quark is tuned using the mass of the $\eta_{s}$ meson Davies et al. (2010b), a fictitious pseudoscalar $s\overline{s}$ state which is not allowed to decay on the lattice. Its properties can be very accurately determined in lattice QCD and we find $M_{\eta_{s}}=$ 0.6885(22) GeV Dowdall et al. (2013b). Values for valence $s$ masses are given in Table 4 and corresponding values of $M_{\eta_{s}}$ in lattice units in Dowdall et al. (2013b). We allow for uncertainties from mistuned valence masses in our determination of the physical results.

III.2 Simulation Results and Error Budget

We describe the 2-point and 3-point correlators used in our analysis in Appendix A. We use the 2-point correlators to extract the decay constants $f_{B_{q}}$, including the $1/m_{b}$ corrections (Eq. (34)). We also combine them with the 3-point correlators to calculate lattice matrix elements of the $O_{n}$. Also in that Appendix, we discuss the Bayesian fits used to extract physics from these correlators. Our final results for $\langle{B}_{q}|O_{n}|\overline{B}_{q}\rangle_{\mathrm{latt}}/(f_{B_{q}}M_{B_{q}})^{2}$ are summarized in Table 8 of Appendix A. As discussed in the Appendix, this was a blind analysis.

We convert the lattice expectation values into ${\overline{\mathrm{MS}}}$ matrix elements using Eq. (II.2) (divided by $(f_{B_{q}}M_{B_{q}})^{2}$). Our results are listed in Table 5. In addition to the statistical errors from the simulation and the (negligible) errors in the $z_{ij}$, we include uncertainties (for each entry in the table) coming from three additional sources:

• •

$\mathcal{O}(\alpha_{s}^{2})$: We estimate this uncertainty to be twice (to be conservative) $\alpha_{s}$ times the magnitude of the $\mathcal{O}(\alpha_{s})$ correction we include for each of our three lattice spacings. (These corrections are correlated between configuration sets with similar lattice spacings.)

• •

$\mathcal{O}({\alpha_{s}\Lambda_{\mathrm{QCD}}/m_{b}})$: $\mathcal{O}(\Lambda_{\mathrm{QCD}}/m_{b})$ corrections have been measured for the temporal axial-vector current and found to be 5% of the leading-order contribution Dowdall et al. (2013a). This suggests that $\mathcal{O}(\Lambda_{\mathrm{QCD}}/m_{b})$ corrections, which are included in our simulations, are 10% for the 4-quark operators. We account for the $\mathcal{O}(\alpha_{s})$ radiative corrections to these terms by adding the following uncertainty to our results:

[TABLE]

where each $c_{\alpha_{s}}^{n,i}=0\pm 0.1$ and

[TABLE]

allows for variation in the coefficients between the lattice spacings. ($\delta_{a}$ is defined to vary from $-1/2$ to $1/2$ over our mass range; see Dowdall et al. (2012a) for more details.)

• •

\mathcal{O}\big{(}\alpha_{s}(a\Lambda_{\mathrm{QCD}})^{2},(a\Lambda_{\mathrm{QCD}})^{4}\big{)}: The NRQCD and HISQ actions we use in the simulation are highly corrected. In particular there are no tree-level $a^{2}$ errors in either. We account for $a^{2}\alpha_{s}$ errors by adding the following uncertainty to our results:

[TABLE]

where each $c_{a^{2}}^{n,i}=0\pm 1$, each $c_{a^{4}}^{n,i}=0\pm 1$, $\Lambda_{\mathrm{QCD}}=0.5$ GeV is the QCD scale, and again the $\delta_{a}$ terms allow for variation between different lattice spacings.

The final entries (“phys.”) in both parts of Table 5 are our final results at the physical values of the light-quark masses after our chiral fit. We use chiral perturbation theory to combine the values obtained on the different configuration sets with different light quark masses; see Appendix B for details. Figure 3 compares our final values for $B_{s}$ mesons with the results from individual configuration sets. These plots show that the dependence on lattice spacing and light-quark mass is negligible compared with our uncertainties. The analogous plot for $B_{d}$ mesons is very similar.

Adding uncertainties to the lattice results to allow for operator normalisation and lattice spacing effects, as we have done above, is equivalent to including them in our fit function with the coefficients treated as fit parameters; see the Appendix of McNeile et al. (2010). The uncertainties that are included are correlated between lattice results on different sets through these coefficients. Figure 3 shows that, for example, the lattice spacing effects that we allow for through Eq. (10) are overestimates of what is seen in the results, since the variation with lattice spacing of the central values is much smaller than the uncertainties shown on the coarser lattices. A further test of this is that omitting the results from the smallest lattice spacing (set 7) shifts our final central values by less than half a standard deviation and often much less.

Finally we convert our final results into bag parameters using Eq. (2). The bag parameters are listed in Table 6. Despite the wide variation in values for $\langle O_{n}\rangle/(fM)^{2}$, the bag parameters are within 30% of 1. This shows that the vacuum saturation approximation can be of some utility.

Figure 4 compares our final results for ratios of bag parameters $B_{B_{s}}^{(n)}/B_{B_{d}}^{(n)}$ with results from the different configuration sets. Results are plotted versus the value of $m_{\pi}^{2}$ used in each simulations. Again there is very little variation with quark mass, with all ratios within 5% of 1. Our final results are shifted by less than half a standard deviation if we omit the data with the largest pion masses, and have errors that are 10-15% larger.

The error budgets for the $B_{s}$ bag parameters are shown in Table 7. The dominant source of error comes from uncalculated terms in perturbation theory ($\alpha_{s}^{2}$ and $\alpha_{s}\Lambda_{\mathrm{QCD}}/m_{b}$ terms). The sensitivity to these terms depends on the operator. For example, it is particularly high for $O_{3}$, because matrix elements for $O_{3}$ are a lot smaller than those of $O_{1}$ (see Eq. 4) which are mixed in by Eq. (II.2). The error budgets for $B_{d}$ mesons are almost identical to those for $B_{s}$, but have twice as large a contribution from statistical uncertainties in the lattice data. Almost all of the uncertainties, and some of the statistical errors, cancel in ratios of $B_{s}$ to $B_{d}$ meson bag parameters.

Matrix elements of the 4-quark mixing operators can be obtained from the ratios in Table 5 given values for the decay constants and masses. Note that the corresponding bag parameters for $O_{2\ldots 5}$ have larger fractional errors than the ratios, and so should not be used for this purpose. The larger errors result from uncertainties due to the factors $\eta^{q}_{i}$ in the bag-parameter definition (see Table 7 and Eq. (2)).

IV Discussion

IV.1 Comparison to previous results

Our results for the bag parameters for all five SM and BSM operators given in Table 6 are more accurate than previous lattice QCD results. This is for a number of reasons:

• •

We work directly with the bag parameters rather than the 4-quark operator matrix elements. The bag parameters are expected from chiral perturbation theory to have little dependence on valence and sea quark masses (see Appendix B). This expectation is borne out in our results and means that we are able easily to combine results at both unphysical and physical light quark masses.

• •

We have results for the physical light quark mass at two values of the lattice spacing improving control of the chiral extrapolation.

• •

The gluon field configurations that we use include the effect of $u$, $d$, $s$ and $c$ quarks in the sea and so we do not have an uncertainty associated with missing flavours of sea quarks (the Fermilab/MILC collaboration include a 2% uncertainty in their 4-quark operator matrix elements from missing $c$ in the sea Bazavov et al. (2016)).

Figure 5 shows a comparison of our bag parameters for the $B_{s}$ meson to those from Bazavov et al. (2016) and Carrasco et al. (2014) (and also, for $O_{1}$, to Gamiz et al. (2009)). The results from Carrasco et al. (2014) include only $u$ and $d$ quarks in the sea and the uncertainty does not include an estimate of the impact of missing $s$ sea quarks. It is therefore not clear whether we should expect agreement between these $n_{f}=2$ results and our $n_{f}=4$ results. The fact that the $n_{f}=2$ purple diamonds from the ETM collaboration are around 20% below our results for $O_{4}$ and $O_{5}$ is reminiscent of what is seen in kaon mixing. ETM use the RI-MOM renormalisation scheme for the purple diamonds and it has been shown in kaon mixing Garron et al. (2016) that the use of the RI-MOM scheme (rather than RI-SMOM) for the equivalent 4-quark operators has large systematic errors that push down the value of the bag parameter. This may then be the main reason (rather than a difference of $n_{f}$) for the discrepancy with our results for $O_{4}$ and $O_{5}$, but more work would be needed to be sure of this.

The $n_{f}=3$ and $n_{f}=4$ results should be comparable because the impact of missing $c$ quarks in the sea on the bag parameters is expected to be very small Bazavov et al. (2016). Our new results agree within $2\sigma$ in each case with Bazavov et al. (2016) but in every case are more accurate. The largest discrepancy is for $B^{(1)}_{B_{s}}$ at $1.9\sigma$.

The weighted average of our $n_{f}=4$ results and the $n_{f}=3$ results from Bazavov et al. (2016) is given by the grey band in the Figure and the value of that average is given in each panel. We assume no correlations, here and subsequently, between our results and those of Bazavov et al. (2016) because they use different actions for both the $b$ quark and the light quarks and different gluon field configurations (with a different sea quark action and generated with a different Monte Carlo updating algorithm).

Figure 6 shows a comparison of the ratio of bag parameters for $B_{s}$ to $B_{d}$ for each operator for our new results and those of Bazavov et al. (2016). Our new results are a lot more accurate, with 2–3% total uncertainty. All of the ratios are very close to 1, but there is a sign of a systematic trend for the ratio for $O_{2}$ and $O_{3}$ to be above 1 and for $O_{4}$ and $O_{5}$ below 1. This is not visible in the results of Bazavov et al. (2016) but does start to emerge with the improved accuracy of our results. This is in general agreement with the results from using sum rules in Grozin et al. (2017); King et al. (2019). We also include in Figure 6 results for $O_{1}$ from HPQCD Gamiz et al. (2009) and RBC/UKQCD Boyle et al. (2018). The RBC/UKQCD result has a 1% uncertainty.

IV.2 Derived quantities

Our results for the bag parameters can be combined with results for the $B$ and $B_{s}$ decay constants to give values for the 4-quark operator matrix elements using Eq. (51) and our results in Table 5. For this we use the most accurate current lattice QCD results obtained on gluon field configurations including $u/d$, $s$ and $c$ quarks in the sea. These have been obtained by the Fermilab/MILC collaboration using the HISQ action for all quarks Bazavov et al. (2018). This ‘heavy-HISQ’ approach, pioneered by HPQCD McNeile et al. (2012a, b), uses pseudoscalar meson 2-point correlators that combine heavy and light quark propagators calculated with multiple heavy quark masses, $am_{h}$, at multiple values of the lattice spacing. $m_{h}$ reaches the $b$ quark mass for $am_{h}<1$ for lattice spacing values $a<0.045$ fm. Since the HISQ action has very small discretisation errors by design, a fit to the $m_{h}$- and $a-$dependence is possible that allows the continuum $m_{h}$-dependence of the decay constant to be reconstructed. It can then be evaluated at the $b$ quark mass to enable the $B$ and $B_{s}$ decay constants to be determined. Note that the correlators can be absolutely normalised in this case and so there is no normalisation uncertainty.

Fermilab/MILC obtain the values $f_{B_{d}}=0.1905(13)$ GeV, $f_{B_{s}}=0.2307(13)$ GeV and $f_{B_{s}}/f_{B_{d}}=1.2109(41)$222Our results obtained on $n_{f}=2+1+1$ gluon field configurations from NRQCD-HISQ calculations Dowdall et al. (2013a); Hughes et al. (2018) agree with these numbers but are less accurate.. Note that we use the decay constant for the neutral $B_{d}$ meson (not the $B_{u}$), which is the appropriate choice here. Our bag parameters are calculated for a light quark $l$ corresponding to the average of $u$ and $d$. Our results show (comparing those for $B_{s}$ with those for $B$) that any difference between bag parameters for $B_{l}$ and $B_{d}$ will be much smaller than our uncertainties. This is not true for the decay constants, where the differences are significant Bazavov et al. (2018).

For the SM phenomenology to be determined from our results for the matrix elements of $\mathcal{O}_{1}$ it is convenient to convert our results from the $\overline{\text{MS}}_{\text{NDR}}$ scheme to the renormalisation-group-invariant quantities $\hat{B}^{(1)}_{B_{q}}$. The conversion is given by

[TABLE]

The matching factor $c_{\text{RGI}}$ is calculated to two-loops in perturbative QCD and we take $c_{\text{RGI}}=1.5158(36)$ Bazavov et al. (2016). This corresponds to the result for $n_{f}=5$ active flavours in the sea and $\overline{\alpha}_{s}(M_{Z})=0.1185(6)$. Our bag parameters are obtained at scale $m_{b}$ for 4 flavours of quarks in the sea. The impact of missing $b$ quarks in the sea, however, should be negligible both for the bag parameters and the resulting 4-quark operator matrix elements. A power-counting estimate of such effects would give a relative contribution of $\alpha_{s}(\Lambda_{\text{QCD}}/2m_{b})^{2}$, which is below 0.1%.

Our results for the RGI bag parameters for $\mathcal{O}_{1}$ are then:

[TABLE]

The ratio of RGI bag parameters is of course the same as that of the $\overline{\text{MS}}$ bag parameters. Combined with the decay constant results from Bazavov et al. (2018) we obtain

[TABLE]

where $\xi$ is the ratio of the two results above it. We form $\xi$ by combining the result for $f_{B_{s}}/f_{B_{d}}$ from Bazavov et al. (2018) with our results for $B^{(1)}_{B_{s}}/B^{(1)}_{B_{d}}$, taking advantage of the correlations that reduce uncertainties in each of these ratios. Note that in combining the decay constant and bag parameter results we add relative uncertainties in quadrature. We expect no significant correlation between the two sets of results because they use a different heavy quark action and, even though both results use $n_{f}=2+1+1$ gluon field configurations, there is little overlap in the ensembles used. The error budgets in the two cases show that the key sources of uncertainty are not the same. The uncertainties in the combinations above are dominated by the uncertainties in our bag parameters and their ratio in Eq. (12) because the decay constant results are now so accurate.

Figure 7 compares our new result for $f_{B_{s}}\sqrt{\hat{B}_{B_{s}}}$ from Eq. (13) to previous lattice QCD results on $n_{f}=3$ gluon field configurations from Fermilab/MILC Bazavov et al. (2016) and HPQCD Gamiz et al. (2009). The Fermilab/MILC results include an uncertainty for missing $c$ in the sea in their calculation. The difference between the central value of our new result and that of Fermilab/MILC is $1.8\sigma$. Because the systematic uncertainties are correlated between our results and those of Gamiz et al. (2009) we do not include the previous HPQCD results in the new lattice QCD $n_{f}=3/n_{f}=4$ average, shown by the grey band in Figure 7. The average value is shown above the grey band.

Figure 8 compares lattice QCD results for the ratio $\xi$ defined in Eq. (13) on $n_{f}=3$ gluon field configurations with our new result here using $n_{f}=4$. There is good agreement between the lattice QCD results with the most recent (including our new result here) having total uncertainties at the level of 1.5%. The result of averaging our new result with that of Bazavov et al. (2016) (both results being obtained at the physical $b$ quark mass) is given by the grey band with the average value quoted above it.

IV.3 $\Delta M$

The phenomenon of neutral $B$-meson oscillations is now well-established experimentally (for recent results see Abe et al. (2005); Aubert et al. (2006); Abazov et al. (2006); Aaij et al. (2016); Abulencia et al. (2006); Aaij et al. (2012, 2013a, 2013b, 2015)), with an oscillation frequency that is set by the mass difference between the two eigenstates. The current experimental average values Tanabashi et al. (2018) for the $B_{s}$ and $B_{d}$ systems are:

[TABLE]

combining statistical and systematic errors in quadrature.

In the SM $\Delta M$ is given by

[TABLE]

Here $S_{0}$ is the Inami-Lim function Inami and Lim (1981) which describes electroweak corrections and has argument $x_{t}=m_{t}^{2}/M_{W}^{2}$. The top quark mass to be used here is in the $\overline{MS}$ scheme, $\overline{m}_{t}(\overline{m}_{t})$ Buchalla et al. (1996). Taking the current average Tanabashi et al. (2018) of direct experimental measurements Khachatryan et al. (2016); Aaboud et al. (2016); Sirunyan et al. (2018) of the top quark mass (172.9(4) GeV) as the pole mass, gives $\overline{m}_{t}(\overline{m}_{t})$ = 163.07(38) GeV using the 4-loop expressions in Marquard et al. (2015). Evaluating the Inami-Lim function then gives: $S_{0}(4.116(19))=2.313(8)$. The QCD correction factor, $\eta_{2B}$, is given at next-to-leading order in Buras et al. (1990). We take $\eta_{2B}$=0.55210(62) Bazavov et al. (2018), again calculated with $n_{f}=5$.

The CKM elements $V_{tq}$ and $V_{tb}$ can be derived in the SM by assuming that the CKM matrix is unitary and determining other CKM elements in the same rows or columns from the comparison of theory and experiment Charles et al. (2005); Bona et al. (2006); CKM (2018); UTf (2018). For Eq. (15) it is important to use values for $V_{tq}$ that did not include $\Delta M_{q}$ itself in their determination. So we use the results from CKMfitter for the case where only tree-level processes were used in the determination. This gives CKM (2018)

[TABLE]

The ratio $\left|V_{td}/V_{ts}\right|_{\text{CKMfitter, tree}}$ is derived from the CKMfitter results for $A$, $\lambda$, $\overline{\rho}$ and $\overline{\eta}$ using the formulae in Charles et al. (2005). The central value differs slightly from the ratio of the two numbers above.

The final terms in Eq. (15) parameterise the hadronic contribution to $\Delta M$ through the matrix element of the appropriate 4-quark operator, $\mathcal{O}_{1}$. Our results for $f_{B_{q}}^{2}\hat{B}^{(1)}_{B_{q}}$ are given in Eq. (13).

Putting all these pieces together we obtain predictions for the mass differences for neutral $B_{s}$ and $B_{d}$ eigenstates of

[TABLE]

where the first error in each case is from the CKM matrix elements and the second error is primarily from the lattice analyses. These results agree well with the experimental values from Eq. (14) — the largest discrepancy is $1.7\sigma$ for the ratio of $\Delta M$ values — but they have much larger uncertainty.

IV.4 $V_{ts}$ and $V_{td}$

Because the experimental values for $\Delta M_{q}$ are so accurate, a better approach to understanding the implications of our improved lattice QCD results for the relevant hadronic matrix elements is to turn the analysis of the previous subsection on its head. That is, to use our results and the experimental values for $\Delta M_{q}$ to determine values for $|V_{ts}|$ and $|V_{td}|$ from Eq. (15) (taking a value for $V_{tb}$ from Eq. (16) CKM (2018)). $|V_{ts}|$ and $|V_{td}|$ obtained this way can then be compared to other determinations that make use of CKM unitarity as a test of that unitarity.

The ratio of $|V_{ts}|$ to $|V_{td}|$ can be obtained more accurately than the separate CKM elements because this makes use of the hadronic parameter $\xi$ (Eq. (13)) in which a lot of the lattice QCD uncertainties cancel (see Section IV.1).

Our results are

[TABLE]

Figure 9 plots the $\pm 1\sigma$ constraints on $|V_{td}|$, $|V_{ts}|$ and their ratio from our results as the dark grey lozenge. Results determined by other lattice QCD calculations Bazavov et al. (2016); Boyle et al. (2018) are also shown along with a recent determination using sum rules King et al. (2019). Also shown as light pink and orange lozenges are results from fits to the CKM unitarity triangle using results from many different processes CKM (2018); UTf (2018). Particularly relevant here is the green lozenge which results from a unitarity triangle fit that includes tree-level processes only CKM (2018), and therefore not $B_{s}/B_{d}$ oscillations. Tension between results derived from $\Delta M_{q}$ (as here) and the results derived from tree-level processes and unitarity would imply the existence of new physics in loop processes.

The Fermilab/MILC results (red lozenge in Figure 9) highlighted an approximately $2.0\sigma$ tension between their values for $V_{ts}$ and $V_{td}$ and those from unitarity fits. See Blanke and Buras (2016); Di Luzio et al. (2018) for examples of the possible implications of this.

Our results show no such tension. Our values for $V_{ts}$ and $V_{td}$ separately agree with the {CKMfitter, tree} results in Eq. (16) within $1\sigma$ and the difference in the ratio amounts to $1.8\sigma$. This limits the scope for new physics in loop-induced processes. However, our ratio for $|V_{td}|/|V_{ts}|$ joins the systematic trend of the previous results shown in Figure 9 in being below that of {CKMFitter, tree}.

IV.5 $B_{q}\rightarrow\mu^{+}\mu^{-}$ decay

The rare decays $B_{q}\rightarrow\mu^{+}\mu^{-}$ have very small branching fractions in the SM since they proceed through $W$ box diagrams and $Z$ penguins and are helicity-suppressed. New physics might then be seen if the experimental and SM branching fractions can be determined to be different to sufficient accuracy.

The hadronic parameter that enters the SM branching fraction is the $B_{q}$ meson decay constant Bobeth et al. (2014) but it appears along with the CKM elements $\left|V^{*}_{tq}V_{tb}\right|$. The uncertainty in the value of the appropriate CKM element is now the largest uncertainty in the value of the SM branching fraction Bazavov et al. (2018).

An alternative method for determining the branching fraction is to take a ratio to $\Delta M$ Buras (2003). In the SM (and extensions with minimal flavour violation) the CKM elements cancel out of this ratio. The decay constant also cancels and the hadronic parameter that remains in the ratio is the bag parameter.

The formula for the time-averaged branching fraction De Bruyn et al. (2012), as measured in the experiment, is then given in the SM by

[TABLE]

Here $C_{A}(\mu_{b})$ includes electroweak and QCD corrections and is given for $\mu_{b}$= 5 GeV in Bobeth et al. (2014). We use $C_{A}(\mu_{b})$ = 0.4694(36) Bazavov et al. (2016). The lifetime $\tau_{B_{q}^{H}}$ that appears in this formula is that of the heavy neutral eigenstate De Bruyn et al. (2012). For the $B_{d}$ this can be taken as the average lifetime, 1.520(4) ps Amhis et al. (2017); hfl (2018) but for the $B_{s}$ there is a measured difference of lifetimes and the heavy eigenstate has the longer lifetime, 1.615(9) ps Amhis et al. (2017); hfl (2018). Values for $S_{0}(x_{t})$ and $\eta_{2B}$ are given in Section IV.4.

Using our results for the bag parameters for $\mathcal{O}_{1}$ given in Eq. (12) and the experimental values for $\Delta M$ in Eq. (14) we obtain the following values for the branching fractions:

[TABLE]

We can also obtain the ratio of branching fractions for $B_{s}$ and $B_{d}$ Buras (2003)

[TABLE]

Here we have dropped the terms in $m^{2}_{\mu}/M^{2}_{B_{q}}$ since they are negligible. There is a lot of cancellation in this ratio, including of systematic errors in the ratio of bag parameters (see our results in Eq. (12)). We obtain the result

[TABLE]

Figure 10 shows our predictions in the SM for the branching fractions from Eqs. (20) and (22) as the grey lozenge. The red lozenge shows lattice QCD predictions Bazavov et al. (2018) for the branching fractions using the direct approach where the hadronic parameter needed is the decay constant and this is combined with input for the CKM elements $V_{tq}$ and $V_{tb}$, along with other factors. The errors in the results from Bazavov et al. (2018) are dominated by uncertainties in the CKM elements, which are taken from a global unitarity triangle fit that includes both tree and loop-induced processes333Note that we include the constraint on the ratio from the July 2019 update of Bazavov et al. (2018)..

Figure 10 shows good agreement between the two lattice QCD predictions. This reflects the fact that, as described in Section IV.4 our results for the bag parameters yield CKM elements $|V_{ts}|$ and $|V_{td}|$ in agreement with CKM unitarity determinations. Our results imply consistency of the CKM matrix (within uncertainties) and hence the two approaches of using the decay constants plus CKM elements or using the bag parameters and $\Delta M_{q}$ will agree.

Note that our results in Eq. (20) include uncertainties in the parameters of Eq. (21). They do not include uncertainties from electromagnetic corrections to the decay process. These are estimated to lead to a reduction of 0.3–1.1% in the muonic branching fractions in Beneke et al. (2018). This is not significant given the current uncertainties in our SM predictions, but will need to be addressed as reduce uncertainties in future.

The blue band in Figure 10 shows the current experimental situation. The decay $B_{d}\rightarrow\mu^{+}\mu^{-}$ has only been seen with $3\sigma$ significance Khachatryan et al. (2015). Recent LHCb Aaij et al. (2017) and ATLAS Aaboud et al. (2019) results give upper bounds to the branching fraction of $3.4\times 10^{-10}$ and $2.1\times 10^{-10}$ respectively. These bounds are outside the range of Figure 10. For the branching fraction for $B_{s}\rightarrow\mu^{+}\mu^{-}$, the Particle Data Group quotes an average value of $3.0(4)\times 10^{-9}$ using results from ATLAS Aaboud et al. (2019), CMS Chatrchyan et al. (2013) and LHCb Aaij et al. (2017). The $\pm 1\sigma$ variation gives the width of the band in Figure 10.

Although no significant tension between experiment and the SM predictions is visible in this figure, it does give encouragement that we are reaching a point where further reductions in uncertainties will start give serious SM constraints, at least for $B_{s}\rightarrow\mu^{+}\mu^{-}$.

IV.6 Contributions to $\Delta\Gamma$

Another physical observable from neutral $B$ meson systems is that of the decay width difference of the eigenstates, $\Delta\Gamma$. This has been measured for the $B_{s}$ at 13% of the average width, but only an upper limit exists for the $B_{d}$ Tanabashi et al. (2018). The prediction for the width differences in the SM is given in Beneke et al. (1999); Lenz and Nierste (2007) in terms of the matrix elements of several 4-quark operators. We give results here for the matrix elements of those labelled $R_{0}$, $R_{1}$ and $\tilde{R}_{1}$ in Lenz and Nierste (2007).

$R_{0}$ Lenz and Nierste (2007) is a combination of $O_{1}$, $O_{2}$ and $O_{3}$ which is a $1/m_{b}$-suppressed operator up to corrections of $\mathcal{O}(\alpha_{s}^{2})$ times a leading order operator.

[TABLE]

evaluating the radiative corrections at $\mu=m_{b}$. We can obtain the matrix elements for this operator through $\mathcal{O}(\alpha_{s})$ from our lattice calculation. $R_{1}$ and $\tilde{R}_{1}$ are proportional to $\mathcal{O}_{4}$ and $\mathcal{O}_{5}$ respectively and will be discussed further below.

The matrix elements for $R_{0}$ in Eq. (23) can be rewritten in terms of our bag parameters using the definition in Eq. (4):

[TABLE]

Writing it in this way makes clear (setting $B^{(n)}_{B_{q}}$ to 1 and $\alpha_{s}$ to zero) the expected cancellation at leading order to leave matrix elements that are $\mathcal{O}(1/m_{b})$. Our evaluation of the term in square brackets above yields

[TABLE]

Note that the uncertainties here include both those for missing $\alpha_{s}^{2}$ terms in the matching of lattice QCD operators to continuum operators and also the effect of missing $\alpha_{s}^{2}$ terms in the definition of $R_{0}$. This latter uncertainty is estimated by calculating the size of the $\alpha_{s}$ corrections in Eq. (23) and multiplying by $\alpha_{s}$. This gives a 35% uncertainty, which dominates the error quoted in Eq. (25). To assist with numerical evaluation we have replaced the square ratio of masses in Eq. (IV.6) with $3\eta^{q}_{3}$; values for this can be found in Table 1. Eq. (25) avoids use of a perhaps somewhat arbitray definition of a bag parameter for $R_{0}$ given in Lenz and Nierste (2007). The numerical factors show clearly that this is a $1/m_{b}$-suppressed operator by being of size $\Lambda_{\text{QCD}}/m_{b}\approx$ 10%.

$R_{1}$ and $\tilde{R}_{1}$ are defined as Lenz and Nierste (2007)

[TABLE]

The matrix elements for $B_{s}$ and $B_{d}$ can then be determined from our results in Table 5. Our bag parameters for $O_{4}$ and $O_{5}$ are given in Table 6. In Lenz and Nierste (2007) bag parameters for $R_{1}$ and $\tilde{R}_{1}$ are defined in such a way as to set the squared mass ratios in Eq. (4) to 1. This means that the bag parameters for $R_{1}$ and $\tilde{R}_{1}$ for the definition in Lenz and Nierste (2007) can be recovered from our bag parameters by multiplying by $3\eta_{4}^{q}/7$ and $3\eta_{5}^{q}/5$ respectively. These factors are larger than 1. Note however that the impact of both $R_{1}$ and $\tilde{R}_{1}$ on $\Delta\Gamma$ is tiny because of the $m_{q}/m_{b}$ factors in their definition.

Matrix elements of the numerically more important $R_{2}$ and $\tilde{R}_{2}$ operators Lenz and Nierste (2007), along with those of $R_{3}$ and $\tilde{R}_{3}$ cannot be directly obtained from our current results because they contain derivatives on the light quark fields inside the 4-quark operator. Results of the calculations of these matrix elements are discussed separately in Davies et al. (2019).

V Conclusions

We give results from the first ‘second-generation’ lattice QCD calculation of the matrix elements that contribute to $B_{s}$ and $B_{d}$ mixing in and beyond the Standard Model. We include $c$ quarks in the sea for the first time and have a range of $u/d$ quark masses (taken to be equal) that go down to the physical value. We use radiatively-improved NRQCD for the $b$ quark action. By calculating the ratio of the matrix elements of the 4-quark operators to the square of the decay constant times mass (proportional to a quantity known as the bag parameter) we obtain results with very little dependence on the $u/d$ quark mass or the lattice spacing. This gives us more accurate results than previous calculations for these ratios, and the associated bag parameters, for all five $\Delta B=2$ operators.

Our key results are given in Tables 5 and 6. Table 5 gives the ratio of matrix elements to $(fM)^{2}$ with our final physical values given in the last row. These are the numbers that should be used to reconstruct the 4-quark operator matrix elements by multiplying by $(f_{B_{q}}M_{B_{q}})^{2}$. Table 6 converts these ratios into bag parameters, defined in Eq. (4). These numbers can be compared to unity, the result expected in the vacuum saturation approximation. Our error budget for the bag parameters is given in Table 7. We have uncertainties of 4–7% for the individual bag parameters. This uncertainty is dominated by missing higher orders in the perturbative matching to the continuum 4-quark operators. The uncertainty is reduced to around 2% in the ratio of bag parameters for $B_{s}$ to $B_{d}$, since this renormalisation cancels. The correlations between results for different operators are given in Table 10 of Appendix C. Our $B_{s}$ to $B_{d}$ ratios are now accurate enough to see that they are above 1 for $O_{1}$ and $O_{3}$ and below 1 for $)_{4}$ and $O_{5}$.

Our results for the key $O_{1}$ bag parameters that appear in SM phenomenology are (repeating Eq. (12))

[TABLE]

where we give the RGI bag parameter as defined in Eq. (11). Multiplying by decay constant values obtained using HPQCD’s approach to $b$-physics with HISQ quarks McNeile et al. (2012a) by the Fermilab/MILC collaboration Bazavov et al. (2018), we also obtain (repeating Eq. (13))

[TABLE]

In Section IV we discuss the phenomenology from our results. We obtain values for $\Delta M$ for $B_{s}$ and $B_{d}$ in Eq. (17) to be compared to experiment.

Alternatively, and more usefully, we can combine our results with experiment to obtain the CKM elements $|V_{ts}|$ and $|V_{td}|$ and their ratio. These values are given in Eq. (18) and Figure 9 shows the constraints they give in the $V_{td}$-$V_{ts}$ plane. Our results are the most accurate determinations of these CKM elements using lattice QCD and show good consistency with determinations from tree-level processes assuming CKM unitarity. This means that we see no signs of new physics in neutral $B$-meson oscillations at this improved level of accuracy.

We derive results, by taking a ratio to $\Delta M$, for the branching fractions for $B_{s}$ and $B_{d}$ to decay to $\mu^{+}\mu^{-}$, a key mode for new physics searches at LHC. Our results are given in Eq. (20) and (22). Figure 10 shows a comparison of our predictions to those from the recent Fermilab/MILC calculation of decay constants Bazavov et al. (2018) along with the current experimental picture. This is encouraging for future tests of new physics contributions to these rare decay processes.

Finally, we give values in Eq. 25 and below Eq. (26) for the matrix elements for the $R_{0}$, $R_{1}$ and $\tilde{R}_{1}$ operators that contribute to the SM prediction for the width difference $\Delta\Gamma$.

To improve accuracy further in future requires improving the matching of the lattice QCD 4-quark operators to those in the continuum; this is the dominant source of uncertainty in our bag parameters. Since lattice QCD perturbation theory is so hard a renormalisation method that can be implemented within the lattice calculation and then matched perturbatively to $\overline{\text{MS}}$ in the continuum could be preferable. A symmetric momentum-subtraction scheme (RI-SMOM) has been found to work well for kaon mixing calculations Garron et al. (2016) but attention must be paid to removing nonperturbative artefacts in these schemes if high accuracy is to be achieved Lytle et al. (2018). Such a method would need to be implemented with a relativistic quark action on lattices with fine enough lattice spacing to allow quark masses close to that of the $b$ for $am\,^{<}_{\sim}1$. This has been a very successful strategy for HISQ quarks for $B$-meson decay constants McNeile et al. (2012a, b); Bazavov et al. (2018), but in that case the decay constants calculated with HISQ do not need any renormalisation. Calculating 4-quark operator matrix elements is much harder, but still feasible. The ETM work with twisted mass quarks Carrasco et al. (2014) is encouraging for this programme (although they used $n_{f}=2$ gluon fields and RI-MOM renormalisation) as is the work by RBC/UKQCD using domain-wall quarks Boyle et al. (2018) (although they have only calculated $B_{s}$ to $B_{d}$ ratios so far). It seems clear that in the next few years improvements in this direction will be possible, pushing uncertainties on bag parameters down to the $\sim$2% level. This will allow $|V_{ts}|$ and $|V_{td}|$ to be determined to 1%.

Acknowledgements We are grateful to the MILC collaboration for the use of their gauge configurations and code. The results described here were obtained using the Darwin Supercomputer of the University of Cambridge High Performance Computing Service as part of STFC’s DiRAC facility. This work was funded by STFC, the Royal Society, the Wolfson Foundation and the US DOE and National Science Foundation. We thank Chris Bouchard, Elvira Gámiz and Alex Lenz for useful discussions. One of the authors (GPL) is grateful to the Department of Applied Mathematics and Theoretical Physics, Cambridge University, for their hospitality during the two visits when much of this analysis was done.

Appendix A Fitting Protocols

We extract mixing amplitudes and decay constants by fitting Monte Carlo data for 2-point and 3-point matrix correlators for each meson:

[TABLE]

where $\mathbf{s}$ is a 3-vector of meson sources, the sums over spatial $\mathbf{x}$ and $\mathbf{y}$ project onto zero 3-momentum, the times satisfy $0<t<T$, and $\beta=1,2\ldots 5$ labels the mixing operator (Figure 11). The sources include a local source, corresponding to

[TABLE]

and two smeared sources: see Dowdall et al. (2013a) for details. We also examine the vector of correlators

[TABLE]

where

[TABLE]

is the leading NRQCD correction to $J_{A_{0}}^{(0)}$. We use the corrected current444Note that the NRQCD-HISQ temporal axial current that we use here is correct through the same order in $\alpha_{s}$ and $\Lambda_{\text{QCD}}/m_{b}$ as that of our 4-quark operators. In Dowdall et al. (2013a) we used a more highly-corrected temporal axial current to determine $f_{B}$. to evaluate the decay constant (in lattice units):

[TABLE]

where coefficients $z_{A_{0}}$ depend on the lattice spacing and were calculated in Dowdall et al. (2013a) using Monahan et al. (2013) (see Table 2 for the values we use here). The values for $\alpha_{s}$ (from Colquhoun et al. (2015)) are given in Table 4.

The 2-point and 3-point correlators are calculated in the standard way. For 3-point correlators this involves combining propagators from a local source (at $t$) into ‘open-meson’ propagators Bazavov et al. (2016) which are then closed off at [math] and $T$ (which cover all time-slices away from $t$) with local or smeared meson and anti-meson operators. We average correlators over 16 values of $t$ for improved statistics. The smearing functions we use are given in Colquhoun et al. (2015).

Fits for the $B_{s}$ proceed in two steps. First the 2-point correlators $\mathbf{G}(t)$ and $\mathbf{G}^{(1)}(t)$ are fit simultaneously. Then the best-fit amplitudes and energies from that fit are used as priors for a simultaneous fit of all the 3-point correlators $\mathbf{G}_{n}(t)$. Having finished the $B_{s}$ fits, we follow the same approach with the $B_{d}$ correlators but constraining (via the priors) the $B_{d}$’s fit parameters to be within 20% of the corresponding values for the $B_{s}$. We discuss all of these fits in what follows.

A.1 Fitting Two-Point Correlators

We fit the 2-point correlators to a formula of the form (in lattice units)

[TABLE]

where the fit parameters $\mathbf{p}$ are comprised of all $\mathbf{c}_{n}$, $\mathbf{c}_{n}^{o}$, $E_{n}$, and $E_{n}^{o}$. Here $\mathbf{c}_{n}$ and $\mathbf{c}_{n}^{o}$ are 3-component vectors, and $j_{n}$ and $j_{n}^{o}$ scalars, where:

[TABLE]

In the exponents, $E_{n}$ and $E_{n}^{o}$ are the energies of the lowest-lying states with zero 3-momentum that couple to the sources. The second (oscillating in time) term in each correlator is due to taste-doubling caused by the staggered-quark HISQ action for the light quarks (see Dowdall et al. (2012b)). $M_{n}$ and $M_{n}^{o}$ are the physical masses corresponding to states $|E_{n}\rangle$ and $|E_{n}^{o}\rangle$, respectively. We keep $N=6$ terms, but fit results are the same for any $N\geq 5$.

We use a Bayesian fit procedure Lepage et al. (2002). Fits for the $B_{s}$ mesons on our coarsest lattices (0.15 fm) use the following Bayesian priors for the energies,

[TABLE]

where $\Delta E_{n}\equiv E_{n}-E_{n-1}$, and the logarithms indicate log-normal priors for energies. Energies are rescaled in proportion to the lattice spacing to make priors for the other lattices. The priors for local and smeared amplitudes are

[TABLE]

on the coarsest lattices. Amplitudes for the local source are rescaled by $a^{3/2}$ for the other lattices. The smeared sources are designed to be lattice-spacing independent and so we use the same prior for the other lattices. Finally the priors for parameters $j_{n}$ and $j_{n}^{o}$ are

[TABLE]

on the coarse lattice, and again scale like $a^{3/2}$ for other lattices.

These central values for these priors were based upon fit values for the ground state $B_{s}$. The uncertainties assigned to the priors are large: for example, the priors are typically 2000–5000 times broader than the final fit errors for the ground state parameters that we need for our analysis. Replacing the central values by random values drawn from the prior distributions leaves our results unchanged within errors.

We need to apply SVD cuts to the data’s correlation matrix because of the large number of correlators being fit. The procedure for determining the SVD cuts is described in Appendix D; typically the cuts affect less than half of the data modes. Fits for $\mathbf{G}(t)$ were constrained to $t$ values between the values for $t_{\mathrm{min}}$ and $t_{\mathrm{max}}$ shown in Table 8; data from larger $t$’s was too noisy to be useful. To keep the number of fit points down (see Appendix D), we restricted the fits for $G^{(1)}(t)$ to the range of $t$s between $(t_{\mathrm{min}}+t_{\mathrm{max}})/2$ and $t_{\mathrm{max}}$.

A.2 Fitting Three-Point Correlators

The fit function for the 3-point correlators is substantially more complicated:

[TABLE]

where the fit parameters $\mathbf{p}$ include all of the 2-point correlator parameters plus the $V_{nm}$, $V_{nm}^{o}$, and $V_{nm}^{oo}$.

We fit the 3-point amplitudes over the range $t_{\mathrm{min}}\leq t\leq T-t_{\mathrm{min}}$ for the values of $T$ shown in Table 8. Parameters $\mathbf{c}_{n}$, $\mathbf{c}_{n}^{o}$, $E_{n}$, and $E_{n}^{o}$ are the same as in the 2-point correlators; we use the results from the fits to the 2-point correlators as priors for these parameters in our 3-point fits. The priors on the coarsest lattices for each of the mixing amplitudes $V(\beta)$, $V^{o}(\beta)$ and $V^{oo}(\beta)$ are:

[TABLE]

these are scaled in proportion to $a^{3}$ for the other lattices. Note that $V_{n,m}(\beta)$ and $V^{oo}_{n,m}(\beta)$ are symmetric under interchange of $n$ and $m$. We are interested in the ground-state value for

[TABLE]

We introduce two simplifications to the analysis that make our fits run 20–100 times faster, without affecting fit results or precision. The first simplification is to replace both the data and the fit function in the fits with their sums over $t$,

[TABLE]

while keeping the same fit parameters Bouchard (private communication). This reduces the number of data points to be fit for our 0.09 fm lattice, for example, from 1050 to 180. Note that the Monte Carlo data for $\mathbf{G}_{n}(t,T)$ do not vary much with $t$, as expected from Eq. (42).

The second simplification is to marginalize all fit parameters other than those associated with the ground state Hornbostel et al. (2012). We do this by splitting the fit function (Eq. (42)) into two parts, one that involves only the ground state (i.e., either the $B_{s}$ or $B_{d}$) and the other with the remaining terms:

[TABLE]

We then replace the fit data $\mathbf{G}^{\mathrm{lat}}_{\beta}$ by

[TABLE]

where the prior values for the fit parameters are used in $\Delta\mathbf{G}_{\beta}^{\mathrm{fit}}$. At the same time, we replace the fit function by just its ground-state term:

[TABLE]

This reduces the number of fit parameters from 450 to 9 (for $N=6$). Marginalization works particularly well here because we have excellent priors for the amplitudes and energies, from the 2-point correlators, and because the mixing parameters enter the fit function linearly. Also the marginalized fit function is $t$-independent, making the first simplification (summing over $t$) quite natural.

Again we need SVD cuts, but the need is greatly reduced by summing over $t$. We used the method outlined in Appendix D; typically the cuts modified around 70% of the data modes.

We tabulate simulation results (using Eqs. (34), (37), (38) and (44)) for the dimensionless ratio

[TABLE]

with $B_{q}=B_{s},B_{d}$ in Table 8. This table also shows sample values of $\chi^{2}$ from the various (2-point and 3-point) correlator fits when we include random SVD and prior noise, as discussed in Appendix D.4. Without noise, the $\chi^{2}$s per degree of freedom are much smaller than 1.0, as expected. Also following Appendix D.4, we tested the uncertainties from our fits using simulated data. Fits to simulated data reproduced the input parameters to within errors.

We verified that marginalization and averaging over $t$ have negligible effect on our results. With our smallest lattice spacing (0.09 fm), for example, undoing both optimizations gives the following values:

[TABLE]

These agree well with the values in Table 8 (for set 7), but took far longer to compute.

In Table 9, we show sample error budgets for these quantities from simulations on the 0.09 fm lattice; others are similar. The dominant source of uncertainty is from the Monte Carlo statistics.

A.3 Blind Analysis

This analysis was blinded by multiplying the 3-point correlators by a random normalization factor. The random factor was removed only after the entire analysis was completed and this paper written.

Appendix B Chiral Fit

Although we have results at physical pion masses we do not rely on these simply for our final value. We include also results at heavier-than-physical pion masses, which are statistically more precise, by using a fit to the dependence on the pion mass based on chiral perturbation theory. This gives the coefficients of the non-analytic ‘chiral logarithms’ in $m_{\pi}^{2}\log(m_{\pi}^{2})$; in addition we include analytic terms to allow both for staggered quark discretisation effects, the unphysically heavy $u/d$ quark masses in the sea and for mistuning of valence quark masses. Performing a fit to results at multiple pion masses then tests the dependence expected from chiral perturbation theory.

Our principal results consist of values for the “reduced” matrix elements of the 4-quark operators,

[TABLE]

on each configuration set (Table 5). We fit the $R^{n}_{q}$ to the form

[TABLE]

where we suppress the indices $n$ and $q$ on each term and each parameter, $p$, for clarity. The functions used in each term are discussed below. $\chi_{\text{log}}$ is given in Eq. (B.1); $J$ in Eq. (B.2); $\delta X$ in Eq. (59); $\delta x_{l}$ in Eq. (61) and $\delta_{v}$ in Eq. (62). Note that this fit is done after applying the additional uncertainties discussed in Section III.2 to allow for matching and discretisation effects.

To derive this form, we make use of the results in Bernard (2013), where Appendix A gives the dependence on light meson masses of the bag parameters at one-loop in heavy meson staggered chiral perturbation theory. This builds on the continuum heavy meson chiral perturbation theory results of Detmold and Lin (2007). There is a lot of cancellation of chiral logarithms between 4-quark operator matrix elements and decay constants so that, as we discuss below, the remaining chiral logarithm terms ($\chi_{\text{log}}$ and $J$ in Eq. B) in the bag parameters (and equivalently in $R$) have small coefficients. This expected very benign dependence on the light quark mass is another reason for working with the bag parameters as we do here, rather than the 4-quark operator matrix elements.

The chiral perturbation theory for the bag parameters is given in Bernard (2013) in the form (using our notation)

[TABLE]

Here $\beta_{n}$ is the low-energy constant (value at zero pion mass) for $R^{n}$ and $\beta_{n}^{\prime}$ is the equivalent term for 4-quark matrix elements between vector heavy-light mesons. For $O_{1}$, $\beta_{1}^{\prime}=\beta_{1}$. $S$ comes from ‘tadpole’ diagrams (with $+$ for $n$ = 1, 2 and 3 and $-$ for $n$ = 4, 5) and $Q$ from ‘sunset’ diagrams that connect pseudoscalar and vector mesons. $\tilde{T}$ and $\tilde{Q}$ are ‘wrong-spin’ tadpole and sunset terms respectively. Below we discuss the content of these functions in terms of the important non-analytic chiral logarithms and the effect on these of the discretisation effects in the staggered quark formalism. This enables us to transcribe Eq. (53) into the simpler Eq. (B) that we will use. We now discuss each of these terms in turn.

B.1 Tadpole diagrams

The results in Bernard (2013) are given in terms of meson masses that include those for pions of different taste that appear in the staggered quark formalism. Thus a term that would be a simple chiral logarithm in the continuum can appear in a number of guises, one of which is as an average over the masses of all tastes of pion. On fine enough lattices this will become a continuum logarithm plus discretisation effects. In fact, for the fully unquenched case that we study here ($m^{\text{val}}=m^{\text{sea}}$), staggered chiral perturbation theory typically arranges itself to cancel taste effects inside chiral logarithms so that non-analyticities in $a^{2}$ cancel as $m_{u/d}\rightarrow 0$ (see Appendix A of Colquhoun et al. (2016)). This also happens here. The chiral logarithm in $R_{d}$ from tadpole diagrams ($S$ in Eq. (53)) appears in the form

[TABLE]

Here $I$ denotes the singlet (largest mass) pion taste. We can compare this function to the corresponding continuum chiral logarithm

[TABLE]

where $P$ denotes the lightest (Goldstone) pion taste. This comparison is shown in Figure 12 using taste-splittings for HISQ pions for the lattice spacing values that we use in this calculation. We give results for our range of $m^{2}_{\pi}\equiv m^{2}_{\pi,P}$ values from the physical point, 0.018 $\mathrm{GeV}^{2}$, to 0.09 $\mathrm{GeV}^{2}$. $\Lambda_{\chi}=4\pi f_{\pi}$ = 1.64 GeV and we take $\mu_{\chi}=1.0$ GeV. The difference between the HISQ and continuum chiral logarithm terms is sufficiently small that we simply allow for that discrepancy in our treatment of discretisation effects. This is included through the $\delta X$ terms in Eq (B) discussed below.

We therefore take the chiral logarithm terms in Eq. (B) with continuum form:

[TABLE]

Here $m^{2}_{\eta}=(2m_{\eta_{s}}^{2}+m_{\pi}^{2})/3$ and we use masses of Goldstone taste $\pi$ and $\eta_{s}$ in this expression. $(\mathrm{phys})$ denotes the value of the previous expression evaluated for physical masses so that the total right-hand side vanishes at that point. $m^{2}_{\eta}$ changes very little as $m_{\pi}$ changes so these terms in the fit do very little. The parameters $p^{q}_{\text{log},n}$ are given priors $+1.0(3)$ for $n=1,2,3$ and $-1.0(3)$ for $n=4,5$ as this logarithm appears with opposite sign for $\mathcal{O}_{4,5}$. The prior width allows for modification of the coefficients from missing higher order terms in chiral perturbation theory.

B.2 Sunset diagrams

We now turn to the term denoted $J$ in Eq. (B) that comes from the sunset diagram term $Q$ in Eq. (53). This would also take the form of a chiral logarithm $m_{\pi}^{2}\log(m^{2}_{\pi}/\mu_{\chi^{2}}^{2})$ in the infinite heavy meson mass limit in the continuum. For finite heavy meson mass, however, $J$ is modified by terms that depend on heavy meson mass differences, because there is a pseudoscalar to vector heavy meson transition inside the sunset diagram. The form of $J$ as a function of pion mass and heavy meson mass difference, $\Delta$, is given in Eq. 6.17 of Bazavov et al. (2012b), which considers chiral perturbation theory terms for the heavy-light meson decay constant. In that case the appropriate value for $\Delta$ can include heavy-strange to heavy-light mass differences as well as vector to pseudoscalar mass differences. Here, when we consider $R_{d}$ and $R_{s}$, that does not happen and we only have to consider the case where $\Delta=M_{B_{(s)}^{*}}-M_{B_{(s)}}$. Then $\Delta$ takes the value 45 MeV for $B^{*}-B$ Tanabashi et al. (2018) and we take the same value for $B^{*}_{s}-B_{s}$ since any differences are expected Dowdall et al. (2012b) and seen to be Tanabashi et al. (2018) small. Figure 13 compares the function $J(m_{\pi},\Delta)/\Lambda^{2}_{\chi}$ with that of the chiral logarithm to which it is equal when $\Delta=0$. Even though $\Delta$ is small, and much smaller than $m_{\pi}$ through most of the range in which we work, we see $\Delta$ does have an impact, so that $J$ has smaller magnitude and gradient in $m^{2}_{\pi}$ than its associated chiral logarithm.

$J$ in Eq. (B) then takes the form

[TABLE]

where $J(m,\Delta)$ is given in Bazavov et al. (2012b). $J$ is multiplied by 3$g^{2}$, where $g$ is the $BB^{*}\pi$ coupling. We take the value of $g$ as 0.5, based on recent lattice QCD calculations Detmold et al. (2012); Bernardoni et al. (2015); Flynn et al. (2016). The uncertainty on $g$, both from the lattice calculations but also from the effect of missing higher order terms in chiral perturbation theory, is absorbed into the coefficient $p_{J}$. $p_{J}$ is the ratio of low-energy constants associated with the bag parameters for vector heavy-light mesons to that for pseudoscalars. For $\mathcal{O}_{1}$ we know that this ratio is 1 Grinstein et al. (1992). For the other operators, $n$=2–5, we do not. We therefore take a prior value and width on $p^{(d,s)}_{J,n}$ of 0(1), allowing either sign. For $p^{(d,s)}_{J,1}$ we take 1.0(3) to allow for uncertainty in $g^{2}$. Note that in the case where $\Delta=0$ the coefficient of the chiral logarithm $(m^{2}_{\pi}/\Lambda^{2}_{\chi}\log(m^{2}_{\pi}/\mu^{2}_{\chi})$ in the chiral perturbation theory for $O_{1}$ is $(3g^{2}-1)/2$ Becirevic et al. (2007), which is small (-0.125) when $g=0.5$.

B.3 Wrong-spin and other effects

The heavy meson staggered chiral perturbation theory analysis of Bernard (2013) showed that $O_{1,2,3}$ and $O_{4,5}$ can mix through ‘wrong-spin’ staggered taste-effects ($\tilde{T}$ and $\tilde{Q}$ in Eq. (53)). The size of these terms depends on the size of the light meson taste-splittings for the staggered action. For the asqtad action used for light quarks in Bazavov et al. (2016) they were of some concern and it was important to include these effects explicitly in a full fit to all 5 operators. For the HISQ action that we use here these effects are much smaller. The question then becomes whether they are distinct in magnitude or form from discretisation effects from other sources that are already included in our analysis.

The wrong-spin contributions from tadpole diagrams for $B_{d}$ involve differences of chiral logarithms for different taste pions (hence cancelling in the absence of taste-splitting effects) along with hairpin terms that have coefficients $a^{2}\delta_{V}^{\prime}$ and $a^{2}\delta_{A}^{\prime}$ that are themselves the size of the unit of taste-splitting Colquhoun et al. (2016). As an illustration of the impact of these terms we examine the terms that are differences of chiral logarithms. These appear in 3 ‘signatures’ : $V-A$, $P-I$ and $-P+2T-I$. Here the letter denotes the pion taste, ordered in increasing mass as: $P$, $A$, $T$, $V$, $I$. Figure 14 illustrates the size and behaviour of these terms for the $-P+2T-I$ example that mixes $O_{2}$ and $O_{3}$, the simplest because there are no additional hairpin corrections. The function plotted is

[TABLE]

We see that $\tilde{X}$ falls rapidly with lattice spacing (approximately as $(a\Lambda)^{4}$ and has a slope with $m^{2}_{\pi}$ that also falls with lattice spacing (approximately as $(a\Lambda)^{2}$). This behaviour is generic for terms that arise from taste-splittings in this way.

The wrong-spin tadpole terms have a variety of coefficients multiplying them that correspond to ratios of 4-quark operator matrix elements (within the groupings ${1,2,3}$ and ${4,5}$). Most of the coefficients for the wrong-spin tadpole terms appearing in the chiral expansion for $O_{y}$ are of the form $\beta_{x}/(4\beta_{y})$ where $\beta_{x}$ is the low energy constant associated with operator $O_{x}$. The exception is $O_{1}$, where the coefficient is $2(\beta_{2}+\beta_{3})/\beta_{1}$. If all the 4-quark operator matrix elements were of the same size then the coefficients would be $\mathcal{O}(1/4)$. $O_{3}$ and $O_{5}$ have smaller matrix elements than the others, however, if we consider the vacuum saturation approximation. Hence $\beta_{x}/\beta_{y}$ can be of size 2 for $O_{5}$ and 5 for $O_{3}$.

We are already including $a^{2}$ and $a^{4}$ errors in Eq. (10), but contributions from wrong-sign tadpole terms differ between $B_{d}$ and $B_{s}$ mesons. The largest contributions are for the $B_{d}$ meson; we allow for them and similar terms that arise from sunset diagrams (and so contain differences of $J(m_{\pi,t},\Delta)$) along with residual right-sign taste-effects by including terms

[TABLE]

in the chiral fit, Eq. (B). We include them with the chiral fit rather than with other discretisation effects in Eq. (10) because they arise from the staggered quark action and hence carry no $am_{b}$ dependence. Here

[TABLE]

allows for $m_{\pi}^{2}$ dependence; it varies between $-0.5$ and $0.5$ over our range of parameters. Similar terms are not needed for $B_{s}$ mesons because the effects are smaller and can be simulated by Eq. (10). We take the priors $p_{a^{4}}$ and $p_{m_{\pi}^{2}a^{2}}$ to be 0(2) since, although this is an unnecessarily broad prior for some $n$, it allows a reasonable size for all the possibilities.

B.4 Analytic terms

The three terms given symbol $\delta$ on the last line of Eq. (B) are simple polynomials to account for mistuning of sea and valence quark masses from their physical values. We use

[TABLE]

where $m_{l}$ and $m_{s}$ are the sea $u/d$ and $s$ quark masses from Table 3. The physical value for the $m_{l}/m_{s}$ ratio we take as 27.18(10) from Bazavov et al. (2018). The factor of 1/10 converts $m_{l}/m_{s}$ to the size of the parameters that appear in chiral perturbation theory as a ratio of meson masses to $\Lambda_{\chi}=4\pi f_{\pi}$. Eq. (B) includes terms in $\delta x_{l}$ and $(\delta x_{l})^{2}$. We take a prior of size 0(1) for each coefficient $p_{l}$ and $p_{l2}$ (for each operator and each $q$). We do not allow for mistuning effects for $c$ quarks in the sea since we expect thse effects to be negligible compared to those from light sea quarks. $\delta_{v}$ accounts for the mistuning of light and strange valence masses appropriate to $B_{d}$ or $B_{s}$. We take

[TABLE]

Again the coefficients for this term, $p^{q}_{v,n}$, have prior 0(1) for each $n$ and $q$. We do not include terms for $b$ quark mass mistuning since the tuning is very accurate and we expect any small mistuning to have negligible effect on bag parameters.

B.5 Finite-volume, Strong Isospin-breaking and QED Effects

Finite-volume effects can be estimated based on chiral perturbation theory. The results in Arndt and Lin (2004) show finite-volume effects of $\mathcal{O}(1\%)$ for the bag parameters of $O_{1}$ in small lattice volumes of size $L=$ 2.5 fm for physical $u/d$ quark masses. On our much larger lattices, with a minimum size of $L=$ 4.6 fm at physical $u/d$ masses for set 3, finite-volume effects will be a lot smaller. We conclude that this is a negligible effect at our current level of uncertainties.

Strong-isospin breaking and electromagnetic effects can also be estimated to be negligible at present. Our bag parameters show very little sensitivity to the $u/d$ quark mass and our ratios for $B_{s}$ to $B_{d}$ differ from 1 by at most 10% (Table 6). This suggests that changing $m_{l}$ to $m_{d}$ should only have a $\mathcal{O}(0.1\%)$ effect. Effects from the fact that the valence quarks have electromagnetic charge are estimated at below 0.1% for the decay constants in Bazavov et al. (2018). They come largely from QED effects on the tuning of quark masses. Since bag parameters are less sensitive to both heavy and light quark masses than decay constants, we conclude that QED effects on the bag parameters will be less than 0.1% and we neglect them. Note that QED effects can still enter $\Delta M_{q}$ or $\text{Br}(B_{q}\rightarrow\mu^{+}\mu^{-})$ through corrections to these processes from adding photons; these effects need to be considered separately.

Appendix C Correlations in Final Results

In this Appendix we describe the correlations between the uncertainties in different final results from our analysis. Our principal results consist of values for the reduced matrix elements of the 4-quark operators, $R^{n}_{q}$, defined in Eq. (51), evaluated for physical quark masses (Table 5). The results for a given meson and different operators are only weakly correlated, as shown in Table 10 for the $B_{s}$ meson. There is more correlation, but still small, between values of the bag parameters (Eq. (2)), because of the normalization factors $\eta_{n}^{s}$ (Eq. (4)). The means and standard deviations for these quantities are collected in Table 11.

Values of $R^{n}_{s}$ are highly correlated with values of $R^{n}_{d}$, for the same $n$, which is why the ratios ${R_{s}^{n}}/{R_{d}^{n}}$ have much smaller uncertainties. Uncertainties in these ratios are almost uncorrelated, however, with those in the $R^{n}_{q}$ (correlations are 0.06 or smaller). Thus correlations for the $B_{d}$ matrix elements $R^{n}_{d}$ can be easily constructed from the results in Table 11 and Table 10 for $R^{n}_{s}$ and $R^{n}_{s}/R^{n}_{d}$. The ratio of bag parameters $B_{B_{s}}^{(n)}/B_{B_{d}}^{(n)}$ is almost equal to $R^{n}_{s}/R^{n}_{d}$ (the difference being only from small quark mass effects) and has the same correlation matrix.

Appendix D SVD Cuts

D.1 The Problem

There are three inputs for our least-square fits to sets of correlators: 1) a collection of $N_{s}$ random (Monte Carlo) samples $\mathbf{G}^{(s)}$, where each sample is packaged as an $N_{G}$-dimensional vector; 2) a (vector) fitting function $\mathbf{G}(\mathbf{p})$ of fit parameters $\mathbf{p}$; and 3) a priori estimates (priors) for the fit parameters.

In the fit, the sample average

[TABLE]

is assumed to be a random sample drawn from a Gaussian distribution with mean $\mathbf{G}(\mathbf{p}^{*})$ for some set $\mathbf{p}^{*}$ of fit parameters, and a covariance matrix given approximately by

[TABLE]

Here $\mathbf{D}$ is the diagonal matrix of standard deviations, $\mathbf{D}_{ij}=\delta_{ij}\,\sigma_{\mathbf{G}_{i}}$, and $\mathbf{M}_{\mathrm{corr}}$ is the correlation matrix. The best-fit parameters are obtained by minimizing

[TABLE]

as a function of the parameters $\mathbf{p}$, where $\lambda_{n}$ and $\mathbf{v}_{n}$ are the eigenvalues and eigenvectors of the correlation matrix:

[TABLE]

Note that

[TABLE]

$\chi^{2}_{\mathrm{prior}}(\mathbf{p})$ is the part of $\chi^{2}(\mathbf{p})$ associated with the Bayesian priors used in the fit.

The approximation for the covariance matrix, Eq. (64), causes problems if the number of samples $N_{s}$ is insufficiently large compared with the number of data points $N_{G}$ Michael (1994); Michael and McKerrell (1995). In particular, the smaller eigenvalues of the correlation matrix are underestimated. Indeed it is obvious from Eq. (64) that there must be $N_{G}-N_{s}$ modes with zero eigenvalue when $N_{s}<N_{G}$. Underestimating eigenvalues exaggerates their importance in $\chi^{2}(\mathbf{p})$ (Eq. (66)), compromising the fit; and $\chi^{2}(\mathbf{p})$ is undefined if there are zero eigenvalues.

The underestimation of small eigenvalues is illustrated in Figure 15, which shows the ratio of $\lambda_{n}^{\mathrm{approx}}/\lambda_{n}^{\mathrm{exact}}$ for approximate eigenvalues estimated from random samples of different sizes drawn from a simulation of a known distribution (so we know the exact eigenvalues). The small eigenvalues are dragged down to zero by the need for zero modes when $N_{s}<N_{G}$. They then increase slowly as new samples are added, until $\lambda_{n}^{\mathrm{approx}}/\lambda_{n}^{\mathrm{exact}}\approx 1$ for all $n$ when $N_{s}\gg N_{G}$. Note that good approximations for all eigenvalues require $N_{s}$ to be 10–100 times larger than $N_{G}$. The figure shows results for $N_{G}=512$ pieces of correlated data; curves for $N_{G}=50$ (or $5000$) would be similar, but with more (or less) noise. The range of values covered by the eigenvalues also has little effect on the overall picture.

D.2 Choosing an SVD Cut

The problematic eigenvalues are those for which

[TABLE]

since, on average, individual terms in $\chi^{2}(\mathbf{p})$ should contribute approximately $1\pm\sqrt{2/N_{G}}$ to the total (based on the width of the $\chi^{2}$ distribution). Following Michael (1994); Michael and McKerrell (1995), we deal with these eigenvalues by introducing a cutoff $\kappa$ such that eigenvalues smaller than $\kappa\lambda_{\mathrm{max}}$ are replaced with $\kappa\lambda_{\mathrm{max}}$, where $\lambda_{\mathrm{max}}$ is the largest eigenvalue:

[TABLE]

Tuning $\kappa$ appropriately, this replacement increases the underestimated eigenvalues to a value that is at least as large as the exact eigenvalue (and probably a lot larger). Unlike in Michael (1994); Michael and McKerrell (1995), we do not renormalize the eigenvalues to preserve the trace of the modified matrix (see below).

We need curves such as those in Figure 15 to set $\kappa$, but we don’t know the exact eigenvalues in real applications. An approximate curve can be generated by comparing the eigenvalues of bootstrapped copies of our simulation results \big{\{}\mathbf{G}^{(s)}\big{\}} with the eigenvalues from the full ensemble. Each bootstrapped copy has $N_{s}$ samples, like the original ensemble. In this analysis, bootstrapped eigenvalues play the role of the approximate eigenvalues above, while eigenvalues computed directly from ensemble \big{\{}\mathbf{G}^{(s)}\big{\}} now play the role of the exact eigenvalues (since they specify the underlying distribution for the bootstrapped copies).

Ratios of these eigenvalues are plotted in Figure 16 (blue points) for examples with $N_{G}=50$ (top panel) and $N_{G}=500$ (bottom panel) data points, and $N_{s}=4N_{G}$ random samples for each data point. The error bars show the spread of values across the different bootstrapped copies. These points give us an approximate curve for $\lambda_{n}/\lambda_{n}^{\mathrm{exact}}$, from which we can determine an SVD cutoff.

The ensembles used in these examples were generated from a known distribution, so in this case we know the correct curve for $\lambda_{n}^{\mathrm{approx}}/\lambda_{n}^{\mathrm{exact}}$ — that is, the ratio of the eigenvalues from the original ensemble to the exact eigenvalues from the underlying distribution. The (solid) red line in the plots shows this curve; it agrees well with the bootstrap estimates.

The vertical (dotted) red lines in each figure show the position where the ratio curves intersect with $1-\sqrt{2/N_{G}}$ (bottom dotted line, see Eq. (69)). We set the SVD cutoff at this point in each case. Fitting these data we find that

[TABLE]

for $N_{G}=50$, showing that the SVD cut has a minimal effect (as expected), while for $N_{G}=500$ we have

[TABLE]

which shows that the SVD cut is essential since $\chi^{2}/N_{G}=1.30$ is much too large for $N_{G}=500$ — it corresponds to a $p$-value of order $3\times 10^{-5}$.

D.3 Conceptual Framework

The nature of the SVD modification can be understood by representing the ensemble-averaged data as a vector of Gaussian random variables,

[TABLE]

where

[TABLE]

and the uncorrelated random variables $z_{n}$ satisfy:

[TABLE]

Here $\delta\mathbf{G}$ represents the uncertainty associated with the ensemble average:

[TABLE]

The effect of the SVD cut is to add more uncertainty, $\delta\mathbf{G}^{\mathrm{SVD}}$:

[TABLE]

where

[TABLE]

and $\tilde{z}_{n}$ are new random variables with zero mean and unit covariance matrix. Then

[TABLE]

is the SVD-modified covariance data. The SVD noise discussed above is a random sample drawn from the distribution described by $\delta\mathbf{G}_{\mathrm{SVD}}$.

These formulas underscore the fact that introducing an SVD cut is a conservative move: it always increases the uncertainties in the data. This would not necessarily be the case if we renormalized the eigenvalues after introducing the SVD cut, as is done in Michael (1994); Michael and McKerrell (1995). In practice, however, the difference between the two approaches is small.

Finally note that another option in an SVD analyses is to discard modes below the cutoff. This corresponds to setting $\lambda_{n}=\infty$ for these modes, which is much larger than is reasonable, much too conservative. We find that fits are more accurate and more stable using the prescription outlined above.

D.4 Goodness of Fit

Note that $\chi^{2}/N_{G}=0.41$ in Eq. (72) is much smaller than expected for $N_{G}=500$: one expects $1.00(6)$ instead. The small value arises because random fluctuations in $\overline{\mathbf{G}}$ are characteristic of the uncertainties in $\delta\mathbf{G}$, but not those in $\delta\mathbf{G}_{\mathrm{SVD}}$. We can demonstrate this by adding a random sample to $\overline{\mathbf{G}}$ drawn from the distribution specified by $\delta\mathbf{G}_{\mathrm{SVD}}$,

[TABLE]

and refitting. In the case of Eq. (72), a typical fit with SVD noise gives $\chi^{2}/N_{G}$ increases to 0.96, which is consistent with expectations.

Parenthetically, we note that overly broad priors — for example, $0\pm 10$ for a set of parameters that are all order 1 — can also result in a small $\chi^{2}$. This situation is addressed in a similar fashion, by replacing the prior distribution $\mathbf{P}$:

[TABLE]

A good fit should have $\chi^{2}/N_{G}\approx 1\pm\sqrt{2/N_{G}}$ when both SVD and prior noise is included, and the fit results should agree (within errors) with the results without noise.

A more direct test of a fitting protocol than adding extra noise is to replace the fit data (Eq. (73)) with simulated data,

[TABLE]

which has the same covariance matrix (from $\delta\mathbf{G}$) as the real data, but whose mean is a random sample drawn from a distribution whose mean is known ($\mathbf{G}(\mathbf{p}_{\mathrm{sim}})$). A good fit to this simulated data should give best-fit results for the parameters that agree with $\mathbf{p}_{\mathrm{sim}}$ to within errors. An obvious choice for the simulation parameters $\mathbf{p}_{\mathrm{sim}}$ are the best-fit results obtained when fitting the real data.