Charmed and $\phi$ meson decay constants from 2+1-flavor lattice QCD

On a lattice with 2+1-flavor dynamical domain-wall fermions at the physical pion mass, we calculate the decay constants of $D_{s}^{(*)}$, $D^{(*)}$ and $\phi$. The lattice size is $48^3\times96$, which corresponds to a spatial extension of $\sim5.5$ fm with the lattice spacing $a\approx 0.114$ fm. For the valence light, strange and charm quarks, we use overlap fermions at several mass points close to their physical values. Our results at the physical point are $f_D=213(5)$ MeV, $f_{D_s}=249(7)$ MeV, $f_{D^*}=234(6)$ MeV, $f_{D_s^*}=274(7)$ MeV, and $f_\phi=241(9)$ MeV. The couplings of $D^*$ and $D_s^*$ to the tensor current ($f_V^T$) can be derived, respectively, from the ratios $f_{D^*}^T/f_{D^*}=0.91(4)$ and $f_{D_s^*}^T/f_{D_s^*}=0.92(4)$, which are the first lattice QCD results. We also obtain the ratios $f_{D^*}/f_D=1.10(3)$ and $f_{D_s^*}/f_{D_s}=1.10(4)$, which reflect the size of heavy quark symmetry breaking in charmed mesons. The ratios $f_{D_s}/f_{D}=1.16(3)$ and $f_{D_s^*}/f_{D^*}=1.17(3)$ can be taken as a measure of SU(3) flavor symmetry breaking.


Introduction
Meson decay constants are important nonperturbative quantities for the study of meson leptonic decays, and their results from lattice Quantum Chromodynamics (QCD) have received much attention. The pseudoscalar meson decay constants (f P ) can be neatly used to determine the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements, if combined with experiment measurements of the corresponding leptonic decays. The newest lattice QCD average of f P can be found in the review by Flavor Lattice Averaging Group (FLAG) [1] In principle vector meson decay constants f V can also be used to determine CKM matrix elements although experimental measurements of leptonic decays of vector mesons are much harder than those of pseudoscalar mesons due to small branching ratios. Nevertheless, the leptonic decay of D * s may be expected to be measured by BES-III or Belle II in the near future for the first time for a vector meson [2]. Then the comparison of f D * s from experiment and theoretical calculation can be used to study the low energy properties of QCD.
Furthermore, decay constants of heavy-light vector mesons can be used to test the accuracy of heavy quark effective theory (HQET). Neglecting terms of O(1/m Q ), where m Q is the heavy quark mass, one has f V /f P = 1 − 2α s (m Q )/(3π) [3] from the leading order QCD calculation, which implies that the ratio f V /f P approaches one since the strong coupling constant α s (m Q ) vanishes in the infinite heavy quark mass limit. We can obtain the corrections from the higher order terms in charmed mesons through the ratio f V /f P from lattice QCD calculations. Also, the ratios f V /f P for charmed mesons are input parameters for QCD factorization studies of charmed nonleptonic B meson decays [4,5]. Another important quantity f T V is the coupling of a vector meson to the tensor current. The nonperturbative determination of the ratio f T V /f V is important in light cone QCD sum rule (LCSR) calculations of form factors in B to vector meson semileptonic decays.
In this paper, we present a lattice calculation of D ( * ) s , D ( * ) and φ meson decay constants in a lattice setup with chiral fermions, which are usually expected to be important when light flavors are involved since chiral symmetry is a fundamental property of QCD. We use overlap fermions for valence quarks and carry out the calculation on 2+1-flavor domain wall fermion gauge configurations generated by the RBC-UKQCD Collaborations. The lattice size is big enough (∼ 5.5 fm) to avoid large finite volume effects. The light sea quark mass is almost at the physical point. There have been four lattice QCD calculations of f D * s in literatures so far. Two of them were performed on 2-flavor gauge ensembles [6,7]. The other two were performed on 2+1-flavor ensembles [2] and 2+1+1-flavor ensembles [8], respectively. An unexpected large quenching effect of the strange quark was observed in f D * s and f D * s /f Ds from the 2-flavor result [6] (confirmed in [9] but with a reduced effect). While the 2-flavor result from [7] shows a much less pronounced effect. We can also check this quenching effect from the strange quark in our study.
The rest of this paper is organized as follows. In Sec. 2 we give our framework of the calculation, including the definitions of the decay constants and the lattice setup. Sec. 3 presents the details of the analyses, the numerical results and discussions. Finally we summarize in Sec. 4.

