Electromagnetic properties of heavy-light mesons

TL;DR

This paper computes electromagnetic form factors for light and heavy-light pseudoscalar mesons using the Bethe-Salpeter framework, achieving good agreement with experimental and lattice QCD data.

Contribution

It introduces a flavour-dependent effective charge in the Bethe-Salpeter approach, improving the description of quark-fluon interactions and mass differences.

Findings

Electromagnetic form factors match experimental data for pion and kaon.

Predicted charge radii of heavy-light mesons agree with lattice QCD.

Nonperturbative quark-photon vertex is effectively determined.

Abstract

Within the Bethe-Salpeter framework, we present a computation of space-like electromagnetic form factors for pseudoscalar mesons, including light and heavy-light systems. Our approach employs a flavour-dependent variation of the standard Taylor effective charge, which contains key contributions from the quark-gluon vertices. This effective interaction is a common ingredient of all relevant dynamical equations, and accommodates the crucial mass differences between the various quark flavours. Particular attention is paid to the nonperturbative determination of the quark-photon vertex. The computed electromagnetic form factors for the pion and the kaon mesons show excellent agreement with experimental determinations. In addition, the predictions for the charge radii of heavy-light systems are in overall good agreement with lattice QCD.

Equations

$$J^{\mu }(p_{\text{av}},q)=2p_{\text{av}}^{\mu }{F}_{\text{s}}(q^2)⟶F_{\text{s}}(q^2)=\frac{J\cdot p_{\text{av}}}{2p_{\text{av}}^2},$$

$$J^{\mu }(p_{\text{av}},q)={\int }_k\bar{\mathcal{A}}(k_f,p_f)S({\ell }_{+}){\Gamma }^{\mu }(k_{+},q)S({\ell }_{-})\mathcal{A}(k_i,p_i)S(k_{-}),$$

$${\int }_k:=(2\pi {)}^{-4}{\int }_{-\infty }^{+\infty }{\mathrm{d}}^{4}k,$$

$$k_{+}=k+\eta p_{\text{av}},k_{-}=k-(1-\eta )p_{\text{av}},{\ell }_{+}=k_{+}+\frac{q}{2},$$

$${\ell }_{-}=k_{+}-\frac{q}{2},s_{+}=k_{-}+\frac{q}{2},s_{-}=k_{-}-\frac{q}{2}.$$

$$F_{\text{s}}(q^2)=e_q{F}_q(q^2)+e_{\overline{q}}{F}_{\overline{q}}(q^2),$$

$${\mathcal{I}}_{f{f}^{\prime }}(q^2)={\tilde{\alpha }}_T(q^2)A_f(q^2)A_{f^{\prime }}(q^2),$$

$$D_{f{f}^{\prime }}^{\mu \nu }(k):=4\pi {\Delta }_0^{\mu \nu }(k){\mathcal{I}}_{f{f}^{\prime }}(k^2),$$

$${\Delta }_0^{\mu \nu }(k)=\left({\delta }^{\mu \nu }-\frac{k^{\mu }{k}^{\nu }}{k^2}\right)\frac{1}{k^2}.$$

$$S_f^{-1}(p)=i/p{A}_f(p^2)+B_f(p^2),$$

$$S_f^{-1}(p)=Z_2(i/p+m_R)+\underset{{\Sigma }_f(p^2)}{\underbrace{C_F{\int }_k{\gamma }_{\alpha }{S}_f(k+p){\gamma }_{\beta }{D}_{ff}^{\alpha \beta }(k)}},$$

$${\mathcal{A}}_{f{f}^{\prime }}(p_1,p_2)=-{\int }_k{\gamma }_{\mu }S(k_1){\mathcal{A}}_{f{f}^{\prime }}(k_1,k_2)S(k_2){\gamma }_{\nu }{D}_{f{f}^{\prime }}^{\mu \nu }(k),$$

$${\mathcal{A}}_{f{f}^{\prime }}(p,P)=\left({\chi }_1^{f{f}^{\prime }}+i{\chi }_2^{f{f}^{\prime }}/P+i{\chi }_3^{f{f}^{\prime }}/p(p\cdot P)+{\chi }_4^{f{f}^{\prime }}[/p,/P]\right){\gamma }_5,$$

$${\Gamma }_{\mu }^f(p_1,p_2)=Z_2{\gamma }_{\mu }-C_F{\int }_k{\gamma }_{\alpha }{S}_f(k_1){\Gamma }_{\mu }^f(k_1,k_2)S_f(k_2){\gamma }_{\beta }{D}_{ff}^{\alpha \beta }(k).$$

$$P^{\mu }{\Gamma }_{5\mu }^{f{f}^{\prime }}(p_1,p_2)=S_f^{-1}(p_1)i{\gamma }_5+i{\gamma }_5{S}_{f^{\prime }}^{-1}(p_2)-i(m_f+m_{f^{\prime }}){\Gamma }_5^{f{f}^{\prime }}(p_1,p_2),$$

