The hidden charm pentaquark states and $\Sigma_c\bar{D}^{(*)}$ interaction in chiral perturbation theory

TL;DR

This paper uses heavy hadron chiral perturbation theory to calculate potentials between $ ext{Sigma}_c$ and $ar{ ext{D}}^{(*)}$, exploring their possible molecular states and relevance to hidden charm pentaquarks.

Contribution

It provides a next-to-leading order calculation of $ ext{Sigma}_car{ ext{D}}^{(*)}$ interactions including heavy quark symmetry violations and explores scenarios for pentaquark states as molecular bound states.

Findings

Reproduces $P_c(4312)$ and $P_c(4440)$ as loosely bound states.

Identifies parameter regions where $P_c(4312)$, $P_c(4440)$, and $P_c(4457)$ coexist as molecular states.

Shows coupled-channel effects can reproduce all three $P_c$ states simultaneously.

Abstract

In this work, we employ the heavy hadron chiral perturbation theory (HHChPT) to calculate the $\Sigma_c\bar{D}^{(*)}$ potentials to the next-to-leading order. The contact, the one-pion exchange and the two-pion exchange interactions are included. Besides, the mass splittings between the heavy quark spin symmetry (HQSS) multiplets are kept in calculations. Our result shows that neglecting the heavy quark symmetry (HQS) violation effect may be misleading to predict the potentials between the charmed hadrons. We perform numerical analysis with three scenarios. In the first scenario, we relate the low-energy constants (LECs) in the contact terms of $\Sigma_c\bar{D}^{(*)}$ to those of nucleon systems, and reproduce the $P_c(4312)$ and $P_c(4440)$ as loosely bound states. In the second scenario, we vary the unknown LECs and find a small parameter region in which $P_c(4312)$, $P_c(4440)$ and $P_c(4457)$ can coexist as molecular states. In the third scenario, we include the coupled-channel effect on the basis of scenario II, and notice that the three $P_c$ states can be reproduced as molecular states simultaneously in a large region of parameters. Our analytical results can be used for the chiral extrapolations in lattice QCD. With the lattice QCD results in the future as inputs, the identification of the $P_c$ states and predictions for other systems would be more reliable.

Tables

Table 1. Table 1: Matrix elements of 𝒍 1 ⋅ 𝒍 2 ⋅ subscript 𝒍 1 subscript 𝒍 2 \bm{l}_{1}\cdot\bm{l}_{2} and 𝑺 1 ⋅ 𝑺 2 ⋅ subscript 𝑺 1 subscript 𝑺 2 \bm{S}_{1}\cdot\bm{S}_{2} operators in the Σ c ( ∗ ) ​ D ¯ ( ∗ ) superscript subscript Σ 𝑐 superscript ¯ 𝐷 \Sigma_{c}^{(*)}\bar{D}^{(*)} wave functions. J 𝐽 J denotes the total spins of the two hadrons.

${\overline{\Sigma }}_c^{(\ast )}{\overline{D}}^{(\ast )};J$	${\Sigma }_c\overline{D};\frac{1}{2}$	${\Sigma }_c{\overline{D}}^{\ast };\frac{1}{2}$	${\Sigma }_c{\overline{D}}^{\ast };\frac{3}{2}$	${\overline{\Sigma }}_c^{\ast }\overline{D};\frac{3}{2}$	${\overline{\Sigma }}_c^{\ast }{\overline{D}}^{\ast };\frac{1}{2}$	${\overline{\Sigma }}_c^{\ast }{\overline{D}}^{\ast };\frac{3}{2}$	${\overline{\Sigma }}_c^{\ast }{\overline{D}}^{\ast };\frac{5}{2}$
$⟨{𝒍}_1\cdot {𝒍}_2⟩$	0	$-\frac{2}{3}$	$\frac{1}{3}$	0	$-\frac{5}{6}$	$-\frac{1}{3}$	$\frac{1}{2}$
$⟨{𝑺}_1\cdot {𝑺}_2⟩$	0	$-1$	$\frac{1}{2}$	0	$-\frac{5}{2}$	$-1$	$\frac{3}{2}$

Table 2. Table 2: The channels we considered in the coupled-channel calculations. The bold ones are our channels of interest.

Channel	1	2	3	4
$J=\frac{1}{2}$	${𝚺}_{𝒄}\overline{𝑫}$	${𝚺}_{𝒄}{\overline{𝑫}}^{\ast }$	${\Sigma }_c^{\ast }{\overline{D}}^{\ast }$	${\Sigma }_c^{\ast }\overline{D}$
$J=\frac{3}{2}$	${𝚺}_{𝒄}{\overline{𝑫}}^{\ast }$	${\Sigma }_c^{\ast }\overline{D}$	${\Sigma }_c^{\ast }{\overline{D}}^{\ast }$

Table 3. Table 3: The numerical results in the coupled-channel calculations in scenario III. The P i subscript 𝑃 𝑖 P_{i} is the proportion of the specific channel.

S-III	Exp.(MeV)	Mass$\left(\text{MeV}\right)$	RMS(fm)	$P_1(\%)$	$P_2(\%)$	$P_3(\%)$
$P_c\left(4312\right)$	$4311.9\pm {0.7}_{-0.6}^{+6.8}$	4305	1.21	99.4	0.5	0.1
$P_c\left(4440\right)$	$4440.3\pm {1.3}_{-4.7}^{+4.1}$	4446	1.22	1.0	98.0	0.9
$P_c\left(4457\right)$	$4457.3\pm {1.3}_{-4.1}^{+0.6}$	4458	1.28	96.8	2.5	0.7

Equations

$$\nu =2L+2-\frac{E_n}{2}+\sum V_i{\Delta }_i,{\Delta }_i=d_i+n_i/2-2,$$

$$\displaystyle\phi=\sqrt{2}\left(\begin{array}[]{cc}\frac{\pi^{0}}{\sqrt{2}}&\pi^{+}\\ \pi^{-}&-\frac{\pi^{0}}{\sqrt{2}}\end{array}\right),~{}U=u^{2}=\exp\left(i\frac{\phi(x)}{F_{0}}\right).$$

$${\Gamma }_{\mu }=\frac{1}{2}[u^{†},{\partial }_{\mu }u],u_{\mu }=\frac{i}{2}\{u^{†},{\partial }_{\mu }u\}.$$