Definitions and lattice setup 2.1 Decay constants of pseudoscalar and vector mesons
The decay constant f P of a pseudoscalar meson P is defined through with p µ being the momentum of the meson. By using the partially conserved axial vector current (PCAC) relation, f P can be obtained from the matrix element of the pseudoscalar density where m 1,2 are quark masses and m P is the pseudoscalar meson mass. For overlap fermions, the quark mass and pseudoscalar densityψ 1 γ 5 ψ 2 renormalization constants cancel each other (Z P = Z −1 m ) due to chiral symmetry. This makes f P obtained from Eq.(2) free of renormalization.
The vector meson decay constant f V is given by the matrix element of the vector current between the vacuum and vector meson V as where µ (p, λ) is the polarization vector of meson V (p, λ) with helicity λ. We use the local vector current on the lattice to compute the above matrix element for convenience. The price to pay is the need of a calculation of the finite renormalization constant for the local current, which was obtained nonperturbatively in Ref. [10] for our lattice setup. Besides f V , vector mesons have another decay constant f T V which is defined through the following matrix element of the tensor current Here in the tensor current σ µν = (i/2)[γ µ , γ ν ]. Since the tensor current has a nonzero anomalous dimension, we will give values of f T V in the commonly used MS scheme and at a scale µ = 2 GeV. The matching factor from the lattice to the continuum MS scheme for the tensor current was presented in Ref. [10].

Lattice setup
Our calculation is carried out on the gauge configurations of N f = 2 + 1 domain wall fermions generated by the RBC-UKQCD Collaborations [11]. We use the gauge ensemble named as 48I with lattice size L 3 × T = 48 3 × 96 and pion mass m (sea) π = 139.2(4) MeV from the sea quarks. The lattice spacing was determined to be a −1 = 1.730(4) GeV [11], thus the spatial extension of the lattice is about La ∼ 5.5 fm. The parameters of the configurations are given in Table 1.
We use overlap fermions for valence quarks to perform a mixed action study. The mismatch of the mixed valence and sea pion masses between the domain-wall fermion and the overlap fermion, measured by ∆ mix , is 0.030(6)(5) GeV 4 [13], which is very small reflecting a small partial quenching effect. The multi-mass algorithm of overlap fermions [14] permits calculations of multiple quark propagators with a reasonable cost. We calculate propagators with a range of masses from the light to charm quark on 45 configurations. The valence quark masses am (val) q (q = l, s, c) in lattice units are given in Table 1. The deflation algorithm is adopted to accelerate the inversion by projecting out the 1000 low eigenvectors (including zero modes) of the overlap Dirac operator, which are calculated explicitly beforehand. (q = l, s, c) are the valence quark mass parameters in lattice units and the corresponding pion masses (in MeV) are from Ref. [12]. The physical charm quark mass am phy c is estimated to be around 0.73 (see below).   Table 1) for the light valence quarks for chiral interpolation. The corresponding pion masses range from 114 MeV to 208 MeV [12]. Two strange quark mass parameters are used to extrapolate to the physical strange quark mass point. The bare charm quark masses that we use are around 0.72 in lattice units, which are not small. Although for chiral lattice fermions the discretization error due to the heavy quark mass starts at O((am c ) 2 ), it could still be large. Thus we shall try to estimate the finite lattice spacing effects in our results for D-mesons.