$${\Gamma }_{5\mu }^{f{f}^{\prime }}(p_1,p_2)={\gamma }_5{\gamma }_{\mu }\underset{G_{5\mu }^{f{f}^{\prime }}(p_1,p_2)}{\underbrace{-{\int }_k{\gamma }_{\alpha }{S}_{f^{\prime }}(k_1){\Gamma }_{5\mu }^{f{f}^{\prime }}(k_1,k_2)S_f(k_2){\gamma }_{\beta }{D}_{f{f}^{\prime }}^{\alpha \beta }(k)}},$$

$${\Gamma }_5^{f{f}^{\prime }}(p_1,p_2)={\gamma }_5\underset{G_5^{f{f}^{\prime }}(p_1,p_2)}{\underbrace{-{\int }_k{\gamma }_{\alpha }{S}_{f^{\prime }}(k_1){\Gamma }_5^{f{f}^{\prime }}(k_1,k_2)S_f(k_2){\gamma }_{\beta }{D}_{f{f}^{\prime }}^{\alpha \beta }(k)}},$$

$$P^{\mu }{\gamma }_5{\gamma }_{\mu }=S_{0,f}^{-1}(p_1)i{\gamma }_5+i{\gamma }_5{S}_{0,f^{\prime }}^{-1}(p)-2i(m_f+m_{f^{\prime }}){\Gamma }_{0,5}^{f{f}^{\prime }}(p_1,p_2),$$

$$P^{\mu }{G}_{5\mu }^{ff}(p_1,p_2)=$$

$$=$$

$$+2i{m}_f{\int }_k{\gamma }_{\alpha }{S}_{f^{\prime }}(k_1){\Gamma }_5^{f{f}^{\prime }}(k_1,k_2)S_f(k_2){\gamma }_{\beta }{D}_{ff}^{\alpha \beta }(k),$$

$$=$$

$$P^{\mu }{G}_{5\mu }^{f{f}^{\prime }}(p_1,p_2)={\tilde{\Sigma }}_f(p_1)i{\gamma }_5+i{\gamma }_5{\tilde{\Sigma }}_{f^{\prime }}(p_2)-i(m_f+m_{f^{\prime }})G_5^{f{f}^{\prime }}(p_1,p_2),$$

$${\tilde{\Sigma }}_f(p_1)=C_F{\int }_k{\gamma }_{\mu }{S}_f(k_1){\gamma }_{\nu }{D}_{f{f}^{\prime }}^{\mu \nu }(k),{\tilde{\Sigma }}_f(p_2)=C_F{\int }_k{\gamma }_{\mu }{S}_f(k_2){\gamma }_{\nu }{D}_{f{f}^{\prime }}^{\mu \nu }(k).$$

$$(\delta A{)}_{u{f}^{\prime }}(p^2)=$$

$$(\delta \mathcal{M}{)}_{u{f}^{\prime }}(p^2)=$$

$$q^{\mu }{\Gamma }_{\mu }^f(p_1,p_2)=i{S}_f^{-1}(p_2)-i{S}_f^{-1}(p_1).$$

$$q^{\mu }{\Gamma }_{\mu }^f(p_1,p_2)=$$

$$=$$

$$i\left[i/p_1+m_f+C_F{\int }_k{\gamma }_{\alpha }{S}_f(k_1){\gamma }_{\beta }{D}_{ff}^{\alpha \beta }(k)\right],$$

$${\Gamma }^{\mu }(p_1,p_2)={\Gamma }_{\mathrm{B}\mathrm{C}}^{\mu }(p_1,p_2)+{\Gamma }_{\mathrm{T}}^{\mu }(p_1,p_2),$$

$${\Gamma }_{\mathrm{B}\mathrm{C}}^{\mu }(q,u)={\lambda }_1{\gamma }^{\mu }+2{\lambda }_2{u}^{\mu }/u+2i{\lambda }_3{u}^{\mu }+i{\lambda }_4[{\gamma }^{\mu },/u],$$

$${\lambda }_1=\frac{A(p_1^2)+A(p_2^2)}{2},{\lambda }_2=\frac{A(p_1^2)-A(p_2^2)}{p_1^2-p_2^2},{\lambda }_3=\frac{B(p_1^2)-B(p_2^2)}{p_1^2-p_2^2},{\lambda }_4=0.$$

$${\Gamma }_{\mathrm{T}}^{\mu }(q,u)=\sum _{i=1}^8{h}_i(q,u){\tau }_i^{\mu },$$

Full text

I Introduction

The electromagnetic form factors (EFF) of light mesons, in particular those of the pion and the kaon, are of key importance for the understanding of the internal structure and dynamics of hadrons in terms of their quark and gluon constituents. Over the past decades, both experimental measurements and theoretical investigations have made significant progress in exploring these form factors across various momentum-transfer regimes.

