---
abstract: 'The formation of charged pion condensate under parallel electromagnetic fields is studied within the two-flavor Nambu–Jona-Lasinio model. The technique of Schwinger proper time method is extended to explore the quantity locating in the off-diagonal flavor space, i.e., charged pion. We obtain the associated effective potential as a function of the strength of the electromagnetic fields and find out that it contains a sextic term which possibly induce weakly first order phase transition. Dependence of pion condensation on model parameters is investigated.'
address:
- 'Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, 730000, China'
- 'Institute of High Energy Physics, Chinese Academy of Sciences, Beijing, 100049, China'
- 'School of Physics Sciences, University of Chinese Academy of Sciences, Beijing 100039, China'
- 'Matrosov Institute for System Dynamics and Control Theory , 664033, Irkutsk, Russia'
- 
author:
- Jingyi Chao
- Mei Huang
- Andrey Radzhabov
title: Charged pion condensation under parallel electromagnetic fields
---

Introduction
============

The phase structure of Quantum Chromodynamics (QCD) at high temperature/density and other extreme conditions has attracted lots of attentions and been a main topic of heavy ion collisions. The perturbative QCD predicts a free gas of quarks and gluons at high temperature limit and a color-flavor-locking phase at very high baryon density but low temperature. However, the QCD vacuum has a rather complicated nonperturbative structure, and the QCD phase diagram is not a simple transition between the hadron phase with non-zero chiral condensate to the weakly coupled quark-gluon plasma as expected long time ago [@Cleymans:1985wb], but instead a rich structure of different phases with corresponding condensates. These phases could include different color superconducting states or inhomogeneous chiral condensates [@Buballa:2003qv; @Buballa:2014tba; @Andersen:2018osr]. Recently, QCD phase structure under strong magnetic fields has drawn great interests  [@Kharzeev:2007jp; @Skokov:2009qp; @Hattori:2016emy; @Andersen:2014xxa; @Miransky:2015ava; @Huang:2015oca]. The strong magnetic fields can be generated with the strength up to $B\sim10^{18\sim 20}$ G in the non-central heavy ion collisions [@Skokov:2009qp; @Deng:2012pc], and is expected to be on the order of $10^{18}$-$10^{20}$ G [@Duncan:1992hi; @Blaschke:2018mqw] in the inner core of magnetars.

Lots of interesting phenomena under strong magnetic fields have been discussed, for example, the magnetic catalysis [@Klevansky:1989vi; @Klimenko:1990rh; @Gusynin:1995nb; @Gusynin:1999pq], inverse magnetic catalysis [@Bali:20111213; @Bali:2012zg; @Bali:2013esa] effect, the chiral magnetic effect (CME) [@Kharzeev:2007jp; @Kharzeev:2007tn; @Fukushima:2008xe] and the vacuum superconductivity [@Chernodub:2010qx; @Chernodub:2011mc]. Moreover, it was pointed out that under the parallel electromagnetic fields, the neutral pion condensation can be formed [@Cao:2015cka; @Wang:2018gmj] due to the connection of field with axial anomaly. If only QCD interaction is included, the axial isospin currents is anomaly free. It turns out that anomaly emerges associated with the coupling of quarks to electromagnetism, where the axial isospin currents is given by $${\partial}_{\mu}j_{5}^{\,\mu 3}=-\frac{e^{2}}{16\pi^{2}}{\varepsilon}^{{\alpha}{\beta}\mu\nu}F_{{\alpha}{\beta}}F_{\mu\nu}\cdot{\hbox{tr}}{\left}[\tau^{3}Q^{2}{\right}].$$ Here $Q$ is the matrix of quark electric charges and $F$ is the field strength. The corresponding process is $\pi_{0}\to {\gamma}{\gamma}$. The decay of a neutral pion into two photons, which had been a puzzle for some time in the 1960s, is the most successful proof of chiral anomaly. Above solution led to the discovery of the Adler–Bell–Jackiw anomaly [@Adler:2004qt].

In the asymmetric flavor space, one can introduce a chiral isospin chemical potential $\mu_{I}^5$ corresponding to the current $\bar{\psi}{\gamma}_{0}{\gamma}_{5}\tau_{3}\psi$, which is similar to the isospin chemical potential $\mu_{I}$ with respect to $\bar{\psi}{\gamma}_{0}\tau_{3}\psi$. It has been a long history of investigating the pion condensation under the isospin asymmetric nuclear matter. In the beginning this effect is discussed for case nuclear matter in neutron-star interiors [@Sawyer:1972cq; @Sawyer:1973fv; @Voskresensky:1980nk] or superdense and supercharged nuclei [@Migdal:1978az]. The pion condensation of charged or neutral pion modes in QCD vacuum are also considered in the frameworks of effective models with quark degrees of freedom [@Son:2000xc; @He:2006tn; @Mao:2014hga; @Khunjua:2017khh] or in lattice calculations [@Brandt:2016zdy; @Brandt:2018bwq].

The degeneracy between $\pi_{0}$ and $\pi_{\pm}$ is destroyed because of the axial isospin chemical potential. It is worth to pursuing the detailed behaviors of charged pions in a strict manner. Hence, in this work, we focus on the possibility of charged pion condensation under the parallel electromagnetic fields in the framework of the ${\mathrm}{SU}(2) \times {\mathrm}{SU}(2)$ NJL model [@Nambu:1961tp; @Nambu:1961fr]. For this purpose, we develop a full routine to derive the mean-field thermodynamical potential of the NJL model with nonzero charged pion condensate ${\langle}\bar{\psi} i \gamma_5 \tau_{\pm} \psi {\rangle}$ in the off-diagonal flavor space under the parallel electromagnetic fields. Calculations are performed with Schwinger proper time method [@Schwinger:1951nm] and the proper time regularization in the NJL model is used. Through the paper we only consider the model at zero temperature and chemical potential and restrict ourselves to the case of the electric field anti-parallel to the magnetic field.

\[intr\]