$$\displaystyle\Sigma_{c}=\left(\begin{array}[]{cc}\Sigma_{c}^{++}&\frac{\Sigma_{c}^{+}}{\sqrt{2}}\\ \frac{\Sigma_{c}^{+}}{\sqrt{2}}&\Sigma_{c}^{0}\end{array}\right),~{}\Sigma_{c}^{*\mu}=\left(\begin{array}[]{cc}\Sigma_{c}^{*++}&\frac{\Sigma_{c}^{*+}}{\sqrt{2}}\\ \frac{\Sigma_{c}^{*+}}{\sqrt{2}}&\Sigma_{c}^{*0}\end{array}\right)^{\mu}.$$

$${\mathcal{L}}_{{\Sigma }_c^{\ast }\varphi }^{(0)}$$

$${\psi }^{\mu }$$

$${\overline{\text{ψ}}}^{\mu }$$

$${\mathcal{L}}_{{\Sigma }_c\varphi }^{(0)}=$$

$$g_1=-\frac{2}{3}{g}_a,g_3=-\sqrt{\frac{1}{3}}{g}_a,g_5=g_a.$$

$$\tilde{H}=({\tilde{P}}_{\mu }^{\ast }{\gamma }^{\mu }+i\tilde{P}{\gamma }_5)\frac{1-/v}{2},$$

$$\overline{\stackrel{~}{H}}=\frac{1-/v}{2}({\tilde{P}}_{\mu }^{\ast †}{\gamma }^{\mu }+i{\tilde{P}}^{†}{\gamma }_5),$$

$$\displaystyle\tilde{P}=\left(\begin{array}[]{c}\bar{D}^{0}\\ \bar{D}^{-}\end{array}\right),\quad\tilde{P}^{*\mu}=\left(\begin{array}[]{c}\bar{D}^{*0}\\ \bar{D}^{*-}\end{array}\right),$$

$$\sqrt{M}{\Phi }^{(\ast )}$$

$$\sqrt{M}{\Phi }^{(\ast )†}$$

$${\mathcal{L}}_{\bar{D}\varphi }^{(0)}=-i⟨\overline{\stackrel{~}{H}}v\cdot D\tilde{H}⟩+g_b⟨\overline{\stackrel{~}{H}}{u}_{\mu }{\gamma }^{\mu }{\gamma }_5\tilde{H}⟩-\frac{{\delta }_b}{8}⟨\overline{\stackrel{~}{H}}{\sigma }^{\mu \nu }\tilde{H}{\sigma }_{\mu \nu }⟩,$$

$${\mathcal{L}}_{\text{contact}}^{(0)}$$

$$m_{\pi }=0.139\text{ GeV},F_{\pi }=0.092\text{ GeV},$$

$${\delta }_a=0.064\text{ GeV},{\delta }_b=0.141\text{ GeV},\lambda =1\text{ GeV},$$

$$g_a=-1.47,g_b=-\mathrm{0.59.}$$

$$\mathcal{V}(\mathbf{q})=-\frac{\mathcal{M}}{\sqrt{2M_{1}2M_{2}2M_{3}2M_4}},$$

