Two-Meson Form Factors in Unitarized Chiral Perturbation Theory

We present a comprehensive analysis of form factors for two light pseudo-scalar mesons induced by scalar, vector, and tensor quark operators. The theoretical framework is based on a combination of unitarized chiral perturbation theory and dispersion relations. The low-energy constants in chiral perturbation theory are fixed by a global fit to the available data of the two-meson scattering phase shifts. Each form factor derived from unitarized chiral perturbation theory is improved by iteratively applying a dispersion relation. This study updates the existing results in the literature and explores those that have not been systematically studied previously, in particular the two-meson tensor form factors within unitarized chiral perturbation theory. We also discuss the applications of these form factors as mandatory inputs for low-energy phenomena, such as the semi-leptonic decays $B_s\to \pi^+\pi^-\ell^+\ell^-$ and the $\tau$ lepton decay $\tau\rightarrow\pi^{-}\pi^{0}\nu_{\tau}$, in searches for physics beyond the Standard Model.


I. INTRODUCTION
The study of multi-meson systems is an interesting problem as they are universal in various physical processes. An example of this is the B → K * (→ Kπ)µ + µ − decay that is induced by the flavor-changing neutral current. Such a process is highly suppressed in the Standard Model (SM), and thus sensitive to physics beyond the Standard Model (BSM). As a result it offers a large number of observables to test the SM ranging from differential decay widths and polarizations to a full analysis of angular distributions of the final-state particles. Recent experimental studies have led to some hints for moderate deviations from the SM [1][2][3]. Note that this process is in fact a four-body decay since the K * meson is reconstructed from the Kπ final state. Therefore, to handle such decay processes, the narrow-width approximation is usually assumed in phenomenological studies. However, this assumption may lead to sizable systematic uncertainties as it captures only part of the Kπ final-state interactions.
To solve this problem, one should use a complete factorization analysis that can systematically separate the low-energy final-state interaction from the short-ranged weak transition. In semileptonic processes like B → M 1 M 2 + − , the final two-meson state decouples from the leptons to a good approximation. Thus, it is guaranteed by the Watson-Madigal theorem [4,5] that the phase of the hadronic transition matrix element is equal to the phase of M 1 M 2 elastic scattering below the first inelastic threshold. More explicitly, as pointed out, e.g., in Ref. [6], the decay matrix element is proportional to a two-meson form factor, where the Dirac matrices Γ = 1, γ µ , σ µν correspond to the scalar, vector, and antisymmetric tensor currents, respectively. The choice depends on whether M 1 and M 2 are in relatively S-wave or Pwave. The relation between vacuum-to-two-meson form factors as in Eq. (1) and those appearing in heavy-meson decays can occasionally be sharpened based on chiral-symmetry relations [7].
One of the standard approaches to calculate these two-meson form factors is using chiral perturbation theory (ChPT), which is a low-energy effective theory of quantum chromodynamics (QCD) that describes the interaction among light mesons and baryons. The next-to-leading-order (NLO) ChPT calculation for the ππ scalar form factor was firstly given in Ref. [8]. Its two-loop representation and some unitarization schmes were discussed in Ref. [9]. After that, Refs. [10,11] performed more complete studies of the scalar form factors in unitarized chiral perturbation theory (uChPT), where the results of the NLO ChPT were extended to a higher energy scale around 1 GeV, which was realized by involving the channel coupling between the ππ and the KK systems to impose unitarity constraints on the form factors. The reconstruction of the scalar ππ and KK form factors based on a Muskhelishvili-Omnès representation, relying on phenomenological phase shift input, has by now a long history [12][13][14][15][16], which includes several dedicated applications in the context of BSM physics searches [17][18][19][20]. Extensions beyond 1 GeV with a new formalism including further inelastic channels were discussed in Ref. [21]. Studies of the πη isovector scalar form factor (also coupled to KK) are much rarer [7,22,23], largely due to far less experimental information on πη scattering. The Kπ, Kη scalar form factors up to 2 GeV were given in Ref. [24][25][26][27] with a coupled-channel dispersive analysis. The two-meson vector form factors for Kπ were first derived in ChPT [28], while it was mostly obtained by fitting to the data of semi-leptonic τ decays in Refs. [29][30][31]. There is a number of works for the pion vector form factor based on the Omnès dispersive representation [32][33][34][35][36][37][38][39][40][41], ChPT calculations [9,[42][43][44], a model based on analyticity [45], and in the large-N c limit [46]. Throughout the present study, we will work in the isospin limit, although also isospin-violating scalar and vector form factors have been studied in uChPT or using dispersive methods, such as in the context of a 0 -f 0 mixing [47] or for studies of second-class currents in τ -decays [48,49].
Unlike the scalar and vector cases, naturally the coupling with the antisymmetric tensor current does not exist in the SM and conventional ChPT (the energy-momentum tensor is symmetric, and has been built into ChPT up to NLO [50,51]). However, in terms of the research for BSM physics, for example in the Standard Model Effective Field Theory (SMEFT), a number of high-dimensional operators including the tensor current are necessary. Besides the conventional ChPT Lagrangian, additional terms with an antisymmetric (antisymmetric is implicit in the following discussion of the tensor part) tensor source was first given in Ref. [52], which is crucial to calculating the tensor form factors. Recently, dispersive analyses of tensor form factors in specific channels (ππ [53], πK [54], and for the nucleon [55]) have been carried out.
In this work, we will perform a study of all three kinds of two-meson form factors based on uChPT and dispersion relations. Section II gives a brief introduction to ChPT and its unitarization, where we will discuss how unitarized meson-meson scattering amplitudes can be obtained by the inverse amplitude method (IAM). The coupled-channel IAM [56] is modified by removing the imaginary parts of the t-and u-channel loops in order to restore unitarity in coupled-channel systems, which is otherwise violated in particular around the ρ-meson region in the isospin-1 sector.
In Section III, we will calculate the two-meson scalar form factors, which are then unitarized by the IAM. There, unphysical sub-threshold singularities, related to the so-called Adler zeros, will show up. To eliminate these defects, an iteration procedure based on dispersion relations is performed for each form factor, so that the improved form factors behave well in a wide energy range 0 . . . 1.2 GeV.
In Sections IV and V, we will apply the same procedure to the calculation of unitarized vector and tensor form factors, respectively. Some of the form factors obtained in this work are compared with those that have appeared in earlier works. For each kind of form factor, we will also briefly introduce their applications in corresponding phenomenological studies. This includes the application of the ππ form factor for the S-wave-dominated decay B s → f 0 (980)(→ π + π − )µ + µ − , the application of two-meson vector form factors in the two-body hadronic decays of a charged lepton l → φφ ν, where φ and φ denote light pseudoscalar mesons, and the application of two-meson tensor form factors for the BSM effects in two-body τ decays τ → φφ ν τ . Finally, Section VI contains a summary.
Various technicalities are relegated to the appendices.

A. Chiral Perturbation Theory and Its Unitarization
ChPT [8,28,57]  . This isomorphism enables one to parametrize any of these quotient elements U ∈ G/H by the eight pNGBs as where F 0 is the pion decay constant in the chiral limit, and contains the pNGB octet. Here, exact isospin symmetry is assumed, which turns off the π 0 -η mixing for simplicity. We use the convention that under SU(3) L × SU(3) R transformations U behaves as U → RU L † , with R ∈ SU(3) R and L ∈ SU(3) L . With U as the building block, the leading-order (LO) effective Lagrangian of ChPT is constructed as where . . . denotes the trace in SU(3) flavor space, χ = 2B 0 (M +s), with M the quark mass matrix, D µ U ≡ ∂ µ U −ir µ U +iU l µ , and s, l µ , r µ are the scalar, the left-handed, and the right-handed external sources. The parameter B 0 is proportional to the QCD quark condensate, 3F 2 0 B 0 = − ūu+dd+ss . Applying Eq. (4) at one-loop produces ultraviolet (UV) divergences that can be regulated using dimensional regularization and then reabsorbed into the low-energy constants (LECs) in the nextto-leading-order (NLO) Lagrangian [8,28]: 6.9 ± 0.7 ---- where L µν and R µν are field-strength tensors of the external sources The UV-finite, scale-dependent renormalized LECs {L r i } are defined as L r i = L i − Γ i λ, with the UV-divergent parts proportional to and the nonzero values for their coefficients Γ i relevant to this work are The scale dependence of these LECs is given by Some details of the loop functions occuring in the one-loop ChPT calculations can be found in Appendix A. Table I collects the numerical results for the L r i at the scale µ = M ρ that were obtained previously. The first column corresponds to the analysis up to O(p 4 ) in ChPT [28,59], and the second refers to the fit with O(p 6 ) corrections [60]. The third column corresponds to the previous fit of meson-meson scattering phase shifts and inelasticities in the coupled-channel IAM [56], which we will discuss further below.
The power counting of ChPT is organized according to the increasing power of the ratio p/Λ χ given in terms of a typical small pNGB momentum p of the order of the pNGB mass and the chiral symmetry breaking scale Λ χ ∼ 4πF π [61], where F π ≈ 92.1 MeV is the physical pion decay constant. Therefore, the perturbative expansion in ChPT is expected to break down when p/Λ χ ∼ 1. Moreover, a perturbative expansion in powers of momenta to any finite order cannot describe the physics of resonances, which are given by poles of the S-matrix on unphysical Riemann sheets. Thus, the masses of the lowest resonances in each meson-meson scattering channel limit the applicability region of ChPT in the corresponding sector.
Unitarization (or resummation) is a systematic prescription intended to extend the applicability of ChPT to higher energies, say 1 GeV, by modifying the perturbative expression such that it satisfies the full instead of only the perturbative unitarity requirement of quantum field theory.
Since unitarity is nonperturbative in its nature, in this way the lowest meson resonances may also be described. Note, however, that unitarization comes with a price: as there are various unitarization schemes, some scheme-dependence is introduced. Also, crossing symmetry is often broken in such approaches.
Let us start with a simple example. In the case of 2 → 2 multichannel scattering with the momenta of initial and final particles as p 1 , p 2 and p 3 , p 4 , respectively, we can define the partialwave amplitude T J (s) with total angular momentum J from the full scattering amplitude T (s, cos θ) where s = (p 1 + p 2 ) 2 = (p 3 + p 4 ) 2 is one of the usual Mandelstam variables, θ is the angle between p 3 and p 1 in the center-of-mass (c.m.) frame, and P J (cos θ) is the Legendre polynomial of order J. If we consider only two-particle intermediate states, then the partial-wave amplitudes should satisfy the following unitarity relation: where time reversal invariance is assumed. The indices i, j, k, and k denote different scattering channels, |p k | is the modulus of the c.m. 3-momentum in the kth channel, and s i th = (M a i + M b i ) 2 is the production threshold of the ith channel particles. Equation (11) can be written in matrix form: which gives exactly the IAM result. The explicit expressions for the scattering amplitudes classified by definite isospin states can be found in Appendix B.
An obvious shortcoming of the IAM formula is that it leads to a peak when the determinant det T (2) − T (4) approaches a minimum. This peak may be unphysical and, in terms of dispersion relations, is due to the failure to incorporate the pole contributions from the so-called Adler zero of the partial wave in the sub-threshold region. This problem can be satisfactorily resolved in the case of the single-channel IAM [64] but not for coupled channels, see Appendix C for a brief explanation of the procedure. In Sect. III C, we will introduce an effective solution based on dispersion relations for coupled-channel systems.

B. Unitarity
Although the uChPT was constructed to fulfill the unitarity relation ImT = T ΣT † , the onestep IAM solution of the partial waves actually satisfies the exact relation only above the highest threshold. In general, the unitarity relation below the highest threshold is broken due to the mixing between the left-hand and right-hand cuts of different matrix elements in T during the process of matrix inversion. This phenomenon occurs since all the particles in the initial and final states of the T (4) -matrix are treated as on shell [65]. Such a problem is well-known and also exists in other methods of unitarization [56,[66][67][68][69]. However, depending on the values of LECs, such unitarity violation is usually very small and would not cause any real problem in practical applications of IAM results. With the LEC values reported in Ref. [56], the only exception is the I = 1, J = 1 channel, see Fig. 2. The imaginary part of the partial-wave amplitude in this channel is peaked at √ s ≈ 0.77 GeV due to the existence of the ρ-resonance, and it turns out that the IAM approach leads to a breaking of the unitarity relation by as much as 20% around the ρ-peak.
As discussed in Ref. [56], this problem can in principle be solved by adopting a multi-step unitarization approach, namely to take the dimension of the T -matrix as a function of s, which changes by one unit every time s crosses a threshold. By doing so one explicitly avoids the mixing of left-hand and right-hand cuts between different matrix elements below the largest thresholds, and thus the unitarity relation is exactly satisfied in all regions. One disadvantage of this approach is that one cannot study the scattering amplitudes (and their associated form factors) below their respective production thresholds because their corresponding matrix elements simply do not exist in the scattering matrix. It is therefore highly desirable to search for a prescription that allows for , which is below the KK threshold, but above the ππ threshold. Similarly, the branch point of the t-channel πη loop occurs at s = 4M 2 K − (M π + M η ) 2 , again above the ππ threshold. If the kaons are off shell, such singularities would not be on the physical Riemann sheet of ππ scattering [65], and would not cause any problem. However, in the IAM treatment, the kaons are on shell, leading to an overlap between such left-hand cuts with the ππ right-hand cut and thus a violation of unitarity. We find that if we remove the imaginary part of the troublesome t-and u-channel loops, unitarity can be exactly maintained, see Fig. 3, where the curves for ImT and T ΣT † coincide. These loops include the t-channel ππ and πη loops in KK → KK (both I = 0 and I = 1), the t-and u-channel ππ loops in ηη → ηη, the t-and u-channel πK loops in KK → ηη, u-channel loops. Here, the LECs take the central values given in Ref. [56]. and the t-channel ππ and u-channel πK loops in Kη → Kη.

C. Global Fit of LECs
After the modification given above, the values of the LECs need to be refixed. We fit to the following data sets using the MINUIT function minimization and error analysis package [83,84]: the ππ scattering phase shifts are taken from the dispersive analysis compiled in Ref. [70] (which is perfectly compatible with the alternative Roy equation analyses of Refs. [85][86][87] at the level of accuracy aimed for with the IAM); 1 the data for the inelasiticity η 0 0 are taken from the analysis of Ref. [81]; the ππ → KK data are from Refs. [79,80] (cf. also Refs. [89,90]); the Kπ phase shifts are taken from Refs. [71][72][73][74][75][76][77][78] (cf. also the corresponding dispersive analyses [89,91,92]); the data for the πη invariant mass distribution are taken from Ref. [82], and the background is extracted from the corresponding curve in that reference.
We notice that in the NLO ChPT amplitudes for ππ → ππ, Kπ → Kπ, and KK → KK, L r 6 1 The table in Appendix D of Ref. [70] gives the ππ scattering phase shifts up to 970 MeV. The δ 0 0 data points above this energy were read off from the band in Fig. 15 of this reference, and those of δ 1 1 were taken from Ref. [88]. The latter reference also provides an analysis of δ 2 0 , which, however, does not match exactly to the data given in Ref. [70]. Thus, for δ 2 0 , we only use those of Ref. [70] up to 970 MeV. data are from Ref. [70]. For the Kπ phase shifts, the data are taken from Refs. [71] (up triangles), [72] (down triangles), [73] (circles), [74] (pentagons), [75] (rectangles), [76] (diamonds), [77] (left triangles), and [78] (right triangles). The ππ → KK data are from Refs. [79] (rectangles) and [80] (circles). The η 0 0 data are from Ref. [81]. The data for dσ πη /d √ s, as well as the corresponding background, are taken from Ref. [82]. and L r 8 always appear as the same linear combination 2L r 6 + L r 8 . Since most of the available data are on these channels, it is difficult to fix L r 6 and L r 8 independently. Thus, we fix L r 6 to the central value given in Ref. [56], and fit the other parameters to the above data. The πη invariant mass distribution is fitted with the following expression [56,93]: where c is a normalization constant to be fitted, q cm is the πη c.m. momentum, and the background is extracted from the experimental analysis [82]. A direct fit to all these data sets leads to a value of χ 2 /dof = 7.76 with the LECs given in the column "Fit 1" in Table I. The large χ 2 /dof value is due to the inconsistency among the data sets. Following Refs. [56,94], we increase the errors of the data points by hand, and find an additional error of 5% (of the central values) to all data points leads to χ 2 /dof = 1.21. The LECs from such a fit are listed in the last column in Table I, labelled as "Fit 2". A comparison of these fits, as well as the results using the central values of LECs in Ref. [56], to the data is shown in Fig. 4. The errors propagated from the data in Fit 2 are plotted as bands, which are rather narrow. One sees that the increase of δ 0 0 around the KK threshold is more abrupt in uChPT than that from the dispersive analysis [70]. Other than that, the data are well described using these different sets of LEC values.
An additional remark is in order. We find that when the LECs take certain values, T (2) − T (4) in the (I, J) = (0, 0) channel can have zeros in the physical region using this modified version of coupled-channel IAM. These are not the Adler zeros in the single-channel scattering amplitudes, and can lead to sharp kinks in phase shifts and other observables at the zeros. Such unphysical singularities also exist in the original coupled-channel IAM. However, in that case, due to the presence of nonvanishing imaginary parts from the unphysical left-hand cuts, the singularities are in the complex s-plane, and thus lead to smoother kinks. 2 Nevertheless, we checked that the bestfit LECs in both Fit 1 and Fit 2, as well as those from Ref. [56], do not have that problem. We will use the central values of these fits (the three last columns in Table I) in the study of form factors in the following to estimate the uncertainty of this method.
In what follows, we shall apply the idea above to calculate the two-meson scalar, vector, and tensor form factors. We shall also consider a dispersion-theoretical improvement that will get rid of the unphysical sub-threshold singularities due to Adler zeros.

III. SCALAR FORM FACTORS
In this section, we will give a systematic calculation for the two-meson scalar form factors in uChPT, where the IAM approach is applied. Note that the unitarization of the two-loop scalar (and vector) pion form factor was already discussed in Ref. [9]. The scalar form factor of a two-meson system is defined by the matrix element where the subscript i is again the channel index, {a i , b i } are the two mesons in channel i, and The unitarity relation of the scalar form factor reads where the subscript "0" at the partial-wave amplitude denotes the J = 0 component. A sketch is shown in Fig. 5, where the imaginary part of the form factor is caused by the on-shell configuration of the intermediate states. Furthermore, time reversal invariance leads to (T 0 ) ij = (T 0 ) ji . Now, we can simplify the expression above by taking Fq q (s) as a column vector in the channel space and write Eq. (21) as a matrix equation This is the exact unitarity relation of the multi-channel scalar form factor. We can also derive the perturbative unitarity relation by expanding Fq q (s) according to the chiral power counting where the superscriptq q and the argument s are suppressed for simplicity. Here we define the leading term of the expansion to be O(p 0 ) because it is not suppressed by the chiral expansion 0 + . . .. Therefore, the perturbative unitarity relation reads We have to first compute the ChPT results for the scalar form factor up to O(p 2 ) as input to the IAM formula. For that purpose, we need to express the scalar current in terms of ChPT fields.
The scalar current in QCD is defined as where q = ( u d s ) T . To obtain the ChPT version of this current, we start with the QCD Lagrangian and promote its quark mass matrix M to a general matrix X q , The scalar current defined in Eq. (25) is then obtained by taking the partial derivative of the Lagrangian with respect to matrix elements of X q , which results in We can now derive the scalar current in ChPT by applying the formula above to the chiral Lagrangian. We can write the scalar current as S ij = S ij is defined as the scalar current derived from the chiral Lagrangian L (n) . The outcome for the scalar currents, up to In particular, the components of S (2) arē The ChPT prediction for scalar form factors is then obtained by calculating the matrix elements for the currents in Eq. (29) with respect to two-meson states up to one loop, as shown in Fig. 6.
The full analytical results can be found in Appendix D.

B. Unitarization
If we restrict ourselves to one single channel, the IAM unitarization formula for form factors (scalar, vector, and tensor) can be derived rigorously from a dispersion relation in complete analogy to the derivation of the single-channel IAM formula for partial waves, see Appendix C for more details (for an early application of the single-channel IAM to scalar and vector form factors, see Ref. [62]). For coupled-channel form factors, a dispersive derivation of the IAM is not available, so here we offer a more empirical derivation of the unitarization formula, which we expect at least to work well above the highest production threshold of the coupled channels considered here, which is sufficient for the applications to most of the interesting processes we mentioned in the Introduction.
From Eq. (22), one notices that a possible solution to this unitarity relation for the scalar form The proof is simple: In the second equality we have used the unitarity relation for T 0 . So the question is how to choose the form of the real vector A such that its expansion reproduces the ChPT result up to O(p 4 ).
The correct choice turns out to be With this choice and a bit of algebra, we obtain the IAM formula for a unitarized scalar form factor Remember that as argued in Sec. II B, the LHCs of T (4) 0 will be transferred into F S . Thus we also need to remove the imaginary part of the troublesome t-and u-channel loops of T Real (solid lines) and imaginary (dashed lines) parts of F n S,ππ from the original IAM calculation. The LECs are taken from the original work [56] (gray line), Fit 1 (red line) and Fit 2 (blue line), respectively. In the sub-threshold region there exist unphysical peaks due to Adler zeros in the scattering amplitude (pointed out by the black arrows) and nonvanishing imaginary parts below the lowest threshold.

C. Improvement by Dispersion Relation
It is well known that the IAM generates spurious structures such as peaks that do not correspond to any physical resonance. This happens in particular in the region below the lightest two-meson production threshold. In fact, since F S = T 0 A is only a possible solution for the unitarity equation (22) above threshold (more rigorously, above the highest two-meson threshold since we are using a one-step unitarization for a coupled-channel problem), it is natural that the outcome can only be trusted above threshold. For example, Fig. 7 shows the scalar nonstrange current form factor from the IAM calculation, where the LECs are taken from Ref. [56] (gray line), Fit 1 (red line), and Fit 2 (blue line), respectively. Obviously it suffers from sub-threshold irregularities, such as unphysical peaks (the tiny peaks near 0 GeV in Fig. 7) and nonvanishing imaginary parts, in all channels.
The unphysical sub-threshold singularities due to the IAM are studied extensively in terms of dispersion relation for the case of single-channel scattering amplitudes [64,95]. There, the existence of spurious poles in the scalar partial wave is identified as a consequence of the failure to include the effect of the so-called Adler zero in the dispersion integral of the inverse amplitude. This problem can be solved by appropriately adding back such contributions in the IAM formula. Unfortunately, a similar solution is not available in the coupled-channel case because there is so far no dispersive derivation of the coupled-channel IAM formula.
On the other hand, the dispersion relations of the form factors themselves are much more straightforward. It simply takes the following form: where = means the principal-value integration, and s th denotes the lowest threshold. Here it is sufficient to employ an unsubtracted dispersion relation because F (s) falls off as 1/s or faster at large s as suggested by perturbative QCD [96]. This integral equation suggests a better way to proceed: use the IAM-predicted imaginary part of the form factor as the input to the dispersion integral, and obtain an improved real part of the form factor. The obtained real part is then used to modify the imaginary part, which will be the input of a subsequent dispersive analysis. Such a procedure can be iterated until the curves of both the real and imaginary parts of the form factor are stable. In the following we shall depict the actual procedure and outline some of the details of such iterations.
First, we use F [n] to denote the form factor after n iterations (to avoid confusion with the chiral order denoted by superscripts with parentheses, here we use square brackets). Obviously, then represents the original IAM result without undergoing any dispersive correction. To start the iteration process, in the first step we set the imaginary part of F [1] as which will be used later as an input to the dispersion integral to obtain ReF [1] . However, we have to apply one extra modification before evaluating the dispersion integral: the IAM result , which certainly does not apply to arbitrarily large s values, fails to reproduce the asymptotic 1/s-behavior. We therefore need to introduce a smooth transition between the IAM-predicted ImF at small s and the expected 1/s at large s by hand. This can be achieved by defining a modified imaginary part ImF [1] as where α [1] is a constant to be determined later, and σ(s) is a monotonically increasing activation function that satisfies σ(−∞) = 0 and σ(+∞) = 1. A simple choice of such a function is This activation function is centered at s 0 and has a width of δs. This means that ImF [1] (s) can smoothly transform from ImF [1] (s) to α [1] /s in the region (s 0 − δs/2, s 0 + δs/2). Now, we should use the modified imaginary part ImF [1] , instead of ImF [1] , in the dispersion integral to ensure the convergence at infinity. The unknown constant α [1] can be fixed by requiring that F (s) reproduces the NLO ChPT result at s = 0, i.e., F [1] (0) = ReF [1] (0) = F ChPT (0).
Once α [1] is fixed, everything in the dispersion integral is known and we can use it to numerically 0.95 1. determine ReF [1] (s). The dispersion relation guarantees that the outcome makes sense both below and above threshold.
The whole procedure above can be iterated until a stable result is obtained. At the second step, we define ImF [2] ≡ Re{T * ΣF [1] } where F [1] = ReF [1] + i ImF [1] , and then modify its UV-behavior by constructing ImF [2] (s). Notice that this construction will involve a new unknown constant α [2] which is in general different from α [1] so that it has to be re-determined. After that, we can plug ImF [2] into the dispersion integral to obtain ReF [2] . This procedure will be iterated for several times so that we can obtain a series of increasingly refined form factors F [3] , F [4] , . . ., which will eventually stabilize. Finally, the unphysical peaks such as those in Fig. 7 are completely wiped out after such a dispersive improvement.
It is worthwhile to stress that the Watson's theorem is still fulfilled perturbatively at the fixed point of the iteration procedure, and this dispersive treatment is applicable to all scalar, vector, and tensor form factors.

D. Numerical Results
In this section we show the numerical results of the unitarized scalar form factors after the dispersive improvement. To show the convergence of the iteration, as an example, Fig. 8 gives the first four iterations for the real and imaginary parts of F n S,ππ , with the LECs taken from Ref. [56]. It can be seen that the iteration successfully removes the kink due to the Adler zeros of the scattering amplitudes, and the imaginary part vanishes below the lowest threshold. The curves of all the scalar form factors after the iteration are plotted in Figs. 9-12 within the plot region 0 GeV < √ s < 1.2 GeV. Here, we use three sets of LECs: the original one from Ref. [56] (gray lines), the one from Fit 1 (red lines), and the one from Fit 2 (blue lines) to plot the form factors. The solid and dashed lines correspond to the real and imaginary parts, respectively. The parameters for the activation function σ(s) are taken to be s 0 = 1.8 GeV 2 and δs = 0.6 GeV 2 , which implies a smooth transition between the IAM result and the 1/s-behavior within the range In the n = (ūu +dd)/ √ 2 channel and thess channel, the real parts of the form factors generate sharp peaks around √ s ∼ 0.99 GeV due to the simultaneous existence of the f 0 (980) resonance and the KK-threshold within a narrow region. Our results for F n S,ππ , Fs s S,ππ , and Fs s S,KK are consistent with those presented in Ref. [97] (barring differences in overall normalization), the latter were obtained by a slightly different version of algebraic unitarization formula, incorporating only the s-channel cuts of the partial waves. However, the f 0 (980) region shown in Fig. 9 has a narrower structure compared with that given in Refs. [16,21]. Other disagreement appears in the F n S,KK form factor. In particular, the outcome of Ref. [97] does not match the NLO ChPT value at s = 0. Our result, on the other hand, guarantees such a matching as it is implemented during the determination of the coefficient α [i] . Also, we present the ηη form factors that were not calculated in that paper.
In theūs channel, the strategy adopted in Ref. [97] therein is computationally involved as one has to first discretize s → {s i } and solve the dispersion relation by inverting a huge rank matrix (in the s-space) to obtain the discretized form factor F (s i ). Our approach is much simpler because we are simply taking the IAM results above threshold as the input of the dispersion integral, and the outcomes quickly stabilize after two or three iterations. In Fig. 13

E. Applications of the Scalar Form Factors
We end this section by discussing applications of the two-meson scalar form factors, especially its s-dependence. Let us consider for example the decay B s → f 0 (980)(→ π + π − )µ + µ − , which is a four-body decay dominated by the S-wave contribution f 0 (980) → π + π − . To study this decay process, the main task is to evaluate the B s → (π + π − ) S transition matrix elements, which are parameterized as where M 2 ππ is the invariant mass square of the two-pion system. These B s → π + π − form factors can be calculated by light-cone sum rules (LCSR) and expressed in terms of the π + π − light-cone distribution amplitudes (LCDAs) [100][101][102][103][104][105][106].
According to the Watson-Migdal theorem, since the B s → π + π − transition totally decouples from the leptonic part at leading order, its amplitude must share the same phase as that of the π + π − scalar form factor Fs s S,ππ (M 2 ππ ) below the lowest inelastic threshold (which is that of KK since the four-pion channel is not considered, and the inelasiticity is known to be negligible below about 1 GeV). Accordingly, in the framework of LCSR, the S-wave π + π − LCDAs are defined almost the same as those of a single scalar meson f 0 but with the normalization factor taken as the scalar ππ form factor [6,[101][102][103][104]. For example, the twist-2 LCDA is defined as where Fs s S,ππ (M 2 ππ )B 0 stands for the original meson decay constant f f 0 in the definition of the f 0 LCDAs. The definition of all the twist-2 and 3 S-wave π + π − LCDAs as well as the explicit form of the resulting form factors [6] can be found in Appendix E. As a result, the form factors take the form Generally, the F i depend on both M 2 ππ and q 2 . However, in the case of B → K * (→ Kπ), as shown by Fig. 3 of Ref. [6], the F i 's have a much weaker dependence on M 2 Kπ than that on q 2 . Such behavior is similar to the case of B s → π + π − . This enables one to approximately suppress the M 2 ππ dependence of F i , which leads to a factorized form so that F i can be described by a suitable parametrization [107][108][109][110]. Practically, one can first fix M ππ = M f 0 (980) to extract F i (q 2 ), and then multiply it with Fs s S,ππ (M 2 ππ ) again in the form of Eq. (40) to recover the complete transition form factor F i (M 2 ππ , q 2 ). The detailed calculation can be found in Ref. [109]. On the other hand, instead of the S-wave π + π − LCDAs, one can firstly use the LCDAs of the f 0 to get the transition form factorsF i (q 2 ) of B s → f 0 , then due to the approximation leading to Eq. (40), one can write the total form factor as Equations (40) and (41) explicitly reflect that the ππ distribution of the S-wave-dominated decay B s → f 0 (980)(→ π + π − )µ + µ − is determined by the distribution of the ππ scalar form factor.

IV. VECTOR FORM FACTORS
Next we discuss the vector form factors. The matrix element of a vector current with respect to a two-meson system can be parametrized as Notice that we label the form factors F V ± as above because they are more commonly defined in the t-channel, where p b i will switch sign. From the equation of motion (EOM) ∂ µ (q γ µ q) = i(m q − m q )q q, one sees that the form factor F V −,i is not independent since it can be expressed in terms of F V +,i and the scalar form factor F S,i according to Therefore, it is sufficient to concentrate only on Fq q V +,i . The unitarity relation of Fq q V +,i is most conveniently expressed in terms of where P ij ≡ |p i |δ ij . One can then straightforwardly express the unitarity relation as Notice thatF is associated with the J = 1 partial-wave scattering amplitude. However, the relation above is rigorously true only above the highest threshold, where all elements of P are real. When s is between the lowest and highest thresholds, a more rigorous form of the unitarity relation is In particular, at the right-hand side of the equation above we have (P −1 ) * instead of P −1 so that the kinematical imaginary part of T * 1 (i.e., the imaginary part due to |p i |) below the highest threshold can be canceled by that of (P −1 ) * .

A. ChPT Result
There are two kinds of vector currents: the SU(3)-octet current V a µ (a = 1, .., 8) and the singlet current V 0 µ due to the SU(3) V and U(1) B (B stands for baryon number) symmetry, respectively. They can be defined as respectively. The easiest way to obtain such currents from the QCD Lagrangian is to first promote SU(3) V and U(1) B to local symmetries by introducing external fields v µ = T a v a µ and v (s) µ : Taking the derivative of the Lagrangian with respect to the external fields gives the currents Again, we can apply the formulae above to obtain vector currents in ChPT. The strict SU (3) symmetry of the ChPT Lagrangian up to O(p 4 ) leads to a vanishing V 0 µ . It should be noted that at O(p 6 ) a certain SU(3) breaking term can be introduced so that V 0 µ no longer vanishes. However, at O(p 4 ) we will not consider this effect. For the octet currents, we have In particular, the components of V (2) aµ in ChPT are (making use of the fact that V 0 µ =ūγ µ u + dγ µ d +sγ µ s = 0 in the meson sector): The one-loop ChPT results for the vector form factors are given in Appendix D.

B. Unitarization, Dispersive Improvement, and Numerical Results
The IAM formula for the unitarized vector form factors can be obtained directly from Eq. (32) by the replacements T 0 → T 1 and F S →F V + : This result is also required to be improved by a dispersion relation. The whole procedure is identical to that of the scalar form factors, except that now the imaginary parts of the vector form factors that enter the dispersion integrals should be taken as where i is the number of iteration following the unitarity relation of the vector form factors.
Our final results for the vector form factors are summarized in Figs. 14-16. In thess channel, we find that the form factor Fs s V +,KK has a peak below the KK threshold, which physically corresponds to the octet part of the φ resonance, but with a mass somewhat too low compared to the physical φ. In principle, thesγ µ s current, being the pure strangeness admixture of octet and singlet vector currents, should create a resonant φ with mass 1.02 GeV, while theūγ µ u +dγ µ d current creates a resonant ω with a much lighter mass of 0.782 GeV. However, without SU(3) symmetry breaking and without an additional singlet current, these two currents are equivalent only up to a minus sign, which means that the particles they create are also identified. In other words, the ω mass and φ mass coincide at a certain value between 0.782 GeV and 1.02 GeV. This may explain why the peak emerging in Fig. 14 is below the physical φ mass. Furthermore, since this peak is below the   The ππ vector form factor is calculated in Ref. [36] through the Omnès representation where P 1 = r 2 /6 with r 2 the pion radius squared (see Refs. [39,117,118] for recent determinations of the pion charge radius from data), and a twice-subtracted dispersion relation is applied.
The Omnès representation relates the form factor with the corresponding phase shift δ 1 1 (s), which can be extracted in a rough approximation from the mass and decay width of the ρ meson using a Breit-Wigner parametrization. The scattering amplitude dominated by the s-channel ρ resonance reads where c is an irrelevant constant. The phase shift is derived as The left plot in Fig. 17 shows a comparison between the ππ vector form factor derived in this work and that from the Omnès representation, and we observe good agreement. We also compare our result for the kaon vector form factor with that from the earlier literature also using uChPT [44].
However, in that paper only the electromagnetic (EM) form factor of the physical K + K − state is given, which is defined as  where i = u, d, s and e i = 2/3, −1/3, −1/3. Generally, due to isospin symmetry, our vector form factors derived above are related to this EM form factor according to Note that due to charge conservation we must have F EM V,K + K − (0) = 1, which can be checked by our results (using the SU(3) symmetry constraintūγ µ u +dγ µ d = −sγ µ s). As explained before, without SU(3) breaking, the last two terms in Eq. (58) will contribute the same peak between the ω mass and the φ mass. On the other hand, in Ref. [44], the two resonances ω and φ are included by hand. Therefore, to perform the comparison we can replace the last two terms in Eq. (58) with the standard Breit-Wigner distributions of ω and φ : where g is the SU(3)-symmetric vector-to-two-pseudoscalars coupling constant, and g φ , g ω refer to the coupling constants of φ and ω to the electromagnetic current. The explicit definitions can be found in Refs. [119,120], where g φ = −12.89, g ω = 17.05 and g = 6.05. The decay widths are taken as Γ φ = 4 MeV and Γ ω = 8.5 MeV. The comparison of this EM form factor was shown in the right diagram of Fig. 17, which also shows nice agreement.

C. Applications of the Vector Form Factors
We end this section by discussing an application of the two-meson vector form factors in the two-body hadronic decays of a charged lepton l → φφ ν. The leading contribution to l → φφ ν is due to a single exchange of a W -boson, which, at low energies, can be approximated by the Fermi interaction. The corresponding amplitude is where G F is Fermi's constant and V qq is the Cabibbo-Kobayashi-Maskawa (CKM) matrix element.
Since the axial vector component in the hadronic matrix element vanishes due to parity, the matrix element can be expressed in terms of the vector form factors defined in Eq. (42). Considering only the spin-averaged decay, the differential decay width is a function of two kinematic variables. They can be chosen as s = (p φ + p φ ) 2 and θ, which is the angle between p φ and p l in the c.m. frame of φφ . After carrying out the phase space integration, one obtains where the φφ subscript andq q superscript in the form factors have been suppressed for notational simplicity, and Two-body hadronic decay of charged leptons is one of the commonly studied processes in the extraction of the CKM matrix elements, for example τ → Kπν τ for V us [121]. Therefore, an improved understanding of the s-dependence in the vector form factors may improve on the V us precision [122] and lead to a better reconciliation of the same quantity measured in other processes such as the kaon leptonic/semileptonic decays.

V. TENSOR FORM FACTORS
Next, we study the tensor form factors of a two-meson system, which were so far only investigated in limited channels. We define the tensor form factors through the following matrix elements: where Λ 2 is an LEC that appears when introducing external tensor sources to the chiral Lagrangian, which we shall discuss later. The unitarity relation obeyed by the tensor form factor is identical with that of the vector form factor F V + [55]: Again, it is associated with the J = 1 partial-wave amplitude.

A. ChPT Result
The derivation of tensor currents in ChPT requires the introduction of an antisymmetric Hermitian tensor sourcet µν into the QCD Lagrangian: The corresponding effective field theory was first investigated in Ref. [52]. The LO chiral Lagrangian coupled to the tensor source scales as O(p 4 ) and is given by where t µν and t µν † are given as and the convention 0123 = 1 has been used for the Levi-Civita tensor.
At NLO, there exist quite a number of corresponding operators, among which the ones contributing to tensor form factors are given as with the covariant derivative defined as They are used to cancel the divergence that occurs in the one-loop corrections to the tensor form factors. The renormalized LECs C r i and divergence coefficients γ T i are defined as: In fact, it turns out that the requirement to cancel all divergences in tensor form factors does not fix all the {γ T i } independently, but only a subset of them, which are determined to be the following constraints: The numerical values of Λ 2 and C r i at a given renormalization scale µ are obviously required to make definite predictions on tensor form factors. Unfortunately, as far as we know, no lattice data are available for these LECs. We therefore make use of the results in Ref. [123] that attempted to evaluate the effective action from first principles (under certain uncontrollable approximations).
However, a critical issue in that approach is that one cannot study the scale dependence of the renormalized LECs, and therefore the issue of how one should match their results to the renormalized LECs with the standard Gasser-Leutwyler subtraction scheme at a given scale, say µ = M ρ , remains ambiguous. To account for this issue, we assume that the LECs in Ref. [123] are given at some unknown scaleμ, and could be run to µ = M ρ by a renormalization group (RG) running which still involves an unknown coefficient ln(M 2 ρ /μ 2 ). To fix this coefficient, we refer to Ref. [124], in which the O(p 4 ) LECs L i are calculated within the same formalism. We perform a RG running of those LECs and compare them to {L r i } at µ = M ρ that are fitted to experimental data [60]. That allows us to get a best-fit value of ln(M 2 ρ /μ 2 ) by minimizing the χ 2 . We find that, with this best-fit value, the changes in the numerical values of {L i } are relatively small. Therefore, for the case of tensor LECs, we shall simply cite the results in Ref. [123], assuming that systematic errors due to the ambiguity in the renormalization scale are much smaller than the other theoretical errors quoted in the paper. Within such a framework, the coefficients C r 36 and C r 37 vanish in the large-N c limit, where N c refers to the number of colors in the QCD Lagrangian, while other nonzero coefficients are collected in Table II (readers should be alerted to certain differences in the definition of operators between Refs. [52,123] that lead to changes in numerical values of LECs and have been properly taken into account in Table II).
The tensor currents are defined as Using the O(p 4 ) and O(p 6 ) chiral Lagrangian with tensor sources, we are able to derive the tensor currents T (4) µν and T (6) µν in ChPT, respectively, which are  The renormalization scale is assumed to be µ = M ρ (see the discussion in the text).
In particular, the components of T (4)µν read The one-loop ChPT results for the tensor form factors are given similarly in Appendix D.

B. Unitarization, Dispersive Improvement, and Numerical Results
The unitarity relation for tensor form factors F T is identical to that of vector form factors, so their IAM formulae should also take the same form T .
The results are given in Figs For the form factor F ud T,ππ , we compare it with that derived in Ref. [53], where F ud T,ππ was obtained using the Omnès representation. Since according to the Watson-Migdal theorem, in the elastic region, the phases of the tensor form factors equal those of the vector form factors, δ T (s) = δ + (s), one can use the dispersion relation to obtain the normalized tensor form factor. The comparison is shown in Fig. 21, and we observe a significant difference. For instance, the sizes of the peak at s = M 2 ρ are quite different in the two calculations, and there exists a zero point in our curve above 1 GeV that does not occur in the phase dispersive representation. The differences are exclusively due to the SU(3)-breaking LECs in the tensor form factors at NLO, and therefore probably a rather large uncertainty should be associated with them. An independent cross-check is therefore highly desirable, and in Appendix F we argue that this is in principle doable through a comparison with future lattice QCD calculations of the tensor charge of the ρ-meson.

C. Applications of the Tensor Form Factors
In the previous section, we have discussed the application of the two-meson vector form factors that characterize the SM contribution to the hadronic charged-lepton decay. At the same time, these processes also provide a suitable platform for searching for the BSM physics due to their large phase space. In the literature, the BSM physics effects are included by introducing a general set of higher-dimensional operators beyond the SM. For instance, the dimension-6 effective operators (with only left-handed neutrinos) responsible for τ → π − π 0 ν τ are given by [53] where the SM Lagrangian is recovered by setting v L = v R = s L = s R = t L = 0. In particular, the t L term contains the tensor interaction. The decay amplitude for τ − (P ) → π − (P π − ) π 0 (P π 0 ) ν τ (P ) with the leptonic and hadronic sectors, respectively, given by and Therefore an improved understanding of the tensor form factors (as well as the scalar and vector form factors) will better constrain the strength of the BSM physics interactions. The same applies to other decay processes such as τ − → (Kπ) − ν τ [54,125,126] and τ − → (πη) − ν τ [127,128]. In this appendix we summarize the relevant loop functions that appear in the one-loop calculations within ChPT. First, from the tadpole diagrams we encounter the following integral: where µ is the renormalization scale and µ i = (M 2 i /32π 2 F 2 0 ) ln(M 2 i /µ 2 ), with i = π, K, η, and F 0 is the pion decay constant in the three-flavor chiral limit. Next, in a loop integral with two propagators we encounter the following two-point function: where P, Q = π, K, and η, s = p 2 , and with For the case of a single mass M P = M Q , theJ-function reads Note that the above integrals have the correct unitarity structure along the right-hand cut, which extends on the real axis from s = (M P + M Q ) 2 to infinity.

Appendix B: Isospin Decomposition of the Scattering Amplitudes and Form Factors
In this appendix we state our conventions in defining one-and two-particle isospin eigenstates, and show how to construct scattering amplitudes of definite isospin.

One-Particle Isospin Eigenstates
The one-particle isospin eigenstate is generically denoted as |φ, I, I 3 . The phases of such states are chosen such that they satisfy standard results when acted on by isospin-raising and -lowering operators:Ĵ For the pion triplet, we choose: For the kaon doublet K + , K 0 , we choose: For the anti-kaon doubletK 0 , K − , we choose: Finally, the η-particle is simply an isospin singlet: |η, 0, 0 = |η .

Two-Particle Isospin Eigenstates
Two-particle isospin eigenstates, denoted generically as |φφ , I, I 3 , are simply obtained by combining one-particle isospin eigenstates with appropriate Clebsch-Gordan coefficients. There is one complication, namely isospin eigenstates that are constructed by two particles in the same isospin multiplet should be multiplied by a factor 1/ √ 2 so that the completeness relation they satisfied is properly normalized as I,I 3 |I, I 3 I, I 3 | = 1. 4 In the following we present the two-particle isospin eigenstates that are relevant to this work:

Form Factors
Finally, we discuss how the two-meson form factors are classified according to the isospin eigenstates. The form factors of interest have the following general form where Γ is any matrix in the non-flavor space. It is obvious that with different choices ofqq and φφ there will be different form factors. However, not all of them are independent because some of them are related via charge conjugation and isospin symmetry. In this appendix, we will extract all the independent form factors in order to minimize the calculation.
Now, with each independent quark bilinear, we just need to compute one matrix element for each independent φφ group; the others are related by Wigner-Eckart theorem. Therefore, the independent form factors can be chosen as the matrix elements of the quark bilinears between the vacuum and two-particle isospin eigenstates given in Table III. such spurious pole in both unitarized partial waves and form factors. One may refer to Ref. [64] for detailed discussions of the topic.
Let us restrict ourselves to a single-channel unitarization of partial waves and form factors. Also, we shall simply use T 2 and T 4 to denote the O(p 2 ) and O(p 4 ) J = 0 partial wave for notational simplicity. First, let us recall the naïve IAM formula for partial waves, The O(p 2 ) amplitude T 2 has a zero at s = s 2 below the production threshold: T 2 (s 2 ) = 0. This is nothing but the Adler zero at the O(p 2 ) level for the S-wave. Meanwhile, the combination T 2 − T 4 has a different zero at s = s . Since in general s = s 2 , the naïve IAM formula (C1) has a spurious pole at s = s . Similarly, the naïve IAM formula for the scalar form factor suffers from a pole at s = s .
From the dispersive point of view, the existence of such spurious poles is due to the negligence of the pole contribution of various inverse amplitudes in the derivation of the IAM formula through a dispersion relation. Therefore, the problem can be resolved by appropriately adding back these contributions. Unfortunately, a dispersive derivation of the multi-channel IAM is still missing, so we can only stick to the single-channel case in this discussion.
To derive the single-channel IAM formula for partial waves, we shall consider the dispersion relation of 1/T , 1/T 2 and T 4 /T 2 2 respectively. Also, here we shall simplify our discussion by considering an unsubtracted dispersion relation of each quantity, the outcome turns out to be equivalent to choosing the subtraction point at the Adler zero of the full amplitude [64]. The dispersion relations we obtain are Here LC and P C denote the left-hand cut and the pole contributions to the dispersion integral, respectively. Since the full and perturbative amplitude satisfies the following unitarity relation on the right-hand cut: − ImT −1 (s) = Im{T 4 (s)/T 2 2 (s)} = the right-hand cut contributions to 1/T and T 4 /T 2 2 simply differ by a sign. Furthermore, one approximates LC(1/T ) ≈ LC(T 4 /T 2 2 ) by arguing that the left-hand cut contribution is weighted at low energies where the usual ChPT expansion is appropriate. With these, we can write where the explicit expressions for the pole contributions are given by , Here s A is the Adler zero of the full partial-wave T -matrix. Of course its exact value is unknown, but we can approximate it by the Adler zero of T 2 +T 4 , which is s A ≈ s 2 +s 4 where s 4 ≈ −T 4 (s 2 )/T 2 (s 2 ).
If one neglects the pole contributions then the naïve IAM formula (C1) is recovered, the inclusion of the pole contributions will eliminate the spurious pole in the sub-threshold region.
The pole subtraction for the scalar form factor follows a similar logic. Let us start by considering the dispersion relations of (F S − F S /T 2 , respectively, Along the right-hand cut, we have Therefore, the right-hand-cut contributions for (F S − F (0) S )/T and F S /T 2 are the same. Furthermore, we assume that the two left-hand-cut contributions are also approximately the same, following an argument similar to the case of the partial-wave scattering amplitude. With this we where the explicit expressions for the pole contributions are given by , .
Again, we may approximate F S and T in the formulae above by their respective ChPT expressions up to NLO.
Notice that if we neglect the pole contributions in Eq. (C9) and substitute T (s) by the naïve IAM formula for the partial-wave T -matrix, then we re-obtain the naïve IAM unitarized scalar form factor. This expression works fine above threshold but suffers from a spurious pole below threshold. On the other hand, if we use Eq. (C9) with the expression of T (s) given in Eq. (C5), then the spurious pole will be smoothly eliminated.

Appendix D: Form Factors in ChPT to One Loop
In this section, we list the NLO ChPT results for scalar, vector, and tensor form factors. Notice that we express our result in terms of the physical pion decay constant F π , which is related to F 0 by (D1)

Scalar Form Factors for the I = 0 System
The two-meson states are chosen to be exact isospin eigenstates. At NLO, calculating the Feynman diagrams shown in Fig. 6 gives 3F n S,ηη (s) = 1 − 3µ π + 4µ K − for the hidden-strangeness form factors.

Scalar Form Factors for the I = 1 System
The scalar form factors for the I = 1 meson-meson systems up to NLO in ChPT are given by Fū d S,πη (s) = 1 − µ π + 2µ K − As mentioned in Section IV A, the SU(3) singlet vector form factor ofūγ µ u +dγ µ d +sγ µ s is zero up to O(p 4 ). Therefore, it is sufficient to present the form factors ofsγ µ s only: µ π + −8M 2 K + 4M 2 π + s 2(M 2 K − M 2 π )

Vector Form Factors for the I = 1 System
The vector form factors for the I = 1 meson-meson systems up to NLO in ChPT are given by For simplicity of notation, we will defineC r i ≡ C r i /Λ 2 . The tensor form factors for the isoscalar systems up to NLO in ChPT are given by then the VMD picture provides an approximate expression of Fq q T,φφ (s) at s ≈ M 2 V , where Γ V is the total decay width of V. An interesting consequence of this formula is that one expects ReFq q T,φφ (s) to vanish and ImFq q T,φφ (s) to peak at s = M 2 V . Furthermore, with future lattice inputs of f T V,q q , the equation above serves as a consistency check of the theoretical result for two-meson tensor form factors at the vector-meson pole.