Tensor form factor of $D \to \pi(K) \ell \nu$ and $D \to \pi(K) \ell \ell$ decays with $N_f=2+1+1$ twisted-mass fermions

We present the first lattice Nf=2+1+1 determination of the tensor form factor $f_T^{D \pi(K)}(q^2)$ corresponding to the semileptonic and rare $D \to \pi(K)$ decays as a function of the squared 4-momentum transfer $q^2$. Together with our recent determination of the vector and scalar form factors we complete the set of hadronic matrix elements regulating the semileptonic and rare $D \to \pi(K)$ transitions within and beyond the Standard Model, when a non-zero tensor coupling is possible. Our analysis is based on the gauge configurations produced by ETMC with Nf=2+1+1 flavors of dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and charm quarks with masses close to their physical values. We simulated at three different values of the lattice spacing and with pion masses as small as 220 MeV. The matrix elements of the tensor current are determined for plenty of kinematical conditions in which parent and child mesons are either moving or at rest. As in the case of the vector and scalar form factors, Lorentz symmetry breaking due to hypercubic effects is clearly observed also in the data for the tensor form factor and included in the decomposition of the current matrix elements in terms of additional form factors. After the extrapolations to the physical pion mass and to the continuum and infinite volume limits we determine the tensor form factor in the whole kinematical region accessible in the experiments. A set of synthetic data points, representing our results for $f_T^{D \pi(K)}(q^2)$ for several selected values of $q^2$, is provided and the corresponding covariance matrix is also available. At zero four-momentum transfer we get $f_T^{D \pi}(0) = 0.506 (79)$ and $f_T^{D K}(0) = 0.687 (54)$, which correspond to $f_T^{D \pi}(0)/f_+^{D \pi}(0) = 0.827 (114)$ and $f_T^{D K}(0)/f_+^{D K}(0)= 0.898 (50)$.

(q 2 ) corresponding to the semileptonic D → π(K) ν and rare D → π(K) decays as a function of the squared four-momentum transfer q 2 . Together with our recent determination of the vector f Dπ(K) + (q 2 ) and scalar f Dπ(K) 0 (q 2 ) form factors we complete the set of hadronic matrix elements regulating the semileptonic D → π(K) ν and rare D → π(K) transitions within and beyond the Standard Model, when a non-zero tensor coupling is possible. Our analysis is based on the gauge configurations produced by the European Twisted Mass Collaboration with N f = 2 + 1 + 1 flavors of dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and charm quarks with masses close to their physical values. We simulated at three different values of the lattice spacing and with pion masses as small as 220 MeV and with the valence heavy quark in the mass range from 0.7 m phys c to 1.2 m phys c . The matrix elements of the tensor current are determined for a plethora of kinematical conditions in which parent and child mesons are either moving or at rest. As in the case of the vector and scalar form factors, Lorentz symmetry breaking due to hypercubic effects is clearly observed also in the data for the tensor form factor and included in the decomposition of the current matrix elements in terms of additional form factors. After the extrapolations to the physical pion mass and to the continuum and infinite volume limits we determine the tensor form factor in the whole kinematical region from q 2 = 0 up to q 2 max = (M D −M π(K) ) 2 accessible in the experiments. A set of synthetic data points, representing our results for f Dπ(K) T (q 2 ) for several selected values of q 2 , is provided and the corresponding covariance matrix is also available. At zero fourmomentum transfer we get f Dπ T (0) = 0.506 (79) and f DK T (0) = 0.687 (54), which correspond to f Dπ T (0)/f Dπ + (0) = 0.827 (114) and f DK