Lagrangian {#NJL}
==========

The Lagrangian of the ${\mathrm}{SU}(2)\times {\mathrm}{SU}(2)$ NJL model is in the form of [@Nambu:1961tp; @Nambu:1961fr; @Volkov:1986zb; @Vogl:1991qt; @Klevansky:1992qe; @Hatsuda:1994pi; @Volkov:2005kw] $$\begin{aligned}
\mathcal{L}_{{\mathrm}{NJL}} = \bar{\psi}\left(i \slashed{D} - m_0 \right) \psi
+ \mathrm{G} \left[ \left(  \bar{\psi} \psi \right)^2+ \left(  \bar{\psi} i\gamma_5 \tau_i \psi \right)^2 \right],\end{aligned}$$ where $\bar{\psi}(x)=(\bar{u}(x),\bar{d}(x))$ are $u$ and $d$ anti-quark fields. The limit of equal current masses for $u,d$, $m_u=m_d\equiv m_{0}$ is considered. $\gamma_i$, $\tau_{i}$ are conventional Dirac and Pauli matrices and $\tau_{0}$ is the unit matrix. $\slashed{D}$ is the covariant derivative and, in the two flavor space, expressed as $$\begin{aligned}
    D_{\mu}={\left}({\partial}_{\mu}-{\mathop{}\!i}QA_{\mu}{\right})\tau_{0}-{\mathop{}\!i}qA_{\mu}\tau_{3},\end{aligned}$$ where $Q=\frac{1}{2}{\left}(q_{u}+q_{d}{\right})$ and $q=\frac{1}{2}{\left}(q_{u}-q_{d}{\right})$.

Introducing auxiliary bosonic fields $\pi$, $\sigma$, with the help of Hubbard-Stratonovich transformation one can integrate over the quark fields, then obtains the following effective Lagrangian: $$\begin{aligned}
\mathcal{L} = \frac{\sigma^2+\vec{\pi}^2}{4\mathrm{G}}-{\mathop{}\!i}{\hbox{Tr}}\ln S^{-1}, \label{eLagrangian}\end{aligned}$$ where $S^{-1}$ is the inverse quark propagator and $$\begin{aligned}
S^{-1}= {\mathop{}\!i}\slashed{D} - M, \quad M=m_{0}\tau_{0}-\sigma\tau_{0}-{\mathop{}\!i}\gamma_5\pi_{i} \tau_i. \label{Propagator}\end{aligned}$$ The auxiliary bosonic fields could have a nonzero vacuum expectation values and therefore it is necessary to shift them as $\sigma = \sigma^\prime-{\langle}\sigma{\rangle}$, $\pi_{i}=\pi_{i}^\prime- {\langle}\pi_{i}{\rangle}$. Equations of motion for mean-fields ${\langle}\sigma{\rangle}$, ${\langle}\pi_{i}{\rangle}$ are obtained from the Lagrangian (\[eLagrangian\]) after elimination from its linear terms, i.e. $$\begin{aligned}
\frac{\delta \mathcal{L} }{\delta {\langle}\sigma{\rangle}}\biggl|_{\substack{ \sigma^\prime =0 \\ \pi_{i}^\prime=0 }}, \quad
\frac{\delta \mathcal{L} }{\delta {\langle}\pi_{i}{\rangle}}\biggl|_{\substack{ \sigma^\prime =0 \\ \pi_{i}^\prime=0 }} =0.\end{aligned}$$ As a result, under different conditions the ${\langle}{\sigma}{\rangle}$, ${\langle}\pi_{i}{\rangle}$ condensates have non-zero values and the non-zero value of scalar condensate leads to a formation of constituent quarks with dynamical quark mass $m=m_{0}-{\langle}{\sigma}{\rangle}$.

Let us denote the second term of effective Lagrangian (\[eLagrangian\]) as $\mathcal{S}_{eff}=-{\mathop{}\!i}{\hbox{Tr}}\ln S^{-1}$. Then the gap equations for ${\langle}{\sigma}{\rangle}$ and ${\langle}\pi_{i}{\rangle}$ takes the form $$\begin{aligned}
m=m_{0}-2\mathrm{G}\frac{\partial\mathcal{S}_{eff}}{\partial {\langle}\sigma{\rangle}},\quad
{\langle}\pi_{i}{\rangle}=-2\mathrm{G}\frac{\partial \mathcal{S}_{eff}}{\partial {\langle}\pi_{i}{\rangle}}.\end{aligned}$$

The calculation of $\mathcal{S}_{eff}$ is presented in the following section.

The effective potential {#main}
=======================

Without loss of generality, one can choose ${\langle}\pi_{i}{\rangle}={\left}(\pi_{1},0,0{\right})$ and therefore “mass” in quark propagator Eq.(\[Propagator\]) is $M=m\tau_{0}+{\mathop{}\!i}\pi_{1}{\gamma}_{5}\tau_{1}$. Since ${\mathrm}{Det}{\left}({\mathop{}\!i}\slashed{D}-M{\right})={\mathrm}{Det}\,{\Gamma}{\left}({\mathop{}\!i}\slashed{D}-M{\right}){\Gamma}$, where ${\Gamma}={\gamma}_{5}\tau_{3}$, the second term of the Lagrangian Eq.(\[eLagrangian\]) is replaced to $$\mathcal{S}_{eff}=-\frac{{\mathop{}\!i}}{2}\ln{\mathrm}{Det}{\left}(\slashed{\mathcal{D}}^{2}+m^{2}+\pi_{1}^{2}{\right}),$$ where $\slashed{\mathcal{D}}^{2}=\slashed{D}^{2}-{\gamma}_{5}{\gamma}^{\mu}\pi_{1}{\left}[\tau_{1},D_{\mu}{\right}]$.

By using the method of proper time, we represent $\mathcal{S}_{eff}$ as following: $$\label{eqn_seff}
    \mathcal{S}_{eff}={\hbox{Tr}}\int\limits_{1/{\Lambda}^{2}}^{\infty}{\mathop{}\!i}\,\frac{{\mathop{}\!d}s}{2s}\int{\hbox{tr}}{\left}{\langle}x\big|{\mathop{}\!e}^{-{\mathop{}\!i}{\left}(\slashed{\mathcal{D}}^{2}+m^{2}+\pi_{1}^{2}{\right})s}\big|x'{\right}{\rangle}{\mathop{}\!d}^{4}x,$$ where the ultraviolet cutoff $1/{\Lambda}^{2}$ has been explicitly introduced, ${\hbox{tr}}$ and ${\hbox{Tr}}$ means the trace taking in the spinor and flavor space, respectively.

From now on, we will work in the Euclidean space. Following notations are introduced: $$\begin{aligned}
\label{AlBeLa}
&   {\alpha}=m^{2}+\pi_{1}^{2}-\frac{1}{2}{\sigma}^{\mu\nu}{\lambda}_{\mu\nu},\quad
    {\beta}_{\nu}=q\pi_{1}{\gamma}_{5}{\gamma}^{\mu}F_{\mu\nu}\tau_{2},\quad\nonumber \\
&    {\lambda}_{\mu\nu}=q_{f}F_{\mu\nu},\end{aligned}$$ where $q_{f}={\mathrm}{Diag}(q_{u},q_{d})$ and ${\sigma}^{\mu\nu}=\frac{{\mathop{}\!i}}{2}{\left}[{\gamma}^{\mu},{\gamma}^{\nu}{\right}]$. In order to obtain $\mathcal{S}_{eff}$, it is then straightforward to look for the solution of $G(x,y;s)$ obeying a second order differential equation ${\left}(\slashed{\mathcal{D}}^{2}+m^{2}+\pi_{1}^{2}{\right})G{\left}(x,y;s{\right})={\delta}{\left}(x,y;s{\right})$. The explicit form is $$\begin{aligned}
\slashed{\mathcal{D}}^{2}+m^{2}+\pi_{1}^{2}&={\partial}^{2}_{x}+\alpha(y)+\beta_{\mu}(y){\left}(x-y{\right})^{\mu}+\nonumber\\
&   +\frac{1}{4}{\lambda}^{2}_{\mu\nu}{\left}(x-y{\right})^{\mu}{\left}(x-y{\right})^{\nu}.\end{aligned}$$ Performing the Fourier transform, one finds, $$\begin{aligned}
\label{eqn_diff_p}
    {\left}(-p^{2}+\alpha-{\mathop{}\!i}\beta_{\mu}\frac{{\partial}}{{\partial}p_{\mu}}-\frac{1}{4}{\lambda}^{2}_{\mu\nu}\frac{{\partial}^{2}}{{\partial}p_{\mu}{\partial}p_{\nu}}{\right})G(p;s)=1.\end{aligned}$$ As suggested in the the reference [@Brown:1975bc] one can solve the equation in the form $$\begin{aligned}
\label{eqn_G_pV2}
    G(p;s)={\mathop{}\!e}^{-\alpha s}{\mathop{}\!e}^{ p\cdot A(s)\cdot p+B(s) \cdot p+C(s)},\end{aligned}$$ whose associated descriptions of matrix $A$, vector $B$ and scalar $C$ are $$\begin{aligned}
\label{eqn_sol_ABC}
&A={\lambda}^{-1}\tan{\lambda}s,\quad
B=-2{\mathop{}\!i}{\beta}\cdot{\lambda}^{-2}{\left}(1-\sec{\lambda}s{\right}),\\
&C=-\frac{1}{2}\,{\hbox{tr}}\ln\cos{\lambda}s-{\beta}\cdot{\lambda}^{-3}{\left}(\tan{\lambda}s-{\lambda}s{\right})\cdot{\beta}\nonumber.\end{aligned}$$ For simplicity here and below indexes are not shown.

Plugging the form of ${\beta}$ in Eq.(\[AlBeLa\]) into vector $B$ and restoring indexes one has $$\begin{aligned}
    B_{\mu}=-2{\mathop{}\!i}q\pi_{1}\tau_{2}{\gamma}_{5}{\gamma}^{\nu}F_{\nu{\alpha}}{\left}[{\lambda}^{-2}{\left}(1-\sec{\lambda}s{\right}){\right}]^{{\alpha}}_{\mu}.\end{aligned}$$ Vector $B$ contains Dirac matrix, not commuting with ${\sigma}^{\mu\nu}$. Therefore, we emphasize that one should be careful with tracing in spinor space and integrating in momentum space. Introducing notations $P_{1}=\frac{1}{2}{\sigma}{\lambda}s$ and $P_{2}=p\cdot A(s)\cdot p+B(s) \cdot p$, one has ${\left}[{\sigma}{\lambda}s, p\cdot A(s)\cdot p{\right}]=0$ and the part with matrices in exponent Eq.(\[eqn\_G\_pV2\]) can be expanded as $$\begin{aligned}
&{\mathop{}\!e}^{P_{1}+P_{2}}\simeq {\mathop{}\!e}^{P_{1}}{\mathop{}\!e}^{P_{2}}{\mathop{}\!e}^{-\frac{1}{2}[P_{1},P_{2}]}=\nonumber \\
&\quad\quad={\mathop{}\!e}^{\frac{1}{2}{\sigma}{\lambda}s}{\mathop{}\!e}^{p\cdot A(s)\cdot p+B(s) \cdot p}{\mathop{}\!e}^{-\frac{1}{4}{\left}[{\sigma}{\lambda}s,B(s) \cdot p{\right}]}.\end{aligned}$$ We denote $-\frac{1}{4}{\left}[{\sigma}{\lambda}s,B(s) \cdot p{\right}]=\frac{1}{2}q\pi_{1} Os$, where $O$ has a structure of the form $O=Q\tau_{2}O_{1}\mathbb{B}_{1}p+q\tau_{1}O_{2}\mathbb{B}_{2}p$ and $\mathbb{B}$ will render in Eq. (\[Bdefinition\]). Shorthand matrix notation is applied, i.e. $\mathbb{F}=F_{\mu}^{\nu}$. To find the eigenvalue of $O$, we square it and get $$\begin{aligned}
O^{2}&=Q^{2}{\left}(\tau_{2}O_{1}\mathbb{B}_{1}p{\right})^{2}+q^{2}{\left}(\tau_{1}O_{2}\mathbb{B}_{2}p{\right})^{2} -\nonumber \\
&-{\mathop{}\!i}qQ\tau_{3}{\left}[O_{1}\mathbb{B}_{1}p,O_{2}\tilde{\mathbb{B}}_{2}p{\right}],\\
&O_{1}={\mathop{}\!i}{\left}[{\sigma}_{\mu\nu},{\gamma}_{5}{\gamma}^{{\alpha}} {\right}]=2{\gamma}_{5}g_{\nu}^{{\alpha}}{\gamma}_{\mu}-2{\gamma}_{5}g_{\mu}^{{\alpha}}{\gamma}_{\nu},\nonumber\\
&O_{2}={\left}\{{\sigma}_{\mu\nu},{\gamma}_{5}{\gamma}^{{\alpha}}{\right}\}=-2{\varepsilon}^{{\alpha}{\beta}}_{\;\mu\nu}{\gamma}_{{\beta}}.\nonumber\end{aligned}$$ With help of relation $\tau_{2}q_{f}\tau_{2}={\mathrm}{Diag}{\left}(q_{d}, q_{u}{\right})=\tilde{q}_{f}$, the $\tilde{\mathbb{B}}{\left}(\mathbb{B}{\right})$ are shown as $$\begin{aligned}
\label{Bdefinition}
&\tilde{\mathbb{B}}_{1}{\left}(\mathbb{B}_{1}{\right})=\frac{1}{\mathsf{q}^{2}}{\left}[1-\sec \mathsf{q}\mathbb{F}s{\right}],\quad \nonumber \\
&\tilde{\mathbb{B}}_{2}{\left}(\mathbb{B}_{2}{\right})=
\frac{\bar{\mathbb{F}}\mathbb{F}}{\mathbb{F}^{2}}
\frac{1}{\mathsf{q}^{2}}
{\left}[1-\sec\mathsf{q}\mathbb{F}s{\right}],\end{aligned}$$ where $\mathsf{q}=\tilde{q}_{f}$ or $q_{f}$ for $\tilde{\mathbb{B}},\mathbb{B}$ respectively; $\mathbb{F}$ and $\bar{\mathbb{F}}$ are field strength tensor $F^{\mu\nu}$ and dual field strength tensor $\bar{F}^{\mu\nu}=\frac{1}{2}{\varepsilon}^{\mu\nu{\alpha}{\beta}}F_{{\alpha}{\beta}}$, in shorthand notations. Moreover, ${\left}(\tau_{2}O_{1}\mathbb{B}_{1}p{\right})^{2}=-16\mathbb{B}_{1}\tilde{\mathbb{B}}_{1}p^{2}$, ${\left}(\tau_{1}O_{2}\mathbb{B}_{2}p{\right})^{2}=16\mathbb{B}_{2}\tilde{\mathbb{B}}_{2}p^{2}$ and $[O_{1}\mathbb{B}_{1}p,O_{2}\tilde{\mathbb{B}}_{2}p]=-32{\gamma}_{5}\mathbb{B}_{1}\tilde{\mathbb{B}}_{2}p^{2}$.

Applying the system that in a Lorentz frame where the electromagnetic field vectors are anti-parallel, e.g., ${\mathbf}{B}=-{\mathbf}{E}=f\hat{z}$, one gets $\mathbb{F}^{2}=f^2\,{\mathrm}{Diag}{\left}(-,+,+,-{\right})$ and $\bar{\mathbb{F}}\mathbb{F}=-f^2{\delta}_{\mu\nu}$ in Euclidean metric $(-,-,-,-)$, hence that $\bar{\mathbb{F}}\mathbb{F}/\,\mathbb{F}^{2}=f^{2}\mathbb{F}^{-2}$. Besides, ${\left}[1-\sec\mathsf{q}\mathbb{F}s{\right}]$ contains even powers of $\mathbb{F}$. It causes $O^{2}=-16Q^{2}\mathsf{p}_{1}^{2}+16q^{2}\mathsf{p}_{2}^{2}+32{\mathop{}\!i}{\gamma}_{5}Qq\mathsf{p}_{1}\cdot\mathsf{p}_{2}$ in a simply manner, where $\mathsf{p}_{1}=p_{\shortparallel}+p_{\perp}$, $\mathsf{p}_{2}=p_{\shortparallel}-p_{\perp}$, $p_{\shortparallel}=b_{\shortparallel}(p_{0},0,0,p_{3})$ and $p_{\perp}=b_{\perp}(0,p_{1},p_{2},0)$. The forms of $b_{\shortparallel}$ and $b_{\perp}$ are taken as $$\label{eqn_b03}
    b_{\shortparallel}=\frac{(1-\sec q_{f}s)^{\frac{1}{2}}(1-\sec\tilde{q}_{f}s)^{\frac{1}{2}}}{q_{f}\tilde{q}_{f}},$$ $$\label{eqn_b12}
    b_{\perp}=\frac{(1-\operatorname{sech}q_{f}s)^{\frac{1}{2}}(1-\operatorname{sech}\tilde{q}_{f}s)^{\frac{1}{2}}}{q_{f}\tilde{q}_{f}}.$$ Here and below in we rescale the integration variable as $s=s^\prime/f$ and omit prime. Because ${\gamma}_{5}^{2}=1$ associated with eigenvalue $\pm 1$, it follows that $O$ has four eigenvalues [@Schwinger:1951nm], written as $$\begin{aligned}
\label{eqn_O_squar}
    \mathcal{O}=\pm 4{\left}({\mathop{}\!i}Q\mathsf{p}_{1}\pm q\tau_{3}\mathsf{p}_{2}{\right}).\end{aligned}$$

Let ${\theta}=q\pi_{1} s/f$, one has $$\begin{aligned}
    {\hbox{tr}}\,{\mathop{}\!e}^{\frac{1}{2}{\theta}O}=\mathsf{T}=\cos{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\cosh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right}),\end{aligned}$$ which follow the method applied in [@dittrich2000probing]. The full statement is that $$\begin{aligned}
\label{eqn_O_exp}
    \exp{\left}[\frac{1}{2}{\theta}O{\right}]=\mathsf{T}+{\mathop{}\!i}{\gamma}_{5}\mathsf{U}+\frac{O\mathsf{V}}{2K^{2}}+\frac{{\mathop{}\!i}{\gamma}_{5}O\mathsf{W}}{2K^{2}},\end{aligned}$$ where $K^{2}=\mathsf{p}_{1}^{2}=\mathsf{p}_{2}^{2}$. $\mathsf{T}, \mathsf{U}, \mathsf{V}$ and $\mathsf{W}$ are scalars. Similarly, $$\begin{aligned}
\label{eqn_F_exp}
    \exp{\left}[q_{f}\frac{{\sigma}Fs}{2f}{\right}]=\mathsf{P}-{\mathop{}\!i}{\gamma}_{5}\mathsf{Q}+\frac{{\sigma}F}{2f}\,\mathsf{R}-\frac{{\mathop{}\!i}{\gamma}_{5}{\sigma}F}{2f}\,\mathsf{S}.\end{aligned}$$ Since $$\begin{aligned}
    {\hbox{tr}}{\left}(O^{2}{\mathop{}\!e}^{\frac{1}{2}{\theta}O}{\right})=\frac{{\partial}^{2}}{{\partial}^{2}{\theta}}{\hbox{tr}}{\left}(4{\mathop{}\!e}^{\frac{1}{2}{\theta}O}{\right})=4\frac{{\partial}^{2}\mathsf{T}}{{\partial}^{2}{\theta}},\end{aligned}$$ apply the identity of [Eq. (\[eqn\_O\_squar\])]{}, it derives that $$\begin{aligned}
    \mathsf{U}=\sin{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\sinh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right}).\end{aligned}$$ Proceeding with the direct differentiation of the exponential function via our basic trick, we get $$\begin{aligned}