From the theoretical point of view, hadronic form factors have been extensively studied within the framework of the Bethe-Salpeter equations (BSEs), a covariant approach for describing bound states in quantum field theory Salpeter and Bethe (1951); Gell-Mann and Low (1951); Bethe and Salpeter (1957); Nakanishi (1969); Jain and Munczek (1993); Munczek (1995). This framework has been adopted in numerous studies of the EFF of light mesons Maris and Tandy (2000); Cloet et al. (2009); El-Bennich et al. (2008); Sanchis-Alepuz et al. (2013); Weil et al. (2017); Eichmann et al. (2019, 2020); Miramontes et al. (2021, 2022a); Xu et al. (2024); Eichmann et al. (2017); Raya et al. (2016); Xu (2024); Miramontes et al. (2025a); Xu et al. (2025a); Miramontes et al. (2025b) and spectroscopy Eichmann et al. (2016); Maris and Tandy (1999); Alkofer et al. (2002); Nicmorus et al. (2009); Hilger et al. (2015a, b); Mojica et al. (2017); Serna et al. (2017); Miramontes et al. (2022b); Gao et al. (2025); Miramontes et al. (2025c), achieving progressively higher levels of precision and sophistication. On the lattice QCD side, a plethora of calculations is available for light mesons Brömmel et al. (2007); Frezzotti et al. (2009); Boyle et al. (2008); Aoki et al. (2009); Brandt et al. (2013); Meyer (2011); Erben et al. (2020).

On the experimental side, the precise determination of the pion and kaon EFF remains a priority for contemporary and future hadronic physics programs. Upcoming experiments at major facilities, such as Jefferson Lab Accardi et al. (2024), the U.S. Electron-Ion Collider (EIC) Abdul Khalek et al. (2022), and the Electron-Ion Collider in China (EicC) Anderle et al. (2021), are expected to provide accurate data on light meson form factors, particularly in the high-$Q^{2}$ region.

Within the BSE formalism, EFF are calculated using the so-called “impulse approximation”, where the interaction of the photon with the quark-antiquark system is mediated via the dressed quark-photon vertex and the Bethe-Salpeter (BS) amplitude of the meson. An important ingredient of this calculation is the quark-photon vertex itself, which we compute here using its full general structure. In particular, we include both the longitudinal components satisfying the Ward-Takahashi identity (WTI) for electromagnetic current conservation, as well as the transverse components that capture dynamical effects beyond the bare vertex. This comprehensive treatment is essential for the accurate description of the EFF, particularly in heavy-light systems, where flavour asymmetries can enhance the impact of transverse contributions.

While the light sector has been thoroughly analyzed Maris and Tandy (2000); Cloet et al. (2009); Eichmann et al. (2008); El-Bennich et al. (2008); Sanchis-Alepuz et al. (2013); Weil et al. (2017); Eichmann et al. (2019, 2020); Miramontes et al. (2021, 2022a); Xu et al. (2024); Eichmann et al. (2017); Raya et al. (2016); Xu (2024); Miramontes et al. (2025a, b), the inclusion of heavy-light systems, such as the $D$, $B$, and their strange counterparts ($D_{s}$, $B_{s}$), presents unique challenges, including the treatment of asymmetric quark masses and singularities in the quark propagators.

In this work, we employ the BSE formalism to compute the EFF of both light and heavy-light pseudoscalar mesons. To account for the nontrivial flavour dynamics in these systems, we implement the flavour-dependent interaction previously developed in Gao et al. (2025). In this approach, the interaction strength acquires a natural dependence on the quark flavour via the dynamically generated quark wave function. By incorporating this flavour sensitivity directly into the BSE and gap equation, we can investigate how mass differences between the constituent quarks influence the momentum dependence of the EFF.

The article is organized as follows. In Section II we present the general framework that will be employed throughout this work. In Section III we present a quantitative study of the amount of symmetry violation stemming from the use of the flavour-dependent interaction in the quark gap equation. In Section IV we carry out a detailed analysis of the quark-photon vertex, and present results for all of its form factors. Section V contains the main results of this work, namely the EFF and charge radii, together with an extensive comparison with experiment, lattice QCD, and a variety of approaches. Finally, in Section VI we discuss our conclusions.

II General framework

In this section we review the standard framework employed for the computation of the EFF, and discuss the key dynamical equations and their main ingredients.

II.1 Electromagnetic form factors in the impulse approximation

Mesonic form factors are extracted from the physical amplitude that describes the interaction between a meson and an electromagnetic current, $j_{\mu}(x)={\overline{\psi}}(x)\gamma_{\mu}\psi(x)$. In the BSE framework, the current is calculated by means of the coupling of an external photon to each of the constituents of the bound state. The corresponding hadronic matrix element is given by $J_{\mu}({p}_{\text{av}},q)=\langle\textbf{s}(p_{f})|j_{\mu}|\textbf{s}(p_{i})\rangle$, where $\textbf{s}({p}_{\text{av}})$ denotes the meson under consideration. This matrix element is related to the electromagnetic form factor $F_{\textbf{s}}$ of the meson through the expression

[TABLE]

where ${p}_{\text{av}}=(p_{f}+p_{i})/2$ is the total momentum of the meson, $p_{i}$ and $p_{f}$ denote the initial and final meson momentum, respectively, and $q=p_{f}-p_{i}$ is the photon momentum.

