$B\to D^{*}\ell{\nu_{\ell}}$ semileptonic form factors from lattice QCD with Möbius domain-wall quarks

Our gauge ensembles are generated for 2+1-flavor lattice QCD, where dynamical up and down quarks are degenerate and dynamical strange quark has a distinct mass. We take three lattice spacings of $a\!\simeq\!0.044,0.055$ and 0.080 fm Colquhoun et al. (2022), which correspond to the lattice cutoffs $a^{-1}\sim 2.453(4)$, 3.610(9) and 4.496(9) GeV, respectively. They are determined using the Yang-Mills gradient flow Borsanyi et al. (2012). The tree-level improved Symanzik gauge action Weisz (1983); Weisz and Wohlert (1984); Luscher and Weisz (1985) and the Möbius domain-wall quark action Brower et al. (2017) are used to control discretization errors and to preserve chiral symmetry to a good precision. We refer the interested reader to Refs. Brower et al. (2017); Colquhoun et al. (2022); Boyle (2015) for the five dimensional formulation of the quark action and our implementation. Its four-dimensional effective Dirac operator is given by

$\displaystyle\frac{1+m_{q}}{2}+\frac{1-m_{q}}{2}\gamma_{5}\,\epsilon\left[H_{M}\right],$

(8)

where $m_{q}$ represents the quark mass. We employ the kernel operator $H_{M}\!=\!2\gamma_{5}D_{W}(-M)/\left\{2+D_{W}(-M)\right\}$, where $D_{W}(-M)$ is the Wilson-Dirac operator with a negative mass with $M\!=\!1$ and with stout smearing Morningstar and Peardon (2004) applied three successive times to the gauge links. With our choice of the staple weight ($\rho\!=\!0.1$), the smearing radius $\sim\!1.5\,a$ is at the scale of the lattice spacing and vanishes in the limit of $a\!\to\!0$. We therefore expect that the link smearing does not change the continuum limit of the $B\!\to\!D^{*}$ form factors. It may induce additional discretization effects, which, however, do not turn out to be large in our continuum and chiral extrapolation in Sec. III. The polar decomposition approximation in Eq. (8) is realized for the sign function $\epsilon\left[H_{M}\right]$ in the five-dimensional domain-wall implementation. With this choice, the residual quark mass, which is a measure of the chiral symmetry violation, is suppressed to the level of 1 MeV at $a\!\simeq\!0.080$ fm, and 0.2 MeV or less on finer lattices, which are significantly smaller than the physical up and down quark masses Colquhoun et al. (2022). With reasonably small values of the fifth-dimensional size $N_{5}\!=\!8$ – 12, the computational cost is largely reduced from that of our previous large-scale simulation Aoki et al. (2008a); Aoki et al. (2012) using the overlap formulation that also preserves chiral symmetry Kaneko et al. (2014).

We employ the Möbius domain-wall action for all relevant quark flavors. Our choices of the degenerate up and down quark mass $m_{ud}$ correspond to pion masses from $M_{\pi}\!\sim\!500$ MeV down to 230 MeV. Chiral symmetry allows us to use heavy meson chiral perturbation theory (HMChPT) Randall and Wise (1993); Savage (2002) for our chiral extrapolation of simulation results to the physical pion mass without introducing terms to describe discretization effects. We take a strange quark mass $m_{s}$ close to its physical value fixed from $M_{\eta_{s}}^{2}\!=\!2M_{K}^{2}-M_{\pi}^{2}$. The charm quark mass $m_{c}$ is set to its physical value fixed from the spin-averaged charmonium mass $(M_{\eta_{c}}+3M_{J/\Psi})/4$ Nakayama et al. (2016). Depending on the lattice spacing $a$, we take three to six bottom quark masses $m_{Q}\!=\!1.25^{n}m_{c}$ $(n\!=\!0,1,\cdots)$ satisfying the condition $m_{Q}\!\leq\!0.7a^{-1}$ to avoid large discretization errors. We note that, with chiral symmetry, the leading discretization error is reduced to $O((am_{Q})^{2})$. As we will see in the next section, $a$ and $m_{Q}$ dependences of the form factors are mild, and the extrapolation to the continuum limit and physical quark masses is under good control. Our simulation parameters are summarized in Tables 1 and 2.

The spatial lattice size $L=N_{s}a$ is chosen to satisfy the condition $M_{\pi}L\!\gtrsim\!4$ to suppress finite volume effects (FVEs). At our smallest $M_{\pi}$ on the coarsest lattice, we also simulate a smaller physical lattice size $M_{\pi}L\!\sim\!3$ to directly examine the FVEs.

The statistics at the lattice cutoffs $a^{-1}\simeq\!2.5$, 3.6 and 4.5 GeV, are 5,000, 2,500 and 2,500 Hybrid Monte Carlo trajectories with a unit trajectory length of 1, 2 and 4 times the molecular dynamics time unit, respectively, which is increased to take account of the longer auto correlation toward the continuum limit. On each ensemble, we take $N_{\rm conf}$ configurations in an equal trajectory interval to calculate correlation functions relevant to the $B\!\to\!D^{*}\ell\nu$ decay. Our choice of $N_{\rm conf}$ is listed in Table 1. More details on our simulation setup can be found in Ref. Colquhoun et al. (2022).

