New Physics in the Visible Final States of $B\to D^{(*)}\tau\nu$

We derive compact expressions for the helicity amplitudes of the many-body $B \to D^{(*)}(\to DY)\tau(\to X\nu)\nu$ decays, specifically for $X = \ell \nu$ or $\pi$ and $Y = \pi$ or $\gamma$. We include contributions from all ten possible new physics four-Fermi operators with arbitrary couplings. Our results capture interference effects in the full phase space of the visible $\tau$ and $D^*$ decay products which are missed in analyses that treat the $\tau$ or $D^*$ or both as stable. The $\tau$ interference effects are sizable, formally of order $m_\tau/m_B$ for the standard model, and may be of order unity in the presence of new physics. Treating interference correctly is essential when considering kinematic distributions of the $\tau$ or $D^*$ decay products, and when including experimentally unavoidable phase space cuts. Our amplitude-level results also allow for efficient exploration of new physics effects in the fully differential phase space, by enabling experiments to perform such studies on fully simulated Monte Carlo datasets via efficient event reweighing. As an example, we explore a class of new physics interactions that can fit the observed $R(D^{(*)})$ ratios, and show that analyses including more differential kinematic information can provide greater discriminating power for new physics, than single kinematic variables alone.


I. INTRODUCTION
Third, experimentally unavoidable phase space cuts, including both missing mass and lepton momentum cuts used to reduce backgrounds, imply that interference effects between the τ and D * spin states affect all pertinent measurements, including dΓ/dq 2 . The experimental acceptances in the presence of NP may therefore differ from the SM ones used to extract R(D ( * ) ).
To properly capture all these effects, one must compute the matrix elements for the full B → Dτ (→ Xν τ )ν τ and B → D * (→ DY )τ (→ Xν τ )ν τ processes, treating both the τ and D * as internal states. Computations of the corresponding full matrix elements for the SM only have long been available and implemented in prevalently used Monte Carlo generators, such as EvtGen [43,44]. Computations for various parts of the full processes with NP are also available [22,[45][46][47][48][49][50][51], variously omitting the coherent D * decays and interference effects, the τ decays and interference effects, the NP interference effects with the SM, or combinations thereof. In this work, we present a set of generalized NP helicity amplitudes, i.e., matrix elements carrying explicit quantum numbers and full differential phase space dependence, for the full B → D ( * ) (→ DY )τ (→ Xν τ )ν τ processes, in particular for X = ν or π and Y = π or γ. We contemplate NP arising from all possible four-Fermi operators withbcντ flavor structure. We include possible CP violating NP, which may introduce additional large interference effects, and right-handed neutrinos, should they be Dirac. (Some of these operators may also be constrained by other flavor-diagonal and flavor-changing processes in the neutrino sector, but the current limits do not significantly constrain the scale of these operators beyond what is probed in B → D ( * ) τ ν τ .) As such, this paper may be considered as an extension of Ref. [14] to include all the effects mentioned above. for approaches that calculate the matrix element squared directly. A software package implementing these results, for use by experimental collaborations, is under preparation. 1 In Sec. II we establish our notation and conventions. After deriving the amplitudes in Sec. III, we proceed to consider example applications of this efficient computational construction. We construct a MC method in Sec. IV, in which MC data samples are reweighted with matrices of weights. This reweighting need only be performed once per sample, and the result can be used to generate data for any new physics model. Post-reweighting, for any set of NP four-Fermi couplings, the distributions of kinematic observables O i in b i bins can be generated by a smaller set of only i b i linear operations. The general problem of reweighting a large MC dataset between different NP theories is thereby reduced to a much smaller set of linear operations. We use this strategy to efficiently generate 1D and 2D distributions in ten kinematic observables, including lepton and pion energies and opening angles, with and without phase space cuts, over a range of NP couplings. To demonstrate the usefulness of efficiently producing multidimensional distributions, we present a sample bivariate analysis that exhibits higher distinguishing power between SM and NP theories, compared to using only single kinematic distributions.

A. Operator basis
In addition to the SM four-Fermi interaction, we consider a complete set of four-Fermi NP operators mediatingb →cτ + ν τ decay, choosing an operator basis Here we have classified each operator according to the Lorentz structure -scalar, vector, or tensor -of the contracted quark and lepton currents,bΓc andν τ Γτ . The CP conjugate operators for b → cτ −ν τ are obtained by complex conjugation. (We are careful to label the tau neutrino inb →cτ + ν τ distinctly from the tau antineutrino in τ →ν τ X, and from the light lepton flavored neutrino for X = ν . Henceforth we drop all other bars and sign superscripts where the meaning is unambiguous.) We use the convention NP couplings to the quark and lepton currents are denoted by α and β, respectively, normalized to g 2 V cb / √ 2 and g 2 / √ 2, where g 2 is the SU (2) electroweak coupling and V cb is the usual CKM element, while the scale of the operator is normalized to the W mass, m W .
If one views each operator as a tree-level exchange of a fictitious particle, then α and β correspond to its quark and lepton current couplings, respectively, and Λ S,V,T corresponds to the mediator mass. The NP couplings may be complex in general, admitting multiple sources of CP violation. We label the chirality of the leptonic β couplings according to the tau neutrino chirality, in order to easily distinguish between contributions involving leftand right-handed neutrinos, and hence contributions that do or do not interfere with the SM operator. Neglecting neutrino masses, β L and β R terms do not interfere. The chirality of the quark couplings α L,R are defined by the chirality of the charm quark. The identity with 0123 = +1, 2 guarantees the absence of α T L β T L or α T R β T R terms, so that there are only two tensor operators. This yields a total of ten independent four-Fermi NP operators. Neutrino flavor-violating effects are GIM-suppressed and may be neglected. Finally, we assume in this paper that τ decays are described by the SM, supported by the good agreement of SM predictions with τ decay data [53]. 2 Our sign conventions imply that Tr[γ µ γ ν γ σ γ ρ γ 5 ] = −4i µνρσ . Fixing instead sign conventions such that Tr[γ µ γ ν γ σ γ ρ γ 5 ] = +4i µνρσ , as done in many places in the literature, changes the sign of eq. (2.2), as well as the sign of g(q 2 ) and a T±,0 (q 2 ) in eqs. (2.8).

B. Form factors
Lorentz symmetry ensures that for the B → D ( * ) transitions, the scalar, pseudoscalar, vector, axial vector and tensor currents have one (zero), zero (one), two (one), zero (three) and one (three) independent form factors, respectively. We define so that q 2 is the only unfixed Lorentz invariant in the B → D ( * ) decay. Note m 2 τ ≤ q 2 ≤ (m B − m D ( * ) ) 2 , and that q µ is equivalently the momentum flowing to the τ ν τ pair. For B → D we adopt the following conventions and definitions for the form factors, The pseudoscalar and axial vector currents D|cγ 5 b |B ≡ 0 and D|cγ µ γ 5 b |B ≡ 0, while the axial tensor current D|cσ µν γ 5 b |B is fixed by the identity (2.2). Under these conventions, at leading order in Λ QCD /m b,c , these form factors are where ξ(w) is the Isgur-Wise function [7,8]. These relations are understood for the value of ). Under CP conjugation, the form factors for the conjugate B →D process are Similarly forB → D * we define The matrix element of the scalar current vanishes, D * |c b |B ≡ 0, while the axial tensor current matrix element D * |cσ µν γ 5 b |B is fixed by the identity (2.2). At leading order in Under CP conjugation, the form factors for the conjugate B →D * process are noting that the pseudoscalar and axial currents do not change sign.
Helicity angle definitions with respect to spatial momenta (bold symbols) in the sequence of particle rest frames. Each subfigure is drawn in the rest frame of the particle denoted in the central grey disk. Transformations between frames are achieved by Euler rotations and Lorentz boosts, denoted by gray arrows (see text for details).

