Angular analysis of new physics operators in polarized τ → 3 ` decays

: In a bottom-up approach we investigate lepton-ﬂavour violating processes τ → 3 ` that are mediated by New Physics encoded in eﬀective-theory operators of dimension six. While the opportunity to scrutinize the underlying operator structure has been investigated before, we explore the beneﬁts of utilising the polarization direction of the initial τ lepton and the angular distribution of the decay. Given the rarity of these events (if observed at all), we focus on integrated observables rather than spectra, such as partial rates and asymmetries. In an eﬀort to estimate the number of events required to extract the coupling coeﬃcients to the eﬀective operators we perform a phenomenological study with virtual experiments.


Introduction
Within the Standard Model (SM) of particle physics, lepton flavour is conserved as long as the neutrino masses are exactly zero. The discovery of massive neutrinos and neutrino oscillations, however, shows that lepton flavour violation (LFV) is in principle allowed. For example the process τ − → µ − γ can occur at the one-loop level if a tau-neutrino oscillates into a muon-neutrino within the loop that is supplemented with a charged vector boson. If this mechanism were the only source of LFV then the branching fractions of the τ − → µ − + − and τ − → e − + − decay channels are non-zero, but tiny -O(10 −45±5 ) or thereabouts -and clearly unobservable. Many New Physics (NP) scenarios, however, predict much higher branching fractions, very roughly O(10 −10 ), which are at the edge of observability (see e.g. the review in [1]). Here and below denotes either a muon or an electron.

JHEP10(2015)082
If such a decay were to be discovered (see e.g. [2]), as is not unrealistic given the hints on deviations from the SM in the lepton sector with the recent measurements 1 of R K = B(B → Kµ + µ − )/B(B → Ke + e − ) by the LHCb collaboration [3] and the LFV Higgs decay B(h → µτ ) = (0.89 +0. 40 −0.37 )% by the CMS experiment [4], it will be highly interesting to investigate the underlying interaction structure in order to disentangle possible NP models. The current limits on the branching fraction are around Br(τ → 3 ) 2 · 10 −8 [5], and the production cross section for tau pairs in e + e − collisions is σ(e + e − → τ + τ − ) = 0.919 nb [6]. If one assumes a branching fraction just below the current limits, say 10 −9 or so, these numbers translate to roughly one ab −1 of integrated luminosity per τ → 3 event. Even with a reduction in efficiency due to our polarization requirement one could expect a few of such events within 50 ab −1 of collected data, which is the stated goal of the world's flavour factories [7,8].
In this paper we employ a bottom-up approach and treat the SM as an effective field theory, i.e. we consider higher-dimensional operators that consist only of SM fields and respect SM symmetries. Our aim is to gain qualitative and quantitative information on the couplings of these operators. It is obvious that such an approach will allow us to only gain direct insight into the couplings at the low scale -set by the tau-lepton massand not at the high scale of the responsible NP mechanisms. The implementation of the renormalization-group running of the couplings, which has been derived at one-loop level in [9][10][11], is beyond the scope of the current paper.
As some of us have shown in a previous publication [12], the τ → 3 decays are mediated by a handful of dimension-6 operators with possibly complex coefficients. In that publication a Dalitz-plot analysis was entertained in order to differentiate between radiative and leptonic operators. Clearly, a reconstruction of a Dalitz distribution requires a large data sample which obviously is hard to obtain for a very rare, lepton number violating decay. Thus, in order to obtain information on the type of operator that mediates the decay and on the helicity structure of the interaction, appropriate observables need to be defined. Obvious candidates are observables such as forward-backward asymmetries, which can be measured even at very low statistics. Our strategy is therefore to define observables of partially integrated phase space, which can be measured by a simple counting experiment.
Additional information on the structure of the interaction can be gained by studying the decay of polarized tau leptons [13]. Such a polarization can be realized at e + e − colliders [14] running close to the τ + τ − threshold with polarized electrons or electrons and positrons. As we shall show in this paper, taking into account the spin direction of the decaying tau lepton allows us to obtain information on the structure of the interaction, even with quite small data samples.
The idea of using polarized tau leptons has also been discussed in [15], however, in the context of a specific NP model. Our focus is a general analysis in order to pin down the structure of the relevant interaction.

