# Exposing the effect of $p$-wave in pion triplet under the strong magnetic
field

Zanbin Xing<EMAIL_ADDRESS>School of Physics, Nankai University,
Tianjin 300071, China  Jingyi Chao<EMAIL_ADDRESS>College of Physics
and Communication Electronics, Jiangxi Normal University, Nanchang 330022,
China  Lei Chang<EMAIL_ADDRESS>School of Physics, Nankai University,
Tianjin 300071, China  Yu-xin Liu<EMAIL_ADDRESS>Department of Physics and
State Key Laboratory of Nuclear Physics and Technology, Peking University,
Beijing 100871, China  Collaborative Innovation Center of Quantum Matter,
Beijing 100871, China  Center for High Energy Physics, Peking University,
Beijing 100871, China

###### Abstract

The static properties, masses and decay constants, of pseudoscalar meson
triplet in a strongly magnetized medium are studied through the Dyson-
Schwinger equation approach treatment of a contact interaction. Complementary
to the usual vector-vector form, a symmetry-preserving formulation of
couplings has been proposed in this work, without modifying the quark
propagator, to control the strength of the $p$-wave component of Bethe-
Salpeter amplitude. It is found that, with the help of flexible auxiliary
interaction, our simple model is able to reproduce the observation in the
lattice QCD simulation, where the spectra of the charged pseudo-scalar meson
shows a non-monotonic behavior as the magnetic field grows. The discovery of
this work implies the strong magnetic field affects the inner structure of
mesons dramatically.

charged pion, weak decay constant, magnetic fields, $p$ wave

##  I introduction

Pion, as the QCD’s Goldstone boson, plays an essential role in low-energy
hadron physics. The deep understanding of pion properties in the vacuum and
under extreme condition might not only provides insight on the emergent hadron
mass but also the phase transition of strong interaction Roberts  _et al._
(2021); Klevansky (1992); Buballa (2005); Brandt  _et al._ (2015); Gao and
Ding (2020). Aside from the extensively explored environments at finite
temperatures and finite baryon densities, magnetic field arises as a new
dimension to analyze the QCD phases and phase transition Agasian and Fedorov
(2008); D’Elia  _et al._ (2010); Mizher  _et al._ (2010); Bali  _et al._
(2012); Andersen  _et al._ (2016). Motivated by the impact on the evolution of
the early universe and the magnetars, where the strength of magnetic field
over its surface is up to $10^{10}\,\mathrm{T}$ Thompson and Beloborodov
(2005); Rea (2013); Kaspi and Beloborodov (2017), and boosted by the novel
topological phenomena emerging in off-central heavy-ion collisions Skokov  _et
al._ (2009); Voronyuk  _et al._ (2011); Deng and Huang (2012), where a
stronger $eB$-field, up to $10^{14}-10^{16}\,\mathrm{T}$, is produced, the
magnetized spectra of pions have been examined in various theoretical
approaches for low energy QCD Fayazbakhsh and Sadooghi (2013); Avancini  _et
al._ (2016); Zhang  _et al._ (2016); Wang and Zhuang (2018); Liu  _et al._
(2018); Mao (2019); Gómez Dumm  _et al._ (2018); Ayala  _et al._ (2018);
Coppola  _et al._ (2019); Sheng  _et al._ (2021); Das and Haque (2020), within
the lattice gauge theory Luschevskaya  _et al._ (2016) and in the chiral
perturbation theory Shushpanov and Smilga (1997); Andersen (2012); Colucci
_et al._ (2014).

While the constituent quarks coupling with the magnetic field, the composite
neutral pion is disturbed by the external magnetic field as a consequence.
However, characterized by the spontaneous chiral symmetry breaking of
$U(1)_{I_{3}}\otimes U(1)_{AI_{3}}\to U(1)_{L+R}$ Chao  _et al._ (2020), it is
expected that the properties of neutral pion is preserved by the Goldstone
nature. Presented by the research groups of Bali  _et al._ (2018a); Ding  _et
al._ (2020), the lattice QCD (LQCD) computations conclude that $\pi_{0}$ is as
light as a pseudo-Goldstone boson. Moreover, an interesting observation is
that the mass of neutral pions monotonically decreases as a function of the
field strength. Some discussions from low-energy approximations of QCD are
given in Ref. Xu  _et al._ (2021). On the contrary, these two non-perturbative
calculations contradict with each other on the energy dispersion of charged
pions. As a tight bound state, the spectra of $\pi_{\pm}$ are perfectly
amenable to the formula of $E^{2}(n,p_{z})=p_{z}^{2}+(2n+1)eB+m_{0}^{2}$ at
weak limit for $eB\leq 0.4\,{\rm GeV}^{2}$. Compared to $\pi_{0}$, the
behaviors of charged mesons are more sensitive to the underlying dynamic of
QCD$\otimes$QED while turning on $eB$ and its point-particle approximation
must become corrupting in large region of the $eB$-field. However, as found in
Ref. Ding  _et al._ (2020), the masses of charged pions and kaons straightly
increase at the beginning, grow slowly, reach to the greatest point, decrease
as the magnetic field increases and end up with a saturation at $eB\gtrsim
2.5\text{GeV}^{2}$. Such a non-monotonic mass curve of charged pions in terms
of the magnetic field is beyond the prediction of low-energy effective
theories Andersen (2012); Colucci  _et al._ (2014); Luschevskaya  _et al._
(2016); Zhang  _et al._ (2016); Wang and Zhuang (2018); Liu  _et al._ (2018);
Mao (2019); Coppola  _et al._ (2019). Meanwhile, Ref. Bali  _et al._ (2012)
shows that, after a same uniform increase and a similarly slower growing, the
mass of charged mesons saturates directly in the ultra-strong limit. Thus,
several theoretical efforts have been done (see, e.g., Refs Xu  _et al._
(2021)) to address the discrepancy between the two first principle results.

Though the pseudoscalar meson is at the state of $J=0^{-+}$, its internal
structure contains a nontrivial $p$-wave component, which has been discussed
in the Bethe–Salpeter equation (BSE) formalism (see, e.g., Refs. Hilger  _et
al._ (2017); Bhagwat  _et al._ (2007)). With regard to this concern, an
intuitive conjecture is that if the higher angular momentum plays an important
role in equipping the Fock states of magnetized pions, it will naturally
reduce the mass of charged mesons. In other words, the isovector meson’s zero-
point energy arising in magnetic fields would be dismissed by the total spin
of two valence quarks like the magnetized vector meson Chernodub (2010). To
demonstrate this idea, in the present paper, we will employ the Dyson-
Schwinger equations (DSE) to study the pions and their properties taking the
advantage that the BSE–DSE includes intrinsically the context of $p$-waves in
the pseudoscalar mesons. As an early investigation, the contact interaction is
adopted. Hence, in order to flexibly control the strength of the $p$-wave, we
blend in an auxiliary term in the interactions, which only manifests itself in
the meson sector. In turn, such an alternative method provides an efficient
way to tune the magnitudes of meson’s scattering kernels.

Furthermore, the scope of this paper will extend to the analysis of the pion-
to-vacuum matrix elements of the axial vector hadronic current based on the
same model approach. Taking into account the energy-momentum tensors of the
magnetic field and numerating all the allowed Lorentz structures, three
leptonic decay constants of charged pions arose. The conventional weak decay
constant $f_{\pi}$ has been analyzed by different points of view Bali  _et
al._ (2018b) and another emerging two have been achieved in the approach of
NJL model Coppola  _et al._ (2019). We note here that we will evaluate only
two decay constants as the functions of the magnetic field since one of them
will vanish in the presented contact model.

This remainder of this paper is organized as follows. In section II, we
describe briefly the full formalism in the sectors of quarks and mesons under
the influence of a uniform magnetic field background and compute the quark
masses, pion masses and weak decay constants with the conventional vector-
vector contact interactions. In section III, we concentrate on the role of
$p$-wave in the meson sector and its impact on the properties of pions and
repeat the calculations containing the auxiliary term. Our conclusions are
summarized in section IV.

##  II general formula in the conventional contact model

###  II.1 Quark Sector

The strong interaction at low energy scale is completely nonperturbative and
extremely difficult by itself and becomes harder to deal with in the presence
of magnetic fields. In this section, we review the model and the formalism
applying in our calculations.

To account for the minimal electromagnetic (EM) coupling along the $z$
direction, the Landau gauge potential $A_{\mu}=\left(0,Bx,0,0\right)$ is used.
For further simplification, a vector-vector contact interaction in coordinate
space is introduced to mimic the strong coupling Gutierrez-Guerrero  _et al._
(2010); Roberts  _et al._ (2011),

$\frac{4}{3m_{G}^{2}}\delta^{4}(a-b)\gamma_{\mu}\otimes\gamma_{\mu},$ (1)

where $a,b$ denote the coordinate $\left(x,y,z,t\right)$. Mass scale of gluon
$m_{G}$ stands for the strength of the instantaneous interaction and is fixed
in order to fit the empirical values. Implementing in the contact interaction,
the quark gap equation develops to

$\displaystyle
S^{-1}\\!\left(a,b\right)=S_{0}^{-1}\\!\left(a,b\right)+\frac{4}{3m_{G}^{2}}\delta^{4}\\!\left(a-b\right)\\!\gamma_{\mu}S(a,b)\gamma_{\mu},$
(2)

where the inverse of the bare quark propagator

$S_{0}^{-1}\\!\left(a,b\right)=\delta^{4}\left(a-b\right)\left(\not{\partial}_{b}+ie_{f}\not{A}+m_{0}\right),$
(3)

and $m_{0}$ is the current quark mass, uniforming in the two-flavor space.
Besides, a general structure of the inverse quark propagator is described as

$S^{-1}(a,b)=\delta^{4}(a-b)\left(\not{\partial}_{b}+ie_{f}\not{A}+M_{f}\right),$
(4)

where the dynamical quark mass $M_{f}$ is a static constant with $f=u,d$. To
obtain the propagator, the usual way is to construct a complete orthogonal
basis and then to form the propagator through it. Belonging to the different
eigenvalues of the Hermitian operator, the magnetized eigenfunctions carry out
a complete set and it allows us apply the operator
$G^{\prime\,2}\left(a,b\right)$ with
$G^{\prime}\left(a,b\right)=\gamma_{5}S^{-1}\left(a,b\right)$, which are
governed by the eigenequations

$\int
d^{4}bG^{\prime\,2}\left(a,b\right)f\left(b\right)=\lambda^{2}f\left(a\right).$
(5)

It is observed that $f(a)$ follows the equation of motion of the simple
harmonic oscillator and the corresponding eigenvalue is double degenerate
excluding at the Lowest Landau Level (LLL).