C. Helicity angles
The helicity amplitudes are most simply expressed in terms of the (θ, φ) helicity angles for each vertex of the B →D ( * ) (→DY )τ + (→ Xν τ )ν τ amplitude. 3 That is, we factorize the phase space of the process into a series of rest frames in the (off-shell) cascade B → D ( * ) (→ DY )W(→ ν τ τ (→ν τ W (→ X))) and so on. Here, for the purpose of defining helicity angles, we treat the τ ν τ pair as originating from a fictitious W particle in the B → D ( * ) transition, with momentum q µ . Similarly we define p µ to be the momentum of the W * in the τ decay, and p 2 ∈ [0, m 2 τ ] neglecting the daughter charged lepton's mass. (Hereafter we always label the momenta of massive particles with the base symbol p and those of massless particles with the base symbol k.) In Fig. 1 we show schematically the helicity angle definitions for B → D ( * ) (→ DY )W(→ ν τ τ (→ν τ W (→ ν ))), with Y = π or γ. Explicit expressions for these helicity angles in terms of Lorentz invariant objects are provided in Appendix A. The polar θ angles in Fig. 1 are well-defined rest frame by rest frame. The orientation of the azimuthal φ angles is, however, defined with respect to an arbitrary direction in the B rest frame, (θ * , φ * ), combined with a sequence of parent-daughter frame transformations. As the B is a spin-0 state, the (θ * , φ * ) angles themselves are unphysical, and vanish from all amplitudes, but we nonetheless keep these angles explicit in Fig. 1. In a parent rest frame with daughter polar coordinates (θ, φ), the parent-daughter frame transformation is defined to be the sequential