Introduction
Precise measurements of hadron weak decays can constrain the Standard Model (SM) and place bounds on New Physics (NP) models. The semileptonic and rare transitions between pseudoscalar (P) mesons can be parametrized, in all extensions of the SM, in terms of three form factors, namely the vector f + , the scalar f 0 and the tensor f T ones. New particles beyond the SM, such as those appearing in NP models with supersymmetry or in a fourth generation or in composite Higgs sectors, can alter the Wilson coefficients of the effective weak Hamiltonian that describes physics below the electroweak scale. Whatever these unknown particles may be, the hadronic physics remains the same.
Recently in Ref. [1] we presented the first N f = 2 + 1 + 1 lattice QCD (LQCD) calculation of the vector and scalar form factors f Dπ(K) + (q 2 ) and f Dπ(K) 0 (q 2 ) governing the semileptonic D → π(K) ν decays within the SM. We employed the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and charm quarks with masses close to their physical values [2,3]. In this work we complete the set of operators relevant for the D → π(K) transitions by analyzing the matrix elements of the weak tensor currentcσ µν d(s) to obtain the tensor form factor f Dπ(K) T (q 2 ). The latter may enter both in the semileptonic D → π(K) ν decays as a contribution beyond the SM and in the rare decays driven by the transition c → u + − , which are loop-suppressed in the SM as they proceed through flavor-changing neutral currents.
Following Ref. [1] we have evaluated the tensor form factor f Dπ(K) T (q 2 ) in the whole range of values of q 2 accessible in the experiments, i.e. from q 2 = 0 up to q 2 max = (M D − M π(K) ) 2 . In our calculations quark momenta are injected on the lattice using non-periodic boundary conditions [4,5] and the matrix elements of the tensor current are determined for plenty of kinematical conditions, in which parent and child mesons are either moving or at rest.
The data coming from different kinematical conditions exhibit a remarkable breaking of Lorentz symmetry due to hypercubic effects for both D → π and D → K form factors. The presence of these effects was already observed in Ref. [1] in the case of the vector and scalar form factors, and there we presented a method to subtract the hypercubic artifacts and to recover the Lorentz-invariant form factors in the continuum limit.
Besides Ref. [1] hypercubic effects were never observed in the context of the D → π(K) transitions. Previous lattice calculations used only a limited number of kinematical conditions (typically the D-meson at rest). This limitation obscures the presence of hypercubic effects in the lattice data. Moreover, in Ref. [1] we found that the hypercubic artifacts strongly depend on the difference between the parent and the child meson masses. This is an important issue, which warrants dedicated investigations. If this is the case, the hypercubic artifacts may play an important role in the determination of the form factors governing the semileptonic B-meson decays and it is therefore crucial to have them under control.
In Ref. [1] the subtraction of the hypercubic effects was achieved by considering a decomposition of the current matrix element which contains, beside the usual Lorentz-covariant part, additional hypercubic structures proportional to a 2 . The form of this decomposition depends on the Dirac structure of the current and, therefore, it is interesting to further validate the method by applying it also to the case of the tensor current.
In this work we present the subtraction of the hypercubic artifacts and the determination of the Lorentz-invariant tensor form factor f Dπ(K) T (q 2 ) after the combined extrapolations to the physical pion mass and to the continuum limit. A set of synthetic data points, representing our results for f Dπ(K) T (q 2 ) for several selected values of q 2 , is provided (see later on Tables 5 and  6). The full covariance matrix, which includes the synthetic data points of the vector and scalar form factors obtained in Ref. [1] and those of the tensor form factor from the present work, is available upon request.
At zero four-momentum transfer our results are and where the errors include both statistical and systematic uncertainties, added in quadrature. The paper is organized as follows. In Section 2 we describe the simulation details. In Section 3 we present the computation of the tensor form factors f Dπ(K) T (q 2 ) using the matrix elements of the weak tensor current relevant for the D → π(K) transition, obtained from two-point and three-point correlation functions. In Section 4 the evidence of Lorentz symmetry breaking in the momentum dependence of the form factors is presented and discussed. In Section 5 we describe the strategy adopted in order to extract the physical, Lorentz invariant, tensor form factors. This is based on a global fit of the data corresponding to all lattice ensembles, studying simultaneously the dependence on q 2 , the light-quark mass m and the lattice spacing a, and using a phenomenological Ansatz to describe the hypercubic effects. In Section 6 the results for f Dπ(K) T (q 2 ) from the global fit, as well as for the ratio f