Employing the basis derived from the Eq. (5), the propagator in
$(x,q_{y},q_{z},q_{t})$ space (we call it as a representation in momentum
space in the rest of the paper even coordinate $x$ remains) is written as

$\displaystyle S\left(\bar{q};x,x^{\prime}\right)=\sum_{n=0}^{\infty}$
$\displaystyle\left\\{P_{+}^{\xi}\frac{-i\not{q}+M_{f}}{q^{2}+M_{n}^{2}}f_{n}\left(x,q_{y}\right)f_{n}\left(x^{\prime},q_{y}\right)-I_{+}^{\xi}\frac{\sqrt{2B_{n}}}{q^{2}+M_{n}^{2}}f_{n-1}\left(x,q_{y}\right)f_{n}\left(x^{\prime},q_{y}\right)\right.$
(6)
$\displaystyle\left.+I_{-}^{\xi}\frac{\sqrt{2B_{n+1}}}{q^{2}+M_{n+1}^{2}}f_{n+1}\left(x,q_{y}\right)f_{n}\left(x^{\prime},q_{y}\right)+P_{-}^{\xi}\frac{-i\not{q}+M_{f}}{q^{2}+M_{n+1}^{2}}f_{n}\left(x,q_{y}\right)f_{n}\left(x^{\prime},q_{y}\right)\right\\},$

where the incomplete-dimension vectors read $\bar{q}=(q_{y},q_{z},q_{t})$ and
$q=(q_{z},q_{t})$, and the full four-dimension vectors are
$\bar{q}_{\mu}=(0,q_{y},q_{z},q_{t})$ and $q_{\mu}=(0,0,q_{z},q_{t})$. Note
that $\xi=\mathrm{sign}(e_{f}B)$, $P_{\pm}^{\xi}=\frac{1\pm
i\xi\gamma_{1}\gamma_{2}}{2}$, $I_{\pm}^{\xi}=\frac{\gamma_{1}\pm
i\xi\gamma_{2}}{2}$, $B_{n}=n\left|e_{f}B\right|$,
$M_{n}^{2}={M_{f}^{2}+2B_{n}}$ and
$f_{n}\left(x,p_{y}\right)=\hat{f}_{n}\left(x+\frac{p_{y}}{e_{f}B}\right)$.
The normalized harmonic oscillator wave functions $\hat{f}_{n}\left(x\right)$
take the forms of

$\displaystyle\hat{f}_{n}\left(x\right)=\frac{\left|e_{f}B\right|^{1/4}}{\sqrt{2^{n}n!\sqrt{\pi}}}H_{n}\left(\sqrt{\left|e_{f}B\right|}x\right)e^{-\left|e_{f}B\right|x^{2}/2},$
(7)

where $H_{n}\left(x\right)$ are the Hermitian polynomials. A detailed
formalism of the propagator $\lim_{a\to a^{\prime}}S(a,a^{\prime})$ can be
found in Appendix A.

Plugging in the magnetized quark propagator, the gap equation of the dynamical
mass transforms to

$\displaystyle
M_{f}=m_{0}+\frac{8M_{f}}{3m_{G}^{2}}\int_{-\infty}^{\infty}\frac{d^{2}q}{\left(2\pi\right)^{2}}\sum_{n=0}^{\infty}\frac{\left|e_{f}B\right|}{2\pi}\left[\frac{1}{q^{2}+M_{n}^{2}}+\frac{1}{q^{2}+M_{n+1}^{2}}\right].$
(8)

The regarded integration and summation are apparently ultraviolet (UV)
divergent. Thus, a regularization scheme must be adopted. In this paper, we
employ the scheme of the proper time regularization Gutierrez-Guerrero  _et
al._ (2010), which is applied as

$\displaystyle\\!\\!\\!\frac{1}{s+M^{2}}=\int_{0}^{\infty}d\tau
e^{-\tau\left(s+M^{2}\right)}\rightarrow\int_{\tau_{uv}^{2}}^{\tau_{ir}^{2}}d\tau
e^{-\tau\left(s+M^{2}\right)}.$ (9)

The UV parameter $\tau_{uv}$ characterizes the dynamical energy scale of all
dimensional quantities in the contact model while as the infrared parameter
$\tau_{ir}$ remains as the description of confinement. The whole used
parameters of the model are

$m_{G}=0.132\,{\rm GeV},\ \tau_{uv}=1/0.905\,{\rm GeV}^{-1},\
\tau_{ir}=1/0.24\,{\rm GeV}^{-1},\ m_{0}=0.0089\,{\rm GeV},$ (10)

which gives $M_{q}=0.37\,{\rm GeV}$ in vacuum. The Landau level is chosen to
be truncated at $n=50$. The numerical results in Fig. 1 demonstrate that
up/down quark mass increases with magnetic field and the quark masses satisfy
a charge relation of $M_{u}(eB)=M_{d}(2eB)$. One can also find out that the
mass plotting presents the well-known magnetic catalysis, which is produced by
the first term of the r.h.s. of the Eq. (8).

Figure 1: Calculated quark masses as the function of magnetic fields in
contact interaction. The solid line denotes the dynamic mass of $u$ quark and
the dashed line represents the mass of $d$ quark.

###  II.2 Meson Sector

Due to the coupling between the magnetic field and charge particles, the
Bethe-Salpeter equation is no longer translation invariant and will be written
in coordinate space. The ladder truncation entails it taking the form

$\displaystyle\Gamma\left(b,b^{\prime}\right)=-\int_{-\infty}^{\infty}d^{4}a\,d^{4}a^{\prime}\,\gamma_{\mu}S_{1}\left(b,a\right)\Gamma\left(a,a^{\prime}\right)S_{2}\left(a^{\prime},b\right)\gamma_{\nu}D_{\mu\nu}(b,b^{\prime}),$
(11)

where
$D_{\mu\nu}(b,b^{\prime})=\frac{4\delta_{\mu\nu}}{3m_{G}^{2}}\delta^{4}(b-b^{\prime})$.
Also, the employed Bethe-Salpeter amplitude (BSA) is
$\Gamma\left(b,b^{\prime}\right)=\Gamma\left(b\right)\delta^{4}(b-b^{\prime})$,
where $\delta^{4}(b-b^{\prime})$ implies a point-like inner structure of
composite particles and $\Gamma\left(b\right)$ describes the motion of
composite degree of freedoms as a whole. Then, the BSE becomes

$\displaystyle\Gamma\left(b\right)=-\frac{4}{3m_{G}^{2}}\int_{-\infty}^{\infty}d^{4}a\,\gamma_{\mu}S_{1}\left(b,a\right)\Gamma\left(a\right)S_{2}\left(a,b\right)\gamma_{\mu}.$
(12)

Converting to momentum space, the BSE is written as

$\displaystyle\Gamma\left(\bar{P};x\right)=-\frac{4}{3m_{G}^{2}}\int_{\bar{q},x^{\prime}}\\!\gamma_{\mu}S_{1}\left(\bar{q},x,x^{\prime}\right)\Gamma\left(\bar{P};x^{\prime}\right)S_{2}\left(\bar{k},x^{\prime},x\right)\gamma_{\mu},$
(13)

where
$\int_{\bar{q},x^{\prime}}=\int_{-\infty}^{\infty}\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}\int_{-\infty}^{\infty}dx^{\prime}$,
$\bar{k}=\bar{q}-\bar{P}$. In general, one has
$S\left(x,x^{\prime}\right)=\phi\left(x,x^{\prime}\right)\tilde{S}\left(\bar{q},x-x^{\prime}\right)$
(see, e.g., Ref. Kojo and Su (2013)), where the position dependence of
magnetic field is attributed to the Schwinger phase factor
$\phi\left(x,x^{\prime}\right)$ and the left $\tilde{S}$ is translation
invariant. Those contributions from the phase factor are trivial for mesons
since only two propagators are involved in the Eq. (13) (a more detailed
explanation can be found in Ref. Chyi  _et al._ (2000)). Hence, Fourier
transforming the translation invariant part and considering the screening mass
along the $z$-direction, it is allowed to use the incomplete-dimension
momentum $\bar{P}=\left(0,P_{z},0\right)$. Finally, the above equation is
turned into an eigenvalue problem in the form

$\displaystyle\lambda(\bar{P}^{2})\Gamma\left(\bar{P};x\right)=-\frac{4}{3m_{G}^{2}}\int_{\bar{q},x^{\prime}}\\!\gamma_{\mu}S_{1}\left(\bar{q},x,x^{\prime}\right)\Gamma\left(\bar{P};x^{\prime}\right)S_{2}\left(\bar{k},x^{\prime},x\right)\gamma_{\mu},$
(14)

and the on-shell condition $m_{\pi}^{2}=-\bar{P}^{2}$ requires

$\lambda(\bar{P}^{2})=1.$ (15)

In contact model, the general structures of pseudoscalar meson’s BSA are given
as

$\Gamma_{\pi}\left(\bar{P};x\right)=\gamma_{5}\left(E_{\pi}\left(\bar{P};x\right)+\frac{-i\not{\bar{P}}}{2M}F_{\pi}\left(\bar{P};x\right)\right),$
(16)

where $M=\frac{M_{1}M_{2}}{M_{1}+M_{2}}$, with $M_{1,2}$ being the effective
mass of propagator $S_{1,2}$. $E$ and $F$ characterize the corresponding
pseudo-scalar and axial-vector components, respectively. Without magnetic
field, the pion mass and leptonic decay constant are
$m_{\pi}(eB=0)=0.157\,{\rm GeV}$, $f_{\pi}(eB=0)=0.102\,{\rm GeV}$ restricting
with both finite $E$ and $F$-terms, while as $m_{\pi}(eB=0)=0.135\,{\rm GeV}$,
$f_{\pi}(eB=0)=0.119\,{\rm GeV}$ with single $E$ component (i.e. setting
$F=0$) in the BSA.

####  II.2.1 Neutral Pion

The neutral pion state
$\left|\pi_{0}\right>=\alpha(B)\left|\pi_{u}\right>-\beta(B)\left|\pi_{d}\right>$
is a mixture of $u,d$ quarks. The associated coefficients satisfy
$\alpha^{2}(B)+\beta^{2}(B)=1$ and
$\alpha(B)=\beta(B)\rightarrow\frac{1}{\sqrt{2}}$ for $B\rightarrow 0$ Bali
_et al._ (2018a). Strictly speaking, the coefficients should be determined by
other quantities which relate to experiment observations of $\pi_{0}$.
However, to keep our discussion simpler, we compute the pure states of
$\left|\pi_{u}\right>$ and $\left|\pi_{d}\right>$ and estimate the real pion
mass locating between them. Here, we will exhibit $\pi_{0}^{u}$ as an example.
We assume that the Bethe-Salpeter amplitude for neutral particles is a delta
function in coordinate space, i.e.
$E_{\pi_{0}^{u}}\left(\bar{P};x\right)=E\delta(x)$ with $E$ being a constant
and so does $F_{\pi_{0}}\left(\bar{P};x\right)$. Since the integration with
respect to $q_{y},\ x$ and $x^{\prime}$ can be done analytically, it reduces
the kernels to