The $B\!\to\!D^{*}$ matrix elements can be extracted from the three- and two-point functions

	$\displaystyle C_{\mathcal{O}_{\Gamma}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf p},{{\bf p}^{\prime}},{\epsilon^{\prime}})$	$\displaystyle=$	$\displaystyle\frac{1}{N_{s}^{3}N_{{t_{\rm src}}}}\sum_{{t_{\rm src}},{{\bf x}_{\rm src}}}\sum_{{\bf x},{{\bf x}^{\prime}}}\left\langle{\mathcal{O}}_{D^{*}}({{\bf x}^{\prime}},{t_{\rm src}}+{\Delta t}+{\Delta t^{\prime}};{\epsilon^{\prime}})\right.$		(9)
			$\displaystyle\hskip 28.45274pt\times\left.{\mathcal{O}}_{\Gamma}({\bf x},{t_{\rm src}}+{\Delta t}){\mathcal{O}}_{B}({{\bf x}_{\rm src}},{t_{\rm src}})^{\dagger}\right\rangle e^{-i{\bf p}({\bf x}-{{\bf x}_{\rm src}})-i{{\bf p}^{\prime}}({{\bf x}^{\prime}}-{\bf x})},\hskip 28.45274pt$
	$\displaystyle C^{P}({\Delta t};{\bf p})$	$\displaystyle=$	$\displaystyle\frac{1}{N_{s}^{3}N_{{t_{\rm src}}}}\sum_{{t_{\rm src}},{{\bf x}_{\rm src}}}\sum_{{\bf x}}\langle{\mathcal{O}}_{P}({\bf x},{t_{\rm src}}+{\Delta t}){\mathcal{O}}_{P}({{\bf x}_{\rm src}},{t_{\rm src}})^{\dagger}\rangle e^{-i{\bf p}({\bf x}-{{\bf x}_{\rm src}})}$		(10)

through their ground state contribution

	$\displaystyle C_{\mathcal{O}_{\Gamma}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf p},{{\bf p}^{\prime}},{\epsilon^{\prime}})$	$\displaystyle\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}$	$\displaystyle\frac{Z_{D^{*}}^{*}({{\bf p}^{\prime}},{\epsilon^{\prime}})\,Z_{B}({\bf p})}{4E_{D^{*}}({{\bf p}^{\prime}})E_{B}({\bf p})}\langle D^{*}({p^{\prime}},{\epsilon^{\prime}})|{\mathcal{O}}_{\Gamma}|B(p)\rangle e^{-E_{D^{*}}({{\bf p}^{\prime}}){\Delta t^{\prime}}-E_{B}({\bf p}){\Delta t}},\hskip 28.45274pt$
	$\displaystyle C^{P}({\Delta t};{\bf p})$	$\displaystyle\xrightarrow[{\Delta t}\to\infty]{}$	$\displaystyle\frac{Z_{P}^{*}({\bf p})\,Z_{P}({\bf p})}{2E_{P}({\bf p})}e^{-E_{P}({\bf p}){\Delta t}},$		(11)

where ${\mathcal{O}}_{\Gamma}$ is the weak vector ($V_{\mu}$) or axial ($A_{\mu}$) current, and we omit the symbol “$\prime$” on the momentum variable and the argument ${\epsilon^{\prime}}$ for the $D^{*}$ two-point function $C^{D^{*}}({\Delta t};{{\bf p}^{\prime}},{\epsilon^{\prime}})$ for simplicity. The meson interpolating field is denoted by ${\mathcal{O}}_{P}$ ($P\!=\!B,D^{*}$), where Gaussian smearing is applied to enhance their overlap with the ground state $Z_{P}({\bf p})\!=\!\langle P|{\mathcal{O}}_{P}^{\dagger}\rangle$. The $B$ meson is at rest (${\bf p}\!=\!{\bf 0}$) throughout this paper. The $w$ dependence of the form factors is studied by varying the three-momentum of $D^{*}$ as $|{{\bf p}^{\prime}}|^{2}\!=\!0,1,2,3,4$ (in this paper, we denote the momentum on the lattice in units of $2\pi/L$). To this end, we also calculate $D^{*}$ two-point functions with the local sink operator ${\mathcal{O}}_{D^{*},\rm lcl}^{\dagger}$

$\displaystyle C^{D^{*}}_{\rm sl}({\Delta t};{{\bf p}^{\prime}},{\epsilon^{\prime}})$

$\displaystyle=$