Two-point correlators
The matrix elements in Eq. (1), (3) and (4), from which the decay constants are defined, can be derived directly from the related two-point functions with the currents being the sink operators. Since the mesons involved in this study are all the ground state hadrons, in order for the matrix elements to be determined precisely, it is desired that the twopoint functions are dominated by the contribution from the ground states. In this work, we adopt the Coulomb wall-source technique. That is to say, we perform the Coulomb gauge fixing to the gauge configurations firstly, and then calculate the two-point functions using the following wall-source operators which are obviously gauge dependent, where ψ f = u, d, s, ... and Γ = γ 5 for pseudoscalar mesons and Γ = γ i (i = 1, 2, 3) for vector mesons. From our experience [12] besides the suppression of excited states, the choice of wall source can also suppress the P -wave scattering states in the vector channels.
For the sink operators, we use spatially extended operators O Γ ( x, t; r) by splitting the quark and anti-quark field with spacial displacement r, namely, O Γ ( x, t; r) ≡ψ f 1 ( x + r, t)Γψ f 2 ( x, t). The operators with the same spatial separation r ≡ | r| are averaged to guarantee the correct quantum number, and also to increase the statistics as a byproduct. Thus the two-point functions we calculate are and where N r is the number of O Γ ( x, t; r)'s with the same | r| = r. The two-point functions C(r, t) with different r can be calculated simultaneously without expensive extra inversions. After the insertion of the intermediate states, the spectral expression of a two-point function reads where Φ n (r) is proportional to the Bethe-Salpeter amplitude 1 Nr | r|=r 0|O Γ ( 0, 0; r)|n for the n-th state. Since the r dependences of Φ n (r) are different for different states in each channel, a proper linear combination of several C(r, t)'s with different r may give an optimal two-point function Obviously, the parameterization of Eq. (9) shows that the spectral weight Φ n (r = 0) is proportional to the matrix element that defines the decay constant of a specific meson state. However, in order to get the decay constant, we need to remove the factor n|O (W ) † Γ |0 , which is the matrix element of the wall-source operator O (W ) † between the vacuum and the meson state and can be derived from the wall-to-wall correlation function 3 Numerical analyses