Simulation details
The gauge ensembles used in this work have been generated by the ETMC with N f = 2 + 1 + 1 dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and the charm quarks [2,3]. The ensembles and the simulations are the same adopted in Ref. [1] for the determination of the vector f Dπ(K) + (q 2 ) and scalar f Dπ(K) 0 (q 2 ) form factors. Here, in Table 1 we recall the basic simulation parameters and the masses of the π, K and D mesons corresponding to each ensemble.
The gauge fields are simulated using the Iwasaki gluon action [7], while sea quarks are implemented with the Wilson Twisted Mass Action at maximal twist [8,9,10]. In order to avoid the mixing of strange and charm quarks induced by lattice artifacts in the unitary twisted-mass formulation we have adopted the non-unitary setup described in Ref. [11], in which the valence strange quarks are regularized as Osterwalder-Seiler (OS) fermions [12], while the valence up and down quarks have the same action as the sea. The use of different lattice regularisations for the valence and sea quarks of the second generation preserves unitarity in the continuum limit    and does not modify the operator renormalization pattern in mass-independent schemes, while producing only a modification of discretization effects. Moreover, since we work at maximal twist, physical observables are guaranteed to be automatically O(a)-improved [10,11].
The QCD simulations have been carried out at three different values of the inverse bare lattice coupling β, to allow for a controlled extrapolation to the continuum limit, and at different lattice volumes. We have simulated quark masses in the range from are the physical values of the average up/down, strange and charm quark masses respectively, as determined in Ref. [6]. The lattice scale is fixed using as input the experimental value of the pion decay constant f π from PDG [13]. The values of the lattice spacing are: a = {0.0885 (36), 0.0815 (30), 0.0619 (18)} fm at β = {1.90, 1.95, 2.10} respectively. The lattice volume goes from 2 to 3 fm and the pion masses range from 220 to 500 MeV. Following Ref. [1] we make use of the input parameters (values of quark masses and lattice spacings) obtained from the eight branches of the analysis carried out in Ref. [6]. The various branches differ by: i) the choice of the scaling variable, which was taken to be either the Sommer parameter r 0 /a [14] or the mass of a fictitious P-meson made of two valence strangelike quarks aM s s ; ii) the fitting procedures, which were based either on Chiral Perturbation Theory (ChPT) or on a polynomial expansion in the light quark mass (for the motivations see the discussion in Section 3.1 of Ref. [6]); iii) the choice between two methods, denoted as M1 and M2 which differ by O(a 2 ) effects (see, e.g., Ref. [15]), used to determine non-perturbatively in the RI -MOM scheme the values of the mass renormalization constant (RC) Z m = 1/Z P [6]. The use of the two sets of RCs lead to the same final results once the continuum limit for the physical quantity of interest is performed.
In what follows we make use of the local version of the tensor operator, which in the twistedmass setup renormalizes multiplicatively and is O(a)-improved [10], namely where q = {d, s}. The tensor current RC, Z T , was computed in the RI -MOM scheme and converted in the M S scheme at a renormalization scale equal to 2 GeV in Ref. [6]. The numerical values of Z T used in our analyses are listed in Table 2. (3) 1.95 0.724(4) 0.711 (2) 2.10 0.774(4) 0.767(2) Throughout this work the results corresponding to the eight branches of the analysis are combined to form our averages and errors according to Eq. (28) of Ref. [6].