$\displaystyle\frac{1}{N_{s}^{3}N_{{t_{\rm src}}}}\sum_{{t_{\rm src}},{{\bf x}_{\rm src}}}\sum_{{\bf x}}\langle{\mathcal{O}}_{D^{*},\rm lcl}({\bf x},{t_{\rm src}}+{\Delta t};{\epsilon^{\prime}}){\mathcal{O}}_{D^{*}}({{\bf x}_{\rm src}},{t_{\rm src}};{\epsilon^{\prime}})^{\dagger}\rangle e^{-i{{\bf p}^{\prime}}({\bf x}-{{\bf x}_{\rm src}})},\hskip 22.76219pt$

(12)

which are used in a correlator ratio (15) below.

The three-point functions are calculated by using the sequential source method, where the total temporal separation ${\Delta t}+{\Delta t^{\prime}}$ between the source and sink operators is fixed and the temporal location of the weak current is varied. In order to control the contamination from the excited states, we repeat our measurement for four values of the source-sink separation ${\Delta t}+{\Delta t^{\prime}}$ listed in Table 1.

The statistical accuracy of the two- and three-point functions are improved by averaging over the spatial location of the source operator ${{\bf x}_{\rm src}}$ as indicated in Eqs. (9) and (10). To this end, we employ a volume source with $Z_{2}$ noise. At $M_{\pi}\!\lesssim\!300$ MeV, we repeat our measurement over two values of the source time-slices ${t_{\rm src}}\!=\!0$ and $N_{t}/2$ to take the average of the correlators. Namely $N_{t_{\rm src}}$ in Eqs. (9) and (10) is 2 at $M_{\pi}\!\lesssim\!300$ MeV, and 1 otherwise. The correlators with non-zero momentum ${\bf p}$ are also averaged over all possible ${\bf p}$’s based on parity and rotational symmetries on the lattice.

The number of measurements on each ensemble is given as $N_{\rm conf}N_{t_{\rm src}}$ in Table 1. As mentioned above, correlation functions for each configuration are averaged over $N_{t_{\rm src}}$ values of the source time-slice ${t_{\rm src}}$. Then, simulation data of $N_{\rm conf}$ configurations are divided into 50 bins: namely, the bin size is two (one) configurations for $\beta\!=\!4.17$ (4.35 and 4.47). The statistical error is estimated by the bootstrap method with 5,000 replicas. Figure 1 shows the bin size dependence of the relative statistical error of a three-point function at our largest cutoff, where the topological charge $Q$ changes much less frequently leading to larger auto-correlation compared to the coarser lattices Colquhoun et al. (2022). Our choice of the bin size on this finest lattice is one configuration, and we do not observe significant increase of the statistical error toward larger bin sizes. This suggests that our bin size is sufficiently large partly due to our choice of larger unit trajectory length towards the continuum limit. We also note that the error estimate is stable when we vary the number of bootstrap replicas as well as when we employ the jackknife method instead.

As observed in our previous simulation in the trivial topological sector Aoki et al. (2008a), the local topological fluctuation are active even if we fix the “global” topological charge $Q$, namely the whole-volume integral of the topological charge density. In Ref. Aoki et al. (2008b), we demonstrated that the topological susceptibility $\chi_{t}$ is calculable within the fixed topology setup. Our result for $\chi_{t}$ shows quark mass dependence consistent with ChPT, and its value extrapolated to the chiral limit is in good agreement with that from our simulation in the $\epsilon$-regime Fukaya et al. (2007a); Fukaya et al. (2007b). References Brower et al. (2003); Aoki et al. (2007) argued that the bias due to the freezing of the global topology can be considered as FVEs suppressed by the inverse space-time volume. In our study of the $B\!\to\!\pi\ell\nu$ decay Colquhoun et al. (2022), we confirmed that such FVEs on the pion effective mass are well suppressed, so that its $|Q|$ dependence is not significant. Figure 2 shows the statistical average of a $B\!\to\!D^{*}$ three-point function calculated for each $|Q|$. We do not observe any significant $|Q|$ dependence of the average at $|Q|\!\leq\!4$. At larger $|Q|$, we do not have enough data for a reliable error estimate, but individual data plotted by crosses are consistent with the averages at $|Q|\!\leq\!4$. We therefore conclude that the topology freezing effect is insignificant within our statistical accuracy.

To precisely extract the form factors, we construct ratios of the correlation functions, in which unnecessary overlap factors and exponential damping factors cancel for the ground-state contribution Hashimoto et al. (1999). The statistical fluctuation is also expected to partly cancel in such a ratio. With $B$ and $D^{*}$ mesons at rest, the vector current matrix element (6) vanishes, and the axial vector one (7) is sensitive only to the axial vector form factor $h_{A_{1}}$ at zero recoil, which is the fundamental input to determine $|V_{cb}|$. In order to precisely determine $h_{A_{1}}(1)$, we employ a double ratio

$\displaystyle R_{1}({\Delta t},{\Delta t^{\prime}})$

$\displaystyle=$