$\displaystyle\left[\begin{array}[]{c}E\\\
F\end{array}\right]=\frac{4}{3m_{G}^{2}}\left[\begin{array}[]{cc}K_{EE}^{\pi_{0}^{u}}&K_{EF}^{\pi_{0}^{u}}\\\
K_{FE}^{\pi_{0}^{u}}&K_{FF}^{\pi_{0}^{u}}\end{array}\right]\left[\begin{array}[]{c}E\\\
F\end{array}\right],$ (17)

where the kernel $K_{EE}^{\pi_{0}^{u}}$ is

$\displaystyle K_{EE}^{\pi_{0}^{u}}$
$\displaystyle=-\hbox{tr}\int_{\bar{q},x}\frac{\gamma_{5}}{4}\gamma_{\mu}S_{u}\left(\bar{q},x,0\right)\gamma_{5}S_{u}\left(\bar{k},0,x\right)\gamma_{\mu}$
(18)
$\displaystyle=\hbox{tr}\int_{\bar{q},x}{\gamma_{5}S_{u}\left(\bar{q},x,0\right)\gamma_{5}S_{u}}\left(\bar{k},0,x\right)$
$\displaystyle=\frac{\left|e_{u}B\right|}{2\pi}2\sum^{\infty}_{n_{u}=0}G_{1,n_{u}}+G_{2,n_{u}}+G_{3,n_{u}}+G_{4,n_{u}}.$

The explicit forms of $G_{j,n_{u}}$ and the other three kernels are listed in
Appendix B.

To manifest the Goldstone nature of neutral pseudoscalar meson in chiral limit
$m_{0}=0$, one writes down the first two elements of the matrix
Eq.$\left(\ref{neutralpionBSE}\right)$, showing as

$\displaystyle
K_{EE}^{\pi_{0}^{u}}=\frac{4}{3m_{G}^{2}}\int\frac{d^{2}q}{\left(2\pi\right)^{2}}\sum^{\infty}_{n_{u}=0}\frac{\left|e_{u}B\right|}{2\pi}\left\\{\frac{2}{q^{2}+M_{n_{u}}^{2}}+\frac{2}{q^{2}+M_{n_{u}+1}^{2}}\right\\},\quad
K_{EF}^{\pi_{0}^{u}}=0.$ (19)

Combining with the gap equation of Eq. (8), it is easy to realize that the
demanded solution of Eq. (15) has been fulfilled automatically by
$m_{\pi}^{2}=-\bar{P}^{2}=0$. Moreover, any mixed state in a two-flavor space
of magnetized neutral pion is massless as required.

To numerically compute the pion mass as a pseudo-Goldstone boson, the scheme
of regularization is applied in the BSE, as well, which is formulated as

$\int\frac{d^{2}q}{(2\pi)^{2}}\frac{q\cdot
k+A}{\left(q^{2}+B\right)\left(k^{2}+C\right)}\
\quad\text{or}\quad\int\frac{d^{2}q}{(2\pi)^{2}}\frac{q\cdot P\ A-k\cdot P\
A^{\prime}}{\left(q^{2}+B\right)\left(k^{2}+C\right)},$ (20)

where $A,A^{\prime},B$ and $C$ are irrelevant to either $q$ or $k=q-P$. The UV
and IR truncation in momentum space are adopted as

$\displaystyle\int\frac{d^{2}q}{(2\pi)^{2}}\frac{q\cdot
k+A}{\left(q^{2}+B\right)\left(k^{2}+C\right)}=\int_{0}^{1}d\alpha\left[\alpha\left(\alpha-1\right)P^{2}+A\right]{\bar{D}}_{1}^{iu}\left[\varsigma\left(\alpha,P^{2},B,C\right)\right]+D^{iu}\left[\varsigma\left(\alpha,P^{2},B,C\right)\right]-D_{1}^{iu}\left[\varsigma\left(\alpha,P^{2},B,C\right)\right],$
(21) $\displaystyle\int\frac{d^{2}q}{(2\pi)^{2}}\frac{q\cdot P\ A-k\cdot P\
A^{\prime}}{\left(q^{2}+B\right)\left(k^{2}+C\right)}=P^{2}\int_{0}^{1}d\alpha\left[\alpha\
A-\left(\alpha-1\right)\
A^{\prime}\right]{\bar{D}}_{1}^{iu}\left[\varsigma\left(\alpha,P^{2},B,C\right)\right],$

where

$\displaystyle
D^{iu}\left[\varsigma\right]=\int_{\tau_{uv}^{2}}^{\tau_{ir}^{2}}d\tau\int\frac{d^{2}q}{(2\pi)^{2}}e^{-\tau\left(q^{2}+\varsigma\right)},\
{\bar{D}}_{1}^{iu}\left[\varsigma\right]=-\frac{d}{d\varsigma}D^{iu}\left[\varsigma\right],\
D_{1}^{iu}\left[\varsigma\right]=\varsigma{\bar{D}}_{1}^{iu}\left[\varsigma\right],$
(22)
$\displaystyle\varsigma\left(\alpha,P^{2},B,C\right)=\alpha\left(1-\alpha\right)P^{2}+\left(1-\alpha\right)B+\alpha
C.$

The obtained neutral pion masses are shown in Fig. 2. While the BSA
incorporates both $E$ and $F$, the mass first decreases and then increases as
the magnetic field grows. In the context of single $E$, the mass first
decreases, but then saturates. Obviously, the component of $F$ plays a non-
negligible role. It is also observed that an individual $E$ term in the BSA
renders a better result to accommodate with the computations of LQCD Ding  _et
al._ (2020); Bali  _et al._ (2018a), which is artificial and a subtler
consideration will show in the next section.

Figure 2: Calculated results of the mass of neutral pion $m_{\pi_{0}^{u}}$ as
the function of the magnetic field in the conventional contact interaction
model. The solid line denotes the mass with complete BSA incorporating both
$E$ and $F$, the dashed line denotes the result extracted from an individual
$E$ term in the BSA.

####  II.2.2 Charged Pion

Unlike the neutral pion, $\pi_{\pm}$ are no longer protected by the Goldstone
theorem and become massive under the influence of external magnetic fields.
Also, the Bethe-Salpeter amplitude is not as simple as a delta function, but
formed by a more complicated Gaussian-like structure. For $\pi_{+}$, the
constituent quarks are $u\bar{d}$ and the Bethe-Salpeter equation is

$\displaystyle\left[\begin{array}[]{c}E_{\pi_{+}}\\\
F_{\pi_{+}}\end{array}\right]=\frac{4}{3m_{G}^{2}}\int
dx^{\prime}\left[\begin{array}[]{cc}K_{EE^{\prime}}^{\pi_{+}}&K_{EF^{\prime}}^{\pi_{+}}\\\
K_{FE^{\prime}}^{\pi_{+}}&K_{FF^{\prime}}^{\pi_{+}}\end{array}\right]\left[\begin{array}[]{c}E^{\prime}_{\pi_{+}}\\\
F^{\prime}_{\pi_{+}}\end{array}\right],$ (23)

where $(E,F)$ and $(E^{\prime},F^{\prime})$ are functions of
$\left(\bar{P};x\right)$ and $\left(\bar{P};x^{\prime}\right)$, respectively.
The detail form of kernel $K_{EE^{\prime}}^{\pi_{+}}$ is

$\displaystyle K_{EE^{\prime}}^{\pi_{+}}$
$\displaystyle=-\hbox{tr}\frac{\gamma_{5}}{4}\int_{\bar{q}}\gamma_{\mu}S_{u}\left(\bar{q},x,x^{\prime}\right)\gamma_{5}S_{d}\left(\bar{k},x^{\prime},x\right)\gamma_{\mu}$
(24)
$\displaystyle=\hbox{tr}\int_{\bar{q}}{\gamma_{5}S_{u}\left(\bar{q},x,x^{\prime}\right)\gamma_{5}S_{d}}\left(\bar{k},x^{\prime},x\right)$
$\displaystyle=2\sum^{\infty}_{n_{u}=0}\sum^{\infty}_{n_{d}=0}G_{1,n_{u},n_{d}}g_{1,n_{u},n_{d}}+G_{2,n_{u},n_{d}}g_{2,n_{u},n_{d}}+G_{3,n_{u},n_{d}}g_{3,n_{u},n_{d}}+G_{4,n_{u},n_{d}}g_{4,n_{u},n_{d}}.$

Here $G_{j,n_{u},n_{d}}$ and $g_{j,n_{u},n_{d}}$ have integrated with respect
to $(q_{z},q_{t})$ and $q_{y}$, respectively. The coordinate dependence of $x$
and $x^{\prime}$ remains in $g_{j,n_{u},n_{d}}$. The UV regularization scheme
has been applied to obtain the functions of $G_{j,n_{u},n_{d}}$ and their
explicit forms are shown in Appendix B.

As shown in Fig. 3, for BSA incorporating both $E$ and $F$, $\pi_{+}$ mass
increases rapidly and the solution of the mesonic bound state disappears for
$eB>0.8\,\mathrm{GeV}^{2}$ in the traditional model of contact interaction.
Such strong enhancement of mass and the disappearance of the bound state are
produced by the incorrect ansatz of the $F$ term. Indeed, for a single $E$
structure adopting in the BSA, the spectra of $\pi_{+}$ is very close to the
result of the lowest Landau level (LLL) approximation. As we mentioned before,
a careful analysis of the $F$-component is necessary and will be addressed in
the next section.

Figure 3: Calculated results of the charged pion mass as the function of the
magnetic field in the conventional contact interaction model. The solid line
denotes the mass of $\pi_{+}$ with complete BSA incorporating both $E$ and
$F$, the dashed line denotes the result extracting from an individual $E$ term
in the BSA and the dotted line represents the LLL approximation.

###  II.3 Weak Decay Constants

According to the equivalence theorem, the Yukawa coupling in pseudo-scalar
channels can be rewritten to pseudo-vector types via the integration by part:

$\pi\bar{\psi}\gamma_{5}\psi\to g^{\mu\nu}\partial_{\mu}\pi\,J_{5\nu},$ (25)