In the impulse approximation Eichmann and Nicmorus (2012); Eichmann et al. (2016), the conserved current $J^{\mu}({p}_{\text{av}},q)$ describes the coupling of a single photon to a quark-antiquark system, see Figure 1. The $J^{\mu}({p}_{\text{av}},q)$ is given by the following integral expression

[TABLE]

where $\mathcal{A}$ is the BS amplitude and $\bar{\mathcal{A}}$ its charge conjugated counterpart, $\Gamma^{\mu}$ represents the quark-photon vertex, and $S$ stands for the quark propagator. The integral measure is denoted by

[TABLE]

where the use of a symmetry-preserving regularization scheme is implicitly assumed. In addition, the relevant kinematic variables are defined as

[TABLE]

Here, the parameter $\eta$ determines the distribution of the total momentum ${p}_{\text{av}}$ between the quark and antiquark. Depending on the quark mass composition, a suitable choice of the value of $\eta$ leads to a tighter control of the singularities in the quark propagator. This, in turn, enhances the numerical stability of the calculation, especially for systems with significant mass disparity, such as heavy-light mesons.

The processes under consideration are elastic, having the same initial and final states, whose common bound state mass is denoted by $M$. Thus, the on-shell momenta satisfy $p_{i}^{2}=p_{f}^{2}=-M^{2}$, and ${p}_{\text{av}}$ may be parametrized as ${p}_{\text{av}}=(0,0,0,i\sqrt{M^{2}+q^{2}/4})$.

It is important to note that 2 accounts for the case where the photon interacts with only one of the valence quarks of the meson; evidently, the contribution of the interaction of the photon with the other quark must be duly added, as shown in Figure 1. Thus, the complete electromagnetic form factor for a pseudoscalar particle reads,

[TABLE]

where $e_{q}$ and $e_{\overline{q}}$ are the corresponding electric charges for the quark and antiquark.

II.2 Dynamical equations

The elements comprising 2, namely $S$, ${\cal A}$, and $\Gamma^{\mu}$ are determined from their own dynamical equations, evaluated within appropriately constructed truncation schemes. Recently, a new framework for the study of heavy-light mesons has been put forth Gao et al. (2025), where the standard one-gluon exchange interaction is complemented by flavour-dependent contributions stemming from the adjacent quark-gluon vertices (blue circles in Figure 2). In this treatment, only the classical form factor of each quark-gluon vertex is retained, and is evaluated in the so-called “symmetric” configuration. The Slavnov-Taylor identity satisfied by the quark-gluon vertex links this form factor to the quark wave-function, $A_{f}(q^{2})$, see 9, which encodes the flavour-dependence. Specifically, the resulting effective interaction, ${\mathcal{I}}_{ff^{\prime}}(q^{2})$, is given by

[TABLE]

where the effective charge ${\tilde{\alpha}}_{T}(q^{2})$ is the so-called “modified Taylor coupling”, introduced in Gao et al. (2025). The difference between ${\tilde{\alpha}}_{T}(q^{2})$ and the standard Taylor coupling, ${\alpha}_{T}(q^{2})$, is the inclusion of certain process-independent contributions extracted from the quark-gluon vertex. A useful fitting function describing ${\tilde{\alpha}}_{T}(q^{2})$ is given by Eq. (20) of Gao et al. (2025).

The effective interaction of 6 is a common ingredient of the three main dynamical equations, namely (10), (11) and (13); in fact, in all of them, it appears in the unique combination

[TABLE]

with

[TABLE]

the tree-level gluon propagator in the Landau gauge.

In particular, we have :

(${\it i}$) The quark propagator, $S_{\!f}^{ab}(p)=i\delta^{ab}S_{\!f}(p)$, where the index $f$ stands for the quark flavour, taking values $f=u,d,s,c,b$. The standard decomposition of $S_{\!f}^{-1}(p)$ is

[TABLE]

where $A_{f}(p^{2})$ and $B_{f}(p^{2})$ are the dressings of the Dirac vector and scalar tensor structures, respectively. The renormalization-group invariant quark mass function, ${\mathcal{M}}_{f}(p^{2})$, is given by ${\mathcal{M}}_{f}(p^{2})=B_{f}(p^{2})/A_{f}(p^{2})$. At tree-level, $S_{0,f}^{-1}(p)=i\not{p}+m_{f}$ , where $m_{f}$ denotes the current quark mass of the flavour $f$. In addition, the self-energy, $\Sigma_{f}(p^{2})$, is defined as $\Sigma_{f}(p^{2})=S_{\!f}^{-1}(p)-S_{0,f}^{-1}(p)$.

The evolution of the components $A_{f}(p^{2})$ and $B_{f}(p^{2})$ is determined from the quark gap equation. In terms of the kernel introduced in Gao et al. (2025), the renormalized gap equation reads

[TABLE]