$\displaystyle\frac{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{\bf 0},{\epsilon^{\prime}})\,C_{A_{1}}^{D^{*}B}({\Delta t},{\Delta t^{\prime}};{\bf 0},{\bf 0},{\epsilon^{\prime}})}{C_{V_{4}}^{BB}({\Delta t},{\Delta t^{\prime}};{\bf 0},{\bf 0})\,C_{V_{4}}^{D^{*}D^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{\bf 0},{\epsilon^{\prime}},{\epsilon^{\prime}})}\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}|h_{A_{1}}(1)|^{2},$

(13)

where the $D^{*}$ polarization is chosen as ${\epsilon^{\prime}}\!=\!(1,0,0,0)$ along the polarization of the current inserted, i.e. $A_{1}$. The left panel of Fig. 3 shows our results for this ratio at our smallest $M_{\pi}$ and $a^{-1}$ with $m_{Q}\!=\!1.25\,m_{c}$. Note that the ratio is symmetrized as $R_{1}({\Delta t},{\Delta t^{\prime}})\!\to\!\left\{R_{1}({\Delta t},{\Delta t^{\prime}})+R_{1}({\Delta t^{\prime}},{\Delta t})\right\}/2$, since we use the same smearing function for both the source and sink operators. We observe a reasonable consistency among data at intermediate values of the source-sink separation ${\Delta t}+{\Delta t^{\prime}}\!=\!12$ (circles) and 16 (squares). The excited state contribution could potentially be significant but not large ($\lesssim\,1$ %) for the smallest ${\Delta t}+{\Delta t^{\prime}}\!=\!8$ (diamonds). Data at the largest ${\Delta t}+{\Delta t^{\prime}}\!=\!20$ (triangles) show a long plateau and are consistent with those at smaller ${\Delta t}+{\Delta t^{\prime}}$ within the large statistical errors. These observations suggest that the data at ${\Delta t}+{\Delta t^{\prime}}\!=\!12$, which is roughly 1 fm, are dominated by the ground state contribution at least at the midpoint ${\Delta t}\!\sim\!{\Delta t^{\prime}}$. The situation is similar in the right panel of the same figure for our largest $a^{-1}$ and with $m_{Q}\!=\!1.25^{4}m_{c}$, where data around midpoint ${\Delta t}\!\sim\!{\Delta t^{\prime}}$ are consistent among all simulated source-sink separations within 2 $\sigma$.

We find that all data are well described by the following fitting form including effects from the first excited states

$\displaystyle R_{1}({\Delta t},{\Delta t^{\prime}})$

$\displaystyle=$

$\displaystyle|h_{A_{1}}(1)|^{2}\left(1+ae^{-\Delta M_{B}{\Delta t}}+be^{-\Delta M_{D^{*}}{\Delta t}}+ae^{-\Delta M_{B}{\Delta t^{\prime}}}+be^{-\Delta M_{D^{*}}{\Delta t^{\prime}}}\right)$

(14)

with $\chi^{2}/{\rm d.o.f.}\!\lesssim\!1$ for a wide range of the values of ${\Delta t}$ and ${\Delta t^{\prime}}$. Here $\Delta M_{B(D^{*})}$ represents the energy difference between the $B(D^{*})$ meson ground state and the first excited state of the same quantum numbers, and is set to the value estimated from a two-exponential fit to the two-point function $C^{B(D^{*})}({\Delta t};{\bf 0})$. The fit range is chosen such that all the data of $R_{1}({\Delta t},{\Delta t^{\prime}})$ with the source-to-current (current-to-sink) separation ${\Delta t}^{(\prime)}$ equal to or larger than a lower cut ${\Delta t}^{(\prime)}_{\rm min}$ are included irrespective of the source-sink separation ${\Delta t}+{\Delta t^{\prime}}$. As shown in the left panel of Fig. 4, our result for $h_{A_{1}}(1)$ is stable against the choice of the lower cuts ${{\Delta t}_{\rm min}}$ and ${{\Delta t^{\prime}}_{\rm min}}$.

The right panels of Figs. 3 and 4 for our largest $a^{-1}$ also show ground state dominance with ${\Delta t}+{\Delta t^{\prime}}\!\gtrsim\!23a\!\sim\!1$ fm and stability of $h_{A_{1}}(1)$ against the choice of the fit range, respectively. The situation is similar for other simulation parameters.

In order to extract form factors at non-zero recoils, the $D^{*}$ momentum and polarization vector need to be chosen appropriately. The axial vector matrix element (7) is sensitive only to $h_{A_{1}}$ with a polarization vector ${\boldsymbol{\epsilon}}^{\prime}\!=\!(1,0,0)$ and a momentum ${{\bf p}^{\prime}_{\perp}}$ that satisfies ${\epsilon^{\prime*}}v\!=\!0$. We note that, here and in the following, three dimensional polarization vector ${\boldsymbol{\epsilon}}$ is accompanied by its temporal component to be fixed from the convention ${\epsilon^{\prime*}}p^{\prime}_{\perp}\!=\!0$, which is assumed for Eqs. (6) and (7). Then, the recoil parameter dependence of $h_{A_{1}}$ may be studied from the following ratio