$$V(r)=\frac{1}{(2\pi {)}^3}\int d^3\mathbf{q}{e}^{i\mathbf{q}\cdot \mathbf{r}}\mathcal{V}(\mathbf{q})\mathcal{F}(\mathbf{q}).$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{X_{1.1}}=-D_1-{\tilde{D}}_1(2{\mathbf{I}}_1\cdot {\mathbf{I}}_2),$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{F_{1.1}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{T_{1.1}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{T_{1.2}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{T_{1.3}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{B_{1.1}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{B_{1.2}}$$

$${\mathcal{V}}_{{\Sigma }_c\bar{D}}^{R_{i,j}}$$

$$q^i{q}^j\rightarrowtail {\delta }^{ij}\frac{1}{d-1}{\mathbf{q}}^2.$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{X_{2.1}}$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{H_{2.1}}$$

$${{\mathbf{ϵ}}^{\ast }}_4\cdot {\mathbf{ϵ}}_2\rightarrowtail 1,q^i{q}^j\rightarrowtail {\delta }^{ij}\frac{1}{d-1}{\mathbf{q}}^2.$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{F_{2.1}}$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{T_{2.1}}$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{T_{2.2}}$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{T_{2.3}}$$

$${\mathcal{V}}_{{\Sigma }_c{\bar{D}}^{\ast }}^{T_{2.4}}$$

Full text

I Introduction

The multiquark state has been a very important topic in hadron physics for a long time (for recent reviews, see Refs. Chen et al. (2016); Esposito et al. (2016); Guo et al. (2018); Olsen et al. (2018); Liu et al. (2019a)). Since 2003, many “XYZ” states have been observed as the candidates of tetraquark states Choi et al. (2003); Aubert et al. (2005); Ablikim et al. (2013); Liu et al. (2013). In 2015, the LHCb Collaboration discovered two pentaquark candidates $P_{c}(4380)$ and $P_{c}(4450)$ in the $J/\psi p$ invariant mass spectrum of $\Lambda_{b}\rightarrow J/\psi Kp$Aaij et al. (2015). Very recently, the LHCb Collaboration reported the new results about pentaquarks Aaij et al. (2019). The previously reported $P_{c}(4450)$ was resolved into two narrow states $P_{c}(4440)$ and $P_{c}(4457)$. Besides, a new state $P_{c}(4312)$ was observed with 7.3$\sigma$ significance, and these three states are measured to be narrow. The masses of these resonances lie below the thresholds of $\Sigma_{c}\bar{D}$ and $\Sigma_{c}\bar{D}^{*}$, respectively. Thus, they are good candidates of the molecular states. Before the discovery of $P_{c}$ states, several groups predicted the existence of the hidden charm molecular states with five quarks Wang et al. (2011); Yang et al. (2012); Wu et al. (2012). Although there are several alternative explanations, like tightly-bound pentaquarks Ali and Parkhomenko (2019); Maiani et al. (2015); Wang and Huang (2016); Weng et al. (2019) and kinetic effect Guo et al. (2015), the molecular state explanation is more favorable. The recent QCD sum rules Chen et al. (2019a) and the one-boson-exchange (OBE) model Chen et al. (2019b); He (2019), the local hidden gauge approach Xiao et al. (2019a) and the quark delocalization color screening model Huang et al. (2019) calculations also support the molecule explanation for the $P_{c}(4312)$, $P_{c}(4440$) and $P_{c}(4457)$ states. The productions and decays of the newly observed pentaquark states were also investigated in Refs. Guo et al. (2019); Xiao et al. (2019b); Cao and Dai (2019).

The OBE model is widely used to study the nuclear force Yukawa (1935); Machleidt et al. (1987). For instance, the deutron was well established as a hadronic molecular state in the framework of OBE model Machleidt et al. (1987). Additionally, in the heavy flavor sector, the molecules with four quarks, five quarks and six quarks were investigated with the OBE model as well Liu et al. (2009); Li and Zhu (2012); Chen et al. (2015); Meng et al. (2017); Yang et al. (2019).

Chiral perturbation theory (ChPT) is used to build the modern the nuclear force Epelbaum et al. (2009); Machleidt and Entem (2011). The idea was proposed by Weinberg Weinberg (1990, 1991). The chiral expansion is performed to obtain the interaction kernel, which is iterated to all orders by solving the Lippmann-Schwinger equation or Schrödinger equation. Compared with the OBE model, the ChPT has the consistent power counting. The potential is calculated order by order, and thus the error is estimable and controllable. The same idea on nuclear force was also applied to the heavy hadron systems. The interactions between two heavy mesons were studied with the ChPT by taking the heavy quark symmetry (HQS) into account Valderrama (2012); Wang and Wang (2013); Baru et al. (2015); Liu et al. (2014); Xu et al. (2019); Wang et al. (2019). The $X(3872)$ and two $Z_{b}$ states were obtained as bound states under this framework.

The HQS was also used to predict the partner states of $P_{c}(4312)$, $P_{c}(4440$), $P_{c}(4457)$ Liu et al. (2019b); Shimizu et al. (2019). The HQS is a good approximation when the heavy quark masses approach the infinity. As we know, the heavy quark spin symmetry (HQSS) violation effect will lead to mass splittings in the charmed sector between the HQSS multiplets. The impacts of the HQS breaking effect on the heavy molecular states are rarely estimated. The molecular states are very shallow bound states and they are subtle to the behaviors of the potentials. Thus, the effect of HQSS violation on the interactions between two heavy hadrons needs to be carefully considered, especially in the charmed sector.

In this work, we derive the potentials between $\Sigma_{c}$ and $\bar{D}^{(*)}$ in the framework of heavy hadron chiral perturbation theory (HHCPT). We try to reproduce the the newly observed $P_{c}(4312)$, $P_{c}(4440$) and $P_{c}(4457)$ as the molecular states. In Sec. II, we discuss the Weinberg’s formalism and construct the Lagrangians. In Sec. III, we perform the Feynman diagrams calculation in the framework of the HHChPT and obtain the analytical results of the effective potentials between $\Sigma_{c}$ and $\bar{D}^{(*)}$ to the next-to-leading order. We give some discussions about the $\Sigma_{c}\bar{D}^{(*)}$ potentials in the heavy quark limit in Sec. IV. In Sec. V, we use three scenarios to give the numerical results. A brief summary is given in Sec. VI. In Appendix A, we list the matrix elements of some operators. In Appendix B, we present the loop integral functions we used.

II Weinberg’s Formalism and Effective Lagrangians

In the framework of the HHChPT, the amplitudes are expanded in powers of $\epsilon=q/\Lambda_{\chi}$, where $q$ is either the momenta of Goldstone bosons or the residual momenta of the matter fields, and $\Lambda_{\chi}$ is the chiral symmetry breaking scale. For the singly heavy hadrons, the mass splitting $\delta$ between the heavy quark multiplets is not vanishing in the chiral limit. Thus, we adopt the small scale expansion in this work Hemmert (1997), where the mass spitting $\delta$ is regarded as another small scale. The amplitudes are also expanded in powers of $\delta/m_{c}$. $m_{c}$ is the heavy quark mass, which is treated as another large scale.

The expansion is organized according to the power counting given by Weinberg Weinberg (1990, 1991). The order of a diagram $\nu$ reads,

[TABLE]

where $L$ and $E_{n}$ are the numbers of loops and external lines of matter fields, respectively. For the $\Sigma_{c}\bar{D}^{(*)}$ potential, $E_{n}=2$. $V_{i}$ is the number of the vertices with order $\Delta_{i}$. $d_{i}$ and $n_{i}$ are the numbers of the derivatives and matter field lines, respectively.

In the Weinberg’s formalism, only two particle irreducible (2PIR) graphs are considered. The amplitudes of the box diagrams would be enhanced by the pinch singularities, which would destroy the power counting in Eq. (1). The pinch singularity originates from the two intermediate on-shell matter fields. Thus, one can recover the power counting by excluding the two particle reducible (2PR) contributions. The amplitudes we get in this way serve as the kernel of Lippmann-Schwinger equation or Schrödinger equation. The tree level one-pion exchange diagrams would be iterated to generate the 2PR contributions automatically by solving the Lippmann-Schwinger equation or Schrödinger equation.

In order to remove the 2PR contributions in the box diagrams, Ordonez et al. adopted the time-ordered perturbation theory Ordonez et al. (1996), while Kaiser et al. removed the contributions from the poles of the intermediate matter fields when performing the loop integrals Kaiser et al. (1997). Here, we use the principal integral method, which is equivalent to the scheme used in Ref. Kaiser et al. (1997). The details can be found in Appendix B.

In order to construct the chiral Lagrangians, we introduce the pion fields,

[TABLE]

The chiral connection $\Gamma_{\mu}$ and the axial vector current $u_{\mu}$ are defined as

[TABLE]

The multiplets of $\Sigma_{c}^{(*)}$ are denoted as

[TABLE]

Their chiral covariant derivative is $D_{\mu}\Sigma^{(*)}_{c}=\partial_{\mu}\Sigma^{(*)}_{c}+\Gamma_{\mu}\Sigma^{(*)}_{c}+\Sigma^{(*)}_{c}\Gamma_{\mu}^{T}$. The leading order chiral Lagrangians for the $\Sigma_{c}^{(*)}$ read,

[TABLE]

where $\text{Tr}[...]$ represents the trace in flavor space. Since the $\Sigma_{c}$ and $\Sigma_{c}^{*}$ are the degenerate states in the heavy quark limit, we can define the superfield as

[TABLE]

where $\mathcal{B}^{*}$ stands for the $\Sigma_{c}^{(*)}$ fields after heavy baryon reduction. The leading order Lagrangian in Eq. (11) can be rewritten as

[TABLE]

The third term in Eq. (13) accounts for the HQS violation effect for the charmed baryons. $\delta_{a}=M_{\Sigma_{c}^{*}}-M_{\Sigma_{c}}$ denotes the mass splitting. Comparing Eq. (13) with Eq. (11), one can easily get

[TABLE]

We introduce the superfield $\tilde{H}$ to denote the $\bar{D}$ and $\bar{D}^{*}$ fields,

[TABLE]

where we use the “$\tilde{X}$” to label the antiparticles. The $\tilde{P}$ and $\tilde{P}^{*}_{\mu}$ are the reduced heavy meson $\bar{D}$ and $\bar{D}^{*}$ fields, respectively, which can be related to the relativistic $D^{(*)}$ fields $\Phi^{(*)}$ by the following relations

[TABLE]

where $M$ is the $D^{(*)}$ meson mass. Their chiral covariant derivative is $D_{\mu}\tilde{P}^{(*)}=\partial_{\mu}\tilde{P}^{(*)}+\Gamma_{\mu}\tilde{P}^{(*)}$. The leading order Lagrangians for $\bar{D}$ and $\bar{D}^{*}$ are

[TABLE]

where $\langle...\rangle$ denotes the trace in spinor space. The third term in Eq. (22) represents the HQS violating effect for the charmed mesons. $\delta_{b}=M_{\bar{D}^{*}}-M_{\bar{D}}$ is the mass splitting between $\bar{D}^{*}$ and $\bar{D}$.

Apart from Eqs. (13) and (22), the leading order Lagrangians also contain the contact terms,

[TABLE]

In this work, the values of the parameters Tanabashi et al. (2018) are taken as

[TABLE]

where $\lambda$ is the large scale we used to perform the small scale expansion. The $g_{a}$ and $g_{b}$ are determined by the partial decay widths of $D^{*}$ and $\Sigma_{c}^{*}$ Ahmed et al. (2001); Meguro et al. (2011); Tanabashi et al. (2018), respectively.

[TABLE]

The signs of $g_{a}$ and $g_{b}$ are determined by the quark model. The LECs in the contact terms will be estimated in Sec. V.

III effective potentials

In order to obtain the effective potential, we calculate the scattering amplitude $\cal{M}$ first. If we use the standard normalization and Feynman rules in the relativistic quantum field theory, the effective potential $\mathcal{V}(\bm{q})$ in the momentum space reads,

[TABLE]

where $M_{1,\dots,4}$ are the masses of the scattering particles. The potential in coordinate space can be obtained by making the Fourier transformation to the $\mathcal{V}(\bm{q})$. In order to regularize the divergence in the Fourier transformation, we introduce the Gauss regulator, ${\cal F}(\bm{q})=\exp(-\bm{q}^{2n}/\Lambda^{2n})$  Ordonez et al. (1996); Machleidt and Entem (2011); Ren et al. (2018). The potential in the coordinate space reads,

[TABLE]

In this work, we set $n=2$ and vary the cutoff $\Lambda$ from $0.5$ GeV to $0.8$ GeV.

III.1 $\Sigma_{c}\bar{D}$ system

The leading order ($\mathcal{O}(\epsilon^{0})$) potential of the $\Sigma_{c}\bar{D}$ system comes from the contact diagram in Fig. 1. There is no one-pion exchange diagram due to the vanishing $\bar{D}\bar{D}\pi$ vertex. This vertex is forbidden by the parity and angular momentum conservation. Thus, the leading order potential of $\Sigma_{c}\bar{D}$ reads,

[TABLE]

where $\bm{I}_{1}\cdot\bm{I}_{2}$ is the isospin-isospin operator. The superscript is the label of the Feynman diagram while the subscript denotes the physical system.

There are large amounts of contact terms contributing to the effective potential at the next-to-leading order ($\mathcal{O}(\epsilon^{2})$) Wang et al. (2019). The LECs of these vertices can be divided into the infinite part and finite part. The infinite part can be used to absorb the divergence in the $\mathcal{O}(\epsilon^{2})$ loop diagrams. The renormalized vertices (the finite part) will contribute to the effective potential. In our calculations, we neglect the finite part of these vertices since we have no experimental data as input to determine the LECs at $\mathcal{O}(\epsilon^{2})$.

There are loop diagrams contributing to the next-to-leading order potential, which can be divided into two types. The first type is the vertex correction and wave function renormalization diagrams in Fig 2. Their contributions can be included when we use the physical values of the parameters in Lagrangians. Another type is the two-pion exchange diagrams in Fig 3, including one football diagram, three triangle diagrams, two box diagrams and two crossed box diagrams. The two-pion vertices in these diagrams stem from the chiral connection terms in Eqs. (13) and (22). The one-pion vertices arise from the axial coupling terms. In the calculation, we keep the mass splittings from the HQS violation effect in the propagators. The analytical results of these diagrams read

[TABLE]

where $J^{T}_{ij}$ are the loop integrals defined in the Appendix B. They are the functions of $m_{\pi}$, $\bm{q^{2}}$ and mass splittings $\delta_{a,b}$. We omit the $m_{\pi}$, $\bm{q^{2}}$ for conciseness and give the specific mass splittings at the end of every expression. The $\bm{q}$ is the transferred three-momentum in the diagrams, and $d$ is the dimension in the dimensional regularization. Since we focus on the $S$-wave interactions, we take the replacement,

[TABLE]

III.2 $\Sigma_{c}\bar{D}^{*}$ system

For the $\Sigma_{c}\bar{D}^{*}$ system, the leading order potential is generated from both the contact terms and the one-pion exchange diagram. The analytical results read,

[TABLE]

where $\bm{\sigma}$ is the Pauli matrix. The $\bm{T}\equiv i\bm{\epsilon^{*}_{4}}\times\bm{\epsilon_{2}}$ is proportional to the spin operator of $\bar{D}^{*}$. Thus, the $\bm{\sigma}\cdot\bm{T}$ terms are the spin-spin interaction. For the $S$-wave potential, we can make the following replacements in the effective potentials,

[TABLE]

Then, only the central terms and spin-spin terms survive in the potential.

For the next-to-leading order potential, we neglect the finite part from the $\mathcal{O}(\epsilon^{2})$ contact terms again. The renormalizations of the vertices, wave functions and masses are included by using the physical values of the coupling constants, decay constants and masses just like we did in the previous subsection.

At $\mathcal{O}(\epsilon^{2})$, there are twelve two-pion exchange diagrams which contribute to the effective potential in Fig. 5. Their analytical results read,

[TABLE]

where the notations are the same as those for the $\Sigma_{c}\bar{D}$ system.

IV The heavy quark symmetry

In Secs. II and III, the LECs can be related to one another by adopting the HQS. The HQS violation effect is introduced through the mass splittings. If we ignore these mass splittings in the loop diagrams, the HQS manifests itself.

The manifestation of the HQS could be clearer at the quark level. In the heavy quark limit, the potential between the $\Sigma_{c}$ and $\bar{D}^{(*)}$ arises from the interactions of their light degrees of freedom. The heavy degrees of freedom are spectators. Their interactions are suppressed by the heavy quark mass. The $S$-wave interactions between the light diquark in the $\Sigma_{c}$ and the light quark in the $\bar{D}^{(*)}$ can be expressed as,

[TABLE]

where $\bm{l}_{1}$ and $\bm{l}_{2}$ are the spin operators of their light degrees of freedom. Since we concentrate on the $S$-wave interactions, only the spin-spin interaction and the central potential exist. We can parameterize the potential at the hadron level as

[TABLE]

With the quark level interaction in the HQS, we can relate the hadron level potentials with each other as follows,

[TABLE]

We give the matrix elements of $\bm{l}_{1}\cdot\bm{l}_{2}$ and $\bm{S}_{1}\cdot\bm{S}_{2}$ in Table 1 and the calculation details in Appendix A. Our analytical results indeed satisfy the above expressions when $d\rightarrow 4$ and $\delta_{a,b}\rightarrow 0$.

With Eq. (50), we can get the potentials of the $\Sigma_{c}^{*}\bar{D}$, $\Sigma_{c}^{*}\bar{D}^{*}$ and even the inelastic channels in the heavy quark limit without calculating the loop diagrams. We first extract the $V_{a}$ and $\tilde{V}_{a}$ from the $\Sigma_{c}\bar{D}^{*}$ potentials as in Eq. (52). Then, we calculate the matrix elements of $\bm{l}_{1}\cdot\bm{l}_{2}$ for the corresponding channels. The potential can then be derived from the quark-level interaction in Eq. (50). In our framework, the leading order potentials satisfy the HQS, i.e., the leading order potentials we derived from above procedures are equal to those from Feynman diagrams.

In the heavy quark limit, the analytical results of the box diagrams become more concise. All the mass splittings between heavy quark multiplets vanish with the HQS, and all the box diagrams with the pinch singularities become the 2PR diagrams. With the expression of $J^{B}_{x}$ in Appendix B, there exists the relation $J^{B}_{x}=-J^{R}_{x}$ when $\delta_{a}=\delta_{b}=0$. The total potentials of the box diagrams and crossed box diagrams are

[TABLE]

The results are similar to those for the nuclear force  Machleidt and Entem (2011).

However, the HQS is still an approximation when the heavy quark mass is not infinite. Whether the HQS is good enough in calculating the heavy hadron potential to obtain the bound states needs to be considered carefully. In our calculation, we keep the mass differences stemming from the HQS breaking effect. For the triangle diagrams, there is no spin-spin interaction. Thus in the heavy quark limit, the $\Sigma_{c}\bar{D}$ and $\Sigma_{c}\bar{D}^{*}$ should have the same potential. We give the potentials of the triangle diagrams for the $I={1\over 2}$ systems in Fig. 6, from which we see that the $\Sigma_{c}\bar{D}^{*}$ potential in the heavy quark limit is very close to its real potential. However, the $\Sigma_{c}\bar{D}$ potential in the heavy quark limit deviates a lot from its real potential. Similar results are obtained for the crossed box diagrams in Fig. 6 as well. For the $\Sigma_{c}\bar{D}^{*}$ system, the potential in the HQS is a good approximation of its real potential. However, the real $\Sigma_{c}\bar{D}$ attractive potential becomes repulsive when the HQS is adopted. In the triangle or crossed box diagrams of the $\Sigma_{c}\bar{D}$ system, ignoring the HQS violation effect will make the potential in the coordinate space change by about 0.02-0.03 GeV. Our numerical results in Sec. V indicate that the minimum of the potential function that generates the loosely bound state is from -0.06 GeV to -0.15 GeV. The correction from the HQS violation is not negligible. Therefore, it may be misleading to adopt the HQS in calculating the charmed hadron potentials, at least for the $\Sigma_{c}\bar{D}$ system.

We notice that the HQS violation effect is more significant for the $\Sigma_{c}\bar{D}$ system than that in the $\Sigma_{c}\bar{D}^{*}$ case. We take the triangle diagrams as an example to illustrate the reason. In Fig. 3, the graphs $(T_{1.1})$ and $(T_{1.3})$ are two diagrams with the HQS violation effect for the $\Sigma_{c}\bar{D}$ system. Their intermediate states are $\Sigma_{c}^{*}$ and $\bar{D}^{*}$, respectively, which are both heavier than the corresponding external particles. The HQS violation effect will deviate the potential in the same direction. Thus, the HQS violation effect from different diagrams are constructive. For the $\Sigma_{c}\bar{D}^{*}$ system, the intermediate state can be either heavier or lighter than its corresponding external field. The HQS breaking (HQSB) effect would cancel with each other. Thus, we can infer that the HQS violation effect from mass splittings is also significant for the $\Sigma_{c}^{*}\bar{D}^{*}$ system.

The HQS violation effect can also be investigated at the quark level. To this end, we present this effect as follows,

[TABLE]

where $V_{c}$, $V_{d}$ and $V_{e}$ are the functions used to parameterize the potential. The first and second terms are the interaction between the light and heavy degrees of freedom, which are suppressed by the $1/m_{c}$. The third term is the interaction between the heavy degrees of freedom, which is the higher order contribution in the heavy quark expansion. For the $\Sigma_{c}\bar{D}$ system, we calculate the matrix elements of three spin-spin operators in Eq. (54),

[TABLE]

The HQS violation effect vanishes, which seems to be contradictory with the conclusion from directly calculating Feynman diagrams. The HQS violation effect in calculating the loop diagrams arises from the mass splittings in the propagators. It is hard to include this effect in the quark model. Quark model can only give the analytical terms which are the polynomials of $m_{\pi}^{2}$ or $\delta$. However, the loop diagrams in ChPT can generate the nonanalytical structures such as the logarithmic terms. The quark level HQS violation effect in Eq. (54) is more likely to appear in the LECs at the hadron level.

The heavy quark limit for the box diagrams is more tricky. When the intermediate states of the box diagrams are the HQS partner states of the external fields, we calculate these diagram directly. They have no pinch singularities due to the existence of the mass splittings. If we decrease the mass splittings to zero, these amplitudes will blow up and these diagrams become two particle reducible. One way to eliminate these singularities is to remove the 2PR contributions. Meanwhile, we shall solve the coupled-channel Schrödinger equation as illustrated in Fig. 7. The one-pion exchange inelastic scattering diagrams at the tree level will be iterated to generate the 2PR contributions automatically. Thus, it is illegitimate to take the mass splittings in our analytical results of the box diagrams to zero directly.

V Numerical results

In our analytical results, there are four unknown LECs from the contact terms, $D_{1}$, $D_{2}$, $\tilde{D}_{1}$ and $\tilde{D}_{2}$. These LECs should be determined by fitting the experimental data. However, there does not exist any $\Sigma_{c}$ and $\bar{D}^{(*)}$ scattering data now. Therefore, we choose three scenarios to present the numerical results in this section. In scenario I, we determine these LECs from the nucleon scattering with the help of quark model. In scenario II, we will assume the three structures observed in LHCb corresponding to three molecular states and check whether they can coexist in our framework. In scenario III, we will include the coupled-channel effect on the basis of scenario II.

V.1 Scenario I

There are two motivations to introduce the contact terms. First, some heavy mesons exchanged between $\Sigma_{c}$ and $\bar{D}^{(*)}$ like $\rho$ and $\omega$, are integrated out, and their contributions are included in the contact terms. Second, the contact terms will absorb the divergence in the loops and remove the scale dependence. The contact terms will cancel the infinity in the loop diagrams. Thus, the values of renormalized contact terms depend on chiral truncation order. Meanwhile, the contact terms will make the potential scale independent. The values of the contact terms will depend on the regularization schemes, types of regulator and values of cutoff. The specific values of the contact terms for the nucleon system in literatures also vary due to above reasons Ordonez et al. (1996); Epelbaum et al. (2005); Machleidt and Entem (2011).

For the nucleon system, the leading order contact Lagrangian and the potential are written as

[TABLE]

There are two independent LECs, $C_{S}$ and $C_{T}$. The isospin-isospin interaction is absorbed through Fierz rearrangement. Since we use the dimensional regularization in calculations, we choose the values determined in Ref. Machleidt and Entem (2011), in which the same regularization scheme was used. We take the values

[TABLE]

In order to relate the contact terms in the nucleon system to those of the $\Sigma_{c}\bar{D}^{(*)}$ systems, we make use of the contact interaction at the quark level,

[TABLE]

where $q=(u,d)^{T}$. $c_{s}$ and $c_{t}$ are the coupling constants at the quark level. The interaction may arise from the heavy meson exchange at the quark level. With the quark model, we can get the values of the LECs for the $\Sigma_{c}\bar{D}^{(*)}$ systems,

[TABLE]

With these LECs, we can solve the Schrödinger equation to obtain some bound states. We vary the cutoff $\Lambda$ from $0.4$ GeV to $0.8$ GeV. The binding solutions are given in Fig. 8. For the $I={1\over 2}$ system, we got binding solutions for the $[\Sigma_{c}\bar{D}]_{J=1/2}$ and $[\Sigma_{c}\bar{D}^{*}]_{J=1/2}$ systems. We reproduce the masses of $P_{c}(4312)$ and $P_{c}(4440)$ when the cutoff is $0.5$ GeV. We present the potentials in Fig. 9. The binding energy and the root mean square radius (RMS) for the $[\Sigma_{c}\bar{D}]_{J=1/2}$ are $-9.21$ MeV and $1.36$ fm, respectively. For the $[\Sigma_{c}\bar{D}^{*}]_{J=1/2}$, the binding energy and the RMS are $-18.93$ GeV and $1.16$ fm, respectively.

However, in this scenario, we can not reproduce the $P_{c}(4457)$. In fact, the potential of $[\Sigma_{c}\bar{D}^{*}]_{J=3/2}^{I=1/2}$ is repulsive as shown in Fig. 9. The scheme we used to determine the LECs is rather rough. Therefore, we can not rule out the possibility of the $P_{c}(4457)$ as a molecular state because of the uncertainty of the LECs. For all three $I={3\over 2}$ systems, there exist loosely bound states when we vary the cutoff $\Lambda$ from $0.4$ GeV to $0.8$ GeV. We give the results in Fig. 8.

V.2 Scenario II

One can use other phenomenological methods to evaluate the LECs, such as heavy meson exchange model. But they also bring large uncertainties. Meanwhile, we drop out the finite contributions from $\mathcal{O}(\epsilon^{2})$ contact terms, which may also influence our final results. In scenario II, we will adopt the general form of contact terms and vary the LECs to search for the bound solutions.

In this scenario, we will focus on the $I=1/2$ systems since the three $P_{c}$ states were all observed in the $J/\psi p$ invariant mass spectrum. Thus, there are only two independent contact terms. We parameterize the contact interaction of $[\Sigma_{c}\bar{D}^{(*)}]^{I=1/2}$ as,

[TABLE]

The isospin-isospin interaction and the $\mathcal{O}(\epsilon^{2})$ contact terms are absorbed into $\mathbb{D}_{1}$ and $\mathbb{D}_{2}$.

We vary $\mathbb{D}_{1}$ and $\mathbb{D}_{2}$ in the range from $-100\text{ GeV}^{-2}$ to $100\text{ GeV}^{-2}$, respectively. We show the parameter regions in which there exist loosely bound states for $[\Sigma_{c}\bar{D}]^{I=1/2}_{J=1/2}$, $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=1/2}$ and $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=3/2}$ in Fig. 10, where we choose $\Lambda=0.5$ GeV. Since the molecules are loosely bound states, we adopt the binding energy $E=-30$ MeV as the lower limit. In Fig. 10, there is a small region, in which three bound states can coexist. In this region, the binding energy ranges for $[\Sigma_{c}\bar{D}]^{I=1/2}_{J=1/2}$, $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=1/2}$ and $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=3/2}$ are $[-30,-25]$, $[-11,0]$ and $[-8,-4]$ MeV, respectively.

We choose one set of parameters in the overlap region of three bands in Fig. 10, $\mathbb{D}_{1}=42\text{ GeV}^{-2}$ and $\mathbb{D}_{2}=-25\text{ GeV}^{-2}$. The potentials are displayed in Fig. 11, where the potentials for all three channels are attractive. For the $[\Sigma_{c}\bar{D}]^{I=1/2}_{J=1/2}$ system, the potential from the contact terms and two-pion exchange are both attractive. The binding energy is also deeper than that of $P_{c}(4312)$ as a $\Sigma_{c}\bar{D}$ bound state. For the $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=1/2}$ system, the interaction of the two-pion exchange are very weak. The attractive one-pion exchange and contact interactions generate a loosely bound state. For the $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=3/2}$ system, both one-pion exchange and two-pion exchange are repulsive. The loosely bound state arises from the very attractive contact interaction. The bound state of $[\Sigma_{c}\bar{D}^{*}]^{I=1/2}_{J=3/2}$ we got is dominated by the short-distance contact interaction, which may arise from the vector meson $\rho$ and $\omega$ exchange in the OBE mode.