JHEP10(2015)082
Considering the possible LFV tau decays into leptons one may classify the six distinct channels as (a) τ − → − − + with = µ, e, or (b) τ − → µ − e − e + , τ − → e − µ − µ + , or (c) τ − → e − e − µ + , τ − → µ − µ − e + . The processes in (a) and (c) involve identical particles in the final state, whereas (b) does not. Radiative operators contribute to classes (a) and (b), but not (c). For definiteness and simplicity we focus our attention solely on case (a), in which two external mass scales suffice, the mass of the tau lepton, m τ , and the mass of the lighter lepton, m.
The paper is organised as follows: in the next section we define the effective Hamiltonian with open coefficients ξ i , where i counts through the various operators of massdimension 6. The setup of our calculation is described in detail, including the polarisation vector of the initial tau and the angles characterizing the position of the polarisation vector relative to the decay plane. As a result the totally differential decay rate is decomposed in terms of trigonometric functions of these angles. In section 3 we integrate the differential decay rate over the relevant parts of the phase space. The results are rather bulky in print, so we have diverted them to the appendix for easier reading. Section 4 consists of a phenomenological study of the two decay channels τ → 3µ and τ → 3e. The main question we try to answer is how well one could determine the Wilson coefficients of the dimension-6 operators under the hypothesis of having discovered a low number of events experimentally. We summarize in section 5.

Operator basis
In the following, we adopt the conventions and notation of [12]. (For more details we also refer the reader to that reference.). The starting point is the most general set of dimensionsix operators respecting the SM gauge symmetries (see [16,17]). After integrating out the weak gauge bosons and the Higgs field after electroweak symmetry breaking, there remain four purely leptonic 4-fermion operators of dimension six. Ordering the operators by the chirality of the involved lepton fields, they can be expressed as . Notice that in order to feed terms of the form (LR)(LR), one would have to include dimensioneight operators in the SM effective field theory. Following [12], we assume that the NP scale Λ is sufficiently large compared to the electroweak scale, so that these contributions can be neglected (for similar reasoning in other context, see e.g. [18,19]). In an explicit UV completion of the SM, the dimensionless couplings g V should be determined from a matching calculation.

JHEP10(2015)082
Similarly, one generates radiative operators, where, at low energies, only couplings to the photon field have to be considered. (Contributions from intermediate W ± or Z 0 bosons are already contained in (2.1), see e.g. the matching calculation in [20].) We are then left with two more terms in the effective Hamiltonian, They contribute to the τ − → − − + amplitude via photon exchange, schematically ¯ |H Here, for simplicity, we used the same notation for lepton fields and on-shell spinors, q µ denotes the momentum flow through the virtual photon, and α em the usual electromagnetic fine-structure constant. In summary, the generic effective Hamiltonian for τ → 3 decays can be written as For future convenience we shall combine the six couplings into a complex-valued vector of mass dimension (-2).

Spin polarization of the tau lepton
The spin polarization of a beam of particles is typically measured in the flight direction of said particles. In the setup of our calculation we will define the z-axis of our lab-frame coordinate system to be the flight direction of the tau lepton. The reference vector for the spin orientation is then chosen as s µ = (0, 0, 0, 1) in the tau lepton's rest frame, with s 2 = −1 and s·p τ = 0. In the calculation of the squared amplitude |M| 2 we will then use the Dirac matrix to project onto tau leptons with "spin-up". In this way the spin vector s µ only appears at most linearly in the squared amplitude. Note that there are then only three linearly independent invariants that one can build from s µ and the lepton momenta in the decay . These can be taken as

JHEP10(2015)082
and therefore the squared amplitude of the spin-up tau decay can be decomposed as Similarly the decay of a "spin-down" polarized tau lepton can be calculated with the help of the projector u ↓ū↓ = ( / p τ +m τ )(1−γ 5 / s)/2. We stress that the direction in which the spin is measured, i.e. the z-axis (or reference vector s µ ) remains fixed. With this point of view 2 the variables defined in (2.6) remain unchanged when describing the decay of a spin-down rather than spin-up tau lepton. In this spin-down case one finds the same functions g 1,t,u,v from above that describe the squared amplitude, albeit in the combination (2.8) Therefore only the part g 1 contributes to the unpolarized decay rate. The information contained in the functions g t,u,v would then be lost and could not be used in the effort to unveil the underlying operator structure. In what follows, we will use a slightly modified version of the decomposition (2.7) based on angles associated with these invariants.