$\displaystyle R_{1}^{\prime}({{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})=\frac{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})\,C^{D^{*}}_{sl}({\Delta t^{\prime}};{\bf 0},{\epsilon^{\prime}})}{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{\bf 0},{\epsilon^{\prime}})\,C^{D^{*}}_{sl}({\Delta t^{\prime}},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})}\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}\frac{w+1}{2}\frac{h_{A_{1}}(w)}{h_{A_{1}}(1)}.$

(15)

Here we use the two-point function with the local sink $C^{D^{*}}_{\rm sl}$ (12) as in our study of the $K\!\to\!\pi\ell{\nu_{\ell}}$ decay Aoki et al. (2017). The vector form factor $h_{V}$ can be extracted from the following ratio through the vector current matrix element (6)

$\displaystyle R_{V}({{\bf p}^{\prime\prime}_{\perp}},{{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})$

$\displaystyle=$

$\displaystyle\frac{C_{V_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime\prime}_{\perp}},{\epsilon^{\prime\prime}})}{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})}\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}\frac{i\varepsilon_{1ij}\epsilon^{\prime\prime*}_{i}v_{\perp j}^{\prime\prime}}{w+1}\frac{h_{V}(w)}{h_{A_{1}}(w)}.$

(16)

Here we use two different polarization vectors, ${\boldsymbol{\epsilon}}^{\prime}\!=\!(1,0,0)$ and ${\boldsymbol{\epsilon}}^{\prime\prime}\!=\!(0,1,0)$, and momenta ${\bf p}^{\prime}_{\perp}$ and ${\bf p}^{\prime\prime}_{\perp}$ that satisfy ${\bf p}^{\prime(\prime\prime)}_{\perp}\!\perp\!{\boldsymbol{\epsilon}}^{\prime(\prime\prime)}$. Note that $|{{\bf p}^{\prime\prime}_{\perp}}|\!=\!|{{\bf p}^{\prime}_{\perp}}|$ to share the same recoil parameter $w$, and $v_{\perp j}^{\prime\prime}\!=\!p_{\perp j}^{\prime\prime}/M_{D^{*}}$. We emphasize that renormalization factors of the weak vector and axial currents cancel thanks to chiral symmetry preserved in our simulations.

Figures 5 and 6 show an example of our results for $R_{1}^{\prime}$ and $R_{V}$ at our largest cutoff $a^{-1}\!\simeq\!4.5$ GeV and $m_{Q}\!=\!1.25^{4}m_{c}$. Around the mid-point ${\Delta t}\!\sim\!{\Delta t^{\prime}}$, we observe good consistency in these ratios among all simulated values of the source-sink separation ${\Delta t}+{\Delta t^{\prime}}\!\gtrsim\!1$ fm. A similar ground state dominance is also observed at other simulation points. We carry out a simultaneous fit using a fitting form that takes account of the first excited state contribution as

$\displaystyle R_{1}^{\prime}({{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})$

$\displaystyle=$

$\displaystyle R_{1}^{\prime}({{\bf p}^{\prime}_{\perp}})\left(1+ae^{-\Delta M_{B}{\Delta t}}+be^{-\Delta M_{D^{*}}{\Delta t^{\prime}}}\right)$

(17)

to extract the form factor ratio $h_{A_{1}}(w)/h_{A_{1}}(1)$ from the ground state contribution $R_{1}^{\prime}({{\bf p}^{\prime}_{\perp}})$, and a similar form for $R_{V}({{\bf p}^{\prime\prime}_{\perp}},{{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})$ to extract $h_{V}(w)/h_{A_{1}}(w)$. Our data are well described by this fitting form with $\chi^{2}/{\rm d.o.f.}\!\lesssim\!1$. From these form factor ratios and $h_{A_{1}}(1)$ from $R_{1}$, we calculate $h_{A_{1}}(w)$ and $h_{V}(w)$ at simulated values of $w$.

There is a 2 $\sigma$ bump of $R_{1}^{\prime}$ for ${\Delta t}\!+\!{\Delta t^{\prime}}\!=\!40$ (black triangles) around ${\Delta t}\!=\!28$ (${\Delta t^{\prime}}\!=\!12$) in Fig. 5. Since the global topological charge does not fluctuate very much at the corresponding cutoff, one may expect bump(s) in the exponential decay of relevant correlators due to instantons frozen at the temporal separation ${\Delta t}^{(\prime)}$ for a certain period of our simulation time. However, Fig. 7 shows that there is no statistically significant bump of the relevant correlators at ${\Delta t}\!\sim\!28$ and ${\Delta t^{\prime}}\!\sim\!12$. As mentioned above, the local topological fluctuation is active even when the global topological charge is fixed. We attribute 1 – 2 $\sigma$ bumps of correlator ratios in Figs. 3, 5 – 6, 8 – 9 to large statistical fluctuations of correlators at the largest temporal separation ${\Delta t}\!+\!{\Delta t^{\prime}}\!=\!40$, and accidental anti correlation of the statistical fluctuations. Since these data have larger statistical error than those at smaller ${\Delta t}\!+\!{\Delta t^{\prime}}$, these bumps do not change results of the fits (14) and (17) to extract the ground state contribution.

The axial vector form factors $h_{A_{2}}$ and $h_{A_{3}}$ can be extracted from the axial vector matrix element (7) with the $D^{*}$ spatial momentum ${{\bf p}^{\prime}_{\not\perp}}$ not perpendicular to the spatial polarization vector. The matrix element of $A_{1}$, for instance, has non-zero sensitivity to $h_{A_{1}}$ and $h_{A_{3}}$. We use the following ratio to extract $h_{A_{3}}$

	$\displaystyle R_{3}({{\bf p}^{\prime}_{\not\perp}},{{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})$	$\displaystyle=$	$\displaystyle r_{Z}\frac{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\not\perp}},{\epsilon^{\prime\prime}})}{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})}$		(18)
		$\displaystyle\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}$	$\displaystyle{\epsilon^{\prime\prime}}^{*}_{1}-\frac{{\epsilon^{\prime\prime}}^{*}_{4}{v^{\prime}}_{\not\perp,1}}{w+1}\frac{h_{A_{3}}(w)}{h_{A_{1}}(w)},$