where the axial current $J_{5\nu}=\bar{\psi}\gamma_{5}\gamma_{\nu}\psi$.
Within EM-field, as a consequence, two more weak decay constants are required
to describe the one-pion-to-vacuum matrix elements by virtue of the couplings
of $\left(F,\tilde{F}\right)^{\mu\nu}\partial_{\mu}\pi\,J_{5\nu}$. Here, those
three decay constants $f_{\pi},\ f^{\prime}_{\pi},\ f^{\prime\prime}_{\pi}$
are defined as:

$f_{\pi}\bar{P}_{\mu}+f^{\prime\prime}_{\pi}eF_{\mu\nu}\bar{P}^{\nu}=N_{c}tr_{D}\int_{x,x^{\prime},\bar{q}}i\gamma_{5}\gamma_{\mu}S(\bar{q};x,x^{\prime})\Gamma_{\pi}(\bar{P};x^{\prime})S(\bar{q}-\bar{P};x^{\prime},x);$
(26)
$f^{\prime}_{\pi}\frac{i}{2}\epsilon_{\mu\nu\rho\sigma}eF^{\nu\rho}\bar{P}^{\sigma}=N_{c}tr_{D}\int_{x,x^{\prime},\bar{q}}-i\gamma_{\mu}S(\bar{q};x,x^{\prime})\Gamma_{\pi}(\bar{P};x^{\prime})S(\bar{q}-\bar{P};x^{\prime},x).$
(27)

Since the Schwinger phase factor doesn’t contribute in our simplified vertex
of $\Gamma_{\pi}$, i.e., two phase factors are involved in the above formulas,
we have replaced the derivative of position $\mathop{}\\!i\partial_{\mu}$ to
$\bar{P}_{\mu}$. In our case, for $\bar{P}_{\mu}=(0,0,P_{z},0)$ and
$F_{12}=-F_{21}=B$, the combination $F_{\mu\nu}\bar{P}_{\nu}$ vanishes
automatically. We therefore examine $f$ and $f^{\prime}$ in the present paper.
Note here that the normalization of BSA is applied by:

$\left[\frac{\partial\ln\lambda(\bar{P}^{2})}{\partial\bar{P}^{2}}\right]^{-1}=-2N_{c}tr_{D}\int_{\bar{q},x,x^{\prime},x^{\prime\prime}}\bar{\Gamma}_{\pi}(-\bar{P};x-x^{\prime\prime})S(\bar{q};x,x^{\prime})\Gamma_{\pi}(\bar{P};x^{\prime})S(\bar{q}-\bar{P};x^{\prime},x).$
(28)

The obtained results of the magnetic field dependence of the charged pion
decay constants are shown in Fig. 4. Constructing the BSA in the contexts of
both $E$ and $F$, $f$ and $f^{\prime}$ increases as the magnetic field
increases. The trend is slower while a single $E$ component is involved.

Figure 4: Calculated results of the charged pion decay constants as the
function of the magnetic field in the conventional contact interaction model.
$f(eB)$ are presented in the upper panel and $f^{\prime}(eB)eB$ are presented
in the lower panel. The solid lines denote the corresponding decay constant
computed by a complete BSA incorporating both $E$\- and $F$-terms; the dashed
lines represent the ones in terms of a single $E$ component.

##  III With an auxiliary term in the contact model

While a single $E$ component of the BSA is applied, the form of the quark-
meson vertex is represented by the bare $\gamma_{5}$ in calculation, which
obviously is a bad approximation. It is known that the Bethe-Salpeter (BS)
wave function is related to the BSA according to

$\chi_{\pi}(P)=S(q_{+})\Gamma(P)S(q_{-}),$ (29)

where $q_{+}=q+\sigma P$ and $q_{-}=q-(1-\sigma)P$ with $0<\sigma<1$. Read
from Eq. (16), the term of $E$ in the BSA contributes as the $s$-wave in the
BS wave function while as the axial-vector component of $F$ is attributed to
both $s$\- and $p$-wave functions, discussed in the referencesHilger  _et al._
(2017); Bhagwat  _et al._ (2007). Therefore, within Dyson-Schwinger formalism,
this allowed us to investigate the internal structures of the meson modes.

Adopting the conventional form of $F$, it has been found that the DSE model
calculation fails to present the tight meson as a bound state starting from
the moderate magnetic field region. A detailed calculation of pseudo-scalar
and vector mesons can be found in the thesis work of Ref. Wang (2013). To
address the failure found in that work, we propose to take full account of the
back reaction on the gluon sector under the influence of the magnetic field
and hence use a more suitable ansatz of $F$-term in BSA. As explored in Ref.
Chao  _et al._ (2014), the Lorentz structures of gluon propagator have been
split into three parts in the magnetized vacuum, which takes the form

$\displaystyle
D^{\mu\nu}(k)=\frac{P_{\shortparallel}(k)}{k^{2}-m^{2}_{G,\shortparallel}}+\frac{P_{\perp}(k)}{k^{2}-m^{2}_{G,\perp}}++\frac{k^{\mu}k^{\nu}/k^{2}-g^{\mu\nu}-P_{\shortparallel}(k)-P_{\perp}(k)}{k^{2}-m^{2}_{G,\mathrm{mix}}},$
(30)

where
$P_{\shortparallel}(k)=k_{\shortparallel}^{\mu}k_{\shortparallel}^{\nu}/k_{\shortparallel}^{2}-g^{\mu\nu}_{\shortparallel}$
and
$P_{\perp}(k)=k_{\perp}^{\mu}k_{\perp}^{\nu}/k_{\perp}^{2}-g^{\mu\nu}_{\perp}$.
Note here that the metric convention $g^{\mu\nu}$ is decomposed into two
orthogonal sub-spaces $g^{\mu\nu}_{\shortparallel}=\left(1,0,0,-1\right)$ and
$g^{\mu\nu}_{\perp}=\left(0,-1,-1,0\right)$. Similar decomposition is adopted
for four dimensional momentum $k_{\mu}=k_{\shortparallel}+k_{\perp}$,
$k_{\shortparallel}=\left(k_{0},0,0,k_{3}\right)$ and
$k_{\perp}=\left(0,k_{1},k_{2},0\right)$. It shows that
$m^{2}_{G,\shortparallel}\neq m^{2}_{G,\perp}\neq m^{2}_{G,\mathrm{mix}}$ at
finite $eB$. As a consequence, the effect of $p$-wave in mesons is modified
regarding the inclusion of magnetic field. To account for the splitting of the
gluonic degrees of freedom, we suggest the following replacement of quark
propagator

$\displaystyle S^{-1}\\!\left(a,b\right)$
$\displaystyle=S_{0}^{-1}\\!\left(a,b\right)+\frac{4}{3m_{G}^{2}}\delta^{4}\\!\left(a-b\right)\\!\gamma_{\mu}S(a,b)\gamma_{\mu}$
(31)
$\displaystyle-\tilde{\eta}^{2}\frac{8}{3m_{G}^{2}}\delta^{4}\\!\left(a-b\right)\\!(\gamma_{\mu}S(a,b)\gamma_{\mu}+\gamma_{5}\gamma_{\mu}S(a,b)\gamma_{5}\gamma_{\mu}),$

where the explicit form of $\tilde{\eta}$ will be written down in the next
context. Specifically, as an early investigation, we require that the second
line in the interaction kernel does not affect the early discussed properties
of quarks in the present work.

To maintain Ward-Takahashi identity (WTI), the corresponding BSE is of the
form

$\displaystyle\Gamma\left(\bar{P};x\right)$
$\displaystyle=-\frac{4}{3m_{G}^{2}}\int_{\bar{q},x^{\prime}}\\!\gamma_{\mu}S_{1}\left(\bar{q},x,x^{\prime}\right)\Gamma\left(\bar{P};x^{\prime}\right)S_{2}\left(\bar{k},x^{\prime},x\right)\gamma_{\mu}$
(32)
$\displaystyle+\tilde{\eta}^{2}\frac{8}{3m_{G}^{2}}\int_{\bar{q},x^{\prime}}\\!(\gamma_{\mu}S_{1}\left(\bar{q},x,x^{\prime}\right)\Gamma\left(\bar{P};x^{\prime}\right)S_{2}\left(\bar{k},x^{\prime},x\right)\gamma_{\mu}+\gamma_{5}\gamma_{\mu}S_{1}\left(\bar{q},x,x^{\prime}\right)\Gamma\left(\bar{P};x^{\prime}\right)S_{2}\left(\bar{k},x^{\prime},x\right)\gamma_{5}\gamma_{\mu}).$

The auxiliary term controls the strength of the $F$-term by the function of
$\tilde{\eta}^{2}$, serving as a tool to exhibit the role of $p$-wave. To
achieve our purpose, the magnitude of $\tilde{\eta}^{2}$ is determined in an
empirical way: (i) since the auxiliary term is introduced as the excitation of
magnetic fields, it is natural to require $\tilde{\eta}^{2}(eB=0)=0$; (ii) to
maintain the spectra of mesons in vacuum, we let
$\tilde{\eta}_{0}^{2}=\frac{1}{4}$ to ensure that $F$-term vanishes at zero
magnetic field, which is proportional to $1-4\tilde{\eta}^{2}$ as shown in the
following. It also means we assume that the contribution of $p$-wave is too
tiny to remain in vacuum; (iii) the form of $\tilde{\eta}^{2}$ is equivalent
to be brought by a one-meson exchange contact interaction and it has to be
proportional to even powers of $\bar{P}$ to avoid affecting the Goldstone
nature in chiral limit; (iv) after an exploration in the parameter space, we
found that, to fully fit the lattice QCD results in Ref. Ding  _et al._
(2020), one has

$\tilde{\eta}^{2}=\tilde{\eta}_{0}^{2}-a*\bar{P}^{2}*\tau_{uv}^{2}*\tanh(b*eB^{2}*\tau_{uv}^{4}),$
(33)

where $a=3.5,\ b=0.153$ and $\tau_{uv}$ is introduced as the energy scale.

Plugging into $\tilde{\eta}^{2}$, the BSA is expressed as

$\displaystyle\left[\begin{array}[]{c}E\\\
F\end{array}\right]=\frac{4}{3m_{G}^{2}}\left[\begin{array}[]{cc}1&0\\\
0&1-4\tilde{\eta}^{2}\end{array}\right]\left[\begin{array}[]{cc}K_{EE}^{\pi_{0}^{u}}&K_{EF}^{\pi_{0}^{u}}\\\
K_{FE}^{\pi_{0}^{u}}&K_{FF}^{\pi_{0}^{u}}\end{array}\right]\left[\begin{array}[]{c}E\\\
F\end{array}\right].$ (34)