D. Phase space
The phase space integration limits are [0, π) and [0, 2π) for each polar and azimuthal helicity angle. In these coordinates, the full phase space measure can be straightforwardly factorized into B → D ( * ) τ ν τ , τ → Xν τ and D * → Dπ, Dγ pieces. These are Here the spatial momentum of the τ ν τ pair in the B rest frame and of the pion in the D * rest frame are, respectively, the usual phase space factor.

III. AMPLITUDES
The helicity amplitudes for the full B → D ( * ) (→ DY )τ (→ Xν τ )ν τ process carry only quantum numbers of external particles (i.e., not the τ and D * spins) corresponding to certain convenient basis choices for external spinors and polarization vectors. For X = ν , these are the spins s ντ , sν τ , s , s ν = −, + that label the helicity amplitudes below, and also the photon helicity κ = ± in the case of D * → Dγ.
The azimuthal helicity angles arise as phases in the helicity amplitudes. These phases are odd under CP , along with those that occur in the NP α or β couplings. In the remainder of this paper, we shall consider explicit expressions for only theb →c process. Results for the CP conjugate b → c process are obtained by conjugation of all these phases, i.e., where s is the set of quantum numbers of all external states, ands the corresponding CP conjugate, obtained by interchanging all spins and helicities with their conjugates.
Since we assume that τ decays are described by the SM, and we can neglect the light charged daughter lepton mass, it is always the case that sν τ = +, s = +, and s ν = −, such that our choice of spinor basis for massless states coincides with the usual helicity basis. We drop these quantum numbers from the amplitude labelling below, with the understanding that all other amplitudes are zero. For the SM, s ντ = − only. However, in the presence of NP currents involving left-(right-)handed ν τ , associated with β L (β R ) couplings, one may further have s ντ = − (s ντ = +) contributions that do (do not) interfere with the SM.

