From deep inelastic scattering to heavy-flavor semi-leptonic decays: Total rates into multi-hadron final states from lattice QCD

TL;DR

This paper introduces a novel lattice QCD technique to compute total decay and transition rates into multi-hadron states, applicable to various processes including semi-leptonic decays and deep inelastic scattering.

Contribution

The authors develop a new method using finite-volume Euclidean correlators and the Backus-Gilbert approach to extract total rates for complex hadronic final states from lattice QCD.

Findings

Demonstrated the method with a numerical example

Achieved reasonable precision for multiple decay channels

Applicable to nucleon structure functions and heavy-flavor decays

Abstract

We present a new technique for extracting decay and transition rates into final states with any number of hadrons. The approach is only sensitive to total rates, in which all out-states with a given set of QCD quantum numbers are included. For processes involving photons or leptons, differential rates with respect to the non-hadronic kinematics may also be extracted. Our method involves constructing a finite-volume Euclidean four-point function, whose corresponding spectral function measures the decay and transition rates in the infinite-volume limit. This requires solving the inverse problem of extracting the spectral function from the correlator and also necessitates a smoothing procedure so that a well-defined infinite-volume limit exists. Both of these steps are accomplished by the Backus-Gilbert method and, as we show with a numerical example, reasonable precision can be expected in cases with multiple open decay channels. Potential applications include nucleon structure functions and the onset of the deep inelastic scattering regime, as well as semi-leptonic $D$ and $B$ decay rates.

Equations

$$⟨E_k(L),\text{p}|\mathcal{J}(0)|\Psi ,\text{P}{⟩}_L=\sum _{\alpha }{\mathcal{C}}_{\alpha }⟨E_k(L),\text{p},\alpha ,\mathrm{out}|\mathcal{J}(0)|\Psi ,\text{P}⟩,$$

$$\Gamma =2\pi \underset{t\to \infty }{\lim }\underset{L\to \infty }{\lim }\sum _k|V_k(L){|}^2{\delta }_{\text{nicefrac}1t}(E_k(L)-M),$$

$${\delta }_{\text{nicefrac}1t}(\omega )=\frac{2}{\pi }\frac{{\sin }^2(\omega t/2)}{{\omega }^{2}t},$$

$$\Gamma$$

$$G(\tau ,L)={\int }_0^{\infty }\frac{d\omega }{2\pi }{e}^{-\omega \tau }\rho (\omega ,L),$$

$$\left[\int d^3\text{x}{\mathcal{H}}_{\mathrm{QCD}}(x)\right]|N,\text{P},\lambda ⟩=|N,\text{P},\lambda ⟩E_N,$$