We draw the binding energies of threes $P_{c}$ states as three solid lines in Fig. 10 if we assume they are molecular states. There are three cross points in which two of three states can coexist. The three cross points are not very close. In other words, if we restrict the binding energies to the experimental values, it is hard to reproduce the three states simultaneously.

V.3 Scenario III

In the above two scenarios, we only consider the potentials of the elastic channels as shown in the upper panel of Fig. 7. We include the HQS partner states of the external lines as intermediate states. Only part of the coupled-channel effects is taken into account. For example, the contributions from the ladder diagrams generated by the inelastic tree diagrams are dropped. In scenario III, we improve our results by including the inelastic channels and solving the coupled-channel Schrödinger equation.

For the $J={1\over 2}$ and $J={3\over 2}$ systems, four and three channels can couple to one another, respectively, which are shown in Table 2. The $[\Sigma_{c}\bar{D}]_{J=1/2}$, $[\Sigma_{c}\bar{D}^{*}]_{J=1/2}$ and $[\Sigma_{c}\bar{D}^{*}]_{J=3/2}$ are the channels we are interested in. For the other channels, we only include their leading order potentials. We can get these potentials either from the tree diagram calculations or the HQS analysis as illustrated in Sec. IV. Both approaches lead to the same results. The matrix elements of the $\bm{l}_{1}\cdot\bm{l}_{2}$ for these channels read