where $C_{\!{\scriptscriptstyle F}}=4/3$ is the Casimir eigenvalue of the fundamental representation, and $m_{\!{\scriptscriptstyle R}}$ is the renormalized current quark mass. The diagrammatic representation of the quark self-energy $\Sigma_{f}(p^{2})$ is shown in Figure 2, panel (${\it i}$). Finally, $Z_{2}$ is the wave-function renormalization of the quark field, determined within the momentum subtraction (MOM) renormalization scheme Celmaster and Gonsalves (1979); Hasenfratz and Hasenfratz (1980); Braaten and Leveille (1981); Athenodorou et al. (2016).

(${\it ii}$) The BS amplitude $\mathcal{A}$ in 2 is given by the homogeneous BSE [see panel (${\it ii}$)] of Figure 2)

[TABLE]

with $k_{i}=k+p_{i},\,\,i=1,2$.

An alternative parametrization of the momenta, used extensively in the related literature, is given by $P=p_{1}-p_{2}$ and $p=(p_{1}+p_{2})/2$. In terms of these variables, the tensorial decomposition of the BS amplitude of a pseudoscalar meson is given by

[TABLE]

where the subamplitudes $\chi_{i}^{{\scriptscriptstyle f}{\scriptscriptstyle f}^{\prime}}:=\chi^{{\scriptscriptstyle f}{\scriptscriptstyle f}^{\prime}}_{i}(p,P)$ are functions of the Lorentz scalars $P^{2}$, $p^{2},p\cdot P$.

(${\it iii}$) The SDE of the quark-photon vertex, depicted diagrammatically in panel (${\it iii}$) of Figure 2, is given by

[TABLE]

The self-consistent treatment of 13 is presented in the next section, where the momentum-dependence of the twelve form-factors comprising $\Gamma^{f}_{\!\mu}(p_{1},p_{2})$ will be determined.

III Axial Ward-Takahashi identity and the effective kernel

Of central importance for the dynamical chiral symmetry breaking are the axial-vector current, $j^{a}_{5\mu}(x)$, and the axial current, $j^{a}_{5}(x)$, see,e.g., ​Itzykson and Zuber (1980) (Ch.11), and Miransky (1994). In terms of the fundamental quark fields, $q(x)$ and $\bar{q}(x)$, they are defined as $j^{a}_{5\mu}(x)=\bar{q}(x)\gamma_{5}\frac{\lambda^{a}}{2}\gamma_{\mu}q(x)$ and $j^{a}_{5}(x)=\bar{q}(x)\gamma_{5}\frac{\lambda^{a}}{2}q(x)$, where $\lambda^{a}$ are the Gell-Mann matrices. Then, the axial-vector vertex, $\Gamma^{a}_{5\mu}(p_{1},p_{2})$, and the axial vertex, $\Gamma^{a}_{5}(p_{1},p_{2})$, are defined as the momentum-space transforms of $\langle 0|T[j^{a}_{5\mu}(0)q(x)\bar{q}(y)]|0\rangle$ and $\langle 0|T[j^{a}_{5}(0)q(x)\bar{q}(y)]|0\rangle$, respectively, where $T$ denotes the standard time-ordering operator.

We next suppress the “isospin” indices, and introduce quark flavour indices, $f$ and $f^{\prime}$. Then, it is well-known that the axial-vector vertex $\Gamma_{5\mu}^{ff^{\prime}}(p_{1},p_{2})$ satisfies the WTI Itzykson and Zuber (1980); Miransky (1994)

[TABLE]

where $P=p_{1}-p_{2}$.

Within the framework of Gao et al. (2025), the vertices $\Gamma_{5\mu}^{ff^{\prime}}(p_{1},p_{2})$ and $\Gamma_{5}^{ff^{\prime}}(p_{1},p_{2})$ satisfy the SDEs shown in Figure 3, namely

[TABLE]

and

[TABLE]

where, for latter convenience, we use the short-hand notation $G^{ff^{\prime}}_{\!5\mu}(p_{1},p_{2})$ and $G^{ff^{\prime}}_{\!5}(p_{1},p_{2})$ to denote the non-trivial parts of the corresponding SDEs.

Evidently, upon contraction with $P^{\mu}$, the tree-level term $\gamma_{5}\gamma_{\mu}$ yields

[TABLE]

where we have used that, at tree-level, $\Gamma_{0,5}^{ff^{\prime}}(p_{1},p_{2})=\gamma_{5}$.

Turning to the quantum part of the calculation, one may directly confirm that, in the case of equal flavours, $f=f^{\prime}$, the use of the kernel $I_{ff}(q^{2})$ in 15 and 16 preserves the validity of 14. In particular, the contraction of $\Gamma_{5\mu}^{ff^{\prime}}(k_{1},k_{2})$ by $P^{\mu}$ triggers the WTI of 14 under the integral sign of $G^{ff}_{\!5\mu}(p_{1},p_{2})$, and one gets

[TABLE]

Then, the sum of 17 and 18 yields precisely the two sides of 14, with $f=f^{\prime}$.

Instead, when $f\neq f^{\prime}$, the use of the kernel $I_{ff^{\prime}}(q^{2})$ leads to a violation of 14. In particular, the repetition of the steps leading to 18 now yields

[TABLE]

where