A. Amplitude factorization and τ spinor basis
It is convenient to express the helicity amplitudes factorized into B → D ( * ) (→ DY )τ ν τ and τ → Xν τ pieces, not only for the sake of presentation, but also in order to enable the B → D ( * ) (→ DY )τ ν τ results to be used modularly with respect to different choices of τ → Xν τ . To obtain the square of the polarized matrix elements, one sums over the internal τ spin, s τ = 1, 2, 4 before squaring, and similarly for B → D * (→ DY )τ (→ Xν τ )ν τ . Here s X (s Y ) is the set of quantum numbers of the X (Y ) external state: For a massive fermion, we label spin states by 1 and 2 (see, e.g., p.48 in Ref. [54]).
empty for X = π (Y = π). The fully differential decay rates are then where we have included the factorized phase space measures (2.10) as well as τ and D * propagators, using the narrow width approximation for both states.
In order to permit extension of the results below to any τ → Xν τ decay, we specify here our choice for the τ spinor basis and phase conventions. Calculation of the helicity amplitudes are achieved by decomposing momenta and spinors (or polarizations) of massive states onto a lightcone basis. For the τ , we choose the ν τ momentum k ντ as a null reference momentum. In the τ rest frame, using phase space coordinates as defined in Fig. 1, the Dirac spinor basis for the τ + is This additional phase factor in the τ + spinors is balanced by a cancelling phase factor e ±iφτ in the corresponding B → D ( * ) (→ DY )τ ν τ amplitudes. We emphasize that this is merely a bookkeeping device, that does not affect the physical phase structure of the full B → D * (→ DY )τ (→ Xν τ )ν τ helicity amplitudes. Under this phase convention the τ → Xν τ helicity amplitudes therefore carry s ντ as a quantum number, even though ν τ itself is not involved in the τ decay.
The quantum numbers in eq. (3.5) need only be matched with those in each of the B → D ( * ) (→ DY )τ ν τ helicity amplitudes below to identify the corresponding τ spinor and phase to be used to compute the τ decay helicity amplitude of interest. We provide below explicit expressions for the τ → ν ν τ and τ → πν τ amplitudes under these conventions.
Let us now proceed to present the helicity amplitudes. For readability, we group terms by form factors. For B → Dτ ν τ , the helicity amplitudes Expressions for the SM helicity amplitudes may be read off taking all α's or all β's to zero. These SM results numerically match the output of EvtGen. In the SM, only A − 1 and A − 2 are non-zero, and contain terms that are all linear or zeroth order in m τ , respectively. Interference effects arising from decay of the s τ = 1, 2 spin states to the same final state therefore enter at O(m τ /m B ) in the SM. When treating the τ as stable, interference terms for operators that respectively couple toν τ L τ L andν τ L τ R , such as the f S f + term between the NP scalar and SM vector operators within A − 1 , are chirally suppressed as expected, entering only at order m τ /m B . However, interference between τ spin states can produce O(1) contributions to these terms, e.g. the f S f + interference term between A − 1 and A − 2 . Similar conclusions follow for B → D * τ ν τ , below.
The decay D * → Dπ proceeds through the operatorĝ π D * µ (π∂ µ D − D∂ µ π), in whichĝ π is the phenomenological couplinĝ withĝ π = (m D * /f π )g π in the notation of Ref. [9]. We define the functions where again r V,S,T ≡ m W /Λ V,S,T . Expressions for the SM helicity amplitudes may be read off taking all α's or all β's to zero. These SM results numerically match the output of EvtGen.
Note that orthogonality of the ∆ and Σ functions permit us to read off from the amplitudes which square and cross-terms contribute under integration over full angular phase space, and which are absent. For instance, the f (q 2 ) g(q 2 ) cross-term integrates to zero.
However, in the presence of angular phase space cuts, such terms do contribute. D * interference terms correspond to cross-terms within or between the ∆ or Σ functions that contain orthogonal θ D or φ D dependence, and are typically O(1).
D. τ → ν ν τ and τ → πν τ Under the conventions of eq. (3.5), the helicity amplitudes Note the quantum number, s ντ , belonging to the τ neutrino in the parent B → D ( * ) τ ν τ process, is a consequence of our spinor phase conventions in eq. (3.5), which ensures that φ τ appears only in the physical combination φ τ − φ W .
For τ → πν τ , we adopt definitions for the helicity angles by replacing the W with a pion in the τ decay within Fig. 1, and replacing (θ W , φ W ) → (θ π , φ π ) and p µ → p µ π . The helicity amplitudes Here f π = 93 MeV is the pion decay constant.

IV. APPLICATIONS
The computation of the NP helicity amplitudes for B → D ( * ) (→ DY )τ (→ Xν τ )ν τ decays permits us to efficiently reweigh large Monte Carlo samples to any theory generated by the NP operators (2.1). We may thereby access the full kinematic structure of the (visible) τ and D * decay products, and explore the NP effects therein. To illustrate the potential usefulness and NP discrimination power of these results, in this section we provide a first exploration of such NP effects for B → D * (→ Dπ)τ (→ ν ν τ )ν τ , focusing on NP scenarios compatible with the B → D ( * ) τ ν τ rate [26]. We include effects of q 2 , missing momentum, and lepton energy cuts in this analysis. However, background modelling, detector simulations, or B → Dτ ν τ pollution, all of which are required for a realistic analysis, are deferred to a future study [52].