&\mathsf{V}=\frac{1}{Q^{2}+q^{2}}
\biggl(Q\mathsf{p}_{1}\sin{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\cosh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right})+\nonumber\\
&\quad\quad\quad+\tau_{3}q\mathsf{p}_{2}\cos{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\sinh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right})\biggr), \nonumber\\
&\mathsf{W}=\frac{1}{Q^{2}+q^{2}}
\biggl(\tau_{3}q\mathsf{p}_{2}\sin{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\cosh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right})-\nonumber\\
&\quad\quad\quad-Q\mathsf{p}_{1}\cos{\left}(2Q{\theta}\mathsf{p}_{1}{\right})\sinh{\left}(2\tau_{3}q{\theta}\mathsf{p}_{2}{\right})\biggr).\end{aligned}$$ Known in [@dittrich2000probing], one has $$\begin{aligned}
&\mathsf{P}=\cos q_{f}s\cosh q_{f}s,\, \mathsf{Q}=\sin q_{f}s\sinh q_{f}s \nonumber\\
&\mathsf{R}=(\sinh q_{f}s\cos q_{f}s+\cosh q_{f}s\sin q_{f}s)/2,\\
&\mathsf{S}=(\sinh q_{f}s\cos q_{f}s-\cosh q_{f}s\sin q_{f}s)/2.\nonumber\end{aligned}$$ Then, we perform an approximate expansion $$\begin{aligned}
\label{eqn_ap2bp}
&{\mathop{}\!e}^{ p\cdot A(s)\cdot p+B(s)\cdot p}\simeq
{\mathop{}\!e}^{p\cdot A(s)\cdot p}{\mathop{}\!e}^{B(s)\cdot p}
\\
&={\mathop{}\!e}^{ p\cdot A(s)\cdot p}{\left}(\cos\varrho
+B(s) \cdot p\frac{\sin\varrho}{\varrho}{\right})\nonumber
 \end{aligned}$$ where $\varrho=2q\pi_{1}k/f$ and $k={\left}(\mathsf{p}_{1}\cdot\mathsf{p}_{2}{\right})^{\frac{1}{2}}$. Now, it is allowed us to integrate with respect to $p$ and take the trace in the spinor space. With help of the [Eq. (\[eqn\_O\_exp\])]{}, [Eq. (\[eqn\_F\_exp\])]{} and [Eq. (\[eqn\_ap2bp\])]{}, one has $$\begin{aligned}