Euler rotations
The three momenta p 1 , p 2 , p 3 of the decay product span the decay plane. In general the spin direction s does not lie in this decay plane, and its orientation can be described with the help of two angles. For example, the angle between the normal of the decay plane and s is the first, and the angle between the momentum p 3 of the antilepton and the projection of s into the decay plane the second angle. For convenience we employ the technique of Euler rotations instead, which also serve to define the orientation of s and the decay plane: we define a reference frame in which the tau lepton is at rest and the z-axis is aligned with the spin vector s = e z . The aim is to then rotate the system into a convenient position, which can either be achieved by rotation of the coordinate system e i while keeping the momenta and spin direction fixed, or by keeping the coordinate system fixed and rotating the physical vectors. Our choice in this paper is the latter. We will hence perform rotations on the vectors p i and s simultaneously until the decay plane coincides with the x − z plane and the antilepton's momentum points in the z direction. This procedure is illustrated and further explained in figure 1. The vectors after rotation are endowed by a prime, to wit with mathematically positive convention: (2.10) 2 Alternatively one may perform the transformation s µ → −s µ to find the differential decay rate of a "spin-down" polarized tau.

JHEP10(2015)082
Hence the polarization direction vector reads s = (cos α sin β, sin α sin β, cos β) . (2.11) Note that there is no dependence on the first rotation angle γ, which reflects the azimuthal symmetry of the problem. Since the rotations are such that the decay plane coincides with the x − z plane we may parameterize where the parameter ϕ = φ 23 has a physical interpretation of the angle between p 2 and p 3 within the decay plane, which is independent of the rotation angles α, β, γ. The lengths of the 3-vectors b 2 = | p 2 | = | p 2 | and similarly b 3 are expressed in terms of the energies E i as While the momentum of the antilepton, p 3 , is well defined, the lepton pair in the final state is indistinguishable if the leptons are of the same flavour, and we must be careful not to double count. For the last Euler rotation, it is sufficient to consider values of α ∈ [0, π). Therefore, in total we restrict the values of the rotation angles as Furthermore it is obvious that either p 1 has a positive x component and p 2 a negative one, or the other way around. The parameter φ 23 is well defined as the angle between p 2 and p 3 if we insist on naming that lepton's momentum p 2 which has a positive 3 x component, i.e. φ 23 ∈ [0, π) .
(2.15) (However, the squared amplitudes, which we are going to calculate, will be symmetric under the exchange of p 1 ↔ p 2 , as they must for identical particles in the final state.) In terms of the angles and energies the invariants in (2.6) are Therefore we may decompose the squared amplitude in (2.7) in terms of the trigonometric functions that appear above instead of the invariants themselves. Specifically we have (2.17) This setup allows us to quite easily pick out the individual contributions from J 1 to J 4 by folding the differential decay rate with appropriate weight functions. 3 With this convention the invariant v = ε αβγδ p α 2 p β 3 p γ τ s δ = ε αβγδ p α 1 p β 2 p γ 3 s δ is negative definite. configuration as seen in one of the rest frames of the tau, in which the spin-direction s is along the z axis. The final-state momenta and spin vector are then rotated by γ around the z-axis. This is the angle between the antilepton's momentum p 3 , projected into the x − y plane, and the (−x) direction. Note that γ depends on our initial choice of orientation of the x-axis, which is arbitrary and of no physical importance. Top right: the next rotation is around the y-axis by β, which is the angle between p 3 and s. Bottom left: the lepton momentum with a negative y component is identified as p 2 . Its projection onto the x − y plane makes an angle α ∈ [0, π) with the x-axis. We finally rotate around the z-axis by α, which puts the decay plane into the x − z plane with p 2 in the +x hemisphere. Bottom right: final result after all three rotations. Note that the angles α and β are now identical to the azimuthal and polar angles of s , respectively.

