An improved calculation of the D ∗ ( s ) D ( s ) V and B ∗ ( s ) B ( s ) V couplings from light-cone sum rules

,


Introduction
The strong couplings among the heavy vector mesons (D * (s) or B * (s) meson), the heavy pseudoscalar mesons (D (s) or B (s) ) and the light vector mesons (V involving ρ, K * , ω and ϕ), represented by g D * DV and g B * BV , offer valuable insights into the non-perturbative aspects of the long-distance dynamics within quantum chromodynamics (QCD).By combining the chiral and heavy quark symmetry of QCD, the heavy meson chiral perturbative theory (HMχPT) [1][2][3][4] provides an effective framework to describe interactions between heavy and light mesons at low energy.The D * (s) D (s) V and B * (s) B (s) V couplings are essential ingredients for the extraction of the effective coupling parameter λ which governing the H * HV interactions in the HMχPT Lagrangian [3,5].
Furthermore, the D * (s) D (s) V and B * (s) B (s) V couplings play a crucial role in determining the residue of the vector and tensor form factors for D → V and B → V transitions at the D * and B * poles, as dictated by the dispersion relation [6][7][8][9].This insight enhances our understandings of the behavior of these form factors, which are vital for extracting the CKM parameters and the lepton flavour universality (LFU) parameters in the D → V ℓν and B → V ℓν semileptonic decays.Moreover, the D * (s) D (s) V couplings are crucial for understanding the final-state interactions (FSIs), which can make significant contributions to the long-distance dynamics and may contribute strong phases, potentially leading to direct CP violation in B meson non-leptonic decays [5].Due to the phase space limit, direct experimental measurements of these couplings are not feasible.Consequently, some couplings, particularly g B * Bρ and g D * DV , have been determined through calculations employing QCD techniques such as SVZ sum rules (QCDSRs) [10][11][12] and light-cone sum rules (LCSRs) [13][14][15].However, these calculations, mainly conducted at leading order [16][17][18][19][20][21][22][23], often come with relatively large uncertainties.Recently, significant progress has been made in calculating the correlation functions using QCD (SCET) factorization approach in light-cone sum rules [24].
Incorporating next-to-leading order (NLO) (α s ) corrections and subleading power (Λ/m Q ) contributions has granted us the ability to more precisely determine the spectral density of the dispersion relation for building up the LCSRs.This approach has been widely used to calculate some essential non-perturbative physical quantities, for instance, the heavy-to-light (B → π, etc.) and heavy-toheavy (B → D, etc.) transition form factors.The calculation of the correlation function based on QCD (SCET) factorization is carried out from two-distinct aspects: the hard-collinear factorization by using the light-cone distribution amplitudes (LCDAs) of light hadrons, such as in [25][26][27][28][29], and a slightly more complex factorization procedure by separating the hard, hard-collinear, collinear and soft scales with the using of the heavy-light hadrons LCDAs, e.g., the calculations in [30][31][32][33][34][35][36][37][38].
The QCD (SCET) factorization-based LCSRs has emerged as a robust technique for investigating strong couplings beyond leading order and power like g H * Hρ [8], and g H * Hπ [9].The higher-order/power contributions can significantly affect the magnitude of strong couplings, potentially altering the overall results by approximately 10% to 20% for NLO and 20% to 30% for subleadingpower effects as shown in [8].However, discrepancies between calculation results in LCSRs [8,[16][17][18] and the extrapolated coupling g B * Bρ from B → ρ form factors [28,32], using dispersion relation, have been observed.These discrepancies may stem from incomplete analyses of subleading power contributions at leading order (LO) and next-to-leading order (NLO) contributions at higher-twist, along with uncertainties arising from different shapes of the duality region at NLO [9].Another potential explanation for this inconsistency is the inclusion of contributions from excited heavy meson states to the dispersion relation [39], which introduces considerable model dependence and raises questions about effectively isolating ground-state contributions within the duality region.
Given the importance of these strong couplings for our understanding of QCD long-distance dynamics, as introduced previously, this paper aims to present a comprehensive investigation of the D * (s) D (s) V and B * (s) B (s) V strong couplings using LCSRs.The essential new ingredients are summarized as follows.
• We improve the calculation accuracy of the established light-cone sum rules from two perspectives.Firstly, we incorporate more complete subleading power corrections1 (up to δ 4  V , δ 1 s , and (δ V δ s ) 1 ) and higher-twist corrections (up to two-particle twist-5 and three-particle twist-4) of the LCDAs of light vector mesons at leading order and the NLO contributions at leading-power.
Secondly, we consider the uncertainties due to the choice of the quark-hadron duality region in the double dispersion relation, an aspect which was not addressed in the previous study [8].
Based on these enhancements, we present the optimized and comprehensive computations for the D * (s) D (s) V and B * (s) B (s) V couplings to date.This enables us to obtain a more accurate prediction for the effective coupling λ of the HMχPT, which still exhibits significant deviation from the estimation provided by the vector meson dominance (VMD) model [40].
• We conduct a detailed investigation into the impact of SU (3) flavour (SU (3) F ) symmetry breaking effects, which stem from modifications in the light quark mass, the decay constants and mass of the vector mesons and heavy mesons, and the parameters related to the vector meson LCDAs.Our analysis indicates that the SU (3) F breaking effects play a non-negligible role in accurately determining the coupling constants g H * s HK * , g H * HsK * , g H * Hω and g H * s Hsϕ .As expected, the most significant breaking effects are observed in the coupling g H * s Hsϕ , primarily due to the involvement of the maximum number of strange quarks in this channel.
• We investigate the dispersion relation, which establishes a connection between the residue of various B (s) → V transition form factors (for V and T 1 ) at the B * (s) pole and the corresponding B * (s) B (s) V strong coupling.The couplings g B * (s) B (s) V , accompanied by systematic uncertainties, are computed from the dispersion relation using form factors obtained from two distinct LCSRs methods [28,32].
We update all input parameters utilized in LCSRs.Specifically, we employ the RunDec package [41] to adopt the four-loop evolution values for the QCD coupling and quark mass within the MS scheme.The decay constants of vector and pseudoscalar heavy-light mesons are determined using two distinct methods: QCD two-point sum rules [9,42] with the same NLO accuracy as the LCSRs, and lattice QCD (LQCD) [43,44].
The structure of this paper is organized as follows: In section 2, we introduce the conventions and definitions for the D * (s) D (s) V and B * (s) B (s) V couplings and establish the corresponding LCSRs using the double dispersion relation.In section 3, we present the calculation details for the LO and NLO double spectral density of the correlation function.Section 4 is dedicated to discussing the uncertainties arising from the quark-hadron duality ansatz for the double dispersion relation and constructing the finalized form of the LCSRs.We focus on determining the various input parameters and analyzing the numerical results in section 5.The final section concludes the work.Essential details regarding the vector meson LCDAs are provided in Appendices A and B, while Appendix C contains expressions for the double spectral densities at NLO.
2 Light-cone Sum Rules for the D * (s) D (s) V and B * (s) B (s) V couplings Our study begins with the introduction of the phenomenological effective Lagrangian [3][4][5] governing the strong interactions between heavy mesons and light vector resonances, with the convention ϵ 0123 = −1, g H * HV represents the strong couplings, and V stands for the nonet matrix of vector mesons, and The hadronic matrix element of the effective Lagrangian, describing the coupling between specific initial and final states, plays a crucial role in determining the strength of the interaction between particles in a given process, where the vector meson V (p, η) is characterized by its mass m V , 4-momentum p, and polarization vector η.These couplings are determined by the charge and flavour quantum numbers of the initial and final hadrons.In this work, we focus on the couplings s ϕ in the charm sector, as well as g B * 0 in the bottom sector.The others can be related to the aforementioned couplings through isospin symmetry, which dictates certain relationships between couplings involving particles with different isospin quantum numbers, D case : In the more general SU (3) flavour symmetry framework, specific relations exist for the D and B cases.These relations govern certain patterns or constraints on the interactions between particles within the D and B meson sectors, based on their shared properties under SU (3) transformations, D case : B case : To establish the sum rules for the coupling g H * HV , we adopt the approach outlined in Ref. [37].
This involves investigating the correlation function between the vacuum and V -meson, utilizing two distinct local interpolating currents for vector and pseudoscalar heavy mesons, denoted by j µ and j 5 .
These currents are employed to effectively probe the relevant properties and dynamics of the heavy meson sector, Here, j µ = q1 γ µ Q and j 5 = (m Q + m q 2 ) Q i γ 5 q 2 , where the heavy quark, denoted as Q, corresponds to either c or b, while q 1,2 represent one of the light quarks, u, d, or s.
Following [45], one can utilize the analytic properties of the amplitude in two distinct invariant variables: (p + q) 2 and q 2 , corresponding to the squared momenta of interpolating currents in the two channels.By incorporating the complete sets of intermediate states with quantum numbers of H and H * into these channels, we derive the double dispersion relation, It's worth noting that we have omitted the continuum part in this derivation.
To establish the LCSRs for the strong coupling constant, we begin by substituting its definition and expressing the matrix elements of interpolating currents in terms of the decay constants f H and f H * .Subsequently, we derive the double dispersion relations for the amplitude F ((p + q) 2 , q 2 ), The relevant matrix elements used to define the heavy meson decay constants are outlined as follows: The leading term in Eq. ( 10) represents the contribution from the ground-state double-pole, which involves the product of the desired strong coupling g H * HV and the corresponding decay constants.The term ρ cont (s 1 , s 2 ) captures the combined spectral density of excited and continuum states, and Σ denotes the duality region occupied by these states in the (s 1 , s 2 )-plane.Additional terms resulting from subtractions, represented by ellipses, vanish after the application of the double Borel transformation.
In QCD, the invariant function F ((p + q) 2 , q 2 ) can be calculated at momentum transfers q 2 and (p + q) 2 much smaller than the heavy quark mass m 2 Q through the light-cone OPE applied to the correlation function (8) near the light-cone x 2 ∼ 0. The result can be factorized in terms of the convolution of a hard kernel and the LCDAs of vector mesons, classified by their twist.The former represents the short-distance perturbative contributions, while the latter parameterizes the longdistance non-perturbative effects.The OPE calculations can be represented in the form of a double dispersion relation: where the involved dual spectral density ρ (OPE) (s 1 , s 2 ) refers to, where Im s 1 and Im s 2 correspond to the sequential extraction of the imaginary part of the corrleation function F (OPE) (s 1 , s 2 ) in the variable s 1 and s 2 .After performing the double Borel transformation with respect to the variables (p + q) 2 → M 2 1 and q 2 → M 2 2 , which is a technique used to enhance the convergence and suppress the contributions from higher-dimensional operators in the operator product expansion (OPE), we arrive at the sum rules for the coupling constant g H * HV , The inclusion of the integration boundary Σ, which is dual to the ground-state contribution to (10), arises from subtracting the continuum contributions using the parton-hadron duality ansatz.
3 Double spectral density of the correlation function

Double spectral density at LO
Near the light-cone1, characterized by x 2 ∼ 0, the OPE for the correlation function given in (8) remains applicable under the condition that both external momenta squared, (p + q) 2 and q 2 , are significantly below the heavy quark threshold m 2 Q .Specifically, to ensure the validity of power counting within the OPE, it suffices to have: Here, τ ≫ Λ QCD is a parameter independent of m Q .The heavy quark propagating in the correlation function is consequently highly virtual and can be expanded near the light-cone.The initial expression in Eq. ( 8) undergoes transformation to: where, S Q (x, 0) consists of the free-quark propagator and a term representing one-gluon emission, to our level of accuracy.Subsequently, we obtain the two-particle (2p) and three-particle (3p) contributions to the correlation function (16), µ,2p ((p + q) 2 , q 2 ) = ϵ µνρσ η ν p ρ q σ F (LO) 2p Upon inserting the DAs into the correlation function, we focus on the key factors discussed below, for two-particle case : for three-particle case : which reveal that for the two-particle case, k ≡ −(q + ūp), with ū = 1 − u and for the three-particle case, k ≡ −(q + ᾱp), with ᾱ = 1−α = 1−(α 1 +vα 3 ).The removal of the 1/px factors is accomplished by applying the replacement rules, as outlined below: for two-particle case : for three-particle case : where the three-particle DAs with hat are defined as: The LO results for the invariant amplitude are obtained by summing the contributions from individual twists up to twist-5 accuracy, For the three-particle case, a parallel expression is obtained by replacing u(ū) with α(ᾱ).Consequently, we arrive at the following compact form through the power expansion: where and The notations include: Our next step is to derive the LO double spectral density, which is given by, We extend the derivation of the dispersion representation within the framework of the tree-level factorization formula (21) by considering the general expression for double spectral densities governing invariant amplitudes in the two-particle case.The expression is given by, where the function ϕ 2p (u) is expressed through a Taylor expansion at u = 0 The three-particle scenario follows a framework similar to that of the two-particle case, where the coefficient c Φ3p ℓ,k and c Φ3p ℓ,k originate from the series expansion of two "effective" DAs, namely Φ 3p (α, ℓ), Φ 3p (α, ℓ), To account for the terms proportional to q 2 in Eq.( 29), we introduce supplementary functions ρ jk,3p (s 1 , s 2 ) as described by Eq.( 33): Performing a variable substitution by replacing α with u, facilitates the combination of contributions from three-particle processes with those from the previously considered two-particle contributions.With the replacement of each twist and multiplicity component in the sum (22) with its corresponding double dispersion form, we derive the LO double spectral density of the correlation function, This expression features an expansion form (32) for each term, where the coefficients c k are readily determined from the polynomial form of the DAs, as explicitly presented in Appendix B. The resulting expression (37) contributes to a new and comprehensive understanding.For instance, we illustrate the contributions to ρ (LO) from the twist-2 and twist-3 DAs, emphasizing their asymptotic behavior,

Double spectral density at NLO
To achieve next-to-leading order (NLO) precision, we express the invariant amplitude associated with the correlation function (8) in the following form, The loop computations in this context closely follow the procedures detailed in [29] for the leading twist.The one-loop diagrams depicted in Figure 2 are computed using the method of regions, initially introduced in Ref. [46].It is evident that only the hard region contributes significantly, while the collinear region yields scaleless integrals within dimensional regularization, resulting in negligible contributions.The one-loop diagrams yield the following results: ,h H Here, we introduce the definitions r 1 ≡ (p + q) 2 /m 2 Q , and µ is presented at the leading twist, The box diagram yields contributions at order O(ϵ) in dimensional regularization, and therefore does not significantly contribute.This observation is consistent with previous findings [37,47], which demonstrate the vanishing of hard-collinear factorization for the one-loop box diagram in the context of hadronic photon correction.Combining the contributions from all one-loop diagrams, the results are in agreement with those presented in a previous study [29]: Previous work [37] has demonstrated that the matrix element of the evanescent SCET operator, denoted as ⟨O E,µ ⟩ (1) , vanishes at order O(α s ) within the NDR scheme of γ 5 when considering hardcollinear factorization.This leads to the factorization formula for twist-2 contributions at next-toleading order (NLO): Now, it is imperative to advance further and obtain the NLO double spectral density: For the NLO terms, we employ asymptotic vector meson DAs.This choice is justified by observing that at LO, the non-asymptotic contributions from Gegenbauer moments in the twist-2 DA (as described in Eq.( 82) constitute only a small percentage in the LCSR, assuming typical values for the moments a 2 and a 4 .Additionally, for the K * particle, terms related to a 1 and a 3 vanish due to the u = 1/2 condition used in the calculation, effectively eliminating odd Gegenbaur terms.The additional O(α s ) factor significantly suppresses these contributions, keeping them well below the uncertainties inherent in the sum rule.

Quark-hadron duality and sum rules
After computing the double spectral density (13), we obtain the expression, where the LO contribution is provided in (37), and the NLO part consists of the twist-2 contributions, as detailed in (97).
Following Ref. [9], we adopt a specific parameterization for the boundaries: We explore three different duality regions: the triangular region where α = 1 and s * = 2s 0 , the concave region with α = 1 2 and s * = 4s 0 , and the convex region where α = 2 and s * = √ 2s 0 .Our findings reveal that the terms containing δ(s 1 − s 2 ) and its derivatives contribute equally across all duality regions.Therefore, only NLO contributions not involving the delta-function or its derivatives are affected by the choice of region.However, most of the LO and the primary part of NLO contributions come from integrating over the diagonal interval, which remains consistent across all regions.Therefore, we opt for the most practical choice: the triangular region, where s 1 + s 2 ≤ 2s 0 .
Returning to the LCSR Eq. ( 14), we then assume equal Borel parameter M 2 1 = M 2 2 = 2M 2 .This allows us to rewrite the sum rule as: Here, we define the integral over the triangular region as follows: Each DA, expressed through its Taylor expansion (32), decomposes the density ρ into individual terms.Importantly, within the context of the triangular duality region, the resulting expressions for the integrals F(M 2 , s 0 ) can be universally formulated for any generic DA.
To achieve this, a variable transformation is applied to the integration variables in (48): and their derivatives over v. Simultaneously, the exponential factor in (48) becomes independent of v.This results in the Taylor expansion of an arbitrary DA ϕ(u) reducing to its value or its derivative at u = 1/2.This yields the following simplified expression: In the subsequent analysis, F 3 stands out as the only contributor to the term in (50).
In conclusion, the LO part of the LCSR in ( 47) is derived as a linear combination of distinct contributions from individual DAs: where After comparing our results with those in Ref. [8], we find agreement in the two-particle components.
However, discrepancies arise in the three-particle counterparts, with numerical values displaying opposite signs.Upon careful examination, we identifiy a confusion in the simplification process observed in Eq. ( 25) of [8] (also referred to in Eqs.(4.10-4.16) of [29]).This confusion stems from overlooking the introduction of additional variables, as demonstrated in Eq. ( 20), when addressing three-particle expressions.Consequently, inaccuracies occurs in integrating when applying the 1/(p•x) replacement rule to the three-particle DAs.Our discovery emphasizes the critical importance of accurately managing higher-order terms and highlights the necessity for validation and verification in theoretical frameworks.
The subsequent stage of our analysis involves deriving the NLO contribution to Eq. ( 47) by incorporating the twist-2 NLO double spectral densities, ρ (tw2,NLO) (s 1 , s 2 ), as explained in Eq. ( 97) of Appendix C. The resulting expression for F (NLO) (M 2 , s 0 ) within the triangular duality region is presented below, capturing the higher-order corrections to the LCSR formulation.
accompanied by the NLO contributions of twist-2, It's important to note that the twist-2 part of the expression precisely matches that derived in [29].When transitioning to the pole-mass scheme for the heavy quark, we incorporate the terms ∆f (tw2) (σ) from (101) to (57).Additionally, we have verified the factorization-scale independence for both twist-2 terms in the LCSR (47) at O(α 2 s ) in the asymptotic limit, which further enhances our confidence in the obtained results.
The derivation of the LCSR expression (47) for the strong H * HV coupling, applicable to both B (s) and D (s) mesons with their respective b or c quark masses, is now complete for the triangular duality region.This comprehensive formulation includes both LO and NLO terms, as described in (52) and (57), respectively.The formulation is now prepared for numerical analysis.

Numerical analysis
This section is dedicated to determining various input parameters and extracting the strong couplings g D * (s) D (s) V and g B * (s) B (s) V from the established LCSRs in Eq. (47).The mass of heavy-light mesons and light vector mesons are taken from the PDG [48].The consideration of decay constants related to pseudoscalar and vector heavy-light mesons requires careful attention.We employ two distinct methods: the first utilizes LQCD values for decay constants of charmed and bottom mesons.Specifically, we incorporate N f = 2 + 1 + 1 results from [49] for heavy pseudoscalar mesons and decay constant ratios from [43] For comparison, the second employs two-point QCD sum rules to determine decay constants f H and f H * , (H = D (s) , B (s) ), as described in [9,42].This approach incorporates NLO corrections to partially cancel the perturbative effects of the LCSRs.Here, we present the decay constant values computed using the two-point QCDSRs Here, we use the RunDec Mathematica package [41] to investigate the evolution of the QCD coupling and quark mass (in the MS scheme) with four-loop precision, as outlined in Table 1.Table 1: QCD parameters used in the LCSRs and two-point QCDSRs.
Our discussion for the vector meson LCDAs primarily focuses on the LO part of the sum rules in Eq.( 52).As we explained in the last section, we need to determine the values of the DAs or their derivatives at the midpoint u = 1/2.This relies on a complete set of coefficients in the conformal expansion, which is shaped by the inherent symmetry properties governed by the renormalization group (RG) equation controlling their scale dependence [50].It differs from the LCSRs for B (s) → V and D (s) → V form factors, where vector meson DAs undergo weighting by the Borel exponent and integration over a duality interval.Additionally, since we consider the NLO for asymptotic twist-2 two-particle DAs, the renormalization of Gegenbauer coefficients should extend up to NLO as well.
This involves a complex multiplicative process, as illustrated by Eq.(84) in Appendix B.
In addition to the decay constants f ⊥,∥ V (µ) and twist-2 parameters a ⊥,∥ n (µ) detailed in Table 2, we present the Gegenbauer moments for the twist-3 and twist-4 DAs of vector mesons in Table 3.
These DAs are developed up to the NLO of conformal expansion in Ref. [51,52].The normalization and non-asymptotic coefficients at µ = 1 GeV, used as input parameters, are determined through computations using QCD sum rules (see Ref. [51][52][53] and citations therein).For the ω meson, the decay constants f ∥ ω = 197(8) MeV and f ⊥ ω = 148 (13) MeV are extracted from the provided reference [28].The remaining parameters for this meson are consistent with those employed for the ρ meson.Subsequently, we specify the renormalization scale of the quark-gluon coupling and quark masses, the factorization scale, the Borel parameter, and the quark-hadron duality thresholds of the LCSRs in Eq. (47).These parameters play crucial role in determining the behavior of the sum rules and extracting the strong couplings.We choose these parameters to ensure that the sum rules exhibit strong convergence and stability, while also minimizing theoretical uncertainties in the final outcomes.5) 220( 5) 215( 5) 215( 5) and the effective threshold s 0 in LCSRs follows established criteria aimed at ensuring the smallness of subleading power contributions in the light-cone OPE while simultaneously suppressing contributions from excited states.These criteria typically involve considerations of convergence and stability in the sum rule framework, and in our analysis, we adhere to them to ensure the reliability and accuracy of our results.These parameters, presented in Table 4, are consistent with prior choices in LCSRs analyses of heavy-to-light decay form factors and couplings [9,54,55].Furthermore, we choose the renormalization scale µ equal to the factorization scale used in the OPE, which is approximately To maintain the perturbative expansion of the correlation function in α s within reasonable bounds, the NLO twist-2 and higher power terms are constrained to be no more than 30% of LO counterparts.Consequently, we adopt default values of µ c = 1.5 GeV for charm and µ b = 3 GeV for bottom and investigate the range of variation within specified intervals (1.0 ∼ 3.0)GeV for charm and (2.5 ∼ 4.5)GeV for bottom to quantify the associated uncertainties.This choice is motivated by considerations of consistency within the theoretical framework and has been adopted in previous analyses within the field [25,26,[54][55][56][57][58].By aligning these scales, we ensure coherence in our calculations and facilitate comparisons with prior studies.decay constant products, respond to changes in the Borel parameter M 2 and the scale µ. Figure 3 depicts the dependence on the Borel parameter M 2 , illustrating a nearly flat behavior within the range considered, indicating minimal sensitivity to variations in M 2 .This observation suggests that the extracted quantities are robust against changes in the Borel parameter, contributing to the reliability of our results.In Figure 4, we observe the dependence of individual contributions on the scale µ.Remarkably, contributions from twist-2 at both LO and NLO display significant but opposite dependencies on µ, almost canceling each other.While other higher-twist contributions exhibit minimal dependence, the overall outcome indicates only slight sensitivity to the scale µ.

Moreover, it's valuable to analyze how
This behavior suggests a delicate balance between different contributions in the LCSRs framework, highlighting the importance of considering various factors when interpreting the results.
To ensure the efficiency of power corrections, we illustrate the contributions of each term in the power correction series, taking the heaviest ϕ meson as an example.From the  M 2 (GeV 2 ) 5.5 (4.5 -6.5) Table 4: The renormalization scale µ, Borel parameter M 2 and M 2 , along with duality threshold s 0 and s0 ( s0 * ), are employed in the LCSR and the two-point QCDSRs to determine the decay constants of H (H * ) mesons for both charmed and bottom particles.
for the D meson, despite the significant value of δ ϕ , the contributions generally decrease in order of power.Contributions at δ 1 ϕ and δ 2 ϕ , compared to δ 0 ϕ , are reduced by an order of magnitude, although δ 1 ϕ and δ 2 ϕ exhibit similar magnitudes.Likewise, δ 3 ϕ and δ 4 ϕ , relative to δ 2 ϕ and δ 3 ϕ , are also suppressed by an order of magnitude, with δ 2 ϕ and δ 3 ϕ displaying comparable magnitudes.This pattern can be explained by noting that in Eq.( 53), higher-power contributions are accompanied by a reduction in the ratio m 2 c /M 2 .In contrast to the D meson case, although the δ ϕ for the B channel is much smaller, the ratio m 2 b /M 2 played a compensatory role, this leads to a similar power expansion behaviour for the B * s B s ϕ coupling as the D * s D s ϕ.We now explore the numerical effects of perturbative QCD and distinct higher-twist corrections  on the products of strong couplings and decay constants computed in Eq. (47).The outcomes of individual twist and NLO effects are summarized in Table 6, along with relevant data uncertainties.
As observed in Table 6, the NLO corrections of twist-2 magnitude exhibit numerical equivalence with the LO contributions of twist-3, constituting roughly 20% of the leading twist for both charm and bottom scenarios.However, it's noteworthy that a sign difference is observed in NLO corrections of twist-2 between charm and bottom.This sign difference may arise from the difference in the scale ratio squared µ 2 /m 2 Q between charm and bottom systems, where µ 2 c /m 2 c >1 for charm and µ 2 b /m 2 b <1 for bottom, leading to distinct behaviors in the higher-order corrections.Additionally, our analysis indicates that the twist-5 contribution in the OPE of the correlation function has negligible impact, as validated through numerical assessments.This observation suggests that higher-order contributions beyond twist-4 have minimal influence on the computed results.Furthermore, the total uncertainty reported in Table 6 should involve variations stemming from changes in the duality region within the NLO aspect.Specifically, for the charm scenario, this variation amounts to around 20%, while for the bottom scenario, it stands at approximately 15%.These variations are reflected in the parameter ∆α, as depicted in Tables 8 and 9, where α is defined in Eq. ( 46).It's suggested that these uncertainties, to some extent, address the "systematic error" inherent in the quark-hadron duality ansatz applied in LCSRs.
In our final analysis, we normalize the couplings using heavy meson decay constants by two methods: the two-point QCDSRs in Eq.( 60) and LQCD results in Eq.( 59).Table 7 highlights the main findings, while Tables 8 and 9 provide more details on the parameters for each set of decay constants.In the D meson sector, the main sources of uncertainty lie in the Gegenbauer moment a ⊥ 2 and tensor decay constants f ⊥ V derived from LQCD.When employing two-point sum rules to calculate decay constants, an additional uncertainty arises from the decay constant itself, which is comparable in magnitude.In the B meson sector, the uncertainties related to the Borel parameter  and threshold somehow increase relative to those in the D meson sector.
The estimation for the D * (s) D (s) V and B * (s) B (s) V strong couplings from LCSRs yields uncertainty ranges of 13% to 18% and 10% to 18%, respectively, when using heavy-meson decay constants from the two-point QCDSRs.These uncertainty intervals are calculated based on variations in the input parameters and methodologies used in the LCSR analysis.While the uncertainty intervals for g D * (s) D (s) V remain consistent regardless of the decay constant choice, only slight agreement is observed in the intervals for g B * (s) B (s) V .The difference in the central values of the B * (s) B (s) V couplings mainly arise from a 20% discrepancy between the lattice-QCD value of f B * (s) and the central value predicted by the two-point sum rule.For consistency, we use the two-point QCDSRs at NLO, with more precise sum rules at NNLO [42], yielding results consistent with LQCD calculations within uncertainties.This ensures that our analysis incorporates the latest advancements in theoretical frameworks and maintains consistency with the most accurate calculations available.
To investigate the SU (3) F symmetry breaking effects, we introduce ∆r H * HV as the parameter characterizing the size of the breaking effects, where g D * Dρ and g B * Bρ are defined in Eq.( 4) and ( 5), respectively.Additionally, takes the ω meson, otherwise, it's set to 1.
In Table 10, we detail the relative magnitudes of SU (3) F symmetry breaking effects for the various strong couplings.Interestingly, both ∆r D * Dω and ∆r B * Bω remain constant at −10%.This stability is attributed to the interplay of specific factors such as f ⊥ ω , f ∥ ω , and m ω , which decrease at rates of 8%, 1.9%, and 0.1%, respectively.The other parameters remain consistent with those associated with the ρ meson, resulting in a cancellation of terms in both the numerator and denominator.Shifting focus to the 9.8% breaking observed in ∆r D * DsK * , analysis using the decay constants from two-point   (47) in units of [GeV −1 ].The negligible uncertainties arising from variations in the remaining input parameters have been considered but are not explicitly enumerated here, but are already taken into account in the determinations of the total errors.QCDSRs reveals four primary sources of this deviation.Firstly, a substantial 12% deviation arises from the LCDAs of the K * meson (8%) and the effect of the light quark mass (4%).Secondly, the decay constants of the D s meson contribute approximately −5%.Thirdly, the decay constant f ⊥ K * explains 9% of the observed breaking.Lastly, the mass of the D s meson contributes roughly −6%.As for ∆r D * s DK * , the contribution from the D s meson's decay constant f D * s is around −12.9%, larger than that of D s .Additionally, a portion of this difference can be attributed to the mass of the D * s meson, resulting in a decrease of about 0.27%.This results in an overall breaking of 7.6%.Substituting decay constants calculated via LQCD for those obtained from two-point QCDSRs notably alters the quantification of symmetry breaking effects, primarily reducing the attributed deviation from decay constants.It's crucial to highlight that the breaking effects observed in ϕ, as determined through LQCD, tends to exceed that calculated using two-point QCDSRs (taking their absolute magnitudes into account).This discrepancy arises from the substantial difference in the decay constant f Ds obtained from the two methods, which causes the negative breaking effects up to 15% from LQCD , larger than 5% from the two-point QCDSRs.However, for f D * s , the situation is reversed but less difference, −13% from two-point QCDSRs compared to −15% from LQCD.Consequently, the   When comparing our numerical outcomes with those of Ref. [8], depicted in Tables 11 and   12, we observed an enhancement of the reported values 3.80 +0.59 −0.45 GeV −1 and 3.89 +0.52 −0.48 GeV −1 for g D * Dρ and g B * Bρ respectively.The primary reason for this difference lies in variations in our chosen input parameters such as scale µ, Borel parameter M 2 , and effective threshold s 0 .The impact of the scale µ is especially pronounced in the twist-2 LO and NLO aspects, although its overall effect decreases.Shifting the scale in the bottom sector from µ = 3 GeV to µ = m b leads to an opposite sign for the NLO terms.Moreover, their specified values of M 2 and s 0 differ from our using. 2We 2 The authors in [8] specified M 2 = 4.5 ± 1.0 GeV 2 for charm and M 2 = 18 ± 3.0 GeV 2 for bottom sector, with differing by a factor of 2. Their choices of s0 = 6.0 ± 0.5 GeV 2 for charm and s0 = 34.0 ± 1.0 GeV 2 for bottom sector contrast with our selections of 7.0 ± 0.5 GeV 2 and 37.5 ± 2.5 GeV  V (GeV −1 ) from several methods.also extend our analysis to include sub-subleading power corrections ranging from twist-3 to twist-5 (considering δ 3 V and δ 4 V ) in Eq.( 53), which were absent in [8].Additionally, our treatment of the threeparticle contributions deviates from that in [8] due to our modifications to the three-particle DAs, resulting in a reversal of signs, although their influence remains negligible compared to the leading twist contribution.Furthermore, it is noteworthy that the LCSR studies conducted in [17,19] only focused on leading-order calculations, while the significant compensatory effects from the twist-2 NLO calculations were absent.As mentioned previously, in the case of charm mesons, the NLO effects at the leading twist effectively counterbalance the corrections from higher-twist contributions.
Consequently, we expect our results to align closely with those from earlier sum rule investigations, within the errors.However, in the case of bottom mesons, the opposite sign of the NLO effects enhance the impact of the leading twist. respectively.

Method
In addition, we introduce an alternative method to determine the B * (s) B (s) V coupling.This method, as extensively described in [8,9], utilizes the hadronic dispersion relation for the B (s) → V vector(tensor) form factor as a foundational element: This relation features the vector-meson B * (s) pole, just below the kinematic threshold q 2 = (m B (s) + m V ) 2 .The residue of the B * pole represents the product of the B * (s) B (s) V coupling and the B * (⊥) (s) decay constant.Multiplying both sides of Eq.( 63) by the denominator of the pole term and taking the limit , we avoid the integral over the hadronic spectral density, leading to: Employing Eq.( 65) and assuming f T B * = f B * (s) , we obtain various numerical values for the coupling g B * (s) B (s) V using the B (s) → V form factors from two distinct LCSRs [28,32].These results are summarized in Table 13.
Our study reveals a discrepancy between our central value and the projections from the transition form factors, which can be attributed to two main factors.Firstly, different theoretical calculations of form factors yield different results, leading to inconsistent conclusions [28,32].In Ref. [32], our Method This work V (q 2 )[28] V (q 2 )[32] T 1 (q 2 )[28] T 1 (q V in units of GeV −1 extracted from the residue of the transition form factors V and T 1 in the framework of LCSR fit, compared with our obtained values.results are consistent with their theoretical predictions when considering uncertainties.This consistency is attributed to the significant magnitude of the uncertainties.However, in Ref. [28], notable deviations persist even after accounting for uncertainties.Moreover, there might be an underestimation of uncertainties stemming from the parameterizations of form factors to address unphysical singularities.Secondly, assuming the validity of these findings, the Borel suppression in the double dispersion relation may not effectively eliminate contributions from isolated excitations, as discussed in Ref. [39].In sum rule calculations, it is common practice to include contributions from higher states in the perturbative estimation of a continuum, often without explicit consideration of isolated excited states.This approach is supported by the stability criterion being met adequately without requiring substantial excitation contributions.However, Ref. [39] proposes a different viewpoint, suggesting that neglecting explicit radial excitation contributions could be the root of the discrepancy.
While acknowledging the Borel suppression effect, they argue that in certain cases, such as the one under consideration, the suppression might be insufficient, leading to the need for additional considerations.Incorporating explicit radial excitation contributions to the hadronic side of the LCSRs provides a plausible explanation for the failure of the LCSR prediction.The relatively weak exponential suppression due to the Borel transformation in the LCSR framework significantly amplifies the impact of radial excitations without necessitating excessively large couplings.This suppression is significantly weaker compared to standard sum rules, primarily due to the additional factor in the exponent and the larger mean value of M 2 within the allowed range.Overall, according to Ref. [39], assuming partial cancellation between contributions of excited and ground states in the dispersion relation may lead to an increased LCSR coupling.
When extending Eq.( 65) to include transition form factors from D to V , we can derive the coupling for D * DV as well.By employing the D to ρ form factor detailed in the Ref. [59], our calculations yielded couplings g D * Dρ of 8.91 +0.77  −0.55 GeV −1 using the heavy meson decay constants from LQCD, and 8.29 +0.71  −0.51 GeV −1 using the decay constants from two-point QCDSRs , which exceed the values of 3.40 +0.49−0.43 GeV −1 and 3.53 +0.61 −0.57GeV −1 obtained from LCSRs a lot, mirroring the trend observed in B → ρ transitions.Nevertheless, since the Ref. [59] only provides uncertainty estimates for V (0) and not for the fitting parameters a 1 and a 2 , the resulting uncertainty remains relatively small.Finally, we also compare our results with other model-dependent studies concerning chiral and heavy quark limits.The coupling parameter λ characterizes the H * Hρ interaction, and it's related to g H * Hρ through the expression, where the parameter g V = m ρ /f π (m ρ is the ρ meson mass and f π denotes the pion decay constant) is determined from the first and second KSRF relation [3,60,61].Our calculations yield λ = 0.21 ± 0.03 GeV −1 with decay constants obtained from two-point sum rule and λ = 0.20 ± 0.03 GeV −1 from LQCD-calculated decay constants at the leading power, based on our determination of g B * Bρ .
These values are compatible with those obtained from Ref. [8,18,19] but significantly smaller than those predicted by the VMD model [40], where λ is estimated to be 0.56 GeV −1 , and the constituent quark model [62], which suggests a value of 0.47 GeV −1 .One possible reason for this discrepancy could be the likelihood of larger uncertainties in the model predictions, favoring the reliability of the LCSRs method.However, uncertainties in the determination of NLO corrections to higher-twists (such as twist-3 NLO) and sub-sub-leading power contributions remain unknown factors contributing to the observed deviation.Advancements in computational techniques, including refinements in sum rule methods and modeling, may offer valuable insights into resolving these discrepancies in the future.

Conclusion
In this paper, we have conducted a comprehensive investigation into the improved calculation of the strong couplings g D * (s) D (s) V and g B * (s) B (s) V , where V = ρ, ω, K * and ϕ, within the framework of LCSR.
This framework is originally developed in Refs.[8,9,29], which primarily focus on the calculation of g D * Dρ and g B * Bρ .
Our study built upon the OPE method, focusing on the correlation function Eq. ( 8) and incorporating vector meson LCDAs with increasing twist.We analytically derived the double spectral representation of the QCD factorization formula, employing the parton-hadron duality ansatz and the double Borel transformation.This allows us to establish LCSRs for both leading and sub-leading contributions, encompassing two-particle DAs up to twist-5 and three-particle DAs up to twist-4.
Exploring the newly determined double spectral densities for subleading twist contributions allows us to perform conduct analytical continuum subtractions.This facilitates the development of highertwist sum rules on the light-cone, improving the precision and scope of our analysis.Furthermore, we derived concise analytical expressions for NLO terms within the MS scheme for heavy quark mass.
This ensures the applicability of our findings to both charmed and bottom mesons.Additionally, we have updated all input parameters for the LCSRs, focusing particularly on the Borel parameter M 2 , scale µ, and effective threshold s 0 .These modifications are based on newly cited data from Ref. [9], reflecting improvements in our understanding compared to the findings presented in Ref. [8].
Examining the numerically calculated LCSRs for the strong couplings g H * HV , we observed that the dominant contribution arises from twist-2.In addition, we found that the NLO corrections closely match the LO contributions at twist-three, constituting nearly 20% of the leading twist.Notably, our analysis revealed opposite signs for the NLO corrections between charm-meson and bottom-meson couplings.Our analysis also indicated that the contributions from subleading twists to charm-meson couplings, particularly from two-particle and three-particle LCDAs, are more significant than those for bottom-meson couplings.
When comparing our optimized LCSRs predictions with the previous studies [8,17,19,21], we found satisfactory agreement in the values obtained for the H * HV couplings within theoretical uncertainties.Have these couplings in hand, we extracted the effective coupling λ in the HMχPT Lagrangian, it is compatible with those obtained in Ref. [8,18,19] but significantly smaller than those predicted by the VMD model [40] and the constituent quark model [62].Furthermore, we fulfilled a detailed exploration to the SU (3) flavour symmetry breaking effects by comparing the H * HV (for The definitions of f ∥ V and f ⊥ V are given by: The three-particle chiral-even DAs [51] are defined by the following matrix elements and the three-particle chiral-odd DAs [51] can be defined as where the formulas for the RG factors E ∥,⊥ n,(N)LO (µ, µ 0 ) and off-diagonal mixing coefficients d k n (µ, µ 0 ) within the MS scheme are explicitly provided in references [47,50].The evolution of a ⊥ n (µ) is closely linked to f ⊥ V (µ) because of the definition of the conformal operator.
The precise definition of ⟨⟨Q (i) ⟩⟩ is given in Ref. [51].ω ⊥ 3K * is a twist-3 parameter.Existing numerical determinations of these parameters from QCD sum rules are far from being precise, so we decide to estimate them using the renormalon model of Ref. [66] instead.
The renormalon model implies the following estimates for the G-conserving twist-4 parameters: • the twist-5 two-particle DAs(asymptotic form) [28]: Recall that ū = 1 − u .The expression discussed is applicable to a V = (q 2 q1 ) meson.In the case of V = (q 1 q2 ), one must replace u with 1 − u and exchange the values of m q 2 and m q 1 .

C Double spectral density and integral at NLO
The expression for the twist-2 component of the Next-to-Leading Order (NLO) double spectral density in the sum rule (47) is obtained from the calculations presented in Ref. [29].The analytical expression for the double spectral density can be collected: with the introduction of two new dimensionless variables carried out through the following definitions:

Figure 3 :
Figure 3: Borel parameter dependence of the products f D* (s) f D (s) g D * (s) D (s) V and f B * (s) f B (s) g B * (s) B (s) Vcalculated from LCSR.Displayed are the total values and separate contributions from twist-2 to twist-5 of the vector meson (ρ, K * , ϕ) DAs, considering central values for all the other input parameters.

Figure 4 :
Figure 4: Scale dependence of the products f D* (s) f D (s) g D * (s) D (s) V and f B * (s) f B (s) g B * (s) B (s)V calculated from LCSR.Displayed are the total values and separate contributions from twist-2 to twist-5 of the vector meson(ρ, K * , ϕ) DAs, considering central values for all other input parameters.
discrepancy attributed to LQCD calculated decay constants amounts to −25%, compared to −13% from two-point QCDSRs.Finally, after accounting for additional factors contributing around 5% to the observed symmetry breaking effects, we arrive at our comprehensive analysis of SU (3) flavour symmetry deviations.The situation with SU (3) symmetry breaking in the B sector resembles that in the D sector, albeit with milder decrease.

Table 2 :
[53]y constants and twist-2 hadronic parameters, initially referenced at µ = 1 GeV[53], are then scaled to µ = 1.5 GeV and µ = 3 GeV using the evolution equations outlined in Appendix B.Since we use the finite heavy quark mass, these parameters are naturally different for the couplings with charmed and bottom mesons.However, the heavy-quark spin symmetry enables certain scale and the Borel parameter to be aligned in the H and H * channels.Detailed default values and intervals for these parameters are provided in Table4, which offers a comprehensive overview of the specific values and intervals for the parameters mentioned above.The determination of the Borel parameter M 2

Table 6 :
Bs g B * BsK * 0.106 +0.034 Summary of the contribution from individual twist and NLO effects on the products of the strong couplings g H * f * HV and decay constants presented with both central values and associated data errors in units of [GeV].

Table 7 :
LCSR results for the strong couplings g H * HV for the two methods of dividing out the decay constants in units of [GeV −1 ]. 53

Table 8 :
Summary of the individual uncertainties for the strong couplings g H * HV predicted from the LCSR

Table 9 :
Summary of the individual theory uncertainties for the strong couplings g H (47)predicted from the LCSR(47)in units of [GeV−1 ].The negligible uncertainties arising from variations in the remaining input parameters have been considered but are not explicitly enumerated here, but are already taken into account in the determinations of the total errors.

Table 10 :
2 , Dω ∆r D * s DK * ∆r D * DsK * ∆r D * Summary of the extent of SU (3) breaking effects as determined by LCSR, with ∆r H * HV

Table 11 :
Numerical values of coupling g D * (s) D (s)

Table 12 :
Numerical values of coupling g B * (s)

Table 13 :
Coupling constant g B * (s) B (s) [63]7,19,21] with g D * Dρ and g B * Bρ .The analysis revealed breaking effects ranging from −25% to −1.8% when employing the decay constants of heavy mesons calculated from LQCD, and from −13% to +9.8% when utilizing the decay constants obtained from two-point QCDSRs.Where, both g H * Hω and g H * s Hsϕ demonstrate negative symmetry breaking effects, however the breaking effects of g H * HsK * and g H * s HK * exhibit opposite signs for the heavy meson decay constants computed from the two methods.The most obvious SU (3) F breaking effects were observed in the coupling g D * s Dsϕ which reach to −25%, this is primarily due to the involvement of the maximum number of strange quarks in this channel.Moreover, we performed a comprehensive investigation of g B * (s) B (s) V and g D * Dρ derived from B → V[28,32]and D → ρ transition form factors by using dispersion relation.Which demonstrated an obvious deviation from our findings in LCSRs and the results obtained in the previous studies[8,17,19,21].Developing the LCSRs for the strong coupling g H * HV beyond our current investigation offers several promising directions.Firstly, it's crucial to calculate NLO QCD corrections to the higher-twist contributions of the vector meson DAs, this is important for clarifying the discrepancies observed in g H * HV between the direct LCSRs calculations and the indirect predictions from the H → V transition form factors. Secondly, updating the non-perturbative parameters in the conformal expansion of vector meson DAs holds significant phenomenological importance and can enhance the accuracy of our LCSRs predictions.As noticed in[9], a twist-2 model for pion LCDA was proposed from the LQCD data[63], which has led to a notable improvement in predictions for B * Bπ.Another potential solution to this problem is to incorporate high excited heavy-meson states in the hadronic components of the sum rules, as suggested in the case of D [39][39].