[TABLE]

Since the off-diagonal terms in the Hamiltonian only arise from the $\bm{l}_{1}\cdot\bm{l}_{2}$ interaction, the forth channel of the $J=1/2$ system does not couple with the other three channels in the leading order potentials.

For our interested channels, we compute their potentials up to the next-to-leading order. In the coupled-channel calculations, we deal with the box diagrams in a different way. As illustrated in Fig. 7, we have included the inelastic tree diagrams at the leading order. These diagrams will be iterated to generate the ladder diagrams. Thus, to avoid the double counting of this contribution, we ignore the mass splittings in the box diagrams and remove the 2PR contributions.

We vary the two unknown LECs $\mathbb{D}_{1}$ and $\mathbb{D}_{2}$ as we did in scenario II. The parameter regions in which the three loosely bound states (with binding energy $-30$ MeV$\sim 0$ MeV) can exist are shown in the second graph of Fig. 10. Because of the coupled-channel effect, the parameter regions are no longer bands with fixed-width. We notice that the overlap region in which three bound states can coexist is much larger than that in scenario II. Meanwhile, the three cross points of the three lines corresponding to $P_{c}(4312)$, $P_{c}(4440)$ and $P_{c}(4457)$ are much closer than those in scenario II. Therefore, the coupled-channel effect is very important to reproduce the three states simultaneously.