$${\mathcal{A}}_{N(\lambda )\to \alpha }(E,\text{p})\equiv ⟨E,\text{p},\alpha ;\mathrm{out}|{\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda ⟩.$$

$$⟨E,\text{p};N\pi ,\text{k};\mathrm{out}|E^{\prime },{\text{p}}^{\prime };N\pi ,{\text{k}}^{\prime };\mathrm{out}⟩=2{\omega }_{N,\text{k}}2{\omega }_{\pi ,\text{p}-\text{k}}(2\pi {)}^6{\delta }^3(\text{k}-{\text{k}}^{\prime }){\delta }^3(\text{p}-\text{k}-{\text{p}}^{\prime }+{\text{k}}^{\prime }),$$

$${\omega }_{N,\text{k}}=\sqrt{M_N^2+{\text{k}}^2},{\omega }_{\pi ,\text{p}-\text{k}}=\sqrt{M_{\pi }^2+(\text{p}-\text{k}{)}^2}.$$

$${\mathcal{A}}_{D\to \alpha }\equiv ⟨E_D,\text{P},\alpha ;\mathrm{out}|{\mathcal{H}}_{\mathcal{Q}}(0)|D,\text{P}⟩.$$

$$M_{k,N(\lambda )\to \mathcal{Q}}(\text{p},L)$$

$$M_{k,D\to \mathcal{Q}}(L)$$

$$d{\Phi }_{\alpha }(k_1,\cdots ,k_{N_{\alpha }})\equiv \frac{d^3{\text{k}}_1}{(2\pi {)}^{3}2{\omega }_{{\text{k}}_1}}\cdots \frac{d^3{\text{k}}_{N_{\alpha }}}{(2\pi {)}^{3}2{\omega }_{{\text{k}}_{N_{\alpha }}}}(2\pi {)}^4{\delta }^4(P-{\sum }_{i=1}^{N_{\alpha }}{k}_i).$$

$${\Gamma }_{D\to \mathcal{Q}}\equiv \frac{1}{2M_D}\sum _{\alpha }\frac{1}{S_{\alpha }}\int d{\Phi }_{\alpha }(k_1,\cdots ,k_{N_{\alpha }})|⟨E_D,\text{P},\alpha ;\mathrm{out}|{\mathcal{H}}_{\mathcal{Q}}(0)|D,\text{P}⟩{|}^2,$$

$${\rho }_{\mathcal{Q},\text{P}}(E,\text{p})\equiv \frac{1}{n_{\lambda }}\sum _{\lambda ,\alpha }\frac{1}{S_{\alpha }}\int d{\Phi }_{\alpha }(k_1,\cdots ,k_{N_{\alpha }})|⟨E,\text{p},\alpha ;\mathrm{out}|{\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda ⟩{|}^2,$$

$$\begin{gathered} {\rho }_{\mathcal{Q},\text{P}}(E,\text{p})=\frac{1}{n_{\lambda }}\sum _{\lambda }\int d^{4}x{e}^{i(E-E_N)t-i(\text{p}-\text{P})\cdot \text{x}}\int \frac{d^3{\text{p}}^{\prime }}{(2\pi {)}^3}\int \frac{d{E}^{\prime }}{2\pi }{e}^{-i(E^{\prime }-E_N)t+i({\text{p}}^{\prime }-\text{P})\cdot \text{x}} \\ \times \sum _{\alpha }\frac{1}{S_{\alpha }}\int d{\Phi }_{\alpha }(k_1,\cdots ,k_{N_{\alpha }})⟨N,\text{P},\lambda |{\mathcal{J}}_{\mathcal{Q}}^{†}(0)|E^{\prime },{\text{p}}^{\prime },\alpha ;\mathrm{out}⟩⟨E^{\prime },{\text{p}}^{\prime },\alpha ;\mathrm{out}|{\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda ⟩. \end{gathered}$$

$$\begin{gathered} {\rho }_{\mathcal{Q},\text{P}}(E,\text{p})=\frac{1}{n_{\lambda }}\sum _{\lambda }\int d^{4}x{e}^{i(E-E_N)t-i(\text{p}-\text{P})\cdot \text{x}}\int \frac{d^3{\text{p}}^{\prime }}{(2\pi {)}^3}\int \frac{d{E}^{\prime }}{2\pi } \\ \times \sum _{\alpha }\frac{1}{S_{\alpha }}\int d{\Phi }_{\alpha }(k_1,\cdots ,k_{N_{\alpha }})⟨N,\text{P},\lambda |{\mathcal{J}}_{\mathcal{Q}}^{†}(x)|E^{\prime },{\text{p}}^{\prime },\alpha ;\mathrm{out}⟩⟨E^{\prime },{\text{p}}^{\prime },\alpha ;\mathrm{out}|{\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda ⟩, \end{gathered}$$

$${\rho }_{\mathcal{Q},\text{P}}(E,\text{p})=\frac{1}{n_{\lambda }}\sum _{\lambda }\int d^{4}x{e}^{i(E-E_N)t-i(\text{p}-\text{P})\cdot \text{x}}⟨N,\text{P},\lambda |{\mathcal{J}}_{\mathcal{Q}}^{†}(x){\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda ⟩.$$

$${\Gamma }_{D\to \mathcal{Q}}=\frac{1}{2M_D}{\rho }_{\mathcal{Q},\text{P}}(E_D,\text{P})=\frac{1}{2M_D}\int d^{4}x⟨D,\text{P}|{\mathcal{H}}_{\mathcal{Q}}(x){\mathcal{H}}_{\mathcal{Q}}(0)|D,\text{P}⟩.$$

$$G_{\mathcal{Q},\text{P}}(\tau ,\text{x},L)\equiv 2E_N{L}^6{e}^{-E_N\tau +i\text{P}\cdot \text{x}}\underset{{\tau }_f\to \infty }{\lim }\underset{{\tau }_i\to -\infty }{\lim }\frac{{\sum }_{\lambda }⟨{\Psi }_{\lambda }({\tau }_f,\text{P}){\mathcal{J}}_{\mathcal{Q}}^{†}(\tau ,\text{x}){\mathcal{J}}_{\mathcal{Q}}(0){\Psi }_{\lambda }^{†}({\tau }_i,\text{P}){⟩}_{\mathrm{conn}}}{{\sum }_{\lambda }⟨{\Psi }_{\lambda }({\tau }_f,\text{P}){\Psi }_{\lambda }^{†}({\tau }_i,\text{P})⟩},$$

$$G_{\mathcal{Q},\text{P}}(\tau ,\text{x},L)=2E_N{L}^3{e}^{-E_N\tau +i\text{P}\cdot \text{x}}\frac{1}{n_{\lambda }}\sum _{\lambda }⟨N,\text{P},\lambda |{\mathcal{J}}_{\mathcal{Q}}^{†}(\tau ,\text{x}){\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda {⟩}_L.$$

$${\tilde{G}}_{\mathcal{Q},\text{P}}(\tau ,\text{p},L)$$

$$=2E_N{L}^3{e}^{-E_N\tau }\int d^3\text{x}{e}^{-i(\text{p}-\text{P})\cdot \text{x}}\frac{1}{n_{\lambda }}\sum _{\lambda }⟨N,\text{P},\lambda |{\mathcal{J}}_{\mathcal{Q}}^{†}(\tau ,\text{x}){\mathcal{J}}_{\mathcal{Q}}(0)|N,\text{P},\lambda {⟩}_L.$$

$${\tilde{G}}_{\mathcal{Q},\text{P}}(\tau ,\text{p},L)$$

$${\tilde{G}}_{\mathcal{Q},\text{P}}(\tau ,\text{p},L)={\int }_0^{\infty }\frac{d\omega }{2\pi }{e}^{-\omega \tau }{\rho }_{\mathcal{Q},\text{P}}(\omega ,\text{p},L),$$

$$\rho_{\mathcal{Q},\textbf{P}}(E,\textbf{p},L)\equiv 2E_{N}\,L^{6}\sum_{k}|M_{k,N\to\mathcal{Q}}(\textbf{p},L)|^{2}\;2\pi\,\delta\big{(}E-E_{k}(L)\big{)}\,,$$

$${\int }_0^{\infty }d\omega {\hat{\delta }}_{\Delta }(\overline{\omega },\omega )=1,\underset{\Delta \to 0}{\lim }{\int }_0^{\infty }d\omega {\hat{\delta }}_{\Delta }(\overline{\omega },\omega )\varphi (\omega )=\varphi (\overline{\omega }),$$

$${\hat{\rho }}_{\mathcal{Q},\text{P}}(\overline{\omega },\text{p},L,\Delta )\equiv {\int }_0^{\infty }d\omega {\hat{\delta }}_{\Delta }(\overline{\omega },\omega ){\rho }_{\mathcal{Q},\text{P}}(\omega ,\text{p},L).$$

$${\rho }_{\mathcal{Q},\text{P}}(E,\text{p})=\underset{\Delta \to 0}{\lim }\underset{L\to \infty }{\lim }{\hat{\rho }}_{\mathcal{Q},\text{P}}(E,\text{p},L,\Delta ),$$

$${\hat{\delta }}_{\Delta }(\overline{\omega },\omega )=\sum _j{C}_j(\overline{\omega },\Delta )e^{-\omega {\tau }_j},$$

$$\Delta ={\int }_0^{\infty }d\omega (\bar{\omega }-\omega {)}^2{\hat{\delta }}_{\Delta }(\overline{\omega },\omega {)}^2,$$

Full text

deep inelastic scattering, weak decays, lattice QCD

††preprint: MITP/17-031

I Introduction

Reliably calculating the low-energy phenomenology of the strong interaction is very challenging. One major source of complication is that the underlying gauge theory, quantum chromodynamics (QCD), is most simply expressed in terms of color-charged quarks and gluons, whereas the low-energy degrees of freedom of the theory are confined color-singlet states, called hadrons. In the past few decades, great progress has been made in extracting the low-energy hadron spectrum and its structure directly from QCD, by numerically estimating the quantum path integral on a discretized, finite-volume Euclidean space-time.

This approach, called lattice QCD, enables one to calculate the discrete spectrum of QCD on a three-dimensional torus, i.e. a cubic spatial volume with periodicity $L$, as well as matrix elements of quark-field operators between finite-volume hamiltonian eigenstates. Taking discretization and finite-time effects to be negligible, it is important to understand how the energies and matrix elements on the torus can be related to the infinite-volume observables of the theory.

For example, the energies of single-particle states with vanishing spatial momentum are known to satisfy $E_{\mathcal{Q}}(L)=M_{\mathcal{Q}}+\mathcal{O}(e^{-M_{\pi}L})$, where $M_{\pi}$ is the physical mass of the lightest degree of freedom, the pion in QCD, and $M_{\mathcal{Q}}$ the mass of the particle with quantum numbers ${\mathcal{Q}}$ Lüscher (1986a). For such states, one aims to calculate $E_{\mathcal{Q}}(L)$ at multiple volumes and extrapolate to the infinite-volume limit. Similarly, local matrix elements with single-hadron states are exponentially close to their infinite-volume counterparts, so that one can hope to estimate the experimental observable by taking $M_{\pi}L$ sufficiently large.

By contrast, for multi-particle states, a naive infinite-volume limit is not useful. For a given box size, the energy of the $n^{\rm th}$ multi-particle state contains information about the interactions of the particles, but as $L\to\infty$ with $n$ fixed, this state flows to the threshold value, $2M_{\mathcal{Q}}$. Instead, the now-standard approach is to use finite-volume as a tool, rather than an unwanted artifact, by applying analytic, field-theoretic relations between finite- and infinite-volume quantities.

This approach was pioneered by Lüscher, who derived a quantization condition relating the finite-volume spectrum to the two-to-two scattering amplitude in the case of identical scalar particles Lüscher (1986b, 1991). The formalism includes all partial waves but assumes that the two-particle states have zero total momentum in the finite-volume frame and also neglects exponentially suppressed terms of the form $e^{-M_{\pi}L}$. More recently, this idea has been generalized to describe systems with nonzero momentum, as well as any number of strongly-coupled two-particle channels, including non-identical and non-degenerate particles as well as intrinsic spin Rummukainen and Gottlieb (1995); h Kim et al. (2005); Christ et al. (2005); Lage et al. (2009); Bernard et al. (2011); Döring et al. (2011); Fu (2012); Göckeler et al. (2012); Hansen and Sharpe (2012); Briceño and Davoudi (2013); Briceño (2014). Extensions to three-particle scattering are also underway Hansen and Sharpe (2014, 2015); Briceño et al. (2017).

A similar approach can be used to relate matrix elements of finite-volume multi-particle states to their infinite-volume counterparts, in order to extract electroweak decay and transition amplitudes. Here the paradigm example is the weak decay $K\to\pi\pi$. The decay amplitude is given at leading order in the weak interaction by the QCD matrix element $\langle\pi\pi,\mathrm{out}|\mathcal{H}(0)|K\rangle$, where $\mathcal{H}(0)$ is the weak hamiltonian density expressed in terms of four-quark operators. In lattice QCD, however, it is only possible to calculate $\langle n,L|\mathcal{H}(0)|K,L\rangle$ where $|n,L\rangle$ is a finite-volume state with the quantum numbers of two pions.

In Ref. Lellouch and Lüscher (2001), Lellouch and Lüscher derived a conversion factor relating these finite- and infinite-volume matrix elements. The factor depends on the box size via a known geometric function and also depends on the derivative of the elastic pion scattering phase shift, $\partial\delta_{\pi\pi}(k)/\partial k$, evaluated at the kaon mass. This relation is being used by the RBC/UKQCD collaboration to perform a full-error-budget calculation of $K\to\pi\pi$ decays, aiming towards a first-principles understanding of the $\Delta I=1/2$ rule and a prediction of $\epsilon^{\prime}/\epsilon$ Blum et al. (2015); Bai et al. (2015).

The Lellouch-Lüscher relation has since been generalized to states with non-zero total momentum in the finite volume and to coupled two-particle channels with non-identical and non-degenerate particles, as well as particles with intrinsic spin h Kim et al. (2005); Christ et al. (2005); Meyer (2011a); Hansen and Sharpe (2012); Briceño and Davoudi (2013); Briceño et al. (2015); Briceño and Hansen (2015). These extensions also accommodate currents that carry angular-momentum, momentum and energy, allowing one to extract time-like form factors as well as semi-leptonic decay amplitudes. In Ref. Briceño and Hansen (2016) the relation was also extended to matrix elements of the form $\langle\pi\pi,\mathrm{out}|\mathcal{J}(0)|\pi\pi,\mathrm{in}\rangle$, providing a rigorous path towards resonance form factors and transition amplitudes. Going beyond single-current insertions, Ref. Christ et al. (2015) used techniques based in the Lellouch-Lüscher approach to analyze long-distance contributions to $K_{L}$-$K_{S}$ mixing from a finite-volume Euclidean four-point function.

Aside from the two-to-two transitions and the neutral meson mixing formalism, all relations of the Lellouch-Lüscher type have the general form

[TABLE]

where $|\Psi,\textbf{P}\rangle$ is a QCD-stable, single-particle state with the indicated three-momentum. Here $\alpha$ is summed over all two-particle channels with the relevant quantum numbers. For example, in the case of unflavored states, the sum includes $\pi\pi$ and $K\overline{K}$. In words, the finite-volume matrix element on the left-hand side of Eq. (1) is equal to a linear combination of all possible infinite-volume matrix elements whose out-states have the appropriate quantum numbers. The coefficients, $\mathcal{C}_{\alpha}$, depend on derivatives of all parameters in this sector of the QCD scattering matrix as well as on the box size.

Using this result to extract decay and transition amplitudes for heavy mesons is extremely challenging. For a system with $N$ open decay channels, one must first use extensions of the Lüscher formalism to extract all QCD scattering parameters, at multiple energies near the target decay energy, in order to estimate the derivatives. Given these, it is possible to calculate the coefficients, $\mathcal{C}_{\alpha}$. In a second step, one must then calculate finite-volume matrix elements at different box sizes but with the same energies, in order to generate multiple equations of the form (1), with only the infinite-volume matrix elements unknown. Given $N$ independent results one can invert the relations to determine the various transition amplitudes.111In certain cases an ambiguity in the overall phase may require an additional constraint beyond the $N$ that are always needed, see also the discussion in Ref. Hansen and Sharpe (2012).

It is also important to note that this approach is currently only available for energies below the production of the lightest state with more than two hadrons. For example, in the case of unflavored final states, one is limited by either the three- or four-pion threshold. The situation is particularly frustrating for the decay of charmed or bottom mesons as these can produce pions copiously, and only by disentangling all open channels can one make a statement, about, say, $D\to\pi\pi,K\overline{K}$. In addition, we note that the number of open channels, $N$, counts not only the species in the asymptotic states but also the angular momentum. For example, in the case of $N\gamma\to N\pi$ transitions, the finite-volume matrix element will be contaminated by multiple angular-momentum states, a problem that worsens as the center-of-mass energy is increased.

Given these complications, it is prudent to investigate whether an alternative approach can be used to extract decay and transition rates from lattice QCD. In this work, we consider cases in which total transition rates into all out-states with given QCD quantum numbers are of interest. This is less information than the individual transition amplitudes of Eq. (1), and it is reasonable to ask whether it can be accessed in a more direct manner.

Our method is inspired by Fermi’s Golden Rule. Imagine a quantum system defined initially on a three-dimensional torus, described by a Hamiltonian $H=H_{0}+V$, with an unperturbed part $H_{0}$ under which a particle $\Psi$ is stable (with mass $M$), and a small perturbation $V$ that allows it to decay. According to the standard derivation of Fermi’s Golden Rule, the decay rate of $\Psi$ in its rest frame in infinite volume can be calculated using the double limit222To be precise, the $k$th term in Eq. (2) is equal to the probability that in a measurement done at time $t$ the system is observed in the unperturbed state $k$, divided by $t$.

[TABLE]

where $V_{k}(L)=\langle k,L|V|\Psi\rangle$ is a finite-volume transition matrix element with unit-norm states and where

[TABLE]

can be interpreted as a regularized delta-function. The sum extends over states whose energy $E_{k}(L)$ lies within a fixed interval centered at $M$. A key observation that we exploit in this paper is that, if one is only interested in $\Gamma$, other forms of regularized delta-functions can be inserted into Eq. (2) without affecting the result, as long as the regularization is removed at the end. Note that the order of limits in Eq. (2) is important.

We can rewrite the particle width as

[TABLE]

where the resolution function $\widehat{\delta}_{\Delta}(M,\omega)$ is a regularized delta-function centered at energy $M$ with a characteristic width given by $\Delta$ and falling off sufficiently rapidly at large $\omega$. The expression in square brackets is the finite-volume spectral function, denoted $\rho(\omega,L)$. Related spectral functions are a central concept in finite-temperature lattice QCD studies (see e.g. Meyer (2011b); Aarts et al. (2015); Brandt et al. (2016); Ghiglieri et al. (2016); Burnier and Rothkopf (2013)333In that context, the projector onto the initial state $|\Psi\rangle\langle\Psi|$ is replaced by the canonical thermal average $\frac{1}{Z}\sum_{n}e^{-\beta E_{n}}|n\rangle\langle n|$.). Bringing these ideas together, in this work we present a formalism for extracting total decay and transition rates by extracting smeared spectral functions from appropriately constructed finite-volume, Euclidean correlation functions. The formalism is valid for final-states with any number of outgoing hadrons and does not require disentangling the exclusive transition amplitudes.

Our method requires solving an inverse problem of the Laplace type

[TABLE]

where $G(\tau,L)$ is a Euclidean correlator computed on the lattice. We envision using the Backus-Gilbert method Backus and Gilbert (1968, 1970); Press et al. (2007); Brandt et al. (2015, 2016) to achieve this aim. This gives an estimate for the finite-volume spectral function, smeared by a resolution function that regulates the delta-functions corresponding to the discrete finite-volume energies. The smeared spectral function has precisely the form appearing in the double limit of Eq. (4). By varying the width of the resolution function and the box size of the calculation, one can estimate the double limit required to extract total decay or transition rates.

We emphasize that, in processes where the hadronic final states are produced with varying energy, such as in semileptonic $D$ or $B$ decays, the smeared spectral function obtained from the lattice can be compared model-independently to the experimental decay rate, at the cost of smearing the experimental information in the exact same way. This comparison may be performed in a model independent way with the Backus-Gilbert method, because the smearing kernel (i.e. the resolution function) is known exactly. In particular, no ad hoc functional ansatz is required for the spectral function.

We note that the present application has two distinct advantages over the finite-temperature studies that use the same method. First, the larger range of Euclidean times in zero-temperature calculations should lead to better constraints on the corresponding spectral function. Second, the Backus-Gilbert method is known to work best for slowly varying spectral functions and to struggle in determining narrow features such as resonance peaks. In the present case, we are interested in the spectral function near the energy of the mother particle, but for differing QCD quantum numbers, i.e. those of the final states after the current-mediated transition. Thus, one can hope that in many instances the spectral function will not exhibit a rapid variation, in which case little information is lost due to our lack of energy-resolution.

We close the introduction by highlighting some additional techniques for studying multi-hadron states, beyond the methods related to that of Lellouch and Lüscher. First, Ref. Meißner et al. (2011) describes an idea for extracting resonance parameters directly from Euclidean two-point functions. The approach is applicable to systems with a well-isolated, low-lying narrow resonance and makes use of a model-independent parameterization that allows one to determine resonance parameters by performing a fit to numerical lattice data. Like our approach, the technique of Ref. Meißner et al. (2011) requires estimating the $L\to\infty$ limit.

Second, Ref. Agadjanov et al. (2016) gives a method for extracting the optical potential from a numerical lattice calculation. This allows one to access a particular scattering channel, for example $K\overline{K}$, above the threshold of another channel, for example $\pi\pi$. The finite-volume optical potential is known to contain an infinite tower of poles, and one must fit to a function of this form, apply an $i\epsilon$ prescription, and finally estimate the limit $L\to\infty$ followed by $\epsilon\to 0$. The $i\epsilon$ provides an effective smearing of the optical potential that makes the infinite-volume limit well defined, analogous to the smearing that we achieve via the Backus-Gilbert method. However, the method of Ref. Agadjanov et al. (2016) differs, for example, in that it relies on extracting finite-volume energy levels and requires twisted boundary conditions to sample the finite-volume optical potential.

Third, Refs. Liu and Dong (1994); Liu (2016, 2017) discuss prospects for studying the hadronic tensor by solving the inverse problem on a Euclidean four-point function in the same way as we describe in this work. In contrast to the present study these references do not advocate the Backus-Gilbert technique and do not discuss finite-volume effects, the need for smoothing and the ordered double limit that plays a central role in the present analysis.

Fourth, and finally, Ref. Hashimoto (2017) proposes to use, and performs an exploratory lattice study of Euclidean four-point functions to study semi-leptonic B decays. Again, the main contrast to our work is that the role of the finite volume is not discussed. Ref. Hashimoto (2017) also advocates avoiding the inverse problem by instead integrating the experimental data against a multi-pole function to extract tailored moments that can be more directly compared to lattice data.

The remainder of this paper is organized as follows. In the following section we detail our formalism for estimating widths and differential rates from Euclidean correlators. In Sec. III we describe two specific examples where our approach may be applied, the transition regime between elastic and deep inelastic scattering as well as semi-leptonic heavy-flavor decays. Next, in Sec. IV we discuss the relation of our approach to the Lellouch-Lüscher formalism of Ref. Lellouch and Lüscher (2001). This is followed by a numerical example of the Backus-Gilbert method applied to a toy system in Sec. V. We close with brief conclusions.

II Formalism

In this section we explain our approach for estimating total decay and transition rates from lattice QCD. We also discuss how the technique may be used to study photo-production processes and semi-leptonic decays with differential rates in the photon or lepton-neutrino invariant mass squared.

We begin with a strongly-interacting quantum field theory, described by the hamiltonian density $\mathcal{H}_{\rm QCD}(x)$, and including a stable single-particle state satisfying

[TABLE]

where $E_{N}=\sqrt{M_{N}^{2}+\textbf{P}^{2}}$ and $M_{N}$ is the physical mass of the particle. Here we have in mind a nucleon state, but our formalism holds for any particle that is stable under the strong interaction. In addition to a flavor label, $N$, and total three-momentum, P, we have included $\lambda$ to denote the azimuthal component of the particle’s intrinsic spin.

We next introduce the infinite-volume matrix element

[TABLE]

Here $\mathcal{J}_{\mathcal{Q}}(x)$ is a local current and $\langle E,\textbf{p},\alpha;\mathrm{out}|$ is a multi-hadron out-state with energy $E$, total momentum p, and all other quantum numbers labeled by the combined index $\alpha$. The multi-particle states have standard relativistic normalization. For example, a two-particle state satisfies

[TABLE]

where we have set the collective index $\alpha$ to represent a two-particle state comprised of a pion and a nucleon, with the nucleon carrying momentum k or $\textbf{k}^{\prime}$, and have also defined

[TABLE]

Throughout this work, we denote the number of hadrons in an asymptotic state by $N_{\alpha}$, e.g. in this case $N_{\alpha}=2$.

In Eq. (7) we have allowed the energy and momentum of the final state to differ from the single-hadron initial state. This type of matrix element is appropriate for describing transitions in which photons or leptons, represented by the current $\mathcal{J}_{\mathcal{Q}}$, inject or carry away some amount of energy and momentum. Another case of interest is when the current represents an insertion of the weak hamiltonian, mediating a decay into a purely hadronic final state. In this case we denote the initial state by $|D,\textbf{P}\rangle$ and the operator by $\mathcal{H}_{\mathcal{Q}}$ and define

[TABLE]

Here we have in mind the $D$-meson of the Standard Model, which can decay into many multi-hadron final states. Again the labeling is only suggestive, as the formalism applies to any QCD-stable states.

In lattice QCD it is only possible to calculate finite-volume matrix elements of the form

[TABLE]

Our convention is such that the state $\mathcal{J}_{\mathcal{Q}}(0)|N,\textbf{P},\lambda\rangle$ has quantum numbers $\mathcal{Q}$, and so the final state must have the same quantum numbers, as indicated. We define the finite-volume states with unit normalization throughout, e.g. $\langle E_{k}(L),\textbf{p},\mathcal{Q}|E_{k}(L),\textbf{p},\mathcal{Q}\rangle=1$.

It is non-trivial to relate the finite-volume matrix elements, $M_{k,N(\lambda)\to\mathcal{Q}}(\textbf{p},L)$ and $M_{k,D\to\mathcal{Q}}(L)$, to the transition amplitudes, $\mathcal{A}_{N(\lambda)\to\alpha}(E,\textbf{p})$ and $\mathcal{A}_{D\to\alpha}$. As was discussed in the introduction, and has been argued in various contexts in Refs. Lellouch and Lüscher (2001); h Kim et al. (2005); Christ et al. (2005); Meyer (2011a); Hansen and Sharpe (2012); Briceño and Davoudi (2013); Briceño et al. (2015); Briceño and Hansen (2015), the finite-volume matrix element is equal to a linear combination of all infinite-volume matrix elements with asymptotic final states that carry the same quantum numbers. In the finite volume it is not possible to define asymptotic states, and therefore one cannot isolate exclusive decay amplitudes. For example, if one considers charm to strange decays, $D\to s+X$, then one must include out-states such as $\overline{K}\pi$, $\overline{K}\pi\pi$, $\overline{K}\pi\pi\pi$, $\overline{K}\pi\pi\pi\pi$, $\overline{K}K\overline{K}$.

In this work we present an alternative approach that allows one to directly calculate transition rates that are integrated over all hadronic kinematics, but are differential rates with respect to non-hadronic degrees of freedom. To define such total rates one requires the standard Lorentz-invariant phase-space measure for an $N_{\alpha}$-particle state

[TABLE]

The phase-space measure can be used, for example, to express total decay widths according to

[TABLE]

where integration runs over all real values of all the three-momenta in the measure. At leading order in the weak interaction, $\Gamma_{D\to\mathcal{Q}}$ gives the total width of the $D$-meson into all open hadronic channels with quantum numbers $\mathcal{Q}$. Note that, in the mother particle’s rest frame, this is an alternative version of Eq. (2) in which the infinite-volume limit has been rewritten as a phase-space integral.

Both differential and total rates can be directly extracted from a more general object that we call the transition spectral function and define as

[TABLE]

where the sum runs over all states in the Hilbert space carrying the four-momentum $(E,\textbf{p})$, and where $S_{\alpha}$ is the symmetry factor that avoids double counting phase-space points related by the exchange of identical particles. Here we use the notation appropriate for the nucleon.

To rewrite $\rho_{\mathcal{Q},\textbf{P}}(E,\textbf{p})$ in a more useful form, we now apply a Fourier transform, together with its inverse, to reach

[TABLE]

The expression can be further simplified by moving the second exponential inside the left-most matrix element and identifying it as a standard translation operator on the current

[TABLE]

where the left-most current is now evaluated at $x=(t,\textbf{x})$.

We now identify the integrals over $E^{\prime}$ and $\textbf{p}^{\prime}$ together with the sum over $\alpha$ and integrals over $d\Phi_{\alpha}$, as a sum over all states with quantum numbers $\mathcal{Q}$ in the QCD Hilbert space. This sum over states, together with the outer product appearing between the factors of $\mathcal{J}_{\mathcal{Q}}$, defines an insertion of the identity. Thus Eq. (17) can be rewritten as

[TABLE]

We deduce that the transition spectral function can be written as the expectation value of a product of field operators in a one-particle external state. We emphasize here that the matrix element on the right-hand side is evaluated in infinite volume, with real Minkowski time coordinates and is not time ordered.

In Sec. III, we describe in specific cases how the transition spectral function can be used to compute decay rates and cross-sections. Here we simply note that the total decay width (14) into hadronic final states can be written as

[TABLE]

The main focus of this work, however, is the more interesting case in which the energy of the outgoing hadrons differs from that of the initial state.

At this point, we have established that the transition spectral function, defined in Eq. (15), gives access to differential transition rates and total decay rates, and have also shown how it may be expressed as a matrix element with two current insertions. However, the discussion thus far, summarized by Eqs. (18) and (19), has relied crucially on the fact that all quantities are defined in an infinite volume and with real, Minkowski-signature time coordinates. Thus, the relations do not seem to be of relevance for calculations in lattice QCD, necessarily restricted to a finite volume and to Euclidean signature.

To bridge this gap, we now consider the finite-volume Euclidean correlator most closely related to that used above

[TABLE]

where $\Psi_{\lambda}^{\dagger}(\tau_{i},\textbf{P})$ is an interpolator for the nucleon with total momentum P and spin component $\lambda$, and the subscript “conn” indicates subtraction of $\langle\mathcal{J}^{\dagger}_{\mathcal{Q}}(\tau,\textbf{x})\mathcal{J}_{\mathcal{Q}}(0)\rangle\langle\Psi_{\lambda}(\tau_{f},\textbf{P})\Psi_{\lambda}^{\dagger}(\tau_{i},\textbf{P})\rangle$. Throughout this work we take $\tau>0$. Evaluating the large time limits we reach

[TABLE]

Next we project the current to definite three-momentum, defining

[TABLE]

Finally, inserting a complete set of finite-volume states into Eq. (23) gives

[TABLE]

where the squared magnitudes on the right-side are given by spin averaging the squared magnitudes of the matrix elements defined in Eq. (11).

Equation (24) can be rewritten as

[TABLE]

where

[TABLE]

is the finite-volume spectral function. Substituting this sum of delta functions into Eq. (25) and evaluating the integral immediately gives back Eq. (24). We emphasize at this stage that, while the infinite-volume spectral function gives direct access to the decay width, in the finite volume we have a sum of delta peaks. Naively sampling $\rho_{\mathcal{Q},\textbf{P}}(E,\textbf{p},L)$ at a specific energy cannot give any useful information as the result will either vanish or diverge.

In order to recover the total decay width, we need to construct a sensible infinite-volume limit of the spectral function. To do so we introduce $\widehat{\delta}_{\Delta}(\overline{\omega},\omega)$ as a regularized delta function, centered at $\overline{\omega}$ with width $\Delta$. We require only that this satisfies

[TABLE]

for a smooth test function $\phi(\omega)$. In our approach, $\delta_{\Delta}(\overline{\omega},\omega)$ will tend to zero exponentially for large $\omega$. We then define

[TABLE]

This smoothing procedure replaces the sum over delta functions with a smooth function that has a well-defined infinite-volume limit.

In particular the smoothing allows us to make contact with the discussion of Fermi’s Golden Rule in the introduction. There we described how the width of the mother particle is given by studying the system in finite volume, and summing over the squared magnitudes of individual finite-volume matrix elements, multiplied by regularized delta functions. That construction is identical to that of the smoothed, finite-volume spectral function, defined in Eq. (28). Thus it follows that differential transition rates can be accessed from the lattice framework via the limits

[TABLE]

where the order of limits is important.

The Backus-Gilbert method applied to the inverse problem of determining $\rho_{\mathcal{Q},\textbf{P}}(E,\textbf{p},L)$ from $\widetilde{G}_{\mathcal{Q},\textbf{P}}(\tau,\textbf{p},L)$ leads precisely to smoothed quantities of the form given in Eq. (28). If the correlation function is known at a discrete set of Euclidean times, $\tau_{j}$, then the ‘resolution function’ $\widehat{\delta}_{\Delta}(\bar{\omega},\omega)$ should be constructed from the Laplace kernel,

[TABLE]

with coefficients $C_{j}$ chosen so as to minimize the width

[TABLE]

under the unit-area constraint of Eq. (27). The Backus-Gilbert method then yields an estimate of the smoothed spectral function

[TABLE]

It is well-known that, due to the numerically ill-posed nature of the inverse Laplace transform, reducing the width $\Delta$ is computationally very costly. In the case of the $K\to\pi\pi$ decay, this method is therefore likely to be less accurate than the Lellouch-Lüscher method. Nevertheless, we note that the positivity property of the spectral function ameliorates the inverse problem. To gain some intuition on this, in Fig. 1 we show the accuracy of estimating the infinite-volume spectral function via the smeared finite-volume result, for various values of $1/(M_{\pi}L)$ and $\Delta/M_{\pi}$, in the case of a constant decay amplitude $A_{K\to\pi\pi}$ and negligible interactions of the outgoing pions. In this example we use a normalized Gaussian for the resolution function, $\widehat{\delta}_{\Delta}$.444 More precisely we define

$\widehat{\delta}_{\Delta}(\overline{\omega},\omega)=\frac{1}{\sqrt{2\pi}\Delta}e^{-(\overline{\omega}-\omega)^{2}/(2\Delta^{2})}\,.$

(33)

We caution that the $\Delta$ used here differs from that defined in Eq. (31) by a factor of $4\sqrt{\pi}$.

While not equivalent to the resolution functions that result from the Backus-Gilbert method, these results give an idea of the optimal trajectory to follow in the $[\Delta/M_{\pi},1/(M_{\pi}L)]$ plane. In particular it is manifest that, if one decreases $\Delta/M_{\pi}$ at a fixed $M_{\pi}L$ then the estimator completely fails at some stage, when the value of width becomes too small. This emphasizes the importance of the order in which the limits are to be taken.

While unlikely to be competitive with the Lellouch-Lüscher approach for a single two-particle channel, our method has the advantage of not requiring a detailed understanding of the connection between the finite-volume and infinite-volume matrix elements, which presumably becomes untractable when many channels are open. As already mentioned in the introduction, for differential transition rates the limited energy resolution means that our approach can only yield predictions in relatively broad ‘energy bins’ of width $\Delta$, but it does not imply an uncontrolled systematic error, given that the resolution function (30) is known exactly. In the following section we discuss various specific examples where the formalism presented here seems particularly promising.

III Example applications

In this section we present two examples of phenomenologically interesting cases where our formalism for extracting total transition rates might be applied. In Sec. III.1, we discuss deep inelastic scattering, and in Sec. III.2, semi-leptonic decays of heavy hadrons.

III.1 Deep inelastic scattering

Deep inelastic scattering (DIS), the collision of high-energy leptons with hadrons via virtual photon exchange, has played a major role in particle and nuclear physics. The deep inelastic scattering experiments in the late 1960’s revealed structure within the nucleon with initially surprising scaling laws, leading to partonic models and eventually to the formulation of quark and gluon degrees of freedom in QCD. In this sub-section we describe the application of our general formalism to studying deep inelastic scattering in lattice QCD.

To make the presentation self-contained, we review some of the basics of DIS. In Fig. 2 we give a schematic of the events of interest. A hadron, usually a nucleon, with mass $M$ and four-momentum $p$ collides with a hard lepton carrying momentum $k$ via the exchange of a virtual, space-like photon with momentum $q$. Defining $k^{\prime}$ as the outgoing lepton momentum, we have $q\equiv k-k^{\prime}$. The final hadronic state produced in the collision is not detected, so that we only have information about the total rates into all allowed states with momentum $p_{x}=p+q=p+k-k^{\prime}$. It is further convenient to define $Q^{2}=-q^{2}$, which is positive for space-like photon momenta. We also introduce the Lorentz-invariant parameters

[TABLE]

Following the discussion of Ref. Manohar (1992), we note that the unpolarized cross-section in the nucleon rest-frame

[TABLE]

is given by contracting the hadronic tensor, $W_{\mu\nu}$, discussed in the following paragraphs, with the leptonic tensor, $\ell^{\mu\nu}$, defined as

[TABLE]

where $m_{\ell}$ is the lepton mass and $\Lambda$ a typical scale of the scattering process. Here $u_{s}$ and $\bar{u}_{s}$ are standard spinors, projected to spin state $s$ and normalized so that $\sum_{s}u_{s}(k)\overline{u}_{s}(k)=\not{k}+m_{\ell}\,$, and the $\gamma^{\mu}$ are standard Dirac matrices. Note that the spin-independent part of $\ell^{\mu\nu}$ is symmetric, while the spin-dependent part is anti-symmetric. Thus an unpolarized lepton beam only probes the symmetric part of the hadronic tensor and the polarization asymmetry of cross sections depends only on the hadronic tensor’s anti-symmetric structure. Here we focus on the unpolarized case.

The spin-averaged hadronic tensor, denoted $W_{\mu\nu}$, is defined as

[TABLE]

where $n_{\lambda}=2$ for the nucleon and $j_{\mu}=\sum_{f}Q_{f}\overline{\psi}_{f}\gamma_{\mu}\psi_{f}$ is the electromagnetic current, expressed as a sum over all flavors $f$ of quark bilinears, weighted by their charge, $Q_{f}$. The notation here is slightly different from the previous section with the incoming nucleon carrying three-momentum p. Following Ref. Manohar (1992), we write

[TABLE]

where $F_{1}$ and $F_{2}$ are the structure functions and depend only on $q^{2}$ and $p\cdot q$. Defining $s_{pq}=M^{2}-{(p\cdot q)^{2}}/{q^{2}}\,$, the individual structure functions can be isolated by taking the linear combinations

[TABLE]

The unpolarized cross-section (35) for the inclusive $e+p\to e^{\prime}+{\rm hadrons}$ process can be written in terms of these as

[TABLE]

where $y=(p\cdot q)/(p\cdot k)$ represents the fractional energy loss of the lepton in the nucleon rest frame.

Combining the definition of $W_{\mu\nu}$ with the discussion of the previous section, the task is to invert

[TABLE]

with the correlator

[TABLE]

where $E_{\textbf{p}}=\sqrt{M^{2}+\textbf{p}^{2}}$. The important point is that the role of $\omega$, in terms of which the spectral representation is written, is played by the energy $p_{x}^{0}$ of the final state.

Applying the Backus-Gilbert method leads to a smoothed spectral function

[TABLE]

from which one can estimate the hadronic tensor via

[TABLE]

This expression makes manifest how the lattice calculation yields a sharply defined outgoing total hadron momentum $\textbf{p}_{x}$, but only limited resolution in the corresponding energy $p_{x}^{0}$. The Lorentz-scalar combinations

[TABLE]

are the more primary quantities in our lattice approach than $F_{1}$ and $F_{2}$, because the kinematic factors $s_{pq}$ and $p\cdot q$ (but not $p^{\mu}$) are affected by the limited resolution of $p_{x}^{0}$.

In a lattice calculation, the four variables555Here we ignore possible violations of O(3) rotation symmetry by the cubic box. $\textbf{p}^{2},\textbf{q}^{2},\textbf{p}\cdot\textbf{q}$ and $p_{x}^{0}$ can be varied independently, with the first three sharply defined, up to statistical uncertainties. Due to Lorentz symmetry however, the nucleon structure functions only depend on two invariants: the photon virtuality $Q^{2}$ and the photon energy in the nucleon rest frame, $\nu$, or alternatively for the latter, the Bjorken variable $x$ [see Eq. (34)]. In general, there is therefore a two-parameter family of lattice kinematic variables that realize a given $(Q^{2},x)$ pair. This redundancy can be exploited to minimize the impact of the limited resolution in $p_{x}^{0}$. For instance, in the deep inelastic regime, the structure functions depend only weakly on $Q^{2}$, and therefore one might wish to reduce the impact of the limited resolution in $p_{x}^{0}$ on the Bjorken variable $x$ and accept a broader ‘bin’ in $Q^{2}$.

Explicitly, the connection between the four lattice variables and the relativistic invariant variables reads

[TABLE]

At fixed p and q, $\nu$ is thus an affine function of the dispersion variable $p_{x}^{0}$. Given these relations, it is straightforward to compute the linear propagation of the uncertainty in the variable $p_{x}^{0}$ on the variables $\nu$, $Q^{2}$ and $x$. Eliminating the ‘fuzzy’ variable $p_{x}^{0}$ from Eqs. (49–50), we obtain

[TABLE]

In other words, fixed values of p and q define an ellipse in the $(\nu,\sqrt{Q^{2}})$ plane, centered at $(-\textbf{q}\cdot\textbf{p}/M,0)$ with long radius $E_{\textbf{p}}|\textbf{q}|/M$ along $\nu$ and short radius $|\textbf{q}|$ along $\sqrt{Q^{2}}$. Thus, while the ellipse is ‘certain’ in a given kinematic setup on the lattice, the finite resolution in the variable $p_{x}^{0}$ results in a spread along the ellipse. Contours of constant $x$ in the $(\nu,\sqrt{Q^{2}})$ plane correspond to square-root functions.

Of particular interest is the point of maximum virtuality $Q^{2}$ for given q, which is reached when $p_{x}^{0}=E_{\textbf{p}}$, implying

[TABLE]

At that point, corresponding to the top of the ellipse, to linear order the uncertainty in $p_{x}^{0}$ only affects the variable $\nu$ and not $Q^{2}$; at the same level of approximation, the spread $\Delta x$ in Bjorken $x$ reads

[TABLE]

In Fig. 3, we plot the ellipses (51) for the case $\textbf{p}_{x}=\textbf{p}+\textbf{q}=\textbf{0}$, corresponding to the rest frame of the outgoing hadronic state, along with the $x=\,$constant curves. In this frame, the maximum value of $Q^{2}$ at fixed $\textbf{q}^{2}$ corresponds to $x=1/2$. For a more general frame defined by $\textbf{p}=-\beta\textbf{q}$, the maximum value of $Q^{2}$ corresponds to

[TABLE]

For $\beta=1/2$, the maximum $Q^{2}$ point then corresponds to probing elastic $e,p$ scattering in the Breit frame, $\textbf{p}=-\textbf{q}/2=-\textbf{p}_{x}$. For $\beta=1$, $\textbf{p}_{x}=0$ and one is probing $x=1/2$. To reach ever smaller values of $x$ while keeping the virtuality large, $\beta$ must be increased further and one approaches the ‘infinite-momentum’ frame of the nucleon. However, it becomes increasingly difficult to maintain a useful resolution in $x$, since the relative spread of $p_{x}^{0}$ gets amplified by a factor $\beta$. Hence $x$ and its spread $\Delta x$ become comparable when $\beta={\mathcal{O}}(p_{x}^{0}/\Delta)$, and for a given relative resolution on $p_{x}^{0}$, effectively the smallest $x$ that can be probed is $x={\mathcal{O}}(\Delta/p_{x}^{0})$.

In summary, while there are kinematic limitations and numerical challenges to a lattice calculation of hadronic structure functions, there does not seem to be a conceptual obstacle to probing structure functions down to $x\approx 1/3$ at ${Q^{2}}\approx 4{\rm GeV}^{2}$. It is well-known to be difficult to achieve a good signal-to-noise ratio even in nucleon two-point functions at large separations. This problem is an algorithmic one. To test the ideas presented here it would be interesting to first calculate the structure functions of the pion, for which the signal-to-noise ratio is more favorable, although it also deteriorates as the pion momentum p increases. For that reason, the most promising opportunity for lattice calculations may be to cover the transition region from real photons to the onset of deep inelastic scattering. The formalism presented may also prove useful at the conceptual level, e.g. in rederiving the theoretical predictions for deep inelastic scattering, in particular the development of the amplitude in a series of twist-two operators or its representation through light-like Wilson lines. The formalism is also flexible enough to accomodate the spin dependent structure functions $g_{1}$ and $g_{2}$, different local currents, for instance those determining neutrino DIS or the interference of $\gamma$ and $Z$ exchange, and off-forward nucleon matrix elements.

We also note that our method is complementary to the approach of Ji, presented in Ref. Ji (2013); Chambers et al. (2017). In that work the author shows how structure functions may be studied via equal-time matrix elements in the large momentum limit. Such equal-time matrix elements can be directly extracted from Euclidean correlation functions without solving the inverse problem, but the large momentum limit is challenging in a realistic lattice calculation and the renormalization of the relevant operators is not yet well understood. It would be interesting to compare the two approaches in more detail.

Finally, we remark that the Euclidean correlation function (44), after applying Euclidean-time-ordering and Fourier-transforming with an imaginary frequency ($e^{\omega\tau}$), gives direct access Ji and Jung (2001) to the nucleon forward Compton amplitude for photon virtuality $Q^{2}=\textbf{q}^{2}-\omega^{2}$, a very worthwhile goal in itself.

III.2 Semileptonic decays: $H_{Q}\to\ell+\bar{\nu}_{\ell}+X$

The matrix elements relevant for deep inelastic scattering, discussed in the previous subsection, have a character similar to those that enter semi-leptonic weak decays of heavy mesons. The process of interest is shown in Fig. 4: an incoming hadron, denoted $H_{Q}$ with four-momentum $p_{H}$, decays into a lepton, carrying momentum $p_{\ell}$, an anti-neutrino, with momentum $p_{\nu}$, and a hadronic state. The subscript $Q$ indicates the heavy quark contained within the incoming hadron. We have in mind mesons with $Q=c$ or $b$, but the approach described here can be used to describe any semi-leptonic decay.

Defining $q\equiv p_{\ell}+p_{\nu}$, the hadronic matrix element needed for $H_{Q}\to\ell+\bar{\nu}_{\ell}+X$ is very similar to that described in the previous subsection: the only distinctions are that $q$ flows away from the vertex and is time-like, and that the current mediating the decay is modified. Following Ref. Blok et al. (1994), we aim to calculate the differential decay rate with respect to the lepton energy, $E_{\ell}$, the lepton-neutrino invariant mass, $q^{2}$, and the lepton-neutrino combined energy, $q^{0}$. The final result can be expressed as

[TABLE]

where $V_{qQ}$ is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element describing the flavor change in the decay, $G_{F}$ is the Fermi decay constant and the $w_{i}$ are hadronic structure functions.

The structure functions result from decomposing

[TABLE]

where $\mathcal{J}_{\mu}=\overline{q}\gamma_{\mu}(1-\gamma_{5})Q$ is the flavor-changing current built from one heavy and one light flavor and $M_{H_{Q}}$ is the mass of the mother particle. Here we have expressed the tensor in terms of the lepton-neutrino four-momentum, $q^{\mu}$, together with the four-velocity of the incoming hadron, $v^{\mu}=p^{\mu}_{H}/M_{H_{Q}}$. Given that the tensor only depends on these two four-vectors, one can show that it must satisfy

[TABLE]

This defines the structure functions appearing in Eq. (55).

Thus, the task is to invert

[TABLE]

with the correlator

[TABLE]

where $E_{\textbf{p}}=\sqrt{M_{H_{Q}}^{2}+\textbf{p}^{2}}$. Applying the Backus-Gilbert method leads to a smoothed spectral function from which one can estimate the hadronic tensor via

[TABLE]

We close this section by commenting that our formalism may also be applied to extract total decay widths for purely hadronic decays. One example where total decay widths are of interest is charm to strange decays, in which a weak vector boson mediates a flavor change, $c\to su\bar{d}$. This transition, relevant for Cabibbo-allowed $D$-meson decays, is given at leading order by a weak hamiltonian of the form Buchalla et al. (1996)

[TABLE]

where $\bar{s}(x)$, $c(x)$, $\bar{u}(x)$, $d(x)$ are Dirac fields in the various flavors.

This weak hamiltonian couples the $D$-meson to many decay channels with two or more hadrons. To study this transition via the Lellouch-Lüscher approach, one would first need to generalize the formalism to accommodate all multi-hadron states, including those with four outgoing particles. One would then aim to disentangle the individual decay rates by calculating finite-volume matrix elements in different volumes, all tuned so that the final finite-volume state has the same energy as the incoming $D$-meson. In our approach the total width, $\Gamma_{D\to s+X}$, can be accessed by calculating the appropriate four-point function and then estimating the corresponding spectral function.

A closely related application is the calculation of the lifetime of charmed baryons. Providing estimates of the lifetime of the doubly charmed baryons $\Xi_{cc}^{+}$ Mattson et al. (2002) and $\Xi_{cc}^{++}$ would help their search at LHCb. The spectroscopy of these heavy baryons has recently been studied on the lattice Padmanath et al. (2015).

We comment that operator (61) mixes under renormalization with a second operator in which color indices are contracted between the two bilinears. Thus, the full calculation requires defining $\mathcal{H}_{\mathcal{Q}}$ as a linear combination with appropriate Wilson coefficients. We direct the reader to Ref. Buchalla et al. (1996) for details.

IV Comparison to the Lellouch-Lüscher method

In this section we study the relation between the formalism presented above and the method of Lellouch and Lüscher, described in Ref. Lellouch and Lüscher (2001). For the case of $K\to\pi\pi$ decays, Lellouch and Lüscher derived a relation between finite- and infinite-volume matrix elements

[TABLE]

where

[TABLE]

Here $\phi(q)$ is a known geometric function and $\delta_{\pi\pi}(p)$ is the s-wave $\pi\pi\to\pi\pi$ scattering phase shift due only to QCD. This relation holds up to neglected, exponentially suppressed corrections of the form $e^{-M_{\pi}L}$. In order to extract a physically meaningful decay amplitude directly, one demands that the final two-pion state has the energy of the initial kaon. In the center-of-mass frame this requires tuning the box size, so that one of the energy levels coincides with the kaon mass, $E_{k}(L)=M_{K}$.

This result has since been generalized to states with non-zero total momentum in the finite-volume frame, multiple coupled two-particle channels with non-identical and non-degenerate particles, to particles with intrinsic spin, and to transitions in which the current carries angular-momentum, momentum and energy, so that these quantum numbers differ between the initial and final states h Kim et al. (2005); Christ et al. (2005); Meyer (2011a); Hansen and Sharpe (2012); Briceño and Davoudi (2013); Briceño et al. (2015); Briceño and Hansen (2015).

In the present work we do not require the precise definition of $\phi(q)$. For our purposes it is sufficient to note that

[TABLE]

where $\nu_{k}$ is the degeneracy of the $k$th state in the non-interacting theory, with the counting appropriate to non-identical particles, i.e. $\nu_{0}=1$, $\nu_{1}=6,\cdots$. To see this, we first recall $4\pi^{2}q_{k}^{2}/\phi^{\prime}(q_{k})=\nu_{k}$ where $q_{k}$ is the dimensionless momentum of the $k$th non-interacting state Lellouch and Lüscher (2001).666In the zero-momentum frame the non-interacting energies take the form $E_{k}(L)=2\sqrt{M_{\pi}^{2}+\textbf{n}^{2}(2\pi/L)^{2}}$ corresponding to $q_{k}^{2}=\textbf{n}^{2}$ where n is a three-vector of integers. Equation (64) then follows from the fact that, as $L\to\infty$, any given finite-volume level will coincide with its non-interacting counterpart. In addition, the second, $\delta_{\pi\pi}$-dependent term is suppressed for fixed $k$ and large $L$ and thus vanishes in the infinite-volume limit.

Substituting Eq. (62) into the definition of the finite-volume spectral function, Eq. (26), we find

[TABLE]

and applying the smearing procedure then gives

[TABLE]

Here we restrict our attention to the case where the incoming kaon and the the outgoing two-pion state both have zero spatial momentum.

The procedure outlined in Sec. II dictates that one must now perform the $L\to\infty$ limit followed by $\Delta\to 0$ in order to recover the infinite-volume spectral function and from this the kaon width. We begin with the infinite-volume limit in isolation

[TABLE]

In the first step we used Eq. (64) and then replaced $\sum_{k}\nu_{k}\to\sum_{\textbf{k}}$. We also used that, for a given state, the difference between $E_{k}(L)$ and the corresponding $2\omega_{\textbf{k}}$ must vanish as $L\to\infty$. In the second step we replaced the sum over finite-volume momenta with an integral. This is justified for smooth integrands and in particular relies on the fact that the integrand depends only on smoothed delta functions with finite width $\Delta$.

At this stage we can send $\Delta\to 0$ and see that we recover the total width as is guaranteed by the general formalism presented in the previous section

[TABLE]

where $A_{K\to\pi\pi}$ indicates the physical decay amplitude, in which the outgoing pion pair carries energy $M_{K}$.

We see that, in contrast to the standard Lellouch-Lüscher approach in which one calculates the finite-volume matrix element and converts it to $A_{K\to\pi\pi}$ via Eq. (62), in this approach one removes the factor of $d\delta_{\pi\pi}(p)/dp$ by sending $L\to\infty$. This is only possible given the utility of the Backus-Gilbert method for directly extracting the smeared spectral function, together with the observation that the latter has a well-defined infinite-volume limit. In the case of a single two-particle decay channel, it is unlikely that this approach will be competitive with the Lellouch-Lüscher method. However, as more channels open, our result continues to provide a viable method for estimating total decay widths.

We have performed this same exercise with the generalization of the Lellouch-Lüscher relation to coupled two-particle channels Hansen and Sharpe (2012); Briceño et al. (2015); Briceño and Hansen (2015), and find that the expected result emerges in a similar manner. In the coupled-channel case the squared finite-volume matrix element can be expressed as a vector product of the infinite-volume matrix elements

[TABLE]

where $A_{D\to\alpha}(E)=\langle E,\,\alpha,\,\mathrm{out}|\mathcal{H}_{\mathcal{Q}}(0)|D\rangle$ and $A_{\beta\to D}(E)=\langle D|\mathcal{H}_{\mathcal{Q}}(0)|E,\,\beta,\,\mathrm{in}\rangle$, with $\alpha$ and $\beta$ representing all quantum numbers of the two-particle asymptotic states. Deriving the corresponding smeared spectral function, and applying the infinite-volume limit, we find that the matrix becomes diagonal with each entry given by the degeneracy of the corresponding non-interacting state together with various kinematic factors. One thus perfectly recovers the integral over all three-momenta and the sum over channels, weighted with the proper factors. Sending $\Delta\to 0$ gives the total decay rate, equal to the individual rates summed over all open two-particle channels.

We close this section by noting that one might combine the Lellouch-Lüscher relations with our approach in order to improve the smoothed spectral function. In particular, if the S-matrix is known in the two-particle sector, then one might fit the first few exponentials in $\widetilde{G}_{\mathcal{Q},\textbf{P}}(\tau,\textbf{p},L)$, subtract these states, and replace them with an integrated contribution of the infinite-volume matrix elements. Extracting $\widehat{\rho}_{\mathcal{Q},\textbf{P}}(E,\textbf{p},L)$ from this modified correlation function would give a flatter extrapolation towards $L\to\infty$ and $\Delta\to 0$, improving the precision of the extracted widths and differential rates.

V A numerical test case

This section is devoted to illustrating, by means of a numerical example, to what extent our procedure is able to reproduce infinite-volume total decay widths. The input data will consist of a Euclidean correlation function $\widetilde{G}_{\mathcal{Q},{\textbf{0}}}(\tau_{i},{\bf 0},L)$ of the type presented in Eq. (20), evaluated at $N_{\tau}$ discrete Euclidean time slices up to a maximum extent $L$, at which point we assume that the signal is lost. In this example we ignore finite-temperature effects, i.e. we take the Euclidean temporal direction to have infinite extent.

We study a toy theory with three scalar particles denoted by, $\pi$, $K$ and $\phi$, with physical masses $M_{\pi},M_{K},M_{\phi}$, respectively, satisfying the hierarchy

[TABLE]

and with interactions given by

[TABLE]

Treating these interactions perturbatively, our goal is, given the Euclidean correlator, to reproduce the complete spectral function $\rho_{\mathcal{Q},\textbf{0}}(\omega,\textbf{0})$, which, when evaluated at $\omega=M_{\phi}$, gives access to the total decay width of the $\phi$ particle to leading order in the dimensionless couplings $\lambda$ and $g$. We will neglect final-state interactions altogether.

V.1 Finite Volume

Before turning to the inverse problem of calculating $\widehat{\rho}_{\mathcal{Q},\textbf{0}}(\overline{\omega},\textbf{0},L,\Delta)$, we need to construct the finite-volume Euclidean correlator that will serve as input data. To proceed, note that the finite-volume matrix element [Eq. (12)] between a $\phi$-state and a two-$K$ state projected to $\bf p=\bf 0$ is given by

[TABLE]

where $\textbf{q}_{k}^{2}=k$ starting with $k=0$, and where $\nu_{k}$ is the number of integer vectors q, satisfying $\textbf{q}^{2}=k$. For example $\nu_{0}=1$ and $\nu_{1}=6$. Following Eq. (24) leads to a correlation function of the form

[TABLE]

where we have divided by $M_{\pi}^{3}$ to give a dimensionless quantity. Analogously, one obtains for the three-$\pi$ amplitude

[TABLE]

where $\textbf{n}_{k},\textbf{m}_{k}\in\mathds{Z}^{3}$ are integer vectors representing the $k$th state and

[TABLE]

For example, the $k=0$ ground state is represented by $\textbf{n}=\textbf{m}=0$ corresponding to $E_{0}(L)=3M_{\pi}$ with degeneracy $1$. The next state can be represented with $\textbf{n}=-\textbf{m}$, with $|\textbf{n}|=|\textbf{m}|=1$, and has degeneracy $\nu_{1}=18$. The Euclidean correlator can then be written as

[TABLE]

The full correlator that serves as input for the Backus-Gilbert procedure is then

[TABLE]

V.2 Infinite Volume

Since we want to compare the outcome of the Backus-Gilbert spectral-function reconstruction to the exact infinite-volume result, we briefly derive the infinite-volume contributions to the total decay width for both the two-$K$ and three-$\pi$ channels. Starting with $\phi\to KK$, the decay amplitude is given by

[TABLE]

The main point here is that the amplitude is energy-independent – the mass $M_{\phi}$ appears only because we have chosen to parametrize the interaction by a dimensionless coupling constant $g$. Integration over the two-particle phase space yields a decay width

[TABLE]

Analogously, the three-$\pi$ decay amplitude is energy-independent,

[TABLE]

Integrating over the three-particle phase space yields

[TABLE]

Here $\mathcal{F}(M_{\phi}/M_{\pi})$, shown in the left pannel of Fig. 5, measures the reduction of phase space relative to the case of $M_{\pi}=0$ where the decay products are massless. The definition is

[TABLE]

where

[TABLE]

Note in particular that $\mathcal{F}(3)=0$ and $\mathcal{F}(\infty)=1$. This function has no simple analytic form, but can easily be calculated numerically to arbitrary precision.

Finally, since we want to reconstruct the full spectral function, it is instructive to write the full infinite-volume result for every $\omega$,

[TABLE]

V.3 Inverse problem and smoothing procedure: Backus Gilbert reconstruction

The aim is to recover the full spectral function

[TABLE]

where in particular for $\overline{\omega}=M_{\phi}$ we get an estimate for the total decay width

[TABLE]

For the numerical application, we choose $M_{K}/M_{\pi}=3.55$ and $M_{\phi}/M_{\pi}=7.30$, $g=1$ and $\lambda=10\sqrt{8}$. We use as input points $\tau_{i}=i\cdot a$, with $aM_{\pi}=0.066$ and $1\leq i\leq N_{\tau}$.

Following the discussion of Sec. II, we construct a family of resolution functions $\widehat{\delta}_{\Delta}(\overline{\omega},\omega)$ by finding the optimal coefficients $C_{i}(\overline{\omega},\Delta)$, $(i=1,\dots,N_{\tau})$, that minimize the width subject to the area constraint of Eq. (27). The width is controlled by the number of points, as can be seen from Fig. 5, where we give examples of resolution functions used in the analysis, centered at ${\omega}=M_{\phi}$. The corresponding estimator $\widehat{\rho}_{\mathcal{Q},\textbf{0}}(\overline{\omega},\textbf{0})$ is calculated via Eq. (32). Since Backus-Gilbert is a linear method, an error estimate on $\widehat{\rho}_{\mathcal{Q},\textbf{0}}(\overline{\omega},\textbf{0})$ has to come through the covariance matrix of the input data $\widetilde{G}_{\mathcal{Q}}(\tau_{i},{\bf 0},L)$. We take a realistic covariance matrix, $S_{ij}$, from a pseudoscalar meson two-point function calculated on an $N_{f}=2$ CLS ensemble, that we have used in the past in a similar context. (See Appendix C of Ref. Brandt et al. (2015).) We scale the covariance matrix to give the same relative uncertainty (about $2\%$) on the toy correlator, $\widetilde{G}_{\mathcal{Q}}(\tau_{i},{\bf 0},L)$, as was observed on the actual lattice data. For $N_{\tau}=64$ we take the entire matrix and for $N_{\tau}=32$ we take the block corresponding to source-sink separations, $\tau$, satisfying $32a\leq\tau<64a$.

In order to determine the optimal coefficients, $C_{i}(\overline{\omega},\Delta)$, a poorly conditioned matrix,

[TABLE]

must be inverted (ideally with high precision). The error on $\widehat{\rho}_{\mathcal{Q},\textbf{0}}(\overline{\omega},\textbf{0})$ is kept under control by making the replacement $W_{ij}(\overline{\omega})\to\lambda_{\text{reg}}W_{ij}(\overline{\omega})+(1-\lambda_{\text{reg}})S_{ij}$. In this way the magnitudes of the coefficients $C_{i}(\overline{\omega},\Delta)$, which otherwise exhibit large oscillations as a function of $i$, are tamed, and the statistical error on $\widehat{\rho}_{\mathcal{Q},\textbf{0}}(\overline{\omega},\textbf{0})$ can be kept under control. In the present example, we aim for a precision of $5-10\%$ and find that this is achieved with $\lambda_{\text{reg}}\sim 10^{-4}$. The regulation parameter, $\lambda_{\text{reg}}$, parametrizes a trade-off between resolving power (smaller $\Delta$) and statistical error. For a more detailed explanation on the choice of $\lambda_{\text{reg}}$ see the discussion in Sec. IV E of Ref. Brandt et al. (2015).

The results of our Backus-Gilbert analysis are summarized in Fig. 6. As can be seen from the figure, the agreement of the Backus-Gilbert result with the infinite-volume spectral function is reasonable and improves for increasing $M_{\pi}L$ as expected. If one instead tries to reduce $\Delta$ at constant $M_{\pi}L$, the result becomes unstable since the finite-volume energy levels start to be resolved individually and the result is eventually dominated by finite-volume effects. This illustrates the importance of the order of the limits $L\to\infty$ and $\Delta\to 0$. Nevertheless, the window of reasonable $\Delta$ for a given $L$ seems to be large enough that it may well be useful in realistic numerical lattice applications. As previously mentioned, the resolution function from the Backus-Gilbert method can be used to smoothen the experimental data in the case of differential decay rates. This procedure removes the uncertainty inherent in the solution of the inverse problem at the cost of comparing a somewhat less differential quantity.

VI Conclusions

In this work we have introduced a new method for directly determining hadronic decay widths and transition rates for semi-leptonic scattering and decay processes. The central advantage of our approach is that it can accommodate final states with any number of hadrons.

As we detail in Sec. II, our idea is to construct a Euclidean correlator such that the corresponding spectral function directly gives the decay width or transition rate of interest. We then propose applying the Backus-Gilbert method, which gives an estimator of the finite-volume spectral function, smeared by a known resolution function $\widehat{\delta}_{\Delta}(\overline{\omega},\omega)$ of width $\Delta$. Taking the limit $L\to\infty$ followed by $\Delta\to 0$ then directly gives the experimental observable. As we discuss in Sec. II, and illustrate with a toy example in Sec. V, in a realistic numerical calculation we expect that reducing $1/L$ and $\Delta$ along a suitable trajectory will provide a good estimate of the target infinite-volume unsmeared spectral function.

Although we have focused on the Backus-Gilbert method, we would like to stress that a vast number of different methods exist for solving the inverse Laplace problem. Any of these may be applied in the formalism presented in this work, provided that one can define a width, $\Delta$, that can be varied together with the box size, $L$, to estimate the double limit. Indeed, even within the Backus-Gilbert approach one has freedom to adjust the regulation parameter, $\lambda_{\text{reg}}$, introduced in Sec. V as well as the minimization condition given in Eq. (31). We have applied the Backus-Gilbert method on vacuum correlators in Refs. Brandt et al. (2016, 2015).

The central advantage of the Backus-Gilbert approach is that it offers a model-independent, unbiased estimator of the smeared, finite-volume spectral function with a precisely known resolution function that is independent of the correlator data. For every value of the final-state energy that one aims to extract, one may vary all inputs into the Backus-Gilbert approach to design the optimal resolution function given the quality of the data available. Nonetheless, it may often be the case that the difficulty of reducing the width $\Delta$ is the limiting factor of the calculation. To this end we re-emphasize that one may also smear experimental or model data with the same resolution function to perform a fully controlled comparison.

In certain cases the resolution function may wash away features that one aims to extract. While this is undesirable, the method does faithfully return the information contained in the Euclidean correlator data about the width or rate of interest. Indeed, in the context of constraining scattering amplitudes via Lüscher’s quantization condition, it is in practice well known that a single correlation function provides limited information on finite-volume energies. The modern approach here is to instead construct a matrix of correlators, with a large operator basis, and diagonalize this in order to reliably extract the excited finite-volume states. Analogous methods in the context of finite-volume spectral functions are under investigation Harris et al. (2016).

Acknowledgements.

We thank all our colleagues in the Mainz lattice group for helpful discussions, encouragement, and support. We additionally acknowledge Steve Sharpe, Keh-Fei Liu and Shoji Hashimoto for helpful comments on the first version of this manuscript. D. R. wishes to specially thank Guy Moore for very interesting discussions, and M. T. H. would like to thank Steve Sharpe and Raúl Briceño for insightful discussions and helpful feedback. This work was supported in part by DFG grant ME 3622/2-2 and by the State of Hesse.