where ${\boldsymbol{\epsilon}}^{\prime}\!=\!(1,0,0)$, ${\boldsymbol{\epsilon}}^{\prime*}{{\bf p}^{\prime}_{\perp}}\!=0$, ${\boldsymbol{\epsilon}}^{\prime\prime*}{{\bf p}^{\prime}_{\not\perp}}\!\neq\!0$, and $|{{\bf p}^{\prime}_{\not\perp}}|\!=\!|{{\bf p}^{\prime}_{\perp}}|$. The factor $r_{Z}\!=\!Z_{D^{*}}({{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})/Z_{D^{*}}({{\bf p}^{\prime}_{\not\perp}},{\epsilon^{\prime\prime}})$ appears, since the overlap factor of $D^{*}$ depends on whether the polarization vector is perpendicular to the momentum, even if $|{{\bf p}^{\prime}_{\not\perp}}|\!=\!|{{\bf p}^{\prime}_{\perp}}|$. While $r_{Z}$ can be estimated from individual fit to the relevant two-point functions, we employ a ratio $r_{Z}\!=\!\sqrt{C^{D^{*}}({\Delta t}_{\rm ref},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})/C^{D^{*}}({\Delta t}_{\rm ref},{{\bf p}^{\prime}_{\not\perp}},{\epsilon^{\prime\prime}})}$, which shows better stability against the temporal separation. The reference temporal separation is chosen as ${\Delta t}_{\rm ref}\!\sim\!T/4$ by inspecting the statistical accuracy and the ground state saturation of $r_{Z}$. A similar ratio but with $A_{4}$ in the numerator is also sensitive to $h_{A_{2}}$

	$\displaystyle R_{2}({{\bf p}^{\prime}_{\not\perp}},{{\bf p}^{\prime}_{\perp}};{\Delta t},{\Delta t^{\prime}})$	$\displaystyle=$	$\displaystyle r_{Z}\frac{C_{A_{4}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\not\perp}},{\epsilon^{\prime\prime}})}{C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}_{\perp}},{\epsilon^{\prime}})}$		(19)
		$\displaystyle\xrightarrow[{\Delta t},{\Delta t^{\prime}}\to\infty]{}$	$\displaystyle{\epsilon^{\prime\prime}}^{*}_{4}\left[1-\frac{1}{w+1}\left\{\frac{h_{A_{2}}(w)}{h_{A_{1}}(w)}+v_{\not\perp,4}\frac{h_{A_{3}}(w)}{h_{A_{1}}(w)}\right\}\right].$

We plot our results for $R_{3}$ and $R_{2}$ on our finest lattice in Figs. 8 and 9, respectively. Similar to other ratios, the excited state contamination is not large. The fit similar to (17) yields the ground state contributions, $R_{3}({{\bf p}^{\prime}_{\not\perp}},{{\bf p}^{\prime}_{\perp}})$ and $R_{2}({{\bf p}^{\prime}_{\not\perp}},{{\bf p}^{\prime}_{\perp}})$, which are stable against the choice of the fit range: we do not observe ${\Delta t}^{(\prime)}_{\rm min}$ dependence beyond 1 $\sigma$ level when ${{\Delta t}_{\rm min}}+{{\Delta t^{\prime}}_{\rm min}}\!\lesssim\!16$, that is our smallest value of the source-sink separation ${\Delta t}+{\Delta t^{\prime}}$. The error rapidly increases at larger values of ${\Delta t}^{(\prime)}_{\rm min}$, because less data are available for the fit. From $R_{2}$, $R_{3}$ and $h_{A_{1}}$ extracted above, we can calculate $h_{A_{2}}$ and $h_{A_{3}}$.