As noted above, a factor of $(1-4\tilde{\eta}^{2})$ arises in the second line
of the kernel matrix. The elements of
$K_{EE}^{\pi_{0}^{u}},K_{EF}^{\pi_{0}^{u}},K_{FE}^{\pi_{0}^{u}}$ and
$K_{FF}^{\pi_{0}^{u}}$ are exactly same as those in Eq.
$(\ref{neutralpionBSE})$.

After solving Eq. (34), we have the spectra of neutral pions as a function of
the strength of the magnetic field as shown in Fig. 5. It is observed that the
mass of $\pi_{0}$ decreases and then tends to saturate, which is roughly in
accordance with both the LQCD simulations in Refs. Ding  _et al._ (2020); Bali
_et al._ (2018a). Such trend of $\pi_{0}$ is originated from that
$1-4\tilde{\eta}^{2}$ limits to zero as $\bar{P}^{2}\sim 0$, which diminishes
the role of $F$-component and ceases the incorrect tendency of the mass
increasing, which was shown in Fig. 2.

Figure 5: Calculated results of neutral pion mass $m_{\pi_{0}^{u}}$ as the
function of the magnetic field in the improved contact interaction model. The
solid line is computed with $\tilde{\eta}^{2}$ in Eq. (33); the dashed line
are lattice computations from Ref. Bali  _et al._ (2018a); the dotted line are
lattice computations from Ref. Ding  _et al._ (2020). Figure 6: Calculated
ratio of $E/F$ as the function of the magnetic field with different model
parameters. The solid line and dashed line present the results of the full BSA
and single $E$, respectively, in the conventional contact model. The dotted
line denotes the results in the improved contact interaction model. Here, the
general form of the ground state in the BSA is $\sim e^{-|eB|x^{2}/2}$ for a
charged pion, and the line is plotted with the strength ratio between the $E$
and $F$ components.

Similarly, we write down the BSE of charged pions after implementing into
$\tilde{\eta}^{2}$,

$\displaystyle\left[\begin{array}[]{c}E_{\pi_{+}}\\\
F_{\pi_{+}}\end{array}\right]=\frac{4}{3m_{G}^{2}}\int
dx^{\prime}\left[\begin{array}[]{cc}1&0\\\
0&1-4\tilde{\eta}^{2}\end{array}\right]\left[\begin{array}[]{cc}K_{EE^{\prime}}^{\pi_{+}}&K_{EF^{\prime}}^{\pi_{+}}\\\
K_{FE^{\prime}}^{\pi_{+}}&K_{FF^{\prime}}^{\pi_{+}}\end{array}\right]\left[\begin{array}[]{c}E^{\prime}_{\pi_{+}}\\\
F^{\prime}_{\pi_{+}}\end{array}\right].$ (35)

Here, the kernels $K_{EE^{\prime}}^{\pi_{+}},\ K_{EF^{\prime}}^{\pi_{+}},\
K_{FE^{\prime}}^{\pi_{+}},\ K_{FF^{\prime}}^{\pi_{+}}$ are exactly the same
with those in Eq. $(\ref{chargedpionBSE})$.

Under the influence of $\tilde{\eta}^{2}$, one finds out that the relative
sign of $E/F$ has been flipped to negative compared with the conventional
approach. It is known that the $p$-wave forms as a triplet-state and those
three behave the same at finite temperatures or densities. However, the energy
states of the charged vector mesons are spin dependent, called the Zeeman
effect, and only one of $s_{z}$ reduces the energy dispersion for spin-one
meson modes. We numerically confirm that the demanded $p$-wave component is
attributed from the negative $F$-term, where the ratio of $E/F$ is shown in
Fig. 6. As discussed in Sec. II, the term of $E$ plays a dominant role in the
weak limit and the contribution of $F$-structure is enhanced in moderate and
strong regimes of magnetic fields.

The final results of charged pion masses are plotted in Fig. 7. Here, we get a
non-monotonic curve as the function of magnetic fields by adopting a non-
trivial $\tilde{\eta}^{2}$. One observes that the spectra of the charged pion
grows and fits the point-particle description in the weak strength regime of
$eB$-field; the line bends at the greatest point, around $eB\sim 0.6~{}{\rm
GeV}^{2}$; and then decreases as the magnetic field increases, which agrees
with the LQCD simulation result in Ref. Ding  _et al._ (2020). We also present
the weak decay constant $f(eB)$ in Fig. 8. It shows a similar non-monotonic
behavior like the mass, which is consistent in the present model approach.
Another axial decay constant $f^{\prime}(eB)eB$ saturates while $eB>0.6~{}{\rm
GeV}^{2}$, which is mainly induced by the crude contact interactions. A
sophisticated investigation is in progress.

Figure 7: Calculated results of charged pion mass as the function of the
magnetic field in the improved contact interaction model. The solid line is
computed with $\tilde{\eta}^{2}$ in Eq. (33); the dotted line denote the
result of LLL approximation and the dashed line is the lattice result, given
in Ref. Ding  _et al._ (2020). Figure 8: Calculated results of the weak decay
constant as the function of the magnetic field in the improved contact
interaction model. $f(eB)$ are presented in the upper panel and
$f^{\prime}(eB)eB$ are presented in the lower panel. The solid line is
computed with $\tilde{\eta}^{2}$ in Eq. (33) and the band is the lattice
result, given in Ref. Bali  _et al._ (2018b)

##  IV Summary

In the present paper, we have systematically studied the pion properties in
the magnetized QCD matter under the framework of the Dyson-Schwinger equation.
Employing the contact model, we have discussed the spectra of dynamical quarks
as well as pion triplets since these are the fundamental properties of hadron
states. Created by the metric $g_{\mu\nu}$ and the additional tensor
structures $F_{\mu\nu}$ and $\tilde{F}_{\mu\nu}$, the relevant leptonic decay
modes of mesons are three. We have calculated the $f$, relating to
$g_{\mu\nu}$ and $f^{\prime}$, coming with $\tilde{F}_{\mu\nu}$, as the
functions of $eB$, while as $f^{\prime\prime}$, constructing with
$F_{\mu\nu}$, vanishes because of the simplified ansatz of $\Gamma_{\pi}$.

In the context of vector-vector contact interaction,
$\gamma_{\mu}\otimes\gamma_{\mu}$, the constituent quark mass increases as the
magnetic field grows and there is no clue that the charged pion spectra will
decrease under the strong magnetic fields. Inspired by the Zeeman splitting in
the energy dispersion of charged vector mesons, we suggested that an enhanced
$p$-wave state in the axial-scalar meson will reduce the pole mass, making it
as light as $\rho^{\pm}$ with spin component $s_{z}=\mp 1$. To demonstrate our
conjecture, we have provided a simple method to stimulate the higher angular
momentum contribution in contact model by implementing an auxiliary term,
$\gamma_{\mu}\otimes\gamma_{\mu}-\gamma_{\mu}\gamma_{5}\otimes\gamma_{\mu}\gamma_{5}$,
into the interactions. In particular, such coupling does not affect the quark
sector, but is highly involved in the Bethe-Salpeter amplitude of pions.

Before our model computations, as comparisons, we have applied a complete
Bethe-Salpeter amplitude to show the quantities of pions without including the
auxiliary term. We have obtained that the spectra of $\pi_{0}$ increases with
magnetic field for $eB>0.5\,\mathrm{GeV}^{2}$ while the charged pion mass
grows too quickly to achieve a bound state for $eB>0.8\,\mathrm{GeV}^{2}$.
These behaviors indicate that, due to the presence of the magnetic field, the
degeneracy of the $p$-wave has been broken and the sign of the $F$-component
has become a significant factor to characterize the spin direction. Meanwhile,
we have found out that both $f$ and $f^{\prime}|eB|$ increase with magnetic
field for the specified BSA. We have also calculated these quantities by
artificially let $F_{\pi}\left(\bar{P};x\right)=0$ and observed that the BSA
reduces to
$\Gamma_{\pi}\left(\bar{P};x\right)=\gamma_{5}E_{\pi}\left(\bar{P};x\right)$,
manifesting as a pure $s$-wave state and very closing to the calculations of
LLL approximation. The resulting neutral pion mass first decreases and then
tends to saturate for $eB>1.5\,\mathrm{GeV}^{2}$. The energy dispersion of
$\pi^{\pm}$ is slightly larger than the point-particle approximation.
Moreover, the two remaining decay constants still grow along the magnetic
field, but much more slowly than the complete BSA case.

In order to control the strength of the $p$-wave state and fit the LQCD
results in Ref. Ding  _et al._ (2020), we have specified the value of
$\tilde{\eta}^{2}$. Not surprising, the neutral pion spectra of present work
are in accordance with the first principle simulations results Bali  _et al._
(2018a); Ding  _et al._ (2020). Besides, we have obtained that the mass of
charged pions bends around $eB\sim 0.8\,\mathrm{GeV}^{2}$ and becomes flatting
in the ultra-strong magnetic fields, which reproduces the curve in Ref. Ding
_et al._ (2020). In our model calculations, the charged pion decay constant
$f_{\pi}(eB)/f_{\pi}(0)$ behaves like the trend of mass and does not deviate
from $1$ too much after $eB>0.3\,\mathrm{GeV}^{2}$ while as
$f^{\prime}_{\pi}(eB)|eB|/f_{\pi}(0)$ increases straightly in the weak limit
and then saturates for $eB>0.3\,\mathrm{GeV}^{2}$. These curves of decay
constants don’t agree well with the lattice results in Ref Bali  _et al._
(2018b). Indeed, at large magnetic field regions, the artificial saturation of
the decay constants is inherited from the employed contact interactions.
Further research is needed to investigate the strongly magnetized medium.

Apparently, the presently employed model is too simple to manifest the full
physics subject to an external magnetic field. The form of $\tilde{\eta}^{2}$
is not derived from an ab initio approach, but adopted through a
phenomenological way. A complete representation of Poincare covariant
components in the Bethe-Salpeter equation hasn’t been fulfilled yet. The main
aim of the current work, as a first attempt, is to serve as a model to exhibit
the effect of the $p$-wave. Therefore, it is still too early to tell the
discrepancy between the two non-perturbatively computations while a
sophisticated description of the higher angular momentum state under the
effect of the magnetic field hasn’t been achieved.

Last but not least, even though the lifetime of magnetic fields in off-central
heavy-ion is too short to result in any observations. Our analysis of weak
decay constants will influence the astrophysical environment and the extension
at finite temperatures will be explored in the future.

###### Acknowledgements.

We thank H. Ding for sharing their pion masses data and B.B. Brandt, G. Endödi
for sharing their pion masses and weak decay constants data. This work is
supported by the National Natural Science Foundation of China (NSFC) Grant
Nos. 11435001, 11775041 and 12135007. The work of JC has been supported by the
start-up funding from Jiangxi Normal University under grant No. 12021211.