L(s)&={\hbox{tr}}\int{\mathop{}\!e}^{\frac{1}{2f}{\sigma}{\lambda}s}{\mathop{}\!e}^{p\cdot A(s)\cdot p+B(s)\cdot p}{\mathop{}\!e}^{-\frac{1}{4f}{\left}[{\sigma}{\lambda}s,B(s) \cdot p{\right}]}{\mathop{}\!d}^{4}p \nonumber\\
&=L_{0}(s)+L_{1}(s)+L_{2}(s).\end{aligned}$$ Here ${\langle}X{\rangle}$ denotes integrating in momentum and tracing in spinor space ${\hbox{tr}}\int X{\mathop{}\!e}^{p\cdot A\cdot p}{\mathop{}\!d}^{4}p$. It gives $$\begin{aligned}
&L_{0}(s)={\left}{\langle}\cos\varrho\mathsf{T}\mathsf{P}{\right}{\rangle}, \, \,
L_{1}(s)={\left}{\langle}\cos\varrho\mathsf{U}\mathsf{Q}{\right}{\rangle}, \, \, L_{2}(s)=\\
&={\left}< \frac{2\tilde{q}_{f}\sin\varrho}{K^{2}k}
{\left}[q_{f}K^{2}{\left}(\mathsf{W}\mathsf{S}-\mathsf{V}\mathsf{R}{\right})+\tilde{q}_{f}k^{2}{\left}(\mathsf{V}\mathsf{S}+\mathsf{W}\mathsf{R}{\right}){\right}]{\right}>. \nonumber\end{aligned}$$ The integration with respect to momentum $p$ is in the Gaussian form, which can be taken easily with result $$\begin{aligned}
&{\langle}1{\rangle}=\mathcal{N}=\pi^{2}{\mathrm}{Det}A^{-\frac{1}{2}},\,A={\mathrm}{Diag}{\left}(a_{\shortparallel}, a_{\perp}, a_{\perp}, a_{\shortparallel}{\right})\nonumber\\
&{\langle}K^{2}{\rangle}=\frac{\mathcal{N}}{2}\,{\hbox{tr}}{\left}(\frac{D_{+}}{A}{\right}),\,
{\langle}k^{2}{\rangle}=\frac{\mathcal{N}}{2}\,{\hbox{tr}}{\left}(\frac{D_{-}}{A}{\right})\end{aligned}$$ The matrices $D_{\pm}={\mathrm}{Diag}{\left}(b^{2}_{\shortparallel}, \pm b^{2}_{\perp}, \pm b^{2}_{\perp}, b^{2}_{\shortparallel}{\right})$, which read from [Eq. (\[eqn\_b03\])]{} and [Eq. (\[eqn\_b12\])]{}. From [Eq. (\[eqn\_sol\_ABC\])]{}, one has $a_{\shortparallel}=\tan q_{f}s/{\left}(q_{f} f{\right})$ and $a_{\perp}=\tanh q_{f}s/{\left}(q_{f} f{\right})$. The higher orders corrections ${\langle}K^{4}{\rangle}$, ${\langle}k^{4}{\rangle}$ and $ {\langle}k^{2}K^{2}{\rangle}$ can be drawn in a similar manner, which are abbreviated here. Since ${\theta}p\sim\pi_{1}ps/f\sim \pi_{1}p/{\Lambda}^{2}\ll 1$ and the integration is exponential suppressed for large $s$, it enables us to approximate $\sin(a{\theta}p)$, $\sinh(a{\theta}p)\sim a{\theta}p$ and $\cos(a{\theta}p)$, $\cosh(a{\theta}p)\sim 1$. Hence, it acquires $\mathsf{T}\sim 1$, $\mathsf{U}\sim k^{2}s$, $\mathsf{V}\sim K^{2}s$ and $\mathsf{W}\sim  K^{2}k^2 s$. Finally, take the integration with respect to $s$ to get $$\begin{aligned}
\label{eqn_eff_potentioal}
&\mathcal{S}_{eff}= {\mathcal}{S}_{eff}^{0}+{\mathcal}{S}_{eff}^{1}+{\mathcal}{S}_{eff}^{2}, \\
&{\mathcal}{S}_{eff}^{i}=\frac{N_c }{4\pi^{2}}{\hbox{Tr}}\int_{f/{\Lambda}^{2}}^{\infty} \frac{{\mathop{}\!d}s}{2s} {\mathop{}\!e}^{-h(s)} S_{eff}^{i}(s),\quad \nonumber \end{aligned}$$ where $-h(s)=-(m^{2}+\pi_{1}^{2})s/f+{C}(s)-\frac{1}{2}\ln{\hbox{tr}}{A}$, and $$\begin{aligned}
&{C}(s)-\frac{\ln{\hbox{tr}}{A}}{2}=
-\ln\frac{\sin q_{f}s\sinh q_{f}s}{q_{f}^{2}f^{2}}\nonumber\\
&\quad\quad\quad-\frac{2q^{2}{\pi_{1}}^{2}}{\tilde{q}_{f}^{3}f}{\left}(2\tilde{q}_{f}s-\tan \tilde{q}_{f}s-\tanh \tilde{q}_{f}s{\right}).\end{aligned}$$ The detailed integrands $S_{eff}^{i}(s)$ are $$\begin{aligned}
&S_{eff}^{0}(s) = \mathsf{P},\quad \nonumber \\
&S_{eff}^{1}(s) =4 \tau_{3} \frac{Q q^3\pi_{1}^2 s^2}{f^2{\mathcal}{N}}\frac{}{}{\left}{\langle}k^2{\right}{\rangle}\mathsf{Q} ,\quad \label{exprS2} \\
&S_{eff}^{2}(s) =\frac{8\tilde{q}_{f}q^2\pi_{1}^2 s}{f^2{\mathcal}{N}}
{\left}(-q_{f}{\left}{\langle}K^2{\right}{\rangle}\mathsf{R}+\tilde{q}_{f}{\left}{\langle}k^2{\right}{\rangle}\mathsf{S}{\right}) 
+\nonumber\\
&\quad\quad+\tau_{3} \frac{32\tilde{q}_{f} Qq^5 \pi_{1}^4 s^3}{3f^4{\mathcal}{N}}
{\left}(q_{f}{\left}{\langle}K^2k^2{\right}{\rangle}\mathsf{S}+\tilde{q}_{f}{\left}{\langle}k^4{\right}{\rangle}\mathsf{R}{\right}).\nonumber  \end{aligned}$$