Meson masses
To extract the meson masses, we apply two fitting strategies. One strategy is applying correlated simultaneous fittings to the correlation functions with different r's using one (for vector mesons) or two (for pseudoscalar mesons) mass terms. The function form used in the simultaneous fits is where T = 96 and the second term in the brackets on the right hand side comes from the propagation of the correlator in the negative time direction. Φ n (r) and m n are fitted with the minimum χ 2 method. We vary the number of mass terms to two or three and check the stability of the fitting results. Within statistical uncertainties the fitted ground state mass m 0 does not depend on the number of mass terms. The upper limit of the fitting range [t min , t max ] is chosen by the following criteria. For the pseudoscalar channel t max is fixed to the maximum value where the relative errors of correlators satisfy δC/C ≤ 5%. For the vector channel t max is chosen by requiring δC/C ≤ 10%. The lower limit of the fitting range is varied in a wide range when doing the fittings and we check the stability of the results. Among all the fittings which have χ 2 /dof ≤ 1.0 and give a consistent ground state mass we then choose the earliest t min to give our final results.
In the left panel of Figure 1 we show the fitted ground state mass M D in lattice units as a function of t min . Here we finally choose the fitting range [11,18] for the D meson. In . M D (the band in the right graph) from fitting range [11,18] is compared with the corresponding effective masses from varies correlators (right panel).
the right panel of Figure 1 the obtained ground state mass M D (the band in the graph) is compared with the corresponding effective masses M eff = log(C(r, t)/C(r, t + 1)) from various correlators with different r. The data points in magenta squares are from the correlator with r = 6.32a ( r = (2, 6, 0) and permutations averaged). The ones in blue triangles are from the local sink correlator with r = 0. The ones in black circles are the effective masses from a combination of two correlators where we can tune the parameter ω and use various C(r, t) to make the effective mass plateau from C(ω, t) appear as early as possible. This leads to our second fitting strategy. Different states with a same quantum number contribute differently to the correlators C Γ (r, t). And these contributions vary as r varies. Thus it is possible to find a large r such that the contribution of the lowest excited state to ωC(r, t) cancels that to C(r = 1, t) and C(ω, t) is dominated by the ground state.  (4)(1) 1.070(3)(1) 1.070(3)(1) 1.071(3)(1) (8)(1) 1.157(8)(1) 1.158(7)(2) 1.160(6)(2) strategy I 1.160(2)(1) 1.160 (2)(1) 1.160 (2)(1) 1.162 (2)(1) strategy II In the right panel of Figure 1, the black circles show a mass plateau which starts much earlier than that from the correlator C(r = 6.32a, t) or C(r = 0, t). Therefore we can fit the combined correlator C(ω, t) easily with a single exponential term. We check that this fitting gives stable and consistent ground state mass as we vary the parameter ω. We also confirm that the results from the above two fitting strategies are in consistency.
The fitting results of am D and am D * from the two strategies with am c = 0.72 are summarized in Table 2 for comparison. The second strategy gives smaller statistical uncertainties since the mass plateau from the combined correlator appears earlier and thus data points with less errors are used in fittings. Similar advantages of the second strategy are observed in the analyses of other meson masses, therefore we adopt strategy II to obtain meson masses in the following.
The results of the pion and kaon masses are shown in Table 3. The pion mass and the combination m 2 ss ≡ 2m 2 K − m 2 π are used to fix the physical up (degenerate with the down quark) and strange quark mass respectively. From Table 3 we can see that 2m 2 K − m 2 π is independent of the pion mass (or equivalently the up/down quark mass) within the statistical uncertainties. This is exactly what we expect from the lowest-order analysis of chiral perturbation theory and it is the reason why we use this combination. The results of the meson masses and decay constants will be interpolated/extrapolated to the physical point where (a 2 m 2 π ) phys = 0.00651(3) and a 2 m 2 ss (phys) ≡ a 2 (2m 2 K − m 2 π ) phys = 0.1565(6) by using m phys π = 139.6 MeV and m phys K = 493.7 MeV [15]. Here the uncertainties come from the error of the lattice spacing. Since these uncertainties are much smaller than our statistical error or the discretization error as we will see later, we ignore them in our estimate of the systematic uncertainty.
For the mass of φ we do a linear extrapolation to the physical point a 2 m 2 ss (phys) = 0.1565 since we only have two data points as given in Table 4. For the corresponding a 2 m 2 ss at each of the two strange quark masses we use the average of the four values in the last column of Table 3. This extrapolation gives m phys φ = 1.018(17) GeV (15) with lattice spacing error included. Both m phys K * and m phys φ are in good agreement with their experiment values. This means that the finite lattice spacing effects in the study of light hadrons are smaller than our current statistical uncertainties.
Vector mesons can decay to two pseudoscalar mesons through P -wave. On our lattice the minimal nonzero momentum is 226 MeV, which is not small. The thresholds of P -wave decays for φ, D * and D * s mesons are not open on our lattice. But K * can decay to Kπ on our lattice. We observed mass plateaus for the K * meson but not for the scattering states of Kπ, which we believe are suppressed by the usage of Coulomb gauge wall source when calculating the 2-point functions [12]. The agreement of m phys K * and m phys φ (from our interpolation/extrapolation) with their experimental values tells us that it is safe to ignore the threshold effects at our current precision.
The masses of D s and D * s mesons are listed in Table. 5. We use the experimental value of D s (together with m 2 ss (phys) in the above) to set the physical charm (and strange) quark mass. With our lattice spacing we have (am Ds ) phys = 1.1378(26) by using m Ds = 1968.34(7) MeV from Particle Data Group (PDG2018) [15]. We use the following function similar to Eq. (13) to interpolate/extrapolate m D * s to the physical strange and charm quark mass point: where ∆m Ds = m Ds − (m Ds ) phys and b 3 is another free parameter. From this we obtain m phys which agrees with the experiment value 2.1122(4) GeV [15]. The interpolation/extrapolation is shown in Figure 2. The function Eq.(16) can describe the data very well.  Table. 6. The following ansatz is used to interpolate/extrapolate our numerical results to the physical quark mass point: Here the term b 2 ∆m 2 ss appears because our lattice results m Ds are not calculated at the physical strange quark mass and m phys Ds is used to set the physical charm quark mass. We get m phys D = 1.873(5) GeV and m phys D * = 2.026(5) GeV (19) for the two mesons respectively after the interpolations/extrapolations. Our D meson mass agrees with the PDG2018 value m D ± = 1.86965(5) GeV within 1σ. However our D * meson mass is heavier than the PDG2018 value m D * ± = 2.01026(5) GeV by about 1%. Thus we estimate the discretization error associated with the large charm quark mass to be about 1% in our results for the charmed meson masses. Table 6: Masses and decay constants of D and D * with statistical uncertainties. The fitting range of correlators for D is t ∈ [11,18]. The range for D * is t ∈ [10,16]. The   [10].

Renormalization constants
Before we go into the data analyses for the meson decay constants, we present first the renormalization constants (RCs) for the local vector current and the tensor current. The RCs of quark bilinear operators for our lattice setup (overlap fermions on domain-wall fermion configurations) were calculated nonperturbatively in Refs. [10,16]. For the 48I ensemble used in this work we employed both the RI/MOM and the RI/SMOM schemes to calculate those constants nonperturbatively [10]. The matching factors to the MS scheme for the local axial vector current Z A and for the tensor current (at scale 2 GeV) are listed in Table. 7. Because we use chiral fermions, we have Z V = Z A which was also confirmed numerically in Ref. [10].