A. Monte Carlo strategy
In accordance with the results of Sec. III, the full B → D * (→ Dπ)τ (→ ν ν τ )ν τ helicity amplitudes may be expressed in the linear form where M v is a vector of amplitudes and the 11-dimensional vector v is The first entry of M v corresponds to the SM contribution. By construction, M v is independent of the particular NP model, but depends only on phase space configuration. Our We shall consider here an MC sample of 10 million events, reweighted once on the full phase space, and once with application of the phase space cuts, motivated by Refs. [2,3], With three neutrinos in the final state, the remaining visible phase space for B → D * (→ Dπ)τ (→ ν ν τ )ν τ is parametrized by seven independent parameters. In the B rest frame we compute an overcomplete set of ten observables, including where cos θ XY is the opening angle between p X and p Y , as well as the normalized triple product and the missing invariant mass, respectively, To generate the B → D * τ ν τ form factors (2.9), we use the ISGW2 parametrization [55,56] for f (q 2 ) as presently implemented in EvtGen [43,44] and obtain the q 2 -dependence of the rest via the leading order HQET relations (2.8).

B. Univariate versus bivariate analyses
Various NP scenarios may produce B → D ( * ) τ ν τ rates commensurate with the central values of current observations. In particular, leptoquark models with couplings In this section, as an example, we focus on the NP model with g T ≡ α T R β T L r 2 T = −0.38. In Figs. 2 and 3, we present the differential distributions for each of the ten kinematic observables (4.4)-(4.5) in the full and cut phase space, respectively, generated by ranging over g T ∈ [−0.76, 0], i.e., over a range spanning twice the best fit g T value. We also show the distributions for g T = −0.38 and the SM. While the q 2 distribution itself has some discriminating power between the SM and the NP along the g T contour, other observables, in particular E , cos θ D , and cos θ π may be just as, if not more, discriminating.
To explore this further, in Fig. 4 we present density plots of the doubly differential decay rates with respect to three pairs of kinematic observables, for the SM (top row), g T = −0.38 (middle row), and their difference (bottom row). In particular, the density plots for the difference of d 2 Γ/dq 2 dE and d 2 Γ/dq 2 d cos θ π have non-trivial level contours, suggesting that an analysis using both of these observables may have significantly more SM-NP discrimination power than q 2 or any other single kinematic observable. (A preliminary multivariate study of all ten observables with a boosted decision tree trained to discriminate the SM and the g T = −0.38 model supports this claim [52].) To roughly quantify the relative discrimination power of single and doubly differential distributions in the q 2 -E space, we proceed to divide the MC sample into two bins -a "2binning" -according to a partitioning in each of the one-dimensional q 2 and E distributions   For each 2-binning, we define a discriminator, where n 1,2 are the two bin entries, T (H) labels the true (hypothesis) theory, and σ 2 is a 2 × 2 covariance matrix. An approximate covariance matrix for the three 2-binnings is constructed based on the distributions presented in Ref. [3], measured in a signal-rich region approximated by the phase space cuts (4.3). We decompose the covariance matrix as where we have suppressed the indices. The first term, σ 2 data , corresponds to the Poisson error of the measured data in each bin, while σ 2 bg corresponds to the error in the normalizations of the main background components, mainly the D * * backgrounds, which are fixed by data in different kinematic regions. Both terms therefore scale with the square root of the luminosity.
Rescaling statistics to a initial benchmark luminosity of 5 ab −1 at Belle II implies σ data 10% and σ bg 14%. While σ data is uncorrelated by construction, we assume σ bg is purely an error in overall normalization, and therefore fully correlated between the two bins. By looking at the systematic error breakdown in Ref. [3], we divide the systematic components into a fully correlated systematic error σ sys and a component σ shape coming from D * * background shape variations of unknown correlation between the two bins. We conservatively assume that systematic errors remain the same in the future, therefore setting σ sys ∼ 4% and σ shape ∼ 3%.
We emphasize that translation of the χ 2 values, obtained from this approximate covariance matrix (4.10), into statistical confidence levels requires a more comprehensive treatment of backgrounds and their correlations than attempted here, beyond the scope of the present work. However, the relative size of χ 2 values for different 2-binnings is less sensitive to background correlation effects, and therefore can be thought of as a proxy for the ratio of the actual χ 2 statistics.
As an example, we now suppose either the SM or the g T = −0.38 model to be the true theory, and consider the space of hypotheses g T = [−0.76, 0.76]. In Fig. 5 we show corresponding χ 2 bands for both theories, generated by ranging over arbitrary correlation for σ shape , with phase space cuts (4.3). We see in Fig. 5 that the two-dimensional 2-binning for the SM (g T = −0.38) true theory excludes the g T = −0.38 (SM) hypothesis with greater confidence than either of the single observable 2-binnings alone. However, for g T hypothesis ranges closer to the true theory values, the lepton energy E 2-binning has greater distinguishing power. An optimized discrimination of these theories using a multivariate analysis will be studied elsewhere. As an example, we have presented a preliminary exploration of kinematical effects in the phase space of B → D ( * ) (→ Dπ)τ (→ ν ν τ )ν τ for a class of theories with a NP antisymmetric tensor current. Our amplitude-level calculation makes it feasible to efficiently compute an event 'weight matrix' in the space of NP couplings, so that reweighting of the MC dataset need be performed only once per data sample. In this way, not only single but also multidimensional distributions can be rapidly computed for any NP theory. We find that bivariate analyses can exhibit greater discriminating power of the SM versus NP models.
Directions for future study include computing the analogous helicity amplitudes for B → D * * τ ν τ using recent form factor results [57], in order to examine the interference effects from the τ and D * * decays. One might also extend the bivariate analysis to consider the hadronic τ → πν mode, given recent results using single kinematic variables [5]. Employment of a boosted decision tree to perform a complete multivariate analysis of the full phase space is also planned. A comprehensive treatment of backgrounds and detector effects will permit estimation of the corresponding statistical confidence levels and future NP exclusion limits achievable with such multivariate analyses at current and upcoming experiments. A software package, Hammer [52], is under development, which can be incorporated into existing software pipelines that account for these background and detector effects.
and for B → D * τ ν τ processes Note that cos θ W is defined with p dependence implicit, so that for τ → πν τ one need only replace θ W → θ π in eq. (A1). In these expressions, the B rest frame energies and the D * rest frame energy E * π = (m 2 D * − m 2 D + m 2 π )/2m D * . For the azimuthal angles, only the combinations φ τ − φ W , φ W − φ and φ D − φ τ appear in the helicity amplitudes. We therefore provide direct expressions for the sine and cosine of these relative twist angles, rather than for the azimuthal helicity angles themselves. To keep expressions short, we express these twist angles iteratively in terms of trigonometric functions of the polar helicity angles, with 0123 = +1, and for B → D * τ ν τ processes Hence, one need only explicitly express half of the helicity amplitudes.
The decay D * → Dγ proceeds via the operator (eµ a /4) µνρσ (∂ µ D * ν − ∂ ν D * µ )F ρσ D, in which, following the notation of Ref. [9], µ a is a magnetic moment such that We define the functions The Ω and Ξ functions play the same role as ∆ and Σ in the D * → Dπ mode above. That is, the s τ = 2 (s τ = 1) helicity amplitudes are linear combinations of the Ω (Ξ) functions exclusively. Each set of Ω and Ξ functions is L 2 (C) orthogonal under integration over the angular phase space dΩ D dΩ τ , while the Ω functions are orthogonal with respect to Σ with the inclusion of an additional e ±iφτ phase in the integration measure, in accordance with our τ spinor phase conventions (3.5). One finds with r V,S,T ≡ m W /Λ V,S,T . The four remaining helicity amplitudes [A γ ] −− sτ and [A γ ] ++ sτ follow immediately from the parity relation (B1).