Phase space
Besides the angles in (2.14) we choose two of the three energies of the final-state leptons, which in the rest frame of the tau satisfy E 1 + E 2 + E 3 = m τ . It is a straightforward exercise to show that the totally differential decay rate is then given by Since |M ↑ | 2 does not depend on the angle γ one can trivially integrate over the allowed range. The phase space for the energies E 2 and E 3 is nearly triangular with edges that are smoothened out by the light-lepton mass m. One way of expressing the phase-space boundaries is, for example, is a function of E 2 that is near 1, except when E 2 is near the endpoint where d 2 rapidly falls to zero. Sometimes it is more advantageous to express these energies in Dalitz-like invariants s ij = (p i + p j ) 2 , for example when addressing the momentum flowing through the virtual photon in the radiative operators in (2.2). In terms of s ij the energies are For later reference we also state the phase-space region in which cos φ 23 ≥ 0, which is Similarly, for cos φ 12 ≥ 0, we have the restrictions

Definition of the observables
In this section we will use the formulae from the last section to define observables, which can be measured even with sparse data samples. All the observables will be defined from the coarsely sliced phase space using the angular variables, so that the observables correspond to partial rates and forward-backward asymmetries with respect to the angles, which also involve the direction of the tau-lepton polarization.

Coupling bilinears
Since lepton-flavour violating processes are rare, we focus on observables where the available phase space in α, β, E 2 , E 3 is at least partially integrated. There are many such observables, and a measurement of them will allow us to draw conclusions on the underlying operator structures. We start by considering the integration of J i over the energies E 2 and E 3 in some region "R" of the phase space. The results are bilinears in the couplings ξ in (2.4), and can be expressed in terms of 6 × 6 hermitian matrices A (R) i that depend on the region R, to wit (3.1) From this notation it is obvious that the doubly differential decay rate in α, β is obtained by integration over the full energy phase space (2.19), 2) from which asymmetries using α, β are readily calculable. For example, the difference in the partial rates Γ(cos β > 0) − Γ(cos β < 0) will only involve A . We may express this difference as the convolution integral over the differential decay rate (3.2) with the weight function [θ(cos β) − θ(− cos β)], where θ is the Heaviside step function. In order to demonstrate that one can pick out each individual matrix using partial rates we note that and list the coefficients c i for a few typical weight functions in table 1. (Note the factor of π that has been absorbed into the prefactor.) We define the following partial rates Γ c ab ,

JHEP10(2015)082
Partial rates combination which will form the basis of our simulated counting experiments in section 4, obtained from disjoint regions of the phase space in α, β, (E 2 , E 3 ). The assignments are Splitting the range of α into more than two regions is necessary for isolating A (R) 4 . Table 1 also states the linear combination of these rates that correspond to the given binning functions. It is clear that any observable constructed from the doubly differential decay rate (3.2) corresponds to a matrix that is such a linear combination of the four A = 0 due to the fact that the two final-state leptons are identical particles and that J 3 and J 4 are antisymmetric in the interchange of E 1 ↔ E 2 . In other words, any weighted integral over (3.2) is just a linear combination of the total decay rate and the above-mentioned asymmetry. Note also that A {full} 4 = 0 is a necessary condition for our discussion on the unpolarized decay rate around equation (2.7) as is evident from the first line in table 1.
If, however, we consider only the part of the phase space (2.22) in which the angle between p 2 and p 3 is between 0 and π/2 (denoted by R= {cos φ 23 > 0}), all four matrices A (cos φ 23 >0) i are non-zero and contribute new information. (Note that R= {cos φ 23 < 0}, which is the complement to the full phase space, would then not yield more information.) Similarly the angle between the two leptons carrying momenta p 1 and p 2 can be utilized. The region R= {cos φ 12 > 0} contributes two more non-zero matrices as again A = 0. We therefore count eight different partial rates from which to construct observables.
Since the vector of coupling constants, ξ, contains six complex parameters, i.e. twelve real unknowns, the above partial rates do not suffice to solve the system. One way out would be to obtain independent information on the coefficients ξ 5 and ξ 6 from the radiative decays τ → µγ. Alternatively, one could, of course, divide the phase space into more (i.e. smaller) regions R, but this would only make sense if sufficient signal events had been measured.

JHEP10(2015)082
For the time being, we will restrict ourselves to a handful of benchmark scenarios that will be defined and analyzed in section 4.