[TABLE]

We note that, in contradistinction to $\Sigma_{f}(p_{1})$ and $\Sigma_{f}(p_{2})$, the $\widetilde{\Sigma}_{f}(p_{1})$ and $\widetilde{\Sigma}_{f}(p_{2})$ in 20 contain a non-diagonal $D_{ff^{\prime}}^{\mu\nu}(k)$, which thwarts the interpretation of these quantities as genuine self-energies. Evidently, this mismatch is the manifestation of the WTI violation associated with the use of the kernel $I_{ff^{\prime}}(q^{2})$ when $f\neq f^{\prime}$.

A way to quantify this violation is to study the quantitative difference between the components of the $S_{f}(p)$ computed using $\Sigma_{f}(p)$, and the same $S_{f}(p)$ computed using $\widetilde{\Sigma}_{f}(p)$ for all possible $I_{ff^{\prime}}(q^{2})$. In particular, if we fix $f=u$, we may first compute $\Sigma_{u}(p)$ using the correct kernel $I_{uu}$; we will denote the resulting quark propagator by $\Sigma_{uu}(p)$, where the additional subscript “$u$” indicates precisely that $I_{uu}$ was used in the computation. Then, we compute $\Sigma_{u}(p)$ using instead $\widetilde{\Sigma}_{u}(p)$, which, depending on the value of $f^{\prime}$, may contain $I_{us}$, $I_{uc}$, or $I_{ub}$; this gives rise to three versions of $\Sigma_{u}(p)$, which we denote by $\Sigma_{us}(p)$, $\Sigma_{uc}(p)$, and $\Sigma_{ub}(p)$, respectively. Employing the standard decomposition given in 9, we will denote the corresponding Dirac components of the above quark propagators by $A_{uu}(p^{2})$, $A_{uf^{\prime}}(p^{2})$, and $B_{uu}(p^{2})$, $B_{uf^{\prime}}(p^{2})$, and the resulting versions of the $u$ quark mass by ${\mathcal{M}}_{uu}(p^{2})$ and ${\mathcal{M}}_{uf^{\prime}}(p^{2})$, with $f^{\prime}=s,c,b$.

With the above ingredients in hand, we next introduce the point-wise percentage errors $(\delta A)_{uf^{\prime}}(p^{2})$ and $(\delta\mathcal{M})_{uf^{\prime}}(p^{2})$, defined as

[TABLE]

for $f^{\prime}=s,c,b$. The curves for $(\delta A)_{uf^{\prime}}(p^{2})$ and $(\delta\mathcal{M})_{uf^{\prime}}(p^{2})$, obtained from the appropriate treatment of the gap equation, are shown in Figure 4. We see that in the worst case, involving the quarks with the largest mass disparity, $f=u$ and $f^{\prime}=b$, the maximum error in the computation of the quark wave-function is slightly over $3\%$ (left panel). The maximum discrepancy in the computation of the constituent quark masses occurs again for the case $f=u$ and $f^{\prime}=b$, and is about $9\%$ (right panel).

As an illustration, we computed the mass of the heavy light meson $D$ using propagators for both the $u$ and $c$ quarks constructed with $\Sigma_{uc}(p)$, i.e., both dressed using the off-diagonal kernel $I_{uc}$. This yields a mass of $1.99~\mathrm{GeV}$, compared to the $1.93~\mathrm{GeV}$ obtained in our previous study using diagonal kernels for the propagators Gao et al. (2025).

IV Quark-Photon vertex

The quark-photon vertex is a crucial ingredient in the study of the electromagnetic interaction of hadrons Sanchis-Alepuz et al. (2013); Eichmann et al. (2016); Xu (2024); Miramontes et al. (2025b), and has been studied in isolation in a series of works Frank (1995); Maris and Tandy (2000); Miramontes and Sanchis-Alepuz (2019); Leutnant and Sternbeck (2018); Tang et al. (2019).

Just as the electron-photon vertex known from QED, the quark-photon vertex $\Gamma^{f}_{\!\mu}(p_{1},p_{2})$ satisfies the fundamental WTI

[TABLE]

It is now straightforward to establish the compatibility of the kernel $I_{ff}(q^{2})$ with the local $U(1)$ symmetry, encoded in the WTI of 22. At the level of the bare vertex SDE ($Z_{2}=1$), one contracts by $q^{\mu}$ both sides of 13, and uses the r.h.s. of 22 under the integral sign of 13, namely

[TABLE]

which, by virtue of 10, is precisely the WTI of 22. The preservation of the $U(1)$ gauge symmetry ensures charge conservation, which is crucial for the veracious calculation of hadronic form factors. Note, in addition, that the requirement that the WTI of 22 be valid after renormalization imposes the well-known relation $Z_{1}=Z_{2}$, where $Z_{1}$ denotes the renormalization constant of the quark-photon vertex; this equality has already been enforced at the level of the vertex SDE in 13, where $Z_{1}$ has been replaced by $Z_{2}$.

The vertex $\Gamma_{\!\mu}(p_{1},p_{2})$ may be naturally decomposed into two parts,