##  Appendix A Quark propagator

In this section, we define $P_{\pm}=P_{\pm}^{\xi=1}$ and
$I_{\pm}=I_{\pm}^{\xi=1}$, and the flavour dependence is suppressed for
convenience. The explicit form of the Hermitian operator
$G_{0}^{\prime\,2}\left(a,b\right)$ are

$\displaystyle
G_{0}^{\prime\,2}\left(a,b\right)=\delta^{4}(a-b)\left[-\partial_{b}^{2}-2ie_{f}A\cdot\partial-
ie_{f}(\partial_{\mu}A_{\nu})\gamma_{\mu}\gamma_{\nu}+e_{f}^{2}A^{2}+M^{2}\right].$
(36)

By solving Eq. (5), we can obtain the eigenfunctions in
$(x,p_{y},p_{z},p_{t})$ space

$\displaystyle
P_{\pm}f_{n,\bar{p}}\left(b\right)=P_{\pm}\hat{f_{n}}\left(b_{x}+\frac{p_{y}}{eB}\right)e^{\left(i\bar{b}\cdot\bar{p}\right)},$
(37)

where $\bar{p}=(p_{y},p_{z},p_{t}),\bar{b}=(b_{y},b_{z},b_{t})$. The Fourier
factor $e^{\left(i\bar{b}\cdot\bar{p}\right)}$ indicates the translation
invariant $(y,z,t)$ directions and $\hat{f}_{n}\left(x\right)$ is the
normalized harmonic oscillator wave function in the Eq. (7). It is convenient
to introduce ladder operators

$\displaystyle\hat{a}_{+}=\frac{\left|eB\right|x-\partial_{x}}{\sqrt{2\left|eB\right|}},\quad\hat{a}_{-}=\frac{\left|eB\right|x+\partial_{x}}{\sqrt{2\left|eB\right|}},$
(38)

where

$\displaystyle\hat{a}_{+}\hat{f}_{n}\left(x\right)=\sqrt{n+1}\hat{f}_{n+1}\left(x\right),\quad\hat{a}_{-}\hat{f}_{n}\left(x\right)=\sqrt{n}\hat{f}_{n-1}\left(x\right).$
(39)

Forming as a complete orthogonal basis, we use the eigenfunctions (37) to
expand the inverse of the dressed propagator (4). The expansion coefficients
are

$\int
d^{4}bS^{-1}\left(a,b\right)P_{\pm}f_{n}\left(b\right)=\left(i\not{p}+M+I_{+}\left(\partial_{\hat{x}}+eB\hat{x}\right)+I_{-}\left(\partial_{\hat{x}}-eB\hat{x}\right)\right)P_{\pm}f_{n}\left(a\right),$
(40)

where ${\hat{x}}=a_{x}+\frac{p_{y}}{eB}$, $p=(0,0,p_{z},p_{t})$,
$I_{\pm}=I_{\pm}^{\xi=1}$ and $I_{\pm}^{\xi}=\frac{\gamma_{1}\pm
i\xi\gamma_{2}}{2}$. $f_{n,\bar{p}}\left(a\right)$ is written as
$f_{n}\left(a\right)$ for convenience. For $eB>0$:

$\int
dbS^{-1}\left(a,b\right)P_{+}f_{n}\left(b\right)=P_{+}\left(i\not{p}+M\right)f_{n}\left(a\right)+I_{+}\sqrt{2B_{n}}f_{n-1}\left(a\right),$
(41) $\int
dbS^{-1}\left(a,b\right)P_{-}f_{n}\left(b\right)=P_{-}\left(i\not{p}+M\right)f_{n}\left(a\right)-I_{-}\sqrt{2B_{n+1}}f_{n+1}\left(a\right).$
(42)

Using $\int
d^{4}bd^{4}cS^{-1}\left(a,b\right)S\left(b,c\right)P_{\pm}f_{n}\left(c\right)=P_{\pm}f_{n}\left(a\right)$,
the expansion coefficients of the dressed propagator are

$\displaystyle\int dbS\left(a,b\right)$ $\displaystyle
P_{+}f_{n}\left(b\right)=\left(P_{+}\left(-i\not{p}+M\right)f_{n}\left(a\right)-I_{+}\sqrt{2B_{n}}f_{n-1}\left(a\right)\right)\frac{1}{p^{2}+M^{2}+2B_{n}},$
(43) $\displaystyle\int dbS\left(a,b\right)$ $\displaystyle
P_{-}f_{n}\left(b\right)=\left(P_{-}\left(-i\not{p}+M\right)f_{n}\left(a\right)+I_{-}\sqrt{2B_{n+1}}f_{n+1}\left(a\right)\right)\frac{1}{p^{2}+M^{2}+2B_{n+1}}.$
(44)

Similarly, we obtain the expansion coefficients of the propagator for $eB<0$.
We summarize the propagator with any $eB$ as

$\displaystyle S\left(\bar{p};x,x^{\prime}\right)=$
$\displaystyle\sum_{n=0}^{\infty}{\left(P_{+}^{\xi}\left(-i\not{p}+M\right)f_{n}\left(x,p_{y}\right)-I_{+}^{\xi}\sqrt{2B_{n}}f_{n-1}\left(x,p_{y}\right)\right)\frac{1}{p^{2}+M^{2}+2B_{n}}}P_{+}^{\xi}f_{n}\left(x^{\prime},p_{y}\right)$
(45) $\displaystyle+$
$\displaystyle\left(P_{-}^{\xi}\left(-i\not{p}+M\right)f_{n}\left(x,p_{y}\right)+I_{-}^{\xi}\sqrt{2B_{n+1}}f_{n+1}\left(x,p_{y}\right)\right)\frac{1}{p^{2}+M^{2}+2B_{n+1}}P_{-}^{\xi}f_{n}\left(x^{\prime},p_{y}\right).$

Now $|eB|$ appears everywhere except in the harmonic wave function
$f_{n}\left(x,p_{y}\right)=\hat{f}_{n}\left(x+\frac{p_{y}}{eB}\right)$.
Finally, for the dynamical quark mass

$\displaystyle M$
$\displaystyle=m_{0}+\frac{4}{3m_{G}^{2}}\gamma_{\mu}S(a,a)\gamma_{\mu},$ (46)

the propagator $S(a,a)$ is explicitly expressed as

$\displaystyle S\left(a,a\right)$
$\displaystyle=\int_{\bar{p}}\sum_{n=0}^{\infty}{\left(\frac{P_{+}^{\xi}\left(-i\not{p}+M\right)}{p^{2}+M^{2}+2B_{n}}f_{n}^{2}\left(a_{x},p_{y}\right)+\frac{P_{-}^{\xi}\left(-i\not{p}+M\right)}{p^{2}+M^{2}+2B_{n+1}}f_{n}^{2}\left(a_{x},p_{y}\right)\right)}e^{i\bar{p}\cdot(\bar{a}-\bar{a})}$
(47)
$\displaystyle=\int_{p}\sum_{n=0}^{\infty}{\frac{\left|eB\right|}{2\pi}\left(\frac{P_{+}^{\xi}\left(-i\not{p}+M\right)}{p^{2}+M^{2}+2B_{n}}+\frac{P_{-}^{\xi}\left(-i\not{p}+M\right)}{p^{2}+M^{2}+2B_{n+1}}\right)}.$

##  Appendix B Meson properties

###  B.1 neutral pion

The kernels in Eq. (17) are formed as

$\displaystyle K_{EE}^{\pi_{0}^{u}}$ $\displaystyle=-tr\int
dx\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{\gamma_{5}}{4}\gamma_{\mu}S_{u}\left(\bar{q},x,0\right)\gamma_{5}S_{u}\left(\bar{q}-\bar{P},0,x\right)\gamma_{\mu}},$
(48) $\displaystyle K_{EF}^{\pi_{0}^{u}}$ $\displaystyle=-tr\int
dx\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{\gamma_{5}}{4}\gamma_{\mu}S_{u}\left(\bar{q},x,0\right)\frac{-i\gamma_{5}\not{\bar{P}}}{M_{u}}S_{u}\left(\bar{q}-\bar{P},0,x\right)\gamma_{\mu}},$
(49) $\displaystyle K_{FE}^{\pi_{0}^{u}}$ $\displaystyle=-tr\int
dx\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{-i\gamma_{5}\not{\bar{P}}M_{u}}{4\bar{P}^{2}}\gamma_{\mu}S_{u}\left(\bar{q},x,0\right)\gamma_{5}S_{u}\left(\bar{q}-\bar{P},0,x\right)\gamma_{\mu}},$
(50) $\displaystyle K_{FF}^{\pi_{0}^{u}}$ $\displaystyle=-tr\int
dx\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{-i\gamma_{5}\not{\bar{P}}M_{u}}{4\bar{P}^{2}}\gamma_{\mu}S_{u}\left(\bar{q},x,0\right)\frac{-i\gamma_{5}\not{\bar{P}}}{M_{u}}S_{u}\left(\bar{q}-\bar{P},0,x\right)\gamma_{\mu}}.$
(51)

For neutral pion, the integration over $q_{y}$ and $x$ can be done
analytically by the orthogonality of the normalized harmonic oscillator wave
function

$\displaystyle g_{1,n_{1},n_{2}}=\int
dx\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{1}}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{1}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(x+\frac{q_{y}}{e_{u}B}\right)=\frac{\left|e_{u}B\right|}{2\pi}\delta_{n_{1},n_{2}},$
(52) $\displaystyle g_{2,n_{1},n_{2}}=\int
dx\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{1}+1}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{1}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}-1}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(x+\frac{q_{y}}{e_{u}B}\right)=\frac{\left|e_{u}B\right|}{2\pi}\delta_{n_{1}+1,n_{2}},$
$\displaystyle g_{3,n_{1},n_{2}}=\int
dx\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{1}-1}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{1}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}+1}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(x+\frac{q_{y}}{e_{u}B}\right)=\frac{\left|e_{u}B\right|}{2\pi}\delta_{n_{1}-1,n_{2}},$
$\displaystyle g_{4,n_{1},n_{2}}=\int
dx\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{1}}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{1}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{2}}\left(x+\frac{q_{y}}{e_{u}B}\right)=\frac{\left|e_{u}B\right|}{2\pi}\delta_{n_{1},n_{2}},$

thus, these kernels can be written as