Eventually, we have the effective potential which takes the following form: $$\begin{aligned}
\Omega = \frac{(m-m_{0})^2+\pi_1^2}{4\mathrm{G}} + \mathcal{S}_{eff}.\end{aligned}$$

Numerical results
=================

The NJL model is nonrenormalizable and therefore the UV cut-off should be employed in order to get reasonable results, where a proper time regularization is applied in the work, i.e., the integration with respect to $s$ start from $f/{\Lambda}^{2}$. We perform calculations of integral expression for $\mathcal{S}_{eff}$ in Eq. (\[eqn\_eff\_potentioal\]) numerically. In the limit of zero field $f$ the expression leads to the usual proper-time regularization scheme of NJL model. Therefore, for numerical estimation we use the model parameterization from Ref. [@Inagaki:2015lma]. Namely, in [@Inagaki:2015lma] there are five sets of model parameters for proper-time regularization scheme which are fitted in favor of observable values of pion mass and weak pion decay constant. For convenience we present them in Table 1. In the set 1 the constituent quark mass is $178$ MeV and for set 5 is $372$ MeV. The constituent quark masses for other parameterizations are in between these two cases. Therefore, one can consider set 1 and set 5 as limited cases for the predictions of the NJL model.

   Set   $m_0$\[MeV\]   $\Lambda$\[MeV\]   G\[GeV$^{-2}$\]   $m$\[MeV\]        
  ----- -------------- ------------------ ----------------- ------------ -- -- --
    1        3.0              1464              1.61            178            
    2        5.0              1097              3.07            204            
    3        8.0              849               5.85            245            
    4        10.0             755               8.13            265            
    5        15.0             645               17.2            372            

  : Parameters of the NJL model in the proper-time regularization taken from [@Inagaki:2015lma].