Lattice calculation of the tensor matrix elements
The matrix element of the renormalized tensor current (3) between an initial D-meson state and a π(K)-meson final state can be written, as required by the Lorentz symmetry, in terms of a single form factor f Dπ(K) T (q 2 ): where P = π(K) can be either the pion or the kaon and the four-momentum transfer q is given by q ≡ p D − p P . The factor 2/(M D + M P ) is conventionally inserted to make the tensor form factor dimensionless.
In order to inject momenta on the lattice we use the same procedure adopted in Refs. [1,16] for the D → π(K) semileptonic decays and the K 3 decays. In particular for the valence quark fields we impose twisted boundary conditions (BC's) [4,5,17] in the spatial directions in order to remove the limitations resulting from the use of periodic BC's and to access the whole kinematical region for momentum-dependent quantities, like the form factors. The sea quarks, on the contrary, have been simulated in Refs. [2,3] with periodic BC's in space. As shown in Refs. [18,19], for physical quantities which do not involve final state interactions (like, e.g., meson masses, decay constants and semileptonic form factors), the use of different BC's for valence and sea quarks produces only finite size effects (FSEs) which are exponentially small. For both valence and sea quarks the BC's in time are anti-periodic. Thus, the valence quark three-momentum is given by where n x,y,z are integers and the θ values are equal in all the three spatial directions, listed in Table 3 for the reader's convenience.  The matrix elements P (p P )| T µν |D(p D ) can be extracted from the large (Euclidean) time distance behavior of a convenient combination of two-point and three-point lattice correlation functions, defined as where t is the time distance between the sink and the source, t is the time distance between the insertion of the tensor current T µν and the source, and P D 5 = icγ 5 u and P π(K) 5 = id(s)γ 5 u are the interpolating fields of the D and π(K) mesons. The Wilson parameters of the two valence quarks in the parent and child mesons are always chosen to have opposite values, i.e. r c = r s = r d = −r u . In this way the squared π(K)-meson mass differs from its continuum counterpart only by terms of order O(a 2 µ (s) Λ QCD ) [10].
The three-point correlation functions C DP Tµν (t, t , p D , p P ) have been simulated by imposing periodic BC's to the spectator valence u quark and twisted BC's (5) to both the initial c and final d(s) quarks. With this choice the π, K and D meson (spatial) three-momenta are given by where θ D(P ) can assume for each gauge ensemble the values of the parameter θ given in Table 3, which have been chosen in order to obtain values of | p D(P ) | equal to ≈ 150 MeV, ≈ 350 MeV and ≈ 650 MeV independently of the lattice spacing and volumes.
Since we have used only momenta distributed democratically in the three spatial directions, the only non-vanishing matrix elements of the tensor current are P (p P )| T 0j |D(p D ) with j = 1, 2, 3, which are equal to each other. Therefore, in order to improve the statistics, in what follows we introduce the tensor operator T sp given by The statistical accuracy of the correlators (6-7) is significantly improved by using the allto-all quark propagators evaluated with the so-called "one-end" stochastic method [20], which includes spatial stochastic sources at a single time slice chosen randomly (see Ref. [21] where the degenerate case of the electromagnetic pion form factor is discussed in details). Statistical errors on the quantities directly extracted from the correlators are always evaluated using the jackknife procedure, while cross-correlations are taken into account by the use of the eight bootstrap samplings (with O(100) events each) corresponding to the eight analyses of Ref. [6] (see Section 2).
In the case of charm quarks we adopt Gaussian-smeared interpolating fields [22] for both the source and the sink in order to suppress faster the contribution of the excited states, leading to an improved projection onto the ground state at relatively small time distances. For the values of the smearing parameters we set k G = 4 and N G = 30. In addition, we apply APE-smearing to the gauge links [23] in the interpolating fields with parameters α AP E = 0.5 and N AP E = 20. The Gaussian smearing is applied as well also for the light and strange quarks.
As is well known, at large time distances the two-point and three-point correlation functions behave as where E D and E P are the energies of the D and P mesons, while Z D and Z P are the matrix elements 0| P D 5 (0) | D( p D ) and 0| P P 5 (0) | P ( p P ) , which depend on the meson momenta p D and p P because of the use of smeared interpolating fields. The matrix elements Z D and Z P can be extracted directly by fitting the large time behavior of the corresponding two-point correlation functions using the exponential behavior given by the r.h.s. of Eq. (10). The time intervals [t min , t max ] adopted for the fit (10) are taken consistently with the one used in Ref. [1] for the calculation of the vector and scalar form factors f , and they are listed in Table 4.
The energies E D(P ) extracted from the fit (10) are consistent (within the statistical errors)  Table 4: Time intervals adopted for the extraction of the P meson energies E D(P ) and the matrix elements Z D(P ) from the two-point correlators in the light ( ), strange (s) and charm (c) sectors. The last column contains the values of the time distance t between the source and the sink adopted for the three-point correlators (7). meson mass extracted from the two-point correlator corresponding to the meson at rest. As in Ref.
[1], we use for the analysis the energy values E disp D(P ) instead of those directly extracted from the fit.
We have evaluated the three-point correlators (7) for various choices of the time distance t between the source and the sink in the case of a representative subset of the ETMC gauge ensembles. In this way the choice of t has been optimized in order to reduce the statistical noise keeping at the same time that the ground-state signals (11) change only within the statistical errors. The values adopted for t are found to be the same as those used for the vector and scalar form factors in Ref. [1] and are shown in the last column of Table 4.
The matrix elements P (p P )| T sp |D(p D ) can be extracted from the time dependence of the ratio R between the two-point and three-point correlation functions given in Eqs. (10)(11), namely where the correlation function C D(P ) 2 is given by which at large time distances behave as i.e. without the backward signal. One has: Taking the square root of R we can get the absolute value of the matrix elements P (p P )| T sp |D(p D ) , while its sign can be easily inferred from that of the correlator C DP Tsp (t, t , p D , p P ) in the relevant time regions.
According to Ref. [10], in order to guarantee the O(a) improvement of the extracted matrix elements we consider the two kinematics with opposite spatial momenta of the parent and child mesons and perform the following average: The quality of the plateau for the matrix elements T Dπ sp imp and T DK sp imp is illustrated in Fig. 1. The time intervals adopted for fitting Eq. (15) are symmetric around t /2 and equal to [t /2 − 2, t /2 + 2]. These values are compatible with the dominance of the π, K and D mesons ground-state observed along the time intervals for the two-point correlation functions. To evaluate possible excited states contaminations, we also checked that different choices in the time interval used to fit Eq. (15), namely [t /2 − 1, t /2 + 1] and [t /2 − 3, t /2 + 3], yield changes in the physical tensor form factor that are smaller than the statistical uncertainties. The standard procedure for determining the tensor form factor f DP T (q 2 ) is to assume the Lorentz-covariant decomposition (4), which for the current (9) implies where In the next Section we present and discuss the result of this determination in which, as anticipated, we find evidence of Lorentz symmetry breaking terms.  Fig. 2, where different markers and colors correspond to different values of the child meson momentum. Were the Lorentz-covariant decomposition (17) adequate to describe the lattice data, the extracted form factor would depend only on the squared four-momentum transfer q 2 (and on the parent and child meson masses). This is not the case and an extra dependence on the value of the child (or parent) meson momentum is clearly visible in Fig. 2. This is due to the fact that the lattice breaks Lorentz symmetry and it  is invariant only under discrete rotations by multiple of 90 • in each direction of the Euclidean space-time. Consequently the form factor may depend also on hypercubic invariants.

Tensor form factor and hypercubic effects
In Fig. 2 finite size effects (FSEs) might in principle play a role. However, by comparing the results corresponding to the gauge ensembles A40.24 and A40.32, which share the same pion mass and lattice spacing, but have different lattice sizes (L = 24a and L = 32a), we found that FSEs are negligible within the current statistical uncertainties.
In Ref. [1], where we observed for the first time the breaking of the Lorentz symmetry in the vector and scalar semileptonic D → π(K) form factors, we have proposed a procedure for the subtraction of these effects, which will be applied to the case of the tensor form factor as discussed in the next Section.
Before closing this Section, we remind that in Ref [1] it was observed that the hypercubic artifacts depend largely on the difference between the parent and the child meson masses. This finding is confirmed also in the present case of the tensor form factor, as shown in Fig. 3. The momentum dependence of the tensor form factor corresponding to a transition between two mass-degenerate P-mesons shows no evidence of hypercubic effects within the statistical uncertainties.

Subtraction of the hypercubic effects
Following the approach developed in Ref. [1] (see there Section 5 for details), we now address the hypercubic effects directly on the matrix elements of the tensor operator.
We start by introducing Euclidean four-momenta defined, in the case of the four-momentum transfer q µ = (q 0 , q), as so that µ q E µ q E µ = −q 2 , and we consider the following decomposition: where T E sp = i T sp . In the r.h.s. of Eq. (20) the first term is the Lorentz-covariant one while P (p P )| T E sp |D(p D ) hyp describes the hypercubic corrections given by with H DP 1 (q 2 ) and H DP 2 (q 2 ) being additional hypercubic form factors and p (22) is the most general structure, up to order O(a 2 ), that transforms properly under hypercubic rotations, is antisymmetric under the exchange of the Dirac indices and is built with four powers of the components of the parent and child momenta p µ D and p µ P . The Lorentz-invariance breaking effects are encoded in the two kinematical structures multiplying the hypercubic form factors H DP j (j = 1, 2). Consequently, we assume that H DP j depends only on q 2 (and on the parent and child meson masses). In particular, we adopt a simple polynomial form where the z variable is defined as [24,25] with t + and t 0 given by In Eq. (23) the coefficients d j 0,1,2 are treated as free parameters and may depend upon the light-quark mass m .
Due to the democratic choice of the spatial components of the meson momenta the hypercubic structure (22) cannot be determined and subtracted separately for each gauge ensemble. Instead it can be fitted by studying simultaneously the data for all the 15 ETMC gauge ensembles used in this work. Thus, we performed a global fit considering the dependences on q 2 , m and a 2 for the tensor form factor f DP T as well as the q 2 and m dependences of the hypercubic form factors H DP j .
For the form factor f DP T (q 2 , a 2 ) we adopt the z-expansion of Ref. [26], modified as in Ref. [27], viz.
where z 0 ≡ z(q 2 = 0) and the kinematical factor (M D + M P )/(M phys D + M phys P ) has been inserted in analogy with the case of the semileptonic K → π transition discussed in Ref. Eq. (26) we consider for the coefficient c DP T (a 2 ) a linear dependence on a 2 , while the pole mass M DP T is treated as a free parameter in the fitting procedure. We have checked that making the coefficient c DP T dependent on the light-quark mass m does not produce significative changes in the fitting procedure. Notice that, in the r.h.s. of Eqs. (26), the term proportional to z 2 is not independent of the other terms in the z-expansion. This constrain comes from the analyticity requirements of f DP T below the annihilation threshold √ t + = M D +M P corresponding to z = −1, as described in Ref. [26]. We have also explicitly checked that, within the current statistical uncertainties, our data are not sensitive to higher order terms in the z-expansion: the inclusion of a free parameter proportional to z 2 , plus the appropriate z 3 −term required by analyticity, produces negligible differences in the final results. At zero four-momentum transfer, the tensor form factor f DP T (0, a 2 ) reduces to For F DP T (0, a 2 ) we use the following Ansatz where ξ ≡ M 2 π /(4πf π ) 2 and the coefficients F DP , B DP and D DP are treated as free parameters in the fitting procedure. The chiral-log coefficient A DP is not known in hard pion SU(2) Chiral Perturbation Theory and therefore it is either left as a free parameter or put equal to zero. The difference between these two Ansatze is used to estimate the systematic uncertainty relative to the chiral extrapolation. An extra term proportional to m 2 in Eq. (28) turns out to be consistent with 0 and, therefore, it is not used for estimating systematic uncertainties.
We have also tried to include an extra term proportional to (aΛ QCD ) 4 . Using a value for Λ QCD equal to 0.35 GeV, we expect that the value of the coefficient of the extra term is in a natural range of order O(1). Therefore, we adopt for the coefficient of the extra term a (conservative) prior distribution equal to 0 ± 3. The results obtained using the above alternative fit show marginal differences with respect to the previous ones, and we use it to estimate the systematic uncertainty relative to the discretization effects.
Using the ingredients described above we fit all the data for the tensor current matrix element T DP sp imp coming from the 15 ETMC gauge ensembles and computed at different parent and child momenta. The quality of these fits, which include a total of 360 data and a number of free parameters between 9 and 12 depending on the specific Ansatz, is quite good with χ 2 /d.o.f. 0.9 for the D → π transition and χ 2 /d.o.f. 0.7 for the D → K one. The inclusion of the Lorentz breaking terms in the fitting procedure is crucial for obtaining a good description of the data: by putting H DP 1 = H DP 2 = 0 the value of χ 2 /d.o.f. goes up to 3.6 for the D → π case and up to 1.8 for the D → K transition.
In Fig. 4 we show the tensor form factors for the same ensemble used in Fig. 2 after the subtraction of the hypercubic contributions determined by the global fit. It can be seen that the tensor form factors depend now only on the 4−momentum transfer q 2 .
Note that an important feature of our analysis is the use of plenty of kinematical conditions corresponding to parent and child mesons either moving or at rest. This makes it possible to observe hypercubic effects in the momentum dependence of the data. As in the case of the vector and scalar form factors [1], by considering only a specific reference frame, for instance the Breit-frame in which p D = − p π(K) or the D−meson at rest, the hypercubic effects are present, but not manifest. This point is illustrated in Fig. 5, which shows the subset of our data for  the D → π (left panel) and D → K (right panel) tensor form factors corresponding only to the D-meson at rest both before and after the subtraction of the hypercubic effects determined in the global fitting procedure. Lorentz-symmetry breaking is not manifest in the limited set of data points with p D = 0, but it is not negligible.
In Ref. [1], where we studied the hypercubic effects in the vector form factor f + (q 2 ), we found that the hypercubic correction is small at q 2 = 0 and that the kinematic with the largest correction at high q 2 is the one corresponding to the child meson at rest. Conversely, in the present case of the tensor form factor f DP T (q 2 ), the comparison of the results shown in Figs. 2 and 4 indicates that mild effects are again present in the low q 2 region, while at high q 2 the data corresponding to child meson momenta different from zero get the largest hypercubic correction. This feature is directly related to the kinematical structures that multiply the hypercubic form factors H DP

Results from the global fit
The momentum dependency of the physical Lorentz-invariant tensor form factor, extrapolated to the physical pion mass and to the continuum and infinite volume limits, is shown in Fig. 6 as a cyan(orange) band for the D → π(K) transition.  In Fig. 6 our results for the tensor form factors are compared with those obtained by choosing only the kinematical configurations corresponding to the D-meson rest frame and by performing the extrapolations to the physical pion mass and to the continuum and infinite volume limits without including the hypercubic terms (22). In this way, for the limited data set corresponding to the D-meson at rest, the continuum extrapolation is based only on the discretization terms contained in Eq. (26), which are unrelated to hypercubic invariants and correspond simply to a 2 or a 2 q 2 terms. This may lead to distortions in the final results at the physical point. From Fig. 6 it can be seen that such distortions are found to be comparable with present global uncertainties within one standard-deviation. They may become more relevant as the precision of the data will be increased in the future.
In Tables 5 and 6 we provide a set of synthetic data points for the tensor D → π and D → K form factors with the corresponding total uncertainties, calculated for eight selected values of q 2 between 0 and q 2 max = (M D − M π(K) ) 2 . We provide also the values of the ratios of the tensor and vector form factors f DP T (q 2 )/f DP + (q 2 ) and of the scalar and vector ones f DP 0 (q 2 )/f DP + (q 2 ), using for the latter the results of Ref. [1].
The errors in Tables 5 and 6 take into account the uncertainties induced by: • the statistical noise and the fitting procedure itself; we stress that this error coming from a multi-combined fit includes also the uncertainty related to the removal of hypercubic effects; • the errors in the determinations of the input parameters of the eight branches of the quark mass analysis of Ref. [6]; • the chiral extrapolation, evaluated by combining the results obtained using the SU(2)inspired Ansatz (28) with a free chiral log (A DP = 0) and without the chiral log (A DP = 0); • the (Lorentz-invariant) discretization effects, calculated by comparing the results obtained either including or excluding in Eq. (26) extra terms proportional to (aΛ QCD ) 4 and adopting for the value of the corresponding parameters a (conservative) prior distribution equal to 0 ± 3.
In order to allow a direct use of the synthetic data points without using our bootstrap samples, we have calculated the covariance matrix among the synthetic data points contained either in Table 5 or in Table 6. Moreover, taking into account also the results of Ref. [1], we have calculated the full covariance matrix corresponding to the sets of synthetic data points corresponding to all the semileptonic form factors f DP + (q 2 ), f DP 0 (q 2 ) and f DP T (q 2 ) for P = π and K (as well as the full covariance matrix corresponding to our data for f DP + (q 2 ), f DP 0 (q 2 )/f DP + (q 2 ) and f DP T (q 2 )/f DP + (q 2 )). The corresponding covariance matrices are available upon request to allow to fit our synthetic data with any functional form, that can be adopted for describing the momentum dependence of the semileptonic form factors.
In Fig. 7 the tensor form factors f D→π(K) T (q 2 ) are compared with the corresponding vector ones f D→π(K) + (q 2 ) extracted from the same ETMC gauge ensembles in Ref. [1].  Table 5: Synthetic data points representing our results for the tensor form factor f Dπ T (q 2 ) and its ratio with the vector one f Dπ T (q 2 )/f Dπ + (q 2 ) (obtained in Ref. [1]), extrapolated to the physical pion point and to the continuum and infinite volume limits for eight selected values of q 2 in the range between q 2 = 0 and q 2 = q 2 max = (M D − M π ) 2 3.0 GeV 2 . The errors correspond to the uncertainties related to the statistical + fitting procedure, the input parameters, the chiral extrapolation and the discretization effects, respectively (see text). In the case of the ratio f Dπ T (q 2 )/f Dπ + (q 2 ) also the error related to finite size effects (present in the vector form factor [1]) is shown. The errors in squared brackets correspond to the combination in quadrature of the statistical and all systematic errors. In the rightmost column the synthetic data points corresponding to the ratio of the scalar and vector form factors f Dπ 0 (q 2 )/f Dπ + (q 2 ), as determined in Ref. [1], are also shown for comparison.  Table 6: The same as in Table 5, but for the D → K transition for eight selected values of q 2 in the range between q 2 = 0 and q 2 = q 2 max = (M D − M K ) 2 1.88 GeV 2 .

Conclusions
We have presented the first lattice N f = 2 + 1 + 1 determination of the tensor form factor f Dπ(K) T (q 2 ) corresponding to the semileptonic(rare) D → π(K) ν ( ) decays as a function of the squared four-momentum transfer q 2 . Together with the vector f Dπ(K) + (q 2 ) and scalar f Dπ(K) 0 (q 2 ) form factors calculated in Ref. [1], the present work completes the set of hadronic matrix elements regulating the semileptonic(rare) D → π(K) ν ( ) transition within and beyond the Standard Model, when a non-zero tensor coupling is possible.
Our analysis is based on the gauge configurations produced by ETMC with N f = 2 + 1 + 1 flavors of dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and charm quarks with masses close to their physical values. The matrix elements of the tensor current are determined for a plethora of kinematical conditions in which parent and child mesons are either moving or at rest. As in the case of the vector and scalar form factors, Lorentz symmetry breaking due to hypercubic effects is clearly observed also in the data for the tensor form factor and included in the decomposition of the current matrix elements in terms of additional form factors.
After the extrapolations to the physical pion mass and to the continuum and infinite volume limits we have determined the tensor form factor in the whole kinematical region from q 2 = 0 up to q 2 max = (M D − M π(K) ) 2 accessible in the experiments. A set of synthetic data points, representing our results for f