Calculation of the matrices
The tree-level calculation of the squared amplitude is straightforward, and we have collected the various parts in appendix A. The mass m of the light leptons in the final state has been kept finite in our calculations. The task is now to integrate the resulting functions J i (E 2 , E 3 ) over (part of) the phase space R to arrive at the matrices A (R) i . Their entries are functions of the mass ratio ratio = m/m τ , which is small. It is tempting to state the results analytically in an expansion of , which can be done using the method of regions [21] (see also [22] for a review). The leading-power contributions involving the leptonic operators are calculable by simply setting m to zero. Interference between radiative operators and leptonic operators give rise to a single logarithm of the form ln after integrating over the phase space, even at leading power. The appearance of this single logarithm is sometimes used to justify a rather pragmatic approach in which the singularity arising from setting m = 0 is regulated by an artificial quantity δ, leading to ln δ after integration, see e.g. [15] and the discussion in the next section.
The most interesting effect happens in the interference of two radiative operators, where a double logarithm (ln 2 ) can appear. This can happen because two of the final-state particles are identical, leading to contributions in the differential decay rate proportional to 1/s 13 (photon propagator in 2 ( 1¯ 3 )|H rad |τ ) and also proportional to 1/s 23  1 , see e.g. |M rad | 2 in (A.3). We give a brief discussion of the strategy of regions on two illustrative examples which can be found in appendix B. For practical purposes, however, it suffices to evaluate the matrices numerically, including the full m dependence without expanding in , and we list the results in appendix C for both τ → 3µ and τ → 3e decays.

Comparison with the literature
The proposal of utilising the polarization vector of the tau for asymmetries [13] was preceded by an analysis within the context of the littlest Higgs model with T parity [15], in which the authors also considered lepton-flavour violating decays of polarized (anti-) taus into light leptons, among other channels. Contrary to our assumption stated after equation (2.1) that (LR)(LR) and (RL)(RL) operators are suppressed by (v 2 /Λ 2 ), Goto et al. included them in their work. Furthermore only the leading-power expressions in (m/m τ ) are considered -massless final-state leptons, in other words. With this approximation there is no interference between operators of different chirality, which requires a mass insertion. In this limit the couplings of the above-mentioned extra operators, called g I Ls and g I Rs enter simply by adding them to g (LL)(LL) V and g (RR)(RR) V (in our notation) in quadrature, see (A8c) of [15].
When neglecting final-state lepton masses one encounters unregulated collinear and soft divergences in the contribution from the radiative operators at the boundaries of phase JHEP10(2015)082 space, where the intermediate photon propagator can become on-shell. Goto et al. introduced a cutoff parameter δ > 0 in order to stay away from this boundary, so that their expressions become functions of δ rather than (m/m τ ). Here, in this present work, we went further in that the final-state lepton masses remain finite, and therefore all interference terms are accounted for and all soft and collinear divergences are naturally regulated.
We have compared some of the terms in the fully differential decay rate stated in [15] with the corresponding leading-power approximations of our results in appendix A and find agreement after accounting for the different setups.

Phenomenology
We start with assuming some values for the LFV couplings, on the basis of which we generate N total events and bin them into N c ab counts congruent with the definition of the partial rates in (3.4). Our goal is then to reconstruct the couplings from a simple and straight-forward least-square fit of the binning counts, i.e. fitting the bin probabilities, to the fraction of events in that bin, N c ab /N total . Note that the bin probabilities are invariant under a simultaneous rescaling of all couplings, ξ −→ ωξ. In order to avoid this flat direction in the χ 2 function we impose a condition that breaks the rescaling invariance, to wit ξ † A {full} 1 ξ ≡ 1 . We stress that a thus fitted result only reflects the relative strength of the couplings to each other, and not the values of the couplings themselves. Those can be inferred from the total decay width of this process and the lifetime of the tau, which is not our main focus here. 4 As we have seen in section 3.1 there are eight independent matrices from which we can calculate observables. One of them, A {cos φ 23 >0} 4 , probes the imaginary parts of the couplings and does not contribute if the couplings are real. One may choose many different sets of observables to fit for the couplings, but here we simply use the bins from which these observables are calculated themselves. The least-square fit is thus performed by minimizing the function where µ is the Lagrange multiplier used to impose the above normalization condition. We use the unsophisticated statistical estimator for the uncertainty of each bin count ∆N c ab = N c ab if N c ab = 0 and 1 otherwise. In the following, as a premature study, we are assuming only real couplings, for simplicity. In this case the above-mentioned matrix A {cos φ 23 >0} 4 does not contribute, and there 4 We stress however, that a measurement of the total decay width itself does depend on the underlying distribution over the energy and angular phase space, and therefore -in case of only a few eventsrequires the information about the relative size of the individual couplings in an essential way. This also affects the experimental procedure to generate bounds on Γ(τ → 3 ). is no need for splitting the range of α into four distinct bins. Hence we will combine the bins with a = 1 and a = 2 into one single bin, as well as the bins a = 3 and a = 4. We will therefore fit 12 observables (seven of which are independent) to 6 real unknowns.