\[TableParameters\]

The important point of calculation is that integrand of $\mathcal{S}_{eff}$ contain singularities and one should specify how to deal with them. The singularities which are generated by trigonometric functions tangent and cotangent of $q_i s$ for quark flavor i are located at real axis and by hyperbolic functions at imaginary axis. We shift $s$ to the complex plane $s+i\epsilon$, see Fig. \[ContTikz\], since we prefer to running a numerical calculation of integral instead of residues summation like what used in [@Inagaki:2003ac; @Ruggieri:2016xww]. In principle, the effective potential at finite $f$ acquires an imaginary part which correspond to pair-production because of Schwinger mechanism [@Schwinger:1951nm; @Tavares:2018poq; @Cao:2015dya]. We figure out that the imaginary part is smaller than the real part in current work. Plus, the subtle effect of Schwinger mechanism is out of the scope of the present paper and will not discuss here.

![Contour on complex $s$-plane. Singularities for a quark of flavor i which are related to tangent are shown by open circles while filled circles correspond to those of cotangent $\cot(q_i s)$.[]{data-label="ContTikz"}](ContTikz2){width="47.00000%"}

In Figs. \[Scan5ff001\], \[Scan5ff010\], \[Scan5ff035\], the behavior of effective potential for Set 5 of model parameters is plotted for field values $f=0.01,0.2,0.450$ GeV$^2$, respectively. We found the following typical behavior for three regions: 1) For small field $f=0.01$ GeV$^2$ as shown in Fig. \[Scan5ff001\], the system is in usual (almost vacuum) chiral symmetry breaking phase with nonzero sigma condensate and zero pion condensate; 2) For moderate field $f=0.2$ GeV$^2$, seen in Fig. \[Scan5ff010\], the additional minima appears in the effective potential and the system takes a chiral rotation in $\sigma -\pi_1$ plane to have a nonzero pion condensate, $\pi_1$; 3) For large field $f=0.450$ GeV$^2$, read from Fig. \[Scan5ff035\], the minimum with $\pi_1=0$ is energetically favorable.