$\displaystyle K_{EE}^{\pi_{0}^{u}}$
$\displaystyle=2\sum^{\infty}_{n_{1}=0}\sum^{\infty}_{n_{2}=0}G_{1,n_{1},n_{2}}g_{1,n_{1},n_{2}}+G_{2,n_{1},n_{2}}g_{2,n_{1},n_{2}}+G_{3,n_{1},n_{2}}g_{3,n_{1},n_{2}}+G_{4,n_{1},n_{2}}g_{4,n_{1},n_{2}}$
(53)
$\displaystyle=\frac{\left|e_{u}B\right|}{2\pi}2\sum^{\infty}_{n_{u}=0}G_{1,n_{u},n_{u}}+G_{2,n_{u},n_{u}+1}+G_{3,n_{u},n_{u}-1}+G_{4,n_{u},n_{u}}$
$\displaystyle\equiv\frac{\left|e_{u}B\right|}{2\pi}2\sum^{\infty}_{n_{u}=0}G_{1,n_{u}}+G_{2,n_{u}}+G_{3,n_{u}}+G_{4,n_{u}},$
$\displaystyle K_{EF}^{\pi_{0}^{u}}$
$\displaystyle=\frac{\left|e_{u}B\right|}{2\pi}2\sum^{\infty}_{n=0}G_{5,n_{u}}+G_{6,n_{u}},$
$\displaystyle K_{FE}^{\pi_{0}^{u}}$
$\displaystyle=\frac{M_{u}^{2}}{2\bar{P}^{2}}K_{EF}^{\pi_{0}^{u}},$
$\displaystyle K_{FF}^{\pi_{0}^{u}}$ $\displaystyle=-2K_{FE}^{\pi_{0}^{u}}.$

In general, the form of $K_{FF}^{\pi_{+}}$ is deduced from WTIs. However, the
WTIs in magnetic field is very complicated, we use the subtraction
($K_{FF}^{\pi_{0}^{u}}=-2K_{FE}^{\pi_{0}^{u}}$), found in the zero magnetic
field, as an approximation Gutierrez-Guerrero  _et al._ (2010). $G_{j,n}$
correspond to the integration over $q_{z},q_{t}$ and should be regularized by
Eq. $(\ref{reg1})$, where

$\displaystyle
G_{1,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
k+M_{u}^{2}}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{u}}^{2})}$ (54) $\displaystyle
G_{2,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{2\sqrt{B_{n_{u}+1}B_{n_{u}+1}}}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{u}+1}^{2})}$
$\displaystyle
G_{3,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{2\sqrt{B_{n_{u}}B_{n_{u}}}}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{u}}^{2})}$
$\displaystyle
G_{4,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
k+M_{u}^{2}}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{u}+1}^{2})}$ $\displaystyle
G_{5,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot P-k\cdot
P}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{u}+1}^{2})}$ $\displaystyle
G_{6,n_{u}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot P-k\cdot
P}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{u}}^{2})}.$

###  B.2 charged pion

The applied kernels in Eq. (23)

$\displaystyle K_{EE^{\prime}}^{\pi_{+}}$
$\displaystyle=-tr\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{\gamma_{5}}{4}\gamma_{\mu}S_{u}\left(\bar{q},x,x^{\prime}\right)\gamma_{5}S_{d}\left(\bar{q}-\bar{P},x^{\prime},x\right)\gamma_{\mu}},$
(55) $\displaystyle K_{EF^{\prime}}^{\pi_{+}}$
$\displaystyle=-tr\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{\gamma_{5}}{4}\gamma_{\mu}S_{u}\left(\bar{q},x,x^{\prime}\right)\frac{-i\gamma_{5}\not{\bar{P}}}{2M}S_{d}\left(\bar{q}-\bar{P},x^{\prime},x\right)\gamma_{\mu}},$
(56) $\displaystyle K_{FE^{\prime}}^{\pi_{+}}$
$\displaystyle=-tr\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{-i\gamma_{5}\not{\bar{P}}2M}{4\bar{P}^{2}}\gamma_{\mu}S_{u}\left(\bar{q},x,x^{\prime}\right)\gamma_{5}S_{d}\left(\bar{q}-\bar{P},x^{\prime},x\right)\gamma_{\mu}},$
(57) $\displaystyle K_{FF^{\prime}}^{\pi_{+}}$
$\displaystyle=-tr\int\frac{d^{3}\bar{q}}{\left(2\pi\right)^{3}}{\frac{-i\gamma_{5}\not{\bar{P}}2M}{4\bar{P}^{2}}\gamma_{\mu}S_{u}\left(\bar{q},x,x^{\prime}\right)\frac{-i\gamma_{5}\not{\bar{P}}}{2M}S_{d}\left(\bar{q}-\bar{P},x^{\prime},x\right)\gamma_{\mu}},$
(58)

and can be rewritten as

$\displaystyle K_{EE^{\prime}}^{\pi_{+}}$
$\displaystyle=2\sum^{\infty}_{n_{u}=0}\sum^{\infty}_{n_{d}=0}G_{1,n_{u},n_{d}}g_{1,n_{u},n_{d}}+G_{2,n_{u},n_{d}}g_{2,n_{u},n_{d}}+G_{3,n_{u},n_{d}}g_{3,n_{u},n_{d}}+G_{4,n_{u},n_{d}}g_{4,n_{u},n_{d}},$
(59) $\displaystyle K_{EF^{\prime}}^{\pi_{+}}$
$\displaystyle=2\sum^{\infty}_{n_{u}=0}\sum^{\infty}_{n_{d}=0}G_{5,n_{u},n_{d}}g_{1,n_{u},n_{d}}+G_{6,n_{u},n_{d}}g_{4,n_{u},n_{d}},$
$\displaystyle K_{FE^{\prime}}^{\pi_{+}}$
$\displaystyle=\frac{2M^{2}}{\bar{P^{2}}}K_{EF^{\prime}}^{\pi_{+}},$
$\displaystyle K_{FF^{\prime}}^{\pi_{+}}$
$\displaystyle=-\frac{M_{u}+M_{d}}{2M}K_{FE^{\prime}}^{\pi_{+}}$

where the form of $K_{FF^{\prime}}^{\pi_{+}}$ is because of the subtraction
from WTIs. $g_{j,n_{u},n_{d}}$ correspond to the integration over $q_{y}$ and
carry the information of $x,x^{\prime}$, where

$\displaystyle
g_{1,n_{u},n_{d}}=\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{u}}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{u}}\left(x^{\prime}+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{d}}\left(x^{\prime}+\frac{q_{y}}{e_{d}B}\right)\hat{f}_{n_{d}}\left(x+\frac{q_{y}}{e_{d}B}\right),$
(60) $\displaystyle
g_{2,n_{u},n_{d}}=\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{u}+1}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{u}}\left(x^{\prime}+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{d}+1}\left(x^{\prime}+\frac{q_{y}}{e_{d}B}\right)\hat{f}_{n_{d}}\left(x+\frac{q_{y}}{e_{d}B}\right),$
$\displaystyle
g_{3,n_{u},n_{d}}=\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{u}-1}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{u}}\left(x^{\prime}+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{d}-1}\left(x^{\prime}+\frac{q_{y}}{e_{d}B}\right)\hat{f}_{n_{d}}\left(x+\frac{q_{y}}{e_{d}B}\right),$
$\displaystyle
g_{4,n_{u},n_{d}}=\int\frac{dq_{y}}{2\pi}\hat{f}_{n_{u}}\left(x+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{u}}\left(x^{\prime}+\frac{q_{y}}{e_{u}B}\right)\hat{f}_{n_{d}}\left(x^{\prime}+\frac{q_{y}}{e_{d}B}\right)\hat{f}_{n_{d}}\left(x+\frac{q_{y}}{e_{d}B}\right).$

$G_{j,n_{u},n_{d}}$ correspond to the integration over $q_{z},q_{t}$ and
should be regularized through Eq. $(\ref{reg1})$, where

$\displaystyle
G_{1,n_{u},n_{d}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
k+M_{u}M_{d}}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{d}+1}^{2})}$ (61)
$\displaystyle
G_{2,n_{u},n_{d}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{-2\sqrt{B_{n_{u}+1}B_{n_{d}+1}}}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{d}+1}^{2})}$
$\displaystyle
G_{3,n_{u},n_{d}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{-2\sqrt{B_{n_{u}}B_{n_{d}}}}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{d}}^{2})}$
$\displaystyle
G_{4,n_{u},n_{d}}=\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
k+M_{u}M_{d}}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{d}}^{2})}$ $\displaystyle
G_{5,n_{u},n_{d}}=\frac{1}{2M}\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
PM_{d}-k\cdot PM_{u}}{(q^{2}+M_{n_{u}+1}^{2})(k^{2}+M_{n_{d}}^{2})}$
$\displaystyle
G_{6,n_{u},n_{d}}=\frac{1}{2M}\int\frac{d^{3}q}{\left(2\pi\right)^{2}}\frac{q\cdot
PM_{d}-k\cdot PM_{u}}{(q^{2}+M_{n_{u}}^{2})(k^{2}+M_{n_{d}+1}^{2})}.$

We also give the detailed charged pion decay constants after including the
auxiliary term. Firstly, the normalization in Eq. (28) can be completed as

$\displaystyle\left[\frac{\partial\ln(\lambda(\bar{P}^{2}))}{\partial\bar{P}^{2}}\right]^{-1}$
$\displaystyle=$
$\displaystyle-2N_{c}\int_{x,x^{\prime},x^{\prime\prime}}\sum_{X=E,F}\left(K_{EX^{\prime}}^{\pi_{+}}X(x^{\prime})E(x-x^{\prime\prime})+\frac{1}{2M}\frac{\bar{P}^{2}}{M}K_{FX^{\prime}}^{\pi_{+}}X(x^{\prime})F(x-x^{\prime\prime})\right)$
(62) $\displaystyle=$
$\displaystyle-2N_{c}\int_{x,x^{\prime\prime}}\frac{3m_{G}^{2}}{4}E(x)E(x-x^{\prime\prime})+\frac{3m_{G}^{2}}{4}\frac{\bar{P}^{2}}{2(1-4\tilde{\eta}^{2})M^{2}}F(x)F(x-x^{\prime\prime})$
$\displaystyle=$
$\displaystyle-2N_{c}\frac{3m_{G}^{2}}{4}\frac{2\pi}{|eB|}(a_{E}^{2}+\frac{\bar{P}^{2}}{2(1-4\tilde{\eta}^{2})M^{2}}a_{F}^{2}).$

The ground states of the charged pion take the form of
$X(x)=a_{X}e^{-|eB|x^{2}/2}$ for $X=E,F$ to complete the integration of the
last second line. Then, the decay constants $f$ and $f^{\prime}$ are given by