Muonic final state
A typical entry in ξ is therefore O(0.1) to which we may compare the errors. From the probabilities we generate N total events that are distributed in the bin counts N c ab . These serve as the output of our virtual experiment. Notice that in this example, all couplings have entries of similar size (or happen to be zero). At this point we pretend that the couplings ξ are unknown and proceed with the fit. The outcome ξ fit will in general deviate from the scenario input ξ in , and we may form the deviation vector ∆ξ = ξ fit − ξ in , not to be confused with the individual fit errors δξ fit . We then repeat this virtual experiment one thousand times and display the entries of the deviation vector ∆ξ in histograms, which peak around zero with a certain width (see figure 3). If a histogram shows a normal distribution then the width (2σ) coincides with twice the mean of the fit errors δξ fit . In table 2 the resulting widths of the histograms (= 2σ i ) are listed for a few numbers of events N total . We calculate σ i as the standard deviation of the data underlying the histogram. The reader should keep in mind that the input values of the couplings are ξ T in ≈ (0.14, 0.04, −0.43, 0.00, 0.26, 0.13) which is to be compared to the uncertainty of the fit results. Regardless of whether we interpret σ i or δξ fit i as the typical uncertainty of the coupling ξ fit i , the conclusion is that with a few tens of events only the couplings of the radiative operators, ξ 5 and ξ 6 , can be extracted in any meaningful way, while the couplings of the leptonic operators, ξ 1 through ξ 4 require hundreds -if not thousands -of events. This pattern is a direct consequence of the magnitude of the entries in the A (R) i matrices, see e.g. (B.14): the largest entries are the log-enhanced diagonal elements of the radiative sector, which is why the observables are quite sensitive to ξ 5 and ξ 6 (unless these couplings happen to be generically suppressed in a particular class of NP models under consideration).
One may ask if the particular choices of parameter values in scenario a) have any influence on the outcome of this sensitivity study. We therefore repeat the above procedure JHEP10(2015)082 Figure 3. Example histograms of the difference ∆ξ 1 = ξ fit 1 − ξ in 1 (left) and ∆ξ 2 (right) for N total = 100 events in scenario a). The standard deviation of the left distribution is σ 1 = 0.288, while the mean of the individual fit errors is δξ fit 1 = 0.256. This example shows a shape that is more peaked than a normal distribution. The right distribution has σ 2 = 0.362 which is less than δξ fit 1 = 0.421, and the shape is more box-like than a normal distribution. with the following twist:

Scenario b)
For each of the one thousand samples the couplings are chosen at random and rescaled to abide by ξ † A We generate a random set of couplings, drawn from a finite interval, which we choose symmetric around zero and universal for all couplings ξ i , and rescale to build a sample of vectors ξ. This procedure removes the preference towards particular values for the couplings, although implicitly it assumes that they are still of the same order of magnitude. Notice that, as a consequence for the condition (4.4), the resulting distributions in the sampled couplings ξ i are not flat. In figure 4 we show the distributions for ξ 1 and ξ 5 , exemplary for couplings to leptonic and radiative operators, respectively. Note that the means for the absolute values are roughly around 0.2 in both cases, in accordance with our previous observation that couplings are of order O(0.1). We present the corresponding fit results in table 3, which may be juxtaposed to the previous scenario. In general the widths σ i have increased due to the fact that one is convoluting the "fixed-coupling results" with the "coupling distributions", leading to more pronounced shoulders in the resulting distributions. However, the mean fit uncertainties δξ fit i are roughly the same, with some entries larger than their fixed-coupling counterparts in table 2, and some entries smaller. Our general conclusion remains unchanged.
We stress again that even with these randomized couplings for each virtual experiment there is still a build-in assumption: that all couplings are of the same order of magnitude with a flat distribution. However, depending on the explicit New Physics model (see e.g. [23][24][25][26][27][28][29][30]) the relative importance of radiative and 4-lepton operators can be quite different. We therefore also considered a scenario in which the radiative operators may be loop induced and are therefore accompanied by Wilson coefficients that could be much smaller compared to those of the leptonic operators. We may name this Scenario b'). To be concrete, here we assume a flat distribution for the couplings of the radiative operators that  Table 3. Histogram widths 2σ i as the error estimator for samples containing N total events, as well as the mean fit uncertainty for Scenario b). have a smaller support interval, by a factor of 1/(4π) 2 . Again the absolute errors quoted in table 3 remain the same, up to small variations. However, the corresponding distributions in the spirit of figure 4 are such that ξ in 1 is of O(1) -representativ for all couplings to leptonic operators -and the ones to radiative operators are of O(1/100). In this case it is the leptonic couplings that can be determined with a few dozen events, while the radiative ones are elusive. If this pattern were to be observed, we could also expect a suppression of the τ → γ decay channel. Ultimately a global fit with all relevant decay channels would be in order.
Using asymmetries to distinguish two scenarios. In this subsection of the analysis the goal is to analyze the discriminating power of our approach to distinguish two different benchmark scenarios that may reflect the dynamical effects of some classes of NP models. Let us say that the dimension-6 operators (2.3) are induced by parity-violating interactions, i.e. one that couples only to left-handed taus or only to right-handed taus,