We choose one set of LECs in the region that three molecular states can coexist to give numerical results, where we take $\mathbb{D}_{1}=55\text{ GeV}^{-2}$ and $\mathbb{D}_{2}=-10\text{ GeV}^{-2}$. The masses, the RMSs and the proportion of each channel are given in Table 3. With this set of LECs, we reproduce the three $P_{c}$ states in the experiment as the molecular states simultaneously. Among them, the coupled-channel effect is more significant for the formation of $P_{c}(4457)$. We plot the potentials of their dominant channels in Fig. 12. We notice that the attractive interaction mainly stems from the contact interaction, which is similar to that in scenario II. As we discussed before, the dominant contact interaction may arise from the vector meson $\rho$ and $\omega$ exchange in the OBE model.

VI Discussion and Summary

In summary, we calculate the effective potentials of the $\Sigma_{c}\bar{D}^{(*)}$ systems in the heavy hadron perturbation theory. We adopt the small scale expansion to keep the mass splittings between the HQS multiplets. We include the contact interaction and one-pion exchange interaction at the leading order. At the next-to-leading order, we take the two-pion exchange interaction into consideration. The renormalizations of vertices and wave functions are included by taking the physical values of the parameters in chiral Lagrangians.

We employ the quark model with heavy quark spin symmetry to get some relations between different systems. Our analytical results are consistent with these relations in the heavy quark limit. Using these relations, we obtain the potentials of the partner channels in the HQS without calculating extra diagrams. We also show that the HQS violation effect is not negligible in calculating the potentials between the charmed hadrons. Since the molecular states are very shallow bound states, their existence is very sensitive to the potentials. The molecular states calculated in the charmed sector in the heavy quark limit might be misleading. For the box diagrams, taking the heavy quark limit will make the original two particle irreducible diagrams reducible. Extra operation to remove the 2PR contributions is needed. One should be cautious about the uncertainty and trap of using the HQS to obtain the potentials between the charmed hadrons.