[TABLE]

where the so-called “Ball-Chiu” part, $\Gamma^{\mu}_{\!\rm{{\scriptscriptstyle B}{\scriptscriptstyle C}}}(p_{1},p_{2})$, saturates the WTI of 22, while the transverse part, $\Gamma^{\mu}_{\!\rm{\scriptscriptstyle T}}(p_{1},p_{2})$, satisfies $q_{\mu}\Gamma^{\mu}_{\!\rm{\scriptscriptstyle T}}(p_{1},p_{2})=0$. According to the construction of Ball and Chiu (1980), the form factors comprising $\Gamma^{\mu}_{\!\rm{{\scriptscriptstyle B}{\scriptscriptstyle C}}}(p_{1},p_{2})$ are expressed entirely in terms of the quark functions $A(p^{2})$ and $B(p^{2})$ appearing in 9. In particular, introducing for convenience the kinematic variables $u=(p_{1}+p_{2})/2$ and $q=p_{1}-p_{2}$, we have

[TABLE]

with the form factors given by

[TABLE]

The transverse component $\Gamma^{\mu}_{\!\rm{\scriptscriptstyle T}}(p_{1},p_{2})$ may be expanded as

[TABLE]

where the $\tau_{i}^{\mu}$ are the elements of an eight-dimensional basis, and the $h_{i}(q,u)$ denote the associated form factors. Specifically, we employ the basis of Eichmann (2014), whose elements are given by

[TABLE]

where we used the projector $t^{\mu\nu}_{ab}=\left(a\cdot b\right)\delta^{\mu\nu}-b^{\mu}a^{\nu}$, and the triple commutator is defined as $[A,B,C]=[A,B]C+[B,C]A+[C,A]B$, with $[A,B,A]=0$. It is understood that all form factors introduced above carry a flavour index “$f$”, which has been suppressed for simplicity.

Note that a key advantage of the decomposition in 24 is that it enables a separation of the transverse components with respect to the total momentum $q$. This is especially relevant for a time-like photon, where the quark-photon vertex develops vector meson bound-state poles (e.g., at $q^{2}=-M_{\rho}^{2}$). These poles reside exclusively in the transverse components, while the time-like domain may also feature non-resonant singularities originating from the analytic structure of the quark propagator Williams (2019); Miramontes and Sanchis-Alepuz (2019).

V Results

In this section we present the numerical results for the BS amplitude and quark-photon vertex. Once these elements have been obtained, the current $J_{\mu}$ can be calculated from 2, and the EFF of pseudoscalar particles may be subsequently extracted.

(${\it i}$) Numerical inputs. The current quark masses have been fixed to the following values:

[TABLE]

at the renormalization point $\mu=4.3$ GeV. Furthermore, the parameter $\eta$ that controls the routing of the total momentum ${p}_{\text{av}}$ in 4 has been fixed to the following values, which depend on the quark content of the meson:

[TABLE]

The resulting quark dressing functions $A(p)$ and $M(p)$ are the same as those computed in Fig. (5) of Ref. Gao et al. (2025). Here, it is important to remark that when employing these fixed values for $\eta_{ff^{\prime}}$ the light mesons, together with the heavy light $D$ and $D_{s}$, can be computed numerically onshell. On the other hand, for the $B,B_{s}$ and $B_{c}$ we employed a parametrization of the quark propagator in terms of complex conjugate poles, as described in Gao et al. (2025).

(${\it ii}$) BS amplitudes. To compute the BS amplitude from 11, we employ standard methods where the integral equation is reformulated as an eigenvalue problem. The physical solutions correspond to mass-shell points given by $P^{2}=-M^{2}$. In our numerical implementation, the BS amplitude is expanded using a basis of 8 Chebyshev polynomials. Figure 5 displays the resulting subamplitudes $\chi_{i}^{ff^{\prime}}$ for the pion and the $D_{s}$ meson, projected onto the leading Chebyshev moment.

(${\it iii}$) Form factors for the quark-photon vertex. For the calculation of spacelike hadron form factors, the quark-photon vertex is required in the region $q^{2}>0$. The solution of 13, unlike the homogeneous case, does not correspond to an eigenvalue problem. Instead, the equation is a linear integral equation driven by the bare vertex, which after discretization reduces to a system of linear equations where the interaction kernel defines the matrix structure. Such systems can be solved using “LU decomposition”, factorizing the kernel matrix into lower and upper triangular parts for efficient inversion Press et al. (1992). In Figure 6 and Figure 7 we present the full structure of the vertex for $u$ and $s$ quarks. As a consistency check, we have verified that our numerical results from the WTI solution from 26 are identical to those obtained from the vertex SDE. On the other hand, the transverse components are consistent with previous studies Miramontes and Sanchis-Alepuz (2019); Sanchis-Alepuz and Williams (2018).

(${\it iv}$) Electromagnetic form factors. Moving to the EFF in Figure 8 and Figure 9, we note that our results are correctly normalized, satisfying $F(0)=1$ for electrically charged particles, and $F(0)=0$ for neutral mesons, as required by the vector WTI.