Scenario c)
Only left-handed taus participate  Table 4. Central values and standard deviations for the asymmetry A β in both scenario c) and d), using 1000 virtual experiments.
There are several asymmetries one can construct from the angles α, β, φ 12 , φ 23 . Some of them have to be combined, for example the asymmetry in cos α: according to table 1 it is proportional to A = 0, so one needs to further cut on the full phase space to gain information. Let us instead look at the asymmetry of the angle β, to wit Given the above two models we can readily calculate the theoretical expectations of A β = 0.202 (Scenario c) and A β = −0.202 (Scenario d). Even with very few total events one can indeed already obtain an impression which scenario to prefer, as is shown in table 4. Here, again, we have repeated the virtual experiment 1000 times to produce a distribution of results from which we estimate the typical error as the standard deviation of the distribution.
The central values of the distributions fluctuate a little around the theoretical expectations due to the finite number of virtual experiments, but even with a rather small N total the likelihood for one scenario over the other is significant. 5 The "Degree of Separation" between the two scenarios is calculated from the overlap of two normal distributions, ρ c (A) and ρ d (A), with the given central values and standard deviations of scenario c) and d), respectively, With this definition the Degree of Separation is zero for completely overlapping distributions and asymptotically approaching 1.0 for distributions that are far apart.

Electronic final state
When the final-state leptons are electrons the operators and Wilson coefficients are different and in general independent from the above-mentioned decay into muons. Since electrons are much lighter than muons the entries of the A (R) i matrices change, and the results can be 5 Interestingly the situation is even better for the corresponding scenarios in which all radiative operators are absent, ξ5 = ξ6 = 0. The central values then shift to A β = ±0.257 and the errors remain the same, so that the separation between the two datasets becomes more pronounced.  Table 5. Histogram widths 2σ i as the error estimator for samples containing N total events, as well as the mean fit uncertainty for Scenario a).  Table 6. Histogram widths 2σ i as the error estimator for samples containing N total events, as well as the mean fit uncertainty for Scenario b).
found in appendix C. Again we start our phenomenological game analogous to the muonic case above with Although we assume here the same relative coupling strengths as in the muonic example above, we remind the reader that there is no relation to the muonic case, and the goal is simply to gain insights into our overall ability to determine the couplings from N total events. However, the normalization condition leads to typical coupling magnitudes that are about half as large as in the muonic case, see (4.3). Comparing the results in table 5 with those in table 2, we observe that the errors on the couplings of leptonic operators are somewhat smaller for electronic final states, but not by half as is the case for the couplings themselves. However, the radiative couplings are now far easier to accurately determine. For example, with only N total = 20 events the typical fit result yields ξ 5 = 0.121 ± 0.039, whereas ξ 1 = 0.065 ± 0.42 does not allow us much insight.