Due to lack of experimental data, we can not determine the four LECs in the contact terms precisely. We use three scenarios to estimate the contact interactions. In the first scenario, we assume a contact interaction at the quark level phenomenologically and then relate the contact terms in nuclear force to those in the $\Sigma_{c}\bar{D}^{(*)}$ systems. With the LECs, we reproduce the $P_{c}(4312)$ and $P_{c}(4440)$ and predict three loosely bound states in the $I={3\over 2}$ channels. We are unable to reproduce $P_{c}(4457)$ due to the large uncertainty of LECs in the first scenario. In the second scenario, we focus on the $I={1\over 2}$ channels. There are only two unknown independent LECs. We vary the two LECs and search for the region that three $P_{c}$ states can coexist as the loosely bound states. We do find a parameter region in which we can reproduce three $P_{c}$ states simultaneously. The region is very small. The solution corresponding to the $P_{c}(4457)$ seems slightly less natural, since the attractions all arise from the short-range contact interactions. In the third scenario, we consider the coupled-channel effect in the leading order on the basis of scenario II. To avoid the double counting, we neglect the mass splittings in the box diagrams and remove their 2PR contributions. Through the coupled-channel calculations in scenario III, we obtain a large parameter region in which three $P_{c}$ states can coexist as molecular states. The attraction mainly comes from the contact interactions.