f P and f V
To obtain decay constants f P and f V we perform simultaneous fits to the wall-to-point (C P/V (r = 0, t)) and wall-to-wall (C W (t)) correlators for a given meson M . These fittings are with two exponentials and the ground state mass is constrained within 10σ to its fitted result from the above strategy II as we determined the meson masses. After removing the matrix element of the source operator 0|O (W ) Γ |M from the spectral weight of C P/V (r = 0, t), we obtain 0|O Γ |M and then the decay constants f P/V by using Eqs.(2,3) and using the fitted meson mass. This fitting and calculation process is repeated for each Jackknife sample to get the statistical uncertainty of f P/V . For f V obtained from the local vector current we need to multiply it with Z V (= Z A ) as discussed in Section 3.2.1. In the following we use a superscript "bare" to indicate decay constants obtained directly from the local vector current.
The bare decay constants f V in lattice units for the φ meson at our two strange quark masses are given in the third column of Table 4. The two center values are almost the same and our statistical uncertainty is big (∼ 6%). Thus it is hard to tell the strange quark mass dependence of af bare as our final result, where the second error comes from the difference between the constant fit and the linear extrapolation and is treated as a systematic error. We use f Ds to estimate our discretization error due to the large charm quark mass since we cannot extrapolate to the continuum limit with only one lattice spacing. The decay constant in lattice units af Ds for all charm and strange quark masses are given in Table 5. One can use the function form given in Eq. (16) (replacing m D * s with f Ds ) to extrapolate/interpolate our lattice results in Table 5 to the physical charm and strange quark mass point. What we find is af phys The difference in the center values of f phys Ds calculated in this work and in our previous work (254(2)(4) MeV) [17] is 5 MeV or 2%. Since our previous result was obtained in the continuum limit, we treat this 2% difference as an estimate of the discretization error and assign it to all our decay constants for the charmed mesons in this work.
The vector meson decay constant af bare D * s from our lattice data is given in the sixth column of Table 5 At each quark mass combination we find the ratio f bare D * s /f Ds as given in the last column of Table 5. The statistical error is from Jackknife by using the Jackknife estimates of f D * s and f Ds . Then the ratio is extrapolated/interpolated to the physical quark mass point by using the function form in Eq.(16) (replacing m D * s with the ratio). What we find is (f bare D * s /f Ds ) phys = 0.999 (24). Multiplying it with Z V and assigning a 2% discretization error, we obtain The decay constants f D and f D * and the ratio f D * /f D from our lattice data are shown in Table 6. Similarly to the above analyses for f Ds and f D * s , we get Here the first error comes from statistics and the interpolation/extrapolation to the physical quark mass point by using Eq. (18) with the replacement of m D ( * ) by the decay constants or their ratio (The example of f D * is shown in the right panel of Fig. 3). For f D * the error of Z V is also included in the first error. The second error is the 2% systematic uncertainty due to the finite lattice spacing. Now we turn to the ratios f Ds /f D and f D * s /f D * which reflect the size of SU(3) flavor symmetry breaking. These ratios can be calculated in two ways. One is using our final results for f D ( * ) (s) at the physical quark mass point. By doing this we get 1.17(4) for both ratios. The other way is first calculating these ratios at our nonphysical quark masses and then interpolating/extrapolating them to the physical point by using the function form in Eq. (18). The second way gives f Ds /f D = 1.163 (14) and f D * s /f D * = 1.17(2) without including the 2% discretization error. Including this error leads to f Ds /f D = 1.163 (14) (23) which we take as our final results for the two ratios. They tell us that SU(3) flavor symmetry breaking effects are of size ∼ 17%. Our value for f Ds /f D agrees with the result from the RBC-UKQCD Collaborations in Ref. [18], which uses unitary lattice setups with eight gauge ensembles including the 48I used in this work.