Scenario b)
For each of the one thousand samples the couplings are chosen at random and rescaled to abide by ξ † A In figure 5 we show the distribution obtained from randomizing the couplings and rescaling. Again the typical size of couplings is of order O(0.1). Just as in the previous case we observe no significant difference from our finding with fixed couplings (compare  table 5 to table 6). Next, we look at the ability to distinguish two scenarios,

Scenario d)
Only right-handed taus participate ξ = ω (0, 1, 0, 1, 1, 0) , The central values for the asymmetry (4.7) in these two scenarios is A β = ±0.387. We show the outcome of 1000 simulations for each N total in table 7. It is notable that even very few events suffice for the distinction of these two scenarios.

Summary
Since the decay τ → 3 is very clean, already a single event of this type would immediately imply New Physics, since the prediction of the Standard Model (extended by including the neutrino masses through the Weinberg operator) is practically zero. However, once such a decay would be observed, the nature of the underlying interaction has to be uncovered. In general this would require to study decay distributions, which is impossible with only a few events.
In this paper we have studied the LFV decays of tau leptons into three muons or electrons, including a polarization of the tau lepton. Such a polarization can be generated from running an e + e − collider in the vicinity of the τ + τ − threshold with polarized electrons JHEP10(2015)082 and/or positrons. Experimentally this could be realized at BES III, once a polarization of the beams would become possible.
We have considered a general, model-independent set-up for the interaction mediating the τ decay, which amounts to parameterizing the effective interaction in terms of a few operators. Any specific model would correspond to specific values for the coupling constants in front of theses operators, and hence even a rough measurement of these couplings could discriminate between different NP models.
However, when determining the couplings in view of very sparse data samples we are required to define proper observables, which we have discussed in this paper. All observables are of the same nature as a forward-backward asymmetry, and therefore can be measured by a simple counting experiment.
On this basis we have performed a feasibility study on how precisely one could assess the values of individual LFV couplings, based on only a small number of total signal events. It turns out that for some simplified cases (i.e. assuming short-distance coefficients to be real, or particular chiral patterns) different NP scenarios could already be distinguished with a quite small number of events. In the case that LFV could be experimentally established, our procedure could be easily extended by refining the binning for the energy phase space and by including independent information on radiative τ → γ decays.

JHEP10(2015)082
Next, And finally Integrating a constant over the full phase space. We integrate the constant 1/m 2 τ over the energies E 2 and E 3 . At leading power in there is no problem and the result is 1/8. However, at subleading power divergences appear at the border of the phase space. We regulate these divergencies by manually introducing and taking the simultaneous limit η 1,2 → 0 in the end. The first integration -over E 3 , say -is straight forward. We then distinguish three regions.

JHEP10(2015)082
• Treating E 2 as O(m τ ) yields (B.2) • When treating E 2 ∼ m one needs to integrate E 2 ∈ [m, ∞]. This gives • E 2 is near its maximum value is akin to treating s 13 ∼ m 2 and integrating s 13 ∈ [4m 2 , ∞]. We find The singular terms cancel in the sum of these contributions, and the dependence on the auxiliary scales µ i drops out as well, resulting in Integrating singular pieces from the radiative contributions.
(B.7) which has a soft E 3 and a collinear s 13 singular structure. We may identify four different regions to solve this integral. We aim to distinguish these regions in a single, one-dimensional integral over a variable that differentiates between the four regions in terms of its power counting. To this end we perform the first integration over E 3 without approximation and are left with a single integral over E 2 . Let us now defineẼ 2 as the energy variable that is smallest for maximal Finally we define the variable ρ =Ẽ 2 /E 2 , which will accomplish our goal: after substituting the E 2 = m 2 ] we find the following in the above four regions and name them 3. ρ ∼ m/m τ , label "soft", 4. ρ ∼ 1, label "hard".
Finally we need to choose a regulator common in all regions. While ρ η η→0 −→ 1 is quite suitable we find that the regulator is more advantageous as it simplifies the practical calculation. We find • Usoft region: expanding the integrand and integrating 0 ≤ ρ < ∞ we find • Soft region: expanding the integrand and integrating 0 ≤ ρ < ∞ we find • Hard region: expanding the integrand and integrating 0 ≤ ρ < ∞ we find • Uhard region: expanding the integrand and integrating 0 ≤ ρ ≤ mτ 2m − 1 − 3m 2mτ we find Armed with the techniques discussed above we are able to state some analytic results. We find in terms of = m/m τ