We have reproduced the three $P_{c}$ states in our calculations. The uncertainties come from either the framework or the LECs. In this work, we only take the $S$-wave into consideration. However, the $S$-$D$ wave mixing plays an important role in reproducing the binding energy of the deutron. We do not consider the HQS violation effect in the LECs. In order to reduce the number of LECs, we relate them to each other through HQS. The approximation will introduce errors, especially for the contact terms. The physical information at the high energy scale is packaged into the contact terms, in which the HQSS tends to be broken. Finally, the coupled-channel effect can be considered more carefully. In this framework, we calculate the potentials of the interested channels to the next-to-leading order. However, for the other inelastic channels, we only calculate their leading order potentials. The numerical results could be improved if one calculates all potentials to the next-to-leading or even higher order.

The main uncertainty comes from the LECs of the contact interactions. Thus the lattice QCD simulation on the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ scattering is called for. Our analytical results can be used to do chiral extrapolation for lattice QCD. With the lattice QCD results in the coming future, the LECs for the contact interaction can be determined more precisely. Then the nature of the $P_{c}$ states in experiment can be identified and more reliable predictions for the other systems can be given.

ACKNOWLEDGMENTS

L. Meng is very grateful to X. Z. Weng, X. L. Chen and W. Z. Deng for very helpful discussions. This project is supported by the National Natural Science Foundation of China under Grants 11575008, 11621131001 and 973 program.

Appendix A Matrix elements

The matrix elements of the isospin-isospin operator is

[TABLE]

where $I$, $I_{1}$ and $I_{2}$ are the total isospin, isospin of $\Sigma_{c}$ and isospin of $\bar{D}^{(*)}$, respectively.

The $\bm{T}$ we defined is proportional to the spin operator $\bm{S_{2}}$ of $\bar{D}^{*}$ as

[TABLE]

The matrix elements of spin-spin operator $\bm{S}_{1}\cdot\bm{S}_{2}$ is

[TABLE]

where $J$ is the total spin of $\Sigma_{c}^{(*)}\bar{D}^{*}$.

The matrix elements of the $\bm{l}_{1}\cdot\bm{l}_{2}$ spin-spin operator in the light degrees of freedom can be calculated via the spin rearrangement. Using the Wigner 9-J symbols, the $\Sigma_{c}^{(*)}\bar{D}^{(*)}$ states can be related to the ones with specific total light spin and heavy spin by the following relations,

[TABLE]

where $l_{i}$ and $h_{i}$ are the light spin and heavy spin for $\Sigma_{c}^{(*)}$ or $\bar{D}^{(*)}$, respectively. $L$ and $H$ are the total light spin and total heavy spin for the two particle states, respectively. Thus, the matrix elements of the $\bm{l}_{1}\cdot\bm{l}_{2}$ can be expressed as

[TABLE]

Appendix B Integrals

B.1 Definitions of integral functions

We will use the “M$x$B$y$” to denote the integrals with $x$ light meson propagators and $y$ heavy hadron propagators in the following.

• •

M1B0

[TABLE]

• •

M2B0

[TABLE]

• •

M1B1

[TABLE]

• •

M2B1

[TABLE]

• •

M2B2

[TABLE]

We use representation $X\vee Y\vee Z...$ to denote the symmetrized tensor structure for concise. For example,

[TABLE]

B.2 Calculations of integral functions

In Ref. Wang et al. (2019), the authors calculated these loop integrals directly. While in this work, we calculate these loop integrals through a different way. We take the $J^{a}_{x}$ as an example. We calculate the $J^{a}_{0}$ directly as in Ref. Wang et al. (2019). Then, we use $v_{\alpha}$ to contract with the $J^{a}_{1x}$ terms, i.e.,

[TABLE]

Here, we get the relation of $J_{11}^{a}$ with the known functions. Third, we can use $g_{\alpha\beta}$ and $v_{\alpha}$ to contact with the $J_{2x}^{a}$ terms,

[TABLE]

[TABLE]

Therefore, following these procedures, we can relate the complicated integrals to the simple ones step by step. The final results we get are equivalent to those in Ref. Wang et al. (2019).

In the following, we give the results of some simple integrals first. The complicated ones can be reexpressed with them,

[TABLE]

where $\bar{\Delta}(z)=m^{2}+q^{2}(z-1)z-i\varepsilon$. When $q^{2}<0$, one get

[TABLE]

where $r=\sqrt{|1-\frac{4m^{2}}{q^{2}}|}$.

[TABLE]

where $\tilde{\Delta}(y)=m^{2}+y^{2}-\omega^{2}-i\varepsilon$. The $J^{T}_{x}$ used in this work reads,

[TABLE]

The results of $J^{B}_{x}(m,\omega_{1},\omega_{2},q)$ depend on whether we remove the 2PR contributions or not. When $\omega_{1}=-\omega_{2}$ there is the pinch singularity. We use the following procedures to remove it,

[TABLE]

where $f(l^{0},\bm{l})$ is the other part of the $J^{B}_{x}$ integrals. In the derivations, the principal integral is used, i.e.,

[TABLE]

This procedure is equivalent to removing the contributions from the poles of the matter fields. When $\omega_{1}\neq-\omega_{2}$, we calculate the $J^{B}_{x}(m,\omega_{1},\omega_{2},q)$ as before. The results read

[TABLE]

where $n$ and $n_{q}$ are the numbers of Lorentz indices and momentum $q$ in the Lorentz structures. $n+n_{q}$ of all the integrals involved in this work are even.