There are two sources to break the chiral symmetry: spontaneous chiral symmetry breaking due to presence of quark condensate ${\langle}\bar{\psi} \psi {\rangle}$ and explicit chiral symmetry breaking due to nonzero current quark mass in the Lagrangian. Therefore, we investigate not only the reality situation but also for $m_0\to 0$. To systematically perform this task, we vary $m_{0}$ and recalculate $m$ while $\mathrm{G}$ and $\Lambda$ have the same values, i.e. we consider $m$ as a function of $m_{0}$ [@Bernard:1992mp]. In the following we denote the physical value of current quark mass as $m_{0}^\star$. The behaviors of $m$ and $\pi_1$ as a function of field $f$ are presented in Fig. \[Set1Set5MassDelta\] for different values of ratio $m_{0}/m_{0}^\star=0.01,0.1,0.5,1.0$. The left and right sides are obtained by model parameter sets 1 and 5, respectively. It is straightforward to figure out that for small current quark mass the system is more preferable to chirally rotate from zero to nonzero value $\pi_1$, leaving the total order parameter of chiral symmetry breaking $|M|=\sqrt{m^2+\pi_1^2}$ unchanged. With increasing of $m_{0}$ the situation becomes more complicated. The phase of pion condensation even never show up for $m_{0}=m_{0}^\star$ in the model parameter Set 1.