f
T V /f V Because of the bad signal-to-noise ratio in C T (r = 0, t) we do not directly determine the decay constant f T V but calculate the ratio f T V /f V from the ratio of two-point functions Úm D s (mc; m 2(phy) ss )" • ·‚rm 2 π §m 2 ss Úm D s ½•Ôn : §5 •ª(J" •Ò, we can see, this ratio approaches a plateau at large t. We do constant fits to R(t) in the range [t min , t max ] to get f T V /f V , where t max is fixed to the maximum value of t with δR/R ≤ 10%. t min is varied to check the stability of the fitting results. The variation ranges of t min are indicated by the red lines in Figure. 4. We make sure all the fittings give consistent results. In this way we get the bare value of f T V /f V at each quark mass point. As an example, the numerical results of this ratio for D * s are presented in Table 8.
Then we use Eq. (16) and Eq.(18) to interpolate/extrapolate our raw data to the physical quark mass point for f T D * s /f D * s and f T D * /f D * respectively. After multiplying the results with the renormalization factor Z T /Z A (2 GeV) = 1.055(31) in the MS scheme and assigning a 2% discretization uncertainty, we find at the scale 2 GeV. Here the first uncertainty includes the errors from statistics and interpolation/extrapolation and the error of Z T /Z A (2 GeV), and is dominated by the error of the renormalization factor. The second uncertainty is from the finite lattice spacing effect.  (7) 274(7) 213(5) 234(6) 241(9) f T V /f V -0.92(4) -0.91(4) -

Summary
We calculated the decay constants f P , f V and f T V /f V of the charmed and light mesons including D ( * ) (s) and φ by using 2+1-flavor domain wall fermion gauge configurations at one lattice spacing. The valence overlap fermion has 4, 2 and 4 mass values respectively for the light, strange and charm quarks. We use the experiment values of m π , m 2 ss ≡ 2m 2 K − m 2 π and m Ds to set the physical light, strange and charm quark masses. The masses of D, D * (s) , φ and K * at the physical point are found by interpolation/extrapolation using the lowest order of Taylor expansion (i.e., a linear interpolation/extrapolation) since our valence quark masses are close to their physical values.
The masses m D , m D * s , m φ and m K * obtained from our lattice calculation are in good agreement with their experiment measurements. The D * mass we found is 1% higher than its experiment value. The center value of f Ds from this calculation is 2% away from our previous lattice QCD calculation extrapolated to the continuum limit [17]. Thus we estimate the discretization uncertainty in this work to be around 2%.
The final results of this work for the decay constants are given in Eqs. (20,(22)(23)(24)(25)(26)28). Quadratically adding together the statistical/fitting uncertainty and the systematic uncertainty, we get the decay constants in Table 9 and some of their ratios in Table 10. For the light vector meson φ the statistical error dominates the uncertainties. While for the heavy mesons the discretization error and the error from Z T /Z A (when needed) are the main sources of uncertainty. We believe our results for f T D * s /f D * s and f T D * /f D * are the first lattice QCD calculations, which can be used as input parameters for LCSR calculations of form factors in B to vector meson semileptonic decays.
Our number f φ = 241(9) MeV is lower than the N f = 2 lattice simulation result in Ref. [19], which gives f φ = 308(29) MeV. This may be due to the dynamical strange quark effects. Note our f φ is in good agreement with that in [20], which is also a 2+1-flavor lattice calculation. The experimental value of f φ can be extracted from Γ(φ → e + e − ) = 1.251(21) keV [15] by using the relation Inputting α em = 1/137.036 and m φ = 1019.461(16) MeV [15], one finds f exp φ = 227(2) MeV. Our result agrees with the experiment value at 1.5σ.
Here the two errors are from the lattice calculation and experiment, respectively. In Figure 5 we compare f D * (s) and the ratio f D * s /f Ds from this work and other lattice QCD calculations [2,6,7,8,9]. The values from 2+1-flavor and 2+1+1-flavor simulations are in consistency. There might be a tension between 2-flavor calculations and the other calculations including dynamical strange quark. More lattice QCD calculations are certainly welcome to clarify this situation.
The ratios of decay constants of charmed mesons in Table 10 show that the size of heavy quark symmetry breaking is about 10%. While the size of SU(3) flavor symmetry breaking is around 17%.
To better control the systematic uncertainty from discretization effects in our work, we need to perform our calculation at more lattice spacings in the future. Also we need to include the quark-line disconnected diagram for the φ meson two-point function.
To accurately estimate the threshold effects of strong decays of vector mesons, further studies on larger volumes are necessary.