For one of the $D^{*}$ momenta ${{\bf p}^{\prime}}\!=\!(1,1,1)$, it is not straightforward to use the above-mentioned correlator ratios, which use $C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}},{\epsilon^{\prime}})$ with ${\boldsymbol{\epsilon}}^{\prime*}{{{\bf p}^{\prime}}}\!=\!0$ as a correlator sensitive only to $h_{A_{1}}(w)$. In this case, we use a linear combination

$\displaystyle C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}},{\epsilon^{\prime}})-C_{A_{1}}^{BD^{*}}({\Delta t},{\Delta t^{\prime}};{\bf 0},{{\bf p}^{\prime}},{\epsilon^{\prime\prime}})$

$\displaystyle\propto$

$\displaystyle(w+1)h_{A_{1}}(w)$

(20)

which only involves $h_{A_{1}}$ with ${\boldsymbol{\epsilon}}^{\prime}\!=\!(1,0,0)$ and ${\boldsymbol{\epsilon}}^{\prime\prime}\!=\!(0,1,0)$.

The double ratio $R_{1}$ enables us to calculate $h_{A_{1}}(1)$ with 1 % statistical error or better, as it only involves correlators with zero momentum and the statistical fluctuation partially cancels between the numerator and denominator. At non-zero recoil, the numerator of $R_{1}^{\prime}$ ($R_{V}$) is exclusively sensitive to $h_{A_{1}}$ ($h_{V}$), and the statistical error of $h_{A_{1}}$ and $h_{V}$ is typically a few %. This is not the case for $h_{A_{2}}$ and $h_{A_{3}}$ with our setup where the $B$ meson is at rest. The numerator of $R_{3}$ involves $h_{A_{1}}$ and $h_{A_{3}}$, and the typical precision for $h_{A_{3}}$ is 10 – 20 %. The statistical error increases to $\gtrsim\!40$ % for $h_{A_{2}}$, since the central value is suppressed by heavy quark symmetry and $R_{2}$ is a linear combination of all axial form factors.

We note that, thanks to chiral symmetry, finite renormalization factors for the weak vector and axial currents cancel in the above correlator ratios (13), (15), (16), (18) and (19). While discretization effects to the wave function renormalization are not small Fahy et al. (2016), these also cancel in the ratios.

Our results for the form factors obtained at various simulation parameters are plotted as a function of the recoil parameter $w$ in Fig. 10. We find that the results with different heavy and light quark masses are all close to each other. It also turns out that the discretization effects are not large.

We extrapolate these results to the continuum limit and physical quark masses based on next-to-leading order heavy meson chiral perturbation theory (NLO HMChPT) Randall and Wise (1993); Savage (2002). The extrapolation form is generically written as

	$\displaystyle h_{X}/\eta_{X}$	$\displaystyle=$	$\displaystyle c+\frac{{g_{D^{*}D\pi}}^{2}}{2\Lambda_{\chi}^{2}}\bar{F}_{\rm log}\!\left({\xi_{\pi}},\Delta_{D},\Lambda_{\chi}\right)+c_{\pi}{\xi_{\pi}}+c_{\eta_{s}}{\xi_{\eta_{s}}}+c_{Q}{\epsilon_{Q}}$		(21)
			$\displaystyle+c_{a}{\xi_{a}}+c_{am_{Q}}{\xi_{am_{Q}}}+c_{w}(w-1)+d_{w}(w-1)^{2}\hskip 14.22636pt(X\!=\!A_{1},A_{2},A_{3}{\rm\ and\ }V),\hskip 14.22636pt$

which comprises the constant term, chiral logarithm at $w\!=\!1$ and leading corrections in small expansion parameters

$\displaystyle{\xi_{\pi}}=\frac{M_{\pi}^{2}}{\Lambda_{\chi}^{2}},\hskip 5.69054pt{\xi_{\eta_{s}}}=\frac{M_{\eta_{s}}^{2}}{\Lambda_{\chi}^{2}},\hskip 5.69054pt{\epsilon_{Q}}=\frac{\bar{\Lambda}}{2m_{Q}},\hskip 5.69054pt{\xi_{a}}=(\Lambda a)^{2},\hskip 5.69054pt{\xi_{am_{Q}}}=(am_{Q})^{2},$

(22)

with $M_{\eta_{s}}^{2}\!=\!2M_{K}^{2}-M_{\pi}^{2}$, $\Lambda_{\chi}\!=\!4\pi f_{\pi}$, and nominal values of $\bar{\Lambda}\!=\!\Lambda\!=\!500$ MeV.