$\displaystyle f_{\pi}$
$\displaystyle=\frac{N_{c}}{\bar{P}^{2}}\int_{x,x^{\prime}}\frac{\bar{P}^{2}}{M}\left[K_{FE^{\prime}}^{\pi_{+}}E(x^{\prime})+K_{FF^{\prime}}^{\pi_{+}}F(x^{\prime})\right]$
(63)
$\displaystyle=\frac{3m_{G}^{2}}{4}\frac{N_{c}}{(1-4\tilde{\eta}^{2})M}\sqrt{\frac{2\pi}{|eB|}}a_{F},$
$f^{\prime}_{\pi}=\frac{N_{c}}{\left|eB\right|\bar{P}^{2}}\int_{x,x^{\prime}}\frac{\bar{P}^{2}}{M}\left[\tilde{K}_{FE^{\prime}}^{\pi_{+}}E(x^{\prime})+\tilde{K}_{FF^{\prime}}^{\pi_{+}}F(x^{\prime})\right],$
(64)

where

$\displaystyle\tilde{K}_{FE^{\prime}}^{\pi_{+}}$
$\displaystyle=\frac{4M^{2}}{\bar{P^{2}}}\sum^{\infty}_{n_{u}=0}\sum^{\infty}_{n_{d}=0}G_{5,n_{u},n_{d}}g_{1,n_{u},n_{d}}-G_{6,n_{u},n_{d}}g_{4,n_{u},n_{d}},$
(65) $\displaystyle\tilde{K}_{FF^{\prime}}^{\pi_{+}}$
$\displaystyle=-\frac{M_{u}+M_{d}}{2M}\tilde{K}_{FE^{\prime}}^{\pi_{+}}.$ (66)

By analogy with $K_{FF^{\prime}}^{\pi_{+}}$, the form of
$\tilde{K}_{FF^{\prime}}^{\pi_{+}}$ also comes from the subtraction of WTIs.

## References

  * Roberts  _et al._ (2021) C. D. Roberts, D. G. Richards, T. Horn, and L. Chang, Prog. Part. Nucl. Phys. 120, 103883 (2021), arXiv:2102.01765 [hep-ph] .
  * Klevansky (1992) S. P. Klevansky, Rev. Mod. Phys. 64, 649 (1992).
  * Buballa (2005) M. Buballa, Phys. Rept. 407, 205 (2005), arXiv:hep-ph/0402234 .
  * Brandt  _et al._ (2015) B. B. Brandt, A. Francis, H. B. Meyer, and D. Robaina, Phys. Rev. D 92, 094510 (2015), arXiv:1506.05732 [hep-lat] .
  * Gao and Ding (2020) F. Gao and M. Ding, Eur. Phys. J. C 80, 1171 (2020), arXiv:2006.05909 [hep-ph] .
  * Agasian and Fedorov (2008) N. Agasian and S. Fedorov, Physics Letters B 663, 445 (2008).
  * D’Elia  _et al._ (2010) M. D’Elia, S. Mukherjee, and F. Sanfilippo, Phys. Rev. D 82, 051501 (2010).
  * Mizher  _et al._ (2010) A. J. Mizher, M. N. Chernodub, and E. S. Fraga, Phys. Rev. D 82, 105016 (2010).
  * Bali  _et al._ (2012) G. S. Bali, F. Bruckmann, G. Endrődi, Z. Fodor, S. D. Katz, S. Krieg, A. Schäfer, and K. K. Szabó, JHEP 02, 044 (2012), arXiv:1111.4956 [hep-lat] .
  * Andersen  _et al._ (2016) J. O. Andersen, W. R. Naylor, and A. Tranberg, Rev. Mod. Phys. 88, 025001 (2016), arXiv:1411.7176 [hep-ph] .
  * Thompson and Beloborodov (2005) C. Thompson and A. M. Beloborodov, The Astrophysical Journal 634, 565 (2005).
  * Rea (2013) N. Rea, in  _Neutron Stars and Pulsars: Challenges and Opportunities after 80 years_, Vol. 291, edited by J. van Leeuwen (2013) pp. 11–18, arXiv:1211.2086 [astro-ph.HE] .
  * Kaspi and Beloborodov (2017) V. M. Kaspi and A. Beloborodov, Annu. Rev. Astron. Astrophys. 55, 261 (2017), arXiv: 1703.00068.
  * Skokov  _et al._ (2009) V. Skokov, A. Illarionov, and V. Toneev, Int. J. Mod. Phys. A 24, 5925 (2009), arXiv:0907.1396 [nucl-th] .
  * Voronyuk  _et al._ (2011) V. Voronyuk, V. Toneev, W. Cassing, E. Bratkovskaya, V. Konchakovski, and S. Voloshin, Phys. Rev. C 83, 054911 (2011), arXiv:1103.4239 [nucl-th] .
  * Deng and Huang (2012) W.-T. Deng and X.-G. Huang, Phys. Rev. C 85, 044907 (2012).
  * Fayazbakhsh and Sadooghi (2013) S. Fayazbakhsh and N. Sadooghi, Phys. Rev. D 88, 065030 (2013), arXiv:1306.2098 [hep-ph] .
  * Avancini  _et al._ (2016) S. S. Avancini, W. R. Tavares, and M. B. Pinto, Phys. Rev. D 93, 014010 (2016), arXiv:1511.06261 [hep-ph] .
  * Zhang  _et al._ (2016) R. Zhang, W.-j. Fu, and Y.-x. Liu, Eur. Phys. J. C 76, 307 (2016), arXiv:1604.08888 [hep-ph] .
  * Wang and Zhuang (2018) Z. Wang and P. Zhuang, Phys. Rev. D 97, 034026 (2018), arXiv:1712.00554 [hep-ph] .
  * Liu  _et al._ (2018) H. Liu, X. Wang, L. Yu, and M. Huang, Phys. Rev. D 97, 076008 (2018), arXiv:1801.02174 [hep-ph] .
  * Mao (2019) S. Mao, Phys. Rev. D 99, 056005 (2019), arXiv:1808.10242 [nucl-th] .
  * Gómez Dumm  _et al._ (2018) D. Gómez Dumm, M. Izzo Villafañe, and N. Scoccola, Phys. Rev. D 97, 034025 (2018), arXiv:1710.08950 [hep-ph] .
  * Ayala  _et al._ (2018) A. Ayala, R. L. Farias, S. Hernández-Ortiz, L. Hernández, D. M. Paret, and R. Zamora, Phys. Rev. D 98, 114008 (2018), arXiv:1809.08312 [hep-ph] .
  * Coppola  _et al._ (2019) M. Coppola, D. Gomez Dumm, S. Noguera, and N. Scoccola, Phys. Rev. D 100, 054014 (2019), arXiv:1907.05840 [hep-ph] .
  * Sheng  _et al._ (2021) B. Sheng, Y. Wang, X. Wang, and L. Yu, Phys. Rev. D 103, 094001 (2021), arXiv:2010.05716 [hep-ph] .
  * Das and Haque (2020) A. Das and N. Haque, Phys. Rev. D 101, 074033 (2020), arXiv:1908.10323 [hep-ph] .
  * Luschevskaya  _et al._ (2016) E. V. Luschevskaya, O. E. Solovjeva, and O. V. Teryaev, Phys. Lett. B 761, 393 (2016), arXiv:1511.09316 [hep-lat] .
  * Shushpanov and Smilga (1997) I. A. Shushpanov and A. V. Smilga, Phys. Lett. B 402, 351 (1997), arXiv:hep-ph/9703201 .
  * Andersen (2012) J. O. Andersen, Phys. Rev. D 86, 025020 (2012).
  * Colucci  _et al._ (2014) G. Colucci, E. S. Fraga, and A. Sedrakian, Phys. Lett. B 728, 19 (2014), arXiv:1310.3742 [nucl-th] .
  * Chao  _et al._ (2020) J. Chao, M. Huang, and A. Radzhabov, Chin. Phys. C 44, 034105 (2020), arXiv:1805.00614 [hep-ph] .
  * Bali  _et al._ (2018a) G. S. Bali, B. B. Brandt, G. Endrődi, and B. Gläßle, Phys. Rev. D 97, 034505 (2018a), arXiv:1707.05600 [hep-lat] .
  * Ding  _et al._ (2020) H.-T. Ding, S.-T. Li, S. Mukherjee, A. Tomiya, and X.-D. Wang, PoS LATTICE2019, 250 (2020), arXiv:2001.05322 [hep-lat] .
  * Xu  _et al._ (2021) K. Xu, J. Chao, and M. Huang, Phys. Rev. D 103, 076015 (2021), arXiv:2007.13122 [hep-ph] .
  * Hilger  _et al._ (2017) T. Hilger, M. Gomez-Rocha, and A. Krassnigg, Eur. Phys. J. C 77, 625 (2017), arXiv:1508.07183 [hep-ph] .
  * Bhagwat  _et al._ (2007) M. S. Bhagwat, A. Krassnigg, P. Maris, and C. D. Roberts, Eur. Phys. J. A 31, 630 (2007), arXiv:nucl-th/0612027 .
  * Chernodub (2010) M. N. Chernodub, Phys. Rev. D 82, 085011 (2010), arXiv:1008.1055 [hep-ph] .
  * Bali  _et al._ (2018b) G. S. Bali, B. B. Brandt, G. Endrődi, and B. Gläßle, Phys. Rev. Lett. 121, 072001 (2018b).
  * Gutierrez-Guerrero  _et al._ (2010) L. Gutierrez-Guerrero, A. Bashir, I. Cloet, and C. Roberts, Phys. Rev. C 81, 065202 (2010), arXiv:1002.1968 [nucl-th] .
  * Roberts  _et al._ (2011) H. L. Roberts, L. Chang, I. C. Cloet, and C. D. Roberts, Few Body Syst. 51, 1 (2011), arXiv:1101.4244 [nucl-th] .
  * Kojo and Su (2013) T. Kojo and N. Su, Phys. Lett. B 726, 839 (2013), arXiv:1305.4510 [hep-ph] .
  * Chyi  _et al._ (2000) T.-K. Chyi, C.-W. Hwang, W. F. Kao, G.-L. Lin, K.-W. Ng, and J.-J. Tseng, Phys. Rev. D 62, 105014 (2000), arXiv:hep-th/9912134 .
  * Wang (2013) K. Wang,  _QCD phase transition and properties of phases with Dyson-Schwinger equations_ , Ph.D. thesis, Peking University (2013), unpublished thesis.
  * Chao  _et al._ (2014) J. Chao, L. Yu, and M. Huang, Phys. Rev. D 90, 045033 (2014), [Erratum: Phys.Rev.D 91, 029903 (2015)], arXiv:1403.0442 [hep-th] .