Conclusions {#con}
===========

In this paper the charged pion condensation under the parallel electromagnetic fields is calculated in the framework of the NJL model by using Schwinger proper-time method. The configuration of field is chosen, the electric field being anti-parallel to the magnetic one, to have a zero first Lorentz invariant, $I_{1}=\mathbf{E}^{2}-\mathbf{B}^{2}$, and a nonzero second Lorentz invariant, $I_2=\mathbf{E} \cdot \mathbf{B}$.

We find that in the chiral limit the system is favorable to form a both nonzero condensation of scalar and charged pion, i.e. rotating in the chiral group. Chiral condensates aligning to pseudo mesons space has been found in [@Cao:2015cka] by the methods of $\chi{\mathrm}{PT}$ and NJL model, where the system is immediately straighten up $\pi_0$ direction in the chiral limit once the second Lorentz invariant $I_2$ turned on. The main difference of charged condensation is that the system will across a weakly first order phase transition to zero pion condensate and then a second order phase transition to chirally symmetric phase as the field strength increasing, while it, characterizing by $\pi_{0}$, is a whole second order phase transition as shown in [@Cao:2015cka]. The underlying mechanism are two folds. One is the obviously coupling between charged pions and electromagnetism. Another reason is that a more complicated influence of anomalous diagrams are implicitly included, not only $\pi_0\rightarrow\gamma\gamma$ but also $\gamma\rightarrow \pi_+\pi_-\pi_0$.

Indeed, if assuming condensation in the neutral channel ${\langle}\sigma{\rangle}$ nears a second order phase transition, its effective potential has the form ${\mathcal}{S}_{eff}^{0}\sim -c_{0}M^{2}+c_{1}M^{4}/f$ according to Ginzburg-Landau theory [@Ginzburg:1950sr]. However, if we include $\pi_{\pm}$ as an additional degree of freedom and non-degenerate with $\pi_{0}$, read from [Eq. (\[eqn\_eff\_potentioal\])]{}, the potential arranges as: ${\mathcal}{S}_{eff}^{2}\sim -\tilde{c}_{1}M^{4}/f+c_{2}M^{6}/f^{2}$. As a result, we have a weakly first order phase transition and effective potential in the form of $$\begin{aligned}
    {\Omega}=\frac{M^{2}}{4 \mathrm{G}}-c_{0}M^{2}+\frac{{\left}(c_{1}-\tilde{c}_{1}{\right})M^{4}}{f}+\frac{c_{2}M^{6}}{f^{2}}.\end{aligned}$$ Our numerical simulations support these arguments, read from Figs. \[Scan5ff001\], \[Scan5ff010\], \[Scan5ff035\]. The mass of current quarks plays an important role and it denies our claim at some regions of the model parameters. It requires a further study via the first principle calculation, such as Dyson-Schwinger equation or functional renormalization group methods.

Application of the charged pion condensation to the case heavy-ion collisions or neutron stars interior need an extension to finite temperature and/or chemical potential. We will explore this extension in future.

Acknowledgments {#Ackn}
===============

We are grateful to Maxim Chernodub, Nikolai Kochelev, Marco Ruggieri and Pengming Zhang for the useful discussions. J.Y.C. is supported by the NSFC under Grant number: 11605254 and Major State Basic Research Development Program in China (No. 2015CB856903). M.H. is supported by the NSFC under Grant No. 11725523, 11735007 and 11261130311(CRC 110 by DFG and NSFC). A.R. is supported by the CAS President’s international fellowship initiative (Grant No. 2017VMA0045), Council for Grants of the President of the Russian Federation (project NSh-8081.2016.9) and numerical calculations are performed on computing cluster “Akademik V.M. Matrosov” (http://hpc.icc.ru).

![image](Slice5ff010a){width="49.00000%"} ![image](SCan5ff010a){width="43.00000%"}

![image](Slice5ff200a){width="49.00000%"} ![image](SCan5ff200wa){width="43.00000%"}

![image](Slice5ff450a){width="49.00000%"} ![image](SCan5ff450wa){width="43.00000%"}

![image](Mult1Set1V3){width="47.00000%"} ![image](Mult1Set5V3){width="47.00000%"}

[10]{}

J. Cleymans, R. V. Gavai, E. Suhonen, Phys. Rept. [**130**]{} (1986) 217. M. Buballa, Phys. Rept. [**407**]{} (2005) 205. M. Buballa, S. Carignano, Prog. Part. Nucl. Phys. [**81**]{} (2015) 39. J. O. Andersen, P. Kneschke, [](http://arxiv.org/abs/1802.01832).

D. E. Kharzeev, L. D. McLerran, H. J. Warringa, Nucl. Phys. A [**803**]{} (2008) 227. V. Skokov, A. [Yu]{}. Illarionov, V. Toneev, Int. J. Mod. Phys. A [**24**]{} (2009) 5925. K. Hattori, X.-G. Huang, Nucl. Sci. Tech. [**28**]{} (2) (2017) 26. J. O. Andersen, W. R. Naylor and A. Tranberg, Rev. Mod. Phys.  [**88**]{} (2016) 025001. V. A. Miransky and I. A. Shovkovy, Phys. Rept.  [**576**]{} (2015) 1. X. G. Huang, Rept. Prog. Phys.  [**79**]{}, no. 7 (2016) 076302. W. -T. Deng and X. -G. Huang, Phys. Rev. C [**85**]{} (2012) 044907. R. C. Duncan, C. Thompson, Astrophys. J. [**392**]{} (1992) L9. D. Blaschke, N. Chamel, [](http://arxiv.org/abs/1803.01836).

S. P. Klevansky and R. H. Lemmer, Phys. Rev. D [**39**]{} (1989) 3478. K. G. Klimenko, Theor. Math. Phys.  [**89**]{} (1992) 1161 \[Teor. Mat. Fiz.  [**89**]{} (1991) 211\]. V. P. Gusynin, V. A. Miransky and I. A. Shovkovy, Nucl. Phys. B [**462**]{} (1996) 249. V. P. Gusynin, V. A. Miransky and I. A. Shovkovy, Nucl. Phys. B [**563**]{} (1999) 361. G. S. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. D. Katz, S. Krieg, A. Schafer and K. K. Szabo, JHEP [**1202**]{} (2012) 044. G. S. Bali, F. Bruckmann, G. Endrodi, Z. Fodor, S. D. Katz and A. Schafer, Phys. Rev. D [**86**]{} (2012) 071502. G. S. Bali, F. Bruckmann, G. Endrodi, F. Gruber and A. Schaefer, JHEP [**1304**]{} (2013) 130. D. Kharzeev and A. Zhitnitsky, Nucl. Phys. A [**797**]{} (2007) 67. K. Fukushima, D. E. Kharzeev and H. J. Warringa, Phys. Rev.  D [**78**]{} (2008) 074033 . M. N. Chernodub, Phys. Rev. D [**82**]{} (2010) 085011. M. N. Chernodub, Phys. Rev. Lett.  [**106**]{} (2011) 142003 . G. Cao, X.-G. Huang, Phys. Lett. B [**757**]{} (2016) 1. L. Wang, G. Cao, X. G. Huang and P. Zhuang, Phys. Lett. B [**780**]{} (2018) 273.

S. L. Adler, [Anomalies to all orders]{}, in: G. ’t Hooft (Ed.), 50 years of Yang-Mills theory, 2005, pp. 187–228. R. F. Sawyer, Phys. Rev. Lett. [**29**]{} (1972) 382. R. F. Sawyer, D. J. Scalapino, Phys. Rev. D [**7**]{} (1973) 953. D. N. Voskresensky, N. [Yu]{}. Anisimov, Sov. Phys. JETP [**51**]{} (1980) 13 \[Zh. Eksp. Teor. Fiz. [**78**]{} (1980) 28\].

A. B. Migdal, Rev. Mod. Phys. [**50**]{} (1978) 107. D. T. Son and M. A. Stephanov, Phys. Rev. Lett.  [**86**]{} (2001) 592. L. He, M. Jin and P. Zhuang, Phys. Rev. D [**74**]{} (2006) 036005. S. Mao, Phys. Rev. D [**89**]{} (2014) 116006. T. G. Khunjua, K. G. Klimenko, R. N. Zhokhov, V. C. Zhukovsky, Phys. Rev. D [**95**]{} (2017) 105010. B. B. Brandt, G. Endrodi, PoS [**LATTICE2016**]{} (2016) 039.

B. B. Brandt, G. Endrodi, E. S. Fraga, M. Hippert, J. Schaffner-Bielich, S. Schmalzbauer, [](http://arxiv.org/abs/1802.06685).

Y. Nambu, G. Jona-Lasinio, Phys. Rev. [**122**]{} (1961) 345. Y. Nambu, G. Jona-Lasinio, Phys. Rev. [**124**]{} (1961) 246. J. S. Schwinger, Phys. Rev. [**82**]{} (1951) 664. M. K. Volkov, Sov. J. Part. Nucl. [**17**]{} (1986) 186, \[Fiz. Elem. Chast. Atom. Yadra [**17**]{} (1986) 433\].

U. Vogl, W. Weise, Prog. Part. Nucl. Phys. [**27**]{} (1991) 195. S. P. Klevansky, Rev. Mod. Phys. [**64**]{} (1992) 649. T. Hatsuda, T. Kunihiro, Phys. Rept. [**247**]{} (1994) 221. M. K. Volkov, A. E. Radzhabov, Phys. Usp. [**49**]{} (2006) 551. M. R. Brown, M. J. Duff, Phys. Rev. D [**11**]{} (1975) 2124. W. Dittrich and H. Gies, [*Probing the quantum vacuum: perturbative effective action approach in quantum electrodynamics and its application*]{} (2000), Springer Science & Business Media.

H. Kohyama, D. Kimura, T. Inagaki, Nucl. Phys. B [**896**]{} (2015) 682. T. Inagaki, D. Kimura, T. Murata, Int. J. Mod. Phys. A [**20**]{} (2005) 4995. M. Ruggieri, Z. Y. Lu, G. X. Peng, Phys. Rev. D [**94**]{} (2016) 116003.

W. R. Tavares, S. S. Avancini, [](http://arxiv.org/abs/1801.10566).

G. Cao, X.-G. Huang, Phys. Rev. D [**93**]{} (2016) 016007. V. Bernard, A. A. Osipov, U. G. Meissner, Phys. Lett. B [**285**]{} (1992) 119. V. L. Ginzburg, L. D. Landau, Zh. Eksp. Teor. Fiz. [**20**]{} (1950) 1064.