We denote the one-loop integral function for the chiral logarithm by $\bar{F}_{\rm log}\!\left({\xi_{\pi}},\Delta_{D},\Lambda_{\chi}\right)$, where $\Delta_{D}\!=\!M_{D^{*}}-M_{D}$ is the $D^{*}$-$D$ mass splitting. While the loop function can be approximated as $\bar{F}_{\rm log}\!\left({\xi_{\pi}},\Delta_{D},\Lambda_{\chi}\right)\!=\!\Delta_{D}^{2}\ln[{\xi_{\pi}}]+O(\Delta_{D}^{3})$, we use its explicit form in Ref. Savage (2002). For the $D^{*}D\pi$ coupling, we take a value of ${g_{D^{*}D\pi}}\!=\!0.53(8)$, which was employed in the Fermilab/MILC paper Bailey et al. (2014) and covers previous lattice estimates Detmold et al. (2012); Becirevic and Sanfilippo (2013); Can et al. (2013); Bernardoni et al. (2015); Flynn et al. (2016). This choice does not lead to a large systematic uncertainty, because the chiral logarithm is suppressed by the small mass splitting squared $\Delta_{D}^{2}$.

In Fig. 11, we plot $h_{A_{1}}$ and $h_{V}$ at our smallest value of $w$ as a function of $M_{\pi}^{2}$. The $M_{\pi}^{2}$ dependence is rather mild, possibly because there is no valence pion in this decay. We note that there is a non-trivial concave structure at small $M_{\pi}^{2}$ due to the opening of the $D^{*}\!\to\!D\pi$ decay. The effect estimated in NLO HMChPT is, however, not large compared to our statistical accuracy.

For our chiral extrapolation, we employ the so-called $\xi$ expansion in ${\xi_{\pi}}$ and ${\xi_{\eta_{s}}}$: namely, we use the measured meson masses squared, $M_{\pi}^{2}$ and $M_{\eta_{s}}^{2}$, instead of quark masses $m_{ud}$ and $m_{s}$, respectively, as well as quark-mass dependent $f_{\pi}$ for $\Lambda_{\chi}$ instead of the decay constant $f$ in the chiral limit. This is motivated by our observation Noaki et al. (2008) that the $\xi$-expansion shows better convergence of the chiral expansion for light meson observables compared to the $x$-expansion with $x\!=\!2Bm_{q}/(4\pi f)^{2}$, where $B$ gives the lowest order relation $M_{\pi}^{2}\!=\!2Bm_{ud}$.

The extrapolation form explicitly includes the one-loop radiative correction $\eta_{X}$ for the matching of the heavy-heavy currents between QCD and HQET Falk et al. (1990); Falk and Grinstein (1990); Neubert (1992). Remaining heavy quark mass dependence of the form factors is then parametrized by polynomials of ${\epsilon_{Q}}\!\propto\!1/m_{Q}$ and ${\xi_{am_{Q}}}\!\propto\!m_{Q}$. The dependence is mild as shown in Fig. 12. Our results are well described by the linear terms $c_{Q}{\epsilon_{Q}}$ and $c_{am_{Q}}{\xi_{am_{Q}}}$ in Eq. (21). We note that, for $h_{A_{1}}$, the term $c_{Q}{\epsilon_{Q}}$ is modified as $c_{Q}(w-1){\epsilon_{Q}}$ to be consistent with Luke’s theorem Luke (1990), and we add the quadratic term $d_{Q}{\epsilon_{Q}}^{2}$ to take account of possible $m_{Q}$ dependence at $w\!=\!1$.

The discretization effects in the form factors also turn out to be mild with our choice of parameters, $a^{-1}\!\gtrsim\!2.5$ GeV and $m_{Q}\!\leq 0.7a^{-1}$. These are described by an expansion, i.e. linear in the $m_{Q}$ independent and dependent parameters ${\xi_{a}}$ and ${\xi_{am_{Q}}}$, respectively. Figure 12 shows that the $m_{Q}$ dependent discretization error is significant in our results. However, the net $m_{Q}$ dependence is described reasonably well with our fitting form (21). We note that the $m_{Q}$ dependence for $h_{A_{2}}$ and $h_{A_{3}}$ is insignificant due to their large relative uncertainty. As we will mention later, fit results do not change significantly if we exclude data at $0.5\!\leq\!am_{Q}\!\leq\!0.7$ in our continuum and chiral extrapolation.

In this study, we first extrapolate the form factors to the continuum limit and physical quark masses, and then parametrize the $w$ dependence by using the model independent BGL parametrization, because it may not be simply combined with the chiral and heavy quark expansions. Concerning the $w$ dependence, therefore, the fit form (21) with an expansion in $w-1$ is used to interpolate our results. Figure 10 shows that our data do not show strong curvature in the simulation region of $1\!\leq\!w\!\lesssim\!1.1$, and can be well described with only up to the quadratic term.

Numerical results of our continuum and chiral extrapolation are summarized in Table 3. The polynomial coefficients turn out to be $O(1)$ or less with our choice of the dimensionless expansion parameters (22), which are normalized by appropriate scales. Many of them are consistent with zero reflecting the mild dependence on the lattice spacing and quark masses mentioned above. As a result, our choice of the fitting forms (21) results in small values of $\chi^{2}/{\rm d.o.f.}\!\lesssim\!0.2$, which, however, are only a rough measure of the goodness of the fits, because we employ approximated estimators of the covariance matrix as discussed in the next subsection.