The left panel of Figure 8 shows that our result for the pion electromagnetic form factor is in excellent agreement with experimental measurements. Similarly, the right panel of Figure 8 shows that our computation for the kaon electromagnetic form factor is consistent with the available experimental data. Furthermore, our results for the $D$ and $D_{s}$ mesons exhibit a similar qualitative behavior to that observed between the pion and kaon, see Figure 9, where the inclusion of a strange quark in the meson’s quark content leads to a slight difference between the form factor of $D(c\bar{d})$ and $D_{s}(c\bar{s})$, reflecting the influence of the quark mass. A similar trend is evident in the comparison between the $B(u\bar{b})$ and $B_{c}(c\bar{b})$ form factors, where the presence of a lighter quark enhances the mass asymmetry, resulting in a more pronounced variation in the form factor behavior. In the case of the $\eta_{c}$ and $\eta_{b}$ mesons, when both quark and antiquark contributions are included, the total electromagnetic form factor vanishes identically due to charge conjugation. Nevertheless, structural information can still be extracted by computing the form factor associated with a single quark contribution from 5. Our results indicate that the electromagnetic form factor for the $\eta_{b}$ exhibits a noticeably slower decrease with increasing momentum transfer compared to the $\eta_{c}$.

(${\it v}$) Charge radii. The charge radius can be determined from the derivative of the electromagnetic form factor evaluated at $q^{2}=0$:

[TABLE]

The computed charge radii for different pseudoscalar mesons are presented in Table 1, where they are compared with experimental measurements, as well as predictions from the Contact Interaction (CI) model Hernández-Pinto et al. (2023), the Algebraic model (AM) Almeida-Zamora et al. (2024), Lattice QCD (LQCD) Wang et al. (2021); Can et al. (2013), Data driven Stamen et al. (2022); Leplumey and Stoffer (2025), Light front Xu et al. (2025b), and the weighted-RL approach of Xu (2024).

Our results for the $\pi^{+}$, $K^{+}$, and $K^{0}$ mesons are in very good agreement with experimental values, showing small deviations of 0.45%, 1.42%, and 2.5%, respectively.

Comparing our results with those from other approaches, we observe some differences. The CI model underestimates the charge radii of the pion and kaon, predicting values of 0.450 fm and 0.420 fm, respectively, which are significantly smaller than both experimental measurements and our own results. This underestimation is also evident in the predictions for heavy-light mesons, where, for instance, the $D_{s}$ meson radius is given as 0.260 fm, notably lower than our 0.368 fm. Nevertheless, the overall qualitative behavior is consistent with our findings.

The Algebraic model yields a charge radius of 0.680 fm for the $D$ meson, which is substantially larger than both our prediction (0.428 fm) and lattice QCD results (0.450 fm). Additionally, this model exhibits a significant difference between the $D$ and $D_{s}$ mesons, with a predicted charge radius of 0.372 fm for the latter, whereas our results suggest that these two mesons should have comparable charge radii.

For heavy-light mesons, our results exhibit a similar overall behavior compared to the weighted-RL investigated in Xu (2024). In particular, our predictions for the $B$ and $B_{s}$ mesons, 0.631 fm and 0.330$i$ fm, respectively, are in close agreement with the corresponding weighted-RL values of 0.619 fm and 0.337$i$ fm. The most noticeable discrepancies arise in the light meson sector, where some variations in the charge radii are observed.

Finally, it is interesting to note that lattice QCD is the only approach that obtains a larger central value for the charge radius of the $D_{s}$ meson than for the $D$ meson. Nevertheless, the error for both particles is rather considerable.

VI Conclusions

In this work, we have applied the Bethe–Salpeter framework to the computation of EFF of heavy–light pseudoscalar mesons by implementing a flavour–dependent interaction kernel. The computed EFF for light mesons ($\pi$ and $K$) shows excellent agreement with experimental data, thereby validating the consistency of the approach. In the heavy–light sector, where experimental and lattice QCD results for form factors are not yet available, we provide predictions. In particular, our estimates of the charge radii are in good agreement with existing theoretical studies, offering valuable benchmarks for forthcoming lattice and experimental efforts.

Furthermore, the present framework can be naturally extended to the study of vector mesons, such as the $\rho$ and $K^{*}$, for which EFF provide complementary insight into internal meson structure beyond spin-zero systems. Beyond mesons, the formalism can also be applicable to baryons, which will be investigated in the future.

Acknowledgements

The authors thank Yin-Zhen Xu for useful communications. The work of A.S.M. and J.P. is funded by the Spanish MICINN grants PID2020-113334GB-I00 and PID2023-151418NB-I00, the Generalitat Valenciana grant CIPROM/2022/66, and CEX2023-001292-S by MCIU/AEI. J.P. is supported in part by the EMMI visiting grant of the ExtreMe Matter Institute EMMI at the GSI, Darmstadt, Germany. J.M.P. is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES Excellence Cluster) and the Collaborative Research Centre SFB 1225 - 273811115 